# Ancient Mars Climate With a Polar Ocean and Ice Sheet Dynamics

Frédéric Schmidt, Michael Way, Aurélien Quiquet, Igor Aleinov, Christophe Dumas

► To cite this version:

Frédéric Schmidt, Michael Way, Aurélien Quiquet, Igor Aleinov, Christophe Dumas. Ancient Mars Climate With a Polar Ocean and Ice Sheet Dynamics. Journal of Geophysical Research. Planets, 2025, 130 (10), 10.1029/2025JE008970. [hal-05296531](https://hal.science/hal-05296531v1)

HAL Id: hal-05296531

https://hal.science/hal-05296531v1

Submitted on 4 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Planets

## RESEARCH ARTICLE

10.1029/2025JE008970

**Key Points:**

* We asynchronously couple a full complexity 3D Global Climate Model with an ice sheet model for Mars at 3 Ga

* The global equilibrium ice sheet mass is evaluated at ~350 m global equivalent layer located in the highest altitude regions

* The timescale to reach a global water cycle equilibrium is ~100 ky, dominated by the ice sheet timescale

**Correspondence to:**

F. Schmidt,
frederic.schmidt@universite-paris-saclay.fr

**Citation:**

Schmidt, F., Way, M. J., Quiquet, A., Aleinov, I., & Dumas, C. (2025). Ancient Mars climate with a polar ocean and ice sheet dynamics. *Journal of Geophysical Research: Planets*, 130, e2025JE008970. https://doi.org/10.1029/2025JE008970

Received 17 FEB 2025

Accepted 9 SEP 2025

**Author Contributions:**

**Conceptualization:** Frédéric Schmidt, Michael J. Way, Aurélien Quiquet, Igor Aleinov, Christophe Dumas

**Formal analysis:** Frédéric Schmidt, Michael J. Way, Aurélien Quiquet, Igor Aleinov, Christophe Dumas

**Funding acquisition:** Michael J. Way, Aurélien Quiquet

**Investigation:** Frédéric Schmidt, Michael J. Way, Aurélien Quiquet, Igor Aleinov, Christophe Dumas

**Methodology:** Frédéric Schmidt, Michael J. Way, Aurélien Quiquet, Igor Aleinov, Christophe Dumas

**Resources:** Michael J. Way, Aurélien Quiquet

**Software:** Frédéric Schmidt, Michael J. Way, Aurélien Quiquet

**Supervision:** Frédéric Schmidt

**Validation:** Frédéric Schmidt, Michael J. Way, Aurélien Quiquet, Igor Aleinov, Christophe Dumas

© 2025. The Author(s).

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

# Ancient Mars Climate With a Polar Ocean and Ice Sheet Dynamics

Frédéric Schmidt¹,² ID, Michael J. Way³,⁴ ID, Aurélien Quiquet⁵ ID, Igor Aleinov³,⁶ ID, and Christophe Dumas⁵ ID

¹Université Paris-Saclay, CNRS, GEOPS, Orsay, France, ²Institut Universitaire de France, Paris, France, ³NASA Goddard Institute for Space Studies, New York, NY, USA, ⁴Theoretical Astrophysics, Department of Physics and Astronomy, Uppsala University, Uppsala, Sweden, ⁵Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France, ⁶Center for Climate Systems Research, Columbia University, New York, NY, USA

**Abstract** In this article, we study the conditions required to maintain a stable ocean on Mars 3 Ga using a new suite of simulations. These simulations couple a 3D Global Climate Model with ocean dynamics and ice sheet flow. The model includes the main processes of the atmosphere/hydrosphere/cryosphere to investigate Mars' ancient climate. The results show that the total water content required to maintain an ocean is ~700 m, global equivalent layer, half in the ocean, half in the ice sheet. This number seems plausible if a significant amount of water has been absorbed by the ground. This could be in the form of mineral alteration, or in a deep porous reservoir. In addition, the results show that the equilibrium mass flux from the ice sheet adjacent to the northern ocean is ~10¹⁵ kg/y with a very low sliding velocity (1 m/y), except for few warm regions in the lowest altitudes that could reach up to 300 m/y. Finally, the global atmosphere/hydrosphere/cryosphere equilibrium should be reached in a few 100 ky. This indicates that the ocean will have a stabilizing feedback on timescales longer than this. An extensive sensitivity study of the ice sheet was performed. This included the effects of a geothermal heat flux, viscosity and basal drag. Finally, we studied the possible effects of planetary obliquity and a reduced ocean extent.

**Plain Language Summary** The ancient climate of Mars remains mysterious since different approaches give inconsistent results. On the one hand, geomorphological clues indicate that an ocean should have been stable 3 Ga, due to indirect tsunami morphologies. On the other hand, some climate model studies indicate that water is more stable in the form of ice on the highlands. Here we propose a new suite of numerical simulations, coupling atmospheric climate with ocean circulation and ice sheet models to investigate the ancient Martian climate. The results quantify the total water required for such a cold, but wet climate: in essence we find as much water in the form of ice sheets as in the ocean. This number is plausible if a significant amount of water has disappeared into the ground, for instance in the form of mineral alteration, or in a deep porous reservoir.

## 1. Introduction

The possibility of an ocean on ancient Mars has been proposed since the 1990s (e.g., Baker et al., 1991) with a lot of controversy. Several reviews of the hypothesis of an ancient northern ocean on Mars have been written in recent years (e.g., Dickeson & Davis, 2020; Kite & Conway, 2024) including insights from the Zhurong rover (Wang & Huang, 2024). These studies point to the episodic presence of an ocean in the early Hesperian to the early Amazonian (about 3.6–2.5 billions years ago).

This hypothesis has been relaunched by the discovery of potential tsunami deposits (Costard et al., 2017; Rodriguez et al., 2016) with at least two impact events. In addition, the Lomonosov crater morphology is coherent with an impact in shallow water, at the very same age of the tsunami deposits (Costard et al., 2019), around 3 Ga, that could be the source of the tsunami. A detailed geological analysis identified similarities between Olympus Mons and other edifices with oceanic islands on Earth (Hildenbrand et al., 2023). Following this theory, the present-day high altitude of the shoreline boundary on volcanoes is due to subsequent isostatic rebound.

In previous work, the long term stability of an ocean in a cold and wet Martian climate seemed impossible in three dimensional General Circulation Models (3D-GCM) (e.g., Forget et al., 2013; Kite et al., 2021; Turbet & Forget, 2019; Turbet et al., 2017; Wordsworth et al., 2013). These studies found that the water tends to accumulate
**Visualization:** Frédéric Schmidt, Aurélien Quiquet

**Writing – original draft:** Frédéric Schmidt

**Writing – review & editing:** Frédéric Schmidt, Michael J. Way, Aurélien Quiquet, Igor Aleinov, Christophe Dumas

in the form of ice in the southern highlands. This view changes when ice sheet processes and ocean circulation are included (Schmidt et al., 2022, hereafter Schmidt2022). In this study, a fully equilibrated water cycle has been proposed with a simplistic ice sheet model. It allows one to estimate the mass balance between the reservoirs but not the actual ice sheet mass. The present article aims at improving those shortcomings.

Previously an ice sheet model was proposed for the Late Noachian period, without 3D-GCM coupling (Fastook & Head, 2015), based on the Wordsworth et al. (2013) model. This previous work focused on the region surrounding Hellas while Schmidt2022 indicated that the major ice accumulation was near the Tharsis region. The main difference is that Wordsworth et al. (2013) did not remove the Polar Layered Deposits, but Schmidt2022 did.

A pioneering work with a coupled ice sheet and 3D-GCM was proposed with a 2 bar atm in the Noachian period (Kamada et al., 2022). The results show that starting with a 500 m Global Equivalent Layer (GEL) of water, all water should be transferred to ice sheets in $10^5$ Mars Years. Here instead, we would like to estimate the total water mass required to keep both the ocean and the ice sheet in equilibrium.

We follow up on Schmidt2022 with improved surface modeling. The aims of this article are to estimate the coupled ice sheet/climate processes, including the following strong bi-directional coupling:

* ice sheet creation depends on snow accumulation, which is driven by the climate.

* the ice sheet transports ice that covers the land and significantly modifies the surface albedo, but also the topography.

* the climate is affected by the surface energy balance, governed by the surface albedo. The topography has a strong influence on the winds. Both effects impact snow accumulation.

To our knowledge, this is the first instance where a 3D-GCM is coupled to a detailed ice sheet model on Mars to estimate the equilibrium water cycle. We use asynchronous coupling, with alternative equilibrium climate and equilibrium ice sheet modeling. This is the standard method used in early Earth to study ice sheet dynamics (Contoux et al., 2015).

This article aims to study the potential distribution of ice, including ice sheet flow, on Mars at 3 Ga and at estimating the water cycle at this time. The results of this numerical study can be extended to Earth-like climate conditions on Mars (1 bar, with stable surface liquid water), that are also foreseen earlier than 3 Ga. An extension to the Noachian is more challenging because the topography may have been completely different due to true polar wander and the formation of the Tharsis bulge.

## 2. Models

The work is here based on 2 models: ROCKE-3D, a three-dimensional General Circulation Model that is used for terrestrial planet climate studies (Way et al., 2017) and GRISLI (Quiquet et al., 2018) a 3D thermo-mechanical ice sheet model.

### 2.1. ROCKE-3D

We use ROCKE-3D version Planet 1.0 (Way et al., 2017) which includes 40 atmospheric layers, a 4 × 5 degree horizontal grid, 6 soil layers, and up to 13 ocean layers. The ocean used is a fully coupled dynamical ocean as described in Way et al. (2017). The atmosphere is 1 bar CO₂ dominated with 22% H₂ for the ocean to be stable. This is a larger amount of H₂ than used in Schmidt2022. We continue to use the same SOCRATES radiative transfer scheme as in Schmidt2022, Edwards and Slingo (1996), Manners et al. (2012), and Amundsen et al. (2016). In Schmidt2022 we used ab-initio calculations of CO₂–H₂ absorption from Wordsworth et al. (2017) where for a 1 bar CO₂ dominated atmosphere, only 10% H₂ was used. In this publication we use the more accurate laboratory measurements of Turbet et al. (2020). We did a series of experiments to determine the amount of H₂ needed with the new tables to obtain a warming that corresponds closely to that using the older ab-initio based tables. We found that 22% of H₂ using the new Turbet et al. (2020) tables is equivalent to 10% using the old Wordsworth et al. (2017) tables in ROCKE-3D with a 1 bar CO₂ dominated atmosphere. To that end we used the following longwave and shortwave SOCRATES spectral files. SOCRATES spectral files can be downloaded from our data portal: https://portal.nccs.nasa.gov/GISS_modelE/ROCKE-3D/spectral_files/. The solar spectrum at 2.9 Ga as defined in the setups for all runs herein sp_lw_16_dsa_pn_hit16, sp_sw_46_dsa_pn_hit16, and solar_spec = sun_2.9gya. This value of 22% of H₂ appears quite large but largely depends on the total
atmospheric pressure. For a 1 bar atm, this value is in agreement with previous studies (e.g. Kamada et al., 2021, Figure 3), while in a 2 bar CO₂ dominated atmosphere the amount of H₂ would be less.

Different scenarios have been proposed to explain a high H₂ content in the early Martian atmosphere. This topic, is still an open area of research and beyond the scope of our article. Here we shortly review and discuss the most relevant scenarios proposed in the literature that may provide high H₂ partial pressure. The first one proposed that CH₄ is trapped in clathrate that can be release in a large event (Chassefiere et al., 2016). In such scenario, a massive release could lead up to 10%–20% H₂ in a 1–2 bar CO₂ dominated atmosphere. A second scenario is the serpentinization of ultramafic rocks (Tutolo & Leong, 2025). This scenario is based on the presence of Fe-rich olivine detected on Mars that is more efficient than Mg-rich olivine to produce H₂ through water alteration. The next scenario is the formation of authigenic magnetite in mudstone under anoxic conditions (Tosca et al., 2018). This scenario is supported by laboratory experiment and the observation of magnetite in Gale Crater. A detailed geochemical model suggests that episodic warm periods of high H₂ content could have lasted ~10⁵ years (Adams et al., 2025) due to crustal alteration. Another scenario is the impact origin as originally proposed in Haberle et al. (2019). They proposed that a very large impactor of diameter >100 km could lead to surface temperatures above freezing for many thousands of years. Hence for a 1 bar CO₂ dominated atmosphere, there are several ways that separately or in combinations could reach 22% H₂ partial pressure. Nevertheless, the residence time for H₂ in such an atmosphere at 3 Ga is still largely unconstrained and is a matter for future investigations.

The paleo-pressure on Mars is not well constrained, especially the pressure 3 Ga (Kite, 2019) and the limited data coming from very different approaches are contradictory. The most recent review proposes multiple climate transitions (Kite & Conway, 2024), where the conditions provided herein are plausible, although as an end-member. The purpose of this article is not to explicitly demonstrate these conditions, but more to assess the implications of such conditions on the climate and the geologic record. An important aspect is the lack of absolute age. The geologic record is only possible using crater counting methods. These have large errors and without radiochronological dating, our progress in that respect is limited. We used the 3 Ga age, as in Schmidt2022 and Costard et al. (2017).

We use an insolation of 452.8 Wm⁻² appropriate for Mars at 3 Ga, approximately 79% of its current value. The obliquity is fixed at 40° (Laskar et al., 2004) and the topography is the same as that used in Schmidt2022 with the exception of changes due to ice sheet accumulation. The physics timestep of the model was set to 1,850 s, while the atmospheric dynamical timestep was set to 230 s and the ocean was set to 450 s. The radiation was called every 5 physics time steps. The model is run until the net radiative balance at the top of the atmosphere between incoming short wave and outgoing long wave is <±0.2 Wm⁻². However, we also desire hydrological equilibrium on the surface, and for this reason we plot a number of additional hydrological quantities as the model moves forward in time to check for stability. These include the mass of water in lakes (quantity defined as MWL in the model output diagnostics), the sea surface height, the ocean and lake ice fraction (OiceFr), the amount of snow deposited on the surface (zsnow or snowdp), and the amount of ground water and ice in the 6 soil layers (GWTR and GICE). If these appear stable/constant in time we consider them to be in equilibrium. For the simulations in this publication we typically find that 200 orbits is sufficient time to reach both types of equilibrium. For simplicity, in this article the time is counted as Earth year (y), except when specified as Martian year. This corresponds to roughly 4 weeks of wall clock time when using 44 cores of 2.4 GHz Intel Xeon Skylake or 2.4 GHz Intel Xeon Cascade CPUs.

### 2.1.1. ROCKE-3D Initialization

As the very first step (before coupling to the icesheet model), we run ROCKE-3D with our 1 bar CO₂ with 10% H₂ simulation from Schmidt2022 with the updated radiation tables described in the previous section. In this new setup our atmosphere is 1 bar CO₂ with 22% H₂. We use a modern Mars topography but with the polar caps removed, but all other aspects of the model setup (i.e., ocean coverage, soil type, albedo, etc.) are the same as those in Schmidt2022. Once this model (noted R for ROCKE-3D) has reached radiative equilibrium (200 orbits) the model results are passed to GRISLI (see next Section).

In the case of a run with an input from GRISLI (for instance RG; G for GRISLI), we input the new topography and albedo maps consistent with the presence of the ice sheet computations from GRISLI.
<table><caption>**Table 1**
GRISLI Model Parameters Used in This Study</caption><thead><tr><th>Variable</th><th>Identifier name</th><th>Value</th></tr></thead><tbody><tr><td>Planet characteristics</td><td></td><td></td></tr><tr><td>Gravitational acceleration</td><td>g</td><td>3.711 m s<sup>-2</sup></td></tr><tr><td>Geothermal heat flux</td><td>ghf</td><td>45 mW m<sup>-2</sup></td></tr><tr><td>Deformation and sliding</td><td></td><td></td></tr><tr><td>Glen flow law enhancement factor</td><td>sf</td><td>5</td></tr><tr><td>Basal drag coefficient</td><td>cf</td><td>0.15 &times; 10<sup>-3</sup> Pa s m<sup>-1</sup></td></tr><tr><td>Surface mass balance</td><td></td><td></td></tr><tr><td>PDD melt coefficient, ice</td><td>C<sub>ice</sub></td><td>0.008 m&deg;C<sup>-1</sup></td></tr><tr><td>PDD melt coefficient, snow</td><td>C<sub>snow</sub></td><td>0.003 m&deg;C<sup>-1</sup></td></tr><tr><td>Temperature standard deviation</td><td>&sigma;</td><td>5&deg;C</td></tr><tr><td>Vertical lapse rate</td><td>&gamma;</td><td>0.004&deg;C m<sup>-1</sup></td></tr></tbody></table>

## 2.2. GRISLI

To simulate ice sheet evolution we use the 3D thermo-mechanically coupled ice sheet model GRISLI (Quiquet et al., 2018). The model uses the standard Glen's flow law constitutive equation in which it is assumed that the ice sheet is a viscous non-Newtonian fluid composed solely of pure incompressible ice. The model uses both the shallow ice and shallow shelf approximations to compute the ice velocity. The shallow ice approximation velocity is driven by vertical shearing and is predominant in cold, slow flowing regions. The shallow shelf approximation neglects vertical deformation and becomes dominant for fast flowing regions, that is, ice streams and ice shelves, where the flow is controlled by lateral shearing. In the model, the two components of the velocities are added, assuming that the basal sliding velocity is given by the shallow shelf approximation component. The vertical profile is given by the shallow shelf approximation component. As in most ice sheet models, a flow enhancement factor (constant and homogeneous, see Table 1) is used to favor deformation due to vertical shear to artificially account for ice anisotropy. For temperate regions, we assume a linear basal drag (see Table 1) reducing the shallow shelf approximation velocity. For cold-based regions, this drag is virtually infinite. The model accounts for glacial isostasy, using an elastic lithosphere—relaxed asthenosphere model (Le Meur & Huybrechts, 1996).

GRISLI has been used extensively to simulate Greenland and Antarctic ice sheet evolution in Earth's past (Crotti et al., 2022; Ladant et al., 2014; Quiquet et al., 2013) and in the future (Le clec'h et al., 2019; Quiquet & Dumas, 2021a, 2021b) as well as for paleo Earth ice sheets over a range of timescales (Benn et al., 2015; Tan et al., 2017). For this work we use relatively standard ice sheet model parameters (listed in Table 1, Quiquet et al. (2021) and van Aalderen et al. (2024)) except that the gravitational acceleration was adapted to Mars. The ice domain is a 40 × 40 km Cartesian grid covering the Tharsis region (236 × 230 GRISLI grid cells in total), projected using a Lambert Azimuthal Equal Area projection (center latitude 5°, center longitude -95°, radius 85°, using Mars 2000 ellipsoid). This domain represents the majority of the snow accumulation over the highest altitude point in Mars. Nevertheless minor accumulation outside this domain, for instance around Hellas and near Elysium, are ignored here but should be investigated further. The main time step is 5 years with a smaller time step for the mass conservation equation. Total run time is 400 ky, which takes approximately 18 hr to complete on a single core of a modern microprocessor.

