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

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 CO$_2$ or H$_2$O clouds. We find that 1.3-to-2.0 bar CO$_2$-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. H$_2$, CH$_4$) or transient heat-injecting phenomena (e.g. asteroid impacts, volcanic eruptions). Moderate amounts (0.1 to 1\%) of CH$_4$ significantly widens the parameter space region in which seasonal thaws are possible.


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 evidence 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 (∼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 single-column or seasonal-latitudinal. Single-column 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).Seasonallatitudinal EBMs add the latitudinal horizontal direction and, as such, are able to estimate the evolution over time of local conditions on a planet surface (North et al. 1981).A simplified horizontal heat diffusion might or might not be present, but only the former class of models can be applied to planets with non-negligible atmospheres and take into consideration some simple climate feedbacks, like the ice-albedo feedback.While historically applied to the study of fast rotating planets, seasonal-latitudinal EBMs have been successfully employed to model also tidally-locked ones (e.g.Kite et al. 2011;Haqq-Misra & Hayworth 2022).Recently, Okuya et al. (2019) proposed an EBM that also considers the longitudinal heat diffusion and is thus capable of producing full surface temperature maps.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 CO 2 -H 2 O atmosphere, but in single-column 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 SO 2 /H 2 S (Halevy et al. 2007;Johnson et al. 2008;Mischna et al. 2013) or H 2 /CH 4 (Ramirez et al. 2014;Wordsworth et al. 2017;Ramirez 2017;Haberle et al. 2019;Turbet et al. 2020a), satisfactorily solving the issues associated with CO 2 -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 CO 2 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 H 2 O 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(Segura et al. , 2008)), volcanic eruptions (Halevy & Head 2014), the destabilization of previously formed CH 4 clathrates due to chaotic obliquity changes (Kite et al. 2017) or episodic injections of H 2 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 SO 2 and H 2 S may have cooled, rather than warmed, the planet (Tian et al. 2010;Kerber et al. 2015); it is uncertain if CH 4 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.Wordsworth et al. 2015).Also, in order to prevent the locking of the hydrosphere in the southern highlands (which act as an efficient cold trap), an even higher global average temperature is required, which puts an extra constraint on the atmospheric chemical composition.The cold-andicy 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 (≳ 200 m) global equivalent layer (GEL) of water, sig-nificant 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 singlecolumn 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 single-column 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. 2020Kamada et al. , 2021)).In this work, we make the case for the employment of seasonal-latitudinal EBMs coupled with a flexible RT code as the natural next step along the path inaugurated by single-column 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 CO 2 , H 2 O and CH 4 ) 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 seasonal-latitudinal EBMs, either non-diffusive (e.g.  1. Orbital and planetary parameters adopted for the Early Mars in this work when not otherwise specified.Armstrong et al. 2004;Kite et al. 2013) or diffusive (Ramirez et al. 2020;Hayworth et al. 2020) is not entirely new, but it is still largely unexplored and promises to yield new results by taking the best of both 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.

