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 (2015) 176–190

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

A 3 Ga old polythermal ice sheet in Isidis Planitia, Mars: Dynamics and thermal regime inferred from numerical modeling Ondˇrej Souˇcek a,∗ , Olivier Bourgeois b,c , Stéphane Pochat b,c , Thomas Guidat b,c a b c

Charles University in Prague, Faculty of Mathematics and Physics, Mathematical Institute, Sokolovská 83, Praha 8, Karlín, CZ 186 75, Czech Republic Université de Nantes, LPG-Nantes, UMR 6112, 2 chemin de la Houssinière, F-44322 Nantes, France CNRS, LPG-Nantes, UMR 6112, F-44322 Nantes, France

a r t i c l e

i n f o

Article history: Received 15 September 2014 Received in revised form 29 May 2015 Accepted 13 June 2015 Available online 8 July 2015 Editor: C. Sotin Keywords: martian climate thumbprint terrain polythermal ice-sheet basal melting Isidis Planitia liquid water

a b s t r a c t Isidis Planitia is a 1350 km wide impact crater located close to the martian equator. To test the hypothesis that the 2.8 to 3.4 Ga old Thumbprint Terrain preserved on the floor of this basin is a glacial landform assemblage, we perform a numerical simulation of glaciation with a thermo-mechanically coupled model of ice sheet dynamics. As model inputs, we use surface temperatures and ice accumulation patterns predicted by a General Circulation Model based on the present-day atmospheric characteristics, and values of the geothermal heat flux provided by a global model of planetary thermal evolution. We find that, under favorable orbital conditions, an ice sheet covering the entire basin can develop in 2 to 5 Ma, with a maximum thickness of 4.9 km. The modeled ice sheet is polythermal: it is permanently coldbased in the periphery and, due to a negative heat-flux anomaly, also in the center, while the pressure melting point is reached in an intermediate ring. Our simulation is consistent with the interpretation that the Thumbprint Terrain is a martian equivalent of terrestrial ribbed moraines and has formed below a wet-based ice sheet. It supports also the interpretation that sinuous ridges and linear valleys observed at the periphery of the basin are parts of a subglacial network of eskers and tunnel valleys that drained the glacial meltwater outwards, across the cold-based periphery of the ice sheet. This work strengthens the hypothesis that glaciers thick as much as several km may have existed on Mars several Ga ago and that glacial basal melting may have contributed to the production and flow of surface liquid water at that time, under an atmosphere no thicker than the present-day one. © 2015 Elsevier B.V. All rights reserved.

1. Introduction On Mars, currently active glaciers comprise two polar ice sheets (e.g. Fishbaugh and Head, 2001; Greve, 2008) and a number of icefilled craters scattered at northern and southern latitudes higher than 70◦ (Conway et al., 2012). Theoretical considerations on the stability of water ice and climate numerical models, however, predict that surface ice accumulation areas may have shifted repeatedly between all latitudes in the past, in response to changes in martian orbital and atmospheric characteristics (Jakosky and Carr, 1985; Mischna et al., 2003; Levrard et al., 2004; Forget et al., 2006, 2013; Madeleine et al., 2009; Wordsworth et al., 2013). In accordance with these theoretical and numerical inferences, abundant observational evidence of Amazonian (<3 Ga

*

Corresponding author. E-mail address: [email protected] (O. Souˇcek).

http://dx.doi.org/10.1016/j.epsl.2015.06.038 0012-821X/© 2015 Elsevier B.V. All rights reserved.

old) glaciations has been described at all latitudes (Squyres, 1978, 1979; Head and Marchant, 2003; Head et al., 2005, 2006a, 2006b, 2010; Shean et al., 2005, 2007; Shean, 2010; Garvin et al., 2006; Milkovich et al., 2006; Levy et al., 2007, 2010; Dickson et al., 2008, 2010; Holt et al., 2008; Kadish et al., 2008; Morgan et al., 2009; Fassett et al., 2010; Souness et al., 2012; Fastook and Head, 2014; Fastook et al., 2014; Hubbard et al., 2014). There is also a growing body of geomorphic evidence in support of earlier glaciations, during the Hesperian (3–3.7 Ga) and perhaps during the Noachian (>3.7 Ga) (Howard, 1981; Kargel and Strom, 1992; Kargel et al., 1995; Head and Pratt, 2001; Milkovich et al., 2002; Ghatan and Head, 2004; Mège and Bourgeois, 2011; Fastook et al., 2012; Gourronc et al., 2014). Some of these ancient glaciers may have been wet-based and, as such, may have contributed to the production and flow of surface liquid water on Mars (Carr and Head, 2003; Gaidos and Marion, 2003; Head and Marchant, 2003; Fassett and Head, 2006). Assessing the scale, dynamics, thermal regime and historical de-

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

177

Fig. 1. a) Location of Isidis Planitia and Thumbprint Terrain on Mars (modified from Kargel et al., 1995). b) Topographic map of Isidis Planitia and surrounding areas derived from MOLA gridded elevation data.

velopment of these glaciers is thus of importance with regards to the age, lifetime, origin, amount and hydrographic environment of liquid water in the history of Mars. Among the geomorphic evidence presented so far in support of ancient wet-based glaciations, a particularly remarkable feature is the so-called Thumbprint Terrain (Grizzaffi and Schultz, 1989; Lockwood et al., 1992; Kargel et al., 1995). This peculiar landform assemblage consists of parallel sets of periodic arcuate ridges and aligned cones, organized in whorl-shaped patterns and generally associated with other distinctive landforms, such as sinuous ridges, linear depressions, and mounds. The Thumbprint Terrain is mostly located at mid-latitudes, in the Hesperian Northern Lowlands, but the largest field of Thumbprint Terrain identified so far is located near the equator, in the Isidis Planitia impact basin (Fig. 1). Building upon previous works (Grizzaffi and Schultz, 1989; Lockwood et al., 1992; Kargel et al., 1995), Guidat et al. (2015) interpreted the landform assemblage in Isidis Planitia as a glacial landsystem inherited from a polythermal, 3 Ga old, glaciation of the entire basin. This interpretation apparently supports the results of General Circulation Models (GCMs), which predict that surface ice may have accumulated near the martian equator in the past

(Mischna et al., 2003; Levrard et al., 2004; Forget et al., 2006; Madeleine et al., 2009; Wordsworth et al., 2013). The aim of the present study is thus to check whether ice accumulation patterns predicted by these climate simulations are effectively able to produce ice sheets consistent with the development of the observed landform assemblages. For that purpose, we simulate the development of an ice sheet in Isidis Planitia with a thermo-mechanical numerical model of ice sheet dynamics constrained by currently available topographic data and by the outputs of recent climatic and geothermal models. Then, we compare the dynamics and thermal regime of the simulated ice sheet with the organization of the landform assemblage preserved in the basin. The current knowledge on the physiography and geology of Isidis Planitia, including a short description of the landform assemblage on its floor, is reviewed in Section 2. The numerical model, the boundary conditions and the simulation set-up are described in Section 3. The modeling results are presented in Section 4. Their significance with regards to landform interpretation, causes of glaciation, ice sheet extent, glacial flow, basal thermal regime and meltwater production, are discussed in Section 5. Our conclusions are summarized in Section 6.

178

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

Fig. 2. Map of the concentric landform assemblage observed on the floor of Isidis Planitia, with interpretation in terms of a polythermal subglacial landsystem (after Guidat et al., 2015). The Thumbprint Terrain (red) comprises a ring where Arcuate Ridges and Aligned Cones (both interpreted as relicts of ribbed moraines formed below wet-based ice) are organized in whorl-shaped patterns centered on an area located in the NW part of the basin. The absence of Arcuate Ridges and Aligned Cones in the most central and most peripheral parts of the basin suggests that cold-based conditions prevailed in these areas. Sinuous Ridges, Linear Depressions and Mounds (blue) are interpreted as relicts of a peripheral radial network of subglacial channels, which drained the meltwater produced within the interior of the ice sheet across its cold-based periphery. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2. Geological setting 2.1. Physiography Isidis Planitia is a circular topographic depression, about 1350 km in diameter, centered at approximately 13◦ N and 87◦ E (Fig. 1). Its flat floor lies at elevations of −3.6 to −3.9 km and dips gently (0.75 m/km) to the South-West. The basin is surrounded by the heavily craterized Noachian Highlands, which comprise two noticeable mountainous regions: Nili Fossae in the North-West and Lybia Montes in the South-East. In the North-East, Isidis Planitia is connected through a gentle topographic saddle to the Hesperian Martian Lowlands of Utopia Planitia. In the South-West, its border and its floor have been partially covered by the Early Hesperian Syrtis Major volcanic complex (Ivanov and Head, 2003; Hiesinger and Head, 2004; Tanaka et al., 2005). 2.2. Crustal thickness and geothermal heat flux Isidis Planitia is amongst the youngest and best preserved multi-ring impact basins of Mars. It formed approximately 4 Ga ago near the end of the Early Noachian (Wilhelms, 1973; Wichman and Schultz, 1989; Schultz and Frey, 1990; Ivanov et al., 2012). Gravimetry measurements suggest that the impact that formed Isidis Planitia may have thinned the underlying crust to as little as 3 km (Zuber et al., 2000; Neumann et al., 2004). The reduction of radiogenic heat sources, related to this substantial crustal thinning, is responsible for a lower surface heat flux in the center of the basin (12–35 mW/m2 ) than in the surrounding high-

