Late Amazonian Ice Survival in Kasei Valles, Mars

High‐obliquity excursions on Mars are hypothesized to have redistributed water from the poles to nourish mid‐latitude glaciers. Evidence of this process is provided by different types of viscous flow features (ice‐rich deposits buried beneath sediment mantle) located there today, including lobate debris aprons (LDAs). During high‐obliquity extremes, ice may have persisted even nearer the equator, as indicated by numerous enigmatic depressions bounded on one side by either isolated mesas or scarps, and on the other by a lava unit. These depressions demarcate the past interaction between flowing lava and ghost LDAs (GLDAs), which have long since disappeared. We term these features GLDA depressions, about which little is known besides their spatial extent. This collection of depressions implies tropical ice loss over an area ∼100,000 km2. To constrain their history in Kasei Valles, we derive model ages for GLDA depressions, mesas, and the lava flow from crater counts. We use a 2D model of glacial ice constrained by the topography of GLDA depressions to approximate the surface and volume of former glacial ice deposits. The model reconstructs former ice surfaces along multiple flow lines orientated normal to GLDA depression boundaries. This reconstruction indicates that 1,400–3,500 km3 of ice—similar to that present in Iceland on Earth—existed at ∼1.3 Ga when the lava was emplaced. Dating shows that the GLDAs survived for up to ∼1 billion years following lava emplacement, before its final demise.

. The Kasei Valles region, three ghost lobate debris apron (GLDA) depressions, and a contemporary LDA. (a) The Kasei Valles study region. Background colors show Mars Orbiter Laser Altimeter gridded topography in a equirectangular projection, overlain by the Murray Lab Context Camera (CTX; Malin et al., 2007) mosaic (Dickson et al., 2018). The dashed line indicates the southern extent of the Amazonian and Hesperian volcanic unit (AHv) lava flow (Tanaka et al., 2014) in the Martian atmosphere (Chamberlain & Boynton, 2007;Head et al., 2003Head et al., , 2005Head et al., , 2006aSchorghofer & Aharonson, 2005). Mars's orbital parameters are chaotic (Laskar et al., 2002(Laskar et al., , 2004 and through modulation of the insolation receipt, their fluctuations cause glaciations at times of high average obliquity (≥35°-40°) when volatiles were redistributed from the polar ice sheets to nourish mid-latitude ice (Head et al., 2006b;Laskar et al., 2004). At the low obliquity (∼25°) of the present day, mid-latitude VFFs represent the relics of more extensive ice sheets and ice caps-formed during phases of high obliquity-the bulk of which has since sublimated. Modeling of Mars's orbital parameters (Laskar et al., 2004) indicates that a major shift in obliquity (from ∼35° to 25°) occurred around 4-6 Ma. Several studies have linked this shift (sometimes referred to as the Last Martian Glacial Maximum; Brough et al., 2015) with landform evidence of recent glacial activity in the last several million years, including the bulk accumulation of the Northern Polar Layered Deposits (NPLD; Laskar et al., 2002;Levrard et al., 2007;Smith et al. , 2016) and the formation of alpine glaciers (Hepburn et al., 2020b) that left behind the glacier-like forms (a sub-classification of VFFs) observed today in the mid-latitudes (Souness et al., 2012).
Climatic modeling indicates obliquity values approaching 40° are not sufficiently high for extensive precipitation at the latitudes of Kasei Valles (Madeleine et al., 2014), and, therefore extremely high obliquity is probably needed to explain the GLDA depressions. Mars's obliquity history beyond 20 Ma is chaotic and can only be bounded, not determined (Head, et al., 2006b;Laskar et al., 2004). Statistical density modeling of 20,000 possible solutions for Mars's orbital history indicates Mars's mean obliquity is ∼37.62° over the last 4 Gyr, but that it may have approached values as high as 90° (Laskar et al., 2004). Although the measurement of elliptical craters suggests this may have only happened a handful of times since the Noachian (Holo et al., 2019), at such high obliquity, volatile redistribution may plausibly have extended beyond the mid-latitudes towards the equatorial region. Indeed, several landforms other than GLDA depressions linked to glaciation at low latitudes have been reported, concentrated around the large volcanoes of the Tharsis province Milkovich et al., 2006;Murray et al., 2005;Scanlon et al., 2014;Shean et al., 2005;Sinha et al., 2017) and elsewhere (Gourronc et al., 2014).
On Earth, glacial action often leaves behind extensive stratigraphic records from which glacial history and dynamics can be reconstructed (e.g., Clark, 2011;Hughes et al., 2014). For Mars, we have limited access to such records. Instead, dating of glacial activity-constrained by counting the size-frequency of supraglacial HEPBURN ET AL.  Hauber et al., 2008). In this model, craters accumulate on the GLDA depression during/after the loss of ice. The model age for the GLDA depression surface reflects this accumulation and pre-GLDA crater counts in (a) (i.e., inheritance). In each panel, we are assuming that the mesa scarp boundary is effectively static.
craters-has mostly been confined to single stratigraphic units, such as large VFFs (e.g., Baker et al., 2010;Berman et al., 2015;Dickson et al., 2008;Head et al., 2010;Levy et al., 2007;Morgan et al., 2009) and the NPLD (e.g., Landis et al., 2016). The target strength of ice modifies the size-frequency of the initial craters produced by impacts and their subsequent preservation, relative to other substrates. This is known to be the case on the NPLD (Landis et al., 2016). However more generally, the effect of icy substrates on the initial expression and long-term preservation of supraglacial craters on VFFs remains poorly constrained (Berman et al., 2015;Hepburn et al., 2020b). Seeking to address this uncertainty regarding VFFs, Fassett et al. (2014) mapped large craters (>2 km in diameter) which appeared modified by, and those which formed following, glacial activity. However, due to the relative scarcity of these craters, the insights derived regarding the effect of icy substrates was limited to broad geographical areas .
The GLDA depressions in our study area present a unique opportunity for reconstructing ice masses at a regional scale and unraveling the glacial stratigraphic record of former low-latitude glaciation. Previous studies have dated the lava unit within Kasei Valles to 1.0-1.6 Ga (Hauber et al., 2008). Mid-latitude VFFs are typically hundreds of Myr old (Berman et al., 2015) and the latitude of GLDA depressions is beyond the known extent of contemporary VFFs. It is not known whether GLDAs in Kasei Valles formed contemporaneously with mid-latitude VFFs or represent an earlier generation (Parsons et al., 2011).
The aim of this study is to constrain the chronology and expression of past glaciation in Kasei Valles, and we use our findings to evaluate competing models of glacier formation and recession for the region. We use a 2D model of glacial ice constrained by the topography of GLDA depressions to approximate the surface and volume of former glacial ice deposits. Crater counting is used to date the GLDA depressions, the mesas they encircle and the lava deposits which mark their boundaries. Our reconstruction provides fresh insights into the chronology and longevity of ice in the low latitudes of Mars.

