Parsons et al. Progress in Earth and Planetary Science
https://doi.org/10.1186/s40645-020-0323-9(2020) 7:13

Progress in Earth andPlanetary Science

RESEARCH ARTICLE

Open Access

Cold-based glaciation of Pavonis Mons,
Mars: evidence for moraine deposition
during glacial advance
Reid A. Parsons1,2* , Tomohiro Kanzaki3, Ryodo Hemmi1 and Hideaki Miyamoto3

Abstract
The three large volcanoes in the Tharsis region of Mars: Arsia, Pavonis, and Ascraeus Montes all have fan-shaped
deposits (FSDs) on their northern or western flanks consisting of a combination of parallel ridges, knobby/
hummocky terrain, and a smooth, viscous flow-like unit. The FSDs are hypothesized to have formed in the
Amazonian during a period of high spin-axis obliquity which redistributed polar ice to the equatorial Tharsis region
resulting in thick (> 2 km), flowing ice deposits. Based on previous ice flow simulations and crater surveys, theridges are interpreted to be recessional drop moraines formed as debris on the ice sheet surface was transported
to the ice margin—forming a long ridge sequence over an extended (∼100 Myr) period of ice sheet retreat. We test
this hypothesis using a high-resolution, thermomechanical ice sheet model assuming a lower ice loss rate (~ 0.5mm/year) than prior work based on new experimental results of ice sublimation below a protective debris layer.
Our ice flow simulation results, when combined with topographic observations from a long sequence of ridges
located interior of the Pavonis FSD, show that the ridged units were more likely deposited during one or more
periods of glacial advance (instead of retreat) when repetitive pulses (approx. 120 kyr periodicity) of ice
accumulation during high obliquity produced kinematic waves which advected a large volume of surface debris tothe ice margin. If ridge deposition does occur during glacial advance, it could explain the cyclic pattern of ridge
spacing and would link the dominant, 120 kyr periodicity in obliquity to the time interval between adjacent ridges.
By measuring the spacing between these ridges and applying this timescale, we constrain the velocity of glacial
margin to be between 0.2 and 4 cm/Earth year—in close agreement with the numerical simulation. This reinterpretation of the FSD ridged unit suggests that the timescale of FSD formation (and perhaps the duration of the
Amazonian high obliquity period) was shorter than previously reported.
Keywords: Thermomechanical ice sheet model, Mars, Amazonian climate, Moraine deposition, Tharsis, Glaciation

Introduction
The fan-shaped deposits (FSDs) on the Tharsis volcanoes (Arsia, Pavonis, and Ascraeus) (Fig. 1a) were firstdescribed in detail using images from the Mariner 9 andViking missions (Zimbelman and Edgett 1992). The
outer margin of FSDs consists of a series of parallelridges which have a lobate shape when viewed fromabove (black and white dotted line, Fig. 1b). A combination of smooth and knobby units is generally found in* Correspondence: rparson4@fitchburgstate.edu
1
University Museum, The University of Tokyo, 7 Chome-3-1 Hongo, Bunkyo,
Tokyo 113-0033, Japan
2
Earth and Geographic Sciences, Fitchburg State University, 160 Pearl St.,
Fitchburg, MA 01420, USA
Full list of author information is available at the end of the article

the interior of the FSDs which are most prominent on
Arsia and Pavonis Mons. Early research efforts invokedvarious processes for FSD formation including gravitydriven landslides (Carr et al. 1977; Scott and Tanaka
1986), glaciation (Lucchitta 1981; Williams 1978), andpyroclastic volcanism (Zimbelman and Edgett 1992).
Later support for glaciation of the Tharsis Montes camefrom the prediction from climate models of ice accumulation during times of high Martian spin-axis obliquity(Forget et al. 2006). During periods of high obliquity, thepoles receive more insolation, and water ice sublimatesin numerical simulations by Jakosky and Carr (1985).
This water vapor is then transported equator-ward
where orographic lifting and adiabatic cooling of air© The Author(s). 2020 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0
International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to
the Creative Commons license, and indicate if changes were made.

Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 2 of 21

Fig. 1 Topography of the Arsia and Pavonis FSDs. a A true color image mosaic from Viking data of the western hemisphere of Mars. b Detailed
elevation (colored regions) of the FSDs (outlined with black/white dotted line) indicate rougher topography within the FSDs and presence of
volcanic rifts at Arsia. White line with labeled distances gives the location of the flowline used in numerical model. c HRSC image of the outermargin of the Arsia FSD where the intersection of moraine-like ridges (arrows) indicates that glacial activity was not capable of eroding previous
ridges during subsequent glacial advancemasses moving over the Tharsis volcanoes results in the
precipitation and accumulation of ice in numericalclimate models (Forget et al. 2006). The geologic observations and climate model results were combined in astudy by (Fastook et al. 2008) to constrain the likely obliquity history which allowed these deposits to accumulate,attain their maximum extent, and retreat in accordancewith crater age dates for the deposits (which indicate an
exposure age of ∼300 Myr).
In this study, we investigate the influence of timevarying obliquity on the ice margin flow velocity, location, and shape using the thermomechanical Parallel Ice
Sheet Model (PISM). Our simulation is compared withthe observed spacing and cross-sectional area of ridges
found within the Pavonis FSD in order to determine
whether moraine deposition is more likely to occur during advance or retreat of the ice sheet. We firstsummarize the prior work regarding the FSD ridged unit
and its geologic interpretation. Next, we describe thesize and location of ridges at Pavonis Mons using CTX
images and DEMs generated from CTX stereo pairs at
both Pavonis and Arsia Montes. In the following section,
we provide a description of the numerical model, its

inputs, and assumptions. The final sections compare thesimulation results with the observations, discuss the process(es) responsible for ridge emplacement, and summarize
our findings.
Prior geologic interpretation

The arcuate, distal portion of the FSDs generally consistsof a series of parallel ridges (∼1 km spacing) with variable heights ranging from 5 to 50 m based on MOLA
laser altimetry observations conducted by Head et al.
(2003) at Arsia Mons. The delicate draping of ridges
over underlying features (craters and lava flows) without
modification (Head and Marchant 2003) is a characteristic shared by cold-based drop moraines in Antarctica(Bockheim 2010; Bockheim 2002; Marchant et al. 1994;
Staiger et al. 2006). The regularity of ridge spacing andlong lateral extent contrasts with landslide deposits
found surrounding Olympus Mons (Francis and Waadge1983) and contrasts with ridges deposited from turbulent pyroclastic flows (Head and Marchant 2003).
Geologic mapping efforts in the Tharsis region havemade use of crater abundances and superposition relationships to conclude that the FSD formed in the late(2020) 7:13

Page 3 of 21

Amazonian and were likely contemporaneous with themost recent volcanic activity on both the volcanic edifices and radial rift zones of Arsia, Pavonis, and Ascraeus
Mons. Landforms such as tuyas and tiered, mounded,leveed and digitate lava flows, which are indicative oflava-ice interaction, have been found in all three FSDs
(Kadish et al. 2008; Scanlon et al. 2014; Scott et al. 1998;
Shean et al. 2005).
The most recent crater density measurements for
Arsia, Pavonis, and Ascraeus give best-fit ages of 210,
125, and 220 Myr, respectively (Kadish et al. 2014) using(Hartmann 2005) isochrons. Kadish et al. (2014) separated the outer ridged unit of the Arsia FSD into three
concentric regions and found a progressive increase in
age with distance away from the Arsia edifice. The authors argued that such a relationship would be expected
for a thick, slowly retreating ice sheet which shields theinterior of the FSD from all but the largest impacts.
Ideally, cosmogenic radionuclide measurements could
provide surface exposure ages to document glacial recession (e.g., Staiger et al. 2006), but such data has not yet
been acquired on Mars.
The ridged units in the Martian FSDs are assumed to
be deposited during glacial recession based on theirsimilarity to terrestrial analogs in the Dry Valleys of Antarctica (Head and Marchant 2003) and Iceland (Lucchita
1984; Williams 1978). Although the crater density measurements made by Kadish et al. (2014) appear to support this hypothesis, they do not rule out a scenario in
which ridges are deposited during glacial advance. If moraines are deposited at the ice margin during advance,then subsequent burial of the ridge by the overriding ice
sheet would shield the ridge from further impacts until
obliquity changes gradually allowed the ice to sublimateback into the atmosphere. Thus, the crater density observations provide a surface exposure age estimate (timesince ice sheet removal) and do not represent the age offeatures formed before or during ice sheet advance. Although the crater density measurements made by Kadish
et al. (2014) help constrain the rate of ice margin retreat,
they do not provide a test for determining whether ridge
deposition happens during ice retreat or advance. Thus,
additional observations and theoretical work needs to bedone in order to test these hypotheses.

