# The internal structure of a debris-covered glacier on Mars revealed by gully incision

F.E.G. Butcher, N.S. Arnold, Susan J. Conway, D.C. Berman, J.M. Davis, M.R. Balme

► To cite this version:

F.E.G. Butcher, N.S. Arnold, Susan J. Conway, D.C. Berman, J.M. Davis, et al. The internal structure of a debris-covered glacier on Mars revealed by gully incision. Icarus, 2024, pp.115717. 10.1016/j.icarus.2023.115717. hal-04262257

HAL Id: hal-04262257

https://hal.science/hal-04262257v1

Submitted on 27 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# The internal structure of a debris-covered glacier on Mars revealed by gully incision

F.E.G. Butcher<sup>a, *</sup>, N.S. Arnold<sup>b</sup>, S.J. Conway<sup>c</sup>, D.C. Berman<sup>d</sup>, J.M. Davis<sup>e</sup>, M.R. Balme<sup>f</sup>

## ARTICLE INFO

**Keywords:**

Mars ice

Mars gullies

Structural glaciology

Ice flow modeling

In situ resource utilization

## ABSTRACT

Viscous flow features (VFFs) in Mars' mid latitudes are analogous to debris-covered glaciers on Earth. They have complex, often curvilinear patterns on their surfaces, which probably record histories of ice flow. As is the case for glaciers on Earth, patterns on the surfaces of VFFs are likely to reflect complexities in their subsurface structure. Until now, orbital observations of VFF-internal structures have remained elusive. We present obser- vations of internal structures within a small, kilometer-scale VFF in the Nereidum Montes region of Mars' southern mid latitudes, using images from the Context Camera (CTX) and High Resolution Imaging Science Experiment (HiRISE) instruments on Mars Reconnaissance Orbiter. The VFF-internal structures are revealed by a gully incision, which extends from the VFF headwall to its terminus and intersects curvilinear undulations and a crevasse field on the VFF surface. Near to the VFF terminus, the curvilinear VFF-surface undulations connect to the VFF-internal layers, which are inclined and extend down to the VFF's deep interior, and possibly all the way to the bed. Similar structures are common near to the termini of glaciers on Earth; they form under ice flow compression where ice thins and slows approaching the ice margin, and ice flow is forced up towards the surface. We performed 3D ice flow modeling which supports this analogy, revealing that the inclined VFF-internal structures, and associated curvilinear structures on the VFF surface, are located in a zone of strong ice flow compression where ice flow deviates upwards away from the bed. The inclined VFF-internal structures we observe could represent up-warped VFF-internal layering transported up to the surface from the VFF's deep interior, or thrust structures representing debris transport pathways between the VFF's bed and its surface. Our observations raise numerous considerations for future surface missions targeting Mars' mid-latitude subsurface ice deposits. Inclined layers formed under flow compression could reduce the requirement for high-cost, high-risk deep drilling to address high-priority science questions. They could allow futures missions to (a) access ice age sequences for palaeoenvironmental reconstruction via shallow sampling along transects of ice surfaces where layers of progressively older age outcrop, and/or (b) access samples of ice/lithics transported to shallow/surface positions from environments of astrobiological interest at/near glacier beds. However, our observations also raise considerations for potential drilling hazards associated with structural complexities and potential dust/debris layers within subsurface ice deposits on Mars. They highlight the importance of characterizing VFF-surface structures, and their relationships to VFF-internal structure and ice flow histories ahead of ice access missions to Mars.

## 1. Introduction

Viscous flow features (VFFs) in Mars' mid-latitudes are widely interpreted as debris-covered glaciers based on geomorphic (e.g., Head

et al., 2005), geophysical (e.g., Holt et al., 2008; Plaut et al., 2009; Petersen et al., 2018), and climatological (e.g., Madeleine et al., 2009) evidence, and direct imaging of ice exhumed by fresh impacts (Dundas et al., 2021). The term VFF covers a family of features that can be

<footnote><sup>a</sup> Department of Geography, University of Sheffield, Winter Street, Sheffield S10 2TN, UK</footnote>

<footnote><sup>b</sup> Scott Polar Research Institute, University of Cambridge, Lensfield Road, Cambridge CB2 1ER, UK</footnote>

<footnote><sup>c</sup> Nantes Université - Université d'Angers - Le Mans Université, CNRS UMR 6112 Laboratoire de Planétologie et Géosciences, France</footnote>

<footnote><sup>d</sup> Planetary Science Institute, 1700 E. Ft. Lowell Rd., Suite 106, Tucson, AZ 85719, USA</footnote>

<footnote><sup>e</sup> Department of Earth Science and Engineering, Imperial College London, London, SW7 2AZ, UK</footnote>

<footnote><sup>f</sup> School of Physical Sciences, The Open University, Walton Hall, Milton Keynes MK7 6AA, UK</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: f.butcher@sheffield.ac.uk (F.E.G. Butcher).</footnote>
grouped into several morphologic subtypes including: small, kilometer-scale glacier-like forms (e.g., Arfstrom and Hartmann, 2005; Hubbard et al., 2011; Souness et al., 2012), and larger (tens of kilometers-scale) lobate debris aprons, lineated valley fills, and concentric crater fills (e. g., Squyres, 1978, 1979). VFFs are key targets for future in situ resource utilization (ISRU) of water by human missions to Mars' surface (Hoffman et al., 2017; Morgan et al., 2021; I-MIM MDT, 2022). Constraining the internal structure and geotechnical properties of Mars' ice deposits is

essential for planning for ISRU of water ice and the acquisition of scientific samples (e.g., ice samples/cores and lithic material) by future surface missions (e.g., Mars Ice Core Working Group, 2021; I-MIM MDT, 2022). Curvilinear structures evocative of viscous deformation on VFF surfaces suggest complex ice flow histories which inevitably manifest in VFF-internal structure, and could record past mass balance processes and climatic fluctuations (e.g., Levy et al., 2021) since VFFs formed millions to hundreds of millions of years ago (e.g., Baker and Head,

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Context map and oblique view of the study site. (A) Mars Orbiter Laser Altimeter (MOLA) elevation map of Mars showing the location of Argyre Planitia and the extent of panel B (black box). (B) MOLA elevation map of Argyre Planitia showing the location of Nereidum Montes, our study site (yellow star; 51.24°W, 42.53°S), and the area in which Berman et al. (2021) obtained a surface age estimate for a population of VFFs. (C) Oblique view of the gully-incised viscous flow feature. Orthorectified High Resolution Imaging Science Experiment (HiRISE) image ESP_051036_1370 overlain on a HiRISE digital terrain model (see Section 3.1). (D) Schematic map overlain on the oblique view in C, with key landscape features colorized. HiRISE data credit: NASA/JPL/University of Arizona/NHM. MOLA data credit: USGS Astrogeology Science Center/Goddard Space Flight Center/NASA. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
2015; Hepburn et al., 2020; Berman et al., 2021). Orbital ground-
penetrating radar observations from the Shallow Radar (SHARAD) in-
strument on Mars Reconnaissance Orbiter (MRO) have identified re-
flectors at/near the beds of some larger (>10 km-scale) VFFs (Holt et al.,
2008; Plaut et al., 2009; Petersen et al., 2018; Berman et al., 2021;
Gallagher et al., 2021), but insights into their internal structure (e.g.,
layering) have remained elusive. This probably results, in part from
resolution limits of existing ground-penetrating radar instruments, as
well as limitations arising from topographic steepness/complexity,
surface roughness, and loss characteristics of target materials (e.g.,
Petersen et al., 2018).

We present new observations and 3D analyses of the internal struc-
ture of a small, kilometer-scale VFF in the Nereidum Montes region of
Mars' southern mid-latitudes (Fig. 1A–B), from images acquired by the
High Resolution Imaging Science Experiment (HiRISE) instrument on
MRO. VFF-internal layers are revealed in the wall of a gully (a hillslope
landform comprising an erosional alcove, a transportation channel, and
a depositional fan) which has incised through the VFF from its headwall
to its terminus (Fig. 1C). The internal layers connect to curvilinear
structures visible on the VFF surface. Near to the VFF terminus, the
layers transition from approximately bed-parallel to inclined structures
that dip up-glacier, extending down from the VFF surface to its deep
interior and possibly all the way to the bed. On Earth such configura-
tions of glacier-internal layers are commonly related to ice flow stress
regime, and often provide pathways for ice and debris transport from
glacier deep interiors/beds to shallower positions at/near glacier sur-
faces (Hambrey et al., 2005; Jennings and Hambrey, 2021). We perform
3D analyses (including dip and dip direction measurements) of the
structures observed within the gully-incised VFF in Nereidum Montes,
and perform 3D ice flow modeling to analyze the relationship of the
structures to ice flow stress regime.

We then explore the implications of the VFF-internal structures for
glacial landform development, and highlight key considerations for
scientific investigations, ice access, and ISRU by future surface missions
targeting mid-latitude ice deposits on Mars. We particularly consider the
possibility that the structures could have transported ice and lithics from
the VFF's deep interior and/or bed to its surface, thus potentially miti-
gating the requirement for deep drilling to access samples from these
scientifically important environments.

## 2. Study location

The incised viscous flow feature (51.24°W, 42.53°S), which was
noted by Berman et al., (2021, their Fig. 7), is located in western Ner-
eidum Montes, in Mars' southern mid-latitudes. Nereidum Montes is a
mountainous region, which forms the northern margin of the ancient
Argyre Planitia impact basin (Fig. 1). It hosts abundant deposits inter-
preted to comprise subsurface water ice, including VFFs and ice-rich
mantling deposits (Milliken et al., 2003; Berman et al., 2021). Subsur-
face ice has been detected in these deposits at sites within and near to
Nereidum Montes, both from subsurface reflections detected by
SHARAD (Berman et al., 2021) and in natural exposures formed by
active mass wasting within a gully (Khuller and Christensen, 2021).
Additionally, ‘ice consistency’ maps from the Subsurface Water Ice
Mapping project, which combine multiple orbital datasets and tech-
niques to derive a measure of the consistency of those datasets with the
presence of subsurface ice across regions, find moderate ice-consistency
values across Nereidum Montes, in both the 1–5 m and > 5 m depth
zones (Putzig et al., 2023).

VFFs across Mars' mid latitudes are geologically young, having
modeled surface ages on the order of millions to hundreds of millions of
years (i.e. mid-to-late Amazonian) (Baker and Head, 2015; Hepburn
et al., 2020; Berman et al., 2021). The VFF at our study site is too small,
and has too few superposed impact craters, to extract a reliable surface
age estimate from impact crater size-frequency statistics on its surface.
However, combined statistics obtained by Berman et al. (2021) over a

