Circumpolar ocean stability on Mars 3 Gy ago.

What was the nature of the Late Hesperian climate, warm and wet or cold and dry? Formulated this way the question leads to an apparent paradox since both options seem implausible. A warm and wet climate would have produced extensive fluvial erosion but few valley networks have been observed at the age of the Late Hesperian. A too cold climate would have kept any northern ocean frozen most of the time. A moderate cold climate would have transferred the water from the ocean to the land in the form of snow and ice. But this would prevent tsunami formation, for which there is some evidence. Here, we provide insights from numerical climate simulations in agreement with surface geological features to demonstrate that the Martian climate could have been both cold and wet. Using an advanced general circulation model (GCM), we demonstrate that an ocean can be stable, even if the Martian mean surface temperature is lower than 0 °C. Rainfall is moderate near the shorelines and in the ocean. The southern plateau is mostly covered by ice with a mean temperature below 0 °C and a glacier return flow back to the ocean. This climate is achieved with a 1-bar CO2-dominated atmosphere with 10% H2 Under this scenario of 3 Ga, the geologic evidence of a shoreline and tsunami deposits along the ocean/land dichotomy are compatible with ice sheets and glacial valleys in the southern highlands.

The possibility of a late Martian ocean is a topic of debate with strong implications on the habitability of the Red Planet (1). A recent review (2) discusses this controversy. There is evidence of Martian paleoshorelines (3) in Deuteronilus Mensae (sometimes noted contact 2) in a geometry closely corresponding to the current equipotential height (4). The Deuteronilus shoreline seems to have formed during the last stage of the true polar wander induced by Tharsis (5). A northern ocean is also supported by specific radar properties (6), smooth surface roughness (7), and a fractal analysis of the topography (8). Crater count dating of the Vastitas Borealis Formation near Deuteronilus Mensae is 3.5 Ga (9) but the ocean may have been more recent. Detailed studies in Kasei Valles imply that such an ocean rose in elevation (∼1,000 m) between ca. 3.6 and 3.2 Ga (10). Along this shoreline, candidate tsunami deposits have been identified at an age of 3 Ga (11, 12) with at least two impact events. The Lomonosov crater morphology is coherent with an impact in shallow water, at the very same age as the tsunami deposits (13), and it is thus the most probable source.

The stability of an ocean in a warm and wet scenario, even transient, at 3 Ga has been questioned since intense and widespread rainfall in such a scenario appeared inconsistent with the few observed dendritic valley networks from this time (14, 15). In previous work, a cold and wet Mars seemed impossible since the long-term stability of an ocean (open or ice covered) in such a scenario has never been achieved by a three-dimensional general circulation model (3D-GCM) (14–19). Mainly the water was found to accumulate in the form of ice in the southern highlands, but no glacier/ice sheet processes were considered in these studies.

Cool and wet scenarios in the case of a faint young sun (at the Noachian-Hesperian boundary at 3.6 Ga) demonstrate that the ocean freezes for pressures below 1 bar and yet predict intense rainfall when the global average temperature is higher than 0 C (20). Nevertheless, liquid brines are possible in this case but only with significant antifreezing properties (21). An investigation with two-dimensional (2D) simulations for early Mars argues for a warm and semiarid early Mars (22) but stresses the need for a coupled ocean/atmosphere model.

The 3D Climate Model

We present fully coupled ocean/atmosphere 3D-GCM simulations based on ROCKE-3D (Resolving Orbital and Climate Keys of Earth and Extraterrestrial Environments with Dynamics, hereafter R3D) (23), which is based upon a parent Earth climate model known as ModelE2 (24). R3D allows us to estimate the interaction between atmosphere and ocean circulation but also encompasses a surface hydrological scheme. We assume the solar luminosity to be ∼79% (25) of its current value. Hence at 3 Ga the flux at Mars was set to 452.8 Wm–2 in the GCM. The total water budget is 150 m global equivalent layer (GEL) to fill the ocean basin in the north and Hellas up to –3,900 m. Our model includes a dynamic fully coupled ocean while sea ice is also interactive and determined by the salinity and corresponding freezing point of water (23). The salinity was set to modern Earth ocean values at model start given the lack of constraints for ancient Mars, but R3D is capable of exploring a range of salinities (26). In addition, we include a glacier flux in R3D to simulate the return flow to the ocean. This scheme allows us to estimate the glacier mass balance, but glacier thickness (including total mass) and flowing path cannot be constrained by our model. In any case, the water budget must be compatible with the deuterium to hydrogen ratio (D/H) measurements that imply 100 to 300 m GEL in the Hesperian (27, 28), and hence we have chosen 150 m GEL.