lands (40–60 mW/m2 ), thus making the center of Isidis Planitia a geothermal “cold spot” with respect to its environment (Ritzer and Hauck, 2009; Grott and Breuer, 2010; Ruiz et al., 2011). 2.3. Surface geology and landforms After its formation, the impact crater served as a sink for sediments (Tanaka et al., 2005; Ritzer and Hauck, 2009; Ivanov et al., 2012). Its current surface is believed to be 2.8 to 3.4 Ga old (Tanaka et al., 2005; Ivanov et al., 2012) and exhibits a distinctive landform assemblage, nicknamed Thumbprint Terrain, which comprises Arcuate Ridges, Aligned Cones and Cone Fields, associated with a radial peripheral network of Sinuous Ridges, Linear Depressions, and Mounds (Figs. 2 and 3). Various geological interpretations have been proposed for all these landforms. Volcanic interpretations include cinder or tuff cones (Plescia, 1980; Bridges et al., 2003), rootless cones (Fagents et al., 2002; Bruno et al., 2004), pyroclastic surge deposits (Ghent et al., 2012) and mud volcanoes (Davis and Tanaka, 1995; Skinner and Mazzini, 2009), while glacial interpretations include terminal recessional moraines, ice-pushed moraines, subglacial ribbed moraines, eskers and subglacial tunnel valleys (Grizzaffi and Schultz, 1989; Lockwood et al., 1992; Kargel et al., 1995; Ivanov et al., 2012). A fluvio-glacial origin has also been proposed for the southern part of the peripheral network of Linear Depressions and Sinuous Ridges (Erkeling et al., 2014). Guidat et al. (2015) recognized a distinctive, concentric, organization of the landforms in Isidis Planitia and interpreted this landform assemblage as a glacial landsystem inherited from the pres-

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

179

ence of a polythermal ice sheet over the entire basin, 2.8 to 3.4 Ga ago (Fig. 2). They argued that wet-based conditions prevailed (and led to the formation of Arcuate Ridges and Aligned Cones) over most parts of basin, while the negative geothermal anomaly due to impact-related crustal thinning (Ritzer and Hauck, 2009; Grott and Breuer, 2010; Ruiz et al., 2011) was responsible for coldbased conditions in its most central part. They interpreted the Sinuous Ridges, Linear Depressions and Mounds at the basin margins as relicts of a radial network of subglacial channels, which drained the glacial meltwater produced within the interior of the ice sheet across its cold-based periphery. 3. Numerical model 3.1. Description of the numerical model

Fig. 3. Characteristic landforms observed on the floor of Isidis Planitia (location of panels indicated by black boxes in Fig. 2) (after Guidat et al., 2015). a) Mosaic of visible CTX images (P03 002400 1887 XN 08N276W//P06 003178

1888 XN 08N276W//B07 012276 1887 XN 08N276W//B17 016311 1903 XN 10N276W), illustrating the typical morphology and the spatial organization of Arcuate Ridges. They are 300 to 600 m wide at their base and 100 to 300 m wide at their top, 300 m to 40 km long (4 to 5 km on average) and 20 to 60 m high. They have more or less pronounced arcuate shapes and are organized in periodic networks, with wavelengths of 1 to 10 km. b) Mosaic of visible CTX images (B19 017181 1978 XN 17N269W//B20 017537 1951 XN 15N268W//B21 017893 1968 XN 16N268W) illustrating the typical morphology and the spatial organization of Aligned Cones. They are 300 to 600 m in diameter at their base, 100 to 300 m in diameter at their top and 20 to 60 m high. They form rectilinear to arcuate lineaments, several km to several tens of km in length, and are organized in periodic networks with wavelengths of 1 to 10 km. c) Mosaic of visible HRSC images (h5162 0000 nd3//h5072 0000 nd3) illustrating the spatial organization of landforms in a portion of the peripheral radial network. Sinuous Ridges, Linear Depressions and Mounds are generally connected with each other and weave across Arcuate Ridges and Aligned Cones.

To test the glacial hypothesis for the formation of the Thumbprint Terrain, we employ a numerical thermo-mechanical threedimensional ice-sheet model that simulates a glaciation in Isidis Planitia. The model provides the solution of the system of evolutionary partial-differential equations in spherical geometry, comprising the balances of mass, linear momentum and energy together with a kinematic condition for the evolution of the upper free surface. Ice is described as an incompressible non-Newtonian fluid with power-law viscosity given by Glen’s flow law (e.g. Paterson, 1981). The thermal and mechanical equations are mutually coupled by (i) strong dependence of viscosity on temperature and (ii) stick–slip basal sliding conditions: ice is assumed to be frozen to the bedrock (no-slip) for local basal temperatures below the pressure-melting point, while sliding initiates when melting is reached at the base. All the model equations and the boundary conditions are summarized in Appendix A. The model is based on the Shallow Ice Approximation (SIA) (Hutter, 1983), which exploits the natural scaling properties of large ice-sheets with small vertical–horizontal aspect ratios. While it is possible in the model to improve the SIA solution by an iterative algorithm, described in Souˇcek and Martinec (2008) and successfully tested in Pattyn et al. (2008), only the first step of the algorithm, corresponding to the SIA solution, is used in the presented simulation. The model is numerically implemented by a finite difference technique (Souˇcek, 2010; Souˇcek and Martinec, 2011). The modeling domain is a spherical rectangle (70◦ –110◦ )E × (0◦ –40◦ )N, divided in 161 × 161 grid points. The corresponding horizontal resolution is thus 0.25◦ × 0.25◦ . For the resolution of the ice-sheet in the vertical direction, we employ the “stretched-coordinates”, i.e. 61 layers mapped equidistantly across the thickness of the ice layer, inducing a curvilinear coordinate system, in which the governing equations are solved. The model is run with a time step of 50 years, initiating from ice-free conditions. The prescribed climatic forcing specified in the next section leads to progressive accumulation of ice in areas with positive mass balance. Gradually an ice-sheet develops and starts to flow by a combination of processes of internal deformation and basal sliding until a (quasi-)steady-state is reached.1 The flow process and the character of the steady-state are mainly controlled by the given material properties (rheology) of ice, and by the prescribed boundary conditions, namely by (i) the topography of the bedrock, (ii) the geothermal heat flux, and the climatic input data consisting of (iii) the surface temperature, and (iv) the net surface mass balance. We now discuss these model inputs in more detail.

1 A true steady-state is, in fact, never reached, because the surface temperatures and surface mass balance are non-stationary due to periodic obliquity variations, specified in the next section. Similarly, the basal stick–slip sliding conditions do not allow steady-state in the wet-based regime, due to a development of basal sliding instabilities. By a steady-state, we thus mean a state of the glacier, when the total mass and average thermal regime do no longer change significantly.

180

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

Fig. 4. a) Global map of the geothermal heat flux on Mars 3 Ga ago. The model corresponds to an interpolation between two end-member simulations of Grott and Breuer (2010) and Grott et al. (2011) – a model with a megaplume beneath Tharsis region and a model with a global melt channel. The black square represents our computational domain in Isidis Planitia. b) The interpolated and regridded detail of the geothermal heat flux in our computational domain.

3.2. Model input data We are interested in modeling glaciation on the surface of Mars approximately 3 Ga ago (the estimated age of the landforms in Isidis Planitia after Ivanov et al., 2012), facing therefore rather poorly constrained input boundary conditions for the numerical model due to a lack of direct observational evidence. Consequently, the climatic and geothermal inputs for our simulation are derived from the results of other numerical simulations while the topographic input is based on the present-day topography. 3.2.1. Bedrock topography Adopting the assumption from Guidat et al. (2015), that the presently observable landforms on the floor of Isidis Planitia are subglacial remnants, we assume that the surface topography of the Isidis region has not changed substantially over the past 3 Ga. We assume the bedrock topography to be fixed throughout the simulation, ignoring thus a possible deformation of the lithosphere due to isostatic compensation of the (modeled) ice sheet. The numerical discretization of the bedrock topography is based on the Mars Orbiter Laser Altimeter (MOLA) data set, smoothed by a Gaussian convolution filter and regridded to a resolution of 0.25◦ × 0.25◦ . 3.2.2. Geothermal heat flux The geothermal heat flux input data are derived from the global model of Grott et al. (2011), 3 Ga ago, corresponding to the estimated age of glacial landforms in Isidis Planitia (Ivanov et al., 2012; Guidat et al., 2015). Grott and Breuer (2010) proposed two scenarios for the global thermal evolution of Mars: one with a plume below Tharsis, the other with a global melt channel. These scenarios provide different values for the geothermal heat flux in the region of Isidis. Both have similar spatial patterns, but differ by a constant shift, the spread between the scenarios being approximately 13 mW m−2 . We pick a model in between the two scenarios; the corresponding global heat flux and the zoom to the computational domain are depicted in Fig. 4. The heat flux value is on the order of 35 mW m−2 in Isidis Planitia, while it reaches values higher than 44 mW m−2 in the surrounding highlands. This negative geothermal anomaly results from the crustal thinning below the basin (Zuber et al., 2000; Neumann et al., 2004) – the missing crust material representing a “negative source” of radiogenic heating, with respect to the surroundings. 3.2.3. Climatic input data – surface temperatures and mass balance It is now widely accepted that the climate of Mars may have changed significantly throughout its history, as a result of the variations of several parameters. These possibly include variations in