population of 23 VFFs in Nereidum Montes (~250 km to the south of our
study site; Fig. 1B) returned a combined VFF model age estimate of 2.7
+/- 0.8 Myr. This suggests that VFFs in Nereidum Montes formed
during the very recent stages of the Late Amazonian epoch. The incised
VFF is one of a population of small, kilometer-scale VFFs (e.g., Milliken
et al., 2003; Arfstrom and Hartmann, 2005; Berman et al., 2021) that are
typically younger than the larger lobate debris apron, lineated valley fill,
and concentric crater fill VFF morphological subtypes (10s km-scale,
100 s Myr old). Many small VFFs, such as the incised VFF we analyze
here, are thought to have formed by downslope flow of ice-rich deposits
(with some component of admixed dust) which form a meters to hun-
dreds of meters-thick mantling deposit over vast areas of Mars' mid-to-
high latitudes (e.g., Mustard et al., 2001; Head et al., 2003; Conway and
Balme, 2014). These ice-rich mantling deposits likely formed by airfall
during multiple cyclical increases in Mars' spin-axis obliquity in the last
few million years, similar to glaciation driven by Milankovitch cycles on
Earth (Head et al., 2003; Laskar et al., 2004). Thicker mantling deposits
and/or those on steep slopes have undergone glacier-like flow, forming
small VFFs (e.g., Berman et al., 2021). Their surficial dust/debris cover
preserves the underlying ice into the present day.

## 3. Materials and methods

### 3.1. Basemap data and digital terrain model

Structures on the surface of the incised VFF are visible in both 6 m/
pixel Context Camera (CTX; Malin et al., 2007) and 25 cm/pixel High
Resolution Imaging Science Experiment (HiRISE; McEwen et al., 2007)
images. The VFF-internal structures exposed in the gully wall are visible
at the limits of CTX resolution, in image
P14_006572_1367_XN_43S051W (Fig. 2). The VFF-internal structures
are clearly visible in panchromatic HiRISE images
(ESP_051036_1370 RED and ESP_015947_1370 RED), and color varia-
tion between layers is visible in 3-band enhanced color (IRB and RGB)
HiRISE images (e.g., ESP_051036_1370_IRB; Fig. 2). The HiRISE false-
color image used herein is the orthorectified image
ESP_051036_1370_IRB_A_01_ORTHO.

Following the procedure of Kirk et al. (2008), we used the USGS
Integrated Software for Images and Spectrometers (ISIS 3) and BAE
Systems SOCET SET to generate a 1 m/pixel digital terrain model (DTM)
and a 25 cm/pixel orthorectified image from HiRISE stereo pair images
ESP_051036_1370 (0.264 m/pixel) and ESP_015947_1370 (0.254 m/
pixel). Based on the approach of Kirk et al. (2008), we estimate the
vertical precision of our HiRISE DTM to be 0.20 m, given a stereo
convergence angle between the input images of 14.8°, and assuming an
RMS pixel matching error of 0.2 pixels.

To minimize distortions for our structural analyses (Section 3.2) and
ice flow modeling (Section 6), we re-projected the HiRISE DTM and
orthorectified image into a sinusoidal projection centered on the VFF-
internal structures (51.25°W).

### 3.2. Structural analysis of VFF-internal and VFF-surface structures

Using ESRI ArcGIS 10.7, we measured the dips and dip directions of
the VFF-internal structures, and the VFF-surface structures to which they
connect, as follows. We first digitized as lines the two structures for
which reliable contacts could be traced clearly from the gully wall
exposure onto the VFF surface (Fig. 3A). While multiple structures have
traces on both the surface and the gully wall (Figs. 2B, 3A), the two
sampled structures had particularly clear contacts with adjacent layers,
allowing reliable, reproducible measurements of dip and dip direction.
The upflow layer is hereafter referred to as $L_{Up}$ and the downflow layer
as $L_{Down}$; we digitized along the contact between these layers and their
adjacent down-flow layers. We also digitized the contact between the
VFF margin and the proglacial surface, continuing around the junction
between the VFF terminus and the gully incision, to infer the slope of the
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** (A) Context Camera image (P14_006572_1367_XN_43S051W) showing the gully-incised VFF and key landsystem components discussed in the text. Black dotted arrows show the approximate inferred ice flow directions for the northern and southern flow units of the VFF (NFU and SFU, respectively). The white dotted line shows the VFF terminus. The location is indicated by the black box in the inset map (top left), showing the regional topographic context in MOLA elevation data. (B) Orthorectified false-color HiRISE image (ESP_051036_1370_IRB_A_01_ORTHO) showing the intersection of curvilinear undulations on the VFF surface with VFF-internal structures visible in the northern wall of the gully chute. Extent shown by black box in panel A. MOLA data credit: USGS Astrogeology Science Center/Goddard Space Flight Center/NASA. CTX image credit: NASA/JPL/MSSS/UfoA. HiRISE data credit: NASA/JPL/University of Arizona.</figcaption></figure>

VFF bed at the ice margin (Fig. 3B). We then converted the lines to points spaced at 1 m and extracted the corresponding elevation values from the HiRISE DTM. Next, we fitted a planar surface through the points with a linear least-squares regression using the ArcGIS 'Trend' tool. We then obtained measurements of the slope and aspect of these surfaces to estimate the dip and dip directions of each of the structures (Fig. 3B). We performed five additional blind repeat measurements for each structure (but with broadly consistent start and end points, to ensure similar numbers of points contributing to each planar fit), taking the average dips and dip directions of all six measurements as the reported values, and their standard deviations as the reported uncertainties. We validated our measurement of the bed slope at the VFF margin by also fitting a linear plane to the exposed proglacial surface immediately beyond the VFF margin (Fig. 3B).

### 3.3. Ground-penetrating radar

Existing orbital ground penetrating radar instruments at Mars are of limited use for studying the specific incised VFF we analyze here. The incised VFF is similar in scale to the 0.3–1 km along-track resolution of the SHARAD instrument on MRO (Seu et al., 2007), and smaller in scale than the 5–10 km along-track resolution of the Mars Advanced Radar for Subsurface and Ionosphere Sounding (MARSIS) instrument on ESA’s Mars Express orbiter (Picardi et al., 2005). The incised VFF is also surrounded by high-relief terrain (which can cause false subsurface reflections, called ‘clutter’, from off-nadir reflections), and has a sub-optimal orientation relative to the existing SHARAD radar track S_00717201 that crosses it.
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** VFF-internal structures exposed in the northern wall of the gully chute. (A) Orthorectified HiRISE false-color image ESP_051036_1370_IRB_A_01_ORTHO of the intersection of the curvilinear undulations with the northern wall of the gully chute. The two curvilinear structures which have clear traces into the gully wall ($L_{Up}$ and $L_{Down}$) are indicated with black arrows. The white arrows indicate examples of structures that appear subparallel to $L_{Up}$ and $L_{Down}$ but have less distinct traces on the VFF surface and/or the gully wall for reliable dip and dip direction measurements. The inset image (extent shown by black box) shows a close-up of $L_{Down}$, which illustrates the relatively red color signature of the VFF internal structures compared to adjacent, relatively blue materials. It also shows the presence of meter-scale clasts on the gully wall. (B) Colorized HiRISE elevation map overlain on orthorectified panchromatic HiRISE image of the VFF margin showing measurements of dip and dip direction (white symbols) for VFF-internal structures $L_{Up}$ and $L_{Down}$, for the VFF bed at the VFF margin, and for the proglacial surface. Black lines show one of six traces per structure, to which we fitted linear planes to obtain the measurements (except for the proglacial surface, for which a linear plane was fitted to all DTM points within the area outlined in black). $L_{Up}$, the bed, and the proglacial surface all dip down towards the West, and $L_{Down}$ dips down towards the North East. Orthorectified HiRISE image ESP_051036_1370 overlain by the colorized HiRISE DTM. HiRISE data credit: NASA/JPL/University of Arizona/NHM. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

# 4. Observations

## 4.1. VFF surface features

The incised VFF is located on the west-facing hillslope of a ~2 km high massif (Fig. 1). The portion of the VFF to the north of the gully incision has a lobate, apron-like form with a surface area of 6.8 km², which extends ~2–4 km E–W from the hillslope over the plain to the west of the massif. The portion of the VFF to the north of the gully incision appears to comprise two units of flow. The northern flow unit (NFU in Fig. 2) is ~4 km long, with flow converging from two source

alcoves within the hillslope to the NE. The southern flow unit (SFU in Fig. 2) is ~2 km long and extends E-W from the face of the hillslope.

The VFF has a relatively smooth surface hosting small-scale (10s m long, meters-wide) linear ridges oriented approximately parallel to the inferred flow direction (Fig. 4A). There is a scattering of meter-scale rocks on the surface, but in general the grain size of the surface materials is below that which can be resolved in HiRISE images (Fig. 4A-B insets). Most striking among the VFF’s surface features are: (a) a set of nested, low-relief, curvilinear undulations near to the lateral and terminal margins of the VFF (Fig. 2A, 3, 4B); and (b) two fields of, quasi-linear ridges and fracture-like troughs, which occur in the upper
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Morphologies on the surface of the incised VFF, as seen in HiRISE images. Dotted arrows show inferred ice flow direction. (A) Small-scale linear ridges oriented approximately parallel to inferred VFF flow. Inset shows a scattering of meter-scale boulders on the surface. (B) A portion of the curvilinear undulations, showing the relatively red color signature of the undulation crests, and the relatively blue signature of the troughs. Inset shows a scattering of meter-scale boulders on the surface. (C) The northern field of quasi-linear flow-transverse ridges; small-scale flow-parallel linear ridges, similar to those in panel A, are also clearly visible here. (D) The southern field of quasi-linear flow-transverse ridges and troughs. Inset map in center shows locations of panels against CTX image P14_006572_1367_XN_43S051W. Panels A, C, and D are orthorectified HiRISE image ESP_051036_1370. Panel B is a false-color HiRISE image ESP_051036_1370_IRB. CTX data credit: NASA/JPL/MSSS/University of Arizona. HiRISE data credit: NASA/JPL/University of Arizona/NHM. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

reaches of both flow units and are oriented transverse to inferred VFF flow (Fig. 2A, 4C-D).

The curvilinear undulations create distinctive banding on the VFF surface (Fig. 2A, 4B), forming low-relief (<1 m) undulations in the HiRISE DTM. They originate close to the lateral margin of the northern flow unit, just above the outlet from its northernmost source alcove (Fig. 2A). Here, the curvilinear features are oriented approximately parallel to the inferred ice flow direction. Upon exiting the alcove, the undulations curve round to the south to form lobate flow-transverse bands perpendicular to the inferred ice flow direction, and parallel to the VFF terminus (Figs. 2A, 4B). They form two distinctive lobes