Four sets of numerical simulations were performed for 0, 20, 40, and 60 obliquity, with 1 bar atmospheric pressure (15, 29). An obliquity of 40 seems to be the most likely (30). We set eccentricity to 0 since it has only a minor effect on mean annual temperature (31, 32). For 20% H2 in a CO2-dominated atmosphere the greenhouse effect is efficient enough to keep most of the planet’s surface temperature above freezing () as expected (20). The global mean surface temperature is always greater than 10 C with a stable ocean (northern ocean and Hellas basin sea) including intense rainfall, corresponding to the warm and wet scenario (33). For 10% H2 the globally averaged surface temperature is below 0 C but the northern ocean surface temperature remains surprisingly around 7 C and thus the ocean is prevented from freezing. The source of H2 may be volcanic outgassing or serpentinization, but such large concentrations are not expected to persist for more than 1 My (34–37). Nevertheless, volcanic activities in Tharsis and Elysium have been identified to cover a very large period of time, from 3.8 Ga to 200 Ma (38). Alternatively, H2 from impact degassing has also been proposed (39). The lower limit of H2 concentration necessary to stabilize the ocean will be a subject of a future investigation.

Fig. 1 shows the simulated surface fields averaged over 10 Martian years for H2 = 10% and obliquity 40. Interestingly the climate changes only slightly with different obliquities () likely because the mean ground albedo does not change significantly (e.g., 23% versus 28% for 0 versus 60 obliquity). The major effects of water, ice, and snow albedo dependence on incidence angles are included (Materials and Methods). The relatively modest effects of obliquity are probably due to the 1-bar atmosphere that more efficiently transports heat than the present-day thin 6-mbar atmosphere.

Fig. 1.

The 3D-GCM output at 40 obliquity and H2 = 10% for the rain precipitation, snow and ice fraction at the surface, snowfall, and sea/ground surface temperature. Black contour lines represent surface elevation level (–2,000, 0, 2,000, and 8,000 m) and the red contour line is the paleoshoreline (–3,900 m).

Despite an average planetary surface temperature below 0 C, the ocean remains above freezing due to its low altitude and low albedo (). In addition, the ability of the ocean to transport heat is responsible for nonnegligible higher local temperatures. Fig. 2, Top shows the net vertical heat flux contributing to the ocean’s stabilization. Ocean gyres are transporting heat poleward, as expected for a fast-rotating planet (40). The importance of ocean heat transport in a dynamic ocean versus a slab ocean was quantified in refs. 26 and 40, and the ability to extract wind-driven ocean upwelling-related nutrient supply for life from R3D simulations has also recently been demonstrated (41). Fig. 2, Bottom presents the surface temperature increase due to the active circulation of the ocean as a function of latitude and obliquity. The warming is globally present at all latitudes but locally higher near the northern polar ocean, from 1 C up to 4.5 C. This effect is larger for extreme obliquities 0 and 60 due to the insolation regime at the poles.

Fig. 2.

(Top) Net outward heat flux transported by the ocean for 40 obliquity and H2 = 10%. The positive value (up to 15 Wm–2 at 60 obliquity) near the North pole indicates that heat goes toward the atmosphere due to ocean circulation. For a slab ocean model, this flux is null. Black contour lines represent surface elevation level and the red contour line is the paleoshoreline. (Bottom) Latitudinal profile of the temperature increase () due to the circulating ocean, as a function of obliquity. In the slab case, the ocean is assumed stratified but noncirculating.

On land, there is a clear boundary at the 0 C isotherm, which corresponds approximately to the Martian dichotomy. The altitude of the boundary varies from –700 to –3,000 m for obliquity from 0 to 60 ().