the activity of the Sun (Gough, 1981), variations in the orbital parameters of the planet (Laskar et al., 2004), and variations in the thickness, composition and dust content of the atmosphere (Kasting, 1991; Haberle, 1998). Many numerical models have been developed (see e.g. reviews in Haberle, 1998; Forget et al., 2013) to constrain the effects of these parameters on the history of the climate and, more specifically, on the surface temperature, temporal and geographic location and form (liquid or solid) of atmospheric precipitation onto the surface. Some GCMs, based on the current atmospheric properties, have been designed to explore dominantly the role of orbital variations in recent times (e.g. Mischna et al., 2003; Madeleine et al., 2009), while others have been designed to explore additionally the role of a fainter young Sun and a different atmosphere earlier in the martian history (e.g. Forget et al., 2013; Wordsworth et al., 2013). All these models consistently predict annual mean surface temperatures below the melting point of water at all latitudes and atmospheric precipitation in the form of solid ice rather than liquid water (Mischna et al., 2003; Madeleine et al., 2009; Forget et al., 2013; Wordsworth et al., 2013). These models suggest that, under a wide range of orbital, surface and atmospheric conditions, net surface ice accumulation may have occurred in the past near the equator and, more specifically, in the region of Isidis Planitia (Mischna et al., 2003; Madeleine et al., 2009; Forget et al., 2013; Wordsworth et al., 2013). In these models, the location of net surface ice accumulation at a given time is highly sensitive, however, to the orbital configuration of the planet (obliquity, eccentricity, latitude of the Sun at perihelion), atmospheric dust content, surface topography, physical properties of the surface materials, and initial geographic location of surface and subsurface volatile reservoirs. To check whether the landform assemblage observed in Isidis Planitia is consistent with the development of an ice sheet in response to such surface ice accumulation, we use the results of the Laboratoire de Météorologie Dynamique (LMD) GCM as input data for the simulation of ice sheet dynamics. Amongst all orbital configurations tested with the LMD GCM, we selected one that favors ice accumulation at low latitudes: the reference simulation of Madeleine et al. (2009). This simulation is based on atmospheric properties similar to the current ones, but with different orbital properties (high obliquity 0 = 35◦ , high eccentricity e = 0.1, aerocentric longitude of the Sun at perihelion Lp = 270◦ ), an initial ice reservoir located on the flanks of the Tharsis volcanoes and a large amount of dust in the atmosphere (dust opacity τdust = 2.5).2 2 The high amount of dust in the atmosphere portrays the effect of frequent dust storms (Haberle et al., 2003) and huge amounts of dust lifting (Newman et al., 2005) at high obliquities (Madeleine et al., 2009).

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

181

Fig. 5. a) Global map of the surface temperature on Mars resulting from the Mars LMD GCM reference simulation of Madeleine et al. (2009) and b) Interpolated and regridded detail of the surface temperature map in our computational domain.

Fig. 6. a) Global map of net surface accumulation rate (annual mass balance), as computed in the reference simulation of Madeleine et al. (2009). The black square represents our computational domain in Isidis Planitia. b) Interpolated and regridded detail of the map of annual surface mass balance in our computational domain.

3.2.4. Surface temperature The global surface temperature map corresponding to the reference simulation in Madeleine et al. (2009) is depicted in Fig. 5 together with a zoom to the computational domain. Because of the relatively large time span of our simulation (5 Ma), it is necessary to take into account changes in the insolation due to variations of martian orbital parameters. As a first approximation, we employ a simplified approach after Greve (2004). Martian obliquity evolution is parameterized by a two-periodic cycle: the mean obliquity value 0 = 35◦ is varied by a short-periodic signal (T main = 125 ka), modulated by a slower cycle (Tmodul =  1.3 Ma). It is given by a formula  (t ) = 0 + ˆ (t ) sin





2π t T main

, where

π t , with ˆ ◦ and ˆ (t ) = ˆmax +2 ˆmin − ˆmax −2 ˆmin cos T 2modul max = 10 ◦ ˆmin = 2.5 . Consequently, the surface temperature data set from

Madeleine et al. (2009) (Fig. 5b) is modulated in the simulation by a time-dependent scaling factor ∼



cos  (t ) cos 0

1 4

, corresponding to an

equilibrium relation between the (approximately) equatorial insolation for the given obliquity and the black-body radiation (Greve, 2004). 3.2.5. Surface mass balance The net surface mass balance is derived from the annual averages of accumulation and sublimation rates, as computed in the reference simulation of Madeleine et al. (2009) (Fig. 6a). This coarse-grid model, interpolated to our model resolution yields positive balance (net annual accumulation, with values as high as 12 kg m−2 a−1 ) along the northwestern margin and over the southern part of Isidis Planitia, whereas negative balance (net annual

ablation) occurs over the northern part of the basin, the southern highlands and the southern part of Utopia Planitia (Fig. 6b). To account for climate changes due to obliquity variations during the simulation, the reference annual balance is varied in time by a temperature-dependent scaling factor, as in Greve 1 1 (2004), ∼ exp( Rλ ( sur )) with λ = 2860 kJ kg−1 and R m = f − sur f m

T0

T

(t )

461.5 J kg−1 K−1 . This factor expresses the temperature dependence of vapor saturation pressure, which is the key factor controlling the kinetics of both the accumulation and sublimation rates. This parameterization of the sublimation rates in particular is rather simplified, since at least two sublimation mechanisms act – buoyancy- and turbulence-driven – and their temperature dependence is more complex than the one proposed (Fastook et al., 2008). While the scaling factor which we use is common to both the accumulation and sublimation rates, and the temperature variations thus cannot change the sign of the original net balance, a more realistic parameterization (Fastook et al., 2008) should lead to slightly asymmetric temperature dependence of sublimation vs. accumulation rates leading to excessive sublimation at very low temperatures. The proposed picture of surface temperatures and surface mass balance is simplified, because it completely ignores the feedback between the glaciation and the atmospheric processes. The development of an ice-sheet in the Isidis Planitia region, and the associated changes in the surface altitude and albedo, would surely have a strong impact on the temperature, atmospheric circulation, and consequently surface mass balance pattern in the region. It may be tempting to partially parameterize these effects e.g. by adopting altitude-dependence of both the accumulation and sublimation rates and specifying adiabatic lapse rate, but however elaborate,

182

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

Fig. 7. a)–f): Snapshots of the evolution of the ice sheet topography and surface flow during the simulation. Colors are used for coding the ice thickness, stippled black contours indicate the elevation of the ice sheet surface, plain black contours indicate the elevation of the bedrock, and arrows show the surface flow velocity. g) Topographic profiles of the simulated ice sheet at different time steps (vertical exaggeration: ×70). Wet-based regions at 5 Ma indicated by red lines at the base of the ice sheet. Location of the profiles indicated by a black line in the inset map. h) Quasi-steady-state (5 Ma) topography of the ice sheet along the same profile as in g) with a smaller vertical exaggeration (×10). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

such parameterizations might still be missing the essential points of the glacier-climate interaction. For example, we would expect that the observed net accumulation pattern on the western margin of Isidis Planitia in the reference simulation of Madeleine et al. (2009) is due to the release of water moisture carried by the north-eastward winds when they hit this topographic barrier. With an ice-sheet growing gradually in the basin, this topographic barrier would be migrating eastwards shifting the accumulation maximum accordingly. None of the parameterizations mentioned above allows to capture such effects without introducing a plethora of additional ad-hoc assumptions and poorly-constrained parameters. This considered, we restrict ourselves to the simple quasi-static cli-

mate setting described above, until we are actually able to perform truly coupled numerical simulation of both the glaciation process and martian climate evolution. 4. Results 4.1. Ice sheet topography and surface flow The evolution of the ice sheet topography and surface flow from the onset of glaciation to the quasi-steady-state and also the corresponding ice-sheet cross-sections along a NW-SE profile are depicted in Fig. 7 (see also videos in the supplementary ma-

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

183

Fig. 8. Evolution of a) – the ice-sheet volume b) – glaciated area. Note the small oscillations superimposed on the main curves due to the imposed two-cycle orbital variations.

Fig. 9. Snapshots of the evolution of basal temperature (relative to the melting point) during the simulation. Red regions correspond to the areas, where the basal ice is at the melting point. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

terial). The ice sheet initiates at the northwestern border of the basin (Fig. 7a), where the net accumulation rates are the highest (Fig. 6b). Then, the ice flows approximately radially from this area and fills progressively the basin as the ice sheet grows (Fig. 7b–f). The entire basin is covered with ice after about 2 Ma. The evolution of ice-sheet volume and area are depicted in Fig. 8. For the first approximately 1 Ma, the ice-sheet volume grows almost linearly (note the sinusoidal variations due to periodic obliquity changes). After this period the growth rate gradually decreases as the ice-sheet flow becomes effective in transporting the ice into zones with negative mass balance, until finally, after about 5 Ma, a quasi-steady-state is reached and the surface topography, the flow lines, the margins and the volume of the ice sheet do not change significantly. At this stage, the thickest part of the ice sheet (ca. 4.9 km) is located in the northwestern part of the basin. Its southwestern, southern and southeastern margins are controlled by the basin topographic borders, while its northeastern part overflows the gentle topographic saddle between Isidis Planitia and Utopia Planitia. In

the North-West, the ice sheet covers the mountains of the Nili Fossae region with as much as 3 km of ice. The surface ice flow lines display a fan-shaped pattern centered on an area located in the northwest: the ice flows southwards along the western border of the basin, southeastwards in its central part and northeastwards (towards Utopia Planitia) in its northeastern part. Surface flow velocities are in the range of 0.5 to 3 ma−1 in most parts of the ice sheet, except for a patch of noticeably slower motion located approximately in the center of the basin and, in contrast, regions of faster flowing ice (up to 8.6 ma−1 ) in the northern region where the ice flow is not hindered by a topographic border. 4.2. Basal thermal regime The evolution of the temperature at the base of the ice sheet from the onset of glaciation to the quasi-steady-state is depicted in Fig. 9 (see also videos in the supplementary material). During the earliest stages of the ice sheet growth, the basal temperature remains below the melting point everywhere. After about 1 Ma, the

