# Seasonal thaws under mid-to-low pressure atmospheres on Early Mars

Paolo Simonetti INAF - Trieste Astronomical Observatory, Via G. B. Tiepolo 11,
34143 Trieste, Italy  Giovanni Vladilo INAF - Trieste Astronomical
Observatory, Via G. B. Tiepolo 11, 34143 Trieste, Italy  Stavro L. Ivanovski
INAF - Trieste Astronomical Observatory, Via G. B. Tiepolo 11, 34143 Trieste,
Italy  Laura Silva INAF - Trieste Astronomical Observatory, Via G. B. Tiepolo
11, 34143 Trieste, Italy  IFPU - Institute for Fundamental Physics of the
Universe, Via Beirut 2, 34014 Trieste, Italy  Lorenzo Biasiotti INAF - Trieste
Astronomical Observatory, Via G. B. Tiepolo 11, 34143 Trieste, Italy
University of Trieste - Dep. of Physics, Via G. B. Tiepolo 11, 34143 Trieste,
Italy  Michele Maris INAF - Trieste Astronomical Observatory, Via G. B.
Tiepolo 11, 34143 Trieste, Italy  Giuseppe Murante  INAF - Trieste
Astronomical Observatory, Via G. B. Tiepolo 11, 34143 Trieste, Italy  IFPU -
Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014
Trieste, Italy  CNR - Institute of Geosciences and Georesources, Via Giuseppe
Moruzzi, 1, 56127 Pisa, Italy  Erica Bisesi  INAF - Trieste Astronomical
Observatory, Via G. B. Tiepolo 11, 34143 Trieste, Italy  IFPU - Institute for
Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy  CNR -
Institute of Geosciences and Georesources, Via Giuseppe Moruzzi, 1, 56127
Pisa, Italy  Sergio Monai  INAF - Trieste Astronomical Observatory, Via G. B.
Tiepolo 11, 34143 Trieste, Italy

###### Abstract

Despite decades of scientific research on the subject, the climate of the
first 1.5 Gyr of Mars history has not been fully understood yet. Especially
challenging is the need to reconcile the presence of liquid water for extended
periods of time on the martian surface with the comparatively low insolation
received by the planet, a problem which is known as the Faint Young Sun (FYS)
Paradox. In this paper we use ESTM, a latitudinal energy balance model with
enhanced prescriptions for meridional heat diffusion, and the radiative
transfer code EOS to investigate how seasonal variations of temperature can
give rise to local conditions which are conductive to liquid water runoffs. We
include the effects of the martian dichotomy, a northern ocean with either 150
or 550 m of Global Equivalent Layer (GEL) and simplified CO2 or H2O clouds. We
find that 1.3-to-2.0 bar CO2-dominated atmospheres can produce seasonal thaws
due to inefficient heat redistribution, provided that the eccentricity and the
obliquity of the planet are sufficiently different from zero. We also studied
the impact of different values for the argument of perihelion. When local
favorable conditions exist, they nearly always persist for $>15\%$ of the
martian year. These results are obtained without the need for additional
greenhouse gases (e.g. H2, CH4) or transient heat-injecting phenomena (e.g.
asteroid impacts, volcanic eruptions). Moderate amounts (0.1 to 1%) of CH4
significantly widens the parameter space region in which seasonal thaws are
possible.

Mars (1007) — Astrobiology (74) — Planetary climates (2184) — Planetary
atmospheres (1244)

##  1 Introduction

The early history of the planet Mars is divided into three periods called
Noachian, Hesperian and Amazonian (see e.g. Hartmann & Neukum, 2001), each of
them linked to a set of geographical regions grouped on the basis of density,
size and superposition of craters. A fourth period, called pre-Noachian,
precedes the most ancient terrains identified so far. The chronological
intervals usually associated to these periods are 4.5 - 4.1 (pre-Noachian),
4.1 - 3.7 (Noachian), 3.7 - 3.0 Gyr ago (Hesperian), with the Amazonian
continuing up to present day (see e.g. Carr & Head, 2010).

A wide range of geomorphological and chemical evidences suggests that, at some
point in the past, liquid water was present on the surface of Mars for
extended (at least $10^{5}$ \- $10^{7}$ yr) periods of time (Hoke et al.,
2011; Balme et al., 2020) and in sufficient quantities to sustain a vigorous
hydrological cycle (e.g. Mangold et al., 2004). Since most of the dendritic
valley networks are mainly associated with Noachian and, to a lesser extent,
Hesperian terrain, the humid period of the martian history is identified with
the Noachian and Hesperian ages ($\sim$4.1 - 3.0 Gyr ago). However, at that
epoch, the power output of the Sun was most probably lower than today by
$20-25\%$ (Gough, 1981), making it difficult to explain how Mars could have
hosted climate conditions sufficiently warm to allow for such an hydrological
cycle to be present. This longstanding problem is known as the Faint Young Sun
(FYS) paradox.

The study of planetary climates is performed using a variety of different
climate models. These models are usually categorized depending on their
complexity, which scales to the number of dimensions assigned to the planet
and its atmosphere. The simplest ones are called Energy Balance Models (EBMs)
and can be either 1D or 2D. 1D EBMs derive the climatology of the entire
planet from a single atmospheric column in radiative-convective equilibrium
(RCE) with fixed surface albedo and thus focus on the radiative transfer (RT)
properties of a given atmospheric composition (Pierrehumbert, 2010). 2D EBMs
add the latitudinal horizontal direction and, as such, are able to produce
seasonal-latitudinal temperature profiles (North et al., 1981). The inclusion
of a simplified horizontal diffusion of heat along the surface allows to take
into consideration some simple climate feedbacks, like the ice-albedo
feedback. Recently, a number of 3D EBMs, that consider also the longitudinal
diffusion of heat, have been proposed for the study of tidally-locked
exoplanets (Haqq-Misra & Hayworth, 2022). At the upper end of the complexity
spectrum there are the General Circulation Models (GCMs), which are coupled
radiative transfer-fluidodynamical models that integrate the Navier-Stokes
equations to reconstruct the full 3D motions of a planetary atmosphere (e.g.
Richardson et al., 2007). The required computational power is proportional to
the model complexity, which makes GCMs unsuitable for parameter space
explorations.

In the last decades, many possible solutions to the FYS paradox have been
proposed, that can be generally classified in two large groups. The first one
searches for an atmospheric composition and thickness capable to warm the
surface up to the required temperature. Early calculations focused on the
viability of a thick CO2-H2O atmosphere, but in 1D EBMs this is both
insufficient to achieve the goal (Kasting, 1991) and difficult to reconcile
with observational data concerning the surface paleopressure of Mars (Lammer
et al., 2013; Kite et al., 2014; Edwards & Ehlmann, 2015, but see also Bultel
et al., 2019). More recent investigations took into consideration the
greenhouse effect produced by other chemical species such as SO2/H2S (Halevy
et al., 2007; Johnson et al., 2008; Mischna et al., 2013) or H2/CH4 (Ramirez
et al., 2014; Wordsworth et al., 2017; Ramirez, 2017; Haberle et al., 2019;
Turbet et al., 2020a), satisfactorily solving the issues associated with
CO2-only atmospheres, but creating new ones related to the availability of
these additional gases (see e.g. Wordsworth, 2016). The role of infrared
backscattering from CO2 clouds to heat the planet has also been investigated
(Forget & Pierrehumbert, 1997), but later calculations cast doubts on its
effectiveness (Colaprete & Toon, 2003; Kitzmann et al., 2013; Forget et al.,
2013), since a vigorous greenhouse effect requires a very specific and narrow
distribution of cloud particle sizes. Similar conclusions were drawn
concerning the radiative forcing of H2O clouds (Wordsworth et al., 2013; Urata
& Toon, 2013; Ramirez & Kasting, 2017). Variations of the martian surface
albedo were also taken into consideration but discarded (Fairén et al., 2012).
The second group of proposed possibilities focuses on the role of transient
warmings as solutions to the FYS paradox, such as asteroid impacts (Segura et
al., 2002, 2008), volcanic eruptions (Halevy & Head, 2014), the
destabilization of previously formed CH4 clathrates due to chaotic obliquity
changes (Kite et al., 2017) or episodic injections of H2 due to crustal
alteration (Wordsworth et al., 2021). None of them, however, are free of
issues. The post-impact rainfall pattern predicted by more modern models
mismatches observations (Turbet et al., 2020b); photochemically driven
formation of highly reflective sulfate aerosols from volcanic SO2 and H2S may
have cooled, rather than warmed, the planet (Tian et al., 2010; Kerber et al.,
2015); it is uncertain if CH4 could have been produced and stored in
sufficient quantities in the martian environment in the first place (Ramirez &
Craddock, 2018).

The two sets of solutions proposed for the FYS paradox give rise to two
possible Early Mars climatologies, called ”warm-and-wet” (e.g. Craddock &
Howard, 2002) and ”cold-and-icy” (e.g. Squyres & Kasting, 1994). The warm-and-
wet scenario is associated with a hydrosphere which rests mostly in a
relatively wide ocean covering a considerable portion of the northern lowlands
(plus the Hellas and Argyre basins in the south), with minor and/or transient
glaciers on the southern highlands. This scenario offers the most
straightforward explanation for valley networks formation and is in accordance
with the existence of a peak in the altitude distribution of fluvial delta-
like formations corresponding to the paleo-shoreline of the northern ocean (di
Achille & Hynek, 2010). This scenario however, has been criticized because it
does not seem to be able to produce the precipitation pattern required to form
the observed valley networks in the right places (see e.g. in Wordsworth et
al., 2015). Also, in order to prevent the locking of the hydrosphere in the
southern highlands (which act as efficient cold traps), an even higher global
average temperature is required, which puts an extra constraint on the
atmospheric chemical composition. The cold-and-icy scenario envisions instead
a situation where most of the water reservoirs are locked in glaciers and
episodic melting events provide the liquid water runoffs necessary to explain
surface erosion. In the case of a large ($\gtrsim 200$ m) global equivalent
layer (GEL) of water, significant basal melting can occur at the bottom of
highland ice sheets, allowing for efficient fluvioglacial erosion (Fastook &
Head, 2014). However, the lack of an extensive glacial erosion in the regions
where these ice sheets should have been existed is a major unsolved issue of
this picture.

The duration and nature of the conditions allowing for liquid water on the
martian surface have important astrobiological implications, too. The possible
survival of martian life in subsurface environments up to the present day,
which is occasionally invoked to explain seasonal variabilities in the martian
atmospheric composition (Webster et al., 2015), hinges on the existence of a
climate state conductive of the emergence of life in the first place. While
the conditions for abiogenesis are still unknown, early Mars is believed to
hold a widespread range of physical and chemical properties potentially
suitable for an independent origin of life, with different possible niches
depending on the prevalent climate state (e.g. Clark et al., 2021).

Most of Early Mars modelings have been performed either via parameter space
explorations with 1D EBMs, or by running very small and specific sets of cases
with 3D GCMs. This is sound considering the respective computational costs of
these models and yielded a vast amount of insights on the surface conditions
during the Noachian and Hesperian periods. Efficient 1D EBMs have been coupled
with a variety of other codes, modeling processes such as the atmospheric
photochemistry (Batalha et al., 2015), the carbonate-silicate cycle (Batalha
et al., 2016) or ecological interactions (Sauterey et al., 2022). On the other
hand, GCMs have been fundamental to cast light on the precipitation pattern
and ice deposition (e.g. Wordsworth et al., 2015; Kamada et al., 2020, 2021).
In this work, we make the case for the employment of 2D EBMs coupled with a
flexible RT code as the natural next step along the path inaugurated by 1D
models by studying the effects of seasonal surface temperature variations,
especially under the influence of specific choices regarding the orbital
elements of Early Mars (eccentricity, obliquity and the argument of
perihelion). In particular, we use EOS-ESTM (Biasiotti et al., 2022) to find
the maximum and minimum latitudinal band-averaged temperatures to see if a
given set of atmospheric compositions (including CO2, H2O and CH4) and orbital
parameters are able to produce seasonal melting of surface ices without
concurrently causing seasonal condensation of the atmosphere at surface. We
also tested different hypotheses in terms of surface albedo and cloud
radiative effects. Using 2D EBMs is not entirely new (see e.g. Hayworth et
al., 2020) but it is still largely unexplored and promises to yield new
results by taking the best of the two worlds, i.e. retaining a high
computational efficiency while allowing a better understanding of the
conditions that might arise on the planetary surface.

The structure of this paper is as follows. In Section 2 we describe our
climate (ESTM) and vertical radiative transfer (EOS) models, together with the
specific choices we made to take into consideration the inhomogeneous water
and altitude distribution of Mars. In Section 3 we show the results of the
EOS-ESTM parameter space explorations. In Section 4 we discuss our results, in
particular comparing them with some bulk observational evidences and in
Section 5 we draw our conclusions.

##  2 Model description

###  2.1 ESTM