In the high-altitude domain, commonly referred to as the “icy highlands,” the surface is mostly frozen and snow precipitation is dominant. The extensive accumulation of snow in the highlands can lead to the formation of significant ice sheets that may flow down to the northern and Hellas basin oceans. Our current model is not able to simulate details regarding the glacier flow, but only a global mass flux from land to the ocean. The simplified mass flux is estimated by restricting snow accumulation to 2 m of H2O. Any excess is immediately redistributed over the ocean to maintain conservation of water in the system. Glacial processes such as accumulation (snow compaction), flow (path, erosion), and ablation (direct melting, fluvio-glacial river) are thus not included. Nevertheless, the global mass flux (Materials and Methods) is estimated to be around kgy–1 and seems to be very robust to obliquity changes and H2 content. For comparison, this flux is currently 1 to kgy–1 on the Earth today (42). In the Last Glacial Maximum, closer to our Martian cold and wet climate, the water discharge was estimated to be 1016 kgy–1 for North America only (43). Given the relative sizes of the planetary bodies, our value seems reasonable.

In the lowest-altitude domain, called the “wet lowlands,” rain, evaporation, and surface runoff are balanced. Wet lowlands represent a minor surface fraction (22% at 40 obliquity) of the planet for H2 = 10% (Fig. 3).

Fig. 3.

Proposed scenario of a cold and wet climate at the Hesperian age (3 Ga). (Top) Map of the climatic zones at 40 obliquity: ocean, wet lowlands, and icy highlands, separated respectively by –3,900 m and the 0 C isotherm (–1,990 m). The topography is modern topography without the North Polar Layered Deposits. Geological evidence of past glaciers with possible paths indicating return flow to the oceans is highlighted. Geological evidence of high drainage density valley networks at 3 Ga is rare but noted above. (Bottom) Simplified scheme of the hydrological cycle using box modeling. Fluxes in parentheses have units of 1015 kgy–1.

The rain occurs mainly over the ocean (∼70% of the total rain), and ∼60% of the ocean evaporation is compensated by rain. The snow precipitation over the icy highlands is intense but ∼80% is compensated by high sublimation rates. Our climate model differs from the “icy highlands model” (17, 31) since the ocean is stable due to glacier return flow.

From these results, one can describe a cold and wet climatic regime (Fig. 3, Bottom) with significant glacier return flow from the cold highlands back to the ocean with moderate surface runoff near the shoreline.

Comparison to Geological Evidence

The proposed scenario is in agreement with the long-term stability of the ocean, with the paleoshoreline and potential tsunami deposits. In addition, this climatic regime is in agreement with independent geological evidence concerning fluvial and glacial processes. Fig. 3 plots the major geological features of a similar late Hesperian/early Amazonian age, not covered by subsequent units ().

Contrary to the Noachian epoch, the modest appearance and development of dendritic valleys (rainfall related) during the late Hesperian/early Amazonian have been tied to major changes in climate. Several valleys at 2.9 Ga have been identified in the North of Alba Patera with a substantial temporal peak in the drainage density (44), located in the dichotomy boundary, in agreement with large rainfall predicted by our model (Fig. 1). In addition, a relatively high-density network has been located in the South of Hellas basin (45). Along the highland/lowland boundary (northwestern Arabia Terra and in Deuteronilus Mensae regions) network-aligned channels are observed displaying streamlined islands of Hesperian age of possible marine (11) or fluvial origin (46). The wet lowlands have been substantially covered by more recent Amazonian processes (including volcanic resurfacing) (46), erasing most of the geological features of this age available at the surface.

Global extensive ice sheets have been proposed in the southern hemisphere. Some deposits in the region of Malea Planum (Pityusa Patera) are possible remnants of an extensive polar ice sheet emplaced during the Hesperian period (47, 48). 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 (49, 50). 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 (51).