184

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

melting temperature is reached in a crescentic region that fringes the western border of the basin, where the ice sheet is thickest. Afterwards, both the northern and southern tips of this crescentic area grow in extent, until the wet-based region eventually forms a ring when the quasi-steady state is reached. This wet-based ring surrounds a permanently frozen circular patch located in the central part of the basin. At the periphery of the ice sheet is a narrow ring of permanently frozen ice that increases in width over the Nili Fossae mountains and in Utopia Planitia. The base of the ice sheet does not remain at the melting temperature neither permanently nor everywhere in the wet-based ring. Rather, in the simulation, this ring is characterized by a banded pattern of “melting waves” – patches of wet-based ice in the cold-based surroundings – that propagate from the northwestern part of the basin, along both the northern and southwestern branches of the wet-based ring (see the videos in the supplementary material). These melting waves have arcuate shapes and are approximately orthogonal to the borders of the wet-based ring. 5. Discussion 5.1. Extent, thickness and volume of ice sheet The simulation produces an ice sheet that covers the entire floor of Isidis Planitia and the mountains around Nili Fossae. This result is consistent with the interpretation that the landform assemblage on the floor of Isidis Planitia is a fossil subglacial landsystem inherited from the former presence of a basinwide ice sheet (Grizzaffi and Schultz, 1989; Lockwood et al., 1992; Kargel et al., 1995; Guidat et al., 2015). By contrast, to our knowledge, a glacial interpretation of the landforms in the mountains around Nili Fossae has never been proposed. On the currently available orbital images, we were not able to find unambiguous geomorphic evidence of glaciation in these mountains. The fact that the simulated ice sheet is cold-based in this region (Fig. 9) might explain this lack of geomorphic evidence. While the ice sheet grows during the simulation, its southern margin remains pinned at the boundary between the floor of Isidis Planitia and the highlands of Lybia Montes and Syrtis Major. This constant position suggests that the location of the margin is controlled mostly by the steep topographic border of the basin in this area. By contrast, in its northeastern part, the simulated ice sheet progressively overlaps the topographic saddle that borders Isidis Planitia, and eventually spreads into the southern part of Utopia Planitia. The location of the ice sheet margin in this area of minor topographic relief is thus controlled by complex relations between the patterns of accumulation, ablation and ice flow on a regional scale. Therefore, the position of the northeastern margin of the simulated ice sheet must be taken with caution. It should be noted, in particular, that our simulation takes into account ice accumulation and ablation patterns predicted by Madeleine et al. (2009)’s reference simulation in our modeling domain only. In this domain, the major area of net accumulation is located between the northwestern border of Isidis Planitia and the mountains of Nili Fossae (Fig. 6b). However, Madeleine et al. (2009)’s reference simulation also predicts the existence, under the same orbital and atmospheric conditions, of major ice accumulation areas further north, in Utopia Planitia and along the martian topographic dichotomy between Nilosyrtis Mensae and Deuteronilus Mensae (Figs. 1 and 6). It has been demonstrated that ice accumulation in these northern areas can lead to the development of glaciers extending down to the martian northern lowlands (Fastook et al., 2011). It is therefore to be expected that, if an ice sheet derived from the accumulation area in the Nili Fossae mountains was once able to cover Isidis Planitia, another ice sheet derived from these northern accumulation areas was also able to cover

large portions of Utopia Planitia at the same time. This hypothesis is consistent with the presence of the Thumbprint Terrain in Utopia Planitia (Lockwood et al., 1992; Kargel et al., 1995; Ivanov et al., 2012, 2014) (Fig. 1). If this happened at some time, both ice sheets might have merged across the topographic saddle between Isidis Planitia and Utopia Planitia. Therefore, the location of the ice sheet northern margin predicted by our simulation is perhaps unrealistic. The volume of the simulated ice sheet is approximately 6 × 106 km3 when the entire basin is covered with ice at 2 Ma and 7 × 106 km3 when the quasi-steady state is reached at 5 Ma. For comparison, the Antarctic and Greenland ice sheets have approximate volumes of 26 × 106 and 3 × 106 km3 respectively (Marshall, 2005), while the present-day Martian North and South Polar Caps have approximate volumes of 1.5 × 106 and 2.5 × 106 km3 respectively (Smith et al., 1999). Volumes on the order of 106 to 107 km3 have been proposed for paleoglaciers on the Tharsis dome, in Valles Marineris, in Dorsa Argentea and in Deuteronilus Mensae (Kite and Hindmarsh, 2007; Fastook et al., 2008, 2011, 2012; Gourronc et al., 2014). The ice sheet simulated here in Isidis Planitia thus stands in the range of volumes suggested for other martian paleo-glaciers, but it is nearly twice as large as the total volume of ice trapped in the present-day polar caps. This difference is consistent with the hypothesis that a fraction of the total ice volume that was available at the surface of Mars in the past has been lost to space, and that another fraction currently resides in the subsurface (Carr and Head, 2015) in the form of ground ice (e.g. Feldman et al., 2004; Mangold et al., 2004; Byrne et al., 2009; Mellon et al., 2009; Mouginot et al., 2010) and sequestered glacial ice (e.g. Gourronc et al., 2014; Levy et al., 2014). 5.2. Direction of ice flow Based on the hypothesis that Arcuate Ridges and Aligned Cones have formed orthogonal to ice flow, Guidat et al. (2015) tentatively reconstructed the pattern of ice flow in Isidis Planitia (Fig. 10c). In their reconstruction, the ice flow lines define a radial pattern centered about an area located in the North-West and they are deflected around a central patch of immobile cold-based ice. They attributed this deflection to the negative geothermal anomaly located below the basin center. The overall flow pattern obtained in our simulation is similar to the reconstruction proposed by Guidat et al. (2015), except in two noticeable regions: (1) in the NW part of the basin, the ice flows eastwards (from the major accumulation area, located close to the mountains of Nili Fossae) in the simulation, while it flows westwards in the reconstruction; (2) in the eastern part of the basin, the ice flows to the northeast in the simulation, while it flows to the southeast in the reconstruction. These local differences suggest that the relations between the direction of ice flow and the direction of landforms might be more complex than those envisaged by Guidat et al. (2015). The deflection of the flow lines around the central patch of cold-based ice is, however, apparent in the simulation, this effect being most pronounced during the initial glaciation phase. A similar control of the subglacial basal thermal regime on ice flow lines and velocities has been described also in terrestrial glaciers (Bourgeois et al., 2000; Kleman and Glasser, 2007). 5.3. Basal thermal regime In the simulation, the central part of the ice sheet and its periphery remain permanently cold-based, while basal melting occurs episodically in the intermediate ring (Figs. 9, 10b). This concentric organization of the basal thermal regime is consistent with the organization of the landforms described by Guidat et al. (2015)

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

185

Fig. 10. Comparison between the reconstructed and the simulated ice sheet dynamics. a) Map of the landform assemblage on the floor of Isidis Planitia, with interpretation in terms of basal thermal regime (after Guidat et al., 2015). b) Map of modeled basal temperature at the quasi-steady-state. Red regions correspond to the areas, where the basal ice is at the melting point. c) Map of ice flow directions derived from the orientation of Arcuate Ridges and Aligned Cones (after Guidat et al., 2015). Flow arrows are red when their senses or directions diverge significantly from those of the modeled velocity vectors. d) Map of the modeled ice flow velocity and surface topography at the quasi-steady-state. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

and with their interpretation that the presence of Arcuate Ridges and Aligned Cones reflects wet-based conditions over large portions of the basin, while their absence at the basin center and periphery reflects cold-based conditions in these areas. The concentric organization of the basal thermal regime obtained in the simulation can be understood in terms of a competition between two processes. First, the basal temperature generally increases with the thickness of the ice sheet from its periphery toward its center, because of the increasing insulating effect of the ice cover. This effect produces a transition from cold-based conditions at the periphery to wet-based conditions in the interior of the ice sheet. This general trend is altered, however, by the existence of a negative geothermal anomaly below the basin: the basal temperature decreases as the geothermal heat flux decreases towards the basin center. This second effect is involved in the development of the cold-based central patch. Uncertainties in the climatic and geothermal boundary conditions are thus critical with respect to the location and width of the simulated wet-based ring. The specific configuration of climatic and geothermal parameters used in the present simulation allows the development of a

wet-based ring that is narrower and laterally offset with respect to the wet-based ring inferred by Guidat et al. (2015). We suspect a better fit could be obtained with slightly different climatic and geothermal inputs for the simulation. In the simulation, a banded pattern of “melting waves” propagates in the wet-based ring (see videos in the supplementary material). They represent variable basal patches of wet-based ice in the cold-based surroundings and they propagate approximately downstream along the flowlines. We suspect that the formation of these melting waves may be associated with the development of either internal creep or basal stick–slip instabilities. In both cases, the instability would arise from the following thermo-mechanical feedback mechanism: increased temperature promotes increased deformation/sliding due to a lowered viscosity/friction. Enhanced deformation/sliding in turn increases the dissipative/frictional heating, which increases the temperature even further. Such process typically leads to a speed-up of the deformation/sliding rate up to a level, when the advection of cold ice or the flow-induced change of the ice-sheet topography cool the system down below the activation threshold for the instability.

186

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