bounding the downslope reaches of the northern and southern flow units. At least two of the curvilinear undulations can be traced continuously across the two lobes, while others appear to pinch-out, or become indistinct in the interlobate zone. The curvilinear undulations have distinct visible roughness and color signatures; the ridges appear rougher and relatively red, and the troughs appear smoother and relatively blue (Figs. 2B, 4B). They are associated with a scattering of meter-scale boulders, which appear to be more abundant here than on the smoother portions of the VFF surface in its middle reaches (Fig. 4A).

The two fields of quasi-linear, flow-transverse ridges and troughs (Fig. 4C-D) occur in the upflow reaches of both of the VFF's constituent
flow units (Fig. 2A). The ridges and troughs in both fields have heights/ depths of ~1 m or less. The northern field (Fig. 4C) comprises only ridges, spaced ~25–50 m apart. It has a down-flow extent of ~300 m, and initiates ~200 m down-slope of the convergence between the two source alcoves of the northern flow unit. The southern field (Fig. 4D) is more extensive and complex. The upslope portion of the field, which initiates just ~50 m from the VFF headwall, comprises mostly ridges. The downslope portion comprises linear, fracture-like troughs. These features become less distinct downslope, but remain visible in HiRISE images up to ~750 m downflow of the headwall. The typical spacing between adjacent ridges/troughs is similar to the northern field.

## 4.2. Gully incising VFF

A gully incises through the southern flow unit of the VFF, from its headwall to its terminus. The gully originates with a ~2 km long, ~1 km wide, ~100 m deep, v-shaped alcove (Fig. 5, profile a–a') incised into the bedrock headwall above the VFF. Where it intersects the VFF downslope, it transitions to a v-shaped chute incised through the VFF (being ~300 m wide and ~40 m deep at profile b–b' in Fig. 5). The chute widens and becomes more flat-floored downslope (being ~600 m wide and ~40 m deep at profiles c–c' and d–d' in Fig. 5). The northern wall of the chute has a scattering of meter-scale boulders on its surface (Fig. 3A inset), but in general the grain size of materials within the gully wall is finer than that which can be resolved in HiRISE images. The chute terminates at the VFF's downslope margin. A depositional fan originates within the gully chute, and extends a further ~3 km beyond the VFF margin to the SW (Fig. 5). The fan is incised by a narrow (~60 m wide, ~10 m deep) sinuous channel, which originates within the chute and extends ~1 km beyond the VFF terminus, terminating in a smaller depositional fan superposed on the main gully fan deposit. The gully system is ~6 km long from the head of the alcove to the distal margin of

the fan deposits. Profile e–e' in Fig. 5 shows a longitudinal profile along the portion of the gully that is covered by the HiRISE DTM.

## 4.3. VFF-internal structures exposed in the gully wall

The northern wall of the gully chute intersects the curvilinear undulations near to the VFF terminus, and the southern field of quasi-linear flow-transverse ridges and troughs further upslope. Both types of structure can be traced from the VFF surface into the wall of the gully chute, thus providing a view of their subsurface manifestations within the VFF.

Several of the curvilinear undulations can be correlated to layers with distinct color and visible roughness signatures in the gully wall (Fig. 2B and Fig. 3A). As on the VFF surface, the layers correlated with the crests of the undulations generally appear rougher and redder, while those correlated to the troughs of the undulations generally appear smoother and bluer. The layers correlated with the undulation crests also appear to have relatively consistent thicknesses throughout (measured in planform, $L_{Down}$ appears ~20–30 m thick; Fig. 3A inset). Two of the curvilinear undulation crests correlate particularly clearly to layers in the gully wall, from which we could reliably obtain dip and dip direction measurements (Fig. 3). The upflow layer, $L_{Up}$, outcrops on the VFF surface ~500 m from the terminus (Fig. 3A). It appears to be part of a sequence of sub-horizontal layers exposed in the wall of the gully chute. $L_{Up}$ dips down towards the west at 6.9 ± 0.2° (dip direction 277 ± 7°, based on six fits to an average of 224 points each). $L_{Up}$ has a particularly rough texture compared to other structures in the gully wall and on the VFF surface (Fig. 3A). The VFF bed at the VFF margin dips down towards the west at ~12 ± 0.05° (dip direction 262 ± 0.2°, based on six fits to an average of 440 points each), demonstrating that $L_{Up}$, and associated subparallel layers, conform approximately to the slope and orientation of the bed (deviating from the slope of the bed by just 5°).

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Elevation profiles of the gully incision through the VFF, from the HiRISE DTM. Elevation profiles a–a’ to d–d’ are cross sections through the gully alcove and chute, and elevation profile e–e’ shows the long-axis of the gully. Note that profile e–e’ is truncated relative to the full extent of the gully because the head of the alcove and the terminus of the fan are not covered by the HiRISE DTM. The locations of the profiles are shown as white lines on CTX image P14_006572_1367_XN_43S051W. CTX image credit: NASA/JPL/MSSS/University of Arizona.</figcaption></figure>
Our measurement of the subglacial surface around the VFF margin is consistent with the slope of a linear plane fitted to the proglacial surface immediately beyond the VFF margin, which also dips down towards the west at ~12° (dip direction 262°, based on 58,317 points; Fig. 3B).

The second structure we measure, $L_{Down}$, occurs further downflow. It outcrops as a curvilinear undulation crest ~350 m up-flow from the VFF terminus (Fig. 3A). $L_{Down}$ appears to be one of a series of layers exposed in the wall of the gully chute, which tilt relatively steeply up to the ice surface from the VFF interior and are spaced tens of meters apart (Fig. 2B and 3A). $L_{Down}$ dips down towards the North East at ~11 ± 2° (dip direction 36 ± 3°, based on six fits to an average of 264 points each). This corresponds to a surface that is angled upwards at approximately 23° relative to the VFF bed (as illustrated schematically in Fig. 6), but we note that the dip-direction of $L_{Down}$ is not bed-parallel (deviating from the bed dip direction by 134°). Thus, we take 23° as a maximum value for the angle between $L_{Down}$ and the bed. Measured in planform within the wall of the gully chute, $L_{Down}$ is a ~20 m thick layer (Fig. 3A inset). It can be traced from the VFF surface near the VFF terminus, into the deep VFF interior (almost all the way to the floor of the gully chute, i.e., the inferred VFF bed). Based on a measurement from the HiRISE DTM, and assuming that the base of the northern wall of the gully chute is the VFF bed, the VFF is ~35 m thick where $L_{Down}$ intersects the surface.

In the upper reaches of the VFF, at least three of the quasi-linear flow-transverse ridges and troughs (two ridges, one trough) can be clearly connected to sub-vertical linear traces in the gully wall (Fig. 7). They are traceable through a vertical thickness of ~5–12 m of VFF materials. Hence, they appear to be shallow structures compared to those connected to the curvilinear undulations further downslope; (Fig. 3).

# 5. Interpretations

## 5.1. VFF surface features

We interpret the curvilinear undulations on the VFF surface as manifestations of VFF-internal layers outcropping beneath the VFF’s surficial dust/debris cover. On Earth, curvilinear banding on glacier surfaces, related to internal layering, commonly reorients from ice-flow parallel to ice-flow transverse approaching the ice margin (e.g., Jen- nings and Hambrey, 2021). The up-flow onset of the curvilinear

undulations on the incised VFF appears to coincide with a location where we can qualitatively predict a zone of lateral (i.e. flow-transverse) ice flow compression near to the northern lateral margin of the VFF, due to the topographic funneling of the northern flow unit out of its northern source alcove. Down-flow of this, the loss of lateral compression by topographic confinement likely gives way to flow-parallel compression induced by thinning (and slowing) of the VFF flow towards its margin. This is consistent with the reorientation of the curvilinear undulations as they approach the ice margin, where they become flow-transverse. We explore this interpretation further with an ice flow model in Section 6.

The true bulk properties (e.g., composition, grain size) of layers outcropping from the VFF to form the curvilinear undulations may not be visible due to the shallow burial of the ice surface that is required to preserve subsurface ice in Mars’ mid latitudes. However, given the observed continuation of the VFF surface structures into the wall of the gully chute, this covering layer is probably only a thin veneer; it is not thick enough to fully obscure meter-to-decimeter-scale textures of the underlying VFF’s structure. The relatively rough texture of the crests of the curvilinear undulations (and their connected VFF-internal struc- tures) could arise from higher debris concentrations and/or larger debris grain sizes within the shallowly buried layers that form them. Their relatively red appearance could arise from preferential trapping of dust where these rougher debris layers outcrop. Some of the meter-scale boulders near to the curvilinear undulations could have been released from the VFF-internal layers (as suggested for boulder bands analyzed by Levy et al., 2021), while others are likely to be supraglacial material that has been transported into this zone over the ice surface.

We interpret the two fields of quasi-linear, flow-transverse ridges and troughs (Fig. 4 C and D) as fields of crevasses (troughs) and supraglacial crevasse fills (ridges). Crevasses are fractures formed in regions of ice flow extension. Supraglacial crevasse fills comprise lithic materials that have infilled pre-existing crevasses. On Mars, this could originate from airfall (i.e., dust), windblown material (e.g., sand), and/or hillslope/ VFF-surface debris. Supraglacial crevasse fills can evolve into supra- glacial ridges where ice-surface lowering removes the bounding mate- rials (e.g., Bennett et al., 2000). The crevasse fills we observe are preferentially located close to the VFF headwall, with open, unfilled crevasses occurring further downslope. This suggests that the delivery of hillslope debris to the upper reaches of the VFF could have been a major contributor to crevasse-filling material.

The small-scale flow-parallel linear ridges, which occur extensively across the VFF surface (Fig. 4A and C), are different in terms of morphology, scale, spacing, and orientation compared to the curvilinear undulations and the quasi-linear, flow-transverse ridges and troughs. Their local orientations appear to deviate with VFF-surface slopes sug- gesting that they reflect the overall VFF flow pattern. We do not inter- pret them as surface manifestations of major ice-internal structure (e.g., layer outcrops); instead we suggest that they could represent the flow- induced fabric of the surface and near-surface materials.

## 5.2. VFF-internal structure

The dissection of the VFF, and its major surface flow deformation structures, by the gully chute (Fig. 2) indicates that gully incision postdates the last major phase of ice flow. We do not find evidence for exposed water ice within the gully chute, and consider it likely that the VFF-interior exposure became mantled by a layer of lithic material. This material could include airfall dust, lithics transported down the wall of the chute from the VFF surface, and perhaps lithics released from the VFF-internal structures themselves. As inferred for the curvilinear VFF- surface structures to which the VFF-internal layers connect (Section 5.1), the VFF-internal structures with relatively red and rough signa- tures could host higher lithic concentrations, larger lithic grain sizes, preferentially trap dust where they outcrop, or a combination of all of these factors.