Reull Vallis is a 1,500-km-long U-shaped valley postdating the Hesperian ridged plains with a long and complex history (52). It has been carved by glacial flow along the channel during a younger resurfacing episode of late Hesperian–early Amazonian age (53). With a length of ∼3,000 km, ∼500 km width, and a depth of ∼3 km, Kasei Valles is the largest outflow channel on Mars. Flooding activity occurred mostly at the Hesperian between 3.6 and 3.8 Ga, with continuing activity possible to 2.5 ± 1 Ga (54). The morphology of the Kasei Valles region is consistent with an origin by ice with a striking analogy to Alaskan glaciers (55) and Antarctic ice streams (56). A glacial interpretation of the Kasei Valles outflow was also proposed to explain the presence of scour features, tunnel valleys, and esker forms (57). Ares Vallis is a 1,500-km-long Hesperian outflow channel displaying multiple catastrophic flooding events. Its floor shows geomorphological evidence of ice-covered flooding (58) and glacial erosion (59) with subsequent thermokarst degradation (60). The presence of thermokarst lakes formed after flooding (60, 61) supports a hypothesis of the thawing of ground ice during the Hesperian (i.e., 3.6 to 3.0 Ga). In addition, evidence of past glaciations (between 1.4 and 3.5 Ga) was reported within Valles Marineris. Morphological evidence shows extensive subglacial erosion of the lower parts of Valles Marineris wallslopes, together with the occurrence of sackung deformations by deglaciation of the upper wallslopes (62, 63).

Finally, one has to note that the extensive outflow channels are contemporaneous with the ocean. An ocean at the Hesperian is also compatible with the timing of peak outflow channel activities (64, 65). Outflows probably formed by aquifer disruption (66) or glacial erosion (55, 59). More abundant signs of glaciation in the icy highlands may have been prevented by less erosive cold-based glaciers (67) and by postdating Amazonian processes, such as volcanic events (46).

One also has to consider the arguments that confront the stable ocean scenario. First, the interpretation and location of the paleo shoreline itself are matters of controversy (68, 69). Some studies claim that the putative shorelines are mutually inconsistent (70) but it mainly concerns contact 1 (Arabia potential shoreline), not contact 2 (younger Deuteronilus shoreline) that we use in our study. Second, mineralogical data seem inconsistent with a Late Hesperian ocean (71, 72). In particular, mineralogy from impact crater excavated terrains (73) and clay occurrence in central peaks (74) are more consistent with a volcanic origin. But widespread sedimentary rock has been identified in the northern plains of Mars (75). Finally, there is a lack of observed carbonates at the surface that implies that the ocean was acidic. However, an acidic ocean seems unfavorable in the presence of mafic rocks (76). Alternative scenarios to an ocean have been proposed implying lava flows (77) or ice-related processes (64, 78) creating the smooth northern plains we observe today. In addition, crater ejecta have been proposed as an alternative to tsunami deposits to form thumbprint terrains (29).

Conclusion

Our results demonstrate that a cold and wet climate could have been stable and is consistent with geomorphological evidence due to glacier return flow. Using fully coupled atmosphere/dynamic ocean modeling, we show that the ocean’s circulation regionally warms the surface up to 4.5 C. In these conditions the ocean is stable even for a global mean planetary temperature below 0 C. This Martian climate may be similar to ancient Earth’s with an active water cycle around the time of the early stages of life’s appearance. Future work should encompass careful analysis of this stable ocean domain and its application to the past Martian climate, especially in the Noachian epoch. In addition, we hope that this work will stimulate more precise ice sheet simulations on Mars at a global scale, in comparison with detailed photo-geomorphological studies. Further in situ analysis of the ocean boundary, for instance by Zhurong and ExoMars Rosalind Franklin rover, should confirm the shoreline and tsunami deposits interpretation. In the longer term, Mars Ice Mapper and particularly its high-resolution radar capability will also be able to discover additional clues.

Materials and Methods

Atmosphere 3D-GCM

