Landforms indicative of regional warm based glaciation, Phlegra Montes, Mars

Viscous low features (VFF) occur in the mid-latitudes of Mars and have characteristics consistent with being glaciers. Climate models suggest that martian glaciers are cold-based systems in which meltwater has never been widely produced. VFF are common in Phlegra Montes, a mountain range in the mid-latitudes of the northern hemisphere of Mars. However, in Phlegra Montes, the presence of an esker associated with an extant Amazonian Period VFF provides evidence that warm-based glacial processes did formerly operate. The problem at the centre of this paper is that the glacial meltwater responsible for this esker could have been produced as a consequence of its setting in a graben, with locally enhanced geothermal heating having been the driver of melt, not systemic heating associated with a regional warm-based regime. Given this uncertainty, this paper aims to determine if there are indicators of more widespread warm-based glacial processes in Phlegra Montes. The paper briely describes the distribution and characteristics of VFF across the region, before focussing on the search for key landforms considered diagnostic of erosion by warm-based glaciers. From our observations, including discrim- inant morphometrics, we conclude that the landscape of Phlegra Montes is indicative of widespread warm-based glacial processes, including subglacial scour, linear abrasion and the incision of subglacial meltwater channels. Our indings have signiicance in constraining the contexts and process environments within which liquid water has been produced during the Amazonian Period on Mars and point to several lines of future research into martian glaciation, climate and landscape evolution.


Introduction
Phlegra Montes are a mountain complex on Mars extending NNE-SSW, over ~1000 km, between latitudes 30 • -52 • N (Fig. 1). The region is dominated by two uplands of Hesperian-Noachian transition age (Tanaka et al., 2014; Fig. 1), composed of impact breccias, sediments and volcanic deposits. The uplands are separated by low-lying Hesperian-Amazonian lood lavas and lava lows (dated locally to 1.6 ± 0.1 Ga by Gallagher and Balme, 2015), and bordered by Hesperian sedimentary units and isolated Noachian highland remnants (Tanaka et al., 2014). The uplands are characterised by rounded, dome-like mountains, with intervening valleys and basins, some of which are occupied by viscous low features (VFF) of different types (e.g. Fig. 2, Appendix S1).
Characteristic ridges, troughs and lobes of the VFF (Milliken et al., 2003) suggest that these are glaciers (Souness et al., 2012;Hubbard et al., 2014), which can be classiied further based on form and topographic context. Commonly, low plains bounding mid-latitude uplands on Mars are characterised by lobate debris aprons (LDA; Squyres, 1979;Carr, 2006), which are accumulations of ice mantled by lithic debris, including dust (Kochel and Peake, 1984;Pierce and Crown, 2003;Holt et al., 2008;Plaut et al., 2009;Parsons et al., 2011;Fastook et al., 2014;Baker and Head, 2015;Petersen et al., 2018). Mid-latitude valleys within uplands on Mars are commonly occupied by lineated valley ills (LVF; Head et al., 2006), which closely resemble valley glaciers on Earth. In Phlegra Montes, several LDA are fed by higher LVF, creating composite systems. Viscous, concentric crater ills (CCF; Dickson et al., 2009) are common in the region. The crater retention ages of martian mid-latitude VFF indicate that they formed over the past 1 Gy Baker and Head, 2015), most recently over the last ~0.6 Gy, in the Late Amazonian epoch (Milliken et al., 2003;Hubbard et al., 2011;Souness et al., 2012;Hubbard et al., 2014). Dating by Gallagher and Balme (2015) of the VFF referred to in this paper as LVF 393 (Fig. 1,  Fig. 2a, Appendix S1) places it also into the Amazonian, with a crater retention age of 150 ± 20 Ma, although this is probably a very approximate age; the few larger craters on LVF 393 might follow an older isochron, still Amazonian, but on the order of 1-1.5 Ga.
As the majority of Amazonian VFF are considered to have been cold/ dry based or minimally polythermal glaciers (Hubbard et al., 2011;Hubbard et al., 2014), the dynamics of these systems are believed to have been dominated by creep (Milliken et al., 2003;Parsons et al., 2011). Warm-based sliding, required to produce the regional-scale landscapes associated with widespread glacial erosion so common on Earth (Sugden, 1974;Glasser and Bennett, 2004), is not predicted to have occurred widely in recent martian environments. Although Conway et al. (2018) concluded that widespread melt of crater wall glaciers occurred 5-10 Ma, observations of landforms, in contextually consistent landsystems (Eyles, 1983) diagnostic of warm/wet based glacial regimes on Mars, especially eskers (Close, 1867), are rare in comparison to the widespread presence of VFF (Souness et al., 2012;Levy et al., 2014). However, Gallagher and Balme (2015) identiied an esker in Phlegra Montes, emerging from the distal deposits of a glacial landsystem morphologically integrated with LVF 393 (Fig. 1). Given the wide consensus on the dominance of cold/dry-based mid-latitude glaciation on Mars, the conclusion was that the Phlegra Montes esker represents the former occurrence of spatially limited warm-based conditions, associated with enhanced geothermal heating along the graben in which the esker is located. Thermal modelling by Butcher et al. (2017) of a glacier-linked esker in a similar setting in northwestern Tempe Terra indicates that elevated geothermal heat lux and viscous strain heating could have permitted restricted wet-based glaciation, even in a general environment in which climate was too cold for more extensive wetbased process environments. These conclusions raise the question at the core of this paper; is the esker in Phlegra Montes an extremely rare exception to a cold/dry-based norm -or did warm/wet-based VFF exist more widely across the region? In answer to this question, the paper irst describes the extant VFF in Phlegra Montes, including a determination of the thickness and composition of one of the largest, with estimates of the properties of several others. The focus of the paper, however, is to determine if the landscape context of these VFF provides evidence of widespread warm-based processes. This is done by reference to geomorphological criteria diagnostic of warm based glacial processes on Earth (Table 1), identiied principally by Glasser and Bennett (2004). The paper describes candidate analogues of these key terrestrial landforms and discusses the implications of their occurrence in Phlegra Montes with respect to former subglacial thermal regime. Contextualization of the extant VFF within this geomorphological framework provides valuable insights to the evolution of glaciation in the region.

Data
Landforms were mapped and measured from mosaicked CTX (Context Camera, NASA Mars Reconnaissance Orbiter; Malin et al., 2007; image credits: NASA/JPL-Caltech/MSSS) visible-wavelength images (resolution ~6 m per pixel), displayed in an equirectangular projection. To preclude projection-induced distortions, all directions were measured relative to true north and all spatial measures were geodesic. CTX images were overlaid in a Geographic Information System (GIS) onto the MOLA-HRSC Blended Global 200 m/pixel digital elevation model (DEM) v2  Tanaka et al., 2014): lNt, late Noachian transition; HNt, Hesperian-Noachian transition; eHt, early Hesperian transition; AHv, Amazonian-Hesperian volcanic. Elevation data Mars Orbiter Laser Altimeter (MOLA)/ USGS 463 m global mosaic. Crosshatch ill = LVF; light grey ill, thin edge = LDA; vertical cross hatch, thin edge = compound LVF-LDA (COMP); light grey ill, thick edge = CCF. Arrow shows location of the esker associated with LVF 393, described by Gallagher and Balme (2015). Labels with numbers refer to the type and surface area (km 2 ) of the largest VFF. Note the lack of VFF in summital areas, e.g. east of LDA 1694 and along the northern mountain spine, east of COMP 1540. MOLA data credit: USGS, NASA, JPL. (Fergason et al., 2018), derived from Mars Orbiter Laser Altimeter (MOLA; Smith et al., 2001) and High-Resolution Stereo Camera (HRSC; Neukum et al., 2004) data. Elevation proiles were plotted from digital elevations models (DEM) based on HRSC elevation data (50 m/pixel) or (where available) High Resolution Imaging Science Experiment (HiRISE; McEwen et al., 2007) elevation data (DEM at 2 m/pixel). Only meso-scale to macro-scale landforms (>20 m to 1 km) are widely considered in the observations that follow, although HiRISE images and anaglyphs were used, where available and relevant, to describe and measure smaller landforms. To conirm that the VFF are composed of ice, rather than being lithic or volcanic valley-ills (cf. Byrne et al., 2013) determinations of VFF ice content, purity and total thickness were attempted using data from the Mars Reconnaissance Orbiter (MRO) Shallow Radar instrument (SHARAD, Seu et al., 2007;Appendix S2.1). Owing to topographic clutter, these determinations were possible for only one VFF. However, the average thickness of this VFF and several others was estimated (Appendix S2.2) and complementary determinations of its surface properties were made (Appendix S2.3).