Parameter | Description | Value
---|---|---
$a$ | semi-major axis | 1.524 AU
$e$ | eccentricity | 0.0934
$I$ | insolation | 439.5 W m-2
$\varepsilon$ | obliquity | 25∘
$\omega$ | argument of perihelion | 286.5∘
$P_{\text{rot}}$ | rotation period | 1.0275 d
$P_{\text{orb}}$ | orbital period | 686.96 d
$g$ | surface gravity acceleration | 3.72 m s-1
$R$ | mean radius | 3389 km
Table 1: Orbital and planetary parameters adopted for the Early Mars in this
work when not otherwise specified.

The Earth-like planet surface temperature model (ESTM) is a largely upgraded
version of a latitudinal and seasonal Energy Balance Model (EBM). In common
with classic EBMs (North et al., 1981), the meridional transport is treated
with a formalism of heat diffusion. In practice, in each latitude zone the
energy balance is described by the equation

$C\frac{\partial T}{\partial t}-\frac{\partial}{\partial
x}\left[D\,(1-x^{2})\,\frac{\partial T}{\partial x}\right]+I=S\,(1-A),$ (1)

where $T$ is the zonal temperature mediated over one rotation period and all
terms are normalized per unit area. Most of the coefficients depend on both
time, $t$, and latitude, $\varphi$, either directly or indirectly, through
their dependence on $T$. The output of ESTM is $T$ as a function of $\varphi$
and the orbital phase. The first term of the equation represents the seasonal
evolution of the zonal heat storage; $C$ is the zonal heat capacity. The
second term represents the amount of heat per unit time leaving the zone along
the meridional direction; $D$ is the “diffusion term” and $x=\sin\varphi$. The
term $I$ is the thermal radiation emitted by the zone, also called Outgoing
Longwave Radiation (OLR). The right hand of the equation represents the zonal
heating due to stellar photons; $S$ is the incoming stellar radiation and $A$
the albedo at the top of the atmosphere (TOA albedo).

At variance with classic EBMs, where $D$ is a constant, ESTM incorporates a
physically-based description of $D$ (Vladilo et al., 2015). If we define
$\Phi$ in such a way that $2\pi R^{2}\Phi\cos\varphi$ is the net rate of
energy transport across a circle of constant $\varphi$ in a planet of radius
$R$ (Pierrehumbert, 2010), the coefficient $D$ can be derived as a function of
physical quantities from the analogy with the equation of heat diffusion, i.e.
$\Phi\equiv-D(\partial T/\partial\varphi)$. By introducing an analytical
description of the longitudinally averaged atmospheric transport validated
with 3D climate models, $D$ is eventually calculated as a function of $R$,
surface atmospheric pressure, $P$, surface gravity acceleration, $g$, angular
rotation rate, $\Omega$, and other planetary properties (Vladilo et al.,
2015).

Another substantial improvement of ESTM with respect to classic EBMs concerns
the description of the TOA albedo, $A$, and of the OLR, $I$. In classic EBMs,
$A$ is a constant and $I$ a linear function of $T$. In ESTM, $A$ and $I$ are
calculated by taking into account the atmospheric radiative transfer in the
short and long wavelength ranges, respectively. In practice, this is achieved
by interpolating on previously calculated lookup tables in which $A$ and $I$
for a specific atmospheric composition and pressure are stored as a function
of the surface temperature $T$ and other relevant quantities. These
calculations are performed for clear-sky atmospheres of given chemical
composition and surface pressure, making use of the EOS code described below.
Clouds are then incorporated in ESTM through a parameterization of their
impact on the albedo and on the OLR. Starting from a realistic description of
the surface albedo, the TOA albedo is calculated by taking into account the
atmospheric transport in the short wavelength range and the simplified
treatment of the cloud albedo. Both the cloud and the surface albedos are
dependent on $\mu=\cos z$, where $z$ is the zenith angle of the incident
light.

Apart for the parameters that we explicitly changed and that we discuss below,
we adopted the same input values as in Biasiotti et al. (2022).

Figure 1: Left panel: band-averaged altitude of modern Mars surface with
respect to the altimetric datum, calculated from the MOLA Mission Experiment
Gridded Data Records (MEGDRs) map. Blue line: original resolution of 0.125°
per band. Red line: ESTM resolution of 3° per band and modern ice caps
removed. Gray and black dashed lines: isolines corresponding to a 150-m GEL
ocean and a 550-m GEL ocean, respectively. Right panel: The average planetary
temperature at the altimetric datum for different values of the surface
pressure. Red, cyan, blue, black and green lines refer respectively to the
A0F0, B0G0, BLG0, BHG0 and C0G0 models described in Table 3. Each point
represents an increase of 33% of the pressure with respect to the preceding
one.

####  2.1.1 Topography

Mars shows a marked topographical dichotomy between the northern and southern
hemisphere. Approximately a third of the martian surface is composed of
relatively young lowland plains, while the rest are highly cratered highlands.
This geographical feature precedes the impact that formed the Hellas basin
(around 4.0 Gyr ago, Fassett & Head, 2011) and has been attributed either to
plate tectonics (Sleep, 1994), to mantle convection (Zhong & Zuber, 2001;
Šrámek & Zhong, 2012) or to a giant impact (Wilhelms & Squyres, 1984; Andrews-
Hanna et al., 2008; Marinova et al., 2008).

The martian dichotomy deeply influences the climate of Mars, especially in
presence of a thicker atmosphere. First, it impacts the global general
circulation, changing how heat is horizontally redistributed on the surface.
Second, it changes the actual surface temperature $T_{s}$ with respect to the
temperature at a given reference altitude (specifically, the altimetric datum,
see below) $T_{0}$ via the atmospheric lapse rate.

ESTM calculates the seasonal-latitudinal map of temperatures $T(t,x)$ at a
single altitude across the entire planet. In this paper, we refer to these
temperature values as $T_{0}$. However, since the martian dichotomy is mostly
a latitudinal feature, we were able to include its effects by adjusting the
ESTM output for the band-averaged lapse rate, thus deriving the actual surface
temperature, $T_{s}$. In practice, we followed the steps described below.
First, we took the altimetric map with a 0.125° horizontal resolution
generated by the Mars Orbiter Laser Altimeter (MOLA) experiment (Smith et al.,
2003) and calculated the band averaged altitude $h$ of the surface. We then
binned the resultant $h$ at a 3° resolution in order to use it in ESTM and
removed the ice caps, defined as the excess altitude above 77° and below -83°
in latitude. The band averaged altitude profiles that we calculated are shown
in Fig. 1, left panel. As in the original MOLA maps, the zero-altitude point
(also called altimetric datum) is the average Mars radius at the equator.
Then, we calculated the lapse rate for each seasonal and latitudinal point of
our runs. The lapse rate of the martian atmosphere changes with temperature
both in the dry case, where it can be ascribed to changes in the specific
heat, and in the moist case, where it is caused mainly by different amounts of
water vapor in the column. In particular, in this work we tested models with
relative humidities ($RH$) equal to 20% and 40% (see Section 2.2.1). The
equations we used to calculate the lapse rate are described in Appendix A. We
find the lapse rate to be limited between 4.2 and 5.2 K km-1 across all
models. This translates to adjustments of the surface temperature that are
$\leq$10 K across most of the martian surface, reaching $\sim$25 K in the most
extreme cases. Finally, we applied these temperature adjustments, obtaining
$T_{s}$ maps in place of the $T_{0}$ maps produced by ESTM. The results
presented in Section 3 and discussed in Section 4 are based on the post-
processed $T_{s}$ maps.

A third effect of the martian altimetry on the surface temperature is caused
by the amount of atmosphere above different latitudinal bands. This changes
the radiative exchanges between the surface and the outer space, decreasing
the OLR and increasing the TOA albedo locally. With respect to the pressure at
the altimetric datum $P_{0}$, the surface pressure $P_{s}$ at the highest
latitudinal band in our topographical model (the South Pole, at 2180 m) is 15%
lower, while its value at the lowest point (the North Pole, at -4740 m) is 40%
higher. While this radiative effect can not be directly included in ESTM, we
estimated its impact by considering how the average temperature of our models
is influenced by $P_{0}$, as done e.g. by Forget et al. (2013). The result is
shown in Fig. 1, right panel, for the five atmospheric compositions and cloud
models considered in this work (see Table 3). The spacing in pressure is
logarithmic and each pressure value represents an increase of 33% with respect
to the preceding one. As such, the temperature difference between a point in
the plot and its two nearest neighbours can be taken as an approximation of
the net radiative effect that is expected to occur at different locations on
the martian surface. This difference is largest in the high methane models
(BHG0 and BHGH), at the second highest tested pressure (4.22 bar) and amounts
to $\sim$13 K, while it is smallest in the dry, CO2-only model at around 2
bar, where it amounts to less than 0.2 K. At low ($<$0.6 bar) pressures, where
correctly accounting for the surface temperature is important to assess
whether the atmosphere collapses or not, the difference is of the order of
$\sim$2 K, which is four to ten times smaller than the effect of the lapse
rate. The actual difference would be smaller due to the redistribution of
heat, and this effect is stronger in denser atmospheres.

The columnar mass of atmosphere is influenced also by the oblateness of Mars
and the latitudinal variation of $g$ due to planetary rotation. While these
minor effects are not accounted for in ESTM, we can try to estimate them. The
tangential velocity at the martian equator is $\sim$240 m s-1, which decreases
the perceived value of $g$ by 0.017 m s-2. With respect to the average $g$ of
Mars, this amounts to a reduction of $\sim$0.45%. The difference between the
polar and equatorial radii (3376 km and 3396 km, respectively, Seidelmann et
al., 2007), causes $g$ to decrease by $\sim$0.4%, thus the combined effect
amounts to $\sim$0.9%. This variation in the columnar mass of the atmosphere
is expected to cause, at maximum, a temperature increase in the equatorial
regions of 0.35 K with respect to our calculations. In most of the tested
cases and at higher latitudes, however, the effect is much smaller.

####  2.1.2 Ocean fraction

A number of independent clues seems to suggest the existence of a water ocean
on the northern hemisphere of Mars at some point in the past, and most
probably during the Noachian and Hesperian periods ($\sim$4.1-3.0 Gyr ago).
Just a few examples of them are the similar height of ancient river outlets
deposits with respect to the datum (di Achille & Hynek, 2010), the sedimetary-
like dielectric properties of the northern terrain (Mouginot et al., 2012) and
the shape of certain northern craters suggestive of an impact in shallow
waters (Costard et al., 2019). On the other hand, the ability to support a
large, open ocean under the FYS conditions is central in the debate between
supporters of a warm-and-wet and a cold-and-icy Early Mars (see e.g.
Wordsworth, 2016).

The global equivalent layer (GEL) of water present on Early Mars is still a
matter of debate, with estimates ranging from 100 m to 1500 m (e.g. Scheller
et al., 2021). The GEL also underwent evolution over time due to sequestration
in the crust (Scheller et al., 2021) and/or photodissociation and atmospheric
escape (Webster et al., 2013; Villanueva et al., 2015), becoming smaller and
smaller during the Noachian and Hesperian ages until it reached the current
estimated value of 10-20 m (Plaut et al., 2007).

In this work we tested three different scenarios in terms of surface water:
(i) no oceans, (ii) oceans produced by a 150-m GEL and (iii) oceans produced
by a 550-m GEL. The first scenario considers the case of a planet where water
covers a negligible fraction of the surface. This includes the case in which
most of the water is locked in glaciers or subsurface lakes (Orosei et al.,
2018). The second scenario considers the long-term water level identified by
Carr & Head (2003), Ivanov et al. (2017) and Salese et al. (2019) and adopted
also in other studies (e.g. Schmidt et al., 2022a), and explores the lower end
of the estimated Early Mars water GEL. Since most of the evidences on which it
rests are linked to Hesperian features, this can be also related more
specifically to the conditions of the Hesperian Mars. The third and last
scenario adopts instead the estimates of di Achille & Hynek (2010) and is
somewhat in the middle of the Early Mars water GEL interval. This is more
relevant for the study of a late Noachian Mars.

ESTM allows us to specify the ocean fraction for each latitudinal band. This
is usually not useful when dealing with exoplanets (see e.g. Silva et al.,
2017, in which the ocean fraction is the same at each latitude). However, the
Noachian and Hesperian martian topography is reasonably constrained, as it is
its impact on the past distribution of water on the planet. Thus, we derived
the latitudinal distribution of the martian paleo-ocean from the same
topographical map employed for estimating the latitudinal altitude profile and
used it in our climate model. Among different paleo-shorelines identified in
literature, we adopted the one at the -3900 m height (Carr & Head, 2003) in
scenario (ii), and the one at the -2540 m height (di Achille & Hynek, 2010) in
scenario (iii). At this point we averaged over all longitudes to obtain the
latitudinal ocean fraction distribution. The curves obtained for scenarios
(ii) and (iii) are shown in Fig. 2, left panel, and correspond to the global
ocean coverage of 0.16 and 0.30, respectively.

####  2.1.3 Surface ices

Figure 2: Left panel: band averaged ocean fraction based on the position of
the martian paleo-shoreline in the models with a 150-m GEL (late Hesperian
ocean, as modeled by Schmidt et al., 2022a) or in the model with a 550-m GEL
(Noachian ocean, as in di Achille & Hynek, 2010). Right panel: CO2 cloud
effects on the average planetary temperature with respect to the clear-sky
case in the Forget et al. (2013, dashed line) and in the A0F0 model tested in
this work (see Tab 3).