The inclined layers (i.e., $L_{Down}$ and associated sub-parallel structures

<figure><img src="image_8.png" /></figure>

**Fig. 6.** Schematic of the gully wall exposure through the VFF, illustrating the inclined nature of $L_{Down}$ relative to $L_{Up}$ and the VFF bed. Black lines represent the structures for which measurements were taken. We note that the dip di- rection of $L_{Down}$ is not truly parallel with the bed, so the angular separation shown in this schematic between $L_{Down}$ and the bed should be considered as a maximum. The dip direction of $L_{Up}$ is almost bed-parallel. Gray lines illustrate schematically the inferred relationship to other structures visible in the wall of the gully chute. Dotted lines represent the potential configurations of the structures at depth, either as quasi-linear connections between the bed and the surface, or as up-warped layers derived from depth.
<figure><img src="image_9.png" /></figure>

**Fig. 7.** Manifestations of the southern field of quasi-linear flow-transverse ridges and troughs where they are intersected by the gully chute. At least three of these structures (black arrows) can be traced to shallow depths within the gully wall. The two eastern (i.e. upflow) structures manifest as low-relief ridges on the VFF surface, while the westernmost (downflow) structure manifests as a fracture-like trough. Orthorectified HiRISE image ESP_051036_1370. HiRISE data credit: NASA/JPL/University of Arizona/NHM.

<figure><img src="image_10.png" /></figure>

**Fig. 8.** Examples of inclined, debris-rich structures formed by ice-flow compression near to the termini of glaciers on Earth. (A) Debris-rich thrusts propagating up to the ice surface from the basal debris layer of Midtre Lovénbreen, Svalbard (image from 1999, person for scale). Photograph by Michael Hambrey, www.glaciers online.net (last accessed 19 January 2023). (B) Aerial view of debris-rich arcuate bands near the terminus of a glacier in Phillips Inlet in NW Ellesmere Island in the Canadian High-Arctic. The arrow shows the approximate perspective of the ground-based image in panel C. (C) Steeply-dipping bands of debris derived from the bed of the glacier in panel B, and connecting to the arcuate bands visible on the surface. Panels B and C are modified from Evans (2009).
in Fig. 2B and 3) connected to the curvilinear structures on the VFF surface bear a striking resemblance to inclined layers of ice and debris formed under ice-flow compression near to the termini of glaciers on Earth (Fig. 8; Hambrey et al., 2005; Evans, 2009; Jennings and Hambrey, 2021). Flow compression as glaciers thin and slow towards their termini forces ice flow up towards the surface, and can tilt internal layering. This occurs within glaciers of all scales, from small valley glaciers to the margins of sheets (e.g., Reeh et al., 1991; Hambrey et al., 2005). It can cause layers to outcrop in age sequences, with older layers outcropping nearer to the ice margin. Flow compression can also induce thrusting within glacial ice, which can form inclined debris bands that dissect the ice and transport lithic material from glacier beds to their surfaces. On Earth, thrusting commonly occurs under polythermal regimes, where water-lubricated warm-based ice converges with slower, frozen, cold-based ice near to glacier margins (e.g., Hambrey et al., 2005). Fig. 8 shows examples of inclined debris thrusts near the margins of High Arctic glaciers that result from ice flow compression approaching their termini.

The downflow transition in VFF-internal layer inclinations that we observe (from approximately bed-parallel layers to inclined structures that are angled upwards from the bed to the surface) supports the formation of $L_{Down}$ and associated subparallel structures approximately in situ, via a dynamic mechanism, such as compressional ice flow deformation. We infer that $L_{Down}$ and the associated subparallel structures have a 'nested spoon' morphology in the subsurface, as is common in glaciers on Earth (e.g., Figs. 18 and 69b in Jennings and Hambrey, 2021). This is consistent with (a) the lobate nature of the curvilinear undulations, and their re-orientation from ice-flow parallel to ice-margin parallel approaching the VFF terminus; and (b) the North East dip direction of $L_{Down}$ where it outcrops in the gully wall, to the south of the apex of the curvilinear surface structure (Fig. 3B). At the apex of the curvilinear surface structure, we expect that the dip direction is towards the East (i.e., up-glacier), at approximately 180° to the dip direction of the bed.

We interpret the approximately bed-conformal nature of $L_{Up}$ as reflective of the VFF's primary stratification. The particularly rough signature of $L_{Up}$ suggests it could have a higher/coarser lithic content than both the bulk VFF and the other VFF-internal structures, but this interpretation is uncertain. Based on its configuration, we interpret $L_{Up}$ as having been emplaced at the surface (i.e., potentially hosting supraglacial debris including hillslope material) before subsequent burial by overlying ice and transport down-flow. The associated subparallel layers visible above $L_{Up}$ within the gully wall do not have such strong color or roughness signatures, nor do they connect to clear structures on the VFF surface. Hence, they may not be similarly debris-rich.

Our interpretation of $L_{Down}$ and the associated inclined structures is that they formed approximately in situ by ice flow compression, and host material transported from the deep interior and/or bed of the VFF. They could represent either layers of ice and intervening debris/dust bands up-warped from depth, or debris-rich thrusts. If they are up-warped layers, they could represent an age sequence outcropping at the surface that includes old layers sourced from depth; in this configuration, layer ages (of both ice and lithic layers) likely increase along a flow-parallel transect of the surface approaching the VFF margin. If the structures are thrusts, they could form a direct connection (and hence debris transport pathway) from the VFF bed to its surface.

Our in situ flow compression hypothesis contrasts with a passive advection model proposed by Levy et al. (2021) to explain boulder-rich bands associated with arcuate structures on the surfaces of some larger (i.e., tens of kilometer-scale) VFFs on Mars. By analogy to debris-covered glaciers in the Antarctic Dry Valleys (Mackay et al., 2014), Levy et al. (2021) propose that the boulder bands represent outcrops of tilted layers of supraglacial debris which formed in hollows on the ice surface near to VFF headwalls, were buried by later ice accumulation, then passively advected down-flow. The Levy et al. (2021) model requires multiple long-term cycles of glacier accumulation and ablation induced by

orbitally-driven climate change. In contrast, our proposed compression mechanism can occur on shorter timescales, within a single glacial cycle. It is, to a large degree, controlled by ice geometry and bed topography and does not require climate dynamism. We do not refute that the Levy et al. (2021) mechanism could have occurred on Mars. Rather, we present an additional mechanism, which seems most consistent with our observations of the incised VFF in Nereidum Montes, and should be considered among potential mechanisms for the formation of inclined VFF-internal layers on Mars. Compression-driven deformation occurs at a range of scales on Earth, from small valley glaciers to large ice sheets (e.g., Reeh et al., 1991; Hambrey et al., 2005; Jennings and Hambrey, 2021), and is likely to have occurred in both small (kilometer-scale) and larger (tens of kilometer-scale) VFFs on Mars.

We do not consider an additional mechanism explored by Petersen et al. (2020) to be a viable analogue to explain the inclined VFF-internal structures we observe (although it may have occurred on Mars in some locations). Petersen et al. (2020) observe inclined debris bands outcropping at the surface of debris-covered ice within a rock glacier on Earth. These debris bands form small, curvilinear surface outcrops with total lengths of tens of meters (Petersen et al., 2020, their Fig. 2). The small scale of these features directly relates to their origins as accumulations of point-source rockfall debris, derived from the glacier headwall, and subsequently buried by seasonal ice accumulation and advected down-flow. This contrasts with the kilometers-long curvilinear structures on the surface of the incised VFF, which bound ice sourced from multiple alcoves, and form quasi-continuous structures over the entire width of the VFF.

It seems probable that multiple mechanisms could have generated inclined layers within individual VFFs in some locations, particularly in topographically complex, high-relief regions with complex histories of ice mass balance and ice flow. However, the incised VFF at our study site is relatively simple, in both surface morphology and topographic setting, compared to many VFFs, and we find no evidence to suggest that more than one mechanism is required to explain the inclined VFF-internal structures in this location.

## 6. Modeled VFF flow dynamics and stress regime

In order to explore our hypothesis that the inclined VFF-internal structures we observe formed approximately in situ as a result of ice-flow compression, we performed 3D ice flow modeling to calculate the VFF flow dynamics and stress regime in the vicinity of the inclined VFF-internal layers.

We used the Ice Sheet and Sea level System Model (ISSM, Larour et al., 2012), a 3-dimensional finite-element higher-order ice dynamics model, to calculate flow dynamics and internal stress fields within the VFF. We used the 1 m/pixel HiRISE DTM (Section 3.1) with a mean distance between model nodes of 50 m for the VFF surface. In the absence of VFF thickness data over the VFF bed, we derived the elevation of the VFF bed by interpolating the elevation of the visible contact of the VFF with the surrounding topography around the VFF terminus, along the base of the gully wall (which appears to incise to the VFF bed), and at the headwall. We supplemented these data with a row of additional points 150 m from the eastern headwall of the VFF, in which the bed elevation is set to the surface elevation of the VFF at the point minus the inferred VFF thickness based on the depth of the gully incision. These additional inferred points were needed to preserve the break of slope between the headwall up-ice of the VFF to the east (visible just upstream of transect b-b' in Fig. 5) and the base of the gully. They also prevented negative calculated ice thickness if the bed elevation was interpolated solely from marginal elevations. Our model domain (Fig. 9A) encompassed the southern flow unit, but did not cover the entirety of the northern flow unit due to the diminishing constraints on bed topography away from the gully incision.

We initialized ISSM by calculating the isothermal steady-state velocity of the VFF using the mean annual present-day surface temperature
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** (A) Reconstructed ice thickness for the southern flow unit of the VFF. Black outline shows the area shown in Figs. 9B and C. (B) Modeled surface velocity. (C) Modeled xx component of deviatoric stress. Black line a-a' in panels B and C shows the location of the vertical flow transect shown in Fig. 9D. (D) Vertical transect of overall velocity (including x, y, and z components), with black feathers (of standard length) showing the orientation of the x, z component.</figcaption></figure>

at ~40°S of 210 K (from Mars Climate Database; Forget et al., 1999; Millour et al., 2018), assuming a frozen bed. ISSM uses Glen's flow law (Glen, 1955), for which we prescribe a stress exponent of $n = 3$, as is typical used for glaciers on Earth (Cuffey and Paterson, 2010). The potential effects of impurities and ice grain size upon the value of $n$ have been explored by various workers, including for Mars, with resulting values typically varying between $n = 2$ (at low shear stresses) and $n = 4$ (at high shear stresses, where grain size does not act as a control; e.g., Cuffey and Paterson, 2010; Parsons et al., 2011). We therefore retain $n = 3$ as a pragmatic mean of these values, noting that calculated shear stresses for an ice thickness similar to the incised VFF are in the higher range explored by Milliken et al. (2003). Given the calculated isothermal velocity as a model input, we then ran ISSM in a steady-state, thermomechanically-coupled mode to derive the internal temperature of the VFF, and hence the possible impact of shear heating on ice velocity (see e.g., Butcher et al., 2017) to give the thermomechanically-coupled ice velocity and stress regime.

Given the absence of known bed elevations for the interior of the VFF, the interpolated bed topography is almost planar, with a mean slope and orientation very similar to the ice-marginal bed. The reconstructed bed topography gives VFF thicknesses around 55–65 m across most of the domain (Fig. 9A). The hollow in the ice surface around the southern field of quasi-linear, flow-transverse ridges and troughs (Fig. 4D) gives an area of thinner ice (~40 m) in the upstream section of the southern flow unit, and the ice thins towards the margin quite quickly from around 500 m upstream of the VFF margin, particularly in the area adjacent to the gully. Given our focus on the possible

mechanism of formation of the curvilinear undulations, we focus on the flow dynamics in their vicinity, from the VFF margin to around 800 m up-flow (black outline in Fig. 9A). We do not analyze the stress regime in the areas containing the features interpreted as crevasses and crevasse-fill ridges because of uncertainty in bed elevation in that area; the surface hollow (Fig. 9A) may reflect an unknown hollow in the bed, which we cannot constrain and introduces significant uncertainty into the calculated ice thickness and modeled stress regime.

Maximum modeled ice flow velocities of around $2 \times 10^{-5}$ ma$^{-1}$ occur in the northern part of the region of interest, in the region of the thickest ice (Fig. 9B). The velocity decreases rapidly towards the western margin as the ice thins sharply. The effect of the gully on the velocity is restricted to the wall of the gully itself, where velocity reduces towards the south as the ice thins into the gully incision. The gully has a minimal effect on the velocity field elsewhere; ice flow follows the local surface slope, with down-slope component of velocity dominating along transect a-a' in Fig. 9B–D. We therefore deemed it unnecessary to artificially infill the gully for the model reconstruction. The improvement of model representation in the area of interest would be minimal, while reconstruction of the former ice surface across the gully (which has significant topographic asymmetry e.g., transect c-c' in Fig. 5) would introduce additional uncertainties.