Extant VFF and topographic contexts
Extant VFF were identiied based on criteria speciied by Milliken et al. (2003), Head et al. (2010), Souness et al. (2012) and Hubbard et al. (2014). The broad topographic contexts of VFF in Phlegra Montes were characterised using imaging and elevation data. This mapping also allowed the spatial distribution of extant VFF to be determined in topographic-geomorphological context. To be considered glacial, VFF must be characterised by: (1) features suggestive of viscous low over or around obstacles, in response to changes in underlying slope; (2) contextual topography indicative of modiication by viscous low over or around obstacles; (3) surface texture and morphology distinct from summital areas, inter-valleys and peripheral terrain; (4) longitudinal foliation, consisting of narrow ridges and furrows, where VFF are affected by local low constriction; (5) extensively lineated surfaces, where VFF are laterally conined (producing LVF), and; (6) moraine-like ridges (MLR) that bound or partition VFF surfaces displaying foliation and/or lineation associated with viscous low. VFF exceeding 50 km 2 (VFF > 50; Appendix S1) were mapped as area features in a GIS, and then classiied by type (LVF/LDA/COMP/CCF), with surface area also measured -e.g. LVF 393 refers to a LVF of 393 km 2 .

Shallow radar identiication of VFF subsurface relectors, and determination of VFF bulk dielectric value, total thickness and ice purity
To determine the ice content, purity and total thickness (t) of GLF in Phlegra Montes, we analyzed SHARAD data along tracks crossing Phlegra Montes (Appendix S2.1). The bulk dielectric value of VFF is measured to constrain t and ice purity. SHARAD signals relect at contrasts in the real part of the relative permittivity for various materials, e. g. between air and rock or ice and rock. Relative permittivity (dielectric constant) varies on Mars from 1 (atmosphere) to ~12 (dense basaltic lava lows). The dielectric value for most non-porous water ice on Mars ranges from ~3 to 3.2 Holt et al., 2008;Plaut et al., 2009;Grima et al., 2009), values greater than 3.2 generally implying some volume fraction of lithic material. Complementary assessments of VFF average thickness, based on surface area measures (after Karlsson et al., 2015) and surface lithic composition (after Baker and Carter, 2019), are provided in Appendix S2.2 and S2.3.
Radar proiles (radargrams) are interpreted in several steps, including clutter mitigation and depth correction. First, each radargram identiied to have potential subsurface interfaces is compared with a simulated radargram, based on MOLA topographic data. Simulated radargrams use topographic facets to predict where relections originating on the surface, both from nadir and to the side, will lie on the observation in time delay. Frequently, across-track surface relections appear in greater delay than the nadir point (thus beneath the surface) due to their longer travel distance. These relections, commonly called clutter, do not indicate a subsurface interface and are excluded from interpretation. Thus, where clutter is predicted and matches a potential subsurface interface, we do not interpret a subsurface relection. On the other hand, when clutter is not predicted for a location, but a relection appears, we interpret that to indicate a subsurface interface and record its position and time delay.
Once clutter echoes have been identiied and excluded, the second step is to use surface topography and knowledge of the surrounding terrain to estimate a subsurface geometry associated with the relection and derive a dielectric constant. In the case of LDA, this step requires an assumption that the surrounding terrain maintains its geometry in locations covered by the LDA (e.g. Holt et al., 2008;Plaut et al., 2009). To determine the bulk dielectric value of LDA 1694, we use an equation that relates distance to velocity and relative permittivity (Seu et al., 2007).
where t is the two-way travel time, d is distance, ε' is dielectric constant and c is the propagation speed of electromagnetic waves through a vacuum. Having found the peak terrain elevation along each radar proile and measured the time delay to calculate the relative permittivity ε', we also performed a test to latten the basal relector to match the surrounding terrain. To do this, radargrams are "depth corrected" by moving all relections that arrive after the nadir surface point towards the surface using the same equation as before. This method is less rigorous but offers a visual way of comparing dielectric values rather than focusing on time and thickness at a single point (Fig. 6).

Table 1
Landforms and values of morphometric parameters considered diagnostic of erosion by warm-based glacial ice (after Glasser and Bennett, 2004;Galofre et al., 2018;Bouquety et al., 2019)

Morphological analogues of erosional landforms produced by warmbased ice
Candidate morphological analogues of terrestrial landforms considered diagnostic of erosion by warm-based glacial ice were identiied by reference to the geomorphological contexts in which the VFF occur and to the insights offered by Glasser and Bennett (2004); the location of candidate analogue landforms is shown in Appendix S2.4. The set of diagnostic landforms (Table 1) is supported by morphometric indexes that have been used to infer the genesis of valleys, valley-heads and channels. Bouquety et al. (2019) applied the valley index (V-index) to discriminate between glacial and luvial valleys on Mars. V-index compares valley cross-sectional area (Ax) against an ideal V-shaped valley cross-section (Av) of the same height and width (Zimmer and Gabet, 2018). Ideal U-shaped (glaciated) valleys have a V-index >0 but ideal Vshaped (non-glaciated) valleys have a V-index <0 (Zimmer and Gabet, 2018). Bouquety et al. (2019) interpreted martian valleys with a V-index > 0.2 as glacial but <0.1 as luvial; we place the luvial threshold at ≤0.1. For Phlegra Montes, V-index was calculated from elevation data extracted from HRSC 50 m/pixel DEM and is deined as.
Cross-sectional area was measured from the deepest point of each valley cross-section to the elevation of the lowest valley shoulder. This is equivalent to the full-valley cross-section measure of Zimmer and Gabet (2018), which had a predictive accuracy of 65% for glaciated valleys, 68% for non-glaciated; the truncated-valley measure of Zimmer and Gabet (2018) is slightly more predictively accurate (66% for glaciated valleys, 70% for non-glaciated) but is more subjective. Bouquety et al. (2019) also applied morphometric criteria developed in the terrestrial context (Barr and Spagnolo, 2015;Evans and Cox, 2017) to identify glacial corries on Mars. Applying the method of Barr et al. (2019), the ratios of corrie width-to-height (W:H) and length-towidth (L:W) were measured for candidate corries in Phlegra Montes. L is the length of a line plotted from the mid-point along the corrie threshold (the crest of the convexity deining the corrie distal margin) to the mid-point along the corrie backwall crest. W was measured normal to the mid-point of L. H was measured as the elevation difference between the corrie loor (under the intersection point of the measurement lines of W and L) and the elevation of the backwall crest. For corries on Earth, Barr and Spagnolo (2015) reported mean morphometric ratios of W:H = 2.72 and L:W = 1.03. Bouquety et al. (2019) interpreted valley heads on Mars as glacial corries/cirques by reference to these morphometric criteria. Galofre et al. (2018) applied an undulation index (ψ) to identify channels of subglacial genesis on Earth. To measure ψ, channel longitudinal proiles were plotted from elevation data extracted from a 2 m/ pixel HiRISE DEM. ψ is deined as.
where h min and h max are, respectively, the local minimum and maximum elevations along a channel reach with positive topographic gradient and ΔH is the total along-channel elevation change. Galofre et al. (2018) concluded that channels with ψ > 0 are unequivocally of subglacial genesis. Glasser and Bennett (2004) classiied landforms of warm-based glacial erosion on Earth at micro-scale (<1 m), meso-scale (1 m to 1 km) and macro-scale (>1 km). However, owing to constraints placed by relying mainly on CTX images, we limited the search to diagnostic mesoto-macro-scale landforms (Table 1). In deciding if the morphological evidence supports an interpretation of morphogenesis in a warm-based subglacial environment, the most conservative palaeoenvironmental interpretation should be assumed, consistent with the general viewpoint that regional warm-based glaciation on Amazonian Mars has been extremely rare. Therefore, strong morphological indicators must be observed to reject the most conservative interpretation in favour of the alternative of warm-based glaciation. In that respect, ψ and the morphometric indexes offer important discriminant metrics.