The martian topography has an important impact also on the distribution of
surface ice. First of all, the effect of the lapse rate may prevent the
freezing of the low-lying water bodies even if the datum temperature is below
the freezing point of water. Second, the southern highlands and the Tharsis
bulge tend to act as cold traps, blocking most of the martian water reservoir
and decreasing the planetary albedo (Wordsworth, 2016).

ESTM lacks the ability to model the deposition and transportation of ice and
instead evaluates the fractional surface covered by ice at a given latitude as
a function of $T_{s}$ averaged over 12 months. In ESTM, the dependence of the
ice coverage on temperature is given by a generalized logistic function with
two different sets of parameter values to describe land and ocean ice coverage
(Eq. 4 in Biasiotti et al., 2022). These parameters have been tuned to
reproduce Earth data. In general, this approach is adequate to describe Earth-
like seasonal deposition of snow on lands and formation on ice on oceans at
high latitudes.

Since we decided to model Early Mars as an essentially dry planet, with a
desert-like average humidity, we self-consistently expect only sporadic
precipitations across the southern highlands. Thus, as far as lands are
concerned, we choose to turn off the temperature-dependent algorithm of ice
formation. To explore the impact of ice over land, we instead calculated the
seasonal temperature variations for a range of average land albedos
$\alpha_{l}$ (see Sect. 3.3) up to values consistent with a permanent ice
sheet of 50% ($\alpha_{l}=0.41$). On the other hand, as far as oceans are
concerned, we kept the temperature-dependent algorithm of ice formation. The
algorithm was fed by the surface temperature of the oceans calculated with the
procedure described in Sect. 2.1.1. With respect to the temperature of the
datum, we find that the surface ocean temperature is $16-20$ K higher in
scenario (ii), and $10-13$ K higher scenario (iii), depending on the specific
lapse rate obtained in a given run. Due to the band-averaged nature of our
calculations, the actual ocean surface temperature is adjusted within $\sim$2
K at most for models using 150 m of water GEL and $\sim$1 K at most for the
model using 550 m of water GEL.

As a minor note, we underline the fact that, since the original value of the
turning point temperature of the logistic for the ice coverage fraction has
been calibrated on Earth’s oceans, it best describes the behaviour of water
with the same average content of salt of Earth water, i.e. 3.5% (Chester &
Jickells, 2012). This corresponds to a freezing temperature of 271.3 K.

####  2.1.4 CO2 and H2O cloud modeling

Cloud modeling is generally the strongest source of uncertainty when working
with climate simulations in conditions that are different from those of the
modern Earth atmosphere (see e.g. Pierrehumbert, 2010). Clouds impact the
climate of a planet both by increasing the albedo, which cools down the
planet, and by reducing the OLR, which heats up the planet. These mechanisms
act together and the net effect depends on a large amount of factors that are
often difficult to determine, such as the height at which clouds typically
form, the distribution of particle sizes inside the cloud and the total amount
of condensate in the column (Ramirez & Kasting, 2017).

In ESTM clouds are parameterized via seven parameters linked to simple
physical or geographical properties, which are summarized in Tab. 2 and
described below. The effect of clouds on the OLR is modeled via a cloud OLR
forcing parameter (Cloud Radiative Effect, $CRE$), which can be fixed or made
dependent on the average planetary temperature. The effect of clouds on the
TOA albedo are modeled via three parameters: the cloud albedo,
$a_{c}(\mu=0.5)$, the slope of the linear dependence between $\mu$ and the
cloud albedo, $m_{c}$, and the cloud shortwave transmittance, $t$. The overall
impact of clouds then scales linearly with the fractional cloud cover over
each latitudinal band. In ESTM we specify different cloud coverage fractions
for each of the three different planetary surfaces considered in the model:
$f_{cw}$ for water, $f_{cl}$ for land and $f_{ci}$ for ice. Thus, the final
cloud coverage at a given latitude is the weighted average between the three
values. The functional dependence of the parameters described in this
paragraph on the main radiative terms considered in ESTM are detailed in
Section 2.5 of Biasiotti et al. (2022).

Input | CO2 ice | H2O Earth-like | ESTM
---|---|---|---
param. | clouds | clouds | calibration
$CRE$ | 15.0 | 26.1 | 5.2-26.1
$a_{c}(0.5)$ | 0.28 | 0.44 | 0.44
$m_{c}$ | -0.67 | -0.67 | -0.67
$t$ | 0.95 | 0.85 | 0.85-0.95
$f_{cw}$ | 0.5 | 0.5 | 0.72
$f_{cl}$ | 0.5 | 0.5 | 0.55
$f_{ci}$ | 0.5 | 0.5 | 0.56
Table 2: Parameters adopted to model CO2-ice clouds (second column), Earth-
like clouds (third column) and adopted in the ESTM reference paper Biasiotti
et al. (2022) (fourth column). Values for $CRE$ and $t$ in the fourth column
are expressed as an interval since they are temperature-dependent.

In this work, we tested two different models for the Early Mars clouds. The
adopted values for the parameters are reported in Table 2. In order to relate
more directly the input choices and results, we deactivated the dependence of
$CRE$ and $t$ on $T_{s}$. The first model tries to capture the impact of
CO2-ice clouds. In practice, we tuned the parameters in ESTM as to obtain an
overall average $T_{s}$ difference with respect to the clear-sky scenario that
is comparable with the results of Forget et al. (2013). In that paper, the
thermal enhancement caused by clouds is explored for a range of surface
pressure comprised between 0.1 and 6.5 bar, however the authors report that
simulations with pressures above 3 bar are out-of-equilibrium because of
atmospheric collapse. The comparison between our results and those of Forget
and collaborators is shown in Fig. 2, right panel.

The second model is instead similar to that adopted in Biasiotti et al. (2022)
to describe the modern Earth clouds. The values chosen for $CRE$ and $t$ are
adequate for clouds that are present over non-icy surfaces (i.e., outside the
Poles). Biasiotti et al. (2022) report that the exact choice of $t$ has only
minor impacts on the final result.

At variance with Biasiotti et al. (2022) we fixed the fraction of clouds over
all surfaces $f_{c}$ to the same value, in both models with CO2-ice clouds and
H2O Earth-like clouds. This value is set to 0.5 when not otherwise specified.
This choice is driven by the fact that the cloud fraction values adopted in
Biasiotti et al. (2022) are meaningful only for the specific geographical
configuration of the present time Earth and cannot be transposed easily to the
conditions of the Early Mars. Using the same fractional cloud coverage across
all models allows us to compare the results of different ESTM runs more
straightforwardly. In Sect. 2.1.4 we explored also the effect of different
cloud coverages in the [0, 1] range.

###  2.2 EOS

The OLR and TOA albedo lookup tables used in ESTM have been provided by the
EOS code (Simonetti et al., 2022). EOS is based on the RT code HELIOS (Malik
et al., 2017, 2019) and the opacity calculator HELIOS-K (Grimm & Heng, 2015;
Grimm et al., 2021), but tailored for the study of terrestrial-type
atmospheres. EOS inherits the feature of being GPU-accelerated. It operates in
the two-stream approximation with an improved non-isotropic scattering
treatment (Heng et al., 2018) and can carry on calculations either line-by-
line or using k-distribution tables.

For this work we adopted the spectral line lists for CO2, H2O and CH4
available on HITRAN2020 (Gordon et al., 2022). Calculations have been
performed using k-distribution tables that we checked against previous EOS
results obtained with the line-by-line procedure. The absorption of each gas
has been calculated on a log-spaced grid in pressure with a 1/3 dex step in
the interval between $10^{-6}$ and $10^{3}$ bar, and a linearly spaced grid in
temperature with a 50 K step between 50 and 450 K. The spectral interval
considered is 0-35000 cm-1. The CO2 absorption lines have been modeled
following Perrin & Hartmann (1989) and cutting the wings at 500 cm-1 from the
line centre, while H2O and CH4 lines have been modeled using a Voigt
distribution cut at 25 cm-1 from the line centre. For H2O, the residual
absorption at distances larger than 25 cm-1 (the so-called Lorentzian
pedestal) has been removed, since it is already accounted for in the
collisional continuum. We then included the CO2-CO2 Collision-Induced
Absorption (CIA) in the 0-750 cm-1 (Gruszka & Borysow, 1997, 1998) and
1000-1800 cm-1 (Baranov & Vigasin, 1999) ranges (the so-called GBB continuum
used in Wordsworth et al., 2010) and the CO2-CH4 CIAs in the 0-2000 cm-1 range
(Wordsworth et al., 2017). Concerning H2O, we used the self- and foreign-
induced continua calculated using the MT_CKD v3.4 (Clough et al., 1989; Mlawer
et al., 2012) model available as a standalone subroutine of the LBLRTM code
(Clough et al., 2005). It should be noted that the foreign-induced H2O
continuum is calculated for a modern Earth-like composition and thus it might
deviate from that produced by a CO2-dominated atmosphere. However, previous
versions of this model have been used in other studies of the Early Mars
atmosphere (e.g. Wordsworth et al., 2013), thus at least guaranteeing some
level of intercomparison. The atmosphere has been modeled as composed by 10
log-spaced layers per order of magnitude in pressure. Since adopting a coarser
or finer pressure grid for the atmosphere impacts the EOS results (Simonetti
et al., 2022, Fig. 2), we chose the total number of layers in order to be as
homogeneous as possible across different cases. The pressure at the upper face
of the uppermost layer ($P_{\text{TOA}}$) is always $10^{-5}$ bar. Rayleigh
scattering cross-sections for CO2 and CH4 (Sneep & Ubachs, 2005) and H2O
(Wagner & Kretzschmar, 2008) are also included. RT calculations are all clear-
sky, since the effect of clouds are already included in ESTM (see Section
2.1.4).

The adopted vertical atmospheric pressure-temperature (PT) profile consists of
a convective troposphere and an isothermal stratosphere. This vertical PT
profile considers the condensation of H2O (if any) in the lower troposphere,
the condensation of CO2 (if possible) in the upper troposphere and an
isothermal stratosphere at 155 K. The non-ideal behaviour of the involved
gases is taken into account (see Appendix A). The resultant lapse rates are
used both for RT calculations and the altimetric adjustment described in
Section 2.1.1.

For each of the tested atmospheric composition and $P_{0}$ in the [0.1, 5.62]
bar range we constructed two lookup tables, which are then used in the ESTM
climate model: one containing the OLR as a function of $T_{0}$, and another
containing the TOA albedo as a function of $T_{0}$, the surface albedo
$\alpha_{s}$ and the stellar zenith angle $z$. These have been obtained by
running EOS for each combination of the relevant parameters. The temperature
grid is linearly spaced, with a resolution of 10 K in the 160-240 K range and
a resolution of 20 K in the 240-360 K range, the zenith angle grid is
$\\{0^{\circ},45^{\circ},\,60^{\circ},\,70^{\circ},\,80^{\circ},\,85^{\circ},\,88^{\circ}\\}$
and the surface albedo grid is $\\{0,\,0.30,\,0.45,\,0.60,\,0.90\\}$. Some of
the $(T_{0},P_{0})$ points are associated with unphysical conditions, i.e.
positions in the pressure-temperature plane that resides below the
condensation curve of CO2. In that case, we substituted the surface pressure
value taken from the pressure grid with the CO2 condensation pressure at that
temperature. As an example, the OLR and TOA albedo in the
$(200\,\text{K},2.0\,\text{bar})$ point has actually been calculated using
1.578 bar. This choice allows us to self consistently take into consideration
the radiative impact of atmospheric condensation at surface. The OLR for
$T_{0}\leq 150$ K is calculated from the Stefan-Boltzmann law, and is linearly
interpolated between 150 K and 160 K. This is reasonable considering that the
CO2 condensation pressure at 150 K is 0.008 bar and thus the expected
atmospheric radiative forcing is small. For $T_{0}<160$ K the TOA albedo
depends on the zenith angle and the surface albedo, but not on $T_{0}$, since
below that temperature a CO2-dominated atmosphere must be thin due to
condensation.

####  2.2.1 Atmospheric chemical composition

The atmospheres tested in this work are all CO2-dominated. This is a standard
choice for the study of Early Mars (e.g. Kasting, 1991) and it is justified by
considerations on the chemical composition of the outgassed fluids (Grott et
al., 2011). At present time, the martian atmosphere includes a 5% by volume of
N2 and some modelists (e.g. Ramirez et al., 2014) scaled this abundance in
their models. Anyway, our choice to not include N2 has only negligible effects
on the radiative properties of the atmosphere, as it is shown in Ramirez
(2017). Fig. 1, right panel, confirms this result: substituting 5% of CO2 with
N2 in our model would reduce $\overline{T}$ by $\sim$2 K at maximum, and by a
much smaller quantity in most other cases.

We also included a variable amount of H2O, by fixing the relative humidity
$RH$ to either 0%, 20% or 40%. Since our RT lookup tables must be
representative of the average atmospheric mixture, we chose low $RH$ values
that are representative of the typical humidities found in modern Earth
deserts (Yang et al., 2020). More specifically, we associated the null $RH$
value to the model with no oceans, the intermediate $RH$ value to the models
with the 150-m GEL ocean and the highest $RH$ value to the model with the
550-m GEL ocean.