ESTM
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 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, φ, either directly or indirectly, through their dependence on T .The output of ESTM is T as a function of φ 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 φ.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 Φ in such a way that N = 2πR 2 Φ cos φ is the net rate of energy transport across a circle of constant φ 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.Φ ≡ −D(∂T /∂φ).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, Ω, 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 µ = 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).

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 bandaveraged 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°r esolution 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, panel (a).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 ≤10 K across most of the martian surface, reaching ∼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 The yearly global average planetary temperature at the altimetric datum ⟨T 0⟩ 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 marker represents an increase of 33% of the pressure with respect to the preceding one.Grey dashed and dotted lines mark the 273 K and the 252 K level, respectively.
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, panel (b), 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 marker 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 ∼13 K, while it is smallest in the dry, CO 2 -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 ∼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 ∼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 ∼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 ∼0.4%, thus the combined effect amounts to ∼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.

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 (∼4.1-3.0Gyr 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, panel (a), and correspond to the global ocean coverage of 0.16 and 0.30, respectively.

Surface ices
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 selfconsistently 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 α l (see Sect. 3.3) up to values consistent with a permanent ice sheet of 50% (α l = 0.41).On the other hand, as far as oceans are concerned, we kept the temperaturedependent 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 ∼2 K at most for models using 150 m of water GEL and ∼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.

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 (µ = 0.5), the slope of the linear dependence between µ 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 depen- dence 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).
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 CO 2 -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, panel (b).
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 CO 2 -ice clouds and H 2 O Earthlike 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.4we explored also the effect of different cloud coverages in the [0, 1] range.

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(Malik et al. , 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 CO 2 , H 2 O and CH 4 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 CO 2 absorption lines have been modeled following Perrin & Hartmann (1989) and cutting the wings at 500 cm −1 from the line centre, while H 2 O and CH 4 lines have been modeled using a Voigt distribution cut at 25 cm −1 from the line centre.For H 2 O, 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 CO 2 -CO 2 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 CO 2 -CH 4 CIAs in the 0-2000 cm −1 range (Wordsworth et al. 2017).Concerning H 2 O, 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 H 2 O continuum is calculated for a modern Earth-like composition and thus it might deviate from that produced by a CO 2 -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 TOA ) is always 10 −5 bar.Rayleigh scattering cross-sections for CO The adopted vertical atmospheric pressuretemperature (PT) profile consists of a convective troposphere and an isothermal stratosphere.This vertical PT profile considers the condensation of H 2 O (if any) in the lower troposphere, the condensation of CO 2 (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 α 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 • , 45 • , 60 • , 70 • , 80 • , 85 • , 88 • } 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 CO 2 .In that case, we substituted the surface pressure value taken from the pressure grid with the CO 2 condensation pressure at that temperature.As an example, the OLR and TOA albedo in the (200 K, 2.0 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 ≤ 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 CO 2 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 CO 2 -dominated atmosphere must be thin due to condensation.

Atmospheric chemical composition
The atmospheres tested in this work are all CO 2dominated.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 N 2 and some modelers (e.g.Ramirez et al. 2014) scaled this abundance in their models.Anyway, our choice to not include N 2 has only negligible effects on the radiative properties of the atmosphere, as it is shown in Ramirez (2017).Fig. 1, panel (b), confirms this result: substituting 5% of CO 2 with N 2 in our model would reduce ⟨T 0 ⟩ by ∼2 K at maximum, and by a much smaller quantity in most other cases.
We also included a variable amount of H 2 O, 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 CH 4 , f CH4 : 0, 0.1% and 1% by volume.There are two reasons behind these specific choices.First, low amounts of CH 4 prevent the formation of organic hazes via photochemical reactions in the upper atmosphere.These hazes form efficiently when the CH 4 /CO 2 ratio increases above ∼ 0.1 (Trainer et al. 2006 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).

2009
).Second, low amounts of CH 4 are possibly easier to explain in light of the known geochemical features of Mars.In general, the amount of CH 4 in the Early Mars atmosphere is understood to be primarily controlled by the balance between serpentinization1 of the olivine crust (Chassefière et al. 2013) and photodissociation (Wordsworth et al. 2017) in the atmosphere.
The photodissociation rate is limited by the availability of photons with λ < 160 nm, most of which are Lyman-α photons, whose flux at the Noachian-Hesperian boundary was about 3×10 15 m −2 s −1 (Wordsworth et al. 2017).
In few bars-level CO 2 dominated atmospheres, this limit is reached for f CH4 ∼ 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% CH 4 concentration would be 3×10 15 m −2 s −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 −2 s −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 es-timate is extremely simplified, it allows to check the plausibility of our choices.