extent and rates of ablation and accumulation but is stationary when the total mass accumulated over the ice
sheet is equal to the mass ablated (steady state). The
conundrum arises from the observation that the moraines are so closely spaced (∼ 1 km) which accounts fora change of less than 1% of the length of the ice deposit
assuming precipitation occurs on the slopes of the volcanic edifices as indicated in GCM results by Forgetet al. (2006). Such small changes in the ice margin position would involve similarly small changes in the average ablation and/or accumulation rates over the ∼ 1 Myr(Kadish et al. 2014) timescale required to deposit an individual ridge. How did the glacial mass balance continuously change by such small amounts in order to form
hundreds of ridges at Arsia (Kadish et al. 2014) when thedominant frequency of obliquity change happens over
120 kyr timescales (Schorghofer 2008; Ward 1992)?
Superposition of remnant ice deposits (smooth unit)
on ridges (Fig. 1d, (Shean et al. 2005) and intersecting
ridge sequences (Fig. 1c) provide geomorphic evidence
for episodes of glacial advance following times of retreat
within the Pavonis and Arsia FSDs. These features suggest several episodes of retreat and advance may havetaken place (Shean et al. 2005) and also show that ridges
can remain preserved even after being overrun by a later
episode of ice sheet advance which may complicate interpretation of the ridge formation sequence.

Parsons et al. Progress in Earth and Planetary Science

Problems with current interpretation

The interpretation of the ridged unit as a series of coldbased glacier moraines deposited over a long duration(∼100 Myr) of retreat of the ice margin poses a conundrum regarding the Amazonian climate. In order to deposit a moraine, the ice margin needs to be roughly
stationary while ice flow continues to advect ice and englacial/supraglacial debris to the ice margin. The location of the ice margin is a sensitive function of the aerial

Potential ridge formation mechanisms

The dry, cold Amazonian climate of Mars limits the scopeof potential ridge-forming mechanisms to those which
would operate under cold-based, dry glacial conditions.
Given the dominance of obliquity on the Amazonian climate (e.g., Tanaka 2004; Schorghofer 2008; Schorghofer2006), we assume obliquity-driven climate change is primarily responsible for ridge formation and that the formation of a surface sublimation lag is the primary debris
source. This lag acts to reduce the rate of ice ablation via
sublimation. Given these assumptions, there are at least
four potential processes for forming ridges:
1) During retreat of the ice margin, relatively warm
periods (during times of low obliquity) may
temporarily enhance ice flow and stall the
migration of the ice margin to result in a ridge.
2) During glacial advance, the margin of the ice
deposit may be slowed during times of highobliquity (cooler ice temperatures) which could
allow even a low sublimation rate to maintain atemporarily fixed ice margin position—allowing
debris to accumulate and a ridge to form.
3) Renewed ice accumulation at the head of a
sublimating ice deposit during high obliquity may
produce a kinematic (velocity) wave which mayParsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 4 of 21collect surface debris and deposit a ridge uponreaching the ice margin.
4) The steepness of the ice margin is expected toincrease along with the ice viscosity during cold
temperatures experienced during high obliquity atthe equator which may cause slope failure of thesurface lag and the formation of a ridge at the ice
margin.
These processes may work in concert or independently
to produce ridges of varying sizes depending on thedominant process at a given point in time. Furthermore,this list is not exhaustive but represents the most likely
mechanisms based on our assumption of cold-based
glacial conditions.
We investigate the process(es) responsible for depositing
ridges by combining results from a thermomechanical ice
sheet simulation with an analysis of the morphology of aseries of ridges north of Pavonis Mons. Our model inputs
assume ice accumulates in accord with GCM results athigh obliquity and ice ablation occurs at a rate associatedwith the formation of a protective lag cover.
Understanding the process of ridge formation is an
essential component to decoding the Amazonian climatechange record contained in the FSDs and for understanding how active the Martian water cycle was in the
recent past.

Methods
Stereo images and observations of moraine
characteristics

In contrast with ridges in the Arsia FSD, the Pavonis
Mons ridges are generally wider, taller, are exposed overa larger portion of the FSD, and exhibit a quasi-periodic
pattern in their spacing which make them an intriguingfeature to study (Figs. 2 and 3). A CTX stereo DEM (18m/pixel resolution) was produced at the distal portion of
the Arsia FSD (Fig. 2a) and the medial section of the
Pavonis FSD (Fig. 2b) using the Ames Stereo Pipeline.
DEM production made use of raw MOLA profiles ascontrol points to improve quality (Beyer et al. 2018).
The vertical precision of the DEMs are ≈ 2.46 m (seeHemmi and Miyamoto 2017 for a DEM generation
workflow description).
The spacing between ∼ 40 ridges are measured using
the CTX global mosaic (5.5 m/px) of the FSD interior(Fig. 3a). Depending on ridge distinctness, three to fivemeasurements are obtained between each pair of ridges
in a direction perpendicular to the ridge crests. A combination of CTX images and a stereo CTX DEM areused to identify ridges with varying degrees of confidence (Fig. 3b). In some cases, ridges are indistinct in
the imagery but resolved in the DEM. Highconfidence ridge identifications are prominent in bothFig. 2 Topography and morphology of the distal portion of theArsia (a) and Pavonis (b) FSDs. CTX imagery and transparent stereo
topography (colored regions) shows ridge sequences of parallelridges together with smooth (labeled with “S”) and knobby (labeledwith “K”) units. The smooth unit (interpreted as remnant debriscovered ice by Shean et al. 2005) superposes ridges in the Pavonis
FSD and indicates that the flow of ice responsible for forming thesmooth unit was not erosive. Ridge spacing and cross-sectional area
measurements using this data are shown in Figs. 3 and 4

the imagery and the DEM (black numbers/lines inFig. 3a, b), while lower-confidence identifications(blue numbers/lines in Fig. 3) are only expressed in
one of the two datasets. Finally, the dotted lines inFig. 3a indicate that the ridge is more prominent at alocation outside the region shown.
The mean spacing between ridges varies in a quasiperiodic fashion in which ridges are densely spaced in
the southernmost portion of the FSD, become more
widely spaced to the north, with this pattern repeating
again (Fig. 3b) before the ridges become extremely faintin the northern portion of the FSD. A few ridges near
the FSD perimeter (along profile B-B′ in Fig. 2b) arevery prominent. Error bars in Fig. 3b are associated withthe standard deviation in the ridge spacing (3 to 5 measurements per ridge pair) and are generally smaller than
the width of the plotted. The line graph in Fig. 3b isshown at the same horizontal scale as Fig. 3a such that

Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 5 of 21

Fig. 3 CTX imagery shown with transparent stereo topography of interior ridge sequence and calculated ridge spacing. a Section of CTX global
mosaic (rotated to correspond to the x-axis of b) illustrating the emergence of parallel ridges from beneath the smooth unit (upper left) withblack lines/numbers indicating ridges which are distinct in both CTX imagery and the stereo DEM (blue ridges are indistinct in one of thesedatasets). Dotted lines indicate poor ridge expression at this location, but, if black, has a more distinct section further to the east. b Mean spacing
between individual ridge pairs (using three to five measurements per ridge pair) as a function of distance along model flowline (Fig. 1b) witherror bars reflecting the standard deviation in the spacing between ridge pairs (generally smaller than the line width). The black line indicatesspacing between high-confidence ridges only whereas blue line includes all ridges

the positions of points along the x-axis of Fig. 3b lie approximately halfway between ridges shown in Fig. 3a.
The amount of material within the ridges provides anadditional constraint on the deposition timescale and debris layer thickness (Kadish et al. 2014) as will be exploredin the “Results and discussion” section. We extracted several topographic profiles (A-A′, B-B′, C-C′, and D-D′ inFig. 2b) in order to determine the ridge cross-sectional
area (shaded regions in Fig. 4a–c). The results are

summarized in Table 1 for measurements of the 26 mostprominent ridges which can be resolved by the CTX stereo DEM. The ridges found at the perimeter of individual
flow units (identified using sets of parallel ridges) aremuch larger than the inner ridges (compare black line in
Fig. 4b with magenta and blue lines in 4c) which arethemselves larger than ridges found at Arsia Mons (redlines in Fig. 4b, c; Head and Marchant 2003). Individual
ridge width and cross-sectional area measurements from

Fig. 4 Ridge topographic profiles and cross-sectional areas measured from the CTX stereo DEM shown in Fig. 2. a This section of the FSD
consists of outer ridge, smooth, and inner ridge units which rest on a regional slope of 0.1°. We, in agreement with prior work, interpret the outer
ridges as terminal moraines (b) which are much larger than the interior ridges (c). The red lines in (b) and (c) indicate that ridges at Arsia aresmaller and more closely spaced than at Pavonis (black, magenta, and blue lines). d The increase in cross-sectional area (y-axis) with ridge width(x-axis) generally follows the trend given by an isosceles triangle with flank slopes of 3° (dashed line) although the smaller, inner ridges tend tohave steeper slopes, but well below the angle of repose

(2020) 7:13

Parsons et al. Progress in Earth and Planetary Science

Page 6 of 21

Table 1 Pavonis ridge statistics
NumberRidge width

X-sectional area

All ridges

26

660 ± 480 m

9000 ± 12,500 m2

Non-terminal ridges

23

500 ± 210 m

5200 ± 3800 m2

Ridge spacing

Extracted from CTX stereo DEM aExtracted from CTX image mosaic
Prominent, non-terminal ridges

30

1.7 ± 1.3 km

Prominent + minor, non-terminal ridges

42

1.3 ± 1.1 kmaGenerated from CTX images D08_030309_1850_XI_05N115W and D08_030375_1850_XI_05N115W

each topographic profile are shown in Fig. 4d. The increase in cross-sectional area with ridge width generally
follows that of an isosceles triangle with flank slopes of ∼3° (dashed line), although many of the smaller ridges havesteeper flanks (Fig. 4d). This slope morphology contrasts
with terrestrial washboard moraines, ribbed moraines, and“push ridges” which have heights of less than 5 m and exhibit steep (> 10°), asymmetric slopes in which the distal
slope is steeper than the proximal (Price 1970).
Both the ridge spacing (which constrains the ice margin velocity if a 120 kyr obliquity-driven formation timescale of ridge pairs is assumed) and cross-sectional area
provide useful information which can be combined with
the ice sheet simulation results.
Numerical simulation

To investigate the influence of obliquity on glacier dynamics and to constrain the ice flow velocity, we performed numerical simulations by adapting the opensource Parallel Ice Sheet Model (PISM) (Winkelmann
et al. 2011) to Martian conditions. PISM is a thermomechanical ice sheet and ice shelf model which accounts
for the influence of a basal heat flux, strain heating, andclimate forcings on the ice temperature. Although themodel is capable of 3D simulations which incorporate
basal sliding, we only require the shallow ice approximation component of the model which simulates deformation and flow of a grounded, cold-based, ice sheet (see,for example, Greve 1997; Marshall et al. 2000; Payne
1999). PISM solves the following pair of time-evolving
partial differential equations describing mass and energy
conservation:*∂H¼ M−∇ Q
∂t

ð1Þ

*
∂T
k ∂2 T− U ∇T þ Σ¼
2∂tρcp ∂zð2Þwhere H is the ice thickness, M is the surface ice mass*balance, and Q is the ice flux (vertical integral of thehorizontal velocity). The change in temperature (T) withtime (Eq. 2) is given by the vertical conduction,