We then tested three different fractional amounts of CH4, $f_{\text{CH4}}$: 0,
0.1% and 1% by volume. There are two reasons behind these specific choices.
First, low amounts of CH4 prevent the formation of organic hazes via
photochemical reactions in the upper atmosphere. These hazes form efficiently
when the CH4/CO2 ratio increases above $\sim 0.1$ (Trainer et al., 2006;
DeWitt et al., 2009). Second, low amounts of CH4 are possibly easier to
explain in light of the known geochemical features of Mars. In general, the
amount of CH4 in the Early Mars atmosphere is understood to be primarily
controlled by the balance between serpentinization111The process in which
mafic and ultramafic rocks are oxidized at relatively low temperatures, i.e.
in the planetary crust. This process sequester water and releases reduced
species, such as CH4. of the olivine crust (Chassefière et al., 2013) and
photo-dissociation (Wordsworth et al., 2017) in the atmosphere. The
photodissociation rate is limited by the availability of photons with
$\lambda<160$ nm, most of which are Lyman-$\alpha$ photons, whose flux at the
Noachian-Hesperian boundary was about $3\times 10^{15}$ m-2s-1 (Wordsworth et
al., 2017). In few bars-level CO2 dominated atmospheres, this limit is reached
for $f_{\text{CH4}}\sim 1\%$ (Zahnle, 1986). If we approximate the
photodissociation rate as linearly dependent on concentration when it is below
this limit, we can make a zeroth-order estimate of the flux needed to maintain
a given atmospheric concentration. On the basis of the photon flux, the source
flux required to support a 1% CH4 concentration would be $3\times 10^{15}$
m-2s-1, or one order of magnitude smaller for a 0.1% concentration. This
quantity can be compared with the flux from regions undergoing
serpentinization, that has been estimated to be in the $10^{16}-10^{18}$
m-2s-1 interval (Etiope et al., 2013). Taking the lower end of this interval,
the methanogenic regions must cover 30% of the martian surface to maintain our
highest tested concentration. While this estimate is extremely simplified, it
allows to check the plausibility of our choices.

##  3 Results

Name | RH | $f_{\text{CH4}}$ | Cloud | GEL
---|---|---|---|---
| [%] | [%] | model | [m]
A0F0 | 0 | 0 | CO2-ice | 0
B0G0 | 20 | 0 | Earth-like | 0
B0GH | 20 | 0 | Earth-like | 150
BLG0 | 20 | 0.1 | Earth-like | 0
BLGH | 20 | 0.1 | Earth-like | 150
BHG0 | 20 | 1 | Earth-like | 0
BHGH | 20 | 1 | Earth-like | 150
C0G0 | 40 | 0 | Earth-like | 0
C0GN | 40 | 0 | Earth-like | 550
Table 3: A summary of the models studied in this paper. The first letter
refers to the relative humidity (”A” for 0%, ”B” for 20%, ”C” for 40%). The
second letter refers to the fractional amount of CH4 (”0” for 0%, ”L” -low-
for 0.1%, ”H” -high- for 1%). The third letter refers to the cloud model (”F”
for CO2-ice, ”G” for Earth-like). The fourth and last letter refers to the
surface water distribution and depth (”0” for no surface water, ”H” for an
Hesperian-like 150-m GEL, ”N” for a Noachian-like 550-m GEL).

Our goal is to explore the early martian climate over a wide range of
planetary and atmospheric parameters. In particular, we are interested in
studying which set of inputs enables localized, seasonal thaws under Noachian
and Hesperian Mars conditions.

The models tested in this paper are summarized in Table 3. For each of them we
varied (i) the datum pressure $P_{0}$ between 0.1 and 5.62 bar in
logarithmically spaced steps of 1/8 of dex and (ii) another parameter chosen
among the obliquity $\varepsilon$, the eccentricity $e$, the albedo of lands
$\alpha_{l}$ and the cloud fraction $f_{c}$. We also explored the interplay
between obliquity, eccentricity and the argument of the perihelion $\omega$,
for a fixed surface pressure value. We refer to a specific combination of
input variables for a given model as a “case” and the total number of cases
run in this study is $\sim$104.

In practice, we followed the steps described below. First, we run ESTM until
convergence is found for a given set of input parameters. The resulting
seasonal-latitudinal temperature map, $T_{0}$, refers to the altimetric datum
and not to the actual surface. Second, we correct the seasonal-latitudinal
$T_{0}$ map of ESTM for the altitude as described in Section 2.1.1, obtaining
$T_{s}$. Third, we search for the highest ($T_{\text{max}}$) and lowest
($T_{\text{min}}$) temperature in the $T_{s}$ map. If $T_{\text{max}}$ is
above the freezing temperature of water, then we calculate the fraction of the
martian year for which this condition is satisfied at that specific latitude.
$T_{\text{max}}$ imposes the most stringent limit on the melting of ices at a
given latitude. Fourth, we check which latitude lasts the most above-freezing
temperatures and save the result as a fraction of martian year. Usually, but
not always, it is the same latitude at which $T_{\text{max}}$ occurs. In
particular, we considered two different threshold temperatures
$T_{\text{thr}}$ as ”above-freezing”: 273 K and 252 K. The first one refers to
pure water, the second one refers to a saturated solution of NaCl (i.e. brine,
Fairén, 2010). These thresholds are reported in Figs. 3-7 as red and yellow
lines, respectively. Finally, using $T_{\text{min}}$, we check the run for CO2
condensation at surface, which signals the potential beginning of an
atmospheric collapse. Notice that we do not compute the duration or magnitude
of this condition: even if condensation occurs at a single latitude point and
for a single time-step interval ($\sim$14 d), that run will still be flagged
as ”condensing”. In Figs. 3-7, the regions of the parameter space in which we
found condensation of the atmospheric CO2 at surface are hatched. Notice that
$T_{\text{max}}$ and $T_{\text{min}}$ generally occurs at different latitudes
and positions along the orbit.

All runs start using an isothermal latitudinal profile at 300 K. Thus, no
ocean ice is present at the initial time of the simulation. Due to climate
bistability, the final climate state might depend on the starting temperature.
However, since: (i) the only source of bistability in ESTM is the ice-albedo
feedback (Murante et al., 2020), (ii) only ocean ices are tracked dynamically
and (iii) the ocean coverage is small, we expect bistability to be generally
negligible.

###  3.1 Dependence on obliquity

Figure 3: The maximum surface temperature $T_{\text{max}}$ as a function of
the surface pressure and obliquity. From top left moving clockwise: maps for
the B0GH, BLGH, BHGH, C0GN models. Yellow and red contour lines highlight the
regions of the parameter space for which pure liquid water and brine can be
liquid, respectively. Hatched regions refer to combinations of parameters that
causes the CO2 in the atmosphere to condense at surface.

We explored the range of obliquity values between 0∘ and 60∘, sampled at 5∘
intervals. These boundaries have been chosen considering the results of Laskar
et al. (2004), whose calculations showed that the probability of an obliquity
excursion up to 60∘ in a span of 250 Myr is around 45%. The results for the
models B0GH, BLGH, BHGH and C0GN are shown in Fig. 3.

All models exhibit a similar structure, with a low pressure ($<1.0$ bar), mid-
to-high obliquity ($>10^{\circ}$) region of the parameter space in which
atmospheric condensation at surface occurs at some moment along the orbit. The
limiting pressure that identifies this region is particularly easy to
identify, since no case with 1.78 bar or more in models B0GH and C0GN and no
case with 1.33 bar or more in models BLGH and BHGH is predicted to have
surface condensation, independently of the obliquity.

In general, $T_{\text{max}}$ is roughly independent of pressure when
$P_{0}<1.0$ bar, and roughly independent of obliquity when $P_{0}>4.0$ bar.
Between these two values there is a transition region in which both parameters
influence the outcome. The low pressure region shows some non-linear
behaviour, since all models exhibit a local maximum of $T_{\text{max}}$ for
$P\sim 0.4$ bar and $\varepsilon=60^{\circ}$.

Both the transition from obliquity-dominated to pressure-dominated
$T_{\text{max}}$ and the low-pressure local maximum are caused by the
interplay between the increasing heat redistribution efficiency and the
increasing greenhouse effect of a thicker atmosphere. Between 0.1 and 0.4 bar,
$T_{\text{max}}$ likely increases because the stronger infrared absorption
overcomes the more efficient redistribution of heat, that tends to smooth out
the seasonal and latitudinal differences. Above 0.4 bar, instead, the
diffusion term becomes more important, leading to lower $T_{\text{max}}$.
Finally, above 1.33 bar, $T_{\text{max}}$ starts increasing again because the
entire planet becomes warmer, while latitudinal variations are suppressed even
at high obliquities.

In model C0GN, there is a rapid transition in the $T_{\text{max}}$ values near
the isotherm at 273 K. This transition is present also in the other models but
it is less evident, and nearly disappears in the BHGH case. This behaviour is
caused by the transition from a completely and continuously glaciated northern
ocean, and a mostly open ocean. Oceans have a very low albedo even compared to
martian regolith ($\sim$0.05222For $z\sim 60^{\circ}$. vs 0.215, see e.g.
Huang et al., 2019, for data on the former), thus they contribute to warm up
the planet. Since the ocean fraction is highest in the C0GN model, their
effect is strongest.

We do not show the results for the A0G0 model because nearly all the tested
obliquity-pressure combinations caused the atmosphere to collapse at the
surface. More specifically, only the cases with $P_{0}=$ 0.13 or 0.18 bar and
$\varepsilon=5^{\circ}$ are stable, and their $T_{\text{max}}$ is equal to
236.9 and 237.9 K respectively. In the A0F0 model the maximum greenhouse limit
is reached at around 2 bar (see Fig. 1, red line), which causes the
disappearance of the conditions conductive to seasonal thaws at higher
pressures.

###  3.2 Dependence on eccentricity

Figure 4: The maximum surface temperature $T_{\text{max}}$ as a function of
the surface pressure and eccentricity. From top left moving clockwise: maps
for the B0GH, BLGH, BHGH, C0GN models.

The second parameter we investigated is the orbital eccentricity, for which we
explored the range of values between 0 and 0.18, in 0.015 steps. Again, the
boundaries were informed by the work of Laskar et al. (2004), which predicted
a potential eccentricity excursion up to $\sim 0.2$ in their 5 Gyr runs. At
variance with the obliquity interval studied in the previous Section,
eccentricity values above 0.14 have a low ($<5$%) occurrence probability in a
250 Myr time span and are thus unrealistic, even though we included them for
completeness.

Eccentricity variations have two effects on the surface temperature of the
planet. First, increasing eccentricity increases the average instellation of
the planet by a factor $(1-e^{2})^{-1/2}$ (Williams & Kasting, 1997). Second,
the higher eccentricities contribute to seasonal variations, with the effect
being strongest when $\omega=90^{\circ}$ or $270^{\circ}$ and weakest when
$\omega=0^{\circ}$ or $180^{\circ}$333If $\omega=90^{\circ}$, the planet is at
periastron during the southern hemisphere summer solstice, meaning that the
seasonal variations due to obliquity and eccentricity compound together. The
same happens if $\omega=270^{\circ}$, but considering the northern hemisphere
summer solstice. Instead, when $\omega=0^{\circ}$ or $180^{\circ}$, the planet
is at equinox when it passes through the periastron.. In the models we are
considering here, the first effect is small since going from 0 to 0.18
increases the martian insolation by 1.7%, while the second is not sufficient
to influence the final results, despite that $\omega=286.5^{\circ}$. This is
evident observing the results in Fig. 4: most of the variation happens on the
vertical axis, i.e. when we changed the surface pressure. At low pressures
($<1$ bar) there is a modest increase of $T_{\text{max}}$ for higher
eccentricities, but the limiting pressure required for thawing, under the more
conservative condition, changes only by $\sim 30\%$.

Seasonal variability due to obliquity is dominant also because it sets the
limiting $P_{0}$ to avoid atmospheric collapse at the surface. $P_{0}$ is
equal to 1.78 bar for models B0GH and C0GN, 1.33 bar for model BLGH and 1 bar
for BHGH. In models B0GH and C0GN this limit is reduced to 1.33 bar for the
highest eccentricity cases, namely those above 0.15.

Once again, we do not show the results for the A0F0 models, since for all the
tested combinations of pressure and eccentricity the atmosphere is unstable
against surface condensation. This is expected, since seasonality is dominated
by obliquity and for $\varepsilon=25^{\circ}$ there are no non-condensing
cases (see the previous Section). The highest $T_{\text{max}}$, just a few
kelvins above the freezing temperature of pure water, are found for $e>0.15$
and $0.13<P_{0}<0.421$, and as for the other models, most of the variations
are caused by pressure changes.

###  3.3 Dependence on land albedo

Figure 5: The maximum surface temperature $T_{\text{max}}$ as a function of
the surface pressure and the land albedo. From top left moving clockwise: maps
for the B0GH, BLGH, BHGH, C0GN models.