The xx component of deviatoric stress (i.e., the stress acting normal to a plane oriented transverse to ice flow) at the surface of the VFF (Fig. 9C) demonstrates that the curvilinear surface features correspond very closely to a strong band of negative (i.e., compressive) stress values where the velocity decreases rapidly towards the VFF terminus.

A vertical transect through the area of interest (Fig. 9D, location shown by black line a–a' in Figs. 9B and C) shows the velocity has a strong depth-dependence in the upper half of the transect (from around 450 m upflow of the ice margin), and flow is oriented parallel to the bed and surface. In the lower section, as the surface velocity decreases rapidly, depth dependence decreases sharply. This causes the orientations of the 3-dimensional flow vectors (short black lines in Fig. 9D) to change to an increasingly bed-divergent orientation, with strong emergent flow occurring at the surface. This divergence peaks around 300–350 m from the terminus of the VFF, in the region of $L_{Down}$. The decrease in surface-emergent flow upstream of this area also corresponds with the decrease in dip between $L_{Down}$ and $L_{Up}$.

For our investigation, the velocity and stress magnitudes modeled in Fig. 9 are secondary in importance to the overall velocity and stress patterns, which are controlled by the geometry of the ice surface and bed topography and show a high degree of correlation to the curvilinear VFF-surface structures. The absolute values of velocity and deviatoric stress are controlled by the mean annual surface temperature, which is set to the present day value. Variations in simulated mean annual surface temperature in our steady-state model change the magnitudes but not the patterns of the velocity and stress fields. The modeled velocities increase by a factor of five with a mean annual temperature of 220 K, and by a factor of 20 with a mean annual temperature of 230 K. The simulated velocities and stresses for contemporary mean annual surface temperature conditions are conservative; mean annual temperature, and hence velocity and stress maxima are likely to have been higher multiple times since the ice was deposited, due to numerous cycles of planetary spin-axis obliquity around the present-day value (Laskar et al., 2004). Given contemporary mean annual surface temperature, velocities are extremely slow, though significant deformation is still possible given the passage of time. The simulated ice flow velocity in the vicinity of the inclined VFF-internal structures approaches $\sim 10 \times 10^{-6}$ ma$^{-1}$. At this velocity, it would take a particle $\sim 4$ Myr to be transported a distance of 40 m (the approximate thickness of the VFF near to the inclined structures). A factor of five increase in velocity given a mean annual surface temperature of 220 K reduces this timescale to $\sim 800$ kyr. This is consistent with the modeled surface age estimate of $2.7+/– 0.8$ Myr for VFFs in Nereidum Montes (Berman et al., 2021), particularly considering that this represents an absolute minimum age of the surface (rather
than the bulk VFF) and is based on craters accumulated since the last period of crater destruction (e.g., by flow or resurfacing). The compression-driven formation of the inclined structures could therefore have occurred within reasonable timescales, even under conservative temperature scenarios.

The observable flow lengths of VFFs (e.g., Souness et al., 2012), and deformation structures on their surfaces, clearly imply displacements over much larger distances (kilometers to tens of kilometers) than the thickness of the incised VFF. This necessitates higher velocities in the past. Indeed, an alternative explanation that the structures were advected from an accumulation zone near to the headwall (following the Levy et al., 2021 model) requires far greater displacements relative to vertical transport through the ice column.

Our modeling suggests that compressive stresses, caused by the sharp drop in velocity as the ice thins towards the margin, could have led to a change in flow orientation from bed/surface-parallel to bed-divergent/ surface-emergent under a cold-based (i.e. frozen bed) regime in the zone of the VFF hosting the inclined structures. They could represent either up-warping of layers of ice and dust/debris from the VFF's deep interior, or ice dissected by debris-rich thrusts sourced from the bed. Inclined, debris-rich thrust structures typically occur at the margins of polythermal glaciers on Earth, where sliding warm-based ice in the core of a glacier converges with thinner, cold-based ice approaching the terminus (Fig. 8; e.g., Hambrey et al., 2005). Previous studies suggest that some VFFs on Mars, of similar scales (and therefore thicknesses) to the gully-incised VFF analyzed here, may have had polythermal regimes in the past (Hubbard et al., 2011; Conway et al., 2018). A polythermal regime, with a meltwater-lubricated core (even with small amounts of meltwater, in microscopic thin films) flowing towards cold-based mar- ginal ice, would likely enhance ice flow velocities, and the magnitude of ice-flow compression in the vicinity of the inclined VFF-internal struc- tures. However, we have limited observations independent of the in- clined structures with which to constrain past variations in the thermal regime of this VFF from an assumed cold-based present-day state. While a polythermal regime is not necessarily a pre-requisite for the debris-rich thrust mechanism, we consider up-warping of deep layers of ice and dust/debris to be the simpler explanation for the observed structures. Nevertheless, a debris-rich thrust origin should not be ruled out.

# 7. Implications

## 7.1. Implications for the internal structure of VFFs

VFFs on Mars host diverse surface structures which likely record complex ice flow histories (e.g., Dickson et al., 2008, 2012; Souness and Hubbard, 2013; Brough et al., 2015; Levy et al., 2021). Our observations demonstrate that some structures on VFF surfaces are the surface out- crops of inclined layers, which connect to deep VFF interiors and potentially all the way to their beds. Based on observations of concen- trated boulder bands on the surfaces of larger (i.e., tens of kilometer- scale) VFFs on Mars, Levy et al. (2021) suggested that VFFs could contain fine-scale tilted debris layers that have not yet been detected by orbital ground-penetrating radar observations from SHARAD; our ob- servations support this inference for a kilometer-scale VFF. While Levy et al. (2021) inferred the presence of inclined subsurface structures on the basis of surface bands of meter-scale boulders, which we also observe for this site, we do not have a strong constraint on the full grain size distribution within such structures, particularly for smaller grain sizes. We also emphasise that we do not have direct observations of lithics within the structures in the incised VFF, rather we infer higher lithic content for $L_{Down}$ and associated subparallel structures, and $L_{Up}$, on the basis of color and roughness signatures. Thus, the true ice-lithic ratio of the materials comprising these structures remains uncertain.

As summarized in an extensive review by Jennings and Hambrey (2021), there are many and various mechanisms by which glaciers and ice sheets on Earth can evolve complex surface and interior structures.

Our observations and modeling suggests that the inclined VFF-internal structures we observe could have formed approximately in situ as a result of ice flow compression approaching the VFF terminus. We do not refute that passive-advection of buried, tilted supraglacial debris layers formed over multiple climate cycles (as proposed by Levy et al., 2021) could also have occurred in some VFFs; both mechanisms could have operated on Mars, at different times, in different locations, or indeed within a single VFF system. The formation of small-scale inclined layers sourced from localised rockfalls on VFF headwalls, as explored on Earth by Petersen et al. (2020), could also have occurred on Mars, but is inconsistent with the scale and configuration of the curvilinear struc- tures observed for the incised VFF examined here.

## 7.2. Potential implications for glacial landform development

On Earth, inclined, debris-rich englacial layers formed by ice-flow compression can evolve into ice-cored moraines called controlled mo- raines when their bounding ice is removed (e.g., Evans, 2009). Debris released from englacial layers covers and protects underlying ice layers from mass loss. On Earth, controlled moraines have low preservation potential due to melt-out and thermal erosion of the ice core (Evans, 2009). However, ice-cored moraines may not experience the same fate under environmental conditions on Mars. Debris-rich internal structures within VFFs could have supplied bands of debris to VFF surfaces, which subsequently evolved into ice-cored moraine ridges when climate con- ditions changed from favoring ice accumulation to favoring sublimation. Moraine-like arcuate ridges are common in association with viscous flow features on Mars (Milliken et al., 2003; Arfstrom and Hartmann, 2005). Indeed, examples of such features are present near to the incised VFF, likely associated with ice flow originating from the hillslope to the south (Fig. 1). Various hypotheses have been proposed for the formation of moraine-like arcuate ridges on Mars, including supraglacial transport of hillslope debris to the ice margin, proglacial bulldozing, and proglacial glaciotectonic deformation (e.g., Arfstrom and Hartmann, 2005; Jawin et al., 2018; Conway et al., 2018). Hubbard et al. (2011) proposed that rectilinear 'moraine mounds' arranged in arcuate ridges bounding a viscous flow feature in Greg crater could have originated (among various potential mechanisms) as thrusts of basal debris transported up from the glacier bed. Our observations suggest that controlled moraines comprising material sourced from englacial debris layers (be they thrusts or up-warped layering) should be more extensively considered among potential origins for arcuate ridges associated with VFFs on Mars, particularly where they occur in zones of past ice-flow compression.

## 7.3. Implications for future scientific sampling of ice deposits on Mars

Mars' mid-latitude ice deposits are key targets for future orbital, robotic, and human exploration of Mars (Mars Ice Core Working Group, 2021; European Space Agency, 2022; National Academies of Science, Engineering, and Medicine, 2022; I-MIM MDT, 2022). Among the numerous scientific opportunities raised by sampling in these environ- ments are the extraction of paleoclimate records from ice samples, and the search for past and/or present microbial life (e.g., Mars Ice Core Working Group, 2021). Glacial environments on Earth are important microbial habitats (Hodson et al., 2008), and Mars' ice deposits are therefore important targets for astrobiological investigations. Subglacial and deep-englacial environments are of particular astrobiological in- terest because of: (a) the availability of water ice, which has potential to generate microscopic interfacial water films at subfreezing temperatures (including at subfreezing temperatures in glaciers that are nominally cold-based; e.g., Cuffey et al., 2000; Dash et al., 1995); and (b) their protection from seasonal/diurnal temperature fluctuations and ionizing radiation at the surface (e.g., Skidmore et al., 2000; Montross et al., 2014). It is unknown whether micro-scale interfacial meltwater film formation has occurred at ice/lithic grain interfaces within VFFs on Mars; however, it is a distinct possibility, particularly considering the
abundance of salts (which can significantly depress the melting point by up to ~70 K (e.g., Hecht et al., 2009; Clark and Kounaves, 2016)). Lithics enhance the habitability of ice, providing energy sources, and encouraging interfacial meltwater film formation (e.g., Montross et al., 2014; Telling et al., 2015; Onstott et al., 2019). Thus, while lithics could pose hazards for drilling/coring on Mars, sampling near to them could bring benefits for the search for microbial life.

While subglacial/deep-englacial environments on Mars are of great scientific interest, accessing them for scientific sampling, at depths of tens to hundreds of meters beneath the present-day surface, would require high cost, high risk drilling/coring. Our observations suggest that ice flow could have transported sub/englacial materials (lithics and/or older ice) up to more accessible positions at, or near to, the surfaces of VFFs on Mars in some locations. This could allow lower cost, lower-risk sampling of ice and/or lithics transported from subglacial and/or deep englacial environments to the shallow subsurface/surface.

Our observations also demonstrate that VFFs do not necessarily comprise simple sub-horizontal layers throughout, and that in some locations (and particularly in zones of ice flow compression), the age sequence of layers may deviate from vertical. On Earth, traditional ice coring accesses palaeoenvironmental records by coring vertically through horizontally-layered ice. However, deep coring is technologically challenging even on Earth. Our observations suggest that future ice coring missions may be able to exploit inclined VFF layering to reduce the requirement for deep coring. Where inclined layering is formed by up-warping of deep layers under flow compression, layer age sequences could be accessed via near-surface sampling along a transect of a VFF surface (with layer age increasing towards the margin). This ‘horizontal ice core’ approach has been demonstrated on Earth, with surface ice samples extracted from the margin of the Greenland Ice Sheets providing palaeoenvironmental records that are in good agreement with those extracted from deep cores in the ice sheet interior (Reeh et al., 1991; Petrenko et al., 2006). Planning for the extraction of paleo-environmental records from ice samples should therefore include detailed characterization of VFF surface structures and their relationships to ice-internal layer configurations and ice flow histories. A better understanding of these characteristics would inform decisions about the balance of drilling/coring hazards and scientific opportunities of sampling at a given site.

### 7.4. Implications for future ice access and ISRU

Evidence for high ice purity (up to ~90% by volume for lobate debris apron-type VFFs; Holt et al., 2008; Plaut et al., 2009; Petersen et al., 2018; Gallagher et al., 2021) makes VFFs attractive targets for ISRU. However, their internal structure, and particularly its manifestation in the shallow subsurface remains a key knowledge gap that will be critical for planning for ice access and ISRU efforts. As predicted by Levy et al. (2021), our observations demonstrate that fine-scale inclined layers, potentially with relatively high concentrations of lithics, do occur within VFFs on Mars. The grain size distributions, and ice-lithic ratios within these layers (which are likely to vary spatially and according to the specific origins and source regions of a given layer) are poorly constrained, and could range from low concentrations of fine-grained materials (e.g., dust, sand) to concentrated bands of meter-scale boulders, or a combination of both. Thus, while these layers may provide opportunities for scientific analyses (Section 7.4), they should also be considered as potential hazards for drilling and ISRU activities on VFFs.

### 7.5. Future orbital remote sensing priorities for structural glaciology on Mars

The SHARAD instrument on MRO has greatly enhanced our understanding of VFFs on Mars, particularly their bulk compositions as debris-covered water ice glaciers (e.g., Holt et al., 2008; Plaut et al., 2009; Petersen et al., 2018; Gallagher et al., 2021). High resolution (decimeter

and meter-scale) imaging, such as that provided by the HiRISE and CTX cameras on MRO, continues to provide important new insights on the surface and subsurface properties of Mars' ice deposits (Dundas et al., 2021; Levy et al., 2021; this study). More recent, spatially extensive color imaging (at meter scales) by the CaSSIS camera (Thomas et al., 2017) on ESA ExoMars Trace Gas Orbiter provides additional opportunities to better understand surface structures on Mars' ice deposits (e.g., Grau Galofre et al., 2022). Meter-to-decimeter-scale elevation models derived from stereo-pair imaging are also highly valuable for understanding the three-dimensional manifestations of structures visible at the surface, and as inputs into ice flow models.

Based on the observations presented here, we consider it to be a high priority that near-term orbital missions preparing for future robotic and human missions to Mars' ice deposits better constrain the internal structure (i.e., 'structural glaciology') of VFFs, both in the shallow subsurface (upper ~10 m) and in their deep interiors (10s to 100s meters depth). We recommend a combination of high-resolution (decimeter to meter-scale) full-color stereo imaging, and geophysical investigations capable of resolving both the shallow and deep interior structure of VFFs, including inclined interior structures. Meter-to-decimeter scale elevation information, targeted to complement geophysical observations, would permit analyses of the three-dimensional manifestations of subsurface structures outcropping at the surface, and serve as valuable input to ice flow models. Spatially extensive, high-resolution (meter-scale) compositional information would also provide important opportunities to analyze potentially heterogeneous spatial distributions, compositions, sources, and alteration histories of materials outcropping at VFF surfaces.

## 8. Conclusions

Until now, observations of the internal structure of viscous flow features (VFFs; thought to be debris-covered glaciers) in Mars' mid latitudes have remained elusive. We present a view of the internal structure of a small VFF on Mars, as revealed by a natural gully incision through the VFF, from its headwall to its terminus. The gully incision reveals the subsurface configurations of VFF-internal layers and crevasses. Our observations suggest that curvilinear undulations on the VFF surface are outcrops of VFF-internal structures. Near to the VFF margin, the VFF-internal structures ($L_{Down}$ and associated subparallel structures) dip down into the deep VFF interior, and potentially all the way to the bed. We performed 3D ice flow modeling, which demonstrates that these structures correlate very closely with a zone of ice flow compression and divergence of ice flow away from the bed. This suggests the inclined structures could have become inclined in situ under ice flow compression. The structures might be akin to up-warped layers or debris-rich thrusts within glaciers on Earth. They could host lithic material in relatively fine-scale layers, compared to adjacent ice layers and the bulk VFF. Up-flow of the inclined structures, the VFF-internal layers ($L_{Up}$ and associated subparallel layers) appear to be approximately conformal to the bed orientation and slope and may be more reflective of horizontal layering from ice accumulation. The small divergence from the bed of $L_{Up}$ could be associated with the first onset of bed-divergent ice flow in this area. $L_{Up}$ manifests as a rough outcropping curvilinear structure on the VFF-surface, and could host supraglacial debris that was buried by subsequent ice accumulation.

Our observations of the inclined VFF-internal structures ($L_{Down}$ and associated subparallel structures), combined with our modeling (which shows a high degree of correlation with a zone of ice flow compression), suggests that the structures could have transported material (ice and/or lithics) from deep within the VFF, and potentially even its bed, up to positions at (or near to) the shallowly-buried surface of the VFF. By extension, it may be possible for future missions to sample materials from VFF-bed or deep interior environments without the requirement for deep drilling, including age sequences represented by outcrops of up-warped layers. Planning for ice access and ISRU of water ice by future
missions to Mars' subsurface ice deposits should consider the potential for spatially variable subsurface layer configurations, including the potential for inclined layers of lithic material near to VFF termini and in other zones of past/present ice flow compression.

## Declaration of Competing Interest

The authors have no competing interests to declare.

## Data availability

Orbital datasets used herein are freely available from the NASA Planetary Data System at https://ode.rsl.wustl.edu/mars/. CTX images are available from the University of Arizona at http://viewer.mars.asu.edu/viewer/ctx, and HiRISE products are available from the HiRISE website at https://www.uahirise.org/. The HiRISE DTM and orthoimage we generated for this study are available at https://doi.org/10.15131/shef.data.23741973.v1. The Ice Sheet and Sea Level System Model (ISSM) used herein is available at: https://issm.jpl.nasa.gov/.

## Acknowledgements

This manuscript was prepared while FEGB was in receipt of a Leverhulme Trust Early Career Fellowship. FEGB also acknowledges financial support from the Région Pays de la Loire GeoloNet project, which facilitated an in-person collaboration visit with SJC in May 2019 (convention No 2016-10982). SJC is funded for her HiRISE work by the French Space Agency CNES. JMD gratefully acknowledges UK Space Agency support (ST/R002355/1; ST/V002678/1; ST/W002566/1; ST/W002566/2). We thank Stefano Nerozzi and Joe Levy for thoughtful reviews which improved the manuscript. We also thank Isaac Smith for insightful discussions on SHARAD radar, and Bryn Hubbard for discussions on structural glaciology. We gratefully acknowledge the hard work of the various teams who brought Mars Reconnaissance Orbiter to fruition, and who continue to deliver fascinating data more than one and a half decades after MRO's arrival at Mars.

## References

Arfstrom, J., Hartmann, W.K., 2005. Martian flow features, moraine-like ridges, and gullies: terrestrial analogs and interrelationships. Icarus 174, 321–335. https://doi.org/10.1016/j.icarus.2004.05.026.

Baker, D.M.H., Head, J.W., 2015. Extensive middle Amazonian mantling of debris aprons and plains in Deuteronilus Mensae. Mars: implications for the record of mid-latitude glaciation. Icarus 260, 269–288. https://doi.org/10.1016/j.icarus.2015.06.036.

Bennett, M.R., Huddart, D., Waller, R.I., 2000. Glaciofluvial crevasse and conduit fills as indicators of supraglacial dewatering during a surge. Skeiðarárjökull, Iceland. J. Glaciol. 46, 25–34. https://doi.org/10.3189/172756500781833232.

Berman, D.C., Chuang, F.C., Smith, I.B., Crown, D.A., 2021. Ice-rich landforms of the southern mid-latitudes of Mars: a case study in Nereidum Montes. Icarus 355, 114170. https://doi.org/10.1016/j.icarus.2020.114170.

Brough, S., Hubbard, B., Souness, C., Grindrod, P.M., Davis, J., 2015. Landscapes of polyphase glaciation: eastern Hellas Planitia, Mars. J. Maps 1–13. https://doi.org/10.1080/17445647.2015.1047907.

Butcher, F.E.G., Balme, M.R., Gallagher, C., Arnold, N.S., Conway, S.J., Hagermann, A., Lewis, S.R., 2017. Recent basal melting of a mid-latitude glacier on Mars. J. Geophys. Res. Planets 122, 2445–2468. https://doi.org/10.1002/2017JE005434.

Clark, B.C., Kounaves, S.P., 2016. Evidence for the distribution of perchlorates on Mars. Int. J. Astrobiol. 15, 311–318. https://doi.org/10.1017/S1473550415000385.

Conway, S.J., Balme, M.R., 2014. Decameter thick remnant glacial ice deposits on Mars. Geophys. Res. Lett. 41, 5402–5409. https://doi.org/10.1002/2014GL060314.

Conway, S.J., Butcher, F.E.G., de Haas, T., Deijns, A.A.J., Grindrod, P.M., Davis, J.M., 2018. Glacial and gully erosion on Mars: a terrestrial perspective. Geomorphology 318, 26–57. https://doi.org/10.1016/j geomorph.2018.05.019.

Cuffey, K.M., Paterson, W.S.B., 2010. The Physics of Glaciers, 4th ed. Elsevier Science, Oxford.

Cuffey, K.M., Conway, H., Gades, A.M., Hallet, B., Lorrain, R., Severinghaus, J.P., Steig, E.J., Vaughn, B., White, J.W.C., 2000. Entrainment at cold glacier beds. Geology 28, 351–354. https://doi.org/10.1130/0091-7613(2000)28<351: EACGB>2.0.CO;2.

Dash, J.G., Fu, H., Wettlaufer, J.S., 1995. The premelting of ice and its environmental consequences. Rep. Prog. Phys. 58, 115. https://doi.org/10.1088/0034-4885/58/1/003.

Dickson, J.L., Head, J.W., Marchant, D.R., 2008. Late Amazonian glaciation at the dichotomy boundary on Mars: evidence for glacial thickness maxima and multiple glacial phases. Geology 36, 411–414. https://doi.org/10.1130/G24382A.1.

Dickson, J.L., Head, J.W., Fassett, C.I., 2012. Patterns of accumulation and flow of ice in the mid-latitudes of Mars during the Amazonian. Icarus 219, 723–732. https://doi.org/10.1016/j.icarus.2012.03.010.

Dundas, C.M., Mellon, M.T., Conway, S.J., Daubar, I.J., Williams, K.E., Ojha, L., Wray, J., J., Bramson, A.M., Byrne, S., McEwen, A.S., Posiolova, L.V., Speth, G., Viola, D., Landis, M.E., Morgan, G.A., Pathare, A.V., 2021. Widespread exposures of extensive clean shallow ice in the Midlatitudes of Mars. J. Geophys. Res. Planets 126, e2020JE006617. https://doi.org/10.1029/2020JE006617.

European Space Agency, 2022. Terrae Novae 2030+ Strategy Roadmap. URL: https://esamultimedia.esa.int/docs/HRE/Terrae_Novae_2030+strategy-roadmap.pdf, last accessed 27/02/2023.

Evans, D.J.A., 2009. Controlled moraines: origins, characteristics and palaeoglaciological implications. Quat. Sci. Rev. 28, 183–208. https://doi.org/10.1016/j.qusacrev.2008.10.024.

Forget, F., Hourdin, F., Fournier, R., Hourdin, C., Talagrand, O., Collins, M., Lewis, S.R., Read, P.L., Huot, J.-P., 1999. Improved general circulation models of the Martian atmosphere from the surface to above 80 km. J. Geophys. Res. 104, 24155–24175. https://doi.org/10.1029/1999JE001025.

Gallagher, C., Butcher, F.E.G., Balme, M., Smith, I., Arnold, N., 2021. Landforms indicative of regional warm based glaciation, Phlegra Montes. Mars. Icarus 355, 114173. https://doi.org/10.1016/j.icarus.2020.114173.

Glen, J.W., 1955. The creep of polycrystalline ice. Proc. Royal Soc. London. Series A. Math. Phys. Sci. 228, 519–538. https://doi.org/10.1098/rspa.1955.0066.

Grau Galofre, A., Serla, J.K., Becerra, P., Noblet, A., Conway, S.J., 2022. Patterns of martian glacial deformation: implications for glacio-geology, internal structure, and regional climate. Planet. Space Sci. 221, 105548 https://doi.org/10.1016/j.pss.2022.105548.

Hambrey, M.J., Murray, T., Glasser, N.F., Hubbard, A., Hubbard, B., Stuart, G., Hansen, S., Kohler, J., 2005. Structure and changing dynamics of a polythermal valley glacier on a centennial timescale: Midre Lovénbreen, Svalbard. J. Geophys. Res. 110, F01006. https://doi.org/10.1029/2004JF000128.

Head, J.W., Mustard, J.F., Kreslavsky, M.A., Milliken, R.E., Marchant, D.R., 2003. Recent ice ages on Mars. Nature 426, 797–802. https://doi.org/10.1038/nature02114.

Head, J.W., Neukum, G., Jaumann, R., Hiesinger, H., Hauber, E., Carr, M., Masson, P., Foing, B., Hoffmann, H., Kreslavsky, M., Werner, S., Milkovich, S., van Gasselt, S., The HRSC Co-Investigator Team, 2005. Tropical to mid-latitude snow and ice accumulation, flow and glaciation on Mars. Nature 434, 346–351. https://doi.org/10.1038/nature03359.

Hecht, M.H., Kounaves, S.P., Quinn, R.C., West, S.J., Young, S.M.M., Ming, D.W., Catling, D.C., Clark, B.C., Boynton, W.V., Hoffman, J., DeFlores, L.P., Gospodinova, K., Kapit, J., Smith, P.H., 2009. Detection of perchlorate and the soluble chemistry of martian soil at the phoenix lander site. Science 325, 64–67. https://doi.org/10.1126/science.1172466.

Hepburn, A.J., Ng, F.S.L., Livingstone, S.J., Holt, T.O., Hubbard, B., 2020. Polyphase mid-latitude glaciation on Mars: chronology of the formation of superposed glacier-like forms from crater-count dating. J. Geophys. Res. Planets 125, e2019JE006102. https://doi.org/10.1029/2019JE006102.

Hodson, A., Anesio, A.M., Tranter, M., Fountain, A., Osborn, M., Priscu, J., Laybourn-Parry, J., Sattler, B., 2008. Glacial Ecosystems. Ecol. Monogr. 78, 41–67. https://doi.org/10.1890/07-0187.1.

Hoffman, S.J., Andrews, A., Joosten, B.K., Watts, K., 2017. A water rich mars surface mission scenario. In: 2017 IEEE Aerospace Conference, pp. 1–21. https://doi.org/10.1109/AERO.2017.7943911.

Holt, J.W., Safaeinili, A., Plaut, J.J., Head, J.W., Phillips, R.J., Seu, R., Kempf, S.D., Choudhary, P., Young, D.A., Putzig, N.E., Biccari, D., Gim, Y., 2008. Radar sounding evidence for buried glaciers in the southern mid-latitudes of Mars. Science 322, 1235–1238. https://doi.org/10.1126/science.1164246.

Hubbard, B., Milliken, R.E., Kargel, J.S., Limaye, A., Souness, C., 2011. Geomorphological characterisation and interpretation of a mid-latitude glacier-like form: Hellas Planitia, Mars. Icarus 211, 330–346. https://doi.org/10.1016/j.icarus.2010.10.021.

I-MIM MDT, 2022. Final Report of the International Mars Ice Mapper Reconnaissance/Science Measurement Definition Team, p. 239.

Jawin, E.R., Head, J.W., Marchant, D.R., 2018. Transient post-glacial processes on Mars: geomorphologic evidence for a paraglacial period. Icarus 309, 187–206. https://doi.org/10.1016/j.icarus.2018.01.026.

Jennings, S.J.A., Hambrey, M.J., 2021. Structures and deformation in glaciers and ice sheets. Rev. Geophys. 59, e2021RG000743 https://doi.org/10.1029/2021RG000743.

Khuller, A.R., Christensen, P.R., 2021. Evidence of exposed dusty water ice within martian gullies. J. Geophys. Res. Planets 126, e2020JE006539. https://doi.org/10.1029/2020JE006539.

Kirk, R.L., Howington-Kraus, E., Rosiek, M.R., Anderson, J.A., Archinal, B.A., Becker, K. J., Cook, D.A., Galuszka, D.M., Geissler, P.E., Hare, T.M., Holmberg, I.M., Keszthelyi, L.P., Redding, B.L., Delamere, W.A., Gallagher, D., Chapel, J.D., Eliason, E.M., King, R., McEwen, A.S., 2008. Ultrahigh resolution topographic mapping of Mars with MRO HiRISE stereo images: meter-scale slopes of candidate Phoenix landing sites. J. Geophys. Res. 113, E00A24. https://doi.org/10.1029/2007JE003000.

Larour, E., Seroussi, H., Morlighem, M., Rignot, E., 2012. Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM). J. Geophys. Res. 117, F01022. https://doi.org/10.1029/2011JF02140.
Laskar, J., Correia, A.C.M., Gastineau, M., Joutel, F., Levrard, B., Robutel, P., 2004. Long term evolution and chaotic diffusion of the insolation quantities of Mars. *Icarus* 170, 343–364. https://doi.org/10.1016/j.icarus.2004.04.005.