To compute the surface mass balance (SMB) we use the standard positive degree day (PDD) melt model of Reeh (1991). This has to be fed by monthly means of near-surface air temperature and precipitation, provided by the ROCKE-3D climate model. Assuming a monthly temperature standard deviation, this model computes the annual PDDs that are then linearly linked via melt coefficients to surface ablation (fresh snow first and then ice). The PDD model does not account for sublimation. The sublimation can be important here since the simulated climate is very cold and dry. For this reason, we additionally impose the GCM computed sublimation rates onto the SMB computed by the PDD model. In addition, we filter out the very low accumulation areas by imposing a
slightly negative SMB (–1 mm/yr) where the accumulation rate is lower than 5 mm/yr. Elevation changes induce a temperature change using a constant and homogeneous vertical lapse rate of 0.004° Cm⁻¹ (See Table 1).

### 2.2.1. GRISLI Initialization

The first GRISLI run (noted RG) is forced with the first ROCKE-3D output (noted R). For this first run, GRISLI starts from ice free conditions and is forced by the field from R. The ice sheet model is run for 400 ky to reach an equilibrium with the climatic forcing fields. We used the temperature and the precipitation monthly fields, averaged over the last 10 martian years from R. The same scheme is repeated each time, with the corresponding updated ROCKE-3D output.