The third parameter tested is the albedo of lands. Currently, its average
value is 0.215, but during the Noachian and Hesperian its value could have
been different. Two factors that can change the land albedo are: (i) large
volcanic eruptions that cover a substantial fraction of the surface with fresh
basalt, whose albedo is lower than that of regolith, and (ii) reflective
permanent glaciers in the southern highlands. In this work, we varied the land
albedo between 0.11 and 0.41 in 0.015 steps. The lower end is meant to
represent the albedo of surface mostly covered in basalt quenched glass (whose
albedo is 0.09, Essack et al., 2020), while the upper end is meant to
represent a surface which is 50% covered in modern Mars regolith and 50%
covered in ice (albedo: 0.6), or an Earth-like dry sand.

The results of this exploration are shown in Fig. 5. The effect of changing
the land albedo is stronger in models B0GH and C0GH and weaker in models BLGH
and BHGH, and decreases when $P_{0}$ increases. In models B0GH and C0GN, high
land albedos prevent seasonal deglaciations both using the more conservative
and the more lenient condition, even at the highest pressure tested. Samely,
they favor the collapse of the atmosphere when higher than 0.32 (for model
B0GH) and 0.38 (for model C0GN). On the other hand, low albedo runs still
require at least 1 bar of CO2 both in B0GH and C0GN in order to prevent
atmospheric collapse. As in the other set of calculations, the C0GN model
shows a faster transition in $T_{\text{max}}$ around the 273 K isoline, which
again underlines the transition from a completely frozen northern ocean to a
mostly open one, and it is more evident than in the other models because this
ocean covers a larger fraction of the surface. This transition is steeper for
high land albedos due to the higher albedo contrast between light lands and
dark oceans.

The models BLGH and BHGH allows instead a stable atmosphere and seasonal thaws
in the northern lowlands even for the highest land albedo value tested,
provided the surface pressure is sufficiently high. In model BLGH the limiting
$P_{0}$ is 3.16 bar and in model BHGH it is 1.33 bar. On the other hand, the
A0F0 model (not shown) can neither avoid surface atmospheric condensation nor
produce seasonal thaws (under the pure water condition), even at the lowest
land albedo tested.

###  3.4 Dependence on the cloud fraction

Figure 6: The maximum surface temperature $T_{\text{max}}$ as a function of
the surface pressure and the cloud fraction. From top left moving clockwise:
maps for the B0GH, BLGH, BHGH, A0F0, C0GN models. We recall that in models
B0GH, BLGH, BHGH and C0GN we use Earth-like H2O clouds, while in model A0F0 we
use CO2 ice clouds.

The fourth and final parameter considered is the cloud cover fraction.
Modeling the cloud coverage and its effects in conditions different from that
of the modern Earth is a notoriously difficult task, and different GCMs can
produce substantially different results for the same input (see e.g. Sergeev
et al., 2022). In absence of a self-consistent model for clouds, we wanted to
explore the implications of the simple parameterization presented in Sect.
2.1.4 by varying the cloud cover fraction between 0% and 100%, in 10% steps.
Alternatively, this analysis can be intended as studying how, for a fixed but
unknown coverage, the climate responds when our simplified cloud properties
are linearly scaled.

The results are shown in Fig. 6. Since we tested two different cloud parameter
sets (see Tab. 2), it is more meaningful to treat the models B0GH, BLGH, BHGH
and C0GN separately from the model A0F0.

####  3.4.1 H2O clouds

Concerning the cases for which we adopted Earth-like H2O clouds, any increase
in cloud fraction causes $T_{\text{max}}$ to increase. In other words, this
type of clouds are net greenhouse contributors across all the tested pressure
range. A direct consequence of this fact is that, in general, we need at least
some cloud coverage to avoid the atmospheric collapse and the presence of
seasonal thaws. The minimum coverage value for a stable atmosphere is between
40% and 50% for B0GH and C0GN and between 10% and 20% for BLGH. The BHGH model
always allows for both non-condensing and thawing conditions, irrespective of
the cloud coverage. All models see their atmospheric stability region widened
when the cloud fraction value is higher. At 100% coverage, all the models with
Earth-like H2O clouds are stable down to 0.32 bar (and BHGH down to 0.13 bar).

As can be seen from Fig. 6, the net greenhouse effect produced by clouds is
also strongly dependent on $P_{0}$, which in turn causes some interesting
secondary effects. First of all, at variance with all the previously tested
parameters, high pressure ($>4$ bar) cases have not always a hotter summer
than lower pressure ones. Apart for the BHGH model, the minimum
$T_{\text{max}}$ value is found at the maximum tested pressure, for a 0%
coverage. In the BHGH model, the minimum $T_{\text{max}}$ is encountered at
0.1 bar (244 K) as in the previous explorations. However, $T_{\text{max}}$
does not monotonically increase with pressure, and after reaching a maximum
value equal to 272 K it decreases again to 247 K at $P_{0}=5.62$ bar. Our
results indicate that at high pressures, the net effects of clouds is the
strongest, contributing up to $\sim$150 K to $T_{\text{max}}$ (in BHGH), while
for the lower pressures, the effect is as small as 7 K (in B0GH). A second
interesting effect is that, in all models other than BHGH, there exist a cloud
fraction value that makes $T_{\text{max}}$ insensitive to $P_{0}$. This value
is around 40% for B0GH and C0GN and around 20% for BLGH. This means that
variability on the vertical axis of all the models shown in this work depends
at least partially on the exact choice of the cloud parameters, and that the
stronger variability shown by the BLGH and BHGH models in Figs. 3, 4 and 5 is
at least partially linked to the specific cloud coverage that we chose.
Decreasing the cloud coverage, i.e. approaching the clear-sky case, reduces
the variability as a function of $P_{0}$ for all models, while the opposite is
true when we increase it.

####  3.4.2 CO2 clouds

The behavior of $T_{\text{max}}$ as a function of cloud coverage in the model
with CO2 ice clouds is somewhat different. First, this type of clouds have a
minor cooling effect at low pressure. As can be seen in the lower center panel
of Fig. 6 (model A0F0), going from 0% to 100% coverage, $T_{\text{max}}$
decreases from 267 K to 262 K at 0.1 bar and are substantially neutral at 0.32
bar. At 5.62 bar, they instead contribute $\sim 50$ K at 100% coverage. The
effect of clouds on $T_{\text{max}}$ is thus different from the one on the
average planetary temperature reported in Fig. 2, right panel. For pressures
above 1.33 bar, the greenhouse contribution of clouds is sufficient to prevent
the atmospheric collapse, provided the coverage is above 70%. However, even at
100%, $T_{\text{max}}$ is not sufficient for allowing surface thaws. The
maximum value of $T_{\text{max}}$ is 268 K for 0.24 bar and 0% cloud coverage,
while a second, local maximum equal to 265 K can be found at 5.62 bar and 100%
coverage.

At variance with the models using H2O Earth-like clouds, in A0F0 there is only
a minor increase in the sensitivity to $P_{0}$ variations when the cloud
fraction changes, and the maximum variability is produced when coverage is 0%,
rather than 100%. There is no coverage value that makes the model insensitive
to pressure changes as in B0GH, BLGH and C0GN and in general, the impact of
cloud fraction is smaller at any given pressure value.

###  3.5 Combined effects of the obliquity, eccentricity and argument of
perihelion

Figure 7: The maximum surface temperature $T_{\text{max}}$ as a function of
the eccentricity and the obliquity, for a fixed surface pressure of 1.33 bar.
Left column: model B0GH. Right column: model C0GN. Top row:
$\omega=90^{\circ}$ (perihelion at southern hemisphere’s summer solstice).
Middle row: $\omega=180^{\circ}$ (perihelion at equinox). Bottom row:
$\omega=270^{\circ}$ (perihelion at northern hemisphere’s summer solstice).

We can now study the combined effect of eccentricity and obliquity, for a
fixed $P_{0}$. We chose the models B0GH and C0GN for this analysis and we
varied $e$ and $\varepsilon$ in the same interval and with the same spacing
adopted in Sections 3.1 and 3.2. The surface pressure has been fixed at 1.33
bar, because it is the highest $P_{0}$ value for which atmospheric
condensation can either occur or not depending on the obliquity. In other
words, it is the $P_{0}$ value for which the output is maximally dependent on
the specific choice of $e$ and $\varepsilon$, and thus the most interesting to
investigate. Each obliquity-eccentricity combination has been run three times
for different values of the argument of perihelion, namely $180^{\circ}$,
$270^{\circ}$ and $90^{\circ}$.

The results are reported in Fig. 7. First of all, this test confirms that
obliquity is the main driver of seasonality, since $T_{\text{max}}$ variations
are more pronounced on the horizontal, rather than on the vertical axis,
except for very low ($<10^{\circ}$) obliquities. In general, the maximum
$T_{\text{max}}$ value is found for the highest tested $e$ and $\varepsilon$
values, even though its value is different when $\omega$ changes. In
particular, for model B0GH, the $T_{\text{max}}$ maxima are 276, 292 and 283 K
for $\omega=180^{\circ}$, $270^{\circ}$ and $90^{\circ}$, respectively. For
model C0GN they are 279, 288 and 286 K. For a limiting scenario as the one
studied here, these differences are sufficient to completely prevent seasonal
thaws.

Similarly, varying $\omega$ changes the regions of the $e-\varepsilon$ plane
that allows for a non-condensing atmosphere. Condensation usually occurs at
the South Pole due to the latitudinal altimetric profile of Mars. Thus, the
stronger southern seasonal temperature excursions produced when
$\omega=90^{\circ}$ make it easier to cross the condensation curve of CO2
during the colder South Pole winters. The opposite is true when
$\omega=270^{\circ}$: it is the northern hemisphere that is exposed to the
strongest seasonality, but due to the combined effect of altimetry, the
atmosphere is stable under a wider range of eccentricity-obliquity
combinations. The shape of the condensation region underlines again the role
$\omega$: in the $180^{\circ}$ and $270^{\circ}$ cases, increasing $e$ shrinks
it, while in the $90^{\circ}$ case widens it, preventing stability at lower
$\varepsilon$.

The structure of these results is similar also for the models BLGH and BHGH
(not shown), but the threshold $P_{0}$ for a non-trivial structure of the
condensation and thawing regions in the $e-\varepsilon$ plane is different.
For BLGH, this value is 1.0 bar, while for BHGH is 0.75 bar and since these
pressures are lower, the changes of $T_{\text{max}}$ are more pronounced. On
the other hand, for $P_{0}=1.33$ bar, most combinations allow for stable
atmospheres (again, $\omega=90^{\circ}$ produces a higher number of runs with
CO2 surface condensation than $\omega=270^{\circ}$) and roughly $\sim 25-35$%
of them allow for seasonal thaws.

###  3.6 Duration of the thaw seasons

Figure 8: The fraction of $P_{\text{orb}}$ during which the surface
temperature remains above 273 K (left column) or above 252 K (right column) in
at least one latitudinal band. All plots refer to the model BLGH. Top row:
dependence on pressure and obliquity. Mid row: dependence on pressure and
eccentricity. Bottom row: dependence on eccentricity and obliquity for
$\omega=270^{\circ}$. Hatched regions refer to combinations of parameters that
causes the CO2 in the atmosphere to condense at surface.

Finally, we analyze the duration of the thaw phases. The existence, in the
latitude-orbital phase plane for a given ESTM run, of point where $T_{s}$ is
above a given threshold (e.g. 273 K) does not ensure that an actual thaw can
happen. The duration of favorable conditions for deglaciation must be taken
into consideration. These conditions last differently at different latitudes.
In order to simplify the problem, we consider only the latitude at which
$T_{\text{max}}$ occurs, since in our EBM it usually corresponds to the
latitude for which the above-freezing season last the longest, and we express
it in terms of fraction of martian year.

The results for the BLGH model are reported in Fig. 8. The other models are
not shown but the general considerations that we draw here apply also on them.
A first distinction can be made between the region of each parameter space in
which deglaciated conditions are permanent from the region in which they are
seasonal. Permanently warm latitudinal bands are possible only above a certain
pressure threshold, independently of the other parameter varied, and are
associate with temperate global conditions (i.e. the ”warm-and-wet” Early Mars
scenario). Wherever the existence of a thaw season is primarily linked to the
$P_{0}$ choice (e.g. in the pressure-eccentricity plane), then the seasonality
region is small: either the parameter set always prevent deglaciation, or they
always allow it. This is not strange considering that, if thaws are driven by
the presence of relatively high pressures, then the more efficient heat
distribution limits the local seasonal variations of temperature.

On the other hand, seasonal thaws are possible at all pressures, even though
in some cases are prevented by atmospheric collapse (see e.g. the upper left
panel in Fig. 8). The duration of the warm season can in principle be very
short, but if it is present, it is usually longer than $\sim 0.15$ of martian
year ($\gtrsim 100$ d). In other words, when thawing can occur, then it is not
limited to negligibly short ($\lesssim 0.1$) fractions444Since we calculate
the $T_{s}$ for 48 points along the martian orbit, the minimum possible
fraction would be $\sim 0.02$.. For example, a typical fraction value along
the border between the non-deglaciating and the deglaciating regions for a 273
K threshold temperature (left column of Fig. 8) is 0.3. This happens because
$T_{\text{max}}$ nearly always occur in the northern hemisphere, where the
surface is covered in water: when ice melts the albedo decreases, making it
harder for ice to form again. Thus, this is basically a manifestation of a
bistability mechanism, albeit on a regional (rather than a global) scale.