RESULTS
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 ε, the eccentricity e, the land albedo α l and the cloud fraction f c .We also explored the interplay between obliquity, eccentricity and the argument of the perihelion ω, 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 ∼10 4 .
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 seasonallatitudinal 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 max ) and lowest (T min ) temperature in the T s map.If T 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 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 max occurs.In particular, we considered two different threshold temperatures T 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 min , we check the run for CO 2 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 (∼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 CO 2 at surface are hatched.Notice that T max and T 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.In Appendix B, we show and discuss the latitudinal OLR, TOA albedo and diffusion profiles for a small number of selected cases.

Dependence on obliquity
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 • ) 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 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 max for P ∼ 0.4 bar and ε = 60 • .
Both the transition from obliquity-dominated to pressure-dominated T 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 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 max .Finally, above 1.33 bar, T 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 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 (∼0.05 2 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 ε = 5 • are stable, and their T max is equal to 2 For z ∼ 60 • .
236.9 and 237.9K respectively.In the A0F0 model the maximum greenhouse limit is reached at around 2 bar (see Fig. 1, panel b, red line), which causes the disappearance of the conditions conductive to seasonal thaws at higher pressures.

Dependence on eccentricity
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 ∼ 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, 0.00 0.03 0.06 0.09 0.12 0.15 0.18  with the effect being strongest when ω = 90 • or 270 • and weakest when ω = 0 • or 180 •3 .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 ω = 286.5 • .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 max for higher eccentricities, but the limiting pressure required for thawing, under the more conservative condition, changes only by ∼ 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 3 If ω = 90 • , 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 ω = 270 • , but considering the northern hemisphere summer solstice.Instead, when ω = 0 • or 180 • , the planet is at equinox when it passes through the periastron.
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 ε = 25 • there are no non-condensing cases (see the previous Section).The highest T 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.

Dependence on land albedo
The third parameter tested is the land albedo.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 0.11 0.17 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 CO 2 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 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.

Dependence on the cloud fraction
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 selfconsistent model for clouds, we wanted to explore the implications of the simple parameterization presented in Sect.2.1.4by varying the cloud cover fraction be-  tween 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.

H2O clouds
Concerning the cases for which we adopted Earth-like H 2 O clouds, any increase in cloud fraction causes T 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 be-tween 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 H 2 O 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 max value is found at the maximum tested pressure, for a 0% coverage.In the BHGH model, the minimum T max is encountered at 0.1 bar (244 K) as in the previous explorations.However, T 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 ∼150 K to T 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 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.

CO2 clouds
The behavior of T max as a function of cloud coverage in the model with CO 2 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 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 ∼ 50 K at 100% coverage.The effect of clouds on T max is thus different from the one on the average planetary temperature reported in Fig. 2, panel (b).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 max is not sufficient for allowing surface thaws.The maximum value of T 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 H 2 O 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.

Combined effects of the obliquity, eccentricity and argument of perihelion
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 ε 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 ε, 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 • , 270 • and 90 • .
The results are reported in Fig. 7. First of all, this test confirms that obliquity is the main driver of seasonality, since T max variations are more pronounced on the horizontal, rather than on the vertical axis, except for very low (< 10 • ) obliquities.In general, the maximum T max value is found for the highest tested e and ε values, even though its value is different when ω changes.In particular, for model B0GH, the T max maxima are 276, 292 and 283 K for ω = 180 • , 270 • and 90 • , 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 ω changes the regions of the e − ε 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 ω = 90 • make it easier to cross the condensation curve of CO 2 during the colder South Pole winters.The opposite is true when ω = 270 • : 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 ω: in the 180 • and 270 • cases, increasing e shrinks it, while in the 90 • case widens it, preventing stability at lower ε.
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 − ε 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 max are more pronounced.On the other hand, for P 0 = 1.33 bar, most combinations allow for stable atmospheres (again, ω = 90 • produces a higher number of runs with CO 2 surface condensation than ω = 270 • ) and roughly ∼ 25 − 35% of them allow for seasonal thaws.

Duration of the thaw seasons
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 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 pressureeccentricity 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.panel a 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 ∼ 0.15 of martian year (≳ 100 d).In other words, when thawing can occur, then it is not limited to negligibly short (≲ 0.1) fractions4 .For example, a typical fraction value along the border between the non-deglaciating and the deglaciating regions for a 273 K threshold temperature (Fig. 8, panels a, c and e) is 0.3.This happens because T 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.

DISCUSSION
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 H 2 or CH 4 ) are not strictly needed, since when obliquity is sufficiently high seasonal thaws can happen even if the atmosphere is made only of CO 2 with a minor contribution from H 2 O (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 H 2 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 H 2 or CH 4 ), 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 long thaws even at 0.1 bar (panels a and b).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 ∼0.5 bar has been sequestered at surface.By combining their results with the upper limit on CO 2 escape (0.15 bar, Lammer et al. 2013), a reasonable maximum pressure of ∼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.

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 • 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 • S and a second one at 5 • 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 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 conditions compatible with liquid water runoffs at low southern latitudes, even for the lower T 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 ≳ 2 bar of pressure to let pure water ice to melt in the 0 • S -45 • 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 ∼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.

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 H 2 O (RH = 0%, 20% and 40%) and CH 4 (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 CO 2 condensation at surface.In total, we run ∼10 4 cases (i.e.different combinations of input parameters).
Our main findings can be summarized as follows: 1.If both the planet obliquity and datum pressure are sufficiently high (ε ≥ 50 • and P 0 ≥ 2.0 bar) then no other additional greenhouse gases other than CO 2 and desert-like levels of H 2 O are needed to produce seasonal thaws in the northern hemisphere.These conditions also guarantee that no CO 2 condensation happens at the surface.How-ever, the same models fail to produce seasonal melting in the southern highlands.
2. Adding a 0.1% of CH 4 lowers the requirements for northern thaws to 1.33 bar (at ε = 45 • ) and simultaneously allows for melting in the 60 • S -90 We found that the minimum cloud coverage to prevent the atmospheric collapse is 0.5 if there is no CH 4 , 0.2 when f CH4 is at 0.1% and are not needed when f CH4 is at 1%. 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. 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. 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 CO 2dominated atmospheres with high (≥ 0.8) CO 2 -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 CO 2 and H 2 O to allow for an hydrological cycle that involves the southern highlands.However, we show that depending on the combinations of orbital (ε, e and ω) and planetary (land albedo, cloud coverage) parameters, this contribution can be significantly lower than usually expected (∼2-3% of H 2 in a 2 bar CO 2 atmosphere, Ramirez 2017; Wordsworth et al. 2017) 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 advan-tages of using seasonal-latitudinal, rather than singlecolumn, 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 and Ramses Ramirez for the useful suggestions.The Authors also thanks the anonymous referee for the careful reading of the manuscript and the valuable comments.We Both the vertical pressure-temperature profile and the altimetric adjustment have been calculated using a modified version of the moist pseudo-adiabatic lapse rate Γ provided by the American Meteorological Society Glossary5 : 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 ϕ is the ratio between the molecular masses of the condensible and the where T i is the temperature of the ith layer and R m is the gas constant of the atmospheric mix (dry plus moist).Fig. 10, panel (a), shows four examples of PT profiles calculated using this procedure.There is a small offset with respect to the CO 2 condensation curve due to the fact that these profiles refer to the case with f 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 CO 2 compressibility.However, that formulation approximate the pressure of the condensible to zero, which is a reasonable assumption only for relatively cold (T s ≤ 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.

B. DIFFUSION, OLR AND ALBEDO PROFILES OF SELECTED CASES
As discussed in the paper, seasonal thaws can occur on the surface of Early Mars because of the imperfect heat redistribution between latitudinal bands.Since the diffusion parameter D in Eq. 1 has a complex dependence on both input parameters and simulation variables, a more relevant quantity to show is the net flux of heat N along the latitudinal direction.We show the yearly averaged value of N , ⟨N ⟩, for a very small number of selected cases in Fig. 10, panel (b).A positive (negative) N means a poleward flux in the northern (southern) hemisphere.We remind that, in all cases, diffusion is calculated before applying the altimetric correction or, in other words, under the hypothesis that it is dominated by processes happening at the same altitude throughout the entire planet.We start considering the cases with ε = 0 • , represented by yellow, green and blue curves, for which N is nearly time-independent (apart for the effect of obliquity) and thus ⟨N ⟩ ∼ N .Solid curves refer to a configuration without oceans (the BLG0 model) and, as such, are symmetric with respect to the equator.In this configuration the role of P 0 can be isolated: it increases by ∼ 50 − 60% at all latitudes P 0 is increased from 0.1 to 1 bar, while it seems constant when P 0 increases again to 5.62 bar.Analyzing the dependence of the peak transport ⟨N ⟩ max on P 0 (not shown) we notice that it reaches a peak at 2.37 bar and then slightly decreases.This is most probably due to the weakening of the dependence between the albedo and z at high pressures, which reduces the difference in the absorbed radiation between the equator and poles.The addition of the ocean (dashed curves) breaks the north-south symmetry and has the effect of enhancing (when glaciated, yellow and green) or suppressing (when ice-free, blue) the latitudinal heat transfer in the northern hemisphere.If we instead consider the cases with ε = 60 • (the black curves), we find quite a different picture.First of all, ⟨N ⟩ has the opposite sign in both hemispheres, which means that the heat flux is equatorward for at least part of the year (i.e. during summer), and the overall effect of this inverted circulation is stronger when averaged on the entire orbit.Checking the yearly average temperatures ⟨T 0 ⟩ (not shown) we find that polar ⟨T 0 ⟩ are slightly higher than the equatorial one (239.4K vs 233.7 K in the BLG0 model).The higher thermal capacity of oceans with respect to the martian surface makes the northern hemisphere a strong contributor of heat to the rest of the planet.
We complete the picture by discussing also the OLR and the TOA albedo profiles for the same cases described above.In Fig. 10, panel (c), we show the yearly averaged value of the OLR, ⟨I⟩.Starting again from cases with ε = 0 • , we underline the flattening of the ⟨I⟩ curve when pressure increases, which corresponds to a flattening of the ⟨T 0 ⟩ curve.The presence of an ocean breaks the north-south symmetry, in particular lowering the northern longwave emission when it is glaciated (in the 0.1 and 1 bar cases) or slightly increasing it (in the 5.62 bar case).The cases with ε = 60 • show instead an inverted curve, with higher average emissions at the poles, especially when an ocean is present.
Finally, in Fig. 10, panel (d), we show the yearly averaged value of the TOA albedo, ⟨A⟩.In cases with ε = 0 • , the curve is always bowl-shaped due to the fact that, at high latitudes, the z angle is larger, light path is longer and thus Rayleigh scattering is more efficient.The surface albedo is also larger due to the slanted insolation.The presence of ice in models with oceans is underlined by the increase of ⟨A⟩ in the northern hemisphere and at the latitudes corresponding to the Hellas basin (cases with P 0 = 0.1 and 1.0 bar), while the darker deglaciated ocean in the P 0 = 5.62 bar slightly decreases the albedo in the same bands.The cases with ε = 0 • have a substantially different behavior.When an ocean is not present, ⟨A⟩ reaches a maximum at φ = 30 • and then decreases.This is due to the fact that A has an effect on the climate (and thus it is calculated) only when a given latitudinal band is illuminated, and for non-zero obliquities there is always a period of "polar night" in each hemisphere during which the Sun remains below the horizon the entire day.This means that ⟨A⟩ is skewed towards summertime A values, which are smaller due to the lower average z.The irregular behavior of the curve near the poles is due to the discrete nature of the calculation.Again, the presence of an ocean increases the albedo in the northern hemisphere, which is somewhat unexpected considering that the curves in panels (b) and (c) seem to suggest it is mostly deglaciated.The apparent discrepancy comes again from the fact that ⟨A⟩ is calculated only considering the part of the year in which the poles are illuminated by the Sun: due to the high albedo of ices, the northern ocean remains glaciated for ∼60% of the illuminated season, which increases ⟨A⟩.Once the ice is melted, the ocean can efficiently accumulate heat due to the low albedo, which is then released during the polar night season.

Figure 1 .
Figure 1.Panel (a): 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 dashed and dotted lines: isolines corresponding to a 550-m GEL ocean and a 150-m GEL ocean, respectively.Panel (b):The yearly global average planetary temperature at the altimetric datum ⟨T 0⟩ 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 Table3.Each marker represents an increase of 33% of the pressure with respect to the preceding one.Grey dashed and dotted lines mark the 273 K and the 252 K level, respectively.

Figure 2 .
Figure 2. Panel (a): 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).Panel (b): 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).
2 and CH 4 (Sneep & Ubachs 2005) and H 2 O (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).