### 2.3. Coupled Model

The climate and the ice sheets are intimately coupled. Indeed, the ice snow accumulation from the climate is a source for the ice sheet but in return, the elevation and the extent of ice (and thus albedo) due to flow may also change the climate. The typical timescale for ROCKE-3D to reach equilibrium is 200 Martians years, but the ice sheet requires a modeling timescale around 10,000 years. Indeed, some processes in ice sheets can be very long, for instance isostasy movements require ~3,000 years, temperature diffusion ~10,000 years, etc. It is therefore impossible to compute both at the same time. Instead, we propose the standard scenario: asynchronous coupling. It consists of:

1. R computes the equilibrium climate using ROCKE-3D

2. RG computes the equilibrium ice sheet using GRISLI using input from R

3. RGR computes the equilibrium climate using ROCKE-3D using input from RG

4. RGRG computes the equilibrium ice sheet using GRISLI using input from RGR

5. RGRGR computes the equilibrium climate using ROCKE-3D using input from RGRG

6. RGRGRG computes the equilibrium ice sheet using GRISLI using input from RGRGR

In parallel, we also conduct experiments to study the effect of shrinking the northern ocean extent:

1. RGR15 computes the equilibrium climate using ROCKE-3D from RG with 15% less water volume in the ocean (38% less ocean surface)

2. RGRG15 computes the equilibrium ice sheet using GRISLI from RGR15

3. RGR30 computes the equilibrium climate using ROCKE-3D from RG with 30% less water volume in the ocean (51% less ocean surface)

4. RGRG30 computes the equilibrium ice sheet using GRISLI from RGR30

5. RGR50 computes the equilibrium climate using ROCKE-3D from RG with 50% less water volume in the ocean (67% less ocean surface)

6. RGRG50 computes the equilibrium ice sheet using GRISLI from RGR50

## 3. Results

### 3.1. Climate

Table 2 presents the same integrated results as in Schmidt2022. We used the same concept to separate the planet into 3 domains: the ocean, the Wet lowlands, with annual mean temperature above 0°C, dominated by rain and the Icy highlands with annual mean temperature below 0°C, dominated by snow. These definitions are justified by the fact that the surface temperature is mainly controlled by altitude. The first part of the results table clearly demonstrates that the climate is getting colder when coupling with the ice sheet model due to the albedo feedback. The Icy Highland surface is increasing due to the decreasing altitude of the 0°C isotherm. The corresponding Wet lowland is shrinking by a factor of 1.5. The thickness of sea ice and its fraction of the total ocean surface is also significantly increasing.

Figure 1 shows the main output field of ROCKE-3D after 3 alternative couplings (step RGRGR): rain precipitation, snow and ice fraction at the surface, snowfall, and sea/land surface temperature. The general results are comparable to those in Schmidt2022 except that due to the albedo feedback the snow fraction tends to accumulate more on the Tharsis bulge.
<caption>**Table 2**</caption>

Table of the Main Climatic Parameters From ROCKE-3D

<table><thead><tr><th>Parameter</th><th>R</th><th>RGR</th><th>RGRGR</th><th>RGR15</th><th>RGR30</th><th>RGR50</th></tr></thead><tbody><tr><td>Altitude of the ocean (meter)</td><td>-3,900</td><td>-3,900</td><td>-3,900</td><td>-4,310</td><td>-4,500</td><td>-4,690</td></tr><tr><td>Altitude of the 0°C isotherm (meter)</td><td>-2,070</td><td>-2,990</td><td>-2,920</td><td>-2,890</td><td>-2,940</td><td>-2,950</td></tr><tr><td>Mean temperature (°C)</td><td>-7.2</td><td>-10.5</td><td>-12.0</td><td>-10.2</td><td>-10.4</td><td>-10.4</td></tr><tr><td>Mean ocean temperature (°C)</td><td>7.1</td><td>3.8</td><td>2.8</td><td>4.4</td><td>3.8</td><td>3.0</td></tr><tr><td>Mean land temperature (°C)</td><td>-10.0</td><td>-13.2</td><td>-14.8</td><td>-11.8</td><td>-11.6</td><td>-11.1</td></tr><tr><td>Surface of Icy highland (%)</td><td>62.4</td><td>71.0</td><td>73.8</td><td>70.4</td><td>71.1</td><td>71.1</td></tr><tr><td>Surface of Wet lowland (%)</td><td>21.9</td><td>14.5</td><td>11.9</td><td>21.1</td><td>22.8</td><td>25.1</td></tr><tr><td>Surface of Ocean (%)</td><td>16.2</td><td>16.2</td><td>16.2</td><td>10.0</td><td>7.7</td><td>5.4</td></tr><tr><td>Mean thickness of sea ice (meter)</td><td>0.04</td><td>0.25</td><td>0.48</td><td>0.18</td><td>0.23</td><td>0.31</td></tr><tr><td>Mean fraction of sea ice (%)</td><td>1.3</td><td>4.7</td><td>11.5</td><td>5.6</td><td>7.4</td><td>10.1</td></tr></tbody></table>

Note. The ocean could be either liquid water or ice. All quantities are average over the last 10 years of the simulations.

The second part (right 3 columns) of Table 2 presents the effect of shrinking the ocean by 15%, 30%, and 50% in volume. As expected the altitude and the surface proportion of the ocean decreases. The 0° isotherm and the Icy Highland proportion remains stable, but the Wet lowland proportion is increasing, due to the low altitude of the Northern plain. The global average temperature appears to be stable, whatever the ocean size, but the fraction and thickness of sea ice increases.

Figure 2 shows the effect of a shrinking ocean on the temperature field and the difference to the RGR simulation taken as a reference. As the surface of the Northern ocean decreases, the ocean becomes warmer and the near shoreline land surface also warms up. Figure 3 shows the effect on the rain precipitation that tends to decrease when shrinking the ocean.

## 3.2. Ice Sheet

Figure 4 shows the evolution of the main parameters of the mass balance: gain (SMB including accumulation and ablation) and losses (bottom melt, ice front melting and calving). Gains and losses compensate each other after

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Main fields from the RGRGR simulation for the rain precipitation (upper left), snow and ice fraction at the surface (upper right), snowfall (lower left), and sea/land surface temperature (lower right). Black contour lines represent surface elevation level (–2,000, 0, 2,000, and 8,000 m) while the red contour line is the paleoshoreline (–3,900 m in this case). The dashed white contour line in the upper right panel represents the domain of the GRISLI ice sheet simulation area, centered in the major snow accumulation area around the Tharsis plateau.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Temperature field at RGR step (a) surface temperature of RGR, taken as a reference; (b) Difference with 15% less water in the Northern ocean RGR15—RGRG; (c) Difference with 30% less water in the Northern ocean RGR30—RGRG; (d) Difference with 50% less water in the Northern ocean RGR50—RGRG. As the surface of the Northern ocean decreases, the ocean gets slightly warmer and the near shoreline land surface also warms up.</figcaption></figure>

~100 ky, showing that the model has reached equilibrium. The values at equilibrium are ~10¹⁵ kg/y, in agreement within an order of magnitude of the global fluxes from ROCKE-3D (See Table 4). These values are not strictly equal because their evaluations are not exactly the same. Estimations from ROCKE-3D are at a global scale while GRISLI is only within its smaller simulation grid. The mass flux estimated in the same overlapping region are 1.6 × 10¹⁵ kg/y at RGR and 1.6 × 10¹⁵ kg/y at RGRGR for ROCKE-3D, in better agreement with GRISLI than the estimate 1.1 × 10¹⁵ kg/y at RGR and 1.4 × 10¹⁵ kg/y at RGRGR. These last discrepancies are due to differences in the surface balance scheme between the two models. The typical timescale to equilibrate is around 100 ky for each of the coupling steps (RG, RGRG, RGRGRG), meaning that the global atmosphere/hydrosphere/cryosphere equilibrium should be reached in a few 100 ky, whatever the conditions.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Precipitation by rain for the (a) RGR simulation; (b) with 15% less water in the Northern ocean RGR15; (c) with 30% less water in the Northern ocean RGR30; (d) with 50% less water in the Northern ocean RGR50. As the surface of the Northern ocean decreases, the rain locally close to the shoreline and in the ocean decreases.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Mass flux gain (surface mass balance) and losses (bottom melt, ice front melting and calving) from the 3 GRISLI simulations. The global surface balance mass flux is $0.5 \times 10^{15}$ kg/year at RG, $1.1 \times 10^{15}$ kg/year at RGRG, and $1.4 \times 10^{15}$ kg/year at RGRGRG. These values are in agreement within an order of magnitude of global fluxes from ROCKE-3D (See Table 4). The typical timescale to equilibrate is around 100 ky for each of the coupling step.</figcaption></figure>