##  4 Discussion

Figure 9: The latitude interval in which the temperature at the datum is
higher than $T_{\text{thr}}$ for at least one time step, as a function of the
datum pressure $P_{0}$. Results refer to the martian southern hemisphere. Top
row: $T_{\text{thr}}$ = 252 K. Bottom row: $T_{\text{thr}}$ = 273 K. From left
to right: B0GH, BLGH, BHGH and C0GN models. Red, green and blue shaded areas
refer to cases with $\varepsilon$ = $30^{\circ}$, $45^{\circ}$ and
$60^{\circ}$, respectively. The vertical dashed lines mark the pressure below
which surface atmospheric condensation occurs for each of the shown
combinations of pressure and $\varepsilon$.

Searching for seasonally, rather than continuously, temperate conditions on
Early Mars relaxes the constraints on the pressure and the general composition
of the atmosphere. This is especially true when the orbital parameters of the
planet are substantially different from their modern values and might help to
reconcile the somewhat paradoxical observational evidences inferred from
martian geological records. This desirable outcome is especially evident in
the results presented in Sections 3.1, 3.5 and 3.6.

From Fig. 3, it appears that additional greenhouse gases (such as H2 or CH4)
are not strictly needed, since when obliquity is sufficiently high seasonal
thaws can happen even if the atmosphere is made only of CO2 with a minor
contribution from H2O (using desert-like RH values). Even minor greenhouse
contributions from other gases make the parameter space region in which
seasonal liquid water conditions are possible substantially larger. As an
example, the concentrations of H2 needed to produce a stable, year-long warm-
and-wet Mars seems difficult to produce using reasonable outgassing and
atmospheric escape rates (Ramirez et al., 2014), but considering the seasonal
thermal maxima, as we did here, would substantially lower such a requirement.

Fig. 7 hints to a possible mechanism, different from adding greenhouse gases,
to locally warm up Early Mars, which is varying together the eccentricity, the
obliquity and the argument of perihelion. While the combinations of these
parameters that allow for thaws are unlikely (at least for an atmosphere
without H2 or CH4), their existence suggest a possible pathway for the valley
networks formation that is not as explored as others.

Fig. 8 shows instead the different pressure range required for a fully vs a
seasonally deglaciated planet. Depending on the exact value of other
parameters, seasonal deglaciation requires at the very least 25% less surface
pressure for a given atmospheric composition, in some cases allowing thaws
even at 0.1 bar. The biggest problem in this regard seems to be the
condensation of the atmosphere at surface. Specific recipes would be required
to track the climate evolution under this kind of conditions. Upper limits on
martian paleopressure are uncertain, but generally point to a relatively thin
atmosphere. Kite et al. (2014) estimates a value of about 2 bar from the
crater diameter distribution at the beginning of the Hesperian age (3.6 Gyr
ago), while Edwards & Ehlmann (2015), by studying known martian carbonate
deposits, found that it is improbable that more than $\sim$0.5 bar has been
sequestered at surface. By combining their results with the upper limit on CO2
escape (0.15 bar, Lammer et al., 2013), a reasonable maximum pressure of
$\sim$0.7 bar during the Hesperian can be obtained. Interestingly, our model
BHGH is both stable and able to produce seasonal thaws at around that pressure
value, albeit for a narrow ranges of obliquity and eccentricity values. On the
other hands, several combinations of atmospheric and planetary parameters
among those tested here are able to produce seasonal liquid water conditions
in the 1-to-2 bar range.

###  4.1 Comparison with the observed distribution of valleys

We can now compare our results with the observed latitudinal distribution of
valley networks on the martian surface. Hynek et al. (2010) reported that a
majority of the Noachian valley networks (which represents 91% of the total
valley networks identified by them) can be found between 45∘S and
$5^{\circ}$N, with a peak a few degrees south of the equator, while Hesperian
networks present a more uniform distribution with a first peak at
$35^{\circ}$S and a second one at $5^{\circ}$N. The altitude distribution of
Noachian network midpoints is nearly gaussian and peaks at 1500 m above the
datum (roughly corresponding to the mean altitude of Noachian terrains), while
Hesperian networks are again distributed more regularly with a mean altitude
consistent with 0 m.

Due to the effect of the lapse rate, the latitude at which $T_{\text{max}}$
occurs is nearly always in the northern hemisphere. However, Noachian and
Hesperian valleys are far more often found in the South. Thus, we can check if
sufficiently warm conditions can occur also in the latitudinal regions
associated with the presence of these valleys. This is especially relevant
since liquid water conditions tend to occur at high obliquity values: in these
cases, higher latitudes are subject to more intense $T_{s}$ variations than
the equator and thus deglaciations are favored at high, rather than low,
latitudes.

As it is possible to see in Fig. 9, models B0GH and C0GN struggle to produce
surface conditions compatible with liquid water runoffs at low southern
latitudes, even for the lower $T_{\text{thr}}$ considered. On the other hand,
models BLGH and BHGH allow brine to stay liquid for at least the summer season
across the entire southern hemisphere when $P_{0}$ is above 1.77 and 1.33 bar,
respectively, while they require $\gtrsim 2$ bar of pressure to let pure water
ice to melt in the $0^{\circ}$S - $45^{\circ}$S region.

Concerning the dependences on parameters other than the obliquity (not shown
in the Figure), the impact of eccentricity and surface albedo is smaller,
while that of the cloud coverage is larger, thus confirming the trends
encountered in Section 3. For model BHGH, thaws are possible across most of
the southern hemisphere for $P_{0}>1$ bar if we consider brine, and $P_{0}>2$
bar if we consider water independently of the chosen eccentricity. Pressure
limits move to 2 bar for brine and 3 bar for water in the BLGH case. C0GN
requires more than 4 bar to allow for water runoffs in the south, while B0GH
never produces datum temperatures at or above 273 K. Land albedo also has a
minor impact on southern summer temperatures, despite the fact that oceans are
mostly absent in that region. Again, the main determinant are the surface
pressure and the composition. For example, in model BHGH, lowering the albedo
from 0.41 to 0.11 only lowers the required pressure from 2.0 to 1.33 bar (from
3 to 2 bar for BLGH). On the other hand, the cloud fraction has a larger
impact due to their important greenhouse contribution in our model. Cloud
cover fraction at or above 0.5 are capable of warming the southern hemisphere
to a sufficient level to allow thaws in all models, albeit above different
$P_{0}$ threshold values.

At variance with Hesperian valley networks, Noachian ones lie at 1500 m of
altitude, as previously noted. This places a tighter constraint on the climate
state of the planet, since at that height the surface temperature is 7-8 K
lower than at the datum. We thus repeated the analysis including this effect,
finding that it increases the pressure threshold for seasonal thaws by
$\sim$33%. Again, the BHGH model is the best suited to explain the
observational evidences, while B0GH and C0GN can produce brine runoffs but not
water ones.

An interesting aspect of searching for seasonal thaws is their possible role
in the formation process of Equatorial Layered Deposits (ELDs, Schmidt et al.,
2021). ELDs are stratified deposits mostly found in craters and basins in the
equatorial region (especially in Terra Arabia) and probably associated with
water activity. A variety of processes have been proposed as their source (see
e.g. Schmidt et al., 2022b, Table 1 and 2). The thinning-thickening sequence
of the deposits might be explained either by a relatively strong seasonal
variability (Schmidt, G., private communication) or by secular changes in the
hydrological cycle forced by obliquity and eccentricity variations, depending
on the time scale associated with this feature and that is currently under
investigation. Morphological evidences found in the Gale Crater evaporitic
basin seems to point toward a sustained wet-dry cycling compatible with
sustained seasonal variations (Rapin et al., 2023) as those found in some of
our simulations.

##  5 Conclusions and future prospects

In this work we have studied the combined impact of varying the surface
pressure and other martian planetary parameters (obliquity, eccentricity,
argument of perihelion, surface albedo and fraction of clouds) on the seasonal
temperature changes in order to identify the conditions conductive to regional
ice melting. We tested five different atmospheric compositions with varying
H2O ($RH=$ 0%, 20% and 40%) and CH4 (0%, 0.1% and 1%). We employed a seasonal-
latitudinal climate model (ESTM) coupled with an up-to-date RT code (EOS),
including the effects of the band-averaged altimetry, a northern ocean with
either 150 or 550 m of GEL and checking for possible CO2 condensation at
surface. In total, we run $\sim$104 cases (i.e. different combinations of
input parameters).

Our main findings can be summarized as follows:

  1. 1.

If both the planet obliquity and datum pressure are sufficiently high
($\varepsilon\geq 50^{\circ}$ and $P_{0}\geq 2.0$ bar) then no other
additional greenhouse gases other than CO2 and desert-like levels of H2O are
needed to produce seasonal thaws in the northern hemisphere. These conditions
also guarantee that no CO2 condensation happens at the surface. However, the
same models fail to produce seasonal melting in the southern highlands.

  2. 2.

Adding a 0.1% of CH4 lowers the requirements for northern thaws to 1.33 bar
(at $\varepsilon=45^{\circ}$) and simultaneously allows for melting in the
$60^{\circ}$S - $90^{\circ}$S region. Increasing the pressure to 3.5 bars or
above melts also the region near the equator. If $f_{\text{CH4}}$ is instead
set to 1%, both northern and southern thaws are possible independently of the
obliquity when $P_{0}\geq 1.33$ bar.

  3. 3.

Increasing the eccentricity allows for thaws and stabilizes the atmosphere
when combined with $\varepsilon\geq 50^{\circ}$ and $\omega\sim 270^{\circ}$
(i.e. when Mars is at perihelion during the northern hemisphere’s summer
solstice). The effect of eccentricity and obliquity alone can change
$T_{\text{max}}$ by 29-52 K, depending on the specific choice of the argument
of perihelion. This effect is important if the planet is near the threshold
conditions for hosting liquid water.

  4. 4.

Changing the albedo of lands has some impact in the models without CH4, while
it is less important when CH4 is present. If $f_{\text{CH4}}$ is set to 1%,
both northern and southern melting are possible at $P_{0}\geq 2$ bars
independently of the land albedo. This is especially interesting when we
consider that the southern hemisphere has a very low ocean fraction.

  5. 5.

Clouds are difficult to treat and have been included in this work by
parameterizing their contributions to OLR and TOA albedo. Clouds calibrated on
the Forget et al. (2013) model (i.e. CO2-ice clouds) give only a modest
contribution in terms of both average temperature and $T_{\text{max}}$ when
coupled with a dry CO2-dominated atmosphere. However, they can prevent the
atmospheric collapse if fractional coverage is higher than 0.8. Clouds
calibrated on the modern Earth model (i.e. H2O clouds) are more efficient at
warming Mars. We found that the minimum cloud coverage to prevent the
atmospheric collapse is 0.5 if there is no CH4, 0.2 when $f_{\text{CH4}}$ is
at 0.1% and are not needed when $f_{\text{CH4}}$ is at 1%.

  6. 6.

If melting is possible, then mild conditions last usually more than 15% of
martian year. While we did not calculate the amount of liquid water that we
can obtain, we are confident that it is not negligible (as if the thawing
season were very short).

  7. 7.

Changing the GEL level has a relatively small impact on the results. Notably,
the transition from the region of the parameter space in which seasonal thaws
are not allowed to that in which they can happen is sharper. This is due to
the fact that a larger ocean causes a larger swing in the average surface
albedo when it melts.

  8. 8.

If we relax the thawing conditions adopting the freezing temperature of brine
instead of that of pure water, the combinations of planetary and atmospheric
parameters that allow melting significantly widen. Seasonal brine runoffs in
the northern hemisphere are allowed also in dry CO2-dominated atmospheres with
high ($\geq 0.8$) CO2-ice cloud fractions for pressures above 1.33 bars.

To conclude, we confirm the previous results indicating the necessity of other
greenhouse gases other than CO2 and H2O to allow for an hydrological cycle
that involves the southern highlands. However, we show that depending on the
combinations of orbital ($\varepsilon$, $e$ and $\omega$) and planetary (land
albedo, cloud coverage) parameters, this contribution can be significantly
lower than usually expected (5-20% of H2, e.g. Ramirez et al., 2014) and
probably easier to explain within the other known geological constraints of
Mars. In particular, to our knowledge, no previous investigations on the
interplay between eccentricity, obliquity and the argument of perihelion for
different atmospheric models were performed. An investigation on the interplay
between eccentricity, obliquity and the argument of perihelion has been
recently performed to explain the formation of gullies in the recent past of
Mars (Dickson et al., 2023), but here we explored a larger portion of the
parameter space and referred to Early Mars conditions. We also underline the
advantages of using seasonal-latitudinal, rather than single-column, climate
models to combine together computational efficiency and predictive power.