Under favorable conditions, the whole process may repeat itself (quasi)-periodically. Such thermo-viscous instabilities have been proposed to account for rapid changes in the flow and volume of terrestrial ice sheets (e.g. MacAyeal, 1992; Payne, 1995) and have been suggested as a possible mechanism for the development of ice streams (Payne and Dongelmans, 1997). It has been reported, nevertheless, that numerical models based on the shallow-ice approximation cannot reliably resolve the dynamics of such instabilities (Hindmarsh, 2004, 2006, 2011) and that the inclusion of membrane stresses in the momentum balance is essential for the robustness of their simulation (Hindmarsh, 2009). On the other hand, Saito et al. (2006) speculate that despite being obscured by certain numerical artifacts, the thermo-mechanical instabilities observed in the SIA-based simulations may still represent a real physical process. Consequently, the “melting waves” in our simulation should be interpreted with caution, in particular as far as their quantitative characteristics are concerned. We hypothesize, however, that the area of the wet-based intermediate ring where the instability appears in our simulation, could have been prone to development of real physical instabilities by the thermo-mechanical coupling mechanism described above. In addition, if we accept the simulated melting waves as real rather than numerical artifacts, it is remarkable to notice that they might be consistent with the analogy proposed by Guidat et al. (2015) between the Thumbprint Terrain in Isidis Planitia and the fields of ribbed moraines on Earth. The formation mechanism of ribbed moraines is still poorly understood, but it probably involves transitions from cold-based to wet-based conditions, associated with a stick–slip behavior of glaciers (Hättestrand and Kleman, 1999; Dunlop and Clark, 2006; Chapwanya et al., 2011). The Thumbprint Terrain on Mars might thus constitute a geomorphic evidence for stick–slip glacial motion in response to the development of englacial or subglacial thermomechanical instabilities within or beneath ancient ice sheets. 5.4. Relation between glaciation, climate, orbital cycles and atmospheric physical properties The simulation presented here demonstrates that, under a climate corresponding to a specific combination of external and internal parameters at least (high obliquity, high eccentricity, favorable longitude of the Sun at perihelion, high amount of dust in the atmosphere, existence of an ice reservoir in the Tharsis region, existence of a negative geothermal anomaly beneath the basin), a polythermal ice sheet can develop in Isidis Planitia under an atmosphere no thicker than the current one. The 2.8 to 3.4 Ga old landform assemblage preserved on the basin floor is consistent with the dynamics and basal thermal regime of the simulated ice sheet. A thicker atmosphere is thus not needed neither to explain the glaciation nor the presence of liquid water in the basin at that time. Several Ma are required to achieve the full glaciation of the basin in the simulation. Such a glaciation would thus require that conditions favorable to ice accumulation persisted during far more than just one martian obliquity cycle (0.1 Ma at the present-day). The martian eccentricity cycle, with its present-day duration of 2 Ma, might therefore have played a major role in sustaining this glaciation. Another possibility is the persistence of an anomalously high amount of dust in the atmosphere for several Ma. Both options are consistent with the inference by Madeleine et al. (2009) that high eccentricity values and high atmospheric dust contents are necessary for the development of large mid-latitude glaciations. If the glacial interpretation of the landforms is correct, their age suggests that a combination of parameters favorable to the de-

velopment of a basin-wide ice sheet in Isidis Planitia was never met again during the last 2.8 Ga. 5.5. Relation between glaciation and production of surface liquid water Numerous landforms on Mars, especially in the region of Isidis Planitia, have been attributed to erosion or deposition by liquid water (Baker, 2001; Ivanov et al., 2012; Erkeling et al., 2012, 2014). The amount, lifetime, age, origin (atmospheric precipitation, melting of surface or underground ice) and hydrographic environment (rivers, lakes, oceans, underground networks, subglacial networks, . . . ) of this liquid water are matters of intense debates (Baker, 2001). Our modeling results support the hypothesis that basal melting below the wet-based portions of an ice sheet was able to produce liquid water on the floor of Isidis Planitia. They also support the hypothesis by Guidat et al. (2015) that Sinuous Ridges, Linear Depressions and Mounds observed at the periphery of the basin are part of a radial network of channels that drained the meltwater outward, across the cold-based periphery of the ice sheet (Fig. 10a). Glacial meltwater is thus a reasonable candidate for the formation of landforms attributed to erosion or deposition by liquid water in and around Isidis Planitia. GCMs based on the present-day composition of the atmosphere (Mischna et al., 2003; Madeleine et al., 2009) or on other physical conditions (fainter young Sun, thicker atmosphere) (Forget et al., 2013; Wordsworth et al., 2013) predict the existence of widespread ice accumulation areas at many other places on Mars. The present analysis of a regional glaciation related to one of these accumulation areas thus supports the hypothesis that, on a global scale, subglacial melting may have contributed to the production and flow of liquid water on Mars (Kargel and Strom, 1992; Kargel et al., 1995; Carr and Head, 2003; Fastook et al., 2012; Gourronc et al., 2014). 6. Conclusions The present modeling study demonstrates that, with atmospheric physical properties similar to the present-day ones and under favorable orbital conditions (obliquity ∼35◦ , eccentricity ∼0.1, aerocentric longitude of the Sun at perihelion ∼270◦ , high amount of dust in the atmosphere), ice accumulation in the northwestern part of Isidis Planitia (Madeleine et al., 2009) may lead within 2 to 5 Ma to the development of an ice sheet with a maximal thickness of 4.9 km over the entire impact basin. The modeled ice sheet is polythermal: its central part and its periphery are permanently cold-based, while the pressure melting point is reached in an intermediate ring. This concentric organization of the basal thermal regime is related to the existence of a negative geothermal anomaly, due to crustal thinning below the basin center (Ritzer and Hauck, 2009; Grott and Breuer, 2010). These modeling results are consistent with the interpretation that the 2.8 to 3.4 Ga old landform assemblage preserved on the floor of Isidis Planitia is a glacial landsystem, with the Thumbprint Terrain corresponding to the wet-based portions of the former ice sheet (Guidat et al., 2015). This work strengthens the hypothesis that glaciers as thick as several km may have existed on Mars several Ga ago and that basal melting below the wet-based portions of these ice sheets may have contributed to the production and flow of liquid water at that time (Kargel and Strom, 1992; Carr and Head, 2003). A thicker atmosphere is not needed to explain neither the development of such ice sheets nor the presence of liquid water at their bases. Acknowledgements We would like to express our deep gratitude to Jean-Baptiste Madeleine for most kindly providing us with the surface tem-

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

peratures and surface mass balance obtained from the LMD GCM (Madeleine et al., 2009), and likewise to Matthias Grott and his coworkers for sharing with us the heat flux data from their (Grott et al., 2011) simulations. We are extremely grateful to Bryn Hubbard and an anonymous reviewer for their thorough and detailed reviews, which helped us remove a number of weaknesses from the initial manuscript and strengthen the argumentation in the present version. O.S. acknowledges the support by project LL1202 in the program ERC-CZ funded by the Ministry of Education, Youth and Sports of the Czech Republic. Part of this work was performed during a research stay of O.S. at the LPG-Nantes, funded by the Université de Nantes. We also acknowledge financial support from the Observatoire des Sciences de l’Univers Nantes Atlantique (OSUNA), through the “Analyse comparative des systèmes morphologiques glaciaires et paraglaciaires islandais et martiens” project, from the Centre National de la Recherche Scientifique (CNRS), through the INSU SYSTER “Analyse Géomorphologique et Modélisation des Réseaux d’Eau de Fonte Sous-Glaciaires” project, and the PICS “Planetary Interiors constrained by Various Observations” project.

rheology (4 deformation mechanisms even for pure-water ice, cf. Durham et al., 2001; Goldsby and Kohlstedt, 2001), or the possible weakening or strengthening effects due to likely presence of contaminants such as dust, CO2 ice or CO2 clathrate-hydrates (e.g. Durham et al., 1997), salt (Fisher et al., 2010), or effects related to grain-size evolution (Kubo et al., 2006). Presence of contaminants would naturally also influence other thermodynamic properties such as the pressure melting point or the thermal conductivity, heat capacity, density, etc. We ignore these effects in this study, being interested mainly in the leading-order features of large-scale ice-sheet dynamics, without any ambitions of resolving fine-scale features and subtleties by the model. A.2. Boundary conditions

• Upper free surface is explicitly parameterized by an elevation function r = f s (ϑ, φ, t ), where ϑ, φ are the colatitude and longitude, respectively and r is the distance from the origin. The evolution of the free surface is governed by the kinematic equation

vφ ∂ f s ∂ f s vϑ ∂ f s + + − vr = as , ∂t r ∂ϑ r sin ϑ ∂φ

Appendix A. Model equations A.1. Field equations – balance equations Mathematically the problem is formulated as the Stokes– Fourier system for an incompressible fluid in Eulerian description:

div v = 0

(A.1)

0 = −∇ p + div σ + ρ g

(A.2)

˙ = div(k()∇) + σ : d , ρ c v ()

(A.3)

where v denotes the ice flow velocity, p is the pressure, ρ is the ice density, g is the gravity acceleration, c v () is the heat capacity,  is the absolute temperature, k() denotes the heat conductivity of ice, d is the symmetric part of the velocity gradient d = 12 (∇ v +

˙ is the material time derivative () ˙ := ∇ T v), ()

∂() ∂t

+ v · ∇(), σ : d =

σ i j di j 3 and σ is the deviatoric part of the Cauchy stress tensor given by the Glen’s flow law (e.g. Paterson, 1981):

1

σ = 2ηd, η = A(, p )−1/n d(II1−n)/n , n = 3 ,

