Polyphase Mid‐Latitude Glaciation on Mars: Chronology of the Formation of Superposed Glacier‐Like Forms from Crater‐Count Dating

Reconstructing Mars's glacial history informs understanding of its physical environment and past climate. The known distribution of viscous flow features (VFFs) containing water ice suggests that its mid‐latitudes were glaciated during the Late Amazonian period (the last several hundred million years). The identification of a subgroup of VFFs—called superposed glacier like forms (SGLFs)—flowing onto other VFFs, indicates multiple glacial phases may have occurred during this time. To explore the history and spatial extent of these glaciations, we record the distribution of SGLFs globally and use impact‐crater counting to date the SGLFs and the VFFs onto which they flow. Our inventory expands the handful of SGLFs reported in earlier literature to include 320 located throughout the mid‐latitudes. Our dating reveals these SGLFs to be much younger than their underlying VFFs, which implies a spatially‐asynchronous glaciation. SGLFs have been forming since ∼65 Ma, and their ages are clustered in two distinct groups around 2–20 and 45–65 Ma, whereas the ages of their underlying VFFs span the last ∼300 Ma diffusely. We discuss these results in the light of well‐known uncertainties with the crater‐dating method and infer that while ice sheets decayed over the Late Amazonian period, alpine glaciers waxed and waned in at least two major cycles before their final demise approximately two million years ago.