Figure 5 shows the evolution of the main parameters: ice sheet extent and ice sheet volume. At the end of the RG step, the total ice sheet volume was $24 \times 10^6$ km³ (152 m GEL liquid water). At the end of the RGRG step, the value reached was $45 \times 10^6$ km³ (283 m GEL liquid water). At the final step of the RGRGRG step, the value reached $59 \times 10^6$ km³ (369 m GEL liquid water). These values are of the same order of magnitude as the ocean volume, meaning that the total water content should be ~700 m GEL to sustain a closed climate cycle. Those values are in agreement with $75 \times 10^6$ km³ found in Kamada et al. (2022). The ice sheet extent reaches $30 \times 10^6$ km².

<table><caption>**Table 3**
Table of the Main Parameters of the Sensitivity Studies on Ice Sheet Model by GRISLI in the RG and the RGRG Stages</caption><thead><tr><th>Exp. name</th><th>Parameter changed</th><th>Value</th></tr></thead><tbody><tr><td>OBL<sub>00</sub></td><td>Obliquity</td><td>0°</td></tr><tr><td>OBL<sub>20</sub></td><td>–</td><td>20°</td></tr><tr><td>OBL<sub>60</sub></td><td>–</td><td>60°</td></tr><tr><td>HIGH_GHF</td><td>Geothermal heat flux</td><td>85 mWm<sup>-2</sup></td></tr><tr><td>HIGH_GRAV</td><td>Gravity acceleration</td><td>9.81 ms<sup>-2</sup></td></tr><tr><td>LOW_SIGMA</td><td>Daily temperature variability</td><td>2°C</td></tr><tr><td>HIGH_SIGMA</td><td>–</td><td>8°C</td></tr><tr><td>LOW_MELT</td><td>Melt PDD factors</td><td>-50%</td></tr><tr><td>HIGH_MELT</td><td>–</td><td>+50%</td></tr><tr><td>VERY-HIGH_MELT</td><td>–</td><td>+150%</td></tr><tr><td>LOW_VISCO</td><td>Flow enhancement factor</td><td>1</td></tr><tr><td>HIGH_VISCO</td><td>–</td><td>10</td></tr><tr><td>HIGH_DRAG</td><td>Basal drag coefficient</td><td>×10</td></tr><tr><td>LOW_DRAG</td><td>–</td><td>/10</td></tr><tr><td>VERY-LOW_DRAG</td><td>–</td><td>/100</td></tr></tbody></table>

At each iteration, both ROCKE-3D and GRISLI reach their own equilibrium but the coupled system should also converge with one equilibrium state. Here we evaluate this coupled convergence state after three iterations of both models because of the extremely heavy computation time. This has required 3 full months of computation time on a parallel super-computer. In order to assess if the final state represents the equilibrium, we could evaluate the difference after each iterations. Figures 4 and 5 both show that the mass flux, the ice sheet surface and ice sheet volume are converging. The last iteration change is less that ~20%. Because this value is not null, one may consider that the mass flux, the ice sheet surface and ice sheet volume are minimum estimates of the final coupled convergence state.

### 3.2.1. Geometry and Dynamics

Figures 6–8 show the ice sheet topography, ice sheet thickness, basal melt and basal velocity computed by GRISLI (step RGRGRG). The ice sheet is up to 4,300 m thick but the flow is relatively limited with 300 m/y at maximum, compared to Earth where it can reach several km/y. This is mainly due to the low gravity. The typical basal melting values are in cm/y, highly correlated to the ice sheet velocity, reaching locally maximum values at 30 cm/y. One interesting point to note is that the ice sheet reaches the ocean in two points in the North East and North West edges, demonstrating that a glacier can flow through the wetland and reach the ocean to potentially produce icebergs. In
<caption>**Table 4**</caption>

Table of the Fluxes (in kg/Earth Year) of the Icy Highland Toward the Other Main Reservoirs From ROCKE-3D

<table><thead><tr><th>Parameter</th><th>R</th><th>RGR</th><th>RGRGR</th><th>RGR15</th><th>RGR30</th><th>RGR50</th></tr></thead><tbody><tr><td>Evaporation</td><td>-1.06e+16</td><td>-1.18e+16</td><td>1.18e+16</td><td>-1.09e+16</td><td>-1.06e+16</td><td>-1.01e+16</td></tr><tr><td>Rain</td><td>5.39e+14</td><td>5.79e+14</td><td>4.68e+14</td><td>4.16e+14</td><td>3.86e+14</td><td>3.40e+14</td></tr><tr><td>Snow</td><td>1.34e+16</td><td>1.42e+16</td><td>1.41e+16</td><td>1.22e+16</td><td>1.17e+16</td><td>1.09e+16</td></tr><tr><td>Runoff to Wet lowland</td><td>-6.54e+14</td><td>-3.65e+14</td><td>-3.16e+14</td><td>-1.89e+14</td><td>-1.89e+14</td><td>-1.91e+14</td></tr><tr><td>Glacier to Ocean</td><td>-2.67e+15</td><td>-2.54e+15</td><td>-2.49e+15</td><td>-1.55e+15</td><td>-1.26e+15</td><td>-9.14e+14</td></tr><tr><td>Balance</td><td>1.00e+00</td><td>-1.00e+00</td><td>0.00e+00</td><td>-5.00e-01</td><td>7.50e-01</td><td>-5.00e-01</td></tr></tbody></table>

*Note.* All quantities are averaged over the last 10 years of the simulations. The mass flux balance in the table is not perfect due to simplification and rounding but the residue is always lower than two orders of magnitude lower than the largest fluxes.

addition, the relative low velocity would prevent the ice from massively eroding the substratum. The isostatic effect can reach up to 800 m.

Figure 9 shows the ice sheet thickness as a function of the ocean level in the step RGRG. The ocean volume has been reduced in the simulation in the step RGR by 0%, 15%, 30% and 50%. The effect on the total thickness is minor but can reach locally 3,000 m near the northern edges of the ice sheet.

### 3.2.2. Sensitivity Studies

The simulated ice sheet geometry is a result of the climatic forcing fields provided by the climate model and parameter choices in the ice sheet model. In addition to our reference experiment we run alternative ice sheet RG and/or RGRG experiments, grouped in four categories:

* Alternative climatic forcing fields. In the reference experiment, the climate model used an obliquity of 40° (Laskar et al., 2004). To assess the importance of this choice in the simulated ice sheet, we also used three climate simulations with 0°, 20°, and 60° obliquity from Schmidt2022. These RG experiments are labeled OBL_00, OBL_20 and OBL_60.

* Alternative planet characteristics. We run GRISLI with a higher geothermal heat flux (HIGH_GHF, 85 mW/m² instead of 45 mW/m² in the reference experiment) and a terrestrial gravitational acceleration (HIGH_-GRAV, 9.81 m/s² instead of 3.711 m/s²). These simulations were run for the RG and RGRG experiments. Changing those characteristics allows one to estimate the effect of the mechanical forcing on the ice sheet.

* Alternative SMB model parameters. We run GRISLI with different snow and ice melt coefficients of the PDD model (LOW_MELT with -50% w.r.t. the reference, HIGH_MELT with +50% and VERY-HIGH_MELT

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Total surface (a, on left) and total volume extent (b, on right) of the ice sheet, computed by GRISLI in the steps RG, RGRG and RGRGRG. The horizontal light blue line represents the total ocean volume 340 m global equivalent layer. The colored lines represents the change in ocean volume experiments in the RGRG. The typical timescale to equilibrate is around 100 ky for each of the coupling steps.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Ice sheet topography computed by GRISLI (step RGRGRG). For this particular simulation, the ice sheet reaches the ocean and could potentially produce icebergs. The altitude here includes the isostasy that can reach up to 800 m. The ice accumulation and ice sheet flow smooth out the surface but the largest peaks are still present in the topography.</figcaption></figure>

with +150%). We also changed the value of the daily temperature variability (LOW rest with 2°C and HIGH rest with 8°C, instead of 5°C). These simulations were run for the RG and RGRG experiments.

* Alternative ice sheet dynamics parameters. Ice flow sensitivity experiments can be easily performed using two tunable parameters: the flow law enhancement factor and the basal drag coefficient. The enhancement factor is a parameter that increases the deformation due to vertical shearing and is mostly influential in slow moving regions. Conversely, the basal drag coefficient is important for temperate-based fast flowing regions. RG and RGRG experiments in this category are labeled LOW_VISCO and HIGH_VISCO (enhancement factor of 1 and 10, respectively, instead of 5), HIGH_DRAG, LOW_DRAG and VERY_LOW_DRAG (basal drag coefficient multiplied by 10 and divided by 10 and 100 with respect to the reference experiment).

Results from Figure 10 show that our results are robust to most of the parametrization, with few exceptions. The effect of obliquity slightly changes the ice sheet extent, but has only minor influence on the total volume. This result may be counter-intuitive but as already discussed in Schmidt2022, the ocean circulation stabilizes the climate and makes it less obliquity dependent. In addition, the relatively modest effects of obliquity in comparison

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Ice sheet thickness computed by GRISLI (step RGRGRG), same as Figure 6. The ice tends to accumulate in the flattest regions near the topographic peaks.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Ice sheet basal melt and basal velocity computed by GRISLI (step RGRGRG), same as Figure 6. The basal melting is highly correlated with the ice sheet velocity, as expected. The maximum basal melting is 30 cm/y but the typical values are in cm/y. The flow is relatively slow in general (1 m/y) but at the lowest altitude can reach (1,000 m/y) due to higher temperatures.</figcaption></figure>