horizontal advection, and strain heating (Σ) terms (given
from left to right) (see Bueler et al. 2007; Winkelmann
et al. 2011 for details) where k is thermal conductivity,
ρ = 910 kg m− 3 is the ice density, cp = 2009 J kg− 1 K− 1 is*

the ice specific heat capacity, and U is the horizontal
velocity.
To calculate the horizontal velocity at elevation (z) thesimulation vertically integrates the shear strain rate
given by Glen’s Law (n = 3, Glen (1955)) which gives:*n

U ¼ −2ðρg Þ ∇h

n−1

Z z





AðT Þðh−ζ Þdζ ∇hð3Þbwhere h is the ice surface elevation and A(T∗) is atemperature-dependent ice softness parameter given by(Paterson and Budd 1982). Equations 1, 2, and 3 aresolved on a 3-D rectangular grid with a horizontal spacing of 500 m and a vertical spacing that decreases from170 to 27 m toward the base of the ice following a quadratic function (PISM authors 2014). The elevation of thebase of the ice and the ice surface are given by b and h,respectively in Eq. 3.
Laboratory experiments by Goldsby and Kohlstedt(2001) and Durham et al. (2010) suggest that the valueof n decreases at low differential stress (< 100 kPa) and isalso dependent on ice grain size despite terrestrial field
evidence for n = 3 behavior even at low differential stress(Budd and Jacka 1989; Cuffey and Paterson 2010). The
flow law proposed by Goldsby and Kohlstedt (2001) generally results in faster ice flow rates than Glen’s Law at
low stress which motivated (Pettit and Waddington
2003) to adopt an additional, n = 1, term to Glen’s Law
in order to model ice deformation at low differential
stress such as encountered near an ice divide. Although
the lower gravity on Mars will tend to reduce the driving
stress, we will be focusing on thick (∼ 2 km) ice sheets
which generally result in differential stresses in excess of
100 kPa (see “Results and discussion” section) where n =3 behavior tends to dominate (Parsons et al. 2011).
However, use of Glen’s Law may underestimate ice flow
velocities on Mars especially for thin or gently sloping
ice deposits.

Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

Application of the PISM model to Mars

We are the first to apply PISM to Mars, and therefore adjusted a number of configurable parameters.
We change the value for gravity to 3.728 m s− 2, disable basal sliding (cold-based glaciation), and apply aspatially and temporally constant basal heat flux of
30 mW m− 2 based on the estimated flux for Tharsis(Plesa et al. 2016), which is somewhat higher than theestimated global average of 20 mW m− 2. Thermal
conductivity of the ice is allowed to change withtemperature due to the colder Martian temperaturesusing the following relationship from Greve andBlatter (2009):
k ðT Þ ¼ 9:98e−0:0057T

ð4Þ

where T is the temperature in Kelvin. This equation results in k = 2.85 W m− 1 K− 1 at T = 220 K compared witha value of 2.1 W m− 1 K− 1 typical of terrestrial studies(PISM authors 2014) resulting in a factor of 1.36 fasterrate of heat conduction within Martian ice deposits
compared with terrestrial ice sheets.
The ice flow simulation was conducted on the 128
pixels per degree Mars Orbiter Laser Altimeter (MOLA)
DEM regridded to 500 m resolution. In order to reduce
computational demand associated with a high spatial resolution, we use PISM in a flowline mode where flow occurs
along a single axis (x-direction) which follows the steepest
gradient down the northwest flank of Pavonis Mons andthrough the FSD (white line in Fig. 1b). For the purposesof determining the flowline trajectory, we apply an averaging window (50 km in radius) to the MOLA topography.
The steepest downslope direction on the smoothed terrain
starting from an initial point near the summit of Pavonis
defines the flowline location. Once the flowline is calculated, the 500 m/pixel resolution MOLA elevation values
along the flow line are adopted as the basal topography.
The smoothing avoids a meandering course or one thatstops in small depressions to produce a downslope trajectory that would more closely resemble that of a thick ice
sheet. Despite this smoothing, a broad depression within
the FSD north of Pavonis Mons, interior to the ridged
unit, causes the downslope profile to terminate. We forced
the profile to continue on a northward trajectory in orderto follow the central axis of the Pavonis FSD (Fig. 1b).
Ice accumulation

The model assumes an ice-free initial condition at an
initial start time of t = − 26 Myr. As time progresses, ice
accumulation occurs over a zone on the flank of Pavoniswith basal topography between 5 and 12 km during
times when obliquity approaches or exceeds 45° basedon high-obliquity GCM results from Forget et al. (2006).
Figure 5 gives the relationship between ice accumulation

Page 7 of 21(y-axis) and obliquity assumed by our model. The shapeand amplitude of this logistic function was chosen after
several iterations of simulating ice flow on Pavonis Mons
in order to produce an ice deposit which remained
within the extent of the FSD but extended beyond theset of interior ridges shown in Fig. 3a.
Snapshots of the model during high (47o at t = − 25.1
Myr) and low (30o at t = − 24.8 Myr) obliquity are shownin Fig. 6 (circle and square symbols, respectively, ininset) to illustrate our mass balance assumptions. The
extent of the flowline model domain is given by thewhite line in Fig. 6a, with dark and light blue lines indicating the extent of ice at t = − 25.1 and − 24.8 Myr, respectively. Topographic profiles and instantaneous massbalance conditions at these times are given in panels b
and c, while panels d and e give ice thickness profiles.
The shaded region in Fig. 6d indicates the extent of theice accumulation zone during high obliquity.
The GCM result from Forget et al. (2006) produces accumulation at Pavonis Mons at an approximate averageaccumulation rate of 35 mm/year (Earth year) at an obliquity of 45° (compared with the blue line in Fig. 6b) assuming a water vapor source at the northern pole,although a south polar water source also results in accumulation at Pavonis. According to Forget et al. (2006), ice
accumulation in the GCM was insensitive to obliquity atvalues larger than 40°, but precipitation may decrease once
polar source regions become depleted over long durationsof obliquity > 45° when ice is unstable at the poles(Mischna et al. 2003). We assume that the accumulation
zone of Pavonis Mons transitions from net ablation to net
accumulation over an obliquity range of 2° centered at 45°
following a logistic function (Fig. 5) under the assumptionthat precipitation is limited to times when polar ice becomes unstable. Moving the obliquity threshold for precipitation to lower values (by even a couple of degrees)
resulted in many more ice accumulation events whichcauses ice to flow beyond the extent of the FSD.
Ice ablation

The ablation rate is calculated assuming that the majority of ice loss occurs while ice is buried beneath a protective lag of 0.5 m thickness. Sublimation rates of
exposed ice at mid-latitudes under current Mars climate
conditions are estimated to be 4.1 mm/year (Earth years)
based on repeat imagery (at 43.6° N) of ice exposed by arecent impact (Dundas and Byrne 2010). Prior ice flow
simulations for the Amazonian Tharsis ice sheets by Fastook et al. (2008) calculated sublimation rates associatedwith buoyancy- and turbulence-driven processes (Ingersoll 1970; Pathare and Paige 2005) which gives ≈ 5 mm/year net ablation rates at elevations just below zones of
net accumulation. Although our calculated sublimation
rate may be an underestimation of ice loss when glacial

Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 8 of 21

Fig. 5 Obliquity-dependent precipitation rates assumed for the accumulation zone. The model assumes precipitation follows a logistic curve withthe transition from ablation to precipitation (35 mm/year) occurring over a 2° window in obliquity (shaded region) centered at 45°. The rate ofablation varies with temperature (dependent on obliquity and surface slope, see Figs. 8 and 9)

ice in the Tharsis region is newly exposed to lowobliquity climate conditions, the ice will become protected by a surface lag as dust entrained in the ice accumulates on the surface as a result of ice loss.
Observational evidence for a reduced ice sublimation
rate comes from the remnant ice protected by surface
debris which makes up the smooth unit (Shean et al.
2005). Assuming that the maximum thickness of ice wasless than ∼ 3 km (based on a perfectly plastic ice rheology (Shean et al. 2005)) and that the age of that ice
must exceed the time since the most recent episode ofhigh obliquity (> 5 Myr; Additional file 1: Figure S1b),then the average sublimation rate must have been lessthan 0.6 mm/year in order for the smooth unit ice topersist until today. This rough estimation is of similar
magnitude to the ∼ 0.1 mm/year long-term (200 kyr)average ablation rate for massive ice buried below a 0.3
m active layer in the permafrost of Beacon Valley,
Antarctica (Liu 2014).
In order to evaluate the timescale of the lag formation
process and to determine the ice ablation rate in themodel, we calculated the sublimation rate (Es) of buried,
non-deforming ice using a water vapor diffusion equation modified from Chevrier et al. (2007):Fig. 6 Model snapshots illustrating assumptions regarding the timing and location of ice accumulation. Inset plot shows obliquity history at the
beginning of the simulations with time and obliquity on x and y axes, respectively. Precipitation occurs at times indicated by the filled blue
regions (see Fig. 6). a MOLA hillshade contour (1 km interval) map indicating location of model flowline (white line) and extent of ice in themodel (transparent blue lines) at − 25.1 (circle in inset) and − 24.8 Myr (square in inset). Topographic profiles of the ice deposit at − 25.1 and −24.8 Myr are shown in b and c, respectively with the blue line (and right y-axis) indicating ice accumulation/ablation rate. Ice thickness for thesetimes are shown in d and e with the shaded region in d indicating extent of ice accumulation

Parsons et al. Progress in Earth and Planetary ScienceEs ¼

DM H2O psat
ð1−RHÞ
LRT s ρ(2020) 7:13

ð5Þ

where D ≈ 5 × 104 m2 s− 1 is the water vapor diffusivity
through the lag which was measured for the JSC-1 Mars
regolith analog (Hudson et al. 2007) at temperatures between 200 and 263 K and corresponding to a regolith poresize of 3 μm (Liu 2014). The molar mass of water is MH2O =0.018 kg/mol, L is the lag thickness in meters, R = 8.31 J/mol
is the ideal gas constant, Ts is the ice surface temperature, ρis the ice density, and RH is the atmospheric relative humidity. The water vapor saturation pressure (psat) in exchangewith ice is given by the Clausius-Clapeyron equation:
H subl: 1 1
psat ¼ pt exp −
−
R
Ts Tt