There are several possible avenues for future developments. First of all, we
limited ourselves to a very narrow range of atmospheric chemical compositions,
all of which in equilibrium. Using the same methodology, it would be possible
to investigate the effects of pulses of e.g. volcanic gases for different
combinations of martian orbital and surface parameters. Second, we took into
consideration a simple parameterization of clouds, rather than including them
into the RT model. Improving the RT modeling of clouds is very important,
since as we have shown, their contribution to planetary warming is crucial
when the concentrations of other greenhouse gases is low. Third, it would be
desirable to calculate the amount of liquid water that is made available
during the seasonal thaws. A precise assessment of precipitation rates
requires 3D GCMs, while an upgraded version of ESTM would be able to provide
approximate latitudinally-averaged estimates that can lend additional
robustness to these findings.

P.S. wants to thank Gene Walter Schmidt for the insightful discussion
concerning the ELDs. We acknowledge support by the Italian Space Agency with
the Life in Space project (ASI N. 2019-3-U.0), OGS and CINECA with the HPC-
TRES program award N. 2022-02 and INAF with the mini-grant titled Radiative
models for the paleoatmospheres of Mars, Earth and Venus (F.O.
1.05.121.04.03).

##  Appendix A Lapse rate calculation

Figure 10: Four examples of vertical PT profiles calculated using the
procedure described in this Section. Black and red lines are calculated
specifying a $T_{s}$ equal to 280 K and 320 K, respectively. Solid refers to
cases with $RH=100\%$ and dotted to cases with $RH=20\%$. The grey dashed line
is the condensation curve for CO2 from the Kasting (1991) formulation.

Both the vertical pressure-temperature profile and the altimetric adjustment
have been calculated using a modified version of the moist pseudo-adiabatic
lapse rate $\Gamma$ provided by the American Meteorological Society
Glossary555https://glossary.ametsoc.org/wiki/:

$\Gamma=g\,\frac{(1+r_{v})(1+\frac{L_{v}r_{v}}{R_{d}T})}{c_{pd}+r_{v}c_{pv}+\frac{L_{v}^{2}r_{v}(\phi+r_{v})}{R_{d}T^{2}}}$
(A1)

where $g$ is the gravitational acceleration, $r_{v}$ is the mass mixing ratio
of the condensible, $L_{v}$ is the latent heat of vaporization of the
condensible, $R_{d}$ is the gas constant of the dry gas, $c_{pd}$ and $c_{pv}$
are the specific isobaric heats of the dry gas and the condensible,
respectively, and $\phi$ is the ratio between the molecular masses of the
condensible and the dry gas. Dry gas properties are the weighted average of
the component gases considered (here, CO2 and CH4). From a theoretical point
of view, this lapse rate is valid only if the $RH$ of the condensible is 100%,
while the dry lapse rate should be used if this condition is not true.
However, we do know that the real lapse rate in the Earth atmosphere is
generally in-between the moist and the dry lapse rate due to inhomogeneities
in the convecting air mass. To represent this fact, the $r_{v}$ term in the
H2O condensing part of the troposphere is scaled by the $RH$ value specified
in Tab. 3.

The condensible is H2O in the lower troposphere and CO2 in the upper
troposphere. In practice, we integrate Eq. A1 on a 100 m-spaced grid in height
(10 m-spaced grid when we calculate the altimetric correction) starting from
the surface and up to the height corresponding to a specified pressure (i.e.
the $P_{\text{TOA}}$) or temperature (i.e. the stratospheric temperature). The
terms $r_{v}$, $L_{v}$, $c_{pd}$ and $c_{pv}$ depend on temperature and are
updated accordingly by iterating on Eq. A1. At each step we check if CO2
condensation occurs and if this is the case, we begin treating CO2 as the
condensible in place of water. Strictly speaking H2O condensation continue to
happen, but since it contributes a negligible amount of latent heat with
respect to CO2 in this region, we neglect it. The altitude-temperature profile
obtained is then converted into a pressure-temperature profile using the
barometric formula for a non-isothermal atmosphere:

$P_{i+1}=P_{i}\,\biggl{(}\frac{T_{i+1}}{T_{i}}\biggl{)}^{-g/(R_{m}\Gamma)}$
(A2)

where $T_{i}$ is the temperature of the $i$th layer and $R_{m}$ is the gas
constant of the atmospheric mix (dry plus moist). Fig. 10, left panel, shows
four examples of PT profiles calculated using this procedure. There is a small
offset with respect to the CO2 condensation curve due to the fact that these
profiles refer to the case with $f_{\text{CH4}}=1\%$.

Using temperature-dependent latent and isobaric specific heats allows us to
capture the main contributions to deviations from an ideal gas lapse rate. In
Simonetti et al. (2022) we adopted the formulation of Kasting (1991), that
includes also the effects of a varying CO2 compressibility. However, that
formulation approximate the pressure of the condensible to zero, which is a
reasonable assumption only for relatively cold ($T_{s}\leq 280$ K) atmospheres
but breaks down at higher temperatures. Low total surface pressures make the
problem worse, and in general the error caused by not considering the changing
gas compressibility are larger than the error caused by neglecting the vapor
pressure, at least in the pressure and temperature ranges that we worked with.