with current 6-mbar martian atmosphere is probably due to more efficient heat transport at higher pressure. The geothermal heat flux and gravity induces a slight reduction in the ice volume, but no significant change in the extent. The melt model parameters have a moderate effect on both volume and extent. Increasing the parameter sigma, which controls the surface temperature variability inside a month, reduces the extent and volume, due to higher surface melting. Lastly, the ice flow parametrization exhibits very low sensitivity on the ice extent, but only slightly on the total ice volume. Overall, in the parameter space tested the results change less than a factor of 2, and with a few 10s of percent for most of the cases.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Ice sheet thickness at the step RGRG (a) with regular ocean level (b) evolution of the thickness with 15% less ocean volume (RGRG15—RGRG) (c) with 30% less ocean volume (RGRG30—RGRG) (d) with 50% less ocean volume (RGRG50—RGRG). The ocean has been reduced in the simulation since the step RGR. The effect of a reducing ocean level on the ice sheet is limited, but the thickness can locally be reduced by 3,000 m near the low altitude edges of the ice sheet, simply by the retreat of the ice border. The maximum isostatic effect is 300 m, near the border of the ice sheet, so these maps in thickness difference can also be interpreted as altitude difference.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Simulated ice sheet volume and ice extent in the sensitivity experiments: climate forcing (blue), planet characteristics (green), melt model parameters (red), ice flow mechanics (purple). The dark colors stand for the RG experiments while the light colors are for the RGRG experiments. Characteristics of the experiments are provided in Table 3. The gray lines stand for the reference experiment: lower one for RG, upper one for RGRG. These sensitivity experiments were performed without filtering for the small accumulation areas as explained in Section 2.2, leading to slight differences with respect to the results shown in Figure 5.</figcaption></figure>

### 3.3. Volatile Cycle

Tables 4–7 present a simplified representation of the water cycle as in Schmidt2022, since non-relevant reservoirs (such as porous soil) are present in ROCKE-3D (Way et al., 2017). Those tables indicate the mass fluxes when alternate coupling with the ice sheet several times and when shrinking the ocean. In general the mass fluxes are equivalent with no significant order of magnitude change. The major changes are:

1. When more alternate coupling with ice sheet, the runoff from Wet lowland to the ocean is reduced by a factor of 3, more than the surface collecting water (1.5 times less surface). The effect of a shrinking ocean seems less important for this parameter.

2. When more alternate coupling with ice sheet, the snow falling inside the ocean is increased by a factor of 12. In the RGRGR state, the snow to the ocean is almost equivalent to the rain.

3. As expected, as the ocean is shrinking, the snow and rain precipitation into the ocean decreases. As well, the ocean evaporation decreases faster than the ocean surface proportion.

4. The glacier flux from the land to the ocean is decreasing only slowly when more alternate coupling with ice sheet but this parameter seems more influenced by the ocean shrinking.

### 4. Comparison With Geological Records

A comparison with the geological record at 3 Ga - in the transition Hesperian/early Amazonian - is difficult because there is a limited number of terrain not covered by subsequent units with this age (Schmidt et al., 2022).

Numerous dendritic valleys have been identified in the Noachian, but only a few during the late Hesperian/early Amazonian. Several valleys at 2.9 Ga have been identified in the North of Alba Patera with a substantial temporal peak in the drainage density (Hynek et al., 2010), located in the dichotomy boundary, in agreement with the
<caption>**Table 5**</caption>

Table of the Fluxes (in kg/Earth Year) of the Wet Lowlands Toward the Other Main Reservoirs From ROCKE-3D

<table><thead><tr><th>Parameter</th><th>R</th><th>RGR</th><th>RGRGR</th><th>RGR15</th><th>RGR30</th><th>RGR50</th></tr></thead><tbody><tr><td>Evaporation</td><td>-3.02e+15</td><td>-1.90e+15</td><td>-1.69e+15</td><td>-2.33e+15</td><td>-2.43e+15</td><td>-2.27e+15</td></tr><tr><td>Rain</td><td>2.44e+15</td><td>1.30e+15</td><td>1.06e+15</td><td>1.45e+15</td><td>1.34e+15</td><td>1.13e+15</td></tr><tr><td>Snow</td><td>1.15e+15</td><td>7.49e+14</td><td>6.93e+14</td><td>1.10e+15</td><td>1.31e+15</td><td>1.26e+15</td></tr><tr><td>Runoff from Icy highland</td><td>6.54e+14</td><td>3.65e+14</td><td>3.16e+14</td><td>1.89e+14</td><td>1.89e+14</td><td>1.91e+14</td></tr><tr><td>Runoff to Ocean</td><td>-1.22e+15</td><td>-5.12e+14</td><td>-4.01e+14</td><td>-4.04e+14</td><td>-3.93e+14</td><td>-3.21e+14</td></tr><tr><td>Balance</td><td>-1.71e+12</td><td>2.08e+12</td><td>-2.07e+13</td><td>9.84e+12</td><td>1.89e+13</td><td>-8.48e+11</td></tr></tbody></table>

*Note.* All quantities are averaged over the last 10 years of the simulations. The mass flux balance in the table is not perfect due to simplification and rounding but the residue is always lower than two orders of magnitude lower than the largest fluxes.

rainfall field predicted by our model (see Figure 1 and Schmidt2022). In addition a relatively high density network has been located in the South of Hellas basin (Bernhardt et al., 2016). In northwestern Arabia Terra and in Deuteronilus Mensae a network of aligned channels are observed to display streamlined islands of Hesperian age of possible marine (Rodriguez et al., 2016) or fluvial origin (Tanaka et al., 2014). The Northern lowlands have been substantially covered by more recent Amazonian material (Tanaka et al., 2014) and cannot be compared.

In the Tharsis area (Figures 6–8), our model predicts a massive ice sheet. The velocity field predicted by the ice sheet is radial, downhill, toward the exterior, as expected. The major ice streamline is toward the North and North-East but it does not follow the present-day valley network. For instance, Kasei valles has been proposed by several authors (Duran & Coulthard, 2020; Hepburn et al., 2020; Lucchitta, 1980, 1995) to be produced by glacial erosion. In our simulation, the low altitude of Kasei prevents the presence of an ice tongue. The same situation appears in Valles Marineris (Gourronc et al., 2014; Mège & Bourgeois, 2011). These discrepancies of geometry could be due to the fact that the topography 3 Ga was significantly different in this region, due to the Tharsis rise. There is no global map of past ice related features but the available maps (Bouquety et al., 2019, 2020; Galofre et al., 2020) shall be complemented in the future. Glacial tongues are in principle possible to model with GRISLI, but this model cannot produce an ice stream running through Kasei valles because the climate is very warm there due to its low altitude. The Tharsis rise ended in the Early Amazonian period (Bouley et al., 2016, 2018), meaning between 2 and 3 Ga depending on the chronology system. We consider in the article the present topography (except the polar caps) as a starting point, but maybe the later stage of Tharsis formation could also be considered. There is a nontrivial effect of load and rotation (Matsuyama & Manga, 2010) and also true-polar wander (Bouley et al., 2016; Perron et al., 2007) that could potentially resolve these discrepancies. The second option would be to reduce the global temperature, while still avoiding the ocean to freeze over. Our condition seems already close to this second option, since we have ~−10°C global temperature and ~5°C ocean temperature (see Table 2) but we could probably lower the global temperature by decreasing the H₂ content. In future work it also would be

<caption>**Table 6**</caption>

Table of the Fluxes (in kg/Earth Year) of the Ocean Toward the Other the Main Reservoirs From ROCKE-3D

<table><thead><tr><th>Parameter</th><th>R</th><th>RGR</th><th>RGRGR</th><th>RGR15</th><th>RGR30</th><th>RGR50</th></tr></thead><tbody><tr><td>Evaporation</td><td>-9.78e+15</td><td>-8.46e+15</td><td>-8.06e+15</td><td>-4.95e+15</td><td>-3.77e+15</td><td>-2.67e+15</td></tr><tr><td>Rain</td><td>5.72e+15</td><td>3.46e+15</td><td>3.03e+15</td><td>1.54e+15</td><td>8.55e+14</td><td>3.97e+14</td></tr><tr><td>Snow</td><td>1.26e+14</td><td>1.41e+15</td><td>1.58e+15</td><td>9.70e+14</td><td>7.68e+14</td><td>5.20e+14</td></tr><tr><td>Glacier</td><td>2.67e+15</td><td>2.54e+15</td><td>2.49e+15</td><td>1.55e+15</td><td>1.26e+15</td><td>9.14e+14</td></tr><tr><td>Runoff</td><td>1.23e+15</td><td>5.66e+14</td><td>4.38e+14</td><td>4.45e+14</td><td>4.08e+14</td><td>3.41e+14</td></tr><tr><td>Sea ice</td><td>-3.28e+10</td><td>-1.00e+11</td><td>-8.28e+10</td><td>-9.73e+10</td><td>-9.81e+10</td><td>-9.71e+10</td></tr><tr><td>Balance</td><td>-4.31e+13</td><td>-4.79e+14</td><td>-5.17e+14</td><td>-4.44e+14</td><td>-4.82e+14</td><td>-4.94e+14</td></tr></tbody></table>

*Note.* All quantities are averaged over the last 10 years of the simulations. The mass flux balance in the table is not perfect due to simplification and rounding but the residue is always lower than one order of magnitude than the largest fluxes. The imbalance is here larger due to neglected fluxes over porous soil.
<caption>**Table 7**</caption>

Table of the Fluxes (in kg/Earth Year) of the Atmosphere Toward the Other the Main Reservoirs From ROCKE-3D