ð6Þ

where pt = 611 Pa and Tt = 273 K give the triple point ofwater and Hsubl. = 5.1058 × 104 J mol− 1 is the latent heatof sublimation. Combining Eqs.. 5 and 6 give the sublimation rate for a given ice surface temperature, burial
depth, and atmospheric humidity. We can now associate
an assumed lag thickness of 50 cm under dry (RH = 1%)
atmospheric conditions with an ice surface temperatureof 215 K (further addressed in Mars paleoclimate inputs
section) to an ablation rate of 0.5 mm/year.
Assuming an initial burial depth of 1 cm (the depth atwhich sublimation becomes diffusion-limited (Chevrier
et al. 2008)), we calculated the lag deposit growth (Fig. 7a)and time-evolving sublimation rate (Fig. 7b) for a combination of RH = 1% or 50% and volumetric dust concentrations of 0.1% or 5% in the ice deposit for Ts = 210 K. The
higher of these dust concentrations is associated with ice
purity estimates from SHARAD radar observation ofGemina Lingula in the north polar layered deposit (Grimaet al. 2009) while mid-latitude ice deposits have dust concentrations limited to < 10 vol% (Holt et al. 2008; Petersen
et al. 2018; Plaut et al. 2009). Ice flow simulations used toconstrain the ice rheology of a southern hemisphere
mid-latitude ice deposit suggested dust concentrationsof < 5 vol% (Parsons and Holt 2016).
The sublimation lag calculation shown in Fig. 7a indicates that a lag thickness of ≈ 50 cm would accumulate in
20 kyr or 1 Myr under dry atmospheric conditions with adust concentration of 5% or 0.1%, respectively. The total
ice lost at a given time in Fig. 7 can be calculated by dividing the lag thickness by the dust concentration. Therefore,
the 50 cm thick surface lag assumed in our PISM simulation is associated with 10 m or 500 m of ice loss for dust
concentrations of 5% and 0.1%, respectively.
Assuming the presence of a 50 cm thick lag, Fig. 8 shows
how the ice sublimation depends on the relative humidity
and ice temperature. We allow the surface temperature tovary in the model (see below) which causes the

Page 9 of 21

sublimation rate to change with time, but the relative humidity is assumed to be dry RH = 1% at all times in thesimulation.
The formation of the moraines themselves requiresthat debris be transported on or within the glacial ice inthe absence of basal sliding. Even if the ice is initially exposed to the atmosphere, ice in the ablation zone willdevelop a sublimation lag due to the presence of englacial dust over a timescale that is short compared withthe timescale of FSD formation. The 1 Myr timescale forlag formation with a dust fraction of 0.1 vol% and dry
atmospheric conditions (RH = 1%) is short comparedwith the few million years it takes ice to flow from theaccumulation zone to the inner ridges of the FSD (see“Results and discussion” section).
Ideally, we would hope to model the temporal and spatial
variation in debris lag thickness within PISM as sublimation
occurs and surface debris is advected (perhaps including
dust deposition/erosion as additional parameters), but thiscomputationally and work-intensive sub-model is left as a
topic for future research. For the purposes of this study,
our focus is to provide PISM with a set of conditions representative of the Martian climate and debris layer during
FSD formation which can provide a basis for further investigation into the moraine deposition process. Our mass balance parameterizations were chosen after exploring several
different maximum accumulation and ablation rates to givean ice extent that reaches the interior set of ridges but does
not flow beyond the extent of the FSD during the last 26
Myr of known obliquity variations.
Mars paleoclimate inputs

Constraining the age of the Pavonis FSD has been a difficult task as early crater counts using coarse (∼ 100 m/
pixel) resolution THEMIS imagery failed to recognizethat the thin, ridged unit of the FSD actually superpose
fresh-looking impacts (Shean et al. 2005). Using CTX
imagery, 49 craters which superpose the Pavonis FSD
were counted by Kadish et al. (2014) to give an age estimate of 125 Myr. Mapping of the Pavonis edifice volcanic units (outside the FSD) using HRSC (15 m/pixel)
resolution or better by Gwinner et al. (2011) suggests
the most recent volcanic activity took place 80 ± 20 Myr
ago and could have been contemporaneous with thepresence of ice within the FSD (Shean et al. 2005). Given
the uncertainty in both the timing of FSD formation and
the Martian obliquity history, we rely on the most recent, well-constrained, obliquity history of Mars as aproxy for the past climate responsible for Pavonis glaciation and FSD formation.
The obliquity of Mars over the past 26 Myr derived by(Laskar et al. 2004) (black line in Fig. 9a) was used as an
input to our model in order to determine both the surface
temperature and ice mass balance. We use this input as a

Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 10 of 21

Fig. 7 Time-evolution of lag thickness and ice sublimation rate of a Martian dirty ice deposit. a Modeled increase in lag thickness over time fordifferent atmospheric humidities (solid vs. dashed lines) and ice deposit dust concentrations (black vs. gray lines) assuming an ice surface
temperature of 210 K. The minor ticks on the x-axis give even increments from the base. The lag thickness in a is given by the buried ice
sublimation rate calculated from Eq. 5 (b) integrated over time and multiplied by the dust concentration. This calculation assumed an initialdebris thickness of 1 cm. The blue line in a gives the assumed debris thickness used to calculate the ablation rate in the simulation

proxy for variations experienced earlier in the Martian pastand do not claim that FSD formation took place within thelast 26 Myr. Prior to 26 Myr, the astronomical calculations
used to determine obliquity begin to diverge when slightchanges are made to spin-axis precession rate (Laskar andRobutel 1993; Additional file 1: Figure S1). Rather thanevaluate various obliquity scenarios (which has already been
done by Fastook et al. 2008 to match ice flow simulations
to the maximum FSD extent), our aim is to simulate theinfluence of short-term climate variability on ice flow dynamics in order to investigate the process by which moraines are deposited. Using the obliquity-accumulation
parameterization shown in Fig. 5 for the accumulation zone
and an ablation rate given by Eq. 5 with L = 0.5 m elsewhere, PISM has the necessary mass balance inputs. Next,
a connection needs to be made between obliquity and thesurface temperature of the Pavonis ice sheet.
 at latitude (ϕ) is
The Mars-year average insolation ( Q)given by Ward (1974) which we have modified to approximate the influence of topographic slope:Q¼

S
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2π 1−e2Z 2π qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1−½ sinðϕ þ δ Þ cosθ− cosðϕ þ δ Þ sinψ 2 dψ
0ð7Þ
where S = 588.6 W m− 2 is the solar constant for theMartian orbital semi-major axis, e is orbital eccentricity,θ is obliquity, and ψ is an integration variable. The
north-south component of the local slope is given by δwith positive values associated with pole-facing slopes.
Adjusting the latitude by the local slope is a reasonable
first-order approximation for slopes < 30° based on
full radiative balance calculations by Aharonson andSchorghofer (2006)) although it tends to slightly overestimate the influence of slope on temperature. We
implement this slope-adjusted calculation in order toavoid more demanding computation that would beneeded to produce temperature estimates at the ∼500 m topographic resolution of interest in this work.
 is then used to deterAssuming radiative balance, Q
mine the equilibrium surface temperature (Teq):
T eq ¼ 1=4
ð1−AÞQσð8Þwhere A is surface albedo and σ is the Stefan-Boltzmann
constant. Variation in obliquity and eccentricity over the
past 26 Myr determined by Laskar et al. (2004) areused as inputs into Eq. 7. Although Mars has an average surface albedo of about 0.25, we use a value of 0.3
associated with the dust-covered Tharsis region(Christensen et al. 2001) for locations in the modelParsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 11 of 21

Fig. 8 Dependence of ice ablation on temperature and relative humidity. Contours give sublimation rate for a given temperature (x-axis) andrelative humidity (RH; y-axis) assuming an ice burial depth of 0.5 m based on Eq. 5. The model simulation assumes RH = 0.01 but allows
temperature to vary with time and location along the flow line depending on obliquity and local slope, respectively

experiencing ablation. For locations where ice is accumulating we use an albedo of 0.4 associated withvalues for the northern residual polar cap (Kieffer
et al. 1976) as done in the previous Mars ice sheet
models (Fastook et al. 2008).

An important note about Eq. 7 is that it is independent of the cyclic, ∼ 50 kyr variation in the argument of perihelion even though the instantaneous
incoming solar flux at a given latitude is affected bythis variation. However, Schorghofer (2008) used an

Fig. 9 Obliquity and associated temperature and mass balance model inputs for the past 26 Myr. a Obliquity (black line) determined by Laskar
et al. (2004) and the average surface temperature in the accumulation (red line) and ablation zones (blue line). If ice is accumulating in thespecified accumulation zone, the surface albedo in the zone is changed from 0.3 to 0.4 resulting in a ~ 10 K drop in temperature. b The ablation
rate changes with local temperature (Eq. 5), and ablation occurs everywhere at obliquity less than roughly 44°. The rate of accumulation depends
on obliquity using the assumed relationship shown in Fig. 5

(2020) 7:13

Page 12 of 21

alternative temperature calculation (the “standard
thermal model”) to show that the obliquity signal
dominates over the influence of the argument of perihelion at latitudes outside of ≈ 50–70°N and S.
Figure 9a shows how the average temperature of theablation (red line) and accumulation (blue line) zones
change as a function of time due to changes in obliquity(black line) and eccentricity. At equatorial latitudes, temperatures decrease during high obliquity while the opposite is true at the poles (e.g., Schorghofer 2008). The∼ 5o pole-facing slope in the accumulation zone resultsin colder temperatures than the ablation zone despite
the ablation zone being at slightly higher latitudes.
During times of high obliquity, the temperature in theaccumulation zone is significantly reduced due to theincrease in albedo assumed to be associated with ice
accumulation. The surface temperature inputs to
PISM vary with time according to the obliquity
changes as shown for the “P000” scenario in Additional file 1: Figure S1b, but the surface slope changes
induced by the accumulation and flow of ice are not
updated within PISM during the simulation. Although
the MOLA-based temperatures will closely match
those at the thick, interior portions of the ice sheet,these prescribed surface temperatures will overestimate the temperature at the steeper, pole-facing
margin of the ice sheet (see “Results and discussion”
section).
The computed ice mass balance model input is shownin Fig. 9b where obliquities between 44 and 46° produce
accumulation rates smaller than the 32.4 kg m− 2 maximum (equivalent to 35 mm/year) based on the shape ofthe logistic curve shown in Fig. 5.