Distribution, topographic contexts and forms
VFF > 50 occur in Phlegra Montes from ~32 • -41 • N ( Fig. 1; Appendix S1) and have a total surface area of ~16,600 km 2 . The VFF have four main topographic contexts. (1) VFF occur in the intra-montane basins of massifs, on the lanks of isolated mountains or small massifs and in long, WNW-ESE trending grabens (Fig. 2a). VFF are not widely present in summital areas ( Fig. 1; see also Appendix S2.4). Intramontane basins occupied by these VFF include both low linkage and outlet to lower levels through short, interconnecting passes, and through longer troughs. VFF in these settings, partially illing topographically constrained elongated negative relief, were mapped as LVF.
(2) Lobate piedmont aprons with viscous-appearing surface forms were mapped as LDA. They occur around the lanks of the northern spine of Phlegra Montes and on piedmonts surrounding isolated mountains or small massifs ( Fig. 2b). (3) Some massifs or more compact uplands are fringed by LDA fed from the uplands through linear corridors by LVF; these compound VFF were mapped as COMP (Appendix S1). (4) Some craters, both upland and lowland, are partially occupied by viscous, concentrically ridged material, mapped as concentric crater ill (CCF; Fig. 2b).

Sharad observations -LDA 1694 subsurface relector and thickness
A regional survey of available SHARAD data found seven relevant radargrams not adversely affected by topographically induced radar clutter (Appendix S2.1). These all cross over or near to LDA 1694 ( Fig. 2b) and exhibit a strong subsurface relection (Fig. 3). For the observations with detections under LDA 1694, we assume a lat basal geometry that extends from the surrounding, approximately lat terrain towards the west and south. Fig. 3 shows that the foreland of LDA 1694 is a nearly lat surface at ~ − 3300 m, and that the highest point of the LDA along a SHARAD orbital ground track (not at the contact with the abutting massif but slightly offset) is at −2802 m, for a maximum potential thickness of ~500 m (Table 2). Applying the time delay (Section 3.2) yields an average relative permittivity (ε') of 3.21 ± 0.11 (Table 2; range 3.07-3.36). This value is consistent with >96% pure water ice (Bramson et al., 2015). Applying a depth correction (Section 3.2), we ind that a value of ε' = 3.15, consistent with pure water ice, does not provide a perfect match with the surrounding terrain, but a value of ε' = 3.6 does a better job. This value is consistent with an ice purity of 88-90% by volume or a stratiied deposit of pure ice (ε' = 3.15) overlain by regolith. Two possible end members are 396 m of ice overlain by 4 m of high-density regolith (ε' = 8.0) or 390 m of ice overlain by 10 m of low-density regolith (ε' = 5.0). Complementary morphometric and lowbased estimates (Appendix S2.2) indicate that average VFF > 50 thickness in the region ranges from <100 m to 200 m (LDA 1694 having the largest average thickness).