<table><thead><tr><th>Parameter</th><th>R</th><th>RGR</th><th>RGRGR</th><th>RGR15</th><th>RGR30</th><th>RGR50</th></tr></thead><tbody><tr><td>Evaporation from Icy highland</td><td>1.06e+16</td><td>1.18e+16</td><td>1.18e+16</td><td>1.09e+16</td><td>1.06e+16</td><td>1.01e+16</td></tr><tr><td>Rain to Icy highland</td><td>-5.39e+14</td><td>-5.79e+14</td><td>-4.68e+14</td><td>-4.16e+14</td><td>-3.86e+14</td><td>-3.40e+14</td></tr><tr><td>Snow to Icy highland</td><td>-1.34e+16</td><td>-1.42e+16</td><td>-1.41e+16</td><td>-1.22e+16</td><td>-1.17e+16</td><td>-1.09e+16</td></tr><tr><td>Evaporation from Wet lowland</td><td>3.02e+15</td><td>1.90e+15</td><td>1.69e+15</td><td>2.33e+15</td><td>2.43e+15</td><td>2.27e+15</td></tr><tr><td>Rain to Wet lowland</td><td>-2.44e+15</td><td>-1.30e+15</td><td>-1.06e+15</td><td>-1.45e+15</td><td>-1.34e+15</td><td>-1.13e+15</td></tr><tr><td>Snow to Wet lowland</td><td>-1.15e+15</td><td>-7.49e+14</td><td>-6.93e+14</td><td>-1.10e+15</td><td>-1.31e+15</td><td>-1.26e+15</td></tr><tr><td>Evaporation from Ocean</td><td>9.78e+15</td><td>8.46e+15</td><td>8.06e+15</td><td>4.95e+15</td><td>3.77e+15</td><td>2.67e+15</td></tr><tr><td>Rain to Ocean</td><td>-5.72e+15</td><td>-3.46e+15</td><td>-3.03e+15</td><td>-1.54e+15</td><td>-8.55e+14</td><td>-3.97e+14</td></tr><tr><td>Snow to Ocean</td><td>-1.26e+14</td><td>-1.41e+15</td><td>-1.58e+15</td><td>-9.70e+14</td><td>-7.68e+14</td><td>-5.20e+14</td></tr><tr><td>Balance</td><td>5.73e+13</td><td>5.30e+14</td><td>5.74e+14</td><td>4.76e+14</td><td>4.78e+14</td><td>5.15e+14</td></tr></tbody></table>

Note. All quantities are averaged over the last 10 years of the simulations. The mass flux balance in the table is not perfect due to simplification and rounding but the residue is always lower than two orders of magnitude lower than the largest fluxes.

interesting to perform higher spatial resolution modeling and using other possible topographic configurations to estimate the temperature and potential ice flow of such valleys.

In the region of Malea Planum, some deposits could be considered as remnants of an extensive polar ice sheet emplaced during the Hesperian period (Head & Pratt, 2001; Leonard & Tanaka, 2001). In addition, a 3 Ga old polythermal ice sheet should have covered the entire basin of Isidis Planitia with a maximum thickness of 4.9 km (Guidat et al., 2015; Souček et al., 2015). Since this basin is close to the 0°C isotherm, a glacier could have flowed from the accumulation zone in the highlands down to the ablation zone in Isidis, as mapped by Jöns (1987). Future work should be dedicated to these regions with local GRISLI modeling.

## 5. Discussion and Conclusion

We simulate the climate of Mars 3 Ga, with an adapted topography (without current polar caps), including the feedback from an ice sheet. We conducted an extensive sensitivity study of the main parameters governing the ice sheets. In the parameter space tested the results change less than an factor of 2 and a few 10s of percent for most of the parameters.

We found that the total volume of water to reach the coupled equilibrium water cycle is ~700 m GEL (340 m for the ocean and the same amount for the ice sheet). This budget is plausible, if a significant amount of water has been removed from the atmosphere/hydrosphere/cryosphere in the last 3 Gy, for instance by chemical alteration (Scheller et al., 2021). A recent analysis of seismological data proposes that current mid-crust porosity is filled by liquid water (Wright et al., 2024) a reservoir representing ~1,700 m GEL. In this scenario, the water should infiltrate to large depth. This scenario may be relevant in the presence of an initial cryosphere preventing the infiltration of water (Clifford & Parker, 2001), or by slow percolation through deep tectonic features. In any case, the global inventory proposed here would imply that this deep reservoir was filled after the wet climate proposed herein at 3 Ga, since atmospheric escape does not appear efficient enough to remove so much water (Alsaeed & Jakosky, 2019; Jakosky & Treiman, 2023). Nevertheless recent work shows higher sputtering than expected (Curry et al., 2025) especially during solar storm and the atmospheric escape may have been larger than previously expected.

From our results, the global atmosphere/hydrosphere/cryosphere equilibrium should be reached in a few 100 ky, dominated by the ice sheet timescale. This implies that the ocean lifetime should be larger than a few 100 ky. For shorter timescales, the global atmosphere/hydrosphere/cryosphere would not have time to equilibrate. In comparison, Kamada et al. (2022) estimated that the ocean will disappear in 100 ky due to evaporation assuming a 500 m GEL ocean at the starting point. The main difference between those two models is the treatment of the water cycle. Assuming an initial global reservoir in the ocean only, would lead to an evaporation dominated climate and large accumulations on the highlands. But by construction the feedback of their ice sheet cannot be
evaluated, except if the initial total water volume is enough to sustain it. Our strategy instead is independent of the initial amount of water and allows an estimate of the equilibrium state of the water cycle. It is very interesting that the sublimation timescale from Kamada et al. (2022) and our equilibrium timescale are similar. This fact could be interpreted (if the initial water is in the ocean only) that the equilibrium state may never reached. On the contrary, if the initial state includes ice in the highlands already, the equilibrium water cycle should be reached. Recent geochemical and atmospherical modeling have shown that H₂ could have warmed up early Mars in transient events, each one lasting ≥100 ky (Adams et al., 2025), which should be enough to reach the equilibrium state.

In addition, we estimated the water cycle mass flux. We found that the numbers estimated in Schmidt2022 are consistent with a realistic ice sheet model on the order of ~10¹⁵ kg/y. The very slow motion of the bottom of the ice sheet (typically 1 m/y) would prevent massive erosion. This would explain why there is little trace of this event in the present-day surface geology. Our results suggest continuing to search for the geomorphological remnant of this ice sheet near the northern border at the lowest elevation where the ice flow is maximum (300 m/y), while noting that this region may be covered by more recent material.

We found ice melting at the base of the ice sheet, typically cm/y, which is highly correlated with the basal velocity. As expected there is a strong positive feedback between higher velocity, lower friction and higher basal melting. Locally, the melting could reach 30 cm/y, near the northern lowest altitude region of the ice sheet.

Our results can be compared to Fastook and Head (2015) even if their target region is surrounding Hellas. The main difference with our work is that they assume a water limited scenario (34 m GEL) without an ocean. We include an ocean, but do not fix the total volume in order to estimate the equilibrium state. For scenarios with the maximum thickness of ice of ~4,000 m, both models are in agreement, with basal melting, even if the maximum velocity is 1 m/y for Fastook and Head (2015) and 1 km/y for us. The main reason is the temperature since they are limited by the high altitude (500 m at minimum), whereas our model predict an ice sheet at low altitude (–3,900 m at minimum).

Due to the use of cartesian coordinates, we focused on the major ice sheet covering the Tharsis region. In future work the missing minor accumulation regions around Hellas and more generally around longitude 90°E would be interesting to investigate. Because of this limitation, the total volume of water in the glacier evaluated in this article are minimum numbers. Since accumulation seems to focus on disjointed domains, we could create several local domains to cover all the separated regions to solve this problem. Another option would be to generalize the GRISLI code in spherical coordinates.

Due to the already very heavy computation of ROCKE-3D, we limit here the coupling with three iterations. The mass flux, the ice sheet surface and ice sheet volume should thus be considered as minimum estimates of the final coupled convergence state. In the future, instead of asynchronous coupling, a more integrated coupling could be proposed to be more efficient to reach the final coupled equilibrium. Nevertheless, this approach will have strong limitations due to the very different timescales for atmosphere and cryosphere processes.

## Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

## Data Availability Statement

All simulation data are available at DOI: 10.14768/dd5d7093-3df3-4791-9ef8-c2d7add6b0c7 and in the ROCKE-3D publications archive at https://portal.nccs.nasa.gov/GISS_modelE/ROCKE-3D/publication-supplements/Schmidt2025JGRP130-Mars_Climate_Polar_Ocean_Ice_Sheet_Dynamics/. The source code of the ROCKE-3D model can be downloaded from the NASA GISS code repository server https://simplex.giss.nasa.gov/snapshots/ (version 1.1).

## References

* Adams, D., Scheucher, M., Hu, R., Ehlmann, B. L., Thomas, T. B., Wordsworth, R., et al. (2025). Episodic warm climates on early Mars primed by crustal hydration. *Nature Geoscience*, 18(2), 133–139. https://doi.org/10.1038/s41561-024-01626-8

* Alsaeed, N. R., & Jakosky, B. M. (2019). Mars water and D/H evolution from 3.3 Ga to present. *Journal of Geophysical Research: Planets*, 124(12), 3344–3353. https://doi.org/10.1029/2019je006066

## Acknowledgments

We acknowledge support from the “Institut National des Sciences de l’Univers” (INSU), the “Centre National de la Recherche Scientifique” (CNRS) and “Centre National d’Etudes Spatiales” (CNES) through the “Programme National de Planétologie.” MJW and AI were supported by an award from the NASA Interdisciplinary Consortia for Astrobiology Research Program from the Planetary Science Division of the Science Mission Directorate. Resources supporting this work were provided by the NASA High-End Computing Program through the NASA Center for Climate Simulation at Goddard Space Flight Center. MJW and IA acknowledge support from the GSFC Sellers Exoplanet Environments Collaboration (SEEC), which is funded by the NASA Planetary Science Division’s Internal Scientist Funding Model.
Amundsen, D. S., Mayne, N. J., Baraffe, I., Manners, J., Tremblin, P., Drummond, B., et al. (2016). The UK Met Office global circulation model with a sophisticated radiation scheme applied to the hot Jupiter HD 209458b. *Astronomy and Astrophysics*, **595**, A36. https://doi.org/10.1051/0004-6361/201629183