Introduction
Mars has substantial volumes of water ice held in its polar ice caps (Laskar et al., 2002;Levrard et al., 2007;Smith et al., 2016). Beyond the polar regions, water ice can only persist if shielded by overlying debris deposits, as exposed ice sublimates quickly. Of interest in this regard are a distinct group of buried ice deposits located throughout the mid-latitudes (Squyres, 1978) called viscous flow features (VFFs: Milliken et al., 2003). Numerous studies have suggested that VFFs are analogous to glacial systems on Earth, positing a similar origin for VFFs whereby ice accumulates following precipitation over several years Dickson et al., 2008;Head et al., 2010;Levy et al., 2007;Morgan et al., 2009). Initially, only small lobate features emerging from alcoves were referred to as VFFs . However, the definition has since been revised and following Souness et al. (2012), we use VFF as an umbrella term encompassing a range of landforms subdivided according to their size and context. The highest order form of VFFs are lineated valley fill (LVF) comprising highly-integrated anastomising flows on valley floors, which extend continuously for hundreds of km (Baker et   system-emanating from alcoves, converging, coalescing, and continuing downslope before terminating in steep convex lobes-as observed in inter-montane glacial settings on Earth (Dickson et al., 2008;Levy et al., 2007Levy et al., , 2014; Morgan et al., 2009). Impact crater morphology (Kress & Head, 2008) and radar-detected stratigraphy (Holt et al., 2008;Plaut et al., 2009) suggest that VFFs consist of a predominantly water-ice substrate covered by a thin surficial debris layer, akin to terrestrial debris-covered glaciers (Baker & Carter, 2019). It is worth also considering the relationship between VFFs and what has been termed latitude dependent mantle (LDM: Kreslavsky & Head, 2002)-a widespread terrain-draping morphological unit found >45 • north and south of the equator that is inferred to contain ice and dust and covers ∼23% of Mars's surface. Gamma-ray spectrometer data suggest LDM contains the bulk of Mars's shallow ground ice (Feldman et al., 2004). However, disagreement persists with respect to the formation of LDM and the extent to which VFF deposits are overlain by and/or sourced from these deposits. Diffusion of atmospheric vapour into regolith is one formation mechanism proposed (Mellon & Jakosky, 1995), but this alone cannot explain observations of ice content in decametre-thick LDM deposits exceeding 45-90% (Conway & Balme, 2014). We follow the interpretation that thick LDM is more consistent with sustained atmospheric deposition of snow/frost (subsequently shielded by debris) during extensive mid-latitude glaciation events (Berman et al., 2015;Conway & Balme, 2014;Head et al., 2003;Milliken et al., 2003;Mustard et al., 2001). Thick LDM deposits with high ice content are expected to deform plastically and therefore may be indistinguishable from some smaller VFFs (e.g., GLFs) in form and appearance (Conway & Balme, 2014). LDM may also contribute to the debris cover on VFFs.
VFFs are distributed within the 30-50 • latitude zones in each hemisphere on Mars (Levy et al., 2007(Levy et al., , 2014Milliken et al., 2003;Souness et al., 2012), inconsistent with ice emplacement under current surface conditions. This finding motivates the so-called "mid-latitude glaciation" hypothesis (Head et al., 2003, 2005, 10.1029/2019JE006102 2006. Mars's orbital parameters are too chaotic to be definitively modelled beyond 20 Ma (Laskar et al., 2004). However, these authors found, through statistical modelling of 20,000 possible orbital solutions, that Mars's obliquity (25.19 • today) has likely fluctuated strongly about its mean value (37.62 • ) over the last 4 Gyr, inferring that, by modulating insolation receipt, these chaotic fluctuations drive extensive zonal redistribution of volatiles in the martian atmosphere. The mid-latitude glaciation hypothesis sees these shifts as driving climatic cycles, causing glaciations at times of high average obliquity (≥35 • -40 • ) when volatiles were redistributed from the polar ice sheets to nourish mid-latitude ice (Head, Nahm, et al., 2006;Laskar et al., 2004). This idea has received support from global climate models that predict precipitation in VFF-rich regions during high obliquity events (Forget et al., 2006;Madeleine et al., 2009), and from glacial flow modelling that predicts kilometres thick ice sheets extending into the plains abutting the present VFF-rich regions in such scenarios (Fastook et al., 2011).
Accordingly, the mid-latitude glaciation hypothesis regards VFFs not as active glaciers, but as the recessed relics of past glaciation armoured by surface debris (Head et al., 2003. Complex extensive ice flow emanating from plateau ice sheets fits the morphology of LVFs observed today Head et al., 2010;Head, Nahm, et al., 2006), and their armouring by debris is thought to have occurred when these ice sheets collapsed (Dickson et al., , 2008. Modelling by Fastook et al. (2014) confirms that during such a deglaciation, VFFs can form as scarps and massifs become exposed and deposit debris onto VFF surfaces, and that over time, the deposition modulates their distal surface slope such that the flow of exposed ice ensures their entire length becomes debris-covered. Sublimation further consolidates surface deposits into a barrier that retards subsequent mass loss, allowing VFFs to persist to the present day (Fastook et al., 2014;Kowalewski et al., 2006;Levy et al., 2007). Similar long-term survival of buried ice has been observed in the debris-covered glaciers of the Antarctic Dry Valleys on Earth (Kowalewski et al., 2006;Sugden et al., 1995). The mid-latitude glaciation hypothesis considers both the bulk ice emplacement and subsequent recession of these regional ice sheets to be orbitally-driven (Head et al., 2003;Head et al., 2005;. Geomorphological evidence for at least one glacial recession since a maximum advance is abundant (e.g., Brough et al., 2016;Dickson et al., 2008;Levy et al., 2007;Morgan et al., 2009). VFFs within the fretted terrain along Mars's northern dichotomy boundary are thought to have eroded this boundary southwards over the last 2 Gyr , leading to the characteristic formation of interconnected valleys and fragmenting of the terrain into isolated mesas further north. At the local scale, approximately a third of GLFs appear recessed or depressed, as shown by distant forefield arcuate moraines, surface lowering relative to putative lateral moraines, and staggered sequences of latero-terminal moraine ridges (Brough et al., 2016). Tropical piedmont-style LDAs also appear to have retreated by up to 10 km from their maximum extent (Shean et al., 2005), and highstands of palaeo ice-surface elevation up to 0.8 km above current LDA/LVF surfaces have been inferred (Dickson et al., 2008). Further, steep scarps observed at the contact between LVF surfaces and abutting talus deposits indicate substantial LVF surface lowering since the bulk emplacement of that talus (Levy et al., 2007).
Age estimates for LDAs/LVFs from crater size-frequency distributions span 10 Ma to 1.2 Ga (Berman et al., 2015), consistent with the broad decline of regional-scale ice sheets and ice caps over the past several hundred million years Dickson et al., 2008;Head et al., 2010;Levy et al., 2007;Morgan et al., 2009). In addition, it is hypothesised that the switch from high (∼35 • ) to low (∼25 • ) planetary obliquity at 4-6 Ma, the only major change in obliquity during the last 20 Ma, caused the bulk accumulation of the polar ice caps which are approximately four million years old (Laskar et al., 2002;Levrard et al., 2007;Smith et al., 2016). Some authors (e.g., Brough et al., 2016) posit that this obliquity shift also marks the time of the last martian glacial maximum (LMGM), or the last time when mid-latitude glaciation was active on Mars.
That smaller and steeper glaciers on Earth are generally more sensitive to climate change has prompted suggestions (Kargel, 2004) that GLFs may be the most responsive VFFs on Mars. Multiple phases of glacial activity can produce a superposition of some VFFs, formed during one glacial phase, on to other VFFs, formed during an earlier glacial phase. Indeed, several individual examples of VFFs (typically GLFs) flowing onto other VFFs have been recorded (Dickson et al., 2008;Head et al., 2005;Levy et al., 2007;Shean et al., 2005). Figures 1 and 2 show examples of such superposition. The GLFs involved in these superpositions debouch from alcoves and are demarcated by a sharply-defined frontal lobate terminus. Here we term this subset of the GLF population "superposed glacier like forms" (SGLFs). In contrast to the anastomising regional VFF deposits, SGLFs terminate abruptly at the superposition boundary at their lower ends, rather than coalescing with the valley VFF deposits below, and they tend to extend short distances ≤5 km from their source (Figures 1 and 2). Where SGLFs intrude far onto pre-existing VFFs, their flow-parallel surface lineations may be orientated obliquely to those on the latter, with the boundary between the two distinguished by a sharp contact (Dickson et al., 2008;Levy et al., 2007) (Figures 1 and 2). This stratigraphic relationship, indicating that SGLFs are younger than their underlying VFFs, has been invoked as evidence for two separate glacial phases: an early phase emplacing ice sheets/caps, and a later, alpine glaciation whose relics are the SGLFs Brough et al., 2015;Dickson et al., 2008;Morgan et al., 2009;Levy et al., 2007). Further, the lack of deflection and apparent freshness of the SGLF contact indicates that the bulk of SGLF development occurred separately from the emplacement of underlying VFF and relatively recently (Dickson et al., 2008;Levy et al., 2007). Based on the principle of geological superposition, SGLFs may therefore provide a means of deducing the chronological relationships between the distinct scales and phases of glacial advance and recession (Dickson et al., 2008;Levy et al., 2007).
This hypothesis of polyphase glaciation-that Mars's mid-latitudes experienced multiple glacial stages-provides the setting for the investigation presented herein. Multiple phases of glacial activity are plausible given Mars's chaotic obliquity history, and one may generally expect variations in the strength of the orbital cycles to have driven successive glacial phases of different intensities and spatial distributions over time. In this paper, we undertake the first systematic global study of SGLFs in an attempt to elucidate these successions. To unravel the temporal patterns of the mid-latitude glaciation, we derive surface ages with crater counting for SGLFs and the VFFs onto which they flow ("underlying VFFs" hereafter), which date when each of the respective ancestral ice masses stopped flowing actively. It is envisaged that each SGLF is the modern form of an earlier glacier, whose transition to dormancy happened when (or shortly after) climatic changes caused net snow/ice accumulation on it to cease; thinning of the ice under negative mass balance then diminished its flow velocity. Because the transition allowed impact craters as well as supraglacial sediments to build up on the SGLF's surface, whereas active ice flow and snow/ice accumulation prior to the transition caused rigorous resurfacing that continually reset/cleaned the surface, the age determined from today's cratering record reflects the formation age of the SGLF which, in turn, constrains the minimum age of its active ancestral glacier. Note that there are also many locations where laterally-confined viscous flow material (in alcoves or tributary valleys) coalesces with LDAs/LVFs without apparent superposition; not being SGLFs, they are not included in the study.

Mapping
We mapped SGLFs within 25-60 • latitudes in both hemispheres ( Figure 3) on a mosaic of 3,694 georeferenced Context Camera (CTX) images (Malin et al., 2007). We identified SGLFs using (i) the defined characteristics of glacier-like forms (Hubbard et al., 2014;Souness et al., 2012; e.g., alcove heads; enclosing higher topography; presence of longitudinal foliations) and (ii) evidence of superposition onto another VFF (sharp frontal contact delineated by the edge/escarpment of a raised lobe or arcuate moraine-like ridge; surface lineations mismatched from those on the underlying VFF; e.g., Figure 2).
The enclosing boundaries of SGLFs and their underlying VFFs were mapped in the software ArcMap 10.4. Key morphometric parameters of each SGLF (length, width, area, mean surface slope, mean elevation, relief, and orientation) were also derived in ArcMap 10.4 and are detailed in Table 1 of the online data repository linked to this paper (Hepburn et al., 2019). When calculating area and slope, care was taken to use suitable map projections to minimise the adverse effects of distortion. All areas were calculated in the sinusoidal projection. There is potential uncertainty in identifying each SGLF/VFF boundary in our mapping. The associated error in each area, estimated by assuming a uniform 1-pixel (∼6 m) misidentification on each boundary, is a few percent. The uncertainty in each area due to pixel loss/replication during reprojection is expected to be less than this. Surface slope was calculated using the Mars Express High-Resolution Stereo Camera (HRSC; ∼75-m grid resolution) digital elevation model. However, since HRSC coverage is not global, we supplement our analysis with Mars Orbiter Laser Altimeter (MOLA; ∼460-m grid resolution) digital elevation models. Where both data sets were available, both were used to verify that while they yielded slope values differing numerically due to the difference in spatial resolution, the results are broadly consistent so this does not affect our conclusions (section 3). The mean slope of a SGLF is the areal average of elevation gradients within its boundary. To compute this quantity from each gridded elevation map, we followed the method of Conway et al. (2017) of calculating the local gradients of elevation in the conformal Mercator projection and then multiplying these by the secant of the local latitude to correct their baseline distances, before reprojecting the gradients into the sinusoidal projection for areal averaging.

Crater size-frequency distribution measurement
To date the SGLFs and their underlying VFFs, we measured the abundance of supraglacial impact craters. Using the ArcGIS extension CraterTools (Kneissl et al., 2011), with CTX imagery, the diameter D of impact craters (except any discernible chains of secondary craters) in each area defined by a SGLF or underlying VFF boundary was recorded and compiled into the crater size-frequency distribution (CSFD). CraterTools has a built-in correction for the distortion of area under the map projection, where all counting areas and craters are individually reprojected into a sinusoidal projection defined by the central meridian of the feature in question.
Most SGLFs are small landforms at length scales of several kilometres. Their small size and the limited occurrence of impact craters on their surface-especially craters visible at CTX-resolution-means that we had to aggregate multiple SGLFs to enable their dating. In aggregation, the counts from neighbouring SGLFs (typically those emerging from the same peak, mesa, or plateau) were combined to give usable crater statistics for their total area; we assume geological homogeneity (e.g., Arvidson et al., 1979) for each counting area as it pertains to the same landform or the same landform type. Table 1 lists the identification numbers of the SGLFs in each of our 35 aggregates, and Table 2 of our online repository gives the aggregated counts binned according to crater diameter. Aggregation necessarily means that the same age is assigned to all SGLFs in each aggregate so their age differences cannot be resolved, but it extends the dating to individual SGLFs for which ages cannot be ascertained due to scarcity of craters. Two of our SGLFs are large and were dated without aggregation (we number them "aggregates" 10 and 19; Table 1), but most other individual SGLFs must be aggregated as they have too few craters for their CSFDs to show a clear decay limb for fitting with confidence to obtain ages (the fitting method and this issue will be elaborated in section 2.3).
We recorded all visible craters on SGLFs as far as discernible from CTX (∼6 m per pixel; Malin et al., 2007) imagery down to diameters (D) ∼20 m. In contrast, we limited the CSFD measurement on the underlying VFFs to craters larger than ∼80 m in diameter to ensure our recording effort was feasible on the largest VFFs. Age determination for underlying VFFs using crater counts extending down to lower D (< 80 m) occurred Note. CSFD = crater size-frequency distribution, SGLFs = superposed glacier-like forms, VFFs = viscous flow features.
Aggregates 5, 6, 7, and 8 contain the same SGLFs but have different underlying VFF ages. This is because the SGLFs flow onto a large VFF with four sub-units that are distinct in terms of surface texture. We dated each sub-unit to ensure robust age estimation. Number of craters in columns 6 and 11 pertains to the fitted range only. For the total crater counts on each CSFD, see Table 2 of the online data repository linked to this paper. Non-integer counts occur when one or more craters are located on the domain boundary.
10.1029/2019JE006102 in two cases (aggregates 28 and 29, Table 1) because the maximum crater size (i.e., the upper limit of their D ranges) was relatively low. Aggregation was not needed for the underlying VFFs, due to their large size.
In terms of crater morphology, a previous study that dated VFFs (Berman et al., 2015) subcategorised counted craters by their variable states of degradation (e.g., fresh, degraded, and filled) in order to resolve any timing difference between the cessation of flow of the ancestral ice mass (i.e., stabilisation of its surface) and development of the surface-debris cover on the VFF, and to date the emplacement of further overlying mantling deposits due to the latitude-dependent mantling (LDM) process. Dating with subcategorised counts is thus useful for examining the history of surface processes on a VFF after the ancestral ice mass transitioned to the VFF. However, although subcategorisation of our crater counts is possible, it would not aid in our goal of dating the initial transition, which is what defines the formation age of each SGLF or underlying VFF. As explained in section 1, the timing of this formation constrains the minimum age of the formation of the ancestral glacier or ice cap. We emphasise, accordingly, that the actual formation age of the ancestral ice glacier/ice cap can only ever be bounded, not determined. Therefore we did not subcategorise craters in the present study; all counted craters were used in estimating the crater retention age of the surfaces. Note that, as mentioned in section 1, we follow the interpretation that SGLFs may be superposed by younger LDM deposits which contribute material to their surface mantle, but that SGLFs themselves do not derive wholly from re-mobilisation of such deposits, even though GLFs/SGLFs are relatively thin. Our approach of using all craters accords with standard practice as stated by Michael and Neukum (2010). Issues concerned with the potential impact of resurfacing on crater retention ages and the use of small craters are explored further in section 3. It would be a mistake to base our age determination on subcategorised crater counts, as this would yield an underestimation of all formation ages.
The occurrence of "ring-mold" craters on some VFFs (Kress & Head, 2008) is relevant to our procedure of crater identification and size determination. Ring-mold craters are observed across D ∼ 10 2 -10 3 m, and they exhibit an unusual morphology (e.g., some with multiple rings, a central plateau) as a result of post-impact modification processes on an icy substrate overlain by a debris layer (Kress & Head, 2008). Of the total 873 craters we counted on SGLFs, ∼54 (6.2%) were ring-mold craters. This number is approximate as the identification becomes more difficult for small craters at scales approaching the CTX imagery resolution limit. Owing to the icy substrate, ring-mold craters are interpreted to form at diameters larger than the corresponding bowl-shaped craters on basaltic substrates-given the same impacts (Kress & Head, 2008). This size difference led Hartmann (2007) to propose downscaling the diameter of ring-mold craters by a factor of two when their counts are compiled into CSFDs for dating, but this correction factor is uncertain and does not account for the actual ice rheology (which is unknown for the SGLFs and their underlying VFFs). We therefore neglected such scaling and instead followed Berman et al. (2015) by treating ring-mold craters as filled craters without diameter correction, and included them in our counting. Note that while the size difference may have biased some ages upwards, enhanced resurfacing on an icy substrate may have removed some craters to cause an opposite (not necessarily equal) effect on age. As noted below, we also did not attempt to account for resurfacing by using a correction factor.
Care was exercised in our crater identification to exclude circular features unrelated to impact processes, such as sublimation and collapse pits, which sometimes occur in localised populations on SGLF/VFF surfaces. The larger sublimation and collapse pits (at scales ∼100 m) are straightforward to distinguish from impact craters. There is higher chance of misidentifying such features as impact craters as the resolution limit of CTX imagery is approached. We acknowledge that some secondary craters and non-impact features may have been falsely counted within our CSFD data, an issue that pertains to all CSFD-based dating on Mars.

Absolute Surface Age Determination
Next, the absolute ages of each pair of SGLF aggregate and underlying VFF were estimated from their CSFDs.
Since the icy composition of SGLFs and VFFs means that their cratering record may have been modified by some resurfacing after the transition, all absolute ages reported herein are "crater retention ages". In fact, all crater-count based ages of martian surfaces are crater retention ages because all such planetary surfaces are influenced by different kinds of resurfacing processes. We interpret these crater retention ages as formation ages of the SGLFs/underlying VFFs because the surfaces were strongly reset up to the transition (see section 1), and there is limited evidence from our results that equilibrium between subsequent resurfacing and impact cratering on those surfaces has been reached for the crater populations used in our dating. . Impact crater size-frequency distributions of (a) SGLF aggregate 21 and its underlying VFF and of (b) all SGLF surfaces and all underlying VFF surfaces. Data points are derived from crater counts on CTX imagery (Table 2 of the online data repository linked to this paper) and grey curves outline the isochron system. All distributions exhibit roll-off at small crater diameter. In (a), the model ages are determined by isochron fitting (red and black lines) to the filled data points; the corresponding crater-diameter ranges are listed in Table 1. The kink on the underlying VFF data trend above the highest filled black point evidences partial resurfacing. In (b), the surfaces of all underlying VFFs combined show systematically higher crater frequency than for SGLFs on the decay flanks of the distributions. No fits were made to these flanks because their slopes deviate markedly from the isochron slopes (discussed in section 3.3). Both plots are outputs from CraterStats 2.0. SGLF = superposed glacier-like form, VFF = viscous flow feature.
Dating of each cratering record was performed on the reverse cumulative CSFD plot (Figure 4a and supporting information, Figures S1-S35), with the crater diameter (D) axis binned logarithmically according to standard practice (Hartmann, 2005;Michael & Neukum, 2010). The isochron system is constructed from the lunar-specific chronology and production functions by accounting for proximity of the planetary target to meteoric sources and its gravity and collisional cross-section. We used the Hartmann (2005) production and chronology functions for Mars (Hartmann, 2005(Hartmann, , 2007 because they contain improved treatment of these complicating factors for Mars and because they amalgamate crater distribution data sets from both the Hartmann and Neukum traditions. The use of alternative functions, which is not undertaken here, would shift our absolute ages systematically, but not alter their overall pattern. Age was determined by fitting an isochron to each distribution's decay slope (Figure 4a and Figures S1-S35) following established procedures (e.g., Arvidson et al., 1979;Michael & Neukum, 2010). Fitting of each CSFD was done in the software CraterStats 2.0 (https://bit.ly/2U1lcJU; see Michael and Neukum, 2010), which yielded the age and corresponding estimates of uncertainty. CraterStats 2.0 assigns error margins when fitting a CSFD to reflect the statistical uncertainty associated with both the CSFD measurement and the counting area. As detailed by Michael and Neukum (2010), the error bars shown on each cumulative CSFD plot (grey vertical lines in Figures 4a and S1-S35) are given by ±[N cum (D)] 1/2 ∕A, where N cum is the cumulative number of craters of diameter D, and A is the counting area. The "±error" reported with each age estimate in Table 1 and Figures S1-S35 is based on the Poisson distribution model of random cratering impacts (Michael & Neukum, 2010), which describes how the observed crater counts on a surface deviate statistically from the number of crater counts expected from the true age of the surface, for given production and chronology functions. The probability density distribution of the true age constructed from this model yields the age error; the counting area is relevant to this statistical model via the impact rate. The error is calculated automatically in CraterStats 2.0 during fitting. Small counting areas in our study produce relatively large error bars (Figure 5a), which are controlled by our aggregation strategy (section 2.2).
Our CSFDs typically show a roll-off at small diameters (Figures 4 and S1-S35) due to undercounting near the image-resolution limits and removal of small craters by resurfacing processes (e.g., sublimation, mass movement, substrate deformation, burial by debris from hillslopes, and mantling by sediments accreting on the surface of ablating ice). We expect most instances of crater misidentification due to image resolution (section 2.2) lie within the roll-off of the CSFDs; in such cases, the counts do not fall within the fitting range in our dating strategy, so those misidentifications do not influence our age estimates. In other cases, specifically for SGLF aggregates 4, 12-15, 18, 19, 22, 24, 28-30, 32, 33 (Table 1), the entire fitting range lies at small diameter (<80 m) and may be more susceptible to this issue. However, since the ages determined for these SGLFs are all <20 Ma, and misidentification could only have increased crater counts to bias an age upward, our conclusion (below) that these SGLFs belong to recession phase R1 is robust.
In addition to roll-off, some CSFDs show a double-kink reflecting partial resurfacing (Michael & Neukum, 2010) as illustrated by the underlying VFF data in Figure 4a. This double-kink reflects the depletion of small craters by discrete episodic resurfacing events occurring in some narrow time-increment, followed by crater accumulation post-event. We associate the kink/roll-off "shoulder" with the largest crater diameter affected by resurfacing, and therefore fit the isochron to the larger craters plotted to the right of the shoulder, minimising the likelihood that the fitted parts of the distributions have been modified by resurfacing (Figures 4a and S1-S35). The fitted diameter range on the limb was chosen to span the data points that plot parallel to the slope of the isochron system (Berman et al., 2015). In most CSFDs, this fit included all other counted craters ≥D min (the minimum diameter of the fitted range). These criteria ensured that we are dating the end of the last major resurfacing, which marks the time of transition from the ancestral glacier to the SGLF/underlying VFF, rather than subsequent resurfacing episodes that erased craters of smaller diameters. The potential influence of slope-driven resurfacing on the counts on the fitted limb for steep SGLFs is discussed in section 3.3. Table 1 lists the fitted diameter ranges in the CSFD plots for our SGLF aggregates and underlying VFFs, which are also illustrated in Figures 4a and S1-S35.
As explained in section 2.2, we date SGLF aggregates because ages for the majority of individual SGLFs cannot be found due to their small size. The issue lies in the difficulty of deciding where to fit the isochron on the CSFDs of individual SGLFs (these are made up by too few craters and do not show well-defined decay limbs), not so much in achieving manageable errors on the model age after fitting. Dating is thus used at the minimum area scale where the CSFDs are well behaved to yield ages. Consequently, we do not have individual SGLF ages against which the aggregate age estimates can be compared (and potentially validated), and our assumption that SGLFs in an aggregate have similar ages (i.e., the homogeneity assumption; section 2.2) may be open to challenge. We raise this as a limitation of the study, while recognising that there are currently no independent methods that would give individual SGLF ages.
Like other studies of Mars that employ crater chronometry, we recognise that dating with CSFDs is subject to uncertainty from diverse sources. There may be other potential systematic biases/errors in our ages besides those already discussed, stemming from assumptions in the cratering model (impact flux, chronology, and production functions) and assumptions about target material properties. However, crater chronometry is the only means currently available to date the SGLFs/VFFs, and we refer readers to Berman et al. (2015), Hartmann (2005Hartmann ( , 2007, Hartmann and Daubar (2017), Michael and Neukum (2010), Warner et al. (2015), and Williams et al. (2018) for extensive discussions of issues. As in some investigations (e.g., Head et al., 2005), no attempt is made here to correct for potential biases. Crater chronology models should therefore be thought of as best effort calibration to absolute ages, and we hereafter refer to these as "model ages". The use of the Hartmann (2005) function (or indeed any other function) is meant to capture the main characteristics of cratering history on Mars, even if calibration to absolute time is subject to considerable uncertainty. In the specific context of dating SGLFs and VFFs, application of a (extremely poorly constrained) correction for the icy nature of our landforms would probably result in all model ages being revised upward (Landis et al., 2016), but since their near-surface material compositions should be similar (debris mantle on the order of 10 1 m thick, underlain by ice; Baker & Carter, 2019), it is unlikely that such correction will change the relative pattern of the modelled ages to invalidate our qualitative findings about the temporal sequence of mid-latitude alpine glaciation phases. In this sense, although our aim is to report model ages, we are at least as concerned with establishing the relative age contrasts between SGLFs and underlying VFFs and testing whether their CSFDs agree with their relative stratigraphy.
Finally, the use of Mars Reconnaissance Orbiter's High Resolution Imaging Science Experiment (HiRISE: McEwen et al. (2007)) imagery at ∼0.25 m per pixel resolution could give improved statistics for small craters than CTX imagery. If HiRISE imagery covers all of our SGLFs completely, it would enable us to compile CSFDs extending to lower values of D. However, only 31 SGLFs were covered by HiRISE imagery, few of which were close enough to other SGLFs for aggregation. In most cases, the HiRISE imagery provided only partial areal coverage of individual SGLFs, which was further compromised by image shadowing. This sub-sampling of SGLF surfaces limits the number of larger craters recorded on any given feature. Even in the few instances where SGLFs were entirely covered by a HiRISE image, crater counts were at the sub-aggregate level and could not be meaningfully compared to the CSFDs derived from CTX imagery, which included the combined area of several aggregated SGLFs. Therefore, no HiRISE-derived crater counts were used in this study.

SGLF Morphometry and Global Distribution
Using CTX imagery, we identified 320 SGLFs, 251 in the northern hemisphere, and 69 in the southern hemisphere ( Figure 3). The resulting inventory of 320 SGLFs (Table 1 of  Morphometrically, SGLFs are similar to GLFs and statistically indistinguishable from them with respect to length, width, and area data (supporting information, Table S1). Our observations support those of Levy et al. (2007) and Dickson et al. (2008) that SGLFs all share a sharp frontal contact demarcating their boundary with their underlying VFFs. However, morphological differences do exist between SGLFs, particularly in terms of the hemisphere. Northern hemisphere SGLFs generally debouch from an alcove and therefore have clearly 10.1029/2019JE006102 marked boundaries with a sharp terminus (e.g., Figures 1, 2a, and 2c); whereas in the southern hemisphere, up-glacier of the clearly demarcated terminus, the boundaries of SGLFs become more difficult to discern and map (e.g., Figure 2d). SGLFs appear to share with GLFs a similar minimum elevation threshold of approximately 3,000 m below datum (Souness et al., 2012).
SGLFs also share similarity with GLFs in terms of geographical distribution. They occupy the same 30-50 • latitude bands in each hemisphere (Figure 3), where all other forms of VFF are found Levy et al., 2007Levy et al., , 2014. However, note that our mapping window extended to 25-60 • in each hemisphere, similar to other VFF studies (e.g., Levy et al., 2007Levy et al., , 2014Milliken et al., 2003;Souness et al., 2012). SGLFs are clustered in eight regions globally, with the highest densities in the contiguous regions of Deuteronilus Mensae, Protonilus Mensae, and Nilli Fossae (Figure 3), which delineates the most prominent section of the northern dichotomy boundary. High densities of VFFs have also been described amongst the fretted terrain within this region (e.g., Levy et al., 2014;Morgan et al., 2009). Outside the fretted terrain of the dichotomy boundary, SGLFs are also found in Utopia Planitia and Phelgra Montes. In the southern hemisphere, examples of fretted terrain are limited, and SGLFs are concentrated along the rims of the Hellas and Argyre impact basins.

Model Ages of SGLF Aggregates and Underlying VFFs
Crater size-frequency distributions were calculated for 162 SGLFs grouped into 35 aggregates, yielding 35 corresponding model ages. One hundred and fifty eight small SGLFs could not be dated, as their distance from others precluded aggregation. We also derived crater size-frequency distributions for the 35 underlying VFFs (without aggregation) to date their surfaces. Table 1 lists all of the age results. As noted above, our crater counts are given in the online data repository. Figure 4a and Figures S1-S35 document the complete set of 35 crater size-frequency distributions, and their isochron fits in the age determination. Typically, the fitted crater-diameter ranges are one to several hundred metres for underlying VFFs and a few decametres to a hundred metres for SGLF aggregates (Table 1). The potential limitation of small craters in dating the latter population is discussed in section 3.3.
The SGLFs have a model age range of 2-96 Ma (Figures 5a and 5b) and mean model age of 21.9 (±23.8) Ma. Most SGLFs are much younger than their underlying VFFs, whose model ages span 3.6-544 Ma with a mean of 121.5 (±121.9) Ma (Figure 5a). The model age deficit of SGLFs confirm their superposition on the underlying VFFs (despite the current lack of radar evidence capable of resolving the internal structure at the contact) and means that the SGLF's ancestral glacier existed for some time even as the underlying VFF ice cap was transitioning to a VFF. However, of the 35 SGLF aggregates dated, four had a model age exceeding that of their underlying VFF. These four anomalous results of SGLF aggregates dating older than their underlying VFFs contradict their observed superposition relationship, but are reconcilable through model age uncertainty indicated by the error bars in Figure 5a; also, geomorphological examination of the SGLFs and VFFs concerned reveals nothing unique about them compared to other SGLFs/VFFs, nor any shared peculiarities that might explain the contradiction (see supporting information, Figure S36). We also undertook an analysis amalgamating the CSFDs of all SGLFs and all VFFs to calculate two global CSFDs for each of these surface types (Figure 4b). Comparison of these two global CSFDs shows systematically different abundances on their decay slopes (Figure 4b), confirming that SGLF surfaces are collectively younger than VFF surfaces. Although each type of surfaces does not have a unique model age, attempts at fitting the decays in Figure 4b generally yield ages of 20-50 Ma for the SGLF surfaces and 100-250 Ma for the underlying VFF surfaces, which are consistent with the mean model ages reported above.
The SGLF model ages cluster at 2-20 Ma and 45-65 Ma, indicating two distinct glacier recession phases "R1" and "R2", respectively (Figure 5b). These imply pre-exisiting alpine glaciations, "A1" and "A2", respectively, stretching back into the periods preceding R1 and R2 when no SGLFs formed. Formal clustering analysis of the ages with the k-means method (e.g., Everitt et al., 2011) supports our identification of two clusters (supporting information, Figure S37a), and additional Monte-Carlo simulations that account for the uncertainty on each age during the k-means optimisation confirms our cluster assignment for each age (Figures S37b  and 37c). To reiterate, our interpretation takes a SGLF's formation age as the minimum age of its ancestral glacier, assuming that the SGLF formed when the glacier stagnated, thinned and developed a debris mantle. "Recession" in our inferred chronology refers to such mass loss, not necessarily involving length reduction. The ancestral glacier must have formed and been existing for some time prior to SGLF formation. Accordingly, in Figure 5b, we indicate each glaciation (A1 or A2) at older times than the corresponding recession.

10.1029/2019JE006102
We do not rule out new glaciers forming during each recession phase. In contrast, VFFs formed continuously over the past ∼500 Myr (Figure 5a, lower plot), consistent with overall ice-sheet decline; resolving the corresponding history would require dating all LDAs and LVFs, not only those underlying SGLFs. These results reveal a polyphase history of overlapping alpine and ice-sheet style glaciations and specifically two growth-recession cycles of alpine glaciation (Figure 5b). The population at a given time may include active glaciers and ice sheets and relict SGLFs and LDAs/LVFs (see supporting information, Movie S1). It has been hypothesised that the original ice sheets had submerged high peaks, mesas, and plateaus before thinning under long-term ablation by sublimation (Dickson et al., 2008;Fastook et al., 2014;Head et al., 2010). In this picture, our oldest SGLF aggregate ages of ∼65 Ma constrain the latest time when total submergence ended to expose hillslopes for alpine glaciation (the age of 96 Ma for aggregate 10 is one of the four anomalous results described above, so it may be overestimated and is not used here). This inference can be made even if other SGLFs had formed at earlier times but their record has been obliterated.

Robustness of Model Ages
Caution has been expressed regarding chronometry based on populations of small diameter craters (decametres or less), as are often encountered when dating young planetary surfaces. While some studies are critical of its application (Warner et al., 2015;Williams et al., 2018), others are supportive (Hartmann & Daubar, 2017). Our dating of many SGLF aggregates relies on low crater numbers in the decametre range. Besides the care needed to exclude obvious chains of secondary craters (scattered "background" secondaries are included in the production function) and non-impact features from crater counts at small scales (section 2.2), another concern is that resurfacing rates on SGLFs and VFFs may strongly influence their model ages.
Here, we examine whether such process could be responsible for the age difference presented above.
To explore this issue, it is necessary to consider the distinction between discrete episodic resurfacing (which causes kinks in CSFD plots; section 2.3) and quasi-continuous resurfacing. The latter refers to the planet-wide background resurfacing across Mars that acts constantly to degrade and remove craters of all sizes (typically, smaller craters at much higher rates). In contrast, discrete resurfacing events reduce the small crater population only up to a certain size, dependent on the magnitude of the event (Michael & Neukum, 2010).
As discussed in section 2.3, discrete resurfacing is already accounted for by our CSFD fitting strategy. The point of contention here is therefore whether the model ages of the SGLF aggregates significantly underestimate their (true) formation ages as a result of quasi-continuous resurfacing. If so, the model age deficits illustrated by Figure 5a may not quantify the formation-time difference between the SGLFs and their underlying VFFs (instead they may reflect differences in the rates of quasi-continuous resurfacing).
Notably, the debris covered nature and underlying icy substrate of our target surfaces could cause enhanced quasi-continuous resurfacing that decreases the model ages from the formation ages. Although the underlying VFF surfaces are gently sloping (≈10 • or less), the surface slopes of steeper SGLFs (Figure 5c and supporting information, Figure S38) approach the dry-granular angle of repose (≈24-42 • , Kleinhans et al., 2011). This may facilitate surface debris movement on SGLFs as well as internal ice deformation, which could serve to remove larger craters (i.e., craters whose diameters lie in the fitted range, right of the shoulder on each CSFD) as well as small craters.
We consider that quasi-continuous resurfacing has a limited impact on our SGLF dating and does not upset our key conclusions above (section 3.2) for the following four reasons: 1. If the influence of quasi-continuous resurfacing extending into the fitted parts of the crater distributions is widespread on the steeper SGLFs, all decays of the fitted SGLF data would display a log-log slope less than that of the isochrons. However, we do not see evidence of this tendency from Figures S1-S35; most decays are closely aligned with the isochrons. As discussed in section 2.3, we selected our fitting range D min to D max (where D min and D max , respectively, denote the minimum and maximum crater diameters of the range) on each CSFD based on the portion that matches the isochron slope. For those CSFDs without a clear kink, the shoulder between the (visible) roll-off and that portion/limb approximately demarcates the D min value. An inherent assumption behind this choice in our fitting procedure is that crater counts at D > D min are negligibly affected by quasi-continuous resurfacing and so are usable for dating, whereas counts at D < D min are affected and were excluded from our analysis. This assumption is supported by an independent line of evidence from our CSFD data. A general expectation is that quasi-continuous resurfacing is effective (rapid) at small length scales. More precisely, Figure 6. Scatter plots of model ages versus D min (a) and D max (b) for our collection of SGLFs (red) and VFFs (black). D min and D max , respectively, denote the minimum and maximum crater diameters of the isochron fitting range used on the crater size-frequency distribution of the SGLF/VFF to determine its absolute surface age. In each panel, straight lines show best-fit regression lines, for which the corresponding equations, r 2 , and p-values are given. As discussed in section 3.3, the positive correlations found here-notably, that D min is a strong predictor of both SGLF and VFF model age in (a)-validates the assumptions behind our isochron-fitting procedure, meaning that our dating of SGLF/VFF formation is not hampered by the undesirable effects of quasi-continuous resurfacing on the landforms.
resurfacing acts on all length scales, but smaller craters resurface more rapidly than larger ones (Michael & Neukum, 2010). For a given surface of a known true age, one should observe a well-defined crater-size threshold separating those craters significantly affected by resurfacing from those insignificantly affected by resurfacing-even if the transition between these crater types is expected to take place over a finite range of D (i.e., not abruptly). Indeed, the roll-offs on many of our CSFDs curve smoothly to join the straight-sloped limbs (Figures 4a and S1-S35). Our choice of D min in each fit essentially demarcates this transition. Now, the older the surface, the more time it has been exposed to resurfacing. Increased exposure time eliminates progressively larger craters, lowering their frequency. Therefore D min should increase with the true age of the surface. Figure 6a plots (model) age versus D min in our CSFD dating fits for both the SGLF aggregates and the underlying VFFs. These data exhibit trends that confirm this expectation, with best-fit regression lines showing strong correlation between the variables with an r 2 of 0.78 and 0.57 for SGLFs and underlying VFFs, respectively. In fact, the data trends (and trendlines) for the SGLFs and VFFs nearly overlap, further confirming that the cumulative action of resurfacing over time affect crater populations in the same manner on both landforms. This analysis means that (i) our visual determination of D min values in the CSFD fitting procedure, and consequently, (ii) our assumption of the fitted limbs not having suffered from resurfacing (thus closely reflecting true surface ages), are sound (i.e., there is self-consistency between the outcomes of the method of fitting/age determination and its underlying assumptions). Finally, we note that the model ages also correlate (positively) moderately strongly with D max (Figure 6b), though not as strongly as with D min . The correlation with D max is not surprising because the likelihood of observing larger craters in a given area should increase with its surface age (as well as with the size of the area). 2. Despite the steepness of some SGLFs, the rate of resurfacing caused by internal deformation within them should be low, because we expect SGLFs to be thin (much thinner than VFFs generally) and this implies limited ice-flow velocity. For instance, a tentative glaciological model (Karlsson et al., 2015), with its optimised plastic yield stress of 22 kPa and typical observed surface slope angles (15-30 • , Figures 5c and supporting information, S38-S39), yields a ball-park estimate for the thickness of SGLFs of only 11-24 m, which is consistent with their often concave surfaces (Hubbard et al., 2014). Although this estimate does not include a calculation of ice-flow velocity (e.g., using Glen's flow law, Cuffey & Paterson, 2010) because the ice composition/rheology and temperature history are both uncertain, it points towards very subdued ice flow and limited reworking of their impact craters.
3. If the SGLFs were much older than their apparent model ages, then their crater distributions would have reached near-equilibrium, reflecting a balance of cratering and resurfacing rates. Thus, their model ages would be systematically (negatively) correlated with slopes; yet this is not observed (r 2 = 0.01: Figure 5c). This absence of correlation holds for both HRSC and MOLA surface slopes ( Figure S39).
4. If our SGLF model ages have been grossly underestimated due to slope-driven resurfacing (because of debris flow and/or ice flow), then it would not be possible to correct them in a physically-consistent manner. Imagine that the supposition above were true, and a sizeable upward correction is applied to all of the SGLF model ages-a correction large enough to eliminate or drastically reduce the age deficits for those SGLFs in recession R1 (Figure 5a)-to disrupt our interpretation. The SGLF model ages in recession R2 that currently lie close to their underlying VFF model ages (eight points in square symbols near the 1:1 line, Figure 5a) would then exceed the underlying VFF model ages after the correction, but this outcome contradicts the observed morphological superposition of SGLFs on their underlying VFFs. In our picture of SGLFs/LDAs as relicts of ancestral glaciers/ice caps, it is implausible for an underlying VFF to form after its SGLF(s): the underlying VFF lies at lower elevation than its SGLF(s), and it is difficult to imagine how its ancestral ice cap/ice sheet could have a zone of net snow/ice accumulation and avoid transition into a VFF when the ancestral glacier (at higher elevation) transitioned to become a SGLF as a result of negative mass balance. Moreover, those SGLF aggregates belonging to recession R1 are not systematically steeper than the eight aggregates belonging to R2 identified by square symbols (Figures 5c and S39), so an alternative supposition that only the SGLF ages in R1 require an upward correction has no basis.
Taken together, these arguments indicate that slope-driven resurfacing has not hampered our dating of the SGLFs, so our interpretation of the model crater-retention ages derived for them as their formation ages is reasonable. The first, third, and fourth arguments above rest directly on evidence coming from our crater counts, irrespective of the chronology model used. It could still be argued that an age correction should be applied to both SGLFs and underlying VFFs in view of their icy composition, which can promote resurfacing to bias model ages downward and/or cause craters to form at larger diameters than assumed in the crater production model to bias model ages upward (Berman et al., 2015;section 2.3). However, such correction would not change the relative pattern of the ages to disrupt our interpretation of (i) SGLFs being generally substantially younger than their underlying VFFs and (ii) SGLF ages clustering into two discrete ranges. We conclude that even if our model age estimates deviate from true formation ages and one thinks only in terms of contrasts in the number of craters observed on VFFs and SGLFs, our inference of two alpine style glaciations and a polyphase mid-latitude glacial history holds.
We note that the decay flanks of the CSFDs in Figure 4b (which relate the total-thus, improved-crater statistics for all SGLF and all underlying VFF surfaces) are markedly shallower than isochrons. Not only does this shallowing prevent specific fits from being made to the flanks to yield unique ages (section 3.2; caption of Figure 4b), but it may also indicate that the fitted limbs of individual CSFDs used to date our SGLF aggregates and underlying VFFs have been pervasively affected by resurfacing that removed craters at D > D min . Specifically, it could be argued that unless the fitted limbs were compromised by resurfacing, the observed shallowing cannot arise, for example, solely from the large spread of ages of SGLF aggregates and underlying VFF. Consider merging the crater counts of two surfaces with different ages and areas, both unaffected by resurfacing; their CSFDs are Hartmann isochrons following the power laws N 1 = c 1 D − and N 2 = c 2 D − where = 3.1 (the log-log slope shown in our CSFD plots). The aggregated CSFD will be N = cD − , with the log-log slope preserved, and with the coefficient c being an area-weighted combination of c 1 and c 2 . However, this argument overlooks the fact that the distributions aggregated to give Figure 4b, and indeed, the CSFDs collected for any real planetary surface, are not functions extending over an infinite range in D. The CSFD data for each SGLF aggregate or underlying VFF have minimum and maximum crater diameters that vary substantially across aggregates and underlying VFFs. There is no reason why aggregating such data should yield the power law deduced above; consider adding two functions N 1 and N 2 , defined to be ∝ D − within two arbitrarily different finite ranges in D and defined to be zero outside those ranges. Consequently, we do not regard the shallowed CSFD decay flanks in Figure 4b as posing a problem to our dating results; they are due to the specific shapes and combinations of the crater distributions on SGLFs and underlying VFFs. It is interesting to ponder whether there are generic causes behind the shallowing. While we do not undertake further analysis into this here, we note that the CSFD of all craters ≥1 km on Mars is found to exhibit a decay slope ( ≈ 1.1) much shallower than the isochron system (Robbins & Hynek, 2012; see their Figure 8a). This global compilation necessarily includes surfaces with different ages and resurfacing histories, as for the CSFDs in Figure 4b.  (Dickson et al., 2018). (b) Conceptual sequence of SGLF formation through two alpine glaciation-recession cycles; note that these schematics are not to scale. In stages 1 to 3, a glacier formed on top of a viscous flow feature in alpine phase A2 transitions into a SGLF during recession phase R2. Protected by its surface debris mantle, this SGLF is later overlain by another glacier in alpine phase A1 (stage 4). Whether a SGLF forms from this new glacier in recession phase R1 (stage 6) or not (stage 5) depends on its level of protection by surface debris and the thickness/extent of the A1 glacier, as determined by local geological and climatic factors. These contrasting outcomes can lead to the coexistence of R1 and R2 SGLFs within the same geographical area, as illustrated in (a).

Broader Implications
The accelerating formation of SGLFs through recession phase R1 (Figure 5b) documents the demise of the last alpine glaciation (A1) on Mars. Its culmination in 2-4 Ma, shortly after the obliquity transition at ∼5 Ma, coincided with the accumulation of the bulk of the north polar ice cap ∼0-4 Ma (Laskar et al., 2002;Levrard et al., 2007;Smith et al., 2016), providing evidence supporting poleward water transfer during mid-latitude deglaciation due to reduced obliquity (Head et al., 2003. However, even if the ancestral glaciers of all SGLFs were as thick as 1 km, they would still contribute just 0.4% of the present-day north polar ice cap volume (∼800,000 km 3 ; Smith et al., 2016), so a significant water contribution to the polar cap must come from other VFFs or reservoirs such as buried ice deposits or LDM at mid-low latitudes. Notwithstanding dating uncertainty, we gather from Figure 5b that alpine phase A1 apparently reached its maximum between ∼20 and ∼10 Ma. Recent investigations of glacial landforms on the interior walls of martian impact craters (Conway et al., 2018;Hartmann et al., 2014) proposed that glaciers eroded some of the crater walls in the mid and low latitudes in a distinct episode dating to ∼5-10 Ma. This timing is compatible with our finding that numerous mid-latitude alpine glaciers (of phase A1) existed at ∼5 Ma and had not transitioned to SGLFs until later (Figure 5b).
We do not know how many glaciers survived the penultimate recession phase R2 into A1 or whether new glaciers in A1 overprinted some SGLFs from R2. The latter scenario is likely given the overlapping spatial occurrence of SGLFs left as relics of A1 and A2 (Figure 3). In this connection, we find several regions where R1 and R2 SGLFs occur in close proximity-within 60 km-of each other (e.g., Figure 7a). Since these SGLF populations often do not differ distinctly in their orientation preference (e.g., Figure 7a) or other apparent geographical factors, it is difficult to imagine how the A1 phase could have selectively glacierised some areas but not others over such short distances, to give the mixture of R1 and R2 SGLFs observed today. These considerations lead us to infer some R2 SGLFs could have been sufficiently preserved to survive through A1 without being overprinted/buried by R1 SGLFs, owing to preservation factors. Figure 7b illustrates this tentative concept, although we note that our ability to test this hypothesis is limited. In our conceptual model, during recession phase R1, a glacier formed in A1 on top of an R2 SGLF (stage 4 in the figure) would transition to a new SGLF to bury the R2 SGLF if adequate debris/sediment supply causes a protective mantle to develop on its surface (stage 6). However, if the glacier lacked sufficient debris supply or was thin or small in extent, it could decay without leaving an R1 SGLF (stage 5); then the R2 SGLF is observed today. Accordingly, a mixture of neighbouring R1 and R2 SGLFs can arise from variations in geological factors (e.g., rock layering, headwall geometry) at mesa or alcove scale, as well as from local glacio-climatic factors that limited or prevented renewed glacial accumulation on R2 SGLFs during A1. By extension, it is also possible for some glaciers in the earlier A2 phase (stage 2, Figure 7b) not to leave behind R2 SGLFs if they had too little debris supply. Thus, a further potential outcome (not illustrated in Figure 7b) is the absence of SGLFs in some alcoves or tributary valleys, with viscous flow material integrating smoothly with the VFF at lower elevation (without superposition), whereas SGLFs occupy nearby alcoves/valleys. Other hypotheses may invoke flow reactivation on select R2 SGLFs without the need for additional mass input (e.g., as a result of localised subglacial warming from below) which reset their cratering record in phase A1 and caused them to adopt an R1 age.
SGLF locations provide additional insights into the two alpine phases identified above. The older SGLFs from A2-R2 are limited to Deuteronilus-Protonilus Mensae and Tempe Terra in the northern hemisphere, whereas those from A1-R1 straddle both hemispheres with broadly-synchronous age distributions. The latter SGLFs are also three times more abundant (Figures 3 and 5 and Movie S1). The orientation of the hemispheric SGLF population from A2-R2 is biased equatorward, but for A1-R1, it is poleward, with a secondary south-pointing lobe for the northern SGLFs (Figure 5d), as found for the GLF population (Souness et al., 2012). It is tempting to interpret these strong orientation differences in terms of contrasting climates during the high-obliquity periods of these glaciations, in particular by considering glacier accumulation and shifts in prevailing wind patterns. However, we caution that A2's antiquity may mean that its SGLF imprint is incomplete, and a sound comparison must account for contrasting preservation (see last paragraph), which is currently poorly understood. Likewise, since Mars's obliquity history before 20 Ma cannot be reconstructed (Laskar et al., 2004), we cannot verify whether A2 was the less intense of the two alpine glaciations, as the SGLF abundance would suggest. Consequently, we think that the more recent A1 offers a more faithful spatial record for modelling studies of alpine glaciation mechanisms.

Conclusions
We have identified and examined 320 SGLFs throughout the mid-latitudes of Mars. These SGLFs-debouching from alcoves and crater walls onto underlying VFFs-have previously been noted in a few locations along the northern dichotomy boundary. Their apparent superposition on underlying VFFs indicates that SGLFs formed more recently in Mars's history than their underlying VFFs, but the formation age differences between these forms have hitherto not been quantified. Using crater-dating, we confirm the relative age difference in terms of crater density and assign absolute model ages to date the formation of 162 SGLFs (in 35 aggregates) and the corresponding underlying VFFs for which crater populations were recorded. Whereas the underlying-VFF ages appear to be spread diffusively over the last 300 Ma, the SGLF ages are bimodal, clustering at ∼2-20 Ma and ∼45-65 Ma. Despite potential limitations of the dating method, the qualitative aspects of these results are considered to be robust. Based on our interpretation of SGLFs as relict features formed from the decline and stagnation of alpine-style glaciers, we interpret the two clusters as evidence of two growth-recession cycles in the martian mid-latitude glaciation record. These growth-recession cycles evidence climatic changes on Mars on timescales of 10 1 to 10 2 Myr, before its period of reconstructible astronomical forcing and before the earliest times recorded in its north polar ice cap. The SGLF age results suggest ∼65 Ma as the latest possible time at which widespread ice sheet submergence ended in the mid-latitudes to expose hillslopes for alpine glaciation. Our first glimpse of the last stages of mid-latitude glaciation highlights a spatiotemporal complexity that is probably common in Mars's cryospheric evolution-similar to Earth's. Future work should extend the tentative mid-latitude glacial chronology presented here by dating the many remaining GLFs as well as LDAs and LVFs to form an integrated reconstruction. performed all mapping, crater-counting measurements, and dating work. A.J.H. and F.S.L.N. designed the study and analysed the results and produced the manuscript with input from S.J.L., T.O.H., and B.H. All authors discussed the results and commented on the writing. All data derived from our mapping (the inventory of SGLFs and their morphometry, and ArcGIS shapefiles of SGLF boundaries and all craters discussed) and all data used in our dating procedure (impact crater counts) are given in an online data repository (Hepburn et al., 2019). The crater size-frequency distribution plots for our 35 SGLF aggregates and their underlying VFFs are provided in the supporting information.