## References

  * Andrews-Hanna et al. (2008) Andrews-Hanna, J. C., Zuber, M. T., & Banerdt, W. B. 2008, Nature, 453, 1212, doi: 10.1038/nature07011
  * Balme et al. (2020) Balme, M. R., Gupta, S., Davis, J. M., et al. 2020, Journal of Geophysical Research (Planets), 125, e06244, doi: 10.1029/2019JE006244
  * Baranov & Vigasin (1999) Baranov, Y. I., & Vigasin, A. A. 1999, Journal of Molecular Spectroscopy, 193, 319, doi: 10.1006/jmsp.1998.7743
  * Batalha et al. (2015) Batalha, N., Domagal-Goldman, S. D., Ramirez, R., & Kasting, J. F. 2015, Icarus, 258, 337, doi: 10.1016/j.icarus.2015.06.016
  * Batalha et al. (2016) Batalha, N. E., Kopparapu, R. K., Haqq-Misra, J., & Kasting, J. F. 2016, Earth and Planetary Science Letters, 455, 7, doi: 10.1016/j.epsl.2016.08.044
  * Biasiotti et al. (2022) Biasiotti, L., Simonetti, P., Vladilo, G., et al. 2022, MNRAS, 514, 5105, doi: 10.1093/mnras/stac1642
  * Bultel et al. (2019) Bultel, B., Viennet, J.-C., Poulet, F., Carter, J., & Werner, S. C. 2019, Journal of Geophysical Research (Planets), 124, 989, doi: 10.1029/2018JE005845
  * Carr & Head (2003) Carr, M. H., & Head, J. W. 2003, Journal of Geophysical Research (Planets), 108, 5042, doi: 10.1029/2002JE001963
  * Carr & Head (2010) —. 2010, Earth and Planetary Science Letters, 294, 185, doi: 10.1016/j.epsl.2009.06.042
  * Chassefière et al. (2013) Chassefière, E., Langlais, B., Quesnel, Y., & Leblanc, F. 2013, Journal of Geophysical Research (Planets), 118, 1123, doi: 10.1002/jgre.20089
  * Chester & Jickells (2012) Chester, R., & Jickells, T. 2012, Chapter 7: Descriptive oceanography: water-column parameters (Chichester, West Sussex, UK: Wiley/Blackwell)
  * Clark et al. (2021) Clark, B., Kolb, V., Steele, A., et al. 2021, Life, 11, doi: 10.3390/life11060539
  * Clough et al. (1989) Clough, S. A., Kneizys, F. X., & Davies, R. W. 1989, Atmospheric Research, 23, 229, doi: 10.1016/0169-8095(89)90020-3
  * Clough et al. (2005) Clough, S. A., Shephard, M. W., Mlawer, E. J., et al. 2005, J. Quant. Spec. Radiat. Transf., 91, 233, doi: 10.1016/j.jqsrt.2004.05.058
  * Colaprete & Toon (2003) Colaprete, A., & Toon, O. B. 2003, Journal of Geophysical Research (Planets), 108, 5025, doi: 10.1029/2002JE001967
  * Costard et al. (2019) Costard, F., Séjourné, A., Lagain, A., et al. 2019, Journal of Geophysical Research (Planets), 124, 1840, doi: 10.1029/2019JE006008
  * Craddock & Howard (2002) Craddock, R. A., & Howard, A. D. 2002, Journal of Geophysical Research (Planets), 107, 5111, doi: 10.1029/2001JE001505
  * DeWitt et al. (2009) DeWitt, H. L., Trainer, M. G., Pavlov, A. A., et al. 2009, Astrobiology, 9, 447, doi: 10.1089/ast.2008.0289
  * di Achille & Hynek (2010) di Achille, G., & Hynek, B. M. 2010, Nature Geoscience, 3, 459, doi: 10.1038/ngeo891
  * Dickson et al. (2023) Dickson, J. L., Palumbo, A. M., Head, J. W., et al. 2023, Science, 380, 1363, doi: 10.1126/science.abk2464
  * Edwards & Ehlmann (2015) Edwards, C. S., & Ehlmann, B. L. 2015, Geology, 43, 863, doi: 10.1130/G36983.1
  * Essack et al. (2020) Essack, Z., Seager, S., & Pajusalu, M. 2020, ApJ, 898, 160, doi: 10.3847/1538-4357/ab9cba
  * Etiope et al. (2013) Etiope, G., Ehlmann, B. L., & Schoell, M. 2013, Icarus, 224, 276, doi: 10.1016/j.icarus.2012.05.009
  * Fairén (2010) Fairén, A. G. 2010, Icarus, 208, 165, doi: 10.1016/j.icarus.2010.01.006
  * Fairén et al. (2012) Fairén, A. G., Haqq-Misra, J. D., & McKay, C. P. 2012, A&A, 540, A13, doi: 10.1051/0004-6361/201118527
  * Fassett & Head (2011) Fassett, C. I., & Head, J. W. 2011, Icarus, 211, 1204, doi: 10.1016/j.icarus.2010.11.014
  * Fastook & Head (2014) Fastook, J. L., & Head, J. W. 2014, Planet. Space Sci., 91, 60, doi: 10.1016/j.pss.2013.12.002
  * Forget & Pierrehumbert (1997) Forget, F., & Pierrehumbert, R. T. 1997, Science, 278, 1273, doi: 10.1126/science.278.5341.1273
  * Forget et al. (2013) Forget, F., Wordsworth, R., Millour, E., et al. 2013, Icarus, 222, 81, doi: 10.1016/j.icarus.2012.10.019
  * Gordon et al. (2022) Gordon, I. E., Rothman, L. S., Hargreaves, R. J., et al. 2022, J. Quant. Spec. Radiat. Transf., 277, 107949, doi: 10.1016/j.jqsrt.2021.107949
  * Gough (1981) Gough, D. O. 1981, Sol. Phys., 74, 21, doi: 10.1007/BF00151270
  * Grimm & Heng (2015) Grimm, S. L., & Heng, K. 2015, ApJ, 808, 182, doi: 10.1088/0004-637X/808/2/182
  * Grimm et al. (2021) Grimm, S. L., Malik, M., Kitzmann, D., et al. 2021, ApJS, 253, 30, doi: 10.3847/1538-4365/abd773
  * Grott et al. (2011) Grott, M., Morschhauser, A., Breuer, D., & Hauber, E. 2011, Earth and Planetary Science Letters, 308, 391, doi: 10.1016/j.epsl.2011.06.014
  * Gruszka & Borysow (1997) Gruszka, M., & Borysow, A. 1997, Icarus, 129, 172
  * Gruszka & Borysow (1998) —. 1998, Molecular Physics, 93, 1007, doi: 10.1080/002689798168709
  * Haberle et al. (2019) Haberle, R. M., Zahnle, K., Barlow, N. G., & Steakley, K. E. 2019, Geophys. Res. Lett., 46, 13,355, doi: 10.1029/2019GL084733
  * Halevy & Head (2014) Halevy, I., & Head, James W., I. 2014, Nature Geoscience, 7, 865, doi: 10.1038/ngeo2293
  * Halevy et al. (2007) Halevy, I., Zuber, M. T., & Schrag, D. P. 2007, Science, 318, 1903, doi: 10.1126/science.1147039
  * Haqq-Misra & Hayworth (2022) Haqq-Misra, J., & Hayworth, B. P. C. 2022, , 3, 32, doi: 10.3847/PSJ/ac49eb
  * Hartmann & Neukum (2001) Hartmann, W. K., & Neukum, G. 2001, Space Sci. Rev., 96, 165, doi: 10.1023/A:1011945222010
  * Hayworth et al. (2020) Hayworth, B. P. C., Kopparapu, R. K., Haqq-Misra, J., et al. 2020, Icarus, 345, 113770, doi: 10.1016/j.icarus.2020.113770
  * Heng et al. (2018) Heng, K., Malik, M., & Kitzmann, D. 2018, ApJS, 237, 29, doi: 10.3847/1538-4365/aad199
  * Hoke et al. (2011) Hoke, M. R. T., Hynek, B. M., & Tucker, G. E. 2011, Earth and Planetary Science Letters, 312, 1, doi: 10.1016/j.epsl.2011.09.053
  * Huang et al. (2019) Huang, C. J., Qiao, F., Chen, S., Xue, Y., & Guo, J. 2019, Journal of Geophysical Research: Oceans, 124, 4480–4491
  * Hynek et al. (2010) Hynek, B. M., Beach, M., & Hoke, M. R. T. 2010, Journal of Geophysical Research (Planets), 115, E09008, doi: 10.1029/2009JE003548
  * Ivanov et al. (2017) Ivanov, M. A., Erkeling, G., Hiesinger, H., Bernhardt, H., & Reiss, D. 2017, Planet. Space Sci., 144, 49, doi: 10.1016/j.pss.2017.05.012
  * Johnson et al. (2008) Johnson, S. S., Mischna, M. A., Grove, T. L., & Zuber, M. T. 2008, Journal of Geophysical Research (Planets), 113, E08005, doi: 10.1029/2007JE002962
  * Kamada et al. (2021) Kamada, A., Kuroda, T., Kasaba, Y., Terada, N., & Nakagawa, H. 2021, Icarus, 368, 114618, doi: 10.1016/j.icarus.2021.114618
  * Kamada et al. (2020) Kamada, A., Kuroda, T., Kasaba, Y., et al. 2020, Icarus, 338, 113567, doi: 10.1016/j.icarus.2019.113567
  * Kasting (1991) Kasting, J. F. 1991, Icarus, 94, 1, doi: 10.1016/0019-1035(91)90137-I
  * Kerber et al. (2015) Kerber, L., Forget, F., & Wordsworth, R. 2015, Icarus, 261, 133, doi: 10.1016/j.icarus.2015.08.011
  * Kite et al. (2017) Kite, E. S., Gao, P., Goldblatt, C., et al. 2017, Nature Geoscience, 10, 737, doi: 10.1038/ngeo3033
  * Kite et al. (2014) Kite, E. S., Williams, J.-P., Lucas, A., & Aharonson, O. 2014, Nature Geoscience, 7, 335, doi: 10.1038/ngeo2137
  * Kitzmann et al. (2013) Kitzmann, D., Patzer, A. B. C., & Rauer, H. 2013, A&A, 557, A6, doi: 10.1051/0004-6361/201220025
  * Lammer et al. (2013) Lammer, H., Chassefière, E., Karatekin, Ö., et al. 2013, Space Sci. Rev., 174, 113, doi: 10.1007/s11214-012-9943-8
  * Laskar et al. (2004) Laskar, J., Correia, A. C. M., Gastineau, M., et al. 2004, Icarus, 170, 343, doi: 10.1016/j.icarus.2004.04.005
  * Malik et al. (2019) Malik, M., Kitzmann, D., Mendonça, J. M., et al. 2019, AJ, 157, 170, doi: 10.3847/1538-3881/ab1084
  * Malik et al. (2017) Malik, M., Grosheintz, L., Mendonça, J. M., et al. 2017, AJ, 153, 56, doi: 10.3847/1538-3881/153/2/56
  * Mangold et al. (2004) Mangold, N., Quantin, C., Ansan, V., Delacourt, C., & Allemand, P. 2004, Science, 305, 78, doi: 10.1126/science.1097549
  * Marinova et al. (2008) Marinova, M. M., Aharonson, O., & Asphaug, E. 2008, Nature, 453, 1216, doi: 10.1038/nature07070
  * Mischna et al. (2013) Mischna, M. A., Baker, V., Milliken, R., Richardson, M., & Lee, C. 2013, Journal of Geophysical Research (Planets), 118, 560, doi: 10.1002/jgre.20054
  * Mlawer et al. (2012) Mlawer, E. J., Payne, V. H., Moncet, J. L., et al. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2520, doi: 10.1098/rsta.2011.0295
  * Mouginot et al. (2012) Mouginot, J., Pommerol, A., Beck, P., Kofman, W., & Clifford, S. M. 2012, Geophys. Res. Lett., 39, L02202, doi: 10.1029/2011GL050286
  * Murante et al. (2020) Murante, G., Provenzale, A., Vladilo, G., et al. 2020, MNRAS, 492, 2638, doi: 10.1093/mnras/stz3529
  * North et al. (1981) North, G. R., Cahalan, R. F., & Coakley, James A., J. 1981, Reviews of Geophysics and Space Physics, 19, 91, doi: 10.1029/RG019i001p00091
  * Orosei et al. (2018) Orosei, R., Lauro, S. E., Pettinelli, E., et al. 2018, Science, 361, 490, doi: 10.1126/science.aar7268
  * Perrin & Hartmann (1989) Perrin, M. Y., & Hartmann, J. M. 1989, J. Quant. Spec. Radiat. Transf., 42, 311, doi: 10.1016/0022-4073(89)90077-0
  * Pierrehumbert (2010) Pierrehumbert, R. T. 2010, Principles of Planetary Climate (Dover)
  * Plaut et al. (2007) Plaut, J. J., Picardi, G., Safaeinili, A., et al. 2007, Science, 316, 92, doi: 10.1126/science.1139672
  * Ramirez (2017) Ramirez, R. M. 2017, Icarus, 297, 71, doi: 10.1016/j.icarus.2017.06.025
  * Ramirez & Craddock (2018) Ramirez, R. M., & Craddock, R. A. 2018, Nature Geoscience, 11, 230, doi: 10.1038/s41561-018-0093-9
  * Ramirez & Kasting (2017) Ramirez, R. M., & Kasting, J. F. 2017, Icarus, 281, 248, doi: 10.1016/j.icarus.2016.08.016
  * Ramirez et al. (2014) Ramirez, R. M., Kopparapu, R., Zugger, M. E., et al. 2014, Nature Geoscience, 7, 59, doi: 10.1038/ngeo2000
  * Rapin et al. (2023) Rapin, W., Dromart, G., Clark, B. C., et al. 2023, Nature, 620, 299, doi: 10.1038/s41586-023-06220-3
  * Richardson et al. (2007) Richardson, M. I., Toigo, A. D., & Newman, C. E. 2007, Journal of Geophysical Research (Planets), 112, E09001, doi: 10.1029/2006JE002825
  * Salese et al. (2019) Salese, F., Pondrelli, M., Neeseman, A., Schmidt, G., & Ori, G. G. 2019, Journal of Geophysical Research (Planets), 124, 374, doi: 10.1029/2018JE005802
  * Sauterey et al. (2022) Sauterey, B., Charnay, B., Affholder, A., Mazevet, S., & Ferrière, R. 2022, Nature Astronomy, 6, 1263, doi: 10.1038/s41550-022-01786-w
  * Scheller et al. (2021) Scheller, E. L., Ehlmann, B. L., Hu, R., Adams, D. J., & Yung, Y. L. 2021, Science, 372, 56, doi: 10.1126/science.abc7717
  * Schmidt et al. (2022a) Schmidt, F., Way, M. J., Costard, F., et al. 2022a, Proceedings of the National Academy of Science, 119, e2112930118, doi: 10.1073/pnas.2112930118
  * Schmidt et al. (2022b) Schmidt, G., Luzzi, E., Rossi, A. P., et al. 2022b, Journal of Geophysical Research (Planets), 127, e07320, doi: 10.1029/2022JE007320
  * Schmidt et al. (2021) Schmidt, G., Pondrelli, M., Salese, F., et al. 2021, Journal of Geophysical Research (Planets), 126, e06974, doi: 10.1029/2021JE006974
  * Segura et al. (2008) Segura, T. L., Toon, O. B., & Colaprete, A. 2008, Journal of Geophysical Research (Planets), 113, E11007, doi: 10.1029/2008JE003147
  * Segura et al. (2002) Segura, T. L., Toon, O. B., Colaprete, A., & Zahnle, K. 2002, Science, 298, 1977, doi: 10.1126/science.1073586
  * Seidelmann et al. (2007) Seidelmann, P. K., Archinal, B. A., A’Hearn, M. F., et al. 2007, Celestial Mechanics and Dynamical Astronomy, 98, 155, doi: 10.1007/s10569-007-9072-y
  * Sergeev et al. (2022) Sergeev, D. E., Fauchez, T. J., Turbet, M., et al. 2022, , 3, 212, doi: 10.3847/PSJ/ac6cf2
  * Silva et al. (2017) Silva, L., Vladilo, G., Murante, G., & Provenzale, A. 2017, MNRAS, 470, 2270, doi: 10.1093/mnras/stx1396
  * Simonetti et al. (2022) Simonetti, P., Vladilo, G., Silva, L., et al. 2022, ApJ, 925, 105, doi: 10.3847/1538-4357/ac32ca
  * Sleep (1994) Sleep, N. H. 1994, J. Geophys. Res., 99, 5639, doi: 10.1029/94JE00216
  * Smith et al. (2003) Smith, D., Neumann, G., Arvidson, R. E., Guinness, E. A., & Slavney, S. 2003, Mars Global Surveyor Laser Altimeter Mission Experiment Gridded Data Record, MGS-M-MOLA-5-MEGDR-L3-V1.0, NASA Planetary Data System
  * Sneep & Ubachs (2005) Sneep, M., & Ubachs, W. 2005, J. Quant. Spec. Radiat. Transf., 92, 293, doi: 10.1016/j.jqsrt.2004.07.025
  * Squyres & Kasting (1994) Squyres, S. W., & Kasting, J. F. 1994, Science, 265, 744, doi: 10.1126/science.265.5173.744
  * Tian et al. (2010) Tian, F., Claire, M. W., Haqq-Misra, J. D., et al. 2010, Earth and Planetary Science Letters, 295, 412, doi: 10.1016/j.epsl.2010.04.016
  * Trainer et al. (2006) Trainer, M. G., Pavlov, A. A., Dewitt, H. L., et al. 2006, Proceedings of the National Academy of Science, 103, 18035, doi: 10.1073/pnas.0608561103
  * Turbet et al. (2020a) Turbet, M., Boulet, C., & Karman, T. 2020a, Icarus, 346, 113762, doi: 10.1016/j.icarus.2020.113762
  * Turbet et al. (2020b) Turbet, M., Gillmann, C., Forget, F., et al. 2020b, Icarus, 335, 113419, doi: 10.1016/j.icarus.2019.113419
  * Urata & Toon (2013) Urata, R. A., & Toon, O. B. 2013, Icarus, 226, 229, doi: 10.1016/j.icarus.2013.05.014
  * Villanueva et al. (2015) Villanueva, G. L., Mumma, M. J., Novak, R. E., et al. 2015, Science, 348, 218, doi: 10.1126/science.aaa3630
  * Vladilo et al. (2015) Vladilo, G., Silva, L., Murante, G., Filippi, L., & Provenzale, A. 2015, ApJ, 804, 50, doi: 10.1088/0004-637X/804/1/50
  * Šrámek & Zhong (2012) Šrámek, O., & Zhong, S. 2012, Journal of Geophysical Research (Planets), 117, E01005, doi: 10.1029/2011JE003867
  * Wagner & Kretzschmar (2008) Wagner, W., & Kretzschmar, H.-J. 2008, International Steam Tables - Properties of Water and Steam based on the Industrial Formulation IAPWS-IF97 (Springer)
  * Webster et al. (2013) Webster, C. R., Mahaffy, P. R., Flesch, G. J., et al. 2013, Science, 341, 260, doi: 10.1126/science.1237961
  * Webster et al. (2015) Webster, C. R., Mahaffy, P. R., Atreya, S. K., et al. 2015, Science, 347, 415, doi: 10.1126/science.1261713
  * Wilhelms & Squyres (1984) Wilhelms, D. E., & Squyres, S. W. 1984, Nature, 309, 138, doi: 10.1038/309138a0
  * Williams & Kasting (1997) Williams, D. M., & Kasting, J. F. 1997, Icarus, 129, 254, doi: 10.1006/icar.1997.5759
  * Wordsworth et al. (2010) Wordsworth, R., Forget, F., & Eymet, V. 2010, Icarus, 210, 992, doi: 10.1016/j.icarus.2010.06.010
  * Wordsworth et al. (2013) Wordsworth, R., Forget, F., Millour, E., et al. 2013, Icarus, 222, 1, doi: 10.1016/j.icarus.2012.09.036
  * Wordsworth et al. (2017) Wordsworth, R., Kalugina, Y., Lokshtanov, S., et al. 2017, Geophys. Res. Lett., 44, 665, doi: 10.1002/2016GL071766
  * Wordsworth et al. (2021) Wordsworth, R., Knoll, A. H., Hurowitz, J., et al. 2021, Nature Geoscience, 14, 127, doi: 10.1038/s41561-021-00701-8
  * Wordsworth (2016) Wordsworth, R. D. 2016, Annual Review of Earth and Planetary Sciences, 44, 381, doi: 10.1146/annurev-earth-060115-012355
  * Wordsworth et al. (2015) Wordsworth, R. D., Kerber, L., Pierrehumbert, R. T., Forget, F., & Head, J. W. 2015, Journal of Geophysical Research (Planets), 120, 1201, doi: 10.1002/2015JE004787
  * Yang et al. (2020) Yang, F., He, Q., Huang, J., et al. 2020, Bulletin of the American Meteorological Society, 102, E1172, doi: 10.1175/BAMS-D-20-0236.1
  * Zahnle (1986) Zahnle, K. J. 1986, J. Geophys. Res., 91, 2819, doi: 10.1029/JD091iD02p02819
  * Zhong & Zuber (2001) Zhong, S., & Zuber, M. T. 2001, Earth and Planetary Science Letters, 189, 75, doi: 10.1016/S0012-821X(01)00345-4