Levy, J.S., Fassett, C.I., Holt, J.W., Parsons, R., Cipolli, W., Goudge, T.A., Tebolt, M., Kuentz, L., Johnson, J., Ishraque, F., Cvijanovich, B., Armstrong, I., 2021. Surface boulder banding indicates Martian debris-covered glaciers formed over multiple glaciations. *PNAS* 118. https://doi.org/10.1073/pnas.2015971118.

Mackay, S.L., Marchant, D.R., Lamp, J.L., Head, J.W., 2014. Cold-based debris-covered glaciers: evaluating their potential as climate archives through studies of ground-penetrating radar and surface morphology. *J. Geophys. Res. Earth Surf.* 119, 2505–2540. https://doi.org/10.1002/2014JF003178.

Madeleine, J.-B., Forget, F., Head, J.W., Levrard, B., Montmessin, F., Millour, E., 2009. Amazonian northern mid-latitude glaciation on Mars: a proposed climate scenario. *Icarus* 203, 390–405. https://doi.org/10.1016/j.icarus.2009.04.037.

Malin, M.C., Bell, J.F., Cantor, B.A., Caplinger, M.A., Calvin, W.M., Clancy, R.T., Edgett, K.S., Edwards, L., Haberle, R.M., James, P.B., Lee, S.W., Ravine, M.A., Thomas, P.C., Wolff, M.J., 2007. Context camera investigation on board the Mars Reconnaissance Orbiter. *J. Geophys. Res.* 112, E05S04 https://doi.org/10.1029/2006JE002808.