Results and discussion
Given our surface mass balance and temperature inputs(Fig. 9), PISM calculates the flow velocities within theice sheet while accounting for the conduction of heat in
the vertical direction, advection of heat in the horizontal
direction, and internal generation of heat from deformation. The simulation produces an ice sheet of varying
extent and thickness (Fig. 10a), but always stays within
the boundaries of the FSD.
Snapshots of the simulation at 25.1 and 24.8 Myr ago
are provided in Figs. 6 and 11 (see Additional files 2, 3, 4
for animated plots of the simulation). The location of
the model flowline is provided by Fig. 6a. The magnitudeof the driving stress, ice temperature, and ice surface velocity are shown in Fig. 11 at these two times. Surface ice
temperatures depressed by the increase in albedo during
times of ice accumulation (right side of Fig. 11c) areconducted to the base of the ice sheet which reduces theice flow rate and requires more ice accumulate (resulting
in elevated stress, Fig. 11a) in order to drive ice flow. In
the ablation zone, an ice thickness of 2.5 km (Fig. 6e) results in a surface to base ice temperatures difference of≈25 K associated with the conduction of the basal heat
flux (Fig. 11d). However, the temperature gradient
within the ice fluctuates by up to 25% due to the propagation of surface temperature variations. Martian regolith thermal modeling by Grott et al. (2007) of 100 kyr
periodicity climate variations with amplitude 10 K give a
10% variation in the Martian surface heat flux (proportional to the temperature gradient) albeit for a regolithwith a lower thermal conductivity than ice. Heat produced from ice deformation is negligible—less than 0.03
mW m− 3 at the base of the ice sheet during/after accumulation events when deformation rates are highest.
Figure 10 provides a visualization of the changes in
conditions at the distal portion (a flowline distance of −100 km is associated with the base of Pavonis Mons) of
the ice sheet over the first 7.5 Myr of the simulation.
The ice thickness increases rapidly during episodes ofaccumulation and slowly thins during prolonged periodsof ablation (Fig. 10a). There are changes to the basal ice
temperature of less than 2 K in amplitude (Fig. 10b)
caused by obliquity-driven temperature variations propagating through the ice deposit. The ice margin velocity
peaks at between 0.1 to 1 m/year during brief periods(Fig. 10c) of ice accumulation but exponentially decreases with time after accumulation ceases due to a reduction in the driving stress as the ice sheet spreads out.
Obliquity-driven temperature fluctuations perturb theice flow slightly during glacial advance (e.g., at − 24Myr). About 2 Myr after the last pulse of ice accumulation, the ice flow rate is less than a millimeter per year
(gray zone in Fig. 10c) and the ice margin begins to retreat as the entire ice sheet continues to sublimate.

Parsons et al. Progress in Earth and Planetary Science

Post-processing of PISM outputs

In order to determine when and where moraine deposition would be most likely from the PISM results, we
saved ice surface velocities from the model every 2 kyr
during the first 4 Myr of our simulation (which starts at− 26 Myr). We then imported these data into a Python
script which tracked the advection of particles placed on
the surface of the ice sheet using a 2D (time and space)linear interpolation.
The location of particle deposition is at the base of
Pavonis with all particles being placed at the − 100 km
point on the flowline (Fig. 1b) after ice advances pastthis point. After the first particle is advected 50 m toward the toe, a new particle is placed at the − 100 km location in a repeating process which mimics thedevelopment of a uniform-thickness sublimation lag.
Once a particle is advected within 1 km of the toe of theice sheet, it is assumed to be deposited and its flowline
location is recorded.

Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 13 of 21

Fig. 10 Model results depicting changes in the ice extent and ice deposit properties in the distal portion of the simulated ice deposit.
Ice thickness (a), basal ice temperature (b), and ice surface flow velocity (c) are indicated by the colored contours plotted against the
distance along the flowline shown in Fig. 6 (left y-axis) and time (x-axis). White regions in these panels contain no ice. The gray region in(c) indicates that ice is present, but the surface flow velocity less than 10− 5 m/year. The black line in each panel gives the obliquity
history over this time period according to the right y-axis with blue filled sections indicating times of ice accumulation in headward
section of the simulated flowline

Surface velocity response to temperature variations

When obliquity is at a local minimum, the surface
temperature reaches a local maximum. Temperature
variations at the base of the ice sheet have a lower amplitude and are out of phase compared with the surface

temperature due to dampening of the propagating, periodic, obliquity-driven thermal wave (e.g., Grott et al.
2007). The thinner ice near the toe of the glacier meansthat the downward propagation of the thermal wave
reaches the base of the ice sooner at the toeParsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 14 of 21