(A.4)

2

(− pI + σ ) · ns = − p atm ns , atm

σ,

(A.5)



where the second invariants dII , σII , are given by dII = di j di j /2,  σII = σ i j σ i j /2, and the Arrhenius rate factor reads A(, p ) =





A exp − RQ with R the universal gas constant, Q the activation energy and  the absolute temperature corrected for the pressure melting point,  =  + C Cl p, C Cl being a Clausius–Clapeyron slope. The temperature field is artificially cut off above the melting temperature throughout the simulations, and melt production is only computed at the base of the wet-based zones and melt transport is not resolved in the model. Wet-based regions are merely recognized as regions where the temperature equals the local melting temperature. It must be noted that the material model assuming pure-water ice composition and only one deformation mechanism (one viscosity) represents also a rather crude, though widely used, simplification. Consequently, in this study we are thus neglecting possible effects related to the complexities of ice

(A.7) s

p is the atmospheric pressure and n is the outer unit normal to the glacier at the upper surface. At the upper surface we also prescribe the surface temperature

 = s .

(A.8)

• Glacier bed is parameterized by a bed elevation function r = f b (ϑ, φ, t ) which in this application is assumed to be fixed, i.e. we do not solve for isostatic adjustment by viscoelastic deformation of the lithosphere. Boundary condition for the glacier movement differs according to the basal temperature, switching between the frozen-bed (no-slip) conditions and basal sliding:



d = A(, p )σ

v=

0

 < M , (A.9)

− ρC Rg |τ b2| τ b  =  M , 2

Nb

where τ b denotes the tangent and N b the normal component of the traction vector at the glacier base, respectively, and C R is the drag coefficient and  M = 0M − C Cl p is the melting temperature. Concerning the energy equation, at the glacier base, we prescribe either the normal component of the heat flux vector q provided the base is below the melting point, or, for a base at the melting point, we evaluate from the jump of the heat fluxes and the frictional heating the corresponding melting rate:



q · nb = −q geo q · nb + q geo + v · τ b = Lab

 ≤ M ,  = M ,

(A.10)

where q geo is the geothermal heat flux, nb is the outer unit normal to the glacier at the base, L is the latent heat of melting of ice and ab is the rate of melting at the base. A.3. Material parameters

3

Einstein’s summation convention is used for repeating indices throughout the paper.

(A.6)

where as is net surface mass balance (accumulation-ablation function) representing both the ice accumulation by ice depositing in accumulation areas (a s ≥ 0) and ablation in melting/sublimation regions (a s ≤ 0). The surface is assumed to be free, which implies

or, inversely, n −1 II

187

All material parameters are listed in Table 1.

188

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

Table 1 Model parameters. Symbol

Variable

Value

Unit

ρ

ice density surface gravity heat capacity of ice

910 3.72

kg m−3 m s−2 J (kg K)−1 J kg−1 K−2 J (kg K)−1 W (m K)−1 W (m K)−1 K−1 s−1 Pa−3 s−1 Pa−3 kJ mol−1 kJ mol−1 a−1 K m−1 J mol−1 K−1 – K

g c v () = c 1  + c 2 c1 c2 k() = k1 exp(−k2 ) k1 k2 A

7.253 146.3 ice thermal conductivity

Pre-exp factor

Q

activation energy

CR C Cl R n

Basal drag coefficient Clausius–Clapeyron gradient Universal gas constant Glen’s flow law exponent Melting point at zero pressure

0M

Appendix B. Supplementary material Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.epsl.2015.06.038. References Baker, V.R., 2001. Water and the Martian landscape. Nature 412, 228–236. Bourgeois, O., Dauteuil, O., Van Vliet-Lanoe, B., 2000. Geothermal control on flow patterns in the Last Glacial Maximum ice sheet of Iceland. Earth Surf. Process. Landf. 24, 1–18. Bridges, J.C., Seabrook, A.M., Rothery, D.A., Kim, J.R., Pillinger, C.T., Sims, M.R., Golombek, M.P., Duxbury, T., Head, J.W., Haldemann, A.F.C., Mitchell, K.L., Muller, J.P., Lewis, S.R., Moncrieff, C., Wright, I.P., Grady, M.M., Morley, J.G., 2003. Selection of the landing site in Isidis Planitia of Mars probe Beagle 2. J. Geophys. Res. 108, 1–17. Bruno, B.C., Fagents, S.A., Thordarson, T., Baloga, S.M., Pilger, E., 2004. Clustering within rootless cone groups on Iceland and Mars: effect of nonrandom processes. J. Geophys. Res. 109, 1–11. Byrne, S., Dundas, C.M., Kennedy, M.R., Mellon, M.T., McEwen, A.S., Cull, S.C., Daubar, I.J., Shean, D.E., Seelos, K.D., Murchie, S.L., Cantor, B.A., Arvidson, R.E., Edgett, K.S., Reufer, A., Thomas, N., Harrison, T.N., Posiolova, L.V., Seelos, F.P., 2009. Distribution of mid-latitude ground ice on Mars from new impact craters. Science 325 (5948), 1674–1676. Carr, M.H., Head, J.W., 2003. Basal melting of snow on early Mars: a possible origin of some valley networks. Geophys. Res. Lett. 30, 2245. Carr, M.H., Head, J.W., 2015. Martian surface/near-surface water inventory: sources, sinks, and changes with time. Geophys. Res. Lett. 42. http://dx.doi.org/10.1002/ 2014GL062464. Chapwanya, M., Clark, C.D., Fowler, A.C., 2011. Numerical computations of a theoretical model of ribbed moraine formation. Earth Surf. Process. Landf. 36, 1105–1112. Conway, S.J., Hovius, N., Barnie, T., Besserer, J., Le Mouélic, S., Orosei, R., Read, N.A., 2012. Climate-driven deposition of water ice and the formation of mounds in craters in Mars’ north polar region. Icarus 220, 174–193. Davis, P.A., Tanaka, K.L., 1995. Curvilinear ridges in Isidis Planitia, Mars – the results of mud volcanism? In: Proc. Lunar Planet. Sci. Conf., vol. 26, pp. 321–322. Dickson, J.L., Head, J.W., Marchant, D.R., 2008. Late Amazonian glaciation at the dichotomy boundary on Mars: evidence for glacial thickness maxima and multiple glacial phases. Geology 36, 411–414. Dickson, J.L., Head, J.W., Marchant, D.R., 2010. Kilometer-thick ice accumulation and glaciation in the northern mid-latitudes of Mars: evidence for crater-filling events in the late Amazonian at the Phlegra Montes. Earth Planet. Sci. Lett. 294, 332–342. Dunlop, P., Clark, C.D., 2006. The morphological characteristics of ribbed moraine. Quat. Sci. Rev. 25, 1668–1691. Durham, W.B., Kirby, S.H., Stern, L.A., 1997. Creep of water ices at planetary conditions: a compilation. J. Geophys. Res., Planets 102, 16293–16302. Durham, W.B., Stern, L.A., Kirby, S.H., 2001. Rheology of ice I at low stress and elevated confining pressure. J. Geophys. Res. 106 (6), 11031–11042. Erkeling, G., Reiss, D., Hiesinger, H., Poulet, F., Carter, J., Ivanov, M.A., Heuber, E., Jaumann, R., 2012. Valleys, paleolakes and possible shorelines at the Libya Montes/Isidis boundary: implications for the hydrologic evolution of Mars. Icarus 219, 393–413. Erkeling, G., Reiss, D., Hiesinger, H., Ivanov, M.A., Hauber, E., Bernhardt, H., 2014. Landscape formation at the Deuteronilus contact in southern Isidis Planitia, Mars: implications for an Isidis Sea? Icarus 242, 329–351.

9.828 0.0057 3.61 ×10−13 1.73 ×103 60 139 105 8.7 × 10−4 8.314 3 273.15

∗ ≤ 263.15 K ∗ > 263.15 K ∗ ≤ 263.15 K ∗ > 263.15 K