Mars Ice Core Working Group, 2021. First Ice Cores from Mars. co-chairs: M.R. Albert and M. Koutnik. URL: https://www.nasa.gov/sites/default/files/atoms/files/mars_ice_core Working_group_report_final_feb2021.pdf, last accessed 27/02/2023.

McEwen, A.S., Eliason, E.M., Bergstrom, J.W., Bridges, N.T., Hansen, C.J., Delamere, W. A., Grant, J.A., Gulick, V.C., Herkenhoff, K.E., Keszthelyi, L., Kirk, R.L., Mellon, M.T., Squyres, S.W., Thomas, N., Weitz, C.M., 2007. Mars reconnaissance orbiter's high resolution imaging science experiment (HiRISE). *J. Geophys. Res.* 112, E05S02. https://doi.org/10.1029/2005JE002605.

Milliken, R.E., Mustard, J.F., Goldsby, D.L., 2003. Viscous flow features on the surface of Mars: observations from high-resolution Mars orbiter camera (MOC) images. *J. Geophys. Res.-Planets* 108, 5057. https://doi.org/10.1029/2002JE002005.

Millour, E., Forget, F., Spiga, A., Vals, M., Zakharov, V., Montabone, L., Lefèvre, F., Montmessin, F., Chaufray, J.-Y., López-Valverde, M.A., González-Galindo, F., Lewis, S.R., Read, P.L., Desjean, M.-C., Cipriani, F., MCD development team, 2018. The Mars Climate Database (Version 5.3). Scientific Workshop: From Mars Express to ExoMars. https://www.cosmos.esa.int/documents/1499429/1583871/Millour_E.pdf. last accessed 24/07/2023.