The 3D-GCM simulations performed herein utilize R3D (23), which is based upon a parent Earth climate model called ModelE2 (24). ModelE2 has been used for the Coupled Model Intercomparison Project Phase 6 (CMIP6) (79). R3D has a number of extensions in areas relative to non-Earth terrestrial worlds (23). In particular, R3D utilizes a separate and highly flexible radiative transfer scheme known as SOCRATES (Suite Of Community RAdiative Transfer codes based on Edwards and Slingo) (80, 81) that allows for non-Earth atmospheric gas constituents and abundances that are outside the capabilities of the native ModelE2 scheme. R3D also extends the ModelE2 calendaring system and allows for worlds with different rotational and orbital parameters. In this project we use 4 horizontal resolution for all components. The atmosphere has 40 layers with the model top at 0.1 hPa. The major components of surface and atmosphere exchange data at a time step of 1,850 s, although shorter time steps are used internally. We set the solar luminosity to be ∼79% (25) of its current value (1,360.67 Wm–2 at Earth) following ref. 82. Hence at 3 Ga, the flux at Mars was set to 452.8 Wm–2. For this work we rely upon CO2 + H2 Collisionally Induced Absorption (CIA) tables from ref. 36 incorporated into R3D. Recent work (83) has demonstrated the possibility that ref. 36 may have overestimated the warming provided in their CO2 + H2 CIA calculations, confirmed by ref. 84. At the same time the R3D SOCRATES pure CO2 absorption is underestimated, allowing ∼7.5 Wm2 additional outgoing long-wave radiation compared with Spectral Mapping Atmospheric Radiative Transfer (SMART) (85) calculations as shown in figure 2 of ref. 86. The overestimated CO2 + H2 CIA R3D uses from ref. 36 tends to cancel out the underestimated CO2 absorption as demonstrated in experiments conducted in ref. 86. A more detailed explanation of these effects can be found in ref. 86. However, even if we are underestimating the absorption, we have “bounded” the problem in this work by providing GCM results for a 20% H2 atmosphere for the four different obliquities in .

H2 provides a powerful greenhouse component in combination with CO2 as a background gas, but other gas combinations involving CH4 or H2S may have an equivalent radiative effect at sufficient mixing ratios even if the motivation for their use is lacking. Recent measurements of collision-induced absorption of CO2 and CH4 show that this possibility is not favored (83). We run simulations with 10% and 20% H2 in a CO2-dominated atmosphere at 1 bar surface pressure to simulate the climate at 3 Ga, as in ref. 15. We run simulations with 0, 20, 40, and 60 obliquity, since this orbital parameter can have large excursions from the mean value in the past (30, 87).

H2O cloud implementation in R3D is described in detail in refs. 23, 24, and 86. The ice particles are assumed to be Mie scatterers while the parameterization of cirrus ice cloud properties follows ref. 88 via our use of the SOCRATES radiation scheme in R3D. R3D does not include CO2 ice clouds that are presently under development. For the relatively warm temperatures of our simulations they are expected to play a minor role in the radiation budget of the atmosphere, although for cooler climates they are known to play an important role (89).

The surface radiative properties of R3D are described in ref. 90. Dry soil albedo is set to 15%. The spectral and solar zenith angle dependence of ocean albedo is based on calculations of Fresnel reflection from wave surface distributions as a function of wind velocity (91). The effects of foam on ocean albedo are also included (92). In addition, R3D also includes the dependence of the cosine of the zenith angle for snow albedo both on ocean ice and on land (93).

The main point of comparison with ref. 20 is the global average temperature for 20% H2. They found 15.65 C (for us 16.2 ∘C for slab ocean) for an obliquity of 25.19 (20 for us) and a solar flux 441.1 Wm–2 (452.8 Wm–2 for our later epoch). Given the relative numbers, we conclude that both models are in agreement.

Ocean Circulation Model

The ocean shoreline is set to –3,900 m in all runs (4, 9). The bathymetry is assumed to be identical to present time except that the northern polar cap was removed. According to modern topography analysis, the Late Hesperian stage of the paleoshoreline is in agreement with current altitude, i.e., without a correction from True Polar Wander due to Tharsis (4, 5, 94, 95). The ocean horizontal resolution is the same as the atmosphere, 4 × 5, and can extend to 13 layers or 4,647 m. The ocean bathymetry is determined by the topography. In Hellas basin the ocean depth reaches 3,194 m (level 12 of 13) while in the northern ocean its deepest extent is 1,294 m (level 10 of 13). The mesoscale diffusivity is fixed at 1,200 ms–1. At model start the ocean is initialized as liquid without any sea ice.

For all simulations herein the dynamic ocean model assumes Earth’s values of salinity (35 practical salinity units) at model start. The ocean includes a number of additional dynamical processes that change as the model moves forward in time such as sea ice growth/shrinkage, salinity change due to fresh water incorporation by surface runoff, precipitation/evaporation, and associated changes in sea height, etc. Additional details on the ocean model can be found in refs. 23, 24, and 26.