LDA 1694 bulk dielectric properties, ice purity and surface variation
The dielectric value of LDA 1694 is 3.2-3.6, higher than the nominal value of water ice (3.15) and higher than the LDA studied by Holt et al. (2008), Plaut et al. (2009) and Petersen et al. (2018). A dielectric value of 3.6 is consistent with a 396 m thick glacier with ~4 m of a dense, rocky material lag on top. If so, the bulk of the LDA must be pure ice and the lag must be too thin for SHARAD to detect it (i.e. <15 m). However, another possibility is that the dielectric value represents a homogeneous mixture of ice and up to ~10% dust. This interpretation would be consistent with LDA in other parts of Mars. Further insights to the clastic surface properties of LDA 1694 are provided in Appendix S2.3. Finally, if the ground under LDA 1694 were not lat but, instead, glacio-isostatically depressed, the dielectric could represent a glacier of pure ice (ε' of 3.15). However, there is no evidence of crustal lexure beneath or immediately surrounding the LDA at this location or for LDA at other sites around Mars.

Candidate landforms interpreted as indicators of warm based erosion
In this section we present the results of a search for examples of candidate analogues of the terrestrial landforms considered diagnostic of erosion by warm-based glacial ice (Table 1). Following the description of each analogue example, an interpretative justiication of their candidate status as evidence of warm-based erosion processes is provided. The candidate landforms all occur within contextual landscapes characterised by the presence of VFF, from the criteria of Milliken et al. (2003), Souness et al. (2012) and Hubbard et al. (2014). Although the interpretations that follow cannot be deinitive, being based only on remote sensing, we believe they are morphologically consilient with diagnostic terrestrial analogues and reasonable in their landscape context. We do not attempt to order the landforms by any implied signiicance; it is important that each landform is viewed as part of an integrated landscape system. Hence, the landforms are dealt with in the order presented by Glasser and Bennett (2004), which also corresponds broadly to their prevalence in Phlegra Montes (Appendix S2.4 shows Table 2 Estimates of maximum potential thickness of LDA 1694. Fig. 3 shows the terrain surrounding LDA 1694 to be nearly lat ~ −3300 m, with the highest point of the LDA along a SHARAD orbital ground track at −2802 m, giving a maximum potential thickness of ~500 m.  Fig. 4a-c show a landscape of smoothed, often elongated hills separated by viscous low components of the LVF 4258 complex ( Fig. 1 and Appendix S1). Several of these elongated, smoothed hills narrow to a rounded summital ridge and some are traversed by subtle bands orthogonal to the long axis of the hill (e.g. Fig. 4b,c). The long axis orientation of these hills is complemented by branching, arcuate or chevron-shaped ridges within the LVF. The branching ridges compartment basins in the low terrain between elongated hills and smoothed platforms (e.g. Fig. 4a). The orientations and junctions of branching ridges are complemented by viscous low lineations on the adjacent LVF surface (Fig. 4b).

Areal scour
4.2.1.1. Interpretation. In this landscape, the narrow, branching LVF ridges are contextual and morphological analogues of lateral and medial moraines (Fig. 4a, labelled MLR). In context, the arcuate and chevronshaped ridges are strikingly like arcuate moraines and barchanoid moraines (Fig. 4b,c, respectively labelled MLR-Ar and MLR-Bc). Elongated hills with narrow, rounded longitudinal ridges are morphological analogues of whalebacks, streamlined landforms resulting from subglacial scour. Bands draped orthogonally across the long axis of streamlined hills and rock platforms are consistent with being recessional moraines (MLR-Rc), deposited by thicker ice that formerly covered the platforms (Fig. 4b-c). In the context of LVF low direction indicators (topography, general slope, surface lineations), MLR-Bc point apex-downlow (their arms extending up-low). These characteristics are consistent with the interpretation of MLR-Bc as morphological and palaeolow analogues of terrestrial barchanoid moraines, which form when pre-existing glacial sediments are overridden during late stage warm based glaciation (Golledge and Stoker, 2006). The presence of viscous material bounded by a pitted margin beyond the MLR-Bc suggest a similar sequence of events, involving small, localized overriding lows across a pre-existing ice-rich substrate. The downlow pointing curvature of the MLR-Rc and the MLR-Bc are mutually consistent, and congruent with small, streamlined forms in the lee-side indentation into hill X and lineations between hill X and the MLR-Bc (Fig. 4c).
Overall, in the context of Fig. 4a-c, the textures of the eroded surfaces and smoothed forms, and the form of the whaleback hills and MLR, are consistent with streamlining glacial scour across a considerable area. An important implication arising from this landscape interpretation is that these forms have become exposed due to signiicant LVF thinning (or, possibly, increasingly deep dissection of troughs and basins between Lines Z-Z', Y-Y ′ and X-X' are cross-sectional proile paths. X-X' and Y-Y ′ are not signiicantly affected by a LVF and each has a V-index > 0.2, consistent with glacial incision. Elevation data: HRSC (data credit ESA/DLR/FU Berlin). CTX image credit: NASA/JPL/MSSS. increasingly exposed intervening highs). A competing, overarching interpretation is that the smoothed forms and ridges are aeolian in origin, but this would ignore the directional variability of these forms, particularly the apex directions of the MLR-Bc, and the close relationship between the coniguration of those forms and the LVF low indicators. Local katabatic winds would not provide a good explanation either, as most MLR-Bc are located downlow of prominences rather than being in valleys.

Glacial troughs
LVF 164 (Fig. 5) occupies valley N, a feature not sharing the same alignment as the major fault-valleys hosting VFF in Phlegra Montes (e.g. Fig. 2a). The valley possibly has a catenary cross-section, but the full cross-sectional form is obscured by the LVF (proile Z-Z', Fig. 5). Basin A is a perched, breached depression, occupied by a viscous ill. Trough B, separated from basin A by ridges, is occupied by a viscous ill, and is conluent with LVF 164 in valley N. Basins C and D are unilled. The backwalls of both C and D are breached by troughs E incised into the summit of the massif. Basin F is a breached crater, partially occupied by a lineated ill. The mouth of valley N is partially enclosed by rock ridge G, beyond which the LVF recurves and terminates at the entrance to the parallel valley S, crossed by proile Y-Y ′ ; valley S has a catenary crosssection (V-index 0.38) and does not appear to have a signiicant ill. Valley S descends south-westward, enclosed by ridge H, which has a maximum relative height of ~175 m. West of ridge H, cross-section X-X' has a catenary proile but is partially occupied by a transluent LVF that exits valley N through pass I (similarly, LVF-occupied pass J links valleys N and S). Long valleys, other than E-W grabens, are rare in Phlegra Montes. However, short, LVF-occupied passes between mountains are very common. In Fig. 6a a LVF pinches-out in a pass between two mountains and the lee-side slope, beyond the pass, is incised by channels. The pass in Fig. 6b is smooth textured, intervenes between rounded, subtly lineated hilltops and overlooks basins occupied by LVF. Both passes have catenary cross-sections: proile W-W ′ , V-index 0.31; proile V-V ′ , V-index 0.28.

Interpretation.
The V-index of valley S (cross-section Y-Y ′ , Fig. 5) is consistent with it being a glacial trough (cf. Bouquety et al., 2019). The V-index, and landform contexts, of the passes in Fig. 6 (crosssections W-W ′ and V-V ′ ) are consistent with them being transluent glacial cols (terrestrial analogues of which often attain catenary form; e. g. Hambrey, 1994). Both Y-Y ′ (Fig. 5) and W-W ′ (Fig. 6a) retain minimal ill but this is typical of terrestrial glacial analogues determined to be catenary. Cross-section proiles Z-Z' and X-X' (Fig. 5) also most likely represent glacial catenaries but with more signiicant viscous ills. The shallow catenary form of trough Y-Y ′ (Fig. 5) is typical of glacial troughs in several settings on Earth, most notably in Patagonia and Antarctica (Hirano and Aniya, 1988). In context, the valley forms in Fig. 5 are consistent with the effects of warm-based glacial erosion. The LVFoccupied summital basin A and trough B, and trough complex E, indicate that the massif north of valley N is glaciatedand was formerly more extensively glaciated. Ridge H (Fig. 5) is interpreted as a MLR, deposited between the main south-westward lowing glacier in valley S and outlet glaciers descending from valley N through transluent cols (e. g. pass I). Hence, ridge H is indicative of previously more extensive valley glaciation, the ridge height indicating former glacial thickness of ~175 m (H appears to have been overridden, and smoothed, by a latestage outlow through transluent col. I, presumably permitted by the withdrawal of the main glacier snout east of H in valley S). The  Fig. 6 also are indicative of former thicker glaciation. The lee-side channels beyond the col. in Fig. 6a possibly are indicative of glacial meltwater discharges by ice at the col. or, formerly, beyond it. A competing interpretation of the trough and col. landscapes is that these forms are glacially occupied grabens, with minimal modiication by glacial erosion. However, to be convincing, this interpretation would have to include all candidate cols with catenary cross-sectionsand this would imply faulting of very high spatial-density and variability in strike. Instead, it is more likely that the troughs and cols represent erosion by glaciers that lowed through pre-existing valleys (including grabens, e.g. Fig. 2a), and became suficiently thick to exploit saddles between neighbouring mountaintops. Judging the relative robustness of the competing interpretations depends on how the within-trough and col-adjacent forms are interpreted and, hence, will be revisited.

Corries
In Fig. 7 a small, perched basin A is occupied by a viscous ill, which is conluent with the small LVF occupying trough B (which descends to LVF 164; Fig. 5 shows wider context). In the same massif as A and B, the open, amphitheatre-shaped basins C and D (Fig. 5, Fig. 7) have lengthto-width ratios of 0.93 (C) and 0.85 (D) and width-to-height ratios of 3.55 (C) and 2.78 (D). Maximum backwall slopes are ~19 • (C) and ~ 21 • (D). Both basins are fronted by a slight convex-up lip. The sloping forelands beyond these lips are incised by parallel sets of lineations. The backwall of basin C is breached by the LVF-occupied trough B; both C and D are also breached by unilled troughs (E). Basin F is a breached  Fig. 7) are consilient with terrestrial corrie analogues, and with forms interpreted as corries/cirques on Mars (Bouquety et al., 2019). These basins it within a continuum of erosional glacial valley forms of varying length, ranging from cols to troughs (cf. Aniya and Welch, 1981). Additionally, the longitudinal proile of basins C and D, including the development of a slight frontal lip, is morphologically consistent with enhanced low resulting from transference of water to the glacier base in crevasse zones (Alley et al., 2019) and, in thinner ice, erosion by rotational basal sliding (Lewis, 1949). Although the lip is less developed in C and D than in analogues described by Bouquety et al. (2019), it is typical of many terrestrial corries (examples in Appendix S2.5). The presence of parallel lineations beyond these lips is consistent with subglacial abrasion having occurred beyond the candidate corries, providing further evidence of formerly more extensive glaciation. Morphologically, basin A containing the viscous ill is similar to basins C and D, with a similar W:L of 0.97, suggesting it could be a glacially-occupied corrie. However, the equidimensional axial form of corries (Barr and Spagnolo, 2015) is closely matched by the cross-sectional form of impact craters, which tend to be circular in planimetric form and parabolic in section. On Earth, Pleistocene corries developed at the head of pre-existing valleys formed in earlier tectonic and luvial contexts. In the Mars context, it is likely (even probable) that most candidate corries, in the sense of being overdeepened glacial basins, developed from pre-existing impact craters (e.g. open basin F, Fig. 7). In the substrate-context of Phlegra Montes, impact craters with no ill tend to have very V-shaped cross-sections, but craters with ill are U-shaped to rectangular in cross-section (examples in Appendix S2.6). Therefore, it is possible that the characteristic morphometric properties of corries are also integral properties of illed impact craters, even those with no outlet that could not have been overdeepened. Consequently, the morphological evolution of the open upland basins C and D in Fig. 7 can be only tentatively ascribed to glacial erosion solely based on their morphometric afinity to corries. The general context of these forms, however, suggests that interpreting them as corries is not unreasonable.

Lee-side rock basins/cavities and streamlined forms
In Fig. 8a, smooth platforms (e.g. A, B, C, D) lack viscous low features. Craters on the platforms are rimless. The LVF, characterised by prominent MLR and low lineations (FL), are constrained and directionally guided by the platforms. However, the platforms are textured by parallel, shallow grooves and lineations (e.g. L* on platform A, Fig. 8a,b, c), and are onlapped by ridges (e.g. MLR(B ′ ) on platform D, Fig. 8a, b, c). The northeast side of platforms A, B and C are indented (features X, Y, Z), the indentation into platform A being especially pronounced (Fig. 8a, b, c feature X). Lineation orientations on platform A are coherent with MLR and low lineations on adjacent LVF surfaces and platform D. Directional statistics presented in Appendix S2.8a show that lineation orientations sampled from platform A are directionally consistent with MLR and low lineations on surrounding LVF surfaces. Fig. 8a are inferred to be rock. The onlapping MLR and shallow grooves/ lineations suggest that the platforms have formerly been overridden and linearly incised. This interpretation is supported by the statistical coherence between the orientation of lineations on platform A and low indicators (MLR and low lineations) on surrounding LVF surfaces (Appendix S2.8a). Rimless craters on the platforms are consistent with this interpretation, suggesting areal scour occurred with linear abrasion. Given the palaeolow context indicated by the pattern of platform grooving and complementary MLR, indentations on the northeast side of platforms A, B and C are lee-side features (Fig. 8b and c). The indentation into platform A is a well-developed lee-side basin, open to the northeast, while the indentations into platforms B and C are less developed lee-side cavities -possibly incipient basins. North-eastward movement of ice and debris, from the lee-side cavity in platform B, is indicated by MLR(B ′ ) on platform (D). MLR(A'), extending from the platform A basin, and converging with the main assemblage of MLR (arrowed), is likewise indicative of northeasterly debris transport. Hence, the basins and cavities appear to represent lee-side entrainment and erosion of platform material by overriding palaeolows.

Interpretation. Lacking viscous low features, the platforms in
In context, the morphologies of the smooth rock platforms (including rimless craters) are consistent with scour by sliding ice. The smooth, rounded surfaces of the platforms, with zones of deeper linear incision, are consistent with streamlining scour with localized abrasion. The leeside locations of the basins and cavities are consistent with plucking, where overriding ice lost close contact with its bed on the lee-side of protuberances. Formation of lee-side basins by plucking to the degree evident in on Platform A (Fig. 8a) is consistent with "quarrying", the process in which large bedrock blocks are eroded and entrained under sliding glaciers (Hallet, 1996). An overarching implication of plucking that has evolved to quarrying is that LVF thickness was formerly greater than at present, implying that quarry-like basins were originally subglacial lee-side cavities. An alternative interpretation is that they formed by backwasting through basin headwall plucking, without glacial overriding of the platforms. Both interpretations involve subglacial meltwater, but their implications differ with respect to former LVF thickness. Given the closeness of the form-low relationships in this landscape, non-glacial erosion alternatives are hard to envisage but could include aeolian erosion and aspect-inluenced surface degradation through loss of volatiles. Further examples of lee-side modiication, ranging from steepening, without excavation, to cavity development associated with parallel grooves incised on lee-side slopes are provided in Appendix S2.7 and S2.8.

Meso-scale rock grooves
Parallel sets of grooves are a notable component of the landscape of Phlegra Montes across a variety of landscape scales. In Fig. 9, the surface of a platform ~1.5 km in down-valley extent (partly obstructing LVF 164) is textured by pits, 80-250 m wide (Fig. 5 provides context, label  G). These coalesce towards the northeast into broad grooves, typically 150-200 m wide, 300-350 m long, separated by narrow inter-grooves. The pit-grooves parallel the local southwest to northeast low lineations on LVF 164.
At a larger landscape scale, the northwestward sloping rock ramp in Fig. 10 (labelled A) extends downslope for ~6 km. This ramp is smoothed, lacking positive relief forms, but scored by long, narrow grooves (B). Simple, irst order, conluences give the grooves the appearance of a channel-like network. The adjacent dome-shaped rock platform C also is smoothed but deeply incised by troughs D and E, which span dome C. Trough F, bisecting A and C, joins the trough occupied by LVF 393, at a hanging valley junction. The interface between trough F and ramp A comprises a subtle depression that converges into the head of the southernmost channel-like groove (at X in Fig. 10).
The grooves in Fig. 9 and Fig. 10 are microcosms of the extensive sets of parallel grooves incised into the rounded, dome-like mountains summits in Fig. 11a,b. Notably, these summit-crossing grooves are not arranged radially with respect to the mountain domes but, rather, produce a landscape with a "raked" appearance ( Fig. 11a, arrows labelled C). Craters on mountain summits are rimless. Grooves become convergent only where the relief of the terrain is centripetal, i.e. where slopes become focussed radially inward to a common point (Fig. 11a, arrows labelled D). Where convergent macro-scale grooves occur, they develop sinuosity downslope, becoming channel-like in both form and network pattern (Fig. 11a,b). However, where terrain is not centripetally focussed, grooves are parallel and continuous for several kilometres, ascending and descending the smoothed, rounded mountains with no topographically induced deviation in alignment. Fig. 9 appear to have developed from individual "chiselled" pits that extended into longer, roughtextured pit-grooves or shallow, elongated basins. The pit-groove alignment, and their elongation parallel to the LVF lineations, suggest that they too record LVF movement towards the northeast. The development of grooves from pits on the cross-valley ramp is consistent with an erosional environment transitional between plucking and linear abrasion. Moreover, referring to platform A in Fig. 8, the assemblage of lineations there suggest a close relationship between linear abrasion and areal scour, leading to streamlined rock forms locally textured by shallow, parallel grooves. In this process continuum, subglacial protrusions become more deeply grooved and abraded over time as an accommodation governed by feedbacks between substrate topography, ice pressure and the through-put of basal debris from cavities to grooves. In Fig. 9, the distal narrowing of the grooves, beyond the foot of the ramp, is consistent with a diminishing lux of luid and debris, as more uniform contact with the substrate was also re-established. In this setting, the development of pits and grooves is likely to relect wet based subglacial conditions, at least as a local consequence of increased strain heating in the vicinity of the cross-valley ramp, which would have been a signiicant low obstruction (Appendix S2.9 shows a historically signiicant terrestrial glacial analogue). Given the form of the platform in Fig. 9, and the parallel arrangement of the pits and grooves on its surface, aeolian erosion could be considered an alternative to the glacial hypothesis. This would ignore the consistency between the directionality of the platform grooves in Fig. 9 and the LVF low indicators but would not fail as a hypothesis for that reason. Implicitly, however, this interpretation would also require non-glacial processes to have driven the morphological evolution of the LVF-occupied trough. While this could be envisaged at a single locality, the general consistency between form and indications of LVF low suggest that glacial erosion is the more likely explanation on a valley-scale.

Interpretation. The grooves in
In Fig. 10, MLR G in troughs D and E suggest that the troughs in this larger landscape context were incised by transluent glaciers that had expanded into trough F. In this context, it is possible that the channellike grooves B were abraded by wet, sediment-rich, slurries at the base of a locally westward lowing, difluent ice lobe that likewise expanded out of trough F, down-wearing the trough-ramp interface (at X, Fig. 10) in the process. Ice appears also to have lowed divergently through trough F, the point of divergence along the trough, between northward and southward low (approximately at F, Fig. 10). Glacial transluence (via eroded valley-head breaches) and difluence (via valley-side overspill) represent a situation in which accumulation of ice (e.g. in trough F) exceeds its drainage (i.e. out of trough F), as a consequence of increased accumulation or diminished outlow (e.g. due to thickening of LVF 393). In the case of trough F, ice suficiently thickened to have risen above the walls of the trough would have lowed to lower elevation across A and C. This scenario, predicated on the former greater thickness of LVF 393, could also explain the depth of incision of the trough hosting LVF 393. Non-glacial explanations for the landscape in Fig. 10 are hard to envisage, particularly given that the grooves (B) incised into platform A traverse a convex-up surface (from X). This makes mass wasting or subaerial luvial erosion unlikely explanations.
At the largest landscape scale, the parallel alignment of mountaincrossing grooves in Fig. 11 (arrows, C) suggests the grooves are not gravity-driven forms, which would have radial, not parallel, fall-lines descending the mountain-slopes. In context, this suggests that the parallel grooves represent glacially driven linear abrasion beneath formerly topography-inundating, sliding ice. The development of channel-like networks from parallel grooves, where topography is centripetally focussed (Fig. 11a, arrows D), suggests that liquid was produced in the subglacial environment and involved in groove incision, given the continuity from parallel grooves to convergent grooves to channels. This morphological continuum, therefore, represents a process continuum modulated by transitions between glacially driven stress and gravitational stress.

Meso-scale subglacial meltwater channels
A candidate subglacial channel is shown in part of the HiRISE anaglyph in Fig. 12, at the interface between the lineated surface near the northern margin of the LVF 2431 complex (Fig. 12a, region A) and a more densely cratered surface unit (Fig. 12a region B) to the north ( Fig. 1 and Appendix S2.4 provide context). Reach S-S ′ is 5.5 km in length, its mean width between channel shoulders is ~115 m and its Vindex ranges from 0.1 to −0.063 (i.e. V-shaped). This channel reach appears irst at the LVF terminus, in a shallow basin, and cuts across the lank of a hill (Fig. 12 at C), which closes the northern rim of the basin. The channel becomes less apparent (Fig. 12 at D) for 3.4 km (although possible traces of the channel remain apparent on the surface), before reemerging along reach N-N ′ (Fig. 12a, b). Reach N-N ′ is 2.3 km long, has a mean width between shoulders of ~75 m and a V-index of −0.063 (i.e. Vshaped).
On the lank of hill C opposite the initial channel emergence at S (Fig. 12a), the terrain is incised by parallel NW-SE grooves (Fig. 12a between labels C and E). Planimetrically, the channel and grooves converge on the downslope centreline of hill C but the channel crosscuts the grooves. Reach N-N ′ (detail in Fig. 12b) crosscuts another set of parallel grooves, these oriented NNE-SSW and, therefore, possibly associated with a viscous overlow from the large crater partly visible in Fig. 12a close to label B. Surface F (Fig. 12a) is subtly lineated, parallel to the NNE-SSW grooves, and is incised by a short channel.
4.2.6.1. Interpretation. The anaglyph suggests that the initial reach of the channel (Fig. 12a at S) is incised along an adverse slope, and that reach N-N ′ crosscuts grooves, also initially along an adverse slope (Fig. 12 b). Arrowed labels X, Y and Z ( Fig. 12a and Fig. 13) point to the locations of inlexions in the channel proile. X and Y are convex-up inlexions (not adverse undulations) in the curve along reach S-S ′ (Fig. 13, proile a). With a total elevation change of 366 m from S-N ′ , the undulation at the head of S-S ′ and, along reach N-N ′ , feature Z both have Ψ = 0.02 (Fig. 13, proile a). However, the incision of the channel into the substrate along reach N-N ′ must have subdued the feature Z undulation (over time luvial erosion would achieve a concave-down proile). Proile b in Fig. 13 was plotted immediately west of reach N-N ′ and shows the topography that water lows would most likely have encountered at the start of the incision process. This proile indicates that water would have had to low uphill for ~1.7 km to incise the channel as we now see it, ascending by 17 m (Ψ = 0.05).
The contours in Fig. 13c (HiRISE 2 m DEM, 20 m contour interval) Fig. 9. Platform transverse to LVF 164 low. The platform is convex-up, with an adverse stoss face sloping gently into the LVF low and a lee face sloping more steeply to the down-low to the northeast. The stoss face is indented by pits but the lee face is incised by longitudinal (i.e. LVF low-parallel) grooves. CTX image G20_026158_2190_XN_39N195W. Image credits: NASA/JPL/MSSS.
are congruent with the longitudinal proiles. Channel reach X-Y (Fig. 12a, Fig. 13a,c) is not contour-orthogonal but the grooves that this reach cross-cuts are (hence these grooves might be older gravity-driven channels similar to the channels at D in Fig. 11a). The reach between the contour closest to feature Z (−2520 m) and the −2540 m contour (Fig. 13b,c) also is not contour-orthogonal; the channel becomes contour-orthogonal downslope from the −2540 m contour. Hence, the contour lines support the indications from the longitudinal elevation proiles (in both contour deviation and proile undulation) that pumped low, rather than only gravitational low, was involved in the channel incision. Altogether, the measured characteristics of the channel (Vindex ≤ 0.1, Ψ > 0, variation in channel-contour alignment) are consistent with hydraulic incision in a subglacial environment. Fig. 14a shows a channel-like form (D) of ~7 km length, sinuosity 1.1, width 0.8-1.2 km and mean depth ~ 30 m (mean w:d 37; d estimated from PEDR data, HRSC DEM data and shadow length cast by the western bank of both channels). Cross-sections V1 and V2 (Fig. 14b) have a V-index of 0.04 and − 0.019 (i.e. V-shaped). The longitudinal path of D (L-L') is convex-up (Fig. 14b, dotted trendline) and dominated by three adverse undulations, each with Ψ > 0: Ψ (Z) = 0.22; Ψ (Y) = 0.05; Ψ (X) = 0.05). The distal reach of D (downslope of X) contains a ill, incised by smaller (~50 m width), sinuous channels. D is part of a group of low sinuosity, high w:d channel-like forms (context, Appendix S2.10), including D* (Fig. 14a), which is largely occupied by a coarse-textured ill. Fig. 14 also shows giant (~1 km wide, ~2 km long), longitudinally asymmetric, positive relief forms (A) beside channel-like form D. The positive relief forms (A) are moulded into the downslopedistal surface of a smoothed platform (B), characterised by rimless craters, and bounded to the east by channel-like form D. The positive relief forms (A) are incised longitudinally by narrow grooves, and bounded around their steep (upslope-proximal) face and lateral margins by morphologically complementary, tapering depressions (C). More widely, the surfaces incised by forms D and A are densely incised by small channel-like forms (E), in particular joining D* with high conluence angles (Fig. 14a,b).

Interpretation.
In the landscape context of Fig. 14a,b, the high w:d channel-like forms (D and D*), small channels (E) and asymmetric depression-bound positive relief forms (A) are all consistent with having been liquid-incised, albeit in different process environments. The Vindex measures from two cross-sections of D (Fig. 14, V1, V2) are ≤0.1 and, therefore, consistent with luvitile incision (Bouquety et al., 2019). However, the convex-up longitudinal proile of D and its Ψ > 0 undulations are not consistent with sub-aerial luvial incision but, instead, with incision in a subglacial environment (Galofre et al., 2018). The small channels (E) suggest late-stage luid lows with distributed sources, as opposed to earlier, larger episodes of erosion responsible for D and D*. The depression-bounded positive-relief forms (A) closely resemble sichelwanne. Described by Kor et al. (1991) and Shaw (2002), sichelwanne are hairpin-shaped shear-forms (s-forms), produced beneath sliding ice by highly erosive lows of subglacial meltwater or meltwatersediment mixtures (slurries; cf. Meehan et al., 1997). S-form location is not topographically controlled; they develop both in stoss-side and leeside situations (Menzies and Shilts, 2002;Menzies and Shilts, 2002). In context with the s-form analogues, the channel-like forms D and D* (Fig. 14) have the scale, form and geometry of subglacial tunnel valleys (Cofaigh and C., 1996;Beaney, 2002;Lelandais et al., 2016;Livingstone and Clark, 2016). Tunnel valleys are elongated depressions, found close to former ice sheet margins, and have three diagnostic characteristics: they begin and end abruptly; their along-path width varies little as a function of down-channel distance, and; their loors are undulatory, with overdeepened reaches and adverse reaches. Generally, they are rectilinear to mildly sinuous in planform, aligned parallel to local ice low indicators. They can be components of dendritic, anastomosed or braided channel networks. Cross-sections are U-shaped or V shaped, with steep sides and high width-to-depth ratios, reaching 50:1. These characteristics together, but particularly the presence of adverse reaches, indicate that tunnel valleys are products of incision by pressurized, subglacial water (Cofaigh and C., 1996;Beaney, 2002). All these characteristics are present in channel-like form D (in D*, the ill precludes reliable measurement of Ψ but the other tunnel valley characteristics are present). Given the contextual presence of the positive relief forms resembling sichelwanne (s-forms), we suggest the channel-like forms are strong tunnel valley candidates, representing a former warm-based subglacial zone. Alternatively, the channel-like forms could be glacially abraded troughs, incised, then illed, by sediment, or icesediment mixes, prior to late-stage incision of the smaller channel assemblage. Even if so, the candidate sichelwanne (s-forms) (and the grooves along them) are strongly suggestive of a wet, high shear, subglacial environment, possibly involving post-shear linear abrasion (with respect to the superimposed longitudinal grooves).

Summary and synthesis of landscape observations
The insights to the composition and structure of LDA 1694 (Fig. 2b) from SHARAD (Section 4.1.4, Fig. 3), together with its viscous low characteristics, point conclusively to its glacial nature. The consilience between the low structures of LDA 1694 and the other VFF in Phlegra Montes is strongly indicative of the other VFF being glaciers, not lithic or volcanic valley ills (cf. Byrne et al., 2013). In that context, LVF low structures and forms, including various MLR forms ( Fig. 4 and Fig. 8), correspond with the streamlining and development of stoss-to-lee asymmetry in the erosion (Fig. 8) of smooth platforms bypassed by extant LVF. The presence of MLR draped across streamlined platforms suggests that the thickness of these LVF was formerly suficient to override the platforms. The consistent correspondence between LVF low indicators and adjacent platform surface form suggests that the platform smoothness is primarily a product of subglacial scour, with local modiication by abrasion, not simply a consequence of mantling. The observations and morphometrics presented in Section 4.2 (in the context of the extensive regional distribution of smoothed, grooved and raked terrain) point to a large proportion of the Phlegra Montes landscape being an assemblage of exposed, subglacially eroded surfaces. The exposure of these subglacial landforms points to a very signiicant deglaciation having taken place. Rounded, smoothed, km-scale mountaintops are common in Phlegra Montes, including many surrounded by extant LVF or LDA. Similarly, many passes between mountain domes, and forelands beyond them, are densely raked with parallel narrow grooves. Where slopes converge on constrictions between mountain domes, or in terminal basins at the foot of mountain-domed escarpments, grooves converging into arborescent networks of channel-like forms are very common (e.g. Fig. 11).

Landform associations and environmental controls
In Phlegra Montes, narrow grooves consistently evolve downslope to be replaced by sinuous, often networked, channel-like forms. Based on Earth analogues, these landforms could relect a genetically associated trend, evolving from linear glacial abrasion to glacioluvial incision. Linear grooves occur as parallel, diverging or converging sets, the arrangement depending on underlying relief and slope; the "raked" grooves in Fig. 11 are continuous for several kilometres, but also exhibit curvilinearity and convergence depending on contextual terrain slope and topographic constraint (see also Appendix S2.11). Parallel sets of raked grooves are consistent with linear abrasion by a liquid or slurry with no degrees of freedom to migrate laterally or fall under gravity, despite the present landscape showing no such constraint. Where parallel sets of grooves incise rounded mountain domes, the expected groove pattern for unconstrained liquid lows would be radial. This evidence for the grooves being produced by laterally constrained luid abrasion is consistent with subglacial linear erosion. This interpretation is strengthened by the presence of the smooth, elongated, axially grooved positive relief forms within conining tapering depressions, including their contextual landform assemblage (Fig. 14a, b); high w: d channels, sinuous channels and channel networks (Fig. 14) are also consistent with having been liquid-incised in a subglacial setting. The high w:d channels (e.g. Fig. 14) have the scale, form and geometry of subglacial tunnel valleys (Livingstone and Clark, 2016) or abraded troughs containing a post-incision ill. The depression-bounded forms closely resemble sichelwanne (hairpin shaped s-forms) produced by highly erosive lows of subglacial meltwater or meltwater-sediment slurries (Kor et al., 1991;Shaw, 2002). If the smoothed mountain domes and grooves represent former subglacial scour and linear abrasion involving basal meltwater or meltwater-sediment slurries (Fig. 10), the channel-like networks (Fig. 11) could represent distal meltwater discharges, or discharges of the meltwater component of abrasive slurries. In that respect, the groove-channel landform assemblage would effectively map the glacial maximum coniguration of a warm-based ice sheet system far larger, and more integrated, than the extant VFF of the region.
The occurrence of col-like passes occupied by VFF, and both morphologically complementary VFF-free passes and difluent/transluent outlets, are consistent with widespread linear erosion by formerly thicker and more extensive glaciers. This suggests that a period of signiicant deglaciation has occurred, following a regional glacial maximum (Baker and Head, 2015;cf. Brough et al., 2016). The distribution and elevations of channel heads is especially consistent with the wastage of montane LVF > 50 but not with LDA > 50. However, the LDA > 50 are generally lower in elevation than LVF > 50 (Appendix S2.2) and, based on adiabatic atmospheric warming, might be predicted to produce larger amounts of meltwater than LVF > 50. This apparently contradictory observation needs to be investigated further. A starting point lies in a comparison with glaciers on Earth that have termini spanning a restricted altitudinal range (i.e. hypsometrically bottomheavy). These should be the most sensitive population to any change in ELA, due to adiabatic warming at lower elevations (Furbish and Andrews, 1984). In this respect, on Mars, on the basis of elevation, LDA (median MGE −2815 m; Appendix S2.2) might be predicted to be more sensitive than LVF (median MGE −2469; Appendix S2.2). However, the dry adiabatic lapse rate of the martian atmosphere is only 4.3 K.km −1 (Leovy, 2001), in contrast to a value of ~10 K.km −1 for Earth. Moreover, atmospheric dust on Mars reduces the effective adiabatic lapse rate to 2.5 K.km −1 (Leovy, 2001). Fastook et al. (2008) determined that the absence of an effective adiabatic lapse rate on Mars means that sublimation diminishes less with increasing altitude than on Earth. Consequently, higher glaciers on Mars would be expected to have a reduced mass balance, due to persistent sublimation (Fastook et al., 2008). In conformity with this prediction, that mass balance on Mars should have an inverse relationship with elevation (positive mass balance at lower elevations but reducing, possibly becoming negative, at higher elevations), Phlegra Montes is characterised by bottom-heavy glacial systems (Appendix S2.2) that are overlooked by relict glacial landscapes at higher elevations. However, evidence of warm-based glaciation is common in these relict landscapes. This suggests that a warm-based phase of glaciation was succeeded by VFF wastage, particularly at higher elevations, dominated over time by sublimation. Other factors contributing to variation in VFF sensitivity could involve altitude variation in moisture-bearing wind direction, and both aspect and insolation-shadowing, particularly of generally low-lying LDA by large massifs. All these factors could contribute to the evolution of subregional climate and microclimate in Phlegra Montes but need to be investigated further. SHARAD results (Fig. 3), and complementary analysis of HiRISE imaging data (Appendix S2.3), indicate that LDA 1694 is an accumulation of nearly pure ice, up to ~390 m thick (average thickness ~ 200 m; Appendix S2.2), covered by a regolith layer (a rock glacier with up to 10% low density regolith by volume or 4% basaltic rock are other possible conigurations). The regolith could have reduced the sensitivity of LDA 1694 to surface melt suficiently for it to persist through multiple obliquity cycles (Bramson et al., 2017). Unfortunately, the radar data (unaffected by topographic clutter) required to determine the debris cover and content (beyond probably thin suricial deposits resolvable in HiRISE imaging data) of LVF in Phlegra Montes are not yet available.

The possibility of subglacial erosion by dry based glaciers
Warm ice with basal meltwater deforms more easily and moves more rapidly than cold ice (without basal meltwater) due to basal sliding. In addition, warm ice entrains more debris and produces erosive meltwater discharges, due to pressure luctuations induced by variable ice coninement. It is, therefore, more erosive than cold ice. Although cold based glaciers can achieve limited erosion (Cuffey et al., 2000;Hambrey and Fitzsimons, 2010), basal sliding owing to the presence of basal meltwater is required for glaciers to be regionally erosive. However, if suficient dynamic basal pressure and/or strain heating is achieved at the base of very cold glaciers, melting can occur, giving rise to polythermal glaciers. These systems are typically frozen and dry-based around their thin margins but warm-based in their interiors or where they impinge on large bedrock obstacles. Polythermal glaciers move by basal sliding or, where their substrates have high porewater pressure, by subglacial deformation. Their drainage systems are mainly supraglacial and englacial, more rarely subglacial (Hambrey and Glasser, 2012). However, meltwater can penetrate cold-based margins at the base (Kavanaugh and Clarke 2001;Bingham et al. 2005 and2006).
In Phlegra Montes, the macro-scale morphological evidence suggests that subglacial abrasion occurred across the region; cols, passes and crater rim breaches are overlooked by rounded, linearly abraded domelike peaks and whalebacks. Sugden (1974) interpreted morphologically analogous landscapes on Earth as being indicative of basal sliding in limited contexts where local conditions allowed basal pressure melting. Although cold-based glaciers in Antarctica have been observed to produce centimetre-scale erosional landforms (Atkins et al., 2002;Bockheim, 2010;Atkins and Dickinson, 2007), the kilometre-scale abrasion landforms and channels in Phlegra Montes dwarf them, making coldbased erosion alone an unlikely explanation of the landscape of Phlegra Montes. Hence, if the extant VFF of Phlegra Montes are presently cold-based, the abrasion landforms and channels visible beyond the VFF margins most likely represent a phase of expanded warm-based or polythermal glaciation.

Summary and conclusions
A total of >16,000 km 2 of VFF exceeding 50 km 2 have been mapped in Phlegra Montes (Fig. 1, Appendix S1). SHARAD data (and area-based thickness estimates in Appendix S2.2), indicate that LDA 1694 is the thickest regional VFF, at ~390 m (average thickness 200 m). Exposed former sub-glacial erosion surfaces and landforms morphologically and contextually analogous to landforms considered diagnostic of warm based glacial erosion on Earth (Glasser and Bennett, 2004) are common in the region. These landforms relect a process continuum bracketed by wet-based areal scour and subglacial channel incision as endmembers. Intermediate processes involved linear abrasion, s-form erosion, groove incision and the development of distal channels from grooves. The extensive presence in Phlegra Montes of landforms considered indicative of warm based subglacial processes on Earth is, therefore, interpreted as evidence that subglacial abrasion occurred widely in Phlegra Montes. The interpretation of glacial wastage following warm based abrasion, supported by both erosional and depositional markers along the walls of troughs occupied by LVF, and both along and across abraded rock forms now exposed, means that the present VFF of Phlegra Montes do not represent the maximum extent of glaciation in the region. The incision of channels has been common in Phlegra Montes; both their geomorphology and context support the hypothesis that they formed in subglacial, ice contact and proximal ice marginal settings, prior to VFF retreat. Gallagher and Balme (2015) concluded that the presence in Phlegra Montes of an esker still physically associated with LVF 393, and its contextual landsystem, represents spatially limited warm-based conditions, associated with enhanced geothermal heating along a narrow but regionally extensive graben. Butcher et al. (2017) conirmed that the potential for subglacial melt induced solely by atmospheric heat, in the absence of suficient geothermal heat lux and/or strain heating, is very unlikely. The Phlegra Montes esker remains an extremely rare landform but this paper has presented evidence of subglacial abrasion and channel incision consistent with the widespread action of warm-based glaciers in the region. Consequently, a better understanding of the thermal environment, both endogenic and exogenic, in which these glaciers operated is required to reconcile the landscape with the palaeoenvironmental contexts it represents. An environmental transition is strongly expressed by the regional geomorphology; a warm-based phase of glaciation, in which subglacial meltwater was produced and involved in distinctive landscape generation, was succeeded by VFF wastage, particularly at higher elevations and, therefore, probably dominated by sublimation (cf. Fastook et al., 2008).
Some individual landforms or landscapes that appear to be promising candidates for warm-based glacial erosion could be interpreted in a non-glacial paradigm. However, any such interpretation would ignore the total assemblage of landforms present in Phlegra Montes. Contextualized by extant VFF, and by reference to diagnostic morphological analogues on Earth (Glasser and Bennett, 2004), the landform assemblage of Phlegra Montes can be interpreted as a complete and compelling set of morphological indicators of warm based glaciation.

Funding acknowledgement
FEGB is part of the PALGLAC research team supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 787263). FEGB also acknowledges funding from the Science and Technology Facilities Council (STFC) grant ST/N50421X/1.

Declaration of Competing Interest
None.