Fagents, S.A., Lanagan, P., Greeley, R., 2002. Rootless cones on Mars: a consequence of lava–ground ice interaction. Geol. Soc. (Lond.) Spec. Publ. 202, 295–317. Fassett, C.I., Head, J.W., 2006. Valleys on Hecates Tholus, Mars: origin by basal melting of summit snowpack. Planet. Space Sci. 54, 370–378. Fassett, C.I., Dickson, J.L., Head, J.W., Levy, J.S., Marchant, D.R., 2010. Supraglacial and proglacial valleys on Amazonian Mars. Icarus 208, 86–100. Fastook, J.L., Head, J.W., Marchant, D.R., Forget, F., 2008. Tropical mountain glaciers on Mars: altitude-dependence of ice accumulation, accumulation conditions, formation times, glacier dynamics, and implications for planetary spin-axis/orbital history. Icarus 198, 305–317. Fastook, J.L., Head, J.W., Forget, F., Madeleine, J.-B., Marchant, D.R., 2011. Evidence for Amazonian northern mid-latitude regional glacial landsystems on Mars: glacial flow models using GCM-driven climate results and comparisons to geological observations. Icarus 216, 23–39. Fastook, J.L., Head, J.W., Marchant, D.R., Forget, F., Madeleine, J.-B., 2012. Early Mars climate near the Noachian–Hesperian boundary: independent evidence for cold conditions from basal melting of the south polar ice sheet (Dorsa Argentea Formation) and implications for valley network formation. Icarus 219, 25–40. Fastook, J.L., Head, J.W., 2014. Amazonian mid- to high-latitude glaciation on Mars: supply-limited ice sources, ice accumulation patterns, and concentric crater fill glacial flow and ice sequestration. Planet. Space Sci. 91, 60–76. Fastook, J.L., Head, J.W., Marchant, D.R., 2014. Formation of lobate debris aprons on Mars: assessment of regional ice sheet collapse and debris-cover armoring. Icarus 228, 54–63. Feldman, W.C., Prettyman, T.H., Maurice, S., Plaut, J.J., Bish, D.L., Vaniman, D.T., Mellon, M.T., Metzger, A.E., Squyres, S.W., Karunatillake, S., Boynton, W.V., Elphic, R.C., Funsten, H.O., Lawrence, D.J., Tokar, R.L., 2004. Global distribution of nearsurface hydrogen on Mars. J. Geophys. Res., Planets 109, E9. Fishbaugh, K.E., Head, J.W., 2001. Comparison of the North and South polar caps of Mars: new observations from MOLA data and discussion of some outstanding questions. Icarus 154, 145–161. Fisher, D.A., Hecht, M.H., Kounaves, S.P., Catling, D.C., 2010. A perchlorate brine lubricated deformable bed facilitating flow of the north polar cap of Mars: possible mechanism for water table recharging. J. Geophys. Res., Planets 115, E2. Forget, F., Haberle, R.M., Montmessin, F., Levrard, B., Head, J.W., 2006. Formation of glaciers on Mars by atmospheric precipitation at high obliquity. Science 311, 368–371. Forget, F., Wordsworth, R., Millour, E., Madeleine, J.B., Kerber, L., Leconte, J., Marcq, E., Haberle, R.M., 2013. 3D modelling of the early martian climate under a denser CO2 atmosphere: temperatures and CO2 ice clouds. Icarus 202, 81–99. Gaidos, E., Marion, G., 2003. Geological and geochemical legacy of a cold early Mars. J. Geophys. Res. 108, 5055. Garvin, J.B., Head, J.W., Marchant, D.R., Kreslavski, M.A., 2006. High-latitude coldbased glacial deposits on Mars: multiple superposed drop moraines in a crater interior at 70◦ N latitude. Meteorit. Planet. Sci. 41, 1659–1674. Ghatan, G.J., Head, J.W., 2004. Regional drainage of meltwater beneath a Hesperianaged south circumpolar ice sheet on Mars. J. Geophys. Res. 109, E7006. Ghent, R.R., Anderson, S.W., Pithawala, T.M., 2012. The formation of small cones in Isidis Planitia, Mars through mobilization of pyroclastic surge deposits. Icarus 217, 169–183. Goldsby, D.L., Kohlstedt, D.L., 2001. Superplastic deformation of ice: experimental observations. J. Geophys. Res. 106 (B6), 11017–11030. Gough, D.O., 1981. Solar interior structure and luminosity variations. Sol. Phys. 74, 21–34. Gourronc, M., Bourgeois, O., Mège, D., Pochat, S., Bultel, B., Massé, M., Le Deit, L., Le Mouélic, S., Mercier, D., 2014. One million cubic kilometers of fossil ice in Valles

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

Marineris: relicts of a 3.5 Gy old glacial landsystem along the Martian equator. Geomorphology 204, 235–255. Greve, R., 2004. Evolution of the north-polar cap of Mars: a modelling study. Planet. Space Sci. 52 (9), 775–787. Greve, R., 2008. The polar ice caps of Mars. Low Temp. Sci. 66, 139–148. Grizzaffi, P., Schultz, P.H., 1989. Isidis basin: site of ancient volatile-rich debris layer. Icarus 77, 358–381. Grott, M., Breuer, D., 2010. On the spatial variability of the Martian elastic lithosphere thickness: evidence for mantle plumes. J. Geophys. Res. 115 (E3), E03005. Grott, M., Morschhauser, A., Breuer, D., Hauber, E., 2011. Volcanic outgassing of CO2 and H2 O on Mars. Earth Planet. Sci. Lett. 308 (3–4), 391–400. Guidat, T., Pochat, S., Bourgeois, O., Souˇcek, O., 2015. Landform assemblage in Isidis Planitia, Mars: evidence for a 3 Ga old polythermal ice sheet. Earth Planet. Sci. Lett. 411, 253–267. Haberle, R.M., 1998. Early Mars climate models. J. Geophys. Res. 103, 28467. Haberle, R.M., Murphy, J.R., Schaeffer, J., 2003. Orbital change experiments with a Mars general circulation model. Icarus 161, 66–89. Hättestrand, C., Kleman, J., 1999. Ribbed moraine formation. Quat. Sci. Rev. 18, 43–61. Head, J.W., Marchant, D.R., 2003. Cold-based mountain glaciers on Mars: western Arsia Mons. Geology 31, 641–644. 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. J. Geophys. Res. 106, 12275–12299. Head, J.W., Neukum, G., Jaumann, R., Hiesinger, H., Hauber, E., Carr, M., Masson, P., Foing, B., Hoffmann, H., Kreslavsky, M., Werner, S., Milkovich, S., van Gasselt, S., Team, H.C.-I., 2005. Tropical to mid-latitude snow and ice accumulation, flow and glaciation on Mars. Nature 434, 346–351. Head, J.W., Marchant, D.R., Agnew, M.C., Fassett, C.I., Kreslavsky, M.A., 2006a. Extensive valley glacier deposits in the northern mid-latitudes of Mars: evidence for late amazonian obliquity-driven climate change. Earth Planet. Sci. Lett. 241, 663–671. Head, J.W., Nahm, A.L., Marchant, D.R., Neukum, G., 2006b. Modification of the dichotomy boundary on Mars by Amazonian mid-latitude regional glaciation. Geophys. Res. Lett. 33, L08S03. Head, J.W., Marchant, D.R., Dickson, J.L., Kress, A.M., Baker, D.M., 2010. Northern mid-latitude glaciation in the late Amazonian period of Mars: criteria for the recognition of debris-covered glacier and valley glacier landsystem deposits. Earth Planet. Sci. Lett. 294, 306–320. Hiesinger, H., Head, J.W., 2004. The Syrtis Major volcanic province, Mars: synthesis from Mars Global Surveyor data. J. Geophys. Res. 109, E01004. Hindmarsh, R.C.A., 2004. Thermoviscous stability of ice-sheet flows. J. Fluid Mech. 502, 17–40. Hindmarsh, R.C.A., 2006. Stress gradient damping of thermoviscous ice flow instabilities. J. Geophys. Res. 111, B12409. Hindmarsh, R.C.A., 2009. Consistent generation of ice-streams via thermo-viscous instabilities modulated by membrane stresses. Geophys. Res. Lett. 36, L06502. Hindmarsh, R.C.A., 2011. Ill-posedness of the shallow-ice approximation when modelling thermoviscous instabilities. J. Glaciol. 57 (206), 1177-8 (Correspondence). Holt, J.W., Safaeinili, A., Plaut, J.J., Head, J.W., Phillips, R.J., Seu, R., Kempf, S.D., Choudhary, P., Young, D.A., Putzig, N.E., Biccari, D., Gim, Y., 2008. Radar sounding evidence for buried glaciers in the southern mid-latitudes of Mars. Science 322, 1235–1238. Howard, A.D., 1981. Etched plains and braided ridges of the south polar region of Mars: features produced by basal melting of ground ice? In: Reports of Planetary Geology Program, pp. 286–288. Hubbard, B., Souness, C., Brough, S., 2014. Glacier-like forms on Mars. Cryosphere 8, 2047–2061. Hutter, K., 1983. Theoretical Glaciology; Material Science of Ice and the Mechanics of Glaciers and Ice Sheets. Reidel, Dordrecht, Netherlands. 510 pp. Ivanov, M.A., Head, J.W., 2003. Syrtis Major and Isidis basin contact: morphological and topographic characteristics of Syrtis Major lava flows and material of the Vastitas Borealis formation. J. Geophys. Res. 108, 5063. Ivanov, M.A., Hiesinger, H., Erkeling, G., Hielscher, F.J., Reiss, D., 2012. Major episodes of geologic history of Isidis Planitia on Mars. Icarus 218, 24–46. Ivanov, M.A., Hiesinger, H., Erkeling, G., Reiss, D., 2014. Mud volcanism and morphology of impact craters in Utopia Planitia on Mars: evidence for the ancient ocean. Icarus 228, 121–140. Jakosky, B.M., Carr, M.H., 1985. Possible precipitation of ice at low latitudes of Mars during periods of high obliquity. Nature 315, 559–561. Kadish, S.J., Head, J.W., Barlow, N.G., Marchant, D.R., 2008. Martian pedestal craters: marginal sublimation pits implicate a climate-related formation mechanism. Geophys. Res. Lett. 35, L16104. Kargel, J.S., Strom, R.G., 1992. Ancient glaciation on Mars. Geology 20, 3–7. Kargel, J.S., Baker, V.R., Beget, J.E., Lockwood, J.F., Pewe, T.L., Shaw, J.S., Strom, R.G., 1995. Evidence of ancient continental glaciation in the Martian northern plains. J. Geophys. Res. 100, 5351–5368. Kasting, J.F., 1991. CO2 condensation and the climate of early Mars. Icarus 94, 1–13. Kite, E.S., Hindmarsh, R.C.A., 2007. Did ice streams shape the largest channels on Mars? Geophys. Res. Lett. 34, 0094–8276.

189