The ocean surface fraction is ∼16%, which is small in comparison to Earth at ∼71%. The ocean is also shallower than the mean depth for Earth. For this reason the time to bring our ancient Mars model ocean and atmosphere into equilibrium is much shorter (approximately hundreds of years) than would be the case for an Earth-like ocean (approximately thousands of years). We assume this equilibrium has been reached when the net radiative balance (the difference between incoming and outgoing fluxes) is less than 0.2 Wm–2 (). Model runs with the dynamic ocean typically require 300 to 400 Martian years to reach equilibrium, although in some cases a checkpoint-restart file from a similar experiment was used to start the simulation in a condition close to a predicted end state. For the latter such simulations may need only 100 y to reach equilibrium. All outputs are averaged over 10 Martian years after radiative equilibrium is attained, and these are used for most of the analysis in this paper. All quantities are expressed in Earth years, noted “y,” and Earth day. When referring about time before present, we used the Earth year ago, noted “a.”

Land Hydrological Cycle

Land surface hydrology (90, 96) is represented by six layers of soil (up to a depth of 3.5 m), a three-layer snow model, and dynamic lakes. The land receives water in the form of precipitation (either rain or snow) and loses it through evaporation/sublimation and via runoff. For the land surface Martian regolith is approximated by sand (with equivalent porosity, hydraulic conductivity, and thermal inertia), initially saturated with 90% liquid water (representing 1 m GEL). The surface runoff depends on the strength of the rain and the level of upper soil layer saturation, while the underground runoff from deeper layers depends on the amount of water in the layer and local slope (in the present simulations the slope was set uniformly to a typical flat desert value). The runoff water goes directly to the lakes, which expand or shrink depending on the amount of available water and exchange water with the lakes in neighboring cells according to a prescribed river routing scheme based on topography. In addition to sublimation, the snow model incorporates the algorithms of melting and meltwater refreezing in deeper layers as well as gravitational compaction of the layers. There is no transport of snow between the cells; i.e., we do not model the glacier dynamics. It is challenging to model glaciers inside a GCM due to their much longer timescales (∼104 to 10 5 y) versus the typical timescales of GCM simulations (∼103 y). In CMIP studies ice sheets are normally calculated offline and incorporated into GCM simulations at different epochs (97). Here we use a simplified approach, also used to model Earth climate, for instance the Last Glacial Maximum, 21,000 y ago (98). We assume that the areas of permanent presence of snow in our simulations represent glacier heads. We assume that in steady state any additional accumulation of snow over these regions is compensated by ice flow to the lower altitudes where it will either melt or eventually reach the ocean. At the typical surface temperatures in our simulations (>250 K), the ice can easily flow on Mars (56, 67, 99). Since we are not modeling the thickness of glaciers, there is no reason to allow snow to accumulate to very high depth. This will require much longer spin-up times but have no noticeable effect on results. So we restrict the snow thickness to 2 m of H2O equivalent (∼3.5 m of compacted snow), which is enough to properly simulate its effect on climate. Any accumulated snow in excess of this amount is immediately removed from the snowpack and is redistributed over the ocean to maintain conservation of water in the system.

Box Modeling

The complex 3D climate system is summarized by box modeling. We divided the surface into three zones: icy highlands, wet lowlands, and ocean, delimited respectively by the mean altitude of the 0 C isotherm and –3,900 m. For each surface box, but also for the atmosphere and the ocean, we compute the sum of all incoming and outgoing fluxes of water, including snow/rain precipitation, evaporation, surface runoff, and simplified glacier return flow (). We use sand for our soil and it is initially saturated with 90% liquid water. This saturation amount will change according to climatic conditions during the course of the simulation (e.g., precipitation or evaporation). The timescale for this process to reach equilibrium can be quite long (approximately thousands of years). Since we are not interested in this process, we estimate and remove this source of water flux in our box model. We check that the hydrological cycle is in balance by examining related diagnostics in the model, but there are unaccounted for fluxes at 5% due to incomplete GCM diagnostics, from salinity changes and numerical rounding. presents the mass fluxes for all model runs.