Baker, V. R., Strom, R. G., Gulick, V. C., Kargel, J. S., Komatsu, G., & Kale, V. S. (1991). Ancient oceans, ice sheets and the hydrological cycle on Mars. *Nature*, **352**(6336), 589–594. https://doi.org/10.1038/352589a0

Benn, D. I., Le Hir, G., Bao, H., Donnadieu, Y., Dumas, C., Fleming, E. J., et al. (2015). Orbitally forced ice sheet fluctuations during the Marinoan Snowball Earth glaciation. *Nature Geoscience*, **8**(9), 704–707. https://doi.org/10.1038/ngeo2502

Bernhardt, H., Hiesinger, H., Ivanov, M., Ruesch, O., Erkeling, G., & Reiss, D. (2016). Photogeologic mapping and the geologic history of the Hellas basin floor, Mars. *Icarus*, **264**, 407–442. https://doi.org/10.1016/j.icarus.2015.09.031

Bouley, S., Baratoux, D., Matsuyama, I., Forget, F., Séjourné, A., Turbet, M., & Costard, F. (2016). Late Tharsis formation and implications for early Mars. *Nature*, **531**(7594), 344–347. https://doi.org/10.1038/nature17171

Bouley, S., Baratoux, D., Paulien, N., Missenard, Y., & Saint-Bézar, B. (2018). The revised tectonic history of Tharsis. *Earth and Planetary Science Letters*, **488**, 126–133. https://doi.org/10.1016/j.epsl.2018.02.019

Bouquety, A., Sejourné, A., Costard, F., Bouley, S., & Leyguarda, E. (2020). Glacial landscape and paleoglaciation in Terra Sabaea: Evidence for a 3.6 Ga polythermal plateau ice cap. *Geomorphology*, **350**, 106858. https://doi.org/10.1016/j.geomorph.2019.106858

Bouquety, A., Sejourné, A., Costard, F., Mercier, D., & Bouley, S. (2019). Morphometric evidence of 3.6 Ga glacial valleys and glacial cirques in martian highlands: South of Terra Sabaea. *Geomorphology*, **334**, 91–111. https://doi.org/10.1016/j geomorph.2019.02.022

Chassefiere, E., Lasue, J., Langlais, B., & Quesnel, Y. (2016). Early Mars serpentinization-derived CH₄ reservoirs, H₂-induced warming and paleopressure evolution. *Meteoritics & Planetary Science*, **51**(11), 2234–2245. https://doi.org/10.1111/maps.12784

Clifford, S. M., & Parker, T. J. (2001). The evolution of the Martian hydrosphere: Implications for the fate of a primordial ocean and the current state of the northern plains. *Icarus*, **154**(1), 40–79. https://doi.org/10.1006/icar.2001.6671

Contoux, C., Dumas, C., Ramstein, G., Jost, A., & Dolan, A. (2015). Modelling Greenland ice sheet inception and sustainability during the Late Pliocene. *Earth and Planetary Science Letters*, **424**, 295–305. https://doi.org/10.1016/j.epsl.2015.05.018

Costard, F., Séjourné, A., Kelfoun, K., Clifford, S., Lavigne, F., Pietro, I. D., & Bouley, S. (2017). Modeling tsunami propagation and the emplacement of thumbprint terrain in an early Mars ocean. *Journal of Geophysical Research: Planets*, **122**(3), 633–649. https://doi.org/10.1002/2016je005230

Costard, F., Séjourné, A., Lagain, A., Ormö, J., Rodriguez, J., Clifford, S., et al. (2019). The Lomonosov crater impact event: A possible mega-tsunami source on Mars. *Journal of Geophysical Research: Planets*, **124**(7), 1840–1851. https://doi.org/10.1029/2019je006008

Crotti, I., Quiquet, A., Landais, A., Stenni, B., Wilson, D. J., Severi, M., et al. (2022). Wilkes subglacial basin ice sheet response to Southern Ocean warming during late Pleistocene interglacials. *Nature Communications*, **13**(1), 5328. https://doi.org/10.1038/s41467-022-32847-3

Curry, S. M., Hara, T., Luhmann, J. G., Leblanc, F., Jolitz, R., Mitchell, D., et al. (2025). First direct observations of atmospheric sputtering at Mars. *Science Advances*, **11**(22), eadt1538. https://doi.org/10.1126/sciadv.adt1538

Dickeson, Z. I., & Davis, J. M. (2020). Martian oceans. *Astronomy and Geophysics*, **61**(3), 311–317. https://doi.org/10.1093/astrogeo/ataa038

Duran, S., & Coulthard, T. J. (2020). The Kasei Valles, Mars: A unified record of episodic channel flows and ancient ocean levels. *Scientific Reports*, **10**(1), 18571. https://doi.org/10.1038/s41598-020-75080-y

Edwards, J. M., & Slingo, A. (1996). Studies with a flexible new radiation code. I: Choosing a configuration for a large-scale model. *Quarterly Journal of the Royal Meteorological Society*, **122**(531), 689–719. https://doi.org/10.1002/qj.49712253107

Fastook, J. L., & Head, J. W. (2015). Glaciation in the Late Noachian Icy Highlands: Ice accumulation, distribution, flow rates, basal melting, and top-down melting rates and patterns. *Planetary and Space Science*, **106**, 82–98. https://doi.org/10.1016/j.pss.2014.11.028

Forget, F., Wordsworth, R., Millour, E., Madeleine, J.-B., Kerber, L., Leconte, J., et al. (2013). 3D modelling of the early martian climate under a denser CO₂ atmosphere: Temperatures and CO₂ ice clouds. *Icarus*, **222**(1), 81–99. https://doi.org/10.1016/j.icarus.2012.10.019

Galofre, A. G., Jellinek, A. M., & Osinski, G. R. (2020). Valley formation on early Mars by subglacial and fluvial erosion. *Nature Geoscience*, **13**(10), 663–668. https://doi.org/10.1038/s41561-020-0618-x

Gourronc, M., Bourgeois, O., Mège, D., Pochat, S., Bultel, B., Massé, M., et al. (2014). One million cubic kilometers of fossil ice in Valles Marineris: Relicts of a 3.5 Gy old glacial land system along the martian equator. *Geomorphology*, **204**, 235–255. https://doi.org/10.1016/j. geomorph.2013.08.009

Guidat, T., Pochat, S., Bourgeois, O., & Souček, O. (2015). Landform assemblage in Isidis Planitia, Mars: Evidence for a 3 Ga old polythermal ice sheet. *Earth and Planetary Science Letters*, **411**, 253–267. https://doi.org/10.1016/j.epsl.2014.12.002

Haberle, R. M., Zahnle, K., Barlow, N. G., & Steakley, K. E. (2019). Impact degassing of H₂ on early Mars and its effect on the climate system. *Geophysical Research Letters*, **46**(22), 13355–13362. https://doi.org/10.1029/2019gl084733

Head, J. W., & Pratt, S. (2001). Extensive Hesperian-aged south polar ice sheet on Mars: Evidence for massive melting and retreat, and lateral flow and ponding of meltwater. *Journal of Geophysical Research*, **106**(E6), 12275–12299. https://doi.org/10.1029/2000je001359

Hepburn, A. J., Ng, F. S. L., Holt, T. O., & Hubbard, B. (2020). Late Amazonian ice survival in Kasei Valles, Mars. *Journal of Geophysical Research: Planets*, **125**(11), e06531. https://doi.org/10.1029/2020JE006531

Hildenbrand, A., Zeyen, H., Schmidt, F., Bouley, S., Costard, F., Gillot, P., et al. (2023). A giant volcanic island in an early Martian Ocean? *Earth and Planetary Science Letters*, **619**, 118302. https://doi.org/10.1016/j.epsl.2023.118302

Hynek, B. M., Beach, M., & Hoke, M. R. T. (2010). Updated global map of Martian valley networks and implications for climate and hydrologic processes. *Journal of Geophysical Research*, **115**(E9), E90008. https://doi.org/10.1029/2009je003548

Jakosky, B. M., & Treiman, A. H. (2023). Mars volatile inventory and outgassing history. *Icarus*, **402**, 115627. https://doi.org/10.1016/j.icarus. 2023.115627

Jöns, H.-P. (1987). Large fossil mud lakes or giant mud sheet floods in Syrtis Major (Isidis Planitia) and Mare Australe, Mars. In *Lunar and planetary science conference* (Vol. 18, p. 470).

Kamada, A., Kuroda, T., Kasaba, Y., Terada, N., & Nakagawa, H. (2021). Global climate and river transport simulations of early Mars around the Noachian and Hesperian boundary. *Icarus*, **368**, 114618. https://doi.org/10.1016/j.icarus.2021.114618

Kamada, A., Kuroda, T., Kodama, T., Kasaba, Y., & Terada, N. (2022). Evolution of ice sheets on early Mars with subglacial river systems. *Icarus*, **385**, 115117. https://doi.org/10.1016/j.icarus.2022.115117

Kite, E. S., Steele, L. J., Mischna, M. A., & Richardson, M. I. (2021). Warm early Mars surface enabled by high-altitude water ice clouds. *Proceedings of the National Academy of Sciences*, **118**(18), e2101959118. https://doi.org/10.1073/pnas.2101959111

Kite, E. S. (2019). Geologic constraints on early Mars climate. *Space Science Reviews*, **215**(1), 10. https://doi.org/10.1007/s11214-018-0575-5