Fig. 11 Simulation results depicting the ice deposit properties and flow velocities at high (47°) and low (30°) obliquities during ice sheet advance.
Panels (a, c, and e give the basal driving stress, temperature, and ice surface velocity, respectively, at high obliquity (− 25.1 Myr) as shown in Fig. 6
a and b. Panels b, c, and f give the same set of properties at low obliquity (− 24.8 Myr)compared with headward portions of the ice sheet.
The right-most dashed vertical line in Fig. 12a helps
to illustrate this effect. In this figure the obliquity signal is given by the black line (right y-axis) with theblue filled portions indicating times of accumulation
at flowline distances less than − 130 km. The dashedline indicates the point in time associated with an obliquity minimum. The wave-like pattern in thetemperature reaches a local maximum at the toe
sooner (closer to the dashed line) than does the headward portion of the ice sheet.
The flow response to this temperature fluctuation
causes the toe of the ice sheet to advance ahead of theheadward portion of the ice sheet, and thin in theprocess. In this case, the increase in velocity is originating at the toe of the ice sheet and propagating in the
headward direction as the low-obliquity temperature
perturbation propagates downward (indicated by theslight rightward skew in the local velocity maximum
relative to the dashed line in Fig. 12b).
Surface velocity response to episodes of ice accumulation

During the series of ice accumulation episodes occurring
during local maxima in obliquity between − 25.6 and −24.9 Myr, accumulation of ice results in a surface velocity
wave propagating from the headward portion of the ice
sheet toward the toe. Once this velocity wave reaches thetoe of the ice sheet, advance occurs rapidly (Fig. 12b).
Particle advection calculation

Computational limitations on the number of advecting
particles we can track (spaced every 50 m) and thenumber of model outputs used (every 2 kyr) pose a challenge to directly simulate ridge formation due to the generation of large file sizes and long computation times.
Although temporal interpolation between outputs andspatial interpolation between model grid points is necessary, the particle-tracking analysis shown in Fig. 13 provides insight into the potential ridge deposition
mechanism.
Plots 13a–f show the density of particles on the ice surface and those that have been deposited on the substrate
as blue and red histograms (500 m bin sizes), respectively.
Between − 25.15(Fig. 13a) to − 18.05 Myr (Fig. 13f) thesimulated ice sheet advances, retreats, and advances again.
The location of ice margin is indicated by the black arrowin each plot. The blue stars indicate debris accumulations
on the glacier surface which, depending on their size, may
result in a ridge once they reach the ice margin. Ridge deposition events associated with a stationary ice margin are
indicated with black stars whereas ridge deposition associated with a kinematic wave generated from an episode ofice deposition is shown with the red star. The black histogram boxes give the deposited particle concentration
using a 50 m bins size but use the same data depicted by
the red histogram.
A consistency check on the kinematic wave deposition
hypothesis comes from estimating the longitudinal distance (Dist) along the ice deposit which surface lag debris of thickness (L) must be collected in order to deposita ridge with a given cross-sectional area (A). Assumingthat the ice margin is stagnant, and that both the laglayer and ridges are ice-free, this “debris collection distance” is:Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 15 of 21

Fig. 12 Increased time resolution model output showing details of the change in ice flow velocity in response to ice accumulation and obliquitydriven temperature change. Highest ice flow velocities follow episodes of ice accumulation. a Basal ice temperature increases with ice thickness(due to thermal conduction of the geothermal flux), while variations over a shorter timescale result from top-down conduction of obliquity
variations. b Surface ice velocity response to episodes of accumulation (blue obliquity peaks) and obliquity-driven temperature changes. Dashedline indicate direction of velocity wave propagation (toeward after accumulation events, e.g., near − 25 Myr, and possibly headward in response
to obliquity-driven temperature change, e.g., at − 24.4 Myr)Dist ¼ A=L:
Using the lag thickness assumed in our model (L =0.5 m) and a mean interior ridge area of 5000 m2(Table 1) gives a distance of 10 km. By comparison, themean width of the interior ridges is 500 m (Table 1).
Therefore, regardless of how much time it takes to deposit a ridge, the ice flow velocity roughly 10 km fromthe ice margin must be 20 times the velocity of ice margin itself in order to deposit the required volume ofmaterial in the space occupied by the ridge. It is also
worth noting, that this discrepancy between the ice
margin velocity and ice flow velocity near the toe must
subsequently disappear in order for ridge deposition tocease. This would happen either by the ice margin velocity increasing (as would happen during ice sheet advance) or by the ice flow velocity decreasing (which

would promote retreat). The latter case represents theformation of terminal moraines, but the former casebetter explains the pattern of interior ridge spacing
resulting from repeating episodes of ice deposition
causing kinematic waves during ice sheet advance.
The changes to the margin of the simulated Pavonis ice deposit using our paleo-climate assumptions
provide several insights into the deposition of ridges
in the FSD. Our model results suggest that ridge formation is a transient process caused by surface velocity responses to quasi-periodic (120 kyr) obliquity
changes. This result is a departure from prior interpretations which suggested ridges were the result of
longer term (~Myr) periods of near steady-state conditions. We revisit the potential ridge-forming mechanisms described in the introduction in the context
of our simulation results:

Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 16 of 21

Fig. 13 Ridge formation modeled by advection of tracer particles placed on the simulated ice sheet. Blue and red histograms give the number oftracer particles on the ice sheet surface and deposited on the substrate, respectively using a 500 m bin size (resolution of PISM simulations).
Particles are individually placed at the − 100 km flow line position and are advected by the flowing ice. Once a particle has been transported 50
m, a new particle is placed at the − 100 km position. The black histogram gives the number of deposited particles using a 50 m bin size andrepresents the same data in the red histogram.

Conditions during ice margin retreat

Ice sheet thinning results in the cessation of ice flow
during ice margin retreat even during times of low obliquity (warm ice conditions) (Fig. 10b, c). Therefore, advection of surface debris has halted during ice margin
retreat and deposition of the surface lag can only occur
passively as sublimation removes ice from below the thindebris layer (as depicted in Fig. 13c). Therefore, using
the low sublimation rate and obliquity-driven accumulation assumptions in this work (Figs. 5, 8, 9), we find itvery unlikely that FSD ridges were formed during glacial
retreat.
Temperature variations during ice margin advance

Ice accumulation and, to a lesser extent, basal
temperature changes affect the ice flow velocity in our
simulations. During advance, the average ice sheet margin velocity is initially ~ 10 cm/year but slows to about
0.5 cm/year after a million years of flow (Fig. 12). Ice
surface temperatures vary due to obliquity changes,
which, after being conducted to the base of the ice sheet,produce temperature changes of less than 2 °C (Fig 10
and 12). These subtle temperature changes have little effect on the ice flow dynamics (Fig. 12b) and cannot produce significant kinematic waves that would perturb the
surface debris layer. However, because our surface
temperature calculation is based on basal topography
and is not updated to reflect the topography of the ice
deposit, the simulation may be underestimating the roleof temperature changes on the ice flow dynamics. Even

more important is the consequence that this assumption
has on the ablation rate (discussed below) which ishighly temperature dependent. Therefore, we cannot
rule out the potential for temperature variations to cause
variations in the ablation rate which would lead to theperiodic deposition of moraines. Our simulations show,
however, that such a process would need to happen during glacial advance in order for ice flow to advect debris
to the ice margin which would be held in place during
times of high ablation.
Renewed accumulation

A pulsating advance of the ice margin results from several episodes of accumulation during short-lived periodsof high obliquity associated with the 120 kyr periodicity
in obliquity (Fig. 12). Discrete episodes of accumulation
generate kinematic waves in the ice deposit which produce thickness variations in the advecting debris layer(blue stars in Fig. 13a). However, during the initial advance at − 26 Myr, the amplitude of these accumulationdriven velocity perturbations is not significant enough toresult in narrow zones of high debris thickness thatwould result in the formation of a ridge upon depositionat the margin. The − 100 km particle deposition location
is likely responsible for this result because the fluctuation in ice velocity is much greater in the accumulation
zone (and would promote greater variation in debris
thickness) than at the base of Pavonis. However, in orderto provide time for a surface lag to develop, we assume aparticle deposition location which is about 30 km

Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

Page 17 of 21

downslope from the edge of the high-obliquity accumulation zone.
Despite this assumption, a thick, narrow pulse of surface debris does form at about − 20 Myr when a stagnant
toe is reactivated by renewed accumulation producing a
propagating velocity “step” (Fig. 10c). In our simple particle advection analysis, this step in velocity causes particles to “pile-up” in front of the high-amplitude
kinematic wave (blue star in Fig. 13d) and results in a
local deposition of surface particles nearly 20 times moreconcentrated than that given by the assumed surface
particle spacing of 50 m (black bars below arrow in
Fig. 13e). The formation of this ridge is contingent upon
the assumption that the re-activated ice sheet does notoverride the stagnant ice, but “flows into” the distal deposit in a way that advects the surface debris in processsimilar to a bulldozer.
Surface debris slope failure at ice margin

The steepness of the margin of flowing ice sheets on
Mars is expected to be greater than their terrestrial
counterparts due to reduced gravity. This conclusion
can be reached by analyzing the shear stress which, atthe base of the ice sheet of thickness (H)—given by:τ ¼ ρgH

∂h
∂x

ð9Þ

for the 2-D, flowline case. In order for Martian ice sheets
to achieve comparable basal shear stresses to terrestrial
ice sheets, Eq. 9 shows that the ice surface slope (∂h/∂x)
must be greater on Mars to compensate for the reduction in g. Identical simulations performed under Earthgravity resulted in a decrease in the ice margin slope
from roughly 15° under Martian gravity to 10° (see Additional files 1, 2, 3, 4).
However, the simulated ice margin slopes are still lessthan the angle of repose for unconsolidated dust on
Mars of 28–35° (Atwood-Stone and McEwen 2011) andthe lateral continuity of FSD ridges suggest slope failureof surface debris is not the primary process responsible
for ridge formation.
A cartoon depicting our supported hypothesis isshown in Fig. 14. Initial sublimation of a dirty ice deposit
results in the formation of a surficial lag layer (Fig. 14a).
During high obliquity, accumulation of ice produces akinematic wave which accumulates surface debris(Fig. 14b) and subsequently deposits it at the ice margin
as the ice sheet advances (Fig. 14c). This process would
then repeat during the next obliquity cycle. The spacing
between ridges would depend on the timing and magnitude of ice accumulation events which dictate how quickly
the ice margin advances. Rapid advance would ensue soon

Fig. 14 Sketch of proposed mechanism for ridge deposition. a
Sublimation lag develops on surface of ice sheet. b Ice accumulationon the flanks of Pavonis produces a velocity wave which collects
surface debris. c Debris is deposited at the ice margin andsubsequently overridden by the advancing ice sheet

after multiple ice accumulation events but would exponentially decrease over time as shown in Fig. 10.
After deposition, the ridges would have become overridden by the advancing ice but would have remained
preserved due to the lack of basal sliding (as indicatedby intersecting ridges such as those in Fig. 1c and thesuperposed smooth unit on ridges in Fig. 1d). We conclude that the simulation results support an alternativeinterpretation of ridge sequences as processional (as opposed to recessional) drop moraines.
Model limitations

The simulation results are subject to our choice ofmodel inputs and the assumptions made within PISM.

Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

Our mass balance and temperature calculations are supported by GCM results (Forget et al. 2006) and insolation calculations (Schorghofer 2008; Ward 1974) for
Amazonian Mars. The choice to use the recent obliquity
history stems from uncertainties in the Martian obliquity
prior to 26 Myr ago, so we are focusing on the best constraints currently available for the Amazonian climate.
Because ice in our simulation does not reach the full extent of the FSD, there was likely a longer episode, ormore frequent episodes, of ice accumulation earlier inthe Amazonian. Alternatively, the ablation rate may havebeen lower than the ~ 0.5 mm/year rate associated witha 0.5 m thick debris layer used in our simulation.
However, our goal is not to reproduce the extent of
the FSD, but to simulate the process by which ridges interior to the FSD were deposited. By focusing the past
26 Myr, we are able to show that known climate variations can plausibly result in ridge deposition, although
the ridges themselves may be older than this recent obliquity history.

Basal heat flux and heat transport

Because the Tharsis region is likely subject to a higher
than average heat flux, the model assumes a value of 30
mW m− 2 based on thermal models presented in Plesa
et al. (2016). The primary influence of the basal heat flux
is on the temperature at the base of the ice sheet whichacts to control the ice sheet thickness. In our simulations, higher heat fluxes result in warmer basal temperatures and more rapid ice deformation which produces amore extensive and thinner ice sheet. Identical Pavonis
simulations run with a fixed ice volume using heat fluxesof 20 and 50 mW m− 2 give ice thicknesses and extents
which are modified by a factor of 0.85 and 1.1, respectively after flow has nearly ceased (after 3 Myr) for thehigher heat flux case. Therefore, the simulation presented here (e.g., Fig. 9) may be biased toward slightly
slower flow and a somewhat thicker ice sheet. Thinner
ice sheets will respond even more dramatically to theobliquity-driven temperature signal because periodic
temperature variations tend to become dampened withdepth (for example, Fig. 10 and Grott et al. 2007).
Therefore, the conclusions regarding the ridge deposition process inferred from the simulation presented
here would still apply to a somewhat thinner ice sheet
resulting from a higher basal heat flux.
Although PISM provides important insights into the
dynamics of Martian ice sheets and captures much of
the essential physics, there are aspects of the Martianenvironment that are not fully represented in the model.
For instance, the model assumes that lateral heat transport is dominated by advection and that conduction in
the horizontal direction is negligible. An ice thermal

Page 18 of 21

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
conduction lengthscale ðd ¼ kt=ðρcp ÞÞ of 2.4 km is associated with an obliquity timescale (t) of 120 kyr. Dividing the lengthscale by the timescale gives an estimate forthe velocity of a thermal pulse through a Martian ice deposit (equal to 2 cm/year) which is similar to modeled
ice flow rates (e.g., Fig. 11b, e). While it is reasonable toneglect heat conduction in the horizontal direction for aterrestrial ice sheet where ice velocities are high (in partdue to sliding of ice across the bed), the combination ofslow flow and a higher ice thermal conductivity (Eq. 4)
at the low temperatures on Mars make the heat flux inthe horizontal direction via conduction comparable withthat of advection based on the rough calculation above.

Ice surface temperature and mass balance

Secondly, the surface temperature model inputs rely on
the slope of the basal topography rather than the timevarying ice surface topography. Therefore, the influenceof the changing ice surface slope on the surface
temperature is not captured by the simulation and may
be important at the toe of the ice sheet where slopes aresteep and pole-facing (in the case of the Pavonis FSD).
In particular, the influence of the surface slope on ice
temperature may help explain the larger size of ridges
found within the Pavonis FSD. As ice advanced northward across the Pavonis FSD, ice temperatures at the toe
may have varied markedly as the angle of incident sunlight changed with obliquity—resulting in more accentuated episodes of advance during high obliquity followed
by a cold, stagnant margin during low obliquity. In thisrespect, PISM may be better suited to simulating the
westward ice flow which occurred in the Arsia FSD
where the margin slope orientation would have less influence on temperature which may account for thesmaller, more closely spaced ridges found there.
The surface temperature also plays a key role in determining the ice sublimation rate (Eqs. 5 and 6). As notedin the ice ablation section, the low (~ 0.5 mm/year) ice
sublimation rates calculated in our numerical simulation
is associated with a debris cover thickness of 0.5 m. If,
however, the ice temperature increased from 215 to 220
K, the debris thickness required to maintain a 0.5 mm/year sublimation rate would be 0.9 m. If, however, thedebris thickness is maintained at 0.5 m, a 5 K temperature
increase (to 220 K) would nearly double the sublimation
rate. Therefore, although the simulation rate is nearly constant, an increase in ice loss is likely during times of highobliquity when the pole-facing ice margin is warm. Whatis uncertain is whether the influence of temperature on
sublimation would outweigh the influence of relative humidity which one might expect to be higher during the redistribution of polar ice at high obliquity. It is possiblethat the increase in sublimation caused by warm

(2020) 7:13

Page 19 of 21

temperatures on the pole-facing ice margin during highobliquity is canceled by the reduction in the rate of sublimation caused by an increase in humidity (Fig. 5b). The
complicated relationship between obliquity and spatial
variations in the mass balance of equatorial ice on Mars is
an important topic for future research. In general, increasing the average sublimation rate in our simulations resultsin fewer ridges which are in close proximity to the distaledge of the accumulation zone (assumed to be 5 km elevation) while lower sublimation rates produce a longer-lived,
more extensive ice deposit which can produce moreridges.
Disregarding other factors (atmospheric deposition of
volcanic ash or dust) that may influence the debris layer
thickness, the factor of ~ 20 increase in particle concentration produced from accumulating surface debris in a
kinematic wave suggests that the original sublimation
lag thickness can be very roughly estimated by dividing
the height of the ridge by 20. However, this “particle
concentration factor” will likely change depending on location—increasing as the ice sheet becomes more extensive and debris can be collected over a larger area. Most
interior ridges have heights of < 40 m suggesting an original surface debris thickness of < 2 m which would takeless than 200 kyr to form from the sublimation of a 5
vol% dusty ice deposit (Fig. 7a). Subsequent to surface
lag formation, renewed ice accumulation resulted in advance of the glacial front and ridge deposition. Note that
our conclusion that ridge deposition occurs during glacial advance does not run counter to the geologic interpretations of Shean et al. (2005), but does indicate that
ice retreat is not needed to form the ridges.

ablation rate we employ results in dramatic glacial advance during times of accumulation at high obliquity
which gradually slows, stops, and slowly retreats duringextended periods of lower obliquity (Fig. 10). During retreat, ice thicknesses are greatly reduced and are notcapable of driving flow resulting in a “dead” ice deposit(Fig. 10d). The surface debris on top this ice would be
passively dropped and is unlikely to form the long setsof parallel ridges which define a former ice margin.
Instead, we propose that ridges are deposited during
glacial advance when precipitation-induced perturbations to the surface ice velocity (Fig. 12b) concentratesurface debris in front of a kinematic wave (Fig. 13c, d,e). Ridges deposited during advance would become buried by the advancing ice sheet but would remain preserved due to the lack of basal sliding of the ice sheet(see “Results and discussion” section). Observational evidence for ridge preservation below an ice sheet comesfrom intersecting ridge sequences at Arsia (Fig. 1c) andPavonis (Shean et al. 2005), as well as the superpositionof remnant ice (the smooth unit) on top of ridges(Fig. 2b). The periodicity in the ridge spacing (Fig. 3a) is
also indicative of an advancing ice sheet which experiences brief periods of ice accumulation causing the ice
margin velocity temporarily increase, then gradually slow
as the ice thins during advance. This mechanism forridge emplacement is contrary to the recessional drop
moraine hypothesis previously ascribed to these ridges(Fastook et al. 2008; Head and Marchant 2003; Kadish
et al. 2014; Shean et al. 2005) and suggests that the formation of individual ridges takes place over a shorter
timescale than previously thought.
We conclude that initial sublimation of the uppermost
portions of the ice deposit likely resulted in the formation of a debris lag which protected ice from further
ablation—allowing it to slowly advance and deposit
moraines as its flow velocity periodically fluctuated withobliquity. This new interpretation of the FSD ridge formation process provides an opportunity to further constrain the Amazonian climate history of Mars.

Parsons et al. Progress in Earth and Planetary Science

Conclusions
Thermomechanical ice sheet simulation of the Amazonian
Pavonis ice sheet using obliquity inputs from the past 26
Myr provides important insights into the formation of aseries of ridges in the interior of the Pavonis FSD. Our
simulation assumes an obliquity-precipitation relationshipgiven in Fig. 5 at elevations between 5 and 12 km and ablation rates near 0.5 mm/year in places where ice is not accumulating. This ablation rate is associated with ice buried bya debris lag ≈ 50 cm thick which would form in less than 1
Myr from the sublimation of ice containing > 0.1 vol% dust
under low humidity conditions (Fig. 7). The model calculates heat transport in the vertical direction associated witha basal heat flux of 30 mW m− 2, deformation-induced
strain heat, and obliquity-driven surface temperature variations while advection accounts for the horizontal heattransport. Horizontal ice flow rates are calculated using
Glen’s Law (Eq. 3) under the shallow ice approximation.
The reduced ablation rate used in our simulation is a
key departure from previous Tharsis ice sheet simulations performed by Fastook et al. (2008). The low

Supplementary information
Supplementary information accompanies this paper at https://doi.org/10.
1186/s40645-020-0323-9.
Additional file 1: Figure S1. Obliquity History Five representative
scenarios for the obliquity of Mars determined by Laskar et al. 2004 overthe past 250 Myrs (a) and 26 Myrs (b). The different colored lines refer todifferent assumed values for the initial precession rate of the rotation axisof Mars which produces chaotically different results over long timescales(b), but nearly identical results in the past 26 Myrs (a). Thick lines give
obliquity over a 1 Myr averaging window. Scenario "P000" is used as ourmodel input.
Additional file 2. Pavonis Simulation 20 mW/m2. Animated plots of theice flow simulation which uses a basal heat flux of 20 mW/m2 showing

Parsons et al. Progress in Earth and Planetary Science(2020) 7:13

the extent, thickness, basal stress, approximate heat flux, ice temperature,
ice effective viscosity, and ice surface velocity.
Additional file 3. Pavonis Simulation Earth Gravity 20 mW/m2.
Animated plots of the ice flow simulation using Earth-like gravity and abasal heat flux of 20 mW/m2 showing the extent, thickness, basal stress,
approximate heat flux, ice temperature, ice effective viscosity, and ice
surface velocity.
Additional file 4. Pavonis Simulation 30 mW/m2. Animated plots of theice flow simulation presented in the text which uses a basal heat flux of
30 mW/m2 showing the extent, thickness, basal stress, approximate heat
flux, ice temperature, ice effective viscosity, and ice surface velocity.
Abbreviations
CTX: Context camera; DEM: Digital elevation model; FSD: Fan-shaped deposit;
MOLA: Mars Orbiter Laser Altimeter; PISM: Parallel ice sheet model
Acknowledgements
RAP is thankful to Fitchburg State University for granting a leave of absence
that made this research possible and to the University Museum, University ofTokyo, for providing computing resources for running simulations. We
deeply appreciate the work undertaken by Andy Aschwanden, Ed Bueler,
Constantine Khroulev, and the other members of the PISM development
team in documenting and sharing PISM. Feedback from an anonymous
reviewer and insightful comments from A. Grau Galofre greatly improved themanuscript.
Authors’ contributions
TK proposed the topic and carried out initial ridge topography observations.
RAP conducted numerical simulations, final ridge topography profiles, andmanuscript drafting. RH constructed the stereo DEM and helped edit themanuscript. HM aided in the interpretation of results. All authors read andapproved the final manuscript.
Authors’ information
RAP conducted this research while on leave from Fitchburg State University.
Funding
Partial support for this research was provided by a grant NASA’s Mars Data
Analysis Program (MDAP 80NSSC18K0206). Development of PISM is
supported by NASA grant NNX17AG65G and NSF grants PLR-1603799 andPLR-1644277.
Availability of data and materials
Simulation results and CTX stereo elevation data is available upon requestfrom RAP.
Competing interests
The authors declare that they have no competing interests.
Author details
1
University Museum, The University of Tokyo, 7 Chome-3-1 Hongo, Bunkyo,
Tokyo 113-0033, Japan. 2Earth and Geographic Sciences, Fitchburg State
University, 160 Pearl St., Fitchburg, MA 01420, USA. 3Department of Systems
Innovation, The University of Tokyo, 7 Chome-3-1 Hongo, Bunkyo, Tokyo
113-0033, Japan.
Received: 17 February 2019 Accepted: 5 February 2020

References
Aharonson O, Schorghofer N (2006) Subsurface ice on Mars with roughtopography. J Geophys Res 111:10
Atwood-Stone C, McEwen AS (2011) Avalanche slope angles in low-gravity
environments from active Martian sand dunes. Geophys Res Lett 40:2929–2934
Beyer RA, Alexandrov O, & McMichael S (2018) The Ames Stereo Pipeline: NASA's
open source software for deriving and processing terrain data. Earth andSpace Science 5:537–548. 10.1029/2018EA000409.
Bockheim JG (2002) Landform and soil development in the McMurdo Dry Valleys,
Antarctica: a regional synthesis. Arctic Antarct Alp Res 34:308–317 Taylor
FrancisPage 20 of 21

Bockheim JG (2010) Soil preservation and ventifact recycling from dry-based
glaciers in Antarctica. Antarct Sci 22:409–417 Cambridge University Press
Budd WF, Jacka TH (1989) A review of ice rheology for ice sheet modelling. ColdReg Sci Technol 16:107–144
Bueler E, Brown J, Lingle C (2007) Exact solutions to the thermomechanicallycoupled shallow-ice approximation: effective tools for verification. J Glaciol 53
Carr MH, Greeley R, Blasius KR, Guest JE, Murray JB (1977) Some Martian volcanicfeatures as viewed from Viking orbiters. J Geophys Res 82:3985–4015
Chevrier V, Ostrowski DR, Sears DWG (2008) Experimental study of thesublimation of ice through an unconsolidated clay layer: implications for thestability of ice on Mars and the possible diurnal variations in atmosphericwater. Icarus 196:459–476
Chevrier V, Sears DWG, Chittenden JD, Roe LA, Ulrich R, Bryson K et al (2007)
Sublimation rate of ice under simulated Mars conditions and the effect oflayers of mock regolith JSC Mars-1. Geophys Res Let 34:10
Cuffey KM, Paterson WSB (2010) The physics of glaciers. Elsevier Science Available
from: https://books.google.co.jp/books?id=Jca2v1u1EKEC. [cited 2010]
Christensen PR, et al. (2001) Mars Global Surveyor Thermal Emission
Spectrometer experiment: Investigation description and surface science
results, J Geophys Res 106(E10):23823–23871, https://doi.org/10.1029/
2000JE001370.
Dundas CM, Byrne S (2010) Modeling sublimation of ice exposed by new
impacts in the Martian mid-latitudes. Icarus 206:716–728
Durham WB, Prieto-Ballesteros O, Goldsby DL, Kargel JS (2010) Rheological andthermal properties of icy materials. Space Sci Rev 153:273–298
Fastook JL, Head JW, Marchant DR, 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
Forget F, Haberle RM, Montmessin F, Levard B, Head JW (2006) Formation of glacierson Mars by atmospheric precipitation at high obliquity. Science 311:368–371
Francis PW, Waadge G (1983) The Olympus Mons Aureole: formation by
gravitational spreading. J Geophys Res 88:8333–8344
Goldsby D, Kohlstedt D (2001) Superplastic deformation of ice: experimental
observations. J Geophys Res 106:11017–11030
Greve R (1997) Application of a polythermal three-dimensional ice sheet modelto the Greenland ice sheet: response to steady-state and transient climate
scenarios. J Clim 10:901–918
Greve R, Blatter H (2009) Dynamics of ice sheets and glaciers. Springer,
Heidelberg
Grott M, Helbert J, Nadalini R (2007) Thermal structure of Martian soil and themeasurability of the planetary heat flow. J Geophys Res 112:09004
Gwinner K, Head JW, Wilson L, Fassett CI, Dissmore S (2011) Surface ages of the
volcanic deposits of Pavonis Mons and implications for the magma supply ofTharsis. Lunar Planet Sci Conf 42:2466
Hartmann WK (2005) Martian cratering 8: isochron refinement and thechronology of Mars. Icarus 174:294–320
Head JW, Marchant DR (2003) Cold-based mountain glaciers on Mars: Western
Arsia Mons. Geology 31
Head JW, Mustard JF, Kreslavsky MA, Milliken RE, Marchant DR (2003) Recent ice
ages on Mars. Nature 426:797–802
Hemmi R, Miyamoto H (2017) Distribution, morphology, and morphometry ofcircular mounds in the elongated basin of northern Terra Sirenum, Mars.
Prog Earth Planet Sci 4:26
Holt JW, Safaeinili A, Plaut JJ, Head JW, Phillips RJ, Seu R et al (2008) Radarsounding evidence for buried glaciers in the southern mid-latitudes of Mars.
Science 322:1235–1238
Hudson TL, Aharonson O, Schorghofer N, Farmer CB, Hecht MH, Bridges NT(2007) Water vapor diffusion in Mars subsurface environments. J Geophys
Res Planets 112:E05016. https://doi.org/10.1029/2006JE002815
Ingersoll AP (1970) Mars: occurrence of liquid water. Science 168:972–973
Jakosky BM, Carr MH (1985) Possible precipitation of ice at low latitudes ofMars during periods of high obliquity. Nature 315:559 Available from:https://doi.org/10.1038/315559a0. Nature Publishing Group SN. [cited
1985 Jun 13]
Kadish SJ, Head JW, Fastook JL, Marchant DR (2014) Middle to late Amazonian
tropical mountain glaciers on Mars: the ages of the Tharsis Montes fansshaped deposits. Planet Space Sci 91:52–59
Kadish SJ, Head JW, Parsons RL, Marchant DR (2008) The Ascraeus Mons fanshaped deposit: volcano–ice interactions and the climatic implications ofcold-based tropical mountain glaciation. Icarus 197:84–109