Figure 3 .
Figure3.The maximum surface temperature Tmax 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.

Figure 4 .
Figure 4.The maximum surface temperature Tmax as a function of the surface pressure and eccentricity.From top left moving clockwise: maps for the B0GH, BLGH, BHGH, C0GN models.

Figure 5 .
Figure 5.The maximum surface temperature Tmax as a function of the surface pressure and the land albedo.From top left moving clockwise: maps for the B0GH, BLGH, BHGH, C0GN models.

Figure 6 .
Figure 6.The maximum surface temperature Tmax 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.

Figure 8 .
Figure 8.The fraction of P 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 ω = 270 • .Hatched regions refer to combinations of parameters that causes the CO2 in the atmosphere to condense at surface.

Figure 9 .
Figure9.The latitude interval in which the temperature at h = 0 m is higher than T thr for at least one time step, as a function of the datum pressure P0.Results refer to the martian southern hemisphere.Top row: T thr = 252 K. Bottom row: T thr = 273 K. From left to right: B0GH, BLGH, BHGH and C0GN models.Red, green and blue shaded areas refer to cases with ε = 30 • , 45 • and 60 • , respectively.The vertical dashed lines mark the pressure below which surface atmospheric condensation occurs for each of the shown combinations of pressure and ε.
acknowledge support by the Italian Space Agency with the Life in Space (ASI N. 2019-3-U.0)and ASTERIA (ASI N. 2023-5-U.0)projects, 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

Table 2 .
Parameters Biasiotti et al. (2022)e clouds (second column), Earth-like clouds (third column) and adopted in the ESTM reference paperBiasiotti et al. (2022)(fourth column).Values for CRE and t in the fourth column are expressed as an interval since they are temperature-dependent.

Table 3 .
; DeWitt et al.A summary of the models studied in this paper.
• S region.Increasing the pressure to 3.5 bars or above melts also the region near the equator.If f CH4 is instead set to 1%, both northern and southern thaws are possible independently of the obliquity when P 0 ≥ 1.33 bar.3.Increasing the eccentricity allows for thaws andstabilizes the atmosphere when combined with ε ≥ 50 • and ω ∼ 270 • (i.e. when Mars is at perihelion during the northern hemisphere's summer solstice).The effect of eccentricity and obliquity alone can change T 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.Changing the land albedo has some impact in the models without CH 4 , while it is less important when CH 4 is present.If f CH4 is set to 1%, both northern and southern melting are possible at P 0 ≥ 2 bars independently of the land albedo.