Montross, S., Skidmore, M., Christner, B., Samyn, D., Tison, J.-L., Lorrain, R., Doyle, S., Fitzsimons, S., 2014. Debris-rich basal ice as a microbial habitat, Taylor glacier, Antarctica. *Geomicrobiol J.* 31, 76–81. https://doi.org/10.1080/01490451.2013.811316.

Morgan, G.A., Putzig, N.E., Perry, M.R., Sizemore, H.G., Bramson, A.M., Petersen, E.I., Bain, Z.M., Baker, D.M.H., Mastrogiuseppe, M., Hoover, R.H., Smith, I.B., Pathare, A., Dundas, C.M., Campbell, B.A., 2021. Availability of subsurface water-ice resources in the northern mid-latitudes of Mars. *Nat. Astron.* 5, 230–236. https://doi.org/10.1038/s41550-020-01290-z.

Mustard, J.F., Cooper, C.D., Rifkin, M.K., 2001. Evidence for recent climate change on Mars from the identification of youthful near-surface ground ice. *Nature* 412, 411–414. https://doi.org/10.1038/35086515.

National Academies of Science, Engineering, and Medicine, 2022. Origins, Worlds, and Life: A Decadal Strategy for Planetary Science and Astrobiology 2023–2032. The National Academies Press, Washington DC. https://doi.org/10.17226/26522.

Onstott, T.C., Ehlmann, B.L., Sapers, H., Coleman, M., Ivarsson, M., Marlow, J.J., Neubeck, A., Niles, P., 2019. Paleo-rock-hosted life on earth and the search on Mars: a review and strategy for exploration. *Astrobiology* 19, 1230–1262. https://doi.org/10.1089/ast.2018.1960.

Parsons, R.A., Nimmo, F., Miyamoto, H., 2011. Constraints on martian lobate debris apron evolution and rheology from numerical modeling of ice flow. *Icarus* 214, 246–257. https://doi.org/10.1016/j.icarus.2011.04.014.

Petersen, E.I., Holt, J.W., Levy, J.S., 2018. High ice purity of martian lobate debris aprons at the regional scale: evidence from an orbital radar sounding survey in deuteronilus and protonilus mensae. *Geophys. Res. Lett.* 45, 11,595–11,604. https://doi.org/10.1029/2018GL079759.

Petersen, E.I., Levy, J.S., Holt, J.W., Stuurman, C.M., 2020. New insights into ice accumulation at Galena Creek Rock Glacier from radar imaging of its internal structure. *J. Glaciol.* 66, 1–10. https://doi.org/10.1017/jog.2019.67.

Petrenko, V.V., Severinghaus, J.P., Brook, E.J., Reeh, N., Schaefer, H., 2006. Gas records from the West Greenland ice margin covering the last glacial termination: a horizontal ice core. *Quat. Sci. Rev.* 25, 865–875. https://doi.org/10.1016/j.qasqrev.2005.09.005.

Picardi, G., Plaut, J.J., Biccari, D., Bombaci, O., Calabrese, D., Cartacci, M., Cicchetti, A., Clifford, S.M., Edenhofer, P., Farrell, W.M., Federico, C., Frigeri, A., Gurnett, D.A., Hagfors, T., Heggy, E., Herique, A., Huff, R.L., Ivanov, B.A., Johnson, W.T.K., Jordan, R.L., Kirchner, D.L., Kofman, W., Leuschen, C.J., Nielsen, E., Orosei, R., Pettinelli, E., Phillips, R.J., Plettemeier, D., Safaeinili, A., Seu, R., Stofan, E.R., Vannaroni, G., Watters, T.R., Zampolini, E., 2005. Radar soundings of the subsurface of Mars. *Science* 310, 1925–1928. https://doi.org/10.1126/science.1122165.

Plaut, J.J., Safaeinili, A., Holt, J.W., Phillips, R.J., Head, J.W., Seu, R., Putzig, N.E., Frigeri, A., 2009. Radar evidence for ice in lobate debris aprons in the mid-northern latitudes of Mars. *Geophys. Res. Lett.* 36, L02203. https://doi.org/10.1029/2008GL036379.

Putzig, N.E., Morgan, G.A., Sizemore, H.G., Hollibaugh Baker, D.M., Petersen, E.I., Pathare, A.V., Dundas, C.M., Bramson, A.M., Courville, S.W., Perry, M.R., Nerozzi, S., Bain, Z.M., Hoover, R.H., Campbell, B.A., Mastrogiuseppe, M., Mellon, M.T., Seu, R., Smith, I.B., 2023. Ice resource mapping on Mars. In: Badescu, V., Zacny, K., Bar-Cohen, Y. (Eds.), *Handbook of Space Resources*. Springer International Publishing, Cham, pp. 583–616. https://doi.org/10.1007/978-3-030-97913-3_16.

Reeh, N., Oerter, H., Letréguilly, A., Miller, H., Hubberten, H.-W., 1991. A new, detailed ice-age oxygen-18 record from the ice-sheet margin in central West Greenland. *Glob. Planet. Chang.* 4, 373–383. https://doi.org/10.1016/0921-8181(91)90003-F.

Seu, R., Phillips, R.J., Biccari, D., Orosei, R., Masdea, A., Picardi, G., Safaeinili, A., Campbell, B.A., Plaut, J.J., Marinangeli, L., Smrekar, S.E., Nunes, D.C., 2007. SHARAD sounding radar on the Mars reconnaissance orbiter. *J. Geophys. Res. Planets* 112. https://doi.org/10.1029/2006JE002745.

Skidmore, M.L., Foght, J.M., Sharp, M.J., 2000. Microbial life beneath a high Arctic glacier. *Appl. Environ. Microbiol.* 66, 3214–3220. https://doi.org/10.1128/AEM.66.8.3214-3220.2000.

Souness, C.J., Hubbard, B., 2013. An alternative interpretation of late Amazonian ice flow: Protonilus Mensae, Mars. *Icarus* 225, 495–505. https://doi.org/10.1016/j. icarus.2013.03.030.

Souness, C., Hubbard, B., Milliken, R.E., Quincey, D., 2012. An inventory and population-scale analysis of martian glacier-like forms. *Icarus* 217, 243–255. https://doi.org/10.1016/j. icarus.2011.10.020.

Squyres, S.W., 1978. Martian fretted terrain: flow of erosional debris. *Icarus* 34, 600–613. https://doi.org/10.1016/0019-1035(78)90048-9.

Squyres, S.W., 1979. The distribution of lobate debris aprons and similar flows on Mars. *J. Geophys. Res.* 84, 8087–8096. https://doi.org/10.1029/JB084iB14p08087.

Telling, J., Boyd, E.S., Bone, N., Jones, E.L., Tranter, M., MacFarlane, J.W., Martin, P.G., Wadham, J.L., Lamarche-Gagnon, G., Skidmore, M.L., Hamilton, T.L., Hill, E., Jackson, M., Hodgson, D.A., 2015. Rock comminution as a source of hydrogen for subglacial ecosystems. *Nat. Geosci.* 8, 851–855. https://doi.org/10.1038/ngeo2533.

Thomas, N., Cremonese, G., Ziethe, R., Gerber, M., Brändli, M., Bruno, G., Erismann, M., Gambicorti, L., Gerber, T., Ghose, K., Gruber, M., Gubler, P., Mischler, H., Jost, J., Piazza, D., Pommerol, A., Rieder, M., Roloff, V., Servonet, A., Trottmann, W., Weigel, T., Viertl, J., De Roux, N., Lochmatter, P., Sutter, G., Casciello, A., Hausner, T., Ficai Veltroni, I., Da Deppo, V., Orleanski, P., Nowosielski, W., Zawistowski, T., Szalai, S., Sodor, B., Tulyakov, S., Troznai, G., Banaskiewicz, M., Bridges, J.C., Byrne, S., Debei, S., El-Maarry, M.R., Hauber, E., Hansen, C.J., Ivanov, A., Keszthelyi, L., Kirk, R., Kuzmin, R., Mangold, N., Marinangeli, L., Markiewicz, W.J., Massironi, M., McEwen, A.S., Okubo, C., Tornabene, L.L., Wajer, P., Wray, J.J., 2017. The Colour and Stereo Surface Imaging System (CaSSIS) for the ExoMars Trace Gas Orbiter. *Space Sci. Rev.* 212, 1897–1944. https://doi.org/10.1007/s11214-017-0421-1.