Kite, E. S., & Conway, S. (2024). Geological evidence for multiple climate transitions on early Mars. *Nature Geoscience*, **17**(1), 10–19. https://doi. org/10.1038/s41561-023-01349-2
Ladant, J.-B., Donnadieu, Y., Lefebvre, V., & Dumas, C. (2014). The respective role of atmospheric carbon dioxide and orbital parameters on ice sheet evolution at the Eocene-Oligocene transition. *Paleoceanography*, 29(8), 810–823. https://doi.org/10.1002/2013PA002593

Laskar, J., Correia, A. C. M., Gastineau, M., Joutel, F., Levrard, B., & Robutel, P. (2004). Long term evolution and chaotic diffusion of the insolation quantities of Mars. *Icarus*, 170(2), 343–364. https://doi.org/10.1016/j.icarus.2004.04.005

Le clec'h, S., Charbit, S., Quiquet, A., Fettweis, X., Dumas, C., Kageyama, M., et al. (2019). Assessment of the Greenland ice sheet–atmosphere feedbacks for the next century with a regional atmospheric model coupled to an ice sheet model. *The Cryosphere*, 13(1), 373–395. https://doi.org/10.5194/tc-13-373-2019

Le Meur, E., & Huybrechts, P. (1996). A comparison of different ways of dealing with isostasy: Examples from modeling the Antarctic ice sheet during the last glacial cycle. *Annals of Glaciology*, 23, 309–317. https://doi.org/10.1017/s0260305500013586

Leonard, G. J., & Tanaka, K. L. (2001). Geologic map of the Hellas region of Mars. In *U.S. Geological Survey Report*. https://doi.org/10.3133/i2694.2694.

Lucchitta, B. K. (1980). Martian outflow channels sculptured by glaciers, II. In *Lunar and planetary science conference* (pp. 634–636).

Lucchitta, B. K. (1995). Kasei Valles, Mars (II): The ice story. In *Lunar and planetary science conference* (Vol. 26, p. 869).

Manners, J., Vosper, S., & Roberts, N. (2012). Radiative transfer over resolved topographic features for high-resolution weather prediction. *Quarterly Journal of the Royal Meteorological Society*, 138(664), 720–733. https://doi.org/10.1002/qj.956

Matsuyama, I., & Manga, M. (2010). Mars without the equilibrium rotational figure, Tharsis, and the remnant rotational figure. *Journal of Geophysical Research*, 115(E12), E12020. https://doi.org/10.1029/2010je003686

Mège, D., & Bourgeois, O. (2011). Equatorial glaciations on Mars revealed by gravitational collapse of Valles Marineris wallslopes. *Earth and Planetary Science Letters*, 310(3–4), 182–191. https://doi.org/10.1016/j.epsl.2011.08.030

Perron, J. T., Mitrovica, J. X., Manga, M., Matsuyama, I., & Richards, M. A. (2007). Evidence for an ancient martian ocean in the topography of deformed shorelines. *Nature*, 447(7146), 840–843. https://doi.org/10.1038/nature05873

Quiquet, A., & Dumas, C. (2021a). The GRISLI-LSCE contribution to the Ice Sheet Model Intercomparison Project for phase 6 of the Coupled Model Intercomparison Project (ISMIP6) – Part 1: Projections of the Greenland ice sheet evolution by the end of the 21st century. *The Cryosphere*, 15(2), 1015–1030. https://doi.org/10.5194/tc-15-1015-2021

Quiquet, A., & Dumas, C. (2021b). The GRISLI-LSCE contribution to the Ice Sheet Model Intercomparison Project for phase 6 of the Coupled Model Intercomparison Project (ISMIP6) – Part 2: Projections of the Antarctic ice sheet evolution by the end of the 21st century. *The Cryosphere*, 15(2), 1031–1052. https://doi.org/10.5194/tc-15-1031-2021

Quiquet, A., Dumas, C., Paillard, D., Ramstein, G., Ritz, C., & Roche, D. M. (2021). Deglacial ice sheet instabilities induced by proglacial lakes. *Geophysical Research Letters*, 48(9), e2020GL092141. https://doi.org/10.1029/2020gl092141

Quiquet, A., Dumas, C., Ritz, C., Peyaud, V., & Roche, D. M. (2018). The GRISLI ice sheet model (version 2.0): Calibration and validation for multi-millennial changes of the Antarctic ice sheet. *Geoscientific Model Development*, 11(12), 5003–5025. https://doi.org/10.5194/gmd-11-5003-2018

Quiquet, A., Ritz, C., Punge, H. J., & Salas y Mélia, D. (2013). Greenland ice sheet contribution to sea level rise during the last interglacial period: A modelling study driven and constrained by ice core data. *Climate of the Past*, 9(1), 353–366. https://doi.org/10.5194/cp-9-353-2013

Reeh, N. (1991). Parameterization of melt rate and surface temperature on the Greenland ice sheet. *Polarforschung*, 59, 113–128.

Rodriguez, J. A. P., Fairén, A. G., Tanaka, K. L., Zarroca, M., Linares, R., Platz, T., et al. (2016). Tsunami waves extensively resurfaced the shorelines of an early Martian ocean. *Scientific Reports*, 6(1), 25106. https://doi.org/10.1038/srep25106

Scheller, E. L., Ehlmann, B. L., Hu, R., Adams, D. J., & Yung, Y. L. (2021). Long-term drying of Mars by sequestration of ocean-scale volumes of water in the crust. *Science*, 372(6537), 56–62. https://doi.org/10.1126/science.abc7717

Schmidt, F., Way, M. J., Costard, F., Bouley, S., Séjourné, A., & Aleinov, I. (2022). Circumpolar ocean stability on Mars 3 Gy ago. *Proceedings of the National Academy of Sciences*, 119(4), e2112930118. https://doi.org/10.1073/pnas.2112930118

Souček, O., Bourgeois, O., Pochat, S., & Guidat, T. (2015). A 3 Ga old polythermal ice sheet in Isidis Planitia, Mars: Dynamics and thermal regime inferred from numerical modeling. *Earth and Planetary Science Letters*, 426, 176–190. https://doi.org/10.1016/j.epsl.2015.06.038

Tan, N., Ramstein, G., Dumas, C., Contoux, C., Ladant, J.-B., Sepulchre, P., et al. (2017). Exploring the MIS M2 glaciation occurring during a warm and high atmospheric CO₂ Pliocene background climate. *Earth and Planetary Science Letters*, 472, 266–276. https://doi.org/10.1016/j.epsl.2017.04.050

Tanaka, K. L., Skinner, J. A., Jr., Dohm, J. M., Irwin, III, R. P., Kolb, E. J., Fortezzo, C. M., et al. (2014). Geologic map of Mars, *U.S. Geological Survey Report*. https://doi.org/10.3133/sim3292.3292.

Tosca, N. J., Ahmed, I. A. M., Tutolo, B. M., Ashpitel, A., & Hurowitz, J. A. (2018). Magnetite authigenesis and the warming of early Mars. *Nature Geoscience*, 11(9), 635–639. https://doi.org/10.1038/s41561-018-0203-8

Turbet, M., Boulet, C., & Karman, T. (2020). Measurements and semi-empirical calculations of CO₂ + CH₄ and CO₂ + H₂ collision-induced absorption across a wide range of wavelengths and temperatures. Application for the prediction of early Mars surface temperature. *Icarus*, 346, 113762. https://doi.org/10.1016/j.icarus.2020.113762

Turbet, M., & Forget, F. (2019). The paradoxes of the Late Hesperian Mars ocean. *Scientific Reports*, 9(1), 5717. https://doi.org/10.1038/s41598-019-42030-2

Turbet, M., Forget, F., Head, J. W., & Wordsworth, R. (2017). 3D modelling of the climatic impact of outflow channel formation events on early Mars. *Icarus*, 288, 10–36. https://doi.org/10.1016/j.icarus.2017.01.024

Tutolo, B. M., & Leong, J. A. (2025). Serpentine solid solutions and hydrogen production on early Earth and Mars. *Elements*, 21(3), 177–183. https://doi.org/10.2138/gselements.21.3.177

van Aalderen, V., Charbit, S., Dumas, C., & Quiquet, A. (2024). Relative importance of the mechanisms triggering the Eurasian ice sheet deglaciation in the GRISLI2.0 ice sheet model. *Climate of the Past*, 20(1), 187–209. https://doi.org/10.5194/cp-20-187-2024

Wang, L., & Huang, J. (2024). Hypothesis of an ancient northern ocean on Mars and insights from the Zhurong rover. *Nature Astronomy*, 8(10), 1220–1229. https://doi.org/10.1038/s41550-024-02343-3

Way, M. J., Aleinov, I., Amundsen, D. S., Chandler, M. A., Clune, T. L., Del Genio, A. D., et al. (2017). Resolving orbital and climate keys of Earth and extraterrestrial environments with dynamics (ROCKE-3D) 1.0: A general circulation model for simulating the climates of rocky planets. *The Astrophysical Journal - Supplement Series*, 231(1), 12. https://doi.org/10.3847/1538-4365/aa7a06

Wordsworth, R., Forget, F., Millour, E., Head, J., Madeleine, J.-B., & Charnay, B. (2013). Global modelling of the early martian climate under a denser CO₂ atmosphere: Water cycle and ice evolution. *Icarus*, 222(1), 1–19. https://doi.org/10.1016/j.icarus.2012.09.036

Wordsworth, R., Kalugina, Y., Lokshtanov, S., Vigasin, A., Ehlmann, B., Head, J., et al. (2017). Transient reducing greenhouse warming on early Mars. *Geophysical Research Letters*, 44(2), 665–671. https://doi.org/10.1002/2016gl071766

Wright, V., Morzfeld, M., & Manga, M. (2024). Liquid water in the Martian mid-crust. *Proceedings of the National Academy of Sciences*, 121(35), e2409983121. https://doi.org/10.1073/pnas.2409983121