Kasei Valles
Kasei Valles is one of the largest outflow channels on Mars, ∼1,500 km long and ∼300 km across at its widest point. It is continuously traceable from the volcanic Tharsis bulge to Chryse Planitia in the northern plains (Rotto & Tanaka, 1995;Scott & Tanaka, 1986;Tanaka & Chapman, 1992). GLDA depressions are found on a large expansive terrace, previously mapped as Amazonian-aged lava flows (Hauber et al., 2008), to the north of the Kasei Valles channel. This volcanic unit, forming part of the Amazonian and Hesperian volcanic unit (AHv: Tanaka et al., 2014) comprises much of the Tharsis and Elysium rises and consists of stacked lobate flows 10-100's of meters thick and cumulatively several kilometers thick in places (Tanaka et al., 2014). Dating of the AHv unit as a whole suggests it formed between 1 and 3 Ga (Dumke et al., 2010;Hauber et al., 2008;Tanaka et al., 2014). To the north of this volcanic unit is a sharp topographic scarp over 1 km high that separates the Tempe Terra plains from the terrace below ( Figure 1a).

Geomorphology Reconstruction
We mapped GLDA depressions using the Context Camera (CTX; Malin et al., 2007) image mosaic from Caltech's Murray Lab, which is a globally seamless mosaic of georeferenced CTX at 5 m/pixel created by calculating the paths of least contrast among overlapping images (see Dickson et al., 2018). The mosaic was inspected for any obvious image artifacts which may affect mapping, none of which were observed. GLDA depressions were identified according to the definition of Hauber et al. (2008): the presence of a depression, sharply distinguishable from the AHv unit, that either (i) surrounds a mesa or (ii) abuts the base of a scarp (e.g., Figures 1b and 1c).
GLDA depressions were identified within the area of Kasei Valles shown in Figure 1a. The mapped lava deposits extended east and west of this region and we mapped the portion of the AHv lava unit which lies east of the extensive N-S faults (Tanaka et al., 2014). As has been done previously (e.g., Tanaka et al., 2014), we mapped the AHv unit as one single unit (>100,000 km 2 ) given the absence of any clear systematic variations in the surface texture at 5 m/pixel resolution. ArcMap 10.5 was used to digitize the boundary demarcating HEPBURN ET AL. 4 of 18 10.1029/2020JE006531 the transition between the GLDA depressions and the lava deposit, around mesas and at the base of scarps. Scarp-parallel GLDA depressions are only partly demarcated by the lava deposit and, in the absence of any ice-related landforms on the Tempe Terra plains (Figure 1), we use the break of slope at the scarp top to digitize the remaining boundary. The uncertainty of each GLDA depression area is less than a few percent, assuming a uniform 1-pixel (∼5 m) misidentification on each boundary. We also investigated the aspect ratio of GLDAs as indicated by the runout length of each GLDA depression and the relief of the mesa/scarp. To measure the runout length of GLDAs, defined here as the distance between the outer GLDA depression boundary with the lava and inner boundary at the base of the scarp, we calculated the nearest Euclidean distance neighbor between two sets of points spaced 100 m along each boundary. The approximate relief of each mesa or scarp within or abutting GLDA depressions was calculated as the standard deviation of elevation within the bounding geometry of each GLDA depression. Elevation was taken from a blended Mars Orbiter Laser Altimeter (MOLA) and High Resolution Stereo Camera (HRSC) digital elevation model (DEM) created using the Ames Stereo Pipeline function dem_mosaic (Beyer et al., 2018) using the D_Mars datum value of 3,396.19 km. The resolution difference between MOLA (∼460 m/pixel) and HRSC (100 m/ pixel) means that although the cell resolution is 100 m/pixel, this is only effectively true where the higher resolution DEM is available ( Figure S1).