(2020) 7:13

Page 21 of 21

Kieffer HH, Chase CC, Martin TZ, Miner ED, Palluconi FD (1976) Martian North
Pole Summer Temperatures: Dirty Water Ice, Science, 194(4271):1341–1344.
Laskar J, Correia ACM, 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
Laskar J, Robutel P (1993) The chaotic obliquity of the planets. Nature 361:608
Available from: doi: 10.1038/361608a0. Nature Publishing Group SN -[cited
1993 Feb 18]
Liu L (2014) Thermal regime and sublimation rates of subsurface ice in Antarctica
and on Mars under current and past environmental conditions: insights frommodeling and ground data. Univ. of Washington, Seattle
Lucchita BK (1984) Ice and debris in the fretted terrain, Mars. Lunar Planet Sci
Conf XIV(suppl):409
Lucchitta BK (1981) Mars and earth: comparison of cold climate features. Icarus
45:264–303
Marchant DR, Denton GH, Bockheim JG, Wilson SC, Kerr AR (1994) Quaternarychanges in level of the upper Taylor glacier, Antarctica: implications forpaleoclimate and East Antarctic ice sheet dynamics. Boreas 23:29–43
Marshall SJ, Tarasov L, Clarke GKC, Peltier WR (2000) Glaciological reconstructionof the Laurentide ice sheet: physical processes and modelling challenges.
Can J Earth Sci 37:769–793
Mischna MA, Richardson MI, Wilson RJ, McCleese DJ (2003) On the orbital forcingof Martian water and CO2 cycles: a general circulation model study withsimplified volatile schemes. J Geophys Res 108:5062. https://doi.org/10.1029/
2003JE002051, E6
Parsons R, Holt J (2016) Constraints on the formation and properties of a Martian
lobate debris apron: insights from high-resolution topography, SHARAD radar
data, and a numerical ice flow model. J Geophys Res 121:432–453
Parsons RA, Nimmo F, Miyamoto H (2011) Constraints on martian lobate debris
apron evolution and rheology from numerical modeling of ice flow. Icarus
214:246–257
Paterson WSB, Budd WF (1982) Flow parameters for ice sheet modeling. ColdReg Sci Tech 6:175–177
Pathare AV, Paige DA (2005) The effects of martian orbital variations upon the
sublimation and relaxation of north polar troughs and scarps. Icarus 174:419–443
Payne A (1999) A thermomechanical model of ice flow in West Antarctica. Clim
Dyn 15:115–125
Petersen EI, Holt JW, Levy JS (2018) High ice purity of Martian lobate debris
aprons at the regional scale: evidence from an orbital radar sounding survey
in Deuteronilus and Protonilus Mensae. Geophys Res Lett 45:11–595
Pettit EC, Waddington ED (2003) Ice flow at low deviatoric stress. J Glaciol
49(166), 359–369. https://doi.org/10.3189/172756503781830584
PISM authors The PISM, a Parallel Ice Sheet Model. 2014. Available from: http://
www.pism-docs.org. [cited 2014]
Plaut JJ, Safaeinili A, JW Holt , 2009. Radar evidence for ice in lobate debris
aprons in the mid-northern latitudes of Mars. 2009. Geophysical; Available
from: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2008GL036379.
[cited 2009]
Plesa A-C, Grott M, Tosi N, Breuer D, Spohn T, Wieczorek MA (2016) How large
are present-day heat flux variations across the surface of Mars? J Geophys
Res 121:2386–2403
Scanlon KE, Head JW, Wilson L, Marchant DR (2014) Volcano-ice interaction in the
Arsia Mons tropical mountain glacier deposits. Icarus 237:315–339
Schorghofer N (2006) Dynamics of ice ages on Mars. Nature 449:192–194
Available from: https://www.nature.com/articles/nature06082. [cited 2006]
Schorghofer N (2008) Temperature response of Mars to Milankovitch cycles.
Geophys Res Lett 35:L18201. https://doi.org/10.1029/2008GL034954
Scott DH, Dohm JM, Zimbelman JR (1998) Geologic map of Pavonis Mons
volcano, Mars
Scott DH, Tanaka KL (1986) Geologic map of the western equatorial region of
Mars
Shean DE, Head JW, Marchant DR (2005) Origin and evolution of a cold-based
tropical mountain glacier on Mars: the Pavonis Mons fan-shaped deposit. J
Geophys Res 110:E05001. https://doi.org/10.1029/2004JE002360
Staiger JW, Marchant DR, Schaefer JM, Oberholzer P, Johnson JV, Lewis AR et al(2006) Plio-Pleistocene history of Ferrar Glacier, Antarctica: implications forclimate and ice sheet stability. Earth Planet Sci Lett 243:489–503
Tanaka KL (2004) Geology and insolation-driven climatic history of Amazonian
north polar materials on Mars. Nature 437:991–994 Available from: https://www.nature.com/articles/nature04065. [cited 2004]

Ward W (1974) Climatic variation on Mars. I. Astronomical theory of insolation. J
Geophys Res 79:3375–3386
Ward WR (1992) Long-term orbital and spin dynamics of Mars. Mars. Univ. of Ariz.
Press, Tuscon, pp 298–320
Williams RS (1978) Geomorphic processes in Iceland and Mars: a comparative
appraisal from orbital images. Geol Soc Am 10:517
Winkelmann R, Martin MA, Haseloff M, Albrecht T, Bueler E, Khroulev C et al(2011) The Potsdam Parallel Ice Sheet Model (PISM-PIK) part 1: model
description. Cryosphere 5:715–726 Available from: http://www.thecryosphere.net/5/715/2011/tc-5-715-2011.pdf. [cited 2011]
Zimbelman JR, Edgett KS (1992) The Tharsis Montes, Mars: comparison of
volcanic and modified landforms. Lunar Planet Sci Conf 22:31–44

Parsons et al. Progress in Earth and Planetary Science

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.