Kleman, J., Glasser, N.F., 2007. The subglacial thermal organisation (STO) of ice sheets. Quat. Sci. Rev. 26 (5), 585–597. Kubo, T., Durham, W.B., Stern, L.A., Kirby, S.H., 2006. Grain-size-sensitive creep in ice II. Science 311, 1267–1269. 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, 343–364. Levrard, B., Forget, F., Montmessin, F., Laskar, J., 2004. Recent ice-rich deposits formed at high latitudes on Mars by sublimation of unstable equatorial ice during low obliquity. Nature 431, 1072–1075. Levy, J.S., Head, J.W., Marchant, D.R., 2007. Lineated valley fill and lobate debris apron stratigraphy in Nilosyrtis Mensae, Mars: evidence for phases of glacial modification of the dichotomy boundary. J. Geophys. Res. 112, E08004. Levy, J., Head, J.W., Marchant, D.R., 2010. Concentric crater fill in the northern midlatitudes of Mars: formation processes and relationships to similar landforms of glacial origin. Icarus 209, 390–404. Levy, J.S., Fassett, C.I., Head, J.W., Schwartz, C., Watters, J.L., 2014. Sequestered glacial ice contribution to the global Martian water budget: geometric constraints on the volume of remnant, midlatitude debris-covered glaciers. J. Geophys. Res., Planets 119 (10), 2188–2196. Lockwood, J.F., Kargel, J.S., Strom, R.B., 1992. Thumbprint Terrain on the Northern plains: a glacial hypothesis. Lunar Planet. Sci. XXIII, 795–796. Abstract. MacAyeal, D.R., 1992. Irregular oscillations of the West Antarctic Ice Sheet. Nature 359, 29–32. Madeleine, J.-B., Forget, F., Head, James W., Levrard, B., Montmessin, F., Millour, E., 2009. Amazonian northern mid-latitude glaciation on Mars: a proposed climate scenario. Icarus 203 (2), 390–405. Mangold, N., Maurice, S., Feldman, W.C., Costard, F., Forget, F., 2004. Spatial relationships between patterned ground and ground ice detected by the Neutron Spectrometer on Mars. J. Geophys. Res. 109, E08001. Marshall, S.J., 2005. Recent advances in understanding ice sheet dynamics. Earth Planet. Sci. Lett. 240, 191–204. Mège, D., Bourgeois, O., 2011. Equatorial glaciations on Mars revealed by gravitational collapse of Valles Marineris wallslopes. Earth Planet. Sci. Lett. 310, 182–191. Mellon, M.T., Arvidson, R.E., Sizemore, H.G., Searls, M.L., Blaney, D.L., Cull, S., Hecht, M.H., Heet, T.L., Keller, H.U., Lemmon, M.T., Markiewicz, W.J., Ming, D.W., Morris, R.V., Pike, W.T., Zent, A.P., 2009. Ground ice at the phoenix landing site: stability state and origin. J. Geophys. Res. 114, E1. Milkovich, S.M., Head, J.W., Pratt, S., 2002. Meltback of Hesperian-aged ice-rich deposits near the south pole of Mars: evidence for drainage channels and lakes. J. Geophys. Res. 107 (E6), 10-1–10-13. Milkovich, S.M., Head, J.W., Marchant, D.R., 2006. Debris-covered piedmont glaciers along the northwest flank of the Olympus Mons scarp: evidence for low-latitude ice accumulation during the late Amazonian of Mars. Icarus 181, 388–407. Mischna, M.A., Richardson, M.I., Wilson, R.J., McCleese, D.J., 2003. On the orbital forcing of Martian water and CO2 cycles: a general circulation model study with simplified volatile schemes. J. Geophys. Res. 108, E65062. Morgan, G.A., Head, J.W., Marchant, D.R., 2009. Lineated valley fill (LVF) and lobate debris aprons (LDA) in the Deuteronilus Mensae northern dichotomy boundary region, Mars: constraints on the extent, age and episodicity of Amazonian glacial events. Icarus 202, 22–38. Mouginot, J., Pommerol, A., Kofman, W., Beck, P., Schmitt, B., Herique, A., Grima, C., Safaeinili, A., Plaut, J., 2010. The 3–5 MHz global reflectivity map of Mars by MARSIS/Mars express: implications for the current inventory of subsurface H2 O. Icarus 210, 612–625. Neumann, G.A., Zuber, M.T., Wieczorek, M.A., McGovern, P.J., Lemoine, F.G., Smith, D.E., 2004. Crustal structure of Mars from gravity and topography. J. Geophys. Res. 109, E08002. Newman, C.E., Lewis, S.R., Read, P.L., 2005. The atmospheric circulation and dust activity in different orbital epochs on Mars. Icarus 174, 135–160. Paterson, W.S.B., 1981. The Physics of Glaciers. Pergamon Press, Oxford. Pattyn, F., Perichon, L., Aschwanden, A., Breuer, B., de Smedt, B., Gagliardini, O., Gudmundsson, G.H., Hindmarsh, R., Hubbard, A., Johnson, J.V., Kleiner, T., Konovalov, Y., Martin, C., Payne, A.J., Pollard, D., Price, S., Rückamp, M., Saito, F., Souˇcek, O., Sugiyama, S., Zwinger, T., 2008. Benchmark experiments for higher-order and full Stokes ice sheet models (ISMIP-HOM). Cryosphere 2 (2), 95–108. Payne, A.J., 1995. Limit cycles in the basal thermal regime of ice sheets. J. Geophys. Res. 100 (B3), 4249–4263. Payne, A.J., Dongelmans, P., 1997. Self-organization in the thermomechanical flow of ice-sheets. J. Geophys. Res. 102 (B6), 12219–12233. Plescia, J.B., 1980. Cinder cones of Isidis and Elysium. NASA Tech. Memo. 82385, 263–265. Ritzer, J.A., Hauck, S.A., 2009. Lithospheric structure and tectonics at Isidis Planitia, Mars. Icarus 201, 528–539. Ruiz, J., McGovern, P.J., Jiménez-Díaz, A., López, V., Williams, J.P., Hahn, B.C., Tejero, R., 2011. The thermal evolution of Mars as constrained by paleo-heat flows. Icarus 215, 508–517. Saito, F., Abe-Ouchi, A., Blatter, H., 2006. European Ice Sheet Modelling Initiative (EISMINT) model intercomparison experiments with first order mechanics. J. Geophys. Res. 111, F02012.

190

O. Souˇcek et al. / Earth and Planetary Science Letters 426 (2015) 176–190

Schultz, R.A., Frey, H.V., 1990. A new survey of multiring impact basins on Mars. J. Geophys. Res. 95 (B9), 14175–14189. Shean, D.E., Head, J.W., Marchant, D.R., 2005. Origin and evolution of a cold-based tropical mountain glacier on Mars: the Pavonis Mons fan-shaped deposit. J. Geophys. Res. 110, E05001. Shean, D.E., Head, J.W., Fastook, J.L., Marchant, D.R., 2007. Recent glaciation at high elevations on Arsia Mons, Mars: implications for the formation and evolution of large tropical mountain glaciers. J. Geophys. Res. 112, E03004. Shean, D.E., 2010. Candidate ice-rich material within equatorial craters on Mars. Geophys. Res. Lett. 37, L24202. Skinner, J.A., Mazzini, A., 2009. Martian mud volcanism: terrestrial analogs and implications for formational scenarios. Mar. Pet. Geol. 26, 1866–1878. Smith, D.E., Sjogren, W.L., Tyler, G.L., Balmino, G., Lemoine, F.G., Konopliv, A.S., 1999. The gravity field of Mars: results from Mars Global Surveyor. Science 286 (5437), 94-7. Souˇcek, O., Martinec, Z., 2008. Iterative improvement of the shallow-ice approximation. J. Glaciol. 54 (188), 812–822. Souˇcek, O., 2010. Numerical modeling of ice sheet dynamics. PhD thesis. Faculty of Mathematics and Physics, Charles University in Prague. Souˇcek, O., Martinec, Z., 2011. ISMIP-HEINO experiment revisited; effect of higherorder approximation and sensitivity study. J. Glaciol. 57 (206), 1158–1170.

Souness, C., Hubbard, B., Milliken, R.E., Quincey, D., 2012. An inventory and population-scale analysis of martian glacier-like forms. Icarus 217, 243–255. Squyres, S.W., 1978. Martian fretted terrain — flow of erosional debris. Icarus 34, 600–613. Squyres, S.W., 1979. The distribution of lobate debris aprons and similar flows on Mars. J. Geophys. Res. 84, 8087–8096. Tanaka, K., Skinner, J., Hare, T., 2005. Geologic map of the northern plains of Mars. U.S. Geol. Surv. Sci. Invest. Ser. Map 2888. Wichman, R.W., Schultz, P.H., 1989. Sequence and mechanisms of deformation around the Hellas and Isidis impact basins on Mars. J. Geophys. Res. 94, 17333–17357. Wilhelms, D.E., 1973. Comparison of martian and lunar multiringed circular basins. J. Geophys. Res. 78, 4084–4095. Wordsworth, R., Forget, F., Millour, E., Head, J.W., Madeleine, J.-B., Charnay, B., 2013. Global modelling of the early Martian climate under a denser CO2 atmosphere: water cycle and ice evolution. Icarus 222, 1–19. Zuber, M.T., Solomon, S.C., Phillips, R.J., Smith, D.E., Tyler, G.L., Aharonson, O., Balmino, G., Banerdt, W.B., Head, J.W., Johnson, C.L., Lemoine, F.G., McGovern, P.J., Neumann, G.A., Rowlands, D.D., Zhong, S., 2000. Internal structure and early thermal evolution of Mars from Mars Global Surveyor topography and gravity. Science 287 (5459), 1788–1793.