Ice Volume Reconstruction
We use a 2D perfect-plasticity model of ice flow to reconstruct the volume of GLDA ice demarcated by the mapped boundaries of GLDA depressions. Similar models based on the assumption that ice deforms plastically have been extensively used on Earth to reconstruct paleo ice volumes (Ng et al., 2010) and more recently have applied to contemporary LDA deposits on Mars (Fastook et al., 2014;Karlsson et al., 2015;Parsons et al., 2011;Schmidt et al., 2019;Weitz et al., 2018). We adapt the numerical method described by Benn and Hulton (2010), driven by parameters from Karlsson et al. (2015). Perfect-plasticity models assume that ice deforms in response to glaciological driving stress (τ D , determined by gravity, ice thickness, and surface slope) when a threshold yield stress τ y is surpassed. A perfectly plastic flow is assumed to adjust continuously such that τ D = τ y (Nye, 1951(Nye, , 1952; thus, where ρ is glacier ice density (ρ = 918 kg m −3 ), g is gravitational acceleration (g = 3.711 m s −2 on Mars), H is the ice thickness, and h is the ice surface elevation. On each flow line, we orient the x-axis to point horizontally up-glacier, with x = 0 located at the mapped GLDA boundary. To reconstruct the ice surface, we integrate Equation 1 by using the finite-difference method at a horizontal space step Δx of 50 m. The discretized equation is (e.g., Schilling & Hollin, 1981) where i = 1,2,3…n is the position index and H * is the representative ice thickness over each step, with the ice thickness given by H = h−B at all positions. Here, the three choices of taking H * to be H i , (H i + H i+1 )/2 and H i+1 amount to forward, central, and backward difference approximations, respectively. As discussed by Benn and Hulton (2010), the forward scheme breaks down at the start of the integration because H = 0 at the ice margin (i = 1); only the central and backward schemes are appropriate there. We adopt the backward approximation, that is, H * = h i+1 −B i+1 , so that substituting this into Equation 2, rearranging, leads to the quadratic equation Solving Equation 3 for its positive root (the negative root is unphysical) yields the surface elevation on the flow line (Figures 3b-3e), h i+1 , from the known quantities at step i and B i+1 , and its successive application reconstructs the surface profile. Because our Δx is relatively small, our numerical results would change negligibly if the central scheme were used (this was the choice of Benn & Hulton, 2010). On some flow lines, the ice can become very thin as the integration marches towards the scarp wall or flank of a mesa (this happens especially when the chosen yield stress is low), and this can lead to unrealistically steep ice surfaces. We therefore ended the integration if H ≤ 25 m and dh/dx exceeded 0.035 ( Figure 3). This stopping criterion did not activate on most of our flow lines as their computed surfaces extended above and over the mesa topography, and it is not used on the initial steep ice surface near x = 0.
The reconstruction requires τ y to be prescribed. Karlsson et al. (2015) estimated yield stresses ranging from 16 to 80 kPa for several LDAs at the mid-latitudes today, by using Equation 1 in an inverse Monte Carlo scheme to best match their surface and basal topographic data, the latter derived from orbiting radar sounding measurements. In comparison, τ y for glaciers and ice sheets on Earth typically span 50-200 kPa (Cuffey & Paterson, 2010;Li et al., 2012). As the ice-flow dynamics of the GLDAs in our study region may be different from those of the contemporary LDAs, we assumed a large range for τ y in our study, varying it between 16 and 200 kPa ( Figure 3e).
To create a 3D ice surface from the surface profiles reconstructed along flow lines, we treat the output network of points derived above as a pseudo point cloud. The Ames Stereo Pipeline was used to produce a simulated DEM of former ice surfaces (Beyer et al., 2018;Figure 3b-d). We used the point2dem program that takes a point cloud and rasterizes it using a user-defined interpolation scheme.

Crater Size-Frequency Distribution Measurement and Absolute Age Determination
To date the GLDA depressions, mesas, and the lava deposit, we measured the abundance of impact craters using the ArcGIS extension CraterTools (Kneissl et al., 2011) and the 5 m/pixel Murray lab CTX mosaic (Dickson et al., 2018). Craters were recorded within areas of geological homogeneity, defined here as the continuous boundary of the three elements of interest. Not all mesas were dated because some of them are small and/or dominated by steep angular morphometry, which we anticipate would limit the preservation of craters through enhanced resurfacing. On mesas that were dated, we limited counts to flat-topped portions and excluded any steep flanks. CraterTools reprojects each crater-delineated by marking three locations at opposing points of the rim-into a sinusoidal projection. The central meridian of each crater's projection is defined by that crater's centroid coordinates and crater diameter is calculated along this unique meridian. This approach limits projection distortion which can deform craters oblately along their longitudinal axis.
We recorded the diameter, D, of craters visible at a scale determined by the area of the landform queried, excluding all identifiable obvious chains of secondary craters. Care was also exercised in our crater identification to exclude circular features unrelated to impact processes (e.g., sublimation/collapse pits). The counts were compiled into a crater size-frequency distribution (CSFD). We varied the scale at which craters were counted to ensure (i) reliable statistics for GLDA depressions (minimum area ≈10 km 2 ; counting scale 1:10,000) and (ii) a feasible counting effort for the regional lava deposit (area ≈74,000 km 2 ). The minimum D of craters recorded on GLDA depressions was ∼20 m and ∼100 m on the lava deposit. A fixed minimum crater counting diameter is not possible because crater visibility depends on factors besides crater size, including lighting conditions and depth. No systematic distinction was made between craters of different morphologies. In the absence of clear variations in texture within count areas at CTX resolution, we made no attempt to correct for unconstrained variations in target properties which may alter the expression of small craters toward the lower diameter limit of our counts (e.g., Dundas et al., 2010).
Age determination from crater count statistics was performed using the the CraterStats 2.0 software (see Michael and Neukum, 2010) on the "differential" style of CSFD plots. For each CSFD, crater diameter was binned using a √2 strategy. An isochron was fitted to the decay slope ( Figure 4) by selecting a span of diameter bins that plot parallel to the slope of an isochron (Berman et al., 2015) following the established procedure (Michael & Neukum, 2010). Isochrons were created using the Hartmann and Daubar (2017) production and chronology functions for Mars. Our CSFDs show roll-off at small diameters due to undercounting as one approaches the resolution-limit of CTX images (5 m/pixel) and due to removal of small craters by resurfacing processes (e.g., mass movement, deflation of GLDA surfaces and redeposition of surface debris, burial by hillslope derived debris, and mantling by eolian deposits). The D at which this roll-off "shoulder" occurs informs the lower end, or D min , of our fitting range, which is always set to the right of this shoulder to ensure that our fitting excludes the biased data of the roll-off. We anticipate that most misidentified craters (due to image-resolution) are small ones falling within the roll-off and therefore will not influence our age estimates. For GLDA depressions and mesas (the craters of which were measured at 1:10,000 scale), D min typically occurs at D ≈ 75 m. On the lava flow (1:50,000 scale), D min occurs at D ≈ 250 m. The Martian isochron system is constructed from the lunar-specific chronology and production functions by accounting for proximity of the planetary target to meteoric source, its gravity, and its collisional cross-section. We used the Hartmann and Daubar (2017) production and the Hartmann (2005) chronology functions for Mars because they contain improved treatment of these complicating factors by amalgamating crater distribution datasets from both the Hartmann and Neukum traditions. At the diameter of craters studied by us, the Hartmann and Daubar (2017) production function and the earlier Hartmann (2005) function are indistinguishable from each other and are effectively interchangeable. The two functions differ only in their treatment of decameter-scale craters (below our GLDA depression roll-off at diameters <100 m). Dating planetary surfaces using CSFD measurements is subject to uncertainty arising from assumptions in the cratering model (e.g., impactor flux, chronology, and production functions), assumptions about the properties of the target material, and potential scale-dependent control on the age result interpreted from crater populations (see Fassett, 2016 andPalucis et al., 2020). Accordingly, all ages reported herein are referred to as "model ages." The use of the Hartmann and Daubar (2017) function is intended to capture the main characteristics of cratering on Mars. Alternative functions, though not used here, would alter our model absolute ages systematically, but not change their overall pattern.
Differential plots are a common means of presenting CSFD data in planetary science (Hartmann, 2005) with craters binned according to their diameter and plotted independently of other bins. Cumulative plots are another common means of presenting CSFD data whereby craters are binned in a reverse cumulative fashion (Crater Analysis Techniques Working Group, 1979). Cumulative plots are useful because the curve rapidly converges to an isochron from which age can be derived (Michael, 2013). However, resurfacingacting to preferentially remove smaller craters (Fassett, 2016;Michael & Neukum, 2010)-is easier to identify on differential plots (Fassett, 2016). Because bins are plotted independently from one another, smaller bin diameters are not offset (upwards) by the cumulative influence of larger diameter bins and downturns in crater frequency are more obvious (Fassett, 2016). Where resurfacing occurs in a discrete time increment and the maximum erasure of craters is limited by the magnitude of that event, it may be possible to estimate HEPBURN ET AL. 8 of 18 10.1029/2020JE006531 several resurfacing ages on a surface with a downturn in crater frequency marking the shift from one isochron to another (Figure 4d; Michael & Neukum, 2010;Michael, 2013).
Surfaces that have undergone successive episodes of lava deposition are an optimal landscape from which to determine multiple resurfacing ages (Fassett, 2016), and the AHv lava deposit mapped here consists of several superposed lobate flow units (Tanaka et al., 2014). Small craters are easier to remove from a surface than large craters (Fassett, 2016;Michael & Neukum, 2010) and successive lava flows would obliterate all craters on a surface up to a maximum diameter determined by the thickness of each lava flow (Fassett, 2016). In this manner, it is possible to derive multiple ages from the crater population on the AHv lava flow (Section 4.3, Figure 4c). For GLDA depressions, which could have undergone a more complex sequence of resurfacing, it is also possible to fit several isochrons to diameter ranges in some cases (Figure 4b). However, not all geological resurfacing processes act over discrete time intervals and may act to remove craters of all diameters from a surface quasi-continuously (Fassett, 2016). It is difficult to disentangle the precise (potentially highly complex) resurfacing histories of GLDA depressions, particularly following the deflation of GLDAs, and we do not know to what extent quasi-continuous resurfacing may have acted to remove craters. Therefore, in this instance, we are only interested the rightmost of these fits (i.e., the range of largest craters analyzed) to approximate the minimum formation age of GLDA depressions. The simple interpretation is that ages derived from this rightmost isochron for GLDA depressions and mesa tops approximate the most recent time at which these surfaces were last ice-free. We anticipate that age(s) derived for the lava deposit approximate the time at which successive lava deposits ceased to flow (see Section 1). We elaborate on a more sensitive interpretation of the model ages for each landform in our results (Section 4.3).
CraterStats 2.0 automatically calculates model age and the associated statistical uncertainty (Michael & Neukum, 2010). Error margins on model ages assigned by CraterStats 2.0 reflect the statistical uncertainty of both CSFD measurements and the counting area. This "±error" (Figures 4b and 4d and Figures S2-S49) is based on the Poission distribution model of random crater impacts (Michael & Neukum, 2010), describing how observed crater counts on a surface deviate statistically from the number of crater counts expected from the true age of a surface for given chronology and production functions. The age error is derived from the probability density function of the true age constructed from this model. We note this statistical uncertainty is independent of the potential sources of uncertainty arising from differences in target properties, resurfacing and neglects systematic errors associated with the chronology and production functions.

GLDA Mapping and Morphometry
In total, we mapped 37 GLDA depressions in Kasei Valles. Of these, 17 abutted the northernmost edge of the lava flow and the remaining 20 were radial GLDA depressions encircling isolated mesas within the lava flow. The total area of the GLDA depressions (including the area of any internal mesas) is 10,689 km 2 . The average runout of GLDA depressions is 3.17 km (s.d. 4 km); the distribution of runout length is heavily skewed toward lower values. The average relief of each scarp/mesa (the standard deviation of elevation within each GLDA depression) is 150 m. Within GLDA depression, there was no clear evidence of landforms related to glacial action such as streamlined flow-parallel lineations, flutes, or moraines. The area, relief, runout length, and center coordinates of each GLDA depression is reported in Table 1 of the online repository linked to this study (Hepburn et al., 2020a).

Glaciological Modeling
We modeled the surface and volume of 35 GLDAs for nine values of τ y from 16 to 200 kPa as shown in Figure  5. At τ y = 22 kPa-the mean τ y for contemporary LDAs (Karlsson et al., 2015)-our estimates predicts a of ice in the Kasei Valles region, with a mean thickness of 84.9 m for the GLDAs. The amount is equivalent to 1660 Gt of ice or a layer of water spread across the surface of Mars 12.6 mm thick. At τ y = 16 kPa, the lowest value of τ y used in this study, 1440 km 3 of ice is predicted, with a mean thickness of 75.8 m. Where τ y = 80 kPa, the largest value for τ y determined by Karlsson et al. (2015), total ice volume is 3450 km 3 , and mean ice thickness is 159 m. Finally, at the extreme upper value of τ y = 200 kPa, we predict a total volume of 5990 km 3 , with a mean thickness of 274 m. At all values of τ y , our model predicts the partial-total submergence of some HEPBURN ET AL. 9 of 18 10.1029/2020JE006531 mesa tops, even when unrealistically steep ice surfaces are removed. When τ y = 16 kPa, one radial GLDAs totally submerged the mesas within its boundaries and the remaining radial GLDAs partially submerged the lower elevation areas of mesas within GLDA depressions.

Dating
The diameter of 47,584 craters across 37 GLDA depressions, 11 mesas, and the Kasei Valles lava flow were recorded. The mean model age is 0.6 Ga for the GLDA depressions and 1.41 Ga for mesas (when each are treated as a collective unit). GLDA depression ages range from 1.9 to 0.2 Ga and mesa ages range from 2.8 to 0.6 Ga. Two model ages were derived for the lava flow ( Figure 4d): 1.3 ± 0.004 Ga (primary) and 2.5 ± 0.3 Ga (secondary). We interpret the younger of these ages, 1.3 Ga, which is consistent with previous estimates by Hauber et al. (2008), to correspond to the formation age of the uppermost lava flow deposits. The distribution of mesa ages extends above as well as below 1.3 Ga (Figures 6 and 7). Figure 6 shows that each GLDA depression is younger than the associated mesa(s) and younger than the primary model age of the Kasei Valles lava flow, with the exception of two GLDA depression surfaces, which are both reconcilable through age uncertainty. It is possible that our mesa and GLDA depression ages-being based on the approach of Michael and Neukum (2010)-may be biased by the size of count area (Fassett, 2016;Palucis et al., 2020;Warner et al., 2015). Although this deterministic approach to crater-based dating is used widely, we note a potential alternative in the probabilistic approach recently HEPBURN ET AL.  proposed by Palucis et al. (2020), but this was not pursued herein. The CSFD statistics are reported in Table  2 of the online repository linked to this study (Hepburn et al., 2020a).

Geomorphology
GLDA depressions delineate the spatial extent of GLDAs at the time the lava deposits cooled and solidified in Kasei Valles. Each GLDA depression is marked by a single clearly defined outer boundary consisting of decameter-scale crenulation, with little variation evident in the morphology of individual depressions. In order to compare the aspect ratio of contemporary mid-latitude LDAs and GLDAs, we measured both the relief and runout length of GLDA depressions. The range of runout length is 0.4 km-22.6 km, with a mean of 3.21 km (s.d.

Glaciological Modeling
Our modeling of GLDAs in Kasei Valles indicates GLDA depressions were occupied by relatively thin ice masses between 75 and 160 m thick with a total volume of 1,440-3,450 km 3 (where τ y = 16-80 kPa). Radar investigation and modeling of contemporary LDAs suggest that they are ∼200-700 m thick (Holt et al., 2008;Karlsson et al., 2015;Parsons & Nimmo, 2009;Plaut et al., 2009). The predicted mean ice thickness of GL-DAs, even where τ y >80 kPa ( Figure 5), are at the lower bounds of contemporary LDA thickness estimates. We interpret the difference in terms of the small planform dimensions of GLDAs relative to contemporary LDAs. However, our model predicts the partial to total submergence of mesa tops even at low values of τ y and when unrealistically steep ice surfaces are removed.

Model Ages
Piecing together the glacial history of the Kasei Valles region from the CSFD ages determined for the different surfaces requires careful consideration of their impact histories. Before interpreting our model ages, we must first consider the inheritance of preserved craters that predate the accumulation of GLDAs and the uppermost lava flow deposits. Typically, it is assumed that model ages correspond to the time since a geologically homogeneous surface was most recently free of craters, be that the formation age of the surface or the time since the last major resurfacing event (Michael & Neukum, 2010). This assumption is important here, as the lava flow, mesa tops, and the GLDA depressions are not bedrock surfaces (Dumke et al., 2010). Every layer in this sequence has been exposed to cratering, and the apparent modeled age of the uppermost surface depends how much it has been modified (or how well reset/cleaned) by successive geological processes. Crater inheritance is an unknown source of uncertainty, and their inclusion in CSFDs would raise the model ages of surfaces. Below, we consider what potential inheritance means for each type of dated surface.
The lava flow is the most expansive unit dated here and preserves the extent of GLDAs at a narrow time increment. Because emplacement of the lava unit drowns the preexisting surface, we anticipate the inheritance of craters predating the uppermost lava deposit (>1.3 Ga) to be minimal, limited to large craters deeper than the thickness of the lava flow. Of the two ages derived, we therefore associate the rightmost limb of Figure 4d (2.5 ± 0.3 Ga), based on 101 craters spanning 1,000-4000 m, with craters too large to be totally eradicated by the uppermost lava deposit, which are instead inherited from older lava deposits. We associate the younger of two model ages (1.3 ± 0.004 Ga), based on 1809 craters spanning 354-707 m, with the formation age of the uppermost lava deposit (a lava ) or the minimum time since the lava solidified and allowed the accumulation of craters. Our modeled a lava corresponds closely to the age estimated by Hauber et al. (2008) and the secondary age is consistent with the older ages (3.4 ± 0.8 Ga) found by Dumke et al. (2010) for the entire unit (extending west beyond our mapped lava deposits) with a large associated error.
Next, we consider the impact histories of the GLDA depressions. In the simplest model described in Figure  2, we envisage a system in which the deflation of a single GLDA comprehensively resets/cleans the surface during an episode of total deglaciation. With zero crater inheritance, the model age, a ghost , or rightmost limb for a GLDA depression corresponds to the cumulative length of time that depression has been exposed to cratering since the recession of the GLDA. The model ages, a ghost , therefore dates the exposure of the GLDA depression or (equivalently) the timing of this former ice loss after complete GLDA recession. With zero crater inheritance, a ghost dates deglacation in Kasei Valles exactly. In contrast, in the likely event of non-zero levels of crater inheritance on GLDA depressions, a ghost would be an overestimate of the time at which the GLDAs receded.
The chances of crater survival beneath and following the demise of GLDAs depends on the dynamics of the former ice masses, much of which is unknown. Our modeling of ice volume indicates each GLDA depression was occupied by a GLDA with a mean thickness of ≥75 m Section (3.2) at the time the lava deposit was formed ∼1.3 Ga (Section 3.3). For the duration of ice occupation, decameter thick ice will have shielded the subglacial surface from the majority of impacts, preventing the formation of new craters. However, GLDAs submerged preglacial craters which may be present as inherited craters on GLDA depressions. There are two conceivable processes by which GLDAs can modify the preglacial inheritance of craters: (i) the basal erosion of craters beneath principally warm and wet based GLDAs and (ii) the burial of craters beneath former supraglacial debris as GLDAs receded. The removal of preglacial craters by these processes is interesting to consider, particularly at smaller crater diameters, but it is unlikely that all preglacial craters were removed by glacial processes. Nonetheless, we reiterate that the inclusion of crater inheritance in our CSFDs does not affect our use of a ghost to date deglaciation in Kasei Valles. Our model ages of GLDA depressions, based on the rightmost limb, yield minimum exposure ages of the depression surfaces irrespective of crater inheritance. Below, we elaborate on our reasoning and consider whether our CSFDs provide additional insights into the thermodynamics of former GLDAs: 1. Warm-based glaciers with a subglacial drainage system on Earth can cause high erosion rates (Ugelvig et al., 2016). If GLDAs were warm-based, subglacial erosion could have removed significant proportions of preglacial craters during glacial occupation. However, our results indicate GLDA depressions are devoid of landforms suggestive of wet-based glaciation (Section 3.1) and glacial activity elsewhere on Mars is expected to be principally cold-based (Fassett et al., 2010). Putative eskers have been linked to warm-based glaciers in areas of especially high geothermal heat flux in Tempe Terra, northwest of Kasei Valles (Butcher et al., 2017), but such landforms have not been identified in our study region. Although it is possible a landform record of warm-based glaciation has been eroded or is otherwise no longer visible, extensive warm-based glaciation during the early-late Amazonian on Mars would require a much warmer climate (Forget et al., 2013). Even during high-obliquity phases, climate projections anticipate only modest temperature increases (∼15 K; Forget et al., 2013) and we therefore anticipate GLDAs were principally cold-based. Submerged craters beneath cold-based GLDAs would likely remain relatively undisturbed during ice occupation. 2. Following the recession of GLDAs, we expect significant modification of subglacial craters by supraglacial debris. Our results indicate that following lava emplacement ∼1.31 Ga, GLDAs persisted for >200 Myr until ∼1.1 Ga ( Figure 7a). As obliquity driven climatic cycles are known to vary on ∼10 1 -10 2 kyr timescales (Laskar et al., 2002), in a time of 200 Ma GLDAs will have experienced many such cycles, including periods unfavorable for the long-term survival of surficial ice in Kasei Valles. To survive through these cycles, GLDAs were likely buried beneath a supraglacial mantling deposit. Contemporary VFFs exist under such conditions, and any (subpolar) surfical ice on Mars would rapidly sublimate today. All contemporary mid-latitude VFFs have a relativity thin (∼10 m) supraglacial debris layer that retards sublimation (Baker & Carter, 2019;Boynton et al., 2002;Feldman, 2004;Plaut et al., 2009). Over time, this material is incorporated englacially and radar observation of LDAs and LVFs indicates their englacial debris content is up to 10% (Holt et al., 2008). Debris on and within mid-latitude VFFs principally consists of headwall-derived deposits and atmospherically deposited mantle (Baker & Carter, 2019). Each GLDA depression surrounds mesas or abuts scarps, similar to those throughout the mid-latitude domain of VFFs, which would provide debris input for analogous supraglacial and englacial deposits. We can estimate the thickness of material deposited on GLDA surfaces by assuming that the debris content and supraglacial debris thickness of GLDAs approximated that of contemporary LDAs. The melting of ice with a 10% debris content would lead to ≈1 m of lag deposits for every 10 m of surface lowering (Baker & Carter, 2019). Accordingly, we estimate at least 7.5 m of sublimation lag (based on an ice thickness of 75 m at τ y = 16 kPa) on GLDAs will have been generated during their recession, supplemented by supraglacial deposits up to 10 m thick. Beneath ∼20 m of debris, many smaller craters would be entirely buried. A key control on the likelihood of burial is crater depth, d, which scales with crater diameter, D (Caprarelli, 2014;Holsapple, 1987;Robbins & Hynek, 2012). Following a theoretical scaling relationship d = 0.084D 1.245 (Robbins & Hynek, 2012), we estimate that craters where D ≤ 80 m would be totally buried beneath ∼20 m of debris generated as GLDAs receded. This scaling relationship is based on fresh, unmodified craters (Robbins & Hynek, 2012). Preglacial crater populations may already have been exposed to modification and infilling, which would reduce the effective depth of craters without altering crater diameter. We anticipate that modified craters with D larger than 80 m would have been buried beneath ∼20 m of debris deposited following GLDA recession (where τ y = 16 kPa), extending into our fitted CSFD range between D min ( ≈62.5 m) and D max ( ≈500 m).
All GLDA depressions (with two accountable exceptions; Section 3.3) have model ages younger than the lava formation age (1.31 Ga; Figure 7a). Any crater inheritance in the GLDA depression crater counts would strengthen this finding. As the supraglacial deposit left by recession of GLDAs is expected to be thin, it is unlikely to have totally reset/cleaned the depressions of preglacial craters. Therefore, we anticipate some crater inheritance in our CSFD counts, especially within larger diameter bins. On multiple CSFD plots (Figures 4a and 4b and Figures S2-S6, S8, S10-S16, S18-S21, S23, S24, S26, S29, S35, and S36), the largest diameter bins appear to be shifted to the right of the fitted isochron, instead tracking a tangential isochron, which may reflect the presence of inherited craters within these bins. However, in most cases, the associated uncertainty (as reflected by vertical error bars) is large and intersects with the fitted isochron, so these bins are included within the fitted diameter range. With no crater inheritance, the survival duration, a ice , of GLDAs (up to their recession) since the time of lava emplacement (1.3 Ga), given by a ice = a lava −a ghost , is a minimum estimate. Our results for individual GLDA depressions (disregarding the two anomalous depressions which appear older than 1.3 Ga) give a ice ≈ 0.2-1.1 Gyr, implying that GLDAs existed for a substantial period of time before their demise. GLDAs also existed for some unknown duration prior to their demise.
Finally, we do not know exactly how and when the mesas in the Kasei Valles region were formed. Certainly they could not have formed after the lava flow 1.3 Ga because we cannot conceive any physically plausible mechanisms for mesas to appear in the middle of our GLDA depressions and grow out of nothing. Indeed, six mesas are dated to be older than 1.5 Ga. However, five other mesas have modeled surfaces younger than 1.3 Ga and one appears to be as young as 0.6 Ga. It therefore appears that post-formation resurfacing must have removed some of their craters. Our modeling indicates the partial-total submergence of many mesa tops by GLDAs, but it is not clear if mesa resurfacing is related to GLDAs themselves. Beneath cold-based ice, craters are unlikely to be significantly altered, and instead, as surfaces were exposed following the loss of ice these younger mesa surfaces were exposed to other mechanisms of resurfacing such as eolian activity.
In any case, all mesas have model ages older than those of the associated GLDA depression (Figure 6), consistent with the landform succession illustrated in Figure 2 and the long-term survival of GLDAs following lava emplacement.

Discussion and Conclusions
Our modeling indicates between 1,440 and 3,450 km 3 of ice existed in Kasei Valles at the time of lava formation (for τ y = 16-80 kPa). Our CSFD dating indicates this lava deposit cooled and solidified ∼1.3 Ga in our study region, contributing new information for the history of volcanic activity in Kasei Valles. Our model shows that when the lava deposit formed the mean ice thickness of contemporaneous GLDAs was 75-160 m. The GLDA ice volume was 1%-2% of the estimated total volume of today's mid-latitude ice on Mars (Karlsson et al., 2015) or 43%-104% of the ice contained within the glaciers and ice caps of Iceland, Earth. However, the ice did not vanish immediately following the lava emplacement event. CSFD dating of the GLDA depressions and the mesas contained therein show that the GLDAs survived for a minimum additional time (a ice ) ranging from 0.2 to 1.1 Gyr (Section 4.3). Implicit to this interpretation of a ice is that each GLDA retreated in one discrete episode of recession, without any glacial reaccumulation/readvance. Nevertheless, it is possible that GLDA depressions hosted multiple episodes of glacial activity over the last ∼1 Ga. By considering this possibility of polyphase glaciation, and the unknown influence of crater inheritance, we outline three alternative models for the history of regional ice loss in Kasei Valles: 1. At first glance, the large spread of a ghost evident from Figure 7a suggests the recession and eventual demise of GLDAs was asynchronous throughout the Kasei Valles region and had occurred relatively continuous over a period of 1 Gyr between 1.1 and 0.2 Ga. However, to explain this asynchronous scenario of regional deglaciation, it is necessary to invoke mechanisms besides climate forcing to drive glacier recession because climate forcing should modulate glacier mass balance in a somewhat uniform way within such a small region as Kasei Valles. 2. If GLDAs were strongly coupled to climate, then they would have receded synchronously. This would reasonably be expected if climate was perturbed by a major orbitally driven shift. In such a scenario, and assuming no crater inheritance on GLDA depressions, each value of a ghost should be identical, which is not indicated by our CSFD dating (Figure 7a). However, the variation in a ghost could reflect spatial discontinuity in the magnitude of crater inheritance, basal erosion, or in the post-formation modification rates of craters. In this interpretation, the true maximum age of the synchronous deglaciation (of all GLDAs) across the region would be given by the smallest value of a ghost (∼200 Ma), yielding a ice >1 Gyr, that is, the survival duration is even longer than our baseline estimate (≥0.2 Gyr). This would indicate that some major climatic shift occurred around ∼0.2 Ga, well beyond the period of reconstructible astronomical forcing (Laskar et al., 2004).
3. Finally, each GLDA depression may have hosted multiple phases of glacial activity since lava was emplaced 1.3 Ga. Obliquity cycles at ∼10 1 -10 2 kyr timescales drive climate variations (Laskar et al., 2002(Laskar et al., , 2004, and depending on the sensitivity of GLDAs to Martian climate, it is possible that the Kasei Valles region has undergone multiple phases of glacial activity over the last 1.31 Ga. In this third scenario of polyphase glaciation, a ghost records the maximum cumulative time that each GLDA depression has been free of ice between cycles of glacier growth and recession, rather than a maximum age estimate of the timing of its final deglaciation. While we cannot disregard scenarios 2 and 3, based on available evidence, we consider scenario 1 to be the most likely model of glacial history in Kasei Valles. There is no clear and obvious record of polyphase glaciation in the Kasei Valles landform assemblage, nor do we consider that crater inheritance exerts a strong control on the variability in model ages (Section 4.3). Although polyphase glaciation is evident elsewhere on Mars (e.g., Baker et al., 2010;Brough et al., 2015;Dickson et al., 2008;Hepburn et al., 2020b;Levy et al., 2007;Morgan et al., 2009), each GLDA depression studied herein has a single sharply defined boundary with the lava flow. Within or beyond GLDA depressions, there are no clearly identifiable features such as moraines or eskers, which may be linked to multiple phases of glaciation. It is possible that such a record has since been removed or otherwise is no longer visible or may not have been formed at all beneath what we anticipate were cold-based GLDAs. However, as with present day VFFs, we anticipate GLDAs persisted beneath relativity thin (∼10 m) supraglacial debris layers (Baker & Carter, 2019). This debris layer would have limited atmospheric mass balance exchanges such that, following an initial period of accumulation, the decoupling of GLDAs from climate would have enabled the long-term survival of buried ice in Kasei Valles, even as Martian climate varied dramatically. Therefore, we anticipate that local variations in the thickness of glacial ice and supraglacial debris layers exert the strongest controls on the longevity of GLDAs.
Our reconstruction of GLDAs raises interesting questions regarding their initial formation. Our prediction that mesas were submerged by palaeo ice may potentially support the idea that GLDAs formed following the decay of a thick, continuous, and regional-scale ice sheet over the Kasei Valles region. Although we do not find any direct evidence for such an ice sheet (which would well predate the emplacement of the AHv lava unit 1.3 Gyr ago), this configuration has been used by Fastook et al. (2014) as an initial condition for their model of the historical evolution of contemporary mid-latitude LDAs.
The GLDAs reconstructed here represent some of the oldest glacial ice masses on Mars and their dating and modeling yield a unique insight into low-latitude ice in the late Amazonian, beyond the known extent of contemporary ice-rich VFFs and before the oldest estimates of mid-latitude glacial activity. We provide an indirect measurement of ice survival times on Mars, indicating GLDAs persisted in the Kasei Valles region for ≥1 Gyr prior to their demise during the late Amazonian on Mars.

Data Availability Statement
All data derived from our mapping (the inventory of GLDAs and their morphometry, and ArcGIS shapefiles of GLDA boundaries and all craters discussed), all data used in our dating procedure (impact crater counts), and the model script are given in the online repository linked to this study (Hepburn et al., 2020a). The crater size-frequency distribution plots for our 37 GLDAs, 11 mesas, and the lava flow are provided in the supporting information.