Recent basal melting of a mid-latitude glacier on Mars Journal Item

Evidence for past basal melting of young (late Amazonian-aged), debris-covered glaciers in Mars’ mid-latitudes is extremely rare. Thus, it is widely thought that these viscous flow features (VFFs) have been perennially frozen to their beds. We identify an instance of recent, localized wet-based mid-latitude glaciation, evidenced by a candidate esker emerging from a VFF in a tectonic rift in Tempe Terra. Eskers are sedimentary ridges deposited in ice-walled meltwater conduits and are indicative of glacial melting. We compare the candidate esker to terrestrial analogues, present a geomorphic map of landforms in the rift, and develop a landsystem model to explain their formation. We propose that the candidate esker formed during a transient phase of wet-based glaciation. We then consider the similarity between the geologic setting of the new candidate esker and that of the only other candidate esker to be identified in association with an existingmid-latitude VFF; both are within tectonic graben/rifts proximal to volcanic provinces. Finally, we calculate potential basal temperatures for a range of VFF thicknesses, driving stresses, mean annual surface temperatures, and geothermal heat fluxes, which unlike previous studies, include the possible role of internal strain heating. Strain heating can form an important additional heat source, especially in flow convergence zones, or where ice is warmer due to elevated surface temperatures or geothermal heat flux. Elevated geothermal heat flux within rifts, perhaps combined with locally-elevated strain heating, may have permitted wet-based glaciation during the late Amazonian, when cold climates precluded more extensive wet-based glaciation on Mars. Plain Language Summary It is widely thought that debris-covered water ice glaciers in Mars’ mid-latitudes have been permanently frozen under cold Martian climates and have not produced liquid water. Using satellite images, we have identified a rare, localized example of past (approximately 110Myr ago) melting of a mid-latitude glacier on Mars. We identify an “esker” emerging from the glacier terminus. Eskers are ridges of sediment deposited bymeltwater flowing through tunnelswithin glaciers. The newglacier-linked esker is located within a tectonic rift. The only other esker to be identified in association with a mid-latitude glacier on Mars to date is located within a similar tectonic trench. On Earth, tectonic rifts are commonly associated with enhanced geothermal heat—heat provided from the interior of the planet. We explore whether enhanced geothermal heat in or near the glaciated rift could have melted the bed of the glacier despite cold recent climates. We calculate theoretical temperatures beneath a glacier of varying thickness under different surface temperature and geothermal heat conditions. We find that when combined with heat generated by internal deformation of flowing ice, geothermal heat could have permitted localized cases of glacial melting while cold climates prevented more extensive melting of glaciers throughout Mars’ mid-latitudes.

To our knowledge, one other example of an esker associated with a late Amazonian-aged (~150 Ma) midlatitude VFF has been documented to date, in a graben in Phlegra Montes, NE of the Elysium volcanic province (Gallagher & Balme, 2015). These young, mid-latitude, VFF-linked candidate eskers are distinct from other candidate eskers on Mars (see section 4.5) in their age, mid-latitude location, and association with an extant parent glacier.
In this contribution, we present a geomorphic map of landforms in the Tempe Terra rift and develop a landsystem model to explain their formation. We analyze the meter-scale three-dimensional morphology of the candidate esker and compare it to terrestrial analogues. We highlight that the rift setting of the newly identified candidate esker in Tempe Terra is similar to the graben in which the Phlegra Montes candidate VFF-linked esker is located (Figure 1c) (Gallagher & Balme, 2015). Terrestrial graben/rifts are commonly associated with elevated geothermal heat flux (e.g., Fernàndez & Banda, 1990;Jordan et al., 2010). Additionally, strong glacier flow convergence in these steep-sided troughs could drive significant viscous strain heating of glacial ice. That both of the candidate mid-latitude VFF-linked eskers found to date occur within graben/rifts suggest a genetic link between the setting and the landform. The lack of identified evidence for basal melting of modern VFFs elsewhere in Mars' mid-latitudes implies that cyclic (~120 Myr) obliquity-induced climate changes (e.g., Madeleine et al., 2009;Schorghofer, 2008) alone were insufficient to permit widespread basal melting of VFF without an additional heat source. Thus, we explore the possibility that late Amazonian magmato-tectonic processes (e.g., Hauber et al., 2011;Neukum et al., 2004a;Tanaka, 1990), could have induced positive geothermal anomalies in the vicinity of these ancient rifts (Hauber et al., 2010) that possibly supplemented by viscous strain heating within the VFF, facilitated glacial melt (Gallagher & Balme, 2015). We undertake some simple calculations of possible basal temperatures for a range of possible VFF thicknesses, mean annual surface temperatures, geothermal heat fluxes, and driving stresses, which unlike previous studies (e.g., Cassanelli et al., 2015) include the possible role of internal strain heating.

Observations and Mapping
We mapped the landforms in the rift (Figures 3 and S1 in the supporting information) by analyzing~6 m/pixel Context Camera (CTX) images (Malin et al., 2007) and~100 m/pixel day and night infrared Thermal Emission Imaging System (THEMIS) mosaics (Christensen et al., 2004;Edwards et al., 2011) (see Table S1 in the supporting information for list of data products). We digitized the map in Environmental Systems Research Institute ArcGIS 10.1 using a sinusoidal projection centered on the long axis of the N-S oriented rift (83.5°W) to minimize cartographic distortion.
Following the method of Kirk et al. (2008), we generated a 1 m/pixel digital elevation model (DEM) and 25 cm/pixel orthorectified images of the sinuous ridge from two 25 cm/pixel High Resolution Imaging Science Experiment (HiRISE) images (Table S1) (McEwen et al., 2007) in U.S. Geological Survey Integrated Software for Images and Spectrometers 3.0 and BAE Systems SOCET SET 5.6.0.

Estimation of Impact Crater Retention Age of VFFs
We measured the size-frequency distribution of impact craters on the surface of the VFFs within the rift to determine the minimum age of the termination of the last major VFF-surface modification phase. We measured impact craters using the CraterTools (Kneissl, van Gasselt, & Neukum, 2011) add-in for ArcGIS, according to the method of Gallagher and Balme (2015). Complete CTX image coverage permitted inclusion of all undeformed crater-like circular features of diameter (D) ≥ 50 m. Pitting and lineation prevented reliable identification of impact craters of D < 50 m. We estimated the crater retention age in CraterStats 2.0 (Michael & Neukum, 2010) based on the Ivanov (2001) production function and Hartmann and Neukum (2001) chronology function. We were mindful of likely underestimation of the true age due to downwasting of VFF surfaces and flow deformation of impact craters. We acknowledge the potential factor of 10 uncertainty in either direction of the derived age arising from reliance upon small D < 250 m impact craters (Hartmann, 2005), and the uncertainties arising from a small count area (1.51 × 10 3 km 2 ) and relatively small sample size (n = 420) (Warner et al., 2015). Thus, the model crater-retention ages returned by our analyses are very approximate. We assessed the likelihood of the temperature at the base of the VFF reaching the pressure melting point by comparing the potential basal heat sources (geothermal heating, plus viscous strain heating within the VFF) with the likely rate of heat loss to the surface through the body of the VFF. We take a conservative approach and do not include the effect of debris cover of the glacial ice within our calculations. This reflects the  (Table S1) (a) central portion of rift (extent in Figure 1b, see Figure S1 for full map) showing locations of (b) sinuous ridge, (c) LDA-proximal transverse ridges, (d) highland mantle, (e) LDA-proximal floor, (f) LDA-distal knobbly ridge (the white box shows extent of THEMIS night IR inset, in which bright areas have higher surface temperatures relative to darker areas), (g) central rift floor, and (h) polygonized wedge (the dashed line delineates polygonized mounds).
Journal of Geophysical Research: Planets 10.1002/2017JE005434 possibility that ubiquitous debris cover was deposited after the termination of the major phase of glacial ice deposition (e.g., Baker & Head, 2015;Fastook et al., 2014). Debris cover of any thickness would reduce the thermal gradient and rate of heat advection toward the surface from the glacier bed (Benn & Evans, 2010), and thus increase the basal temperature. Therefore, the calculated basal temperatures assuming nondebris-covered ice are likely to be minimum temperatures for the given parameters. Terrestrial glaciers often show nonlinear temperature/depth profiles (Cuffey & Paterson, 2010), but the very slow flow of VFFs on Mars and the lack of a strong adiabatic lapse rate (such that any ice advected downglacier from higher altitudes does not distort the vertical temperature profile toward lower temperatures) suggest that like Cassanelli et al. (2015), we can assume a linear temperature profile within the VFF in the first instance. We can then calculate the temperature difference (dT, K) between the bed and surface of the VFF: where H is ice thickness (m), k T is thermal conductivity of the ice (assumed to be 2.5 W m À1 K À1 ) (Cassanelli et al., 2015), and Q is basal heating (W m À2 ), equal to geothermal heat flux, G, plus any internal strain heating. Volumetric strain heating (P, W m À3 ) at depth can be calculated as (e.g., Lüthi et al., 2015) where A T is the temperature-dependent flow rate factor for ice (s À1 Pa À3 ) and σ xz is shear stress, calculated as ρ g h tan α, where ρ is ice density, g is gravity, h is the depth within the VFF, and α is ice surface slope. Given the strong nonlinearity of P with depth within glaciers, it has been assumed that all of the strain heating over a column of ice can affect the base (e.g., Hindmarsh, 1990), but we conservatively assume that only the heating in the lowest 10% of the ice thickness will affect the basal temperature and integrate P over this depth to give basal strain heating in W m À2 .
The key unknowns for estimating basal temperature are geothermal heat flux, surface temperature, ice thickness, and surface slope. Following Cassanelli et al. (2015), we estimate basal temperatures over a range of possible geothermal heat flux values (10 to 100 m W m À2 ), mean annual surface temperatures (190-230 K) and ice thicknesses (500 to 1,500 m). The minimum ice thickness corresponds to our estimate for the maximum thickness of existing VFF within the rift (e.g., northern VFF in Figure 4), while the maximum ice thickness corresponds to the depth of the rift relative to its northern and southern margins. We use a range of plausible ice surface slopes (from~0.1 to 3.5°) based on observations of Martian VFFs (including the VFF discussed in this study) with the range of ice thicknesses discussed in order to generate a range of driving (shear) stresses (σ xz, from 10 to 100 kPa) covering the range of values reported for Martian VFFs (e.g., Karlsson et al., 2015), and terrestrial ice masses. The surface slopes of the existing VFFs within the rift are approximately 2-2.5°.
To calculate the basal temperature for each set of values, we first calculate basal temperature neglecting strain heating, using equation (1). We use this initial temperature to estimate A T using the values in Cuffey and Paterson (2010, p. 75) and use this value in equation (2) to estimate P in the lowest 10% of the ice. We add this to the geothermal heat to calculate a new basal temperature from equation (1), and then a new value of A T . We iterate this process until the calculated basal temperature changes by less than 0.1 K or reaches the pressure melting point (assumed to be 273 K). Except for the very warmest, thickest, and highest driving stress ice, this process takes fewer than five iterations. It should be noted that A T as given by Cuffey and Paterson (2010) is applicable to clean terrestrial ice. There remains significant uncertainty over both the nature and quantity of impurities (e.g., salts, dust, and debris) within Martian VFFs, and whether they have a net softening or hardening effect upon VFF rheology relative to clean ice (e.g., Parsons et al., 2011). Thus, given present uncertainties, we consider A T as described by Cuffey and Paterson (2010) to be an appropriate middle-ground for current estimations of VFF rheology, but we do consider the effect of variations in A T upon our model results.

Observations
The Tempe Terra rift is~110 km long,~24.5 km wide, and up to~2,300 m deep and is bounded to the east and west by normal fault scarps ( Figure 4) (see also Hauber et al., 2010 Figure 2b). Its northern and southern margins rise~1,550 m above the rift floor, forming a closed basin (Figure 4). Pitted and lineated lobes, consistent with the lobate debris apron (LDA) subtype of Martian mid-latitude VFFs (e.g., Head et al., 2010; Journal of Geophysical Research: Planets 10.1002/2017JE005434 Squyres, 1979), extend~15 km down into its center from all marginal walls (Figures 3a-3c and 4). They converge laterally and their termini confine a region of~225 km 2 of the rift floor. Their upper margins are gradational with highland mantle that drapes topography and infills topographic lows (Figures 3d and 4), similar to that described by van Gasselt et al. (2011).

Sinuous Ridge
An S-N oriented sinuous ridge (Figures 2, 3b and 4-6) emerges from the terminus of the southern LDA at 46.17°N, 83.06°W. It originates~8 km into the LDA tongue and extends a further 7 km north into the foreland. It has a sinuosity (ratio of crest length to shortest end-to-end distance) of 1.12. We identify four morphological subzones along the ridge (Figures 5 and 6a).  (Table S1) (Neukum et al., 2004b;Jaumann et al., 2007). Depths and profiles of subsurface contacts are estimates. The sinuous ridge originates as an 80-150 m wide textural band in the LDA surface (subzone I, 0-5 km), 700 m from its lateral margin. It transitions between 60 and 160 m wide ridges and troughs (sub-DEM resolution) (Figure 6b). The ridge crest becomes topographically distinct from the LDA surface toward the end of this subzone ( Figure 5).  (Table S1) (a) CTX image (Table S1) of sinuous ridge (the white arrows show textured band) showing locations of (b) ridge and trough within textural band (subzone I), (c) multi-crested ridge atop LDA-mantled topographic rise (subzone II), (d) steep-sided flank troughs (subzone II), (e) longitudinal flank fractures (subzone III), (f) pits and polygons (subzone III), (g) polygonally fractured impact crater (subzone III), (h) transverse ridge (white arrows) that crosses the sinuous ridge (subzone III, track shown by white dotted line in Figure 6a), (i) pits on ridge crest (subzone IV), (j) tributary A (white arrows) with textured band, and (k) tributary B (white arrows) with smooth depression at head.

10.1002/2017JE005434
In subzone II (5-7.6 km), the ridge transitions into an~150 m wide,~5-10 m high multi-crested ridge with two sharp, parallel crests, atop a broader,~500-1.5 km wide,~30 m high topographic rise (Figures 5 and  6c). The central~500-800 m wide portion of the topographic rise is exposed above the LDA surface. Its margins are dissected by several steep-sided troughs (~20 m wide,~25-150 m long, and~3-5 m deep; Figure 6d). The flanks of the topographic rise in subzone II are superposed by LDA materials. This limits scrutiny of whether this topographic rise represents an emergent undulation of the bed upon which the central multi-crested ridge is superposed, or a lateral continuation of the ridge that has been incompletely exposed from the basal ice. We favor the latter on the basis of its elongation along the axis of the ridge ( Figure 5). Therefore, draping of the ridge flanks by LDA materials means that our estimate of total ridge height in this subzone (~40 m) is likely an underestimate of the true height.
As it emerges from the LDA terminus, the main ridge widens into a broad (up to 1,300 m wide),~10-55 m high diamond-shaped plateau (subzone III, 7.6-10.7 km; Figures 5 and 6a). In subzone III, the flanks of the ridge are fractured subparallel to its long axis (Figure 6e). Its surface hosts shallow~10 m diameter pits and~15-30 m diameter polygonal fractures (Figure 6f), which cut through impact craters (Figure 6g). At the northernmost end of subzone III, the ridge is superposed by a complex of~5 m high,~30-100 m wide ridges oriented transverse to the sinuous ridge (Figures 6a and 6h). The transverse ridge complex originates at the terminus of the LDA~1 km east of the sinuous ridge and crosses its crest in subzone III, before tracking SW along its NW facing margin (Figure 6a).
At the northernmost end of subzone III, 2 km beyond the LDA terminus, the sinuous ridge transitions to a narrow (150-200 m wide),~5-15 m high, sharp-crested ridge (subzone IV, 10.7-15 km; Figure 5). Pits are densely clustered along its crest ( Figure 6i). The ridge gradually fades into the LDA-proximal floor at the northernmost end of subzone IV.

Rift Floor Units and Landforms
At its emergence from the LDA terminus, the sinuous ridge crosses into the LDA-proximal floor unit (Figure 3e), which forms the LDA foreland. The LDA-proximal floor comprises smoothed, pitted, and fractured materials ( Figure 3e). Fractures commonly intersect to form~100 m diameter polygons. This unit extends 2-10 km from the LDA termini toward the center of the rift, separating the LDA termini from an annulus of LDA-distal "knobbly" ridges ( Figure 3f). The LDA-distal knobbly ridges are oriented oblique to the SSW-NNE oriented floor fault scarps, which cut through them, but correlate well with the termini of the surrounding LDAs, being subparallel to them. A distinct population of arcuate LDA-proximal transverse ridges ( Figure 3c (Figure 3h). Its THEMIS IR signature is intermediate between those units. The~100 m diameter polygons have wider (tens of meters) marginal troughs than those within the LDA-proximal floor unit. They superpose a chain of mounds that is morphologically similar to the parallel LDA-distal knobbly ridge to the south.

10.1002/2017JE005434
owing to the small sample size (n = 5), we do not treat this result as significant. Unavoidable dependence upon small impact craters could result in a factor of 10 uncertainty in these model ages (Hartmann, 2005). Age estimates at the upper extreme of this uncertainty are~300 Ma for D(50-80 m) impact craters, and 1.1 Ga for D(130-150 m) impact craters.

Origin of the Sinuous Ridge
Given the glaciated tectonic setting of the sinuous ridge in Tempe Terra, we consider four possible origins: (i) an esker, (ii) an ice-confined subglacial lava flow, (iii) a medial moraine ridge, or (iv) an inverted fluvial channel, later buried and partially reexposed by LDA emplacement and degradation. The production of liquid water by Amazonian melt of the LDA is implicit in the first two hypotheses; channelized meltwater drainage is a prerequisite for formation of eskers and esker-like lava flows (e.g., Mathews, 1958), and meltwater is an inevitable by-product of subglacial volcanism. Meltwater is not required for medial moraine formation via supraglacial accumulation of rockfall debris. Liquid water is required for hypothesis (iv), but in this case water could have occurred much earlier in Mars' history, divorced temporally from the glacial landscape described here.

Esker Hypothesis
The sinuous ridge emerges from the terminus of an LDA, which leads us to consider a genetic link between these landforms. The sinuous ridge is similar in mean length and sinuosity to terrestrial esker systems in Canada (5,932 measured, mean length 15.6 km, and mean sinuosity 1.08) (Storrar et al., 2014). The longitudinal transitions in crest morphology between subzones I and IV are remarkably similar, in both morphometry and sequence, to those observed along North American eskers that formed within pressurized water-filled subglacial conduits (Perkins et al., 2016;Shreve, 1985). Sharp-to round-crested elements of these eskers (analogous to subzone IV, 5-15 m high, 150-200 m wide) are typically~3-50 m high and~150 m wide (e.g., Perkins et al., 2016 ;Shreve, 1985). Broad-crested elements of the Katahdin esker, Maine, USA (analogous to subzone III, 10-55 m high, 1,300 m wide) are typically~10 m high and 2 km wide (Shreve, 1985). Multicrested reaches of North American eskers (analogous to subzone II, 5-10 m high, 150 m wide, atop~30 m high, 500-1.5 km wide topographic rise) have typical heights of 10-15 m (Perkins et al., 2016;Shreve, 1985), and widths ranging from hundreds of meters (Perkins et al., 2016) to kilometers (Shreve, 1985).
The sequence of morphological transitions between subzones I and IV of the sinuous ridge is the same as those identified along the Katahdin esker system by Shreve (1985). Shreve attributed these transitions to changes in the shape of the meltwater conduit due to bed slope and ice pressure gradient-controlled modification of melt dynamics at the conduit roofs and walls. Under this model, sharp-crested elements form on flat or downslope reaches and broad-crested elements form on upslope reaches. Multi-crested reaches occur on intermediate positive bed slopes between sharp-and broad-crested elements. Obscuration of the LDA bed by existing LDA deposits prevents detailed tests for such topographic relationships; however, the broad-crested subzone (III) does cross upslope reaches of the LDA-proximal floor and transitions back to a sharp-crested morphology downslope ( Figure 5), as predicted under the Shreve (1985) model. As is also predicted by the Shreve (1985) model, the sharp-crested ridge in subzone IV progressively reduces in height as it ascends a gentle upslope reach of the bed at its VFF-distal end ( Figure 5).
Under the esker hypothesis, an alternative explanation for the broad morphology in subzone III is that it represents an outwash fan, which was deposited over a central esker ridge following retreat of the LDA terminus to its present position. Such outwash fans overlie the central ridge of the Vars-Winchester esker in Canada (Cummings et al., 2011). However, outwash fans typically form subaqueously in proglacial lakes, for which we identify no strong evidence within the rift. Continuous meltwater production throughout retreat of the LDA terminus to its present position would be expected to form a continuous sequence of onlapping fans (Cummings et al., 2011) from the northernmost margin of the sinuous ridge to the present LDA terminus, but this is not observed. Therefore, fan formation at this position along the sinuous ridge would require a second discrete phase of meltwater production to that which formed the sinuous ridge in subzone IV. We consider this to be less likely than formation of all subzones during the same phase of meltwater production. Finally, the margins of the ridge in subzone III progressively narrow toward subzone IV. This is more readily explained by a narrowing subglacial meltwater conduit toward subzone IV than by narrowing of an outwash fan in its ice-distal reaches.

10.1002/2017JE005434
Tributary ridges merging with the main ridge are consistent with capture of minor subglacial drainage by a major meltwater conduit (Shreve, 1972). Disruption of the LDA surface at the heads of the tributaries, and troughs in subzone I of the main ridge, is consistent with collapse of the LDA surface under conduit enlargement and evacuation (see Benn et al., 2017). The textural band at the head of tributary A could indicate an underlying (possibly collapsed or infilled) subglacial conduit. This conduit may have formed by lateral escape of overpressurized meltwater from the central conduit, a process that drives formation of anatomozing branches along terrestrial eskers (Gorrell & Shaw, 1991). The broad depression at the head of tributary B is more consistent with capture of water from a subglacial pond or cavity.
Several smaller-scale features of the sinuous ridge can also be readily explained under the esker hypothesis. Pits in its surface can be explained by collapse into void spaces evacuated by sublimation of ice that was trapped during esker deposition. Subsequent exposure from the LDA could have triggered destabilization and sublimation of ice to the atmosphere. Similarly, steep-sided troughs (subzone II) could be explained by sublimation of dead ice that was trapped within the esker surface and subsequently sublimated upon exposure from the LDA.
Polygonal fractures crosscutting impact craters superposed upon subzone III of the sinuous ridge provide strong evidence that fractures are not primary structures associated with ridge formation but formed due to secondary modification of its surface. They can be explained by thermal contraction cracking of volatilerich glaciofluvial sediments comprising an esker. Their size and morphology are consistent with contraction crack polygons observed in ice-rich terrains elsewhere on Mars and on Earth (e.g., Mangold et al., 2004;Marchant & Head, 2007). Longitudinal fractures along the flanks of the ridge in subzone III can be explained by relaxation of the ridge flanks following loss of lateral support by confining walls of a meltwater conduit (e.g., Burke et al., 2008;Burke et al., 2012).
We interpret the small transverse ridge complex that crosses the crest of the sinuous ridge at the northern end of subzone III (Figures 6a and 6h) to be a moraine ridge formed by transient ice margin fluctuation following exposure of subzone IV by LDA retreat, similar to <5 m high transverse ridges crossing eskers on the floors of Spitsbergen fjords (Dowdeswell & Ottesen, 2016).

Ice-Confined Subglacial Lava Flow Hypothesis
Given its glaciated tectonic setting, we consider whether an ice-confined subglacial lava flow could form a viable alternative origin for the sinuous ridge. We find that our observations are inconsistent with such an origin. Terrestrial examples of "esker-like" lava flow ridges have been associated with eruptions of mafic lava into subglacial environments (e.g., Hungerford et al., 2014;Lescinsky & Fink, 2000;Mathews, 1958;Smellie & Skilling, 1994). Their esker-like morphologies have been attributed to exploitation of preexisting subglacial meltwater channels by lava flows. Thus, they may adopt very similar planforms to eskers.
The surfaces of subglacial lava flows are typically characterized by fractures forming columnar joints arising from rapid quenching of lava flow margins upon contact with ice (Lescinsky & Fink, 2000). While polygonal fractures are observed on the surface of the sinuous ridge in Tempe Terra, they have diameters 1 to 2 orders of magnitude greater than those within terrestrial ice-contact lava flows (tens of centimeter to 1 m) (e.g., Hungerford et al., 2014;Lescinsky & Fink, 2000). Furthermore, we assert that the observed polygonal fractures cannot be explained by rapid cooling because they crosscut superposed impact craters, which postdate ridge formation (Figure 6g). Thus, fracture formation was not coincident with ridge emplacement but occurred due to secondary modification of ridge surface materials. Thermal contraction of ice-rich sediments under the esker hypothesis can more adequately explain both the scale and relative time of formation of the polygonal fractures.
Furthermore, there exists insufficient evidence to invoke a subglacial volcanic eruption to explain the origin of the sinuous ridge. While it is conceivable that a subglacial volcanic eruption could occur in such an extensional tectonic setting, no large-scale disruption of the LDA surface indicative of an underlying eruptive edifice, such as an ice cauldron, is observed, and no lava flows or volcanic edifices are identified in the surrounding region. Moreover, origination of tributary B from a distinct point to the main ridge and tributary A is better explained by scavenging of water over wide areas of the glacier bed than by the origin of tributary B from a distinct volcanic edifice. Thus, we argue that a subglacial volcanic eruption is inconsistent with the geomorphic evidence documented here.

Medial Moraine Hypothesis
Medial moraines are bands of debris oriented parallel to glacier flow direction and positioned distal from the lateral margins of a glacier. The insulating effect of supraglacial medial debris bands can promote differential ablation of the glacier surface, resulting in preservation of a debris-covered ice core as a ridge within the ablation zone, and occasionally the proglacial zone. These characteristics lead us to consider a medial moraine as a third alternative explanation for the origin of the sinuous ridge.
Whereas the morphological transitions along the sinuous ridge can be readily reconciled with a subglacial esker origin, we find that they are harder to explain under a medial moraine origin. Medial moraines typically have a relatively simple planform morphology. Where they can be traced to their source, medial moraines originate as broad debris accumulations, becoming concentrated into narrow, low-sinuosity ridges down-glacier (e.g., Eyles & Rogerson, 1978;Goodsell et al., 2005;Vere & Benn, 1989). They rarely undergo abrupt variations in width but may spread laterally toward the glacier terminus. Lateral expansion of medial moraines is typically accompanied by an overall lowering of ridge height (e.g., Eyles & Rogerson, 1978;Vere & Benn, 1989). This contrasts with the abrupt variations in width observed along the sinuous ridge, and the apparent increase in height that accompanies ridge widenings in some locations.
Medial moraines can form via "ice stream interaction," "avalanche-type," or "ablation-dominant" mechanisms, as described by Eyles and Rogerson (1978). "Ice stream interaction-type" medial moraines form by the convergence of two lateral moraines at a glacier confluence (Eyles & Rogerson, 1978). We exclude such an origin for the sinuous ridge, since the textural band in subzone I is not spatially associated with a longitudinal confluence between VFF lobes. The sinuous ridge is clearly superposed by the LDA and emerges from beneath LDA materials in subzones I and II. Based on this superposition relationship, we can also exclude an origin as an avalanche-type medial moraine: rare, often transient bands of debris transported supraglacially away from a subaerial debris source such as a rockfall or avalanche at a headwall or nunatak (Eyles & Rogerson, 1978).
Ablation-dominant (AD) medial moraine formation involves entrainment of subaerially derived debris into englacial transport paths and its reexposure downglacier by ablation of overlying ice (Eyles & Rogerson, 1978). Eyles and Rogerson (1978) identified three subtypes (AD1, AD2, and AD3) based on the mechanism of debris entrainment, which we now consider in sequence. Debris comprising AD1 medial moraines is entrained into shallow englacial transport paths via crevasses below the firn line (Eyles & Rogerson, 1978). A lack of visible crevasses or crevasse traces on the VFF surface leads us to exclude this mechanism for formation of the sinuous ridge.
Debris comprising AD2 medial moraines is entrained into englacial transport paths via burial by snowfall or avalanches above the firn line (Eyles & Rogerson, 1978). Ice flow patterns within terrestrial glaciers commonly dictate that supraglacial debris deposited at the foot of the headwall is transported along englacial flow paths that are close to the bed (Benn & Evans, 2010). During an earlier phase of glaciation within the rift, the bedrock outcrops to the south east of the sinuous ridge were likely overridden by the LDA to form subglacial bedrock knobs, as evidenced by their rounded, planed, and smoothed morphologies. This would likely have resulted in direct ice flow from the southern rift wall to the LDA terminus. Under this configuration, the bedrock knobs could have formed a topographic obstruction to near-bed englacial transport of a debris septum sourced at the southern rift wall. If englacial debris had been transported along a flow path above the bedrock knobs, evidence for let-down of a moraine ridge over the present-day bedrock outcrops would also be expected, but is not observed. Thus, we find insufficient evidence that the sinuous ridge is an AD2 medial moraine sourced from the southern rift wall.
The AD3 formation mechanism for medial moraines is rarer than the AD1 and AD2 subtypes. Under this mechanism, debris is derived directly from subglacial bedrock knobs (Eyles & Rogerson, 1978;Vere & Benn, 1989). As discussed above, the bedrock outcrops to the south east of the sinuous ridge may have formed subglacial bedrock knobs during an earlier stage of glaciation. In the absence of data for the bed topography of the LDA, including the continuation of the bedrock outcrop beneath the LDA, we are unable to exclude conclusively the possibility that these bedrock knobs could have provided a debris source to the sinuous ridge. However, the orientation of the sinuous ridge, lateral to the nearest bedrock outcrop~700 m to the south east, seems inconsistent with debris supply from this location. Meltwater is typically required as an agent for erosion and entrainment of such large volumes of subglacial bedrock material into glacial ice.

Inverted Fluvial Channel Hypothesis
It is possible that the sinuous ridge originated as an ancient subaerial fluvial channel, which was topographically inverted by extensive erosion of the surrounding rift floor and subsequently buried by and reexposed from, later LDA deposits. However, the sinuous ridge does not follow the steepest topographic gradient, and even ascends~5-10 m high topographic undulations in some reaches ( Figure 5). This is inconsistent with subaerial gravity-driven flow, but it can be explained by flow under hydraulic pressure in a subglacial meltwater conduit according to the esker hypothesis (Shreve, 1972). Second, we identify no other evidence for subaerial fluvial erosion within the rift, or inversion of topography in the broader region. We exclude an inverted floodplain for the origin of the broad zone in subzone III. The upper surface of subzone III is tilted, and its crest is higher in elevation than the adjacent portions of the sinuous ridge, which is inconsistent with an origin as a sediment sink. We also observe no evidence for meanders or scroll bars anywhere along the ridge, as are observed along inverted channels and within inverted floodplains in Aeolis Planum (Burr et al., 2009). Third, evidence for disruption of the LDA surface (subzone I; Figure 6b) is more consistent with ridge formation via a process that actively influenced the structure of the overlying LDA (i.e., opening and subsequent evacuation of a meltwater conduit) than with passive advance and retreat of LDA over a preexisting inverted channel. Finally, the inverted fluvial channel hypothesis invokes multiple unrelated processes: aggradation of sediment in a fluvial channel, armoring of channel deposits, differential erosion of surrounding plains, then burial and subsequent exposure by a cold-based glacier. This contrasts with the simpler process explanation under the esker hypothesis, which invokes a genetic relationship between the sinuous ridge and the existing VFF from which it emerges.

Working Hypothesis
Of the four hypotheses considered here, we find a subglacial esker origin to be most consistent with our observations of the sinuous ridge on the basis of (1) its association with an existing putative debris-covered glacier, (2) crossing of topography and longitudinal transitions in cross-sectional morphology that are strikingly similar to those of terrestrial eskers formed under hydraulic pressure in subglacial meltwater conduits, and (3) meter-scale morphologies indicative of a volatile-rich composition and formation within an iceconfined subglacial environment.

Origin of Pro-LDA Landforms and Units
We interpret the LDA-proximal transverse ridges as moraines based on their similarity to other moraine-like ridges on Mars (e.g., Arfstrom & Hartmann, 2005;Hubbard et al., 2011). The oblique orientation of the LDAdistal knobbly ridges to structural faults on the rift floor is inconsistent with tectonic control. We interpret them as moraines based on their correlation in orientation to adjacent LDA fronts. Their knobbly appearance and muted topographic expression could imply derivation from supraglacial or englacial debris released from a glacier terminus, similar to marginal aprons at the termini of cold-based glaciers in the Antarctic Dry Valleys (Cuffey et al., 2000).
The high thermal inertia signature of the central rift floor is consistent with a consolidated lithology of the surface materials relative to the LDA-proximal floor (Putzig & Mellon, 2007). The coherence in strike of the lineations here with local faults suggests that they are structural features of the faulted bedrock, rather than streamlined glacial landforms. We therefore interpret the central rift floor as an unglaciated surface, possibly comprising bedrock, which also outcrops in the rift-marginal fault scarps (Figure 4). The lower thermal inertia of the LDA-proximal floor is consistent with the presence of deposits of less-consolidated (e.g., glacigenic) material (Putzig & Mellon, 2007). Loss of a volatile component of the material can explain the pitted textures, while thermal contraction cracking of an ice-rich sedimentary substrate (Levy, Head, & Marchant, 2010) can explain the observed polygonal fracturing. We therefore interpret the LDA-proximal floor unit as a deglaciated zone across which the LDAs retreated from the terminal knobbly ridges. The lack of evidence for glacial scour, streamlined forms, or channels implies that the LDAs were cold-based as they crossed this unit. A low sediment supply under such a regime explains the apparent lack of thick sedimentary accumulations in the LDA-proximal floor unit and the muted topographic expression of the LDA-distal knobbly ridges.
We interpret the polygonized wedge as relicts of an older, more extensive glacial advance within the rift because (1) it is superposed by the knobbly ridge complex but superposes a chain of mounds that is morphologically similar to it (Figure 3h), implying stratigraphically-separated episodes of terminal apron deposition; (2) wider polygon-marginal troughs than those within the LDA-proximal floor unit are consistent with more advanced

Crater Retention Ages
Our crater counts ( Figure S2) are in agreement with those by van Gasselt et al. (2011), who identify two dominant crater retention ages in their analyses of LDAs in Tempe Terra. They interpret older crater retention ages (up to~200 Ma) as early traces of modification of LDA surfaces and younger (10-50 Ma) isochrons to represent episodes of emplacement and modification of an ice-rich mantle on top of the LDAs. Our best estimate for the minimum age of the LDAs in the Tempe Terra rift (~110 Ma) places their minimum age of stagnation, and thus formation of the candidate esker, during the Late Amazonian epoch (Hartmann & Neukum, 2001), at a similar time to the candidate VFF-linked candidate esker in Phlegra Montes (~150 Ma) (Gallagher & Balme, 2015). Our upper estimate of 1.1 Ga for the age of the bulk LDA-which considers the potential for a factor of 10 uncertainty (see section 3.3)-is consistent with the 100 Ma-1 Ga age range estimates of Berman et al. (2012) for LDAs in Deuteronilus Mensae and East Hellas. VFFs are generally interpreted to have formed during cyclic (~120 kyr) excursions to intermediate (~35°) planetary obliquity during the late Amazonian (Laskar et al., 2004), which permitted prolonged ice accumulation in Mars' mid-latitudes (Madeleine et al., 2009). , we invoke deposition of an ice-rich mantle deposit superposing the LDAs to explain the younger~30 Ma isochron in our crater counts. This can explain the gradational nature of the contact between the upper margins of the LDAs and the highland mantle unit. The origin of this mantle can most readily be explained by the up to tens of meters-thick (Conway & Balme, 2014) latitude-dependent mantle (LDM) identified throughout Mars' mid-latitude to high-latitude regions (e.g., Mustard, Cooper, & Rifkin, 2001). LDM deposition is generally attributed to more recent (0.4-2.1 Ma) Schon, Head, & Fassett, 2012) increases in planetary obliquity than those during which LDAs formed. We attribute formation of the LDA-proximal transverse ridges to moraine deposition by glacier-like flows within this younger mantling deposit, rather than the bulk LDA because the ridges are very similar to moraine-like ridges associated with 0.1-10 Ma old glacier-like lobes expressed within mantle materials across Mars' mid-latitudes (e.g., Arfstrom & Hartmann, 2005;Hartmann et al., 2014;Milliken, Mustard, & Goldsby, 2003;Souness et al., 2012). Furthermore, differential flow evidenced by lineation in portions of the LDA surface bounded by the LDA-proximal transverse ridges is more easily explained by nonuniform deposition of a surface mantle than by reactivation of isolated lobes within the LDA following termination of bulk flow.

Landsystem Model
The candidate esker is the only landform within the mapped assemblage that is indicative of glacial meltwater production. All other landforms within the rift can be explained without invoking glacial melting. Hence, we propose that the candidate esker likely represents a phase of meltwater production resulting from a transient perturbation of a longer-term cold-based regime. The position of the sinuous ridge relative to terminal moraine deposits (knobbly ridges) and the parent LDA implies that it formed during partial deglaciation of the rift according to the following landsystem model. First, LDAs advanced across the LDA-proximal floor unit toward a period of glacial maximum, merging laterally. They advanced under cold-based regimes, which precluded significant reworking of the bed. Thus, terminal deposits (LDA-distal knobbly ridges), fed primarily by supraglacial debris, were poorly developed, and the proglacial zone (central rift floor) was unmodified. Transient basal melting, perhaps associated with the onset of deglaciation, permitted esker formation beneath the southern LDA. Following esker deposition, the parent LDA returned to a cold-based thermal regime and stagnated~110 Myr ago. The esker was preserved in the basal ice and was partially exhumed from the LDA tongue by sublimation. The upper reaches of the esker remain buried within the LDA and may extend further into the LDA than the visible textural band in subzone I. The highlands surrounding the rift and the LDA were subsequently mantled by a thin ice-rich deposit 30 Myr ago (or periodically deposited and ablated since LDA stagnation, resulting in a model crater retention age of~30 Ma), which underwent viscous deformation, forming flow lobes atop the eastern LDA. These flows mobilized sediment and deposited it in moraine ridges at their termini, forming the LDA-proximal transverse ridges.

Eskers on Mars
To our knowledge, one other example of a late Amazonian-aged (~150 Ma) mid-latitude VFF-linked esker has been documented to date, in a graben in Phlegra Montes, NE of the Elysium volcanic province (Gallagher & Balme, 2015). Sinuous ridges interpreted as eskers in the south circumpolar Dorsa Argentea Formation (Butcher et al., 2016;e.g., Head & Pratt, 2001;Kress & Head, 2015) and southern Argyre Planitia (Banks et al., 2009;Bernhardt et al., 2013) are late Noachian-early Hesperian in age (>3.3 Ga). Their formation has been attributed to basal melting of ancient high-latitude ice sheets under warmer climates (Scanlon et al., 2018) and higher global average geothermal heat flux of the late Noachian and early Hesperian (Fastook et al., 2012). The parent ice sheets of these candidate eskers have long since retreated. Thus, these candidate eskers are distinct from the mid-latitude VFF-linked eskers discussed here in their ancient age, high-latitude locations, and lack of preservation of their parent glaciers.
An esker interpretation of a ridge within late Amazonian aged (~210 Ma) (Kadish et al., 2014) deposits on the lower flank of the low-latitude Arsia Mons volcano (Forget et al., 2006;Shean et al., 2005) has been proposed (Scanlon et al., 2014;Scanlon et al., 2015). The landform assemblage at Arsia Mons, including this candidate esker, is interpreted by Scanlon et al. (2014Scanlon et al. ( , 2015 to have formed under localized wet-based conditions induced by a late Amazonian volcanic eruption beneath a hypothesized ice sheet that has since completely ablated. Unlike the mid-latitude VFF-linked eskers discussed here, the candidate esker at Arsia Mons lacks a preserved parent glacier, although the evidence of a glacial origin for the landforms here is well explored at a broad scale (Forget et al., 2006;Kadish et al., 2014;Scanlon et al., 2014Scanlon et al., , 2015Shean et al., 2005). Absence of a clear association with an existing parent glacier does, however, reduce confidence in interpretations of sinuous ridges as eskers, as exemplified by decades of debate surrounding the origin of the candidate eskers in the Dorsa Argentea Formation (Butcher et al., 2016;Head & Pratt, 2001;Kargel, 1993;Parker et al., 1986;Ruff & Greeley, 1990;Tanaka et al., 2014;Tanaka & Kolb, 2001;Tanaka & Scott, 1987). Hence, while there is some evidence for a handful of other candidate eskers on Mars at a range of latitudes, thus far, the only two documented examples of eskers emerging from existing parent glaciers are seen at mid-latitude and are of late Amazonian age.
Amazonian glaciation of the low latitudes is widely thought to have occurred during periods of high (>45°) planetary obliquity that were climatically distinct from the periods of intermediate (~35°) obliquity during which mid-latitude glaciers advanced (Fastook et al., 2008;Madeleine et al., 2009). Thus, the paleoenvironmental implications of possible glacial melt-production on Arsia Mons are different to those of the mid-latitude VFF-linked eskers discussed in the present study. Nonetheless, it is interesting that the rift-occupying tectonic settings of the candidate eskers in Tempe Terra and Phlegra Montes (Gallagher & Balme, 2015) are similar and that Scanlon et al. (2014Scanlon et al. ( , 2015 observe a possible esker of similar age on the flanks of a volcano. Gallagher & Balme (2015) suggest an enhanced geothermal heat flow in the rift setting of the candidate esker in Phlegra Montes, and Scanlon et al. (2014Scanlon et al. ( , 2015 invoke a subglacial magmatic intrusion to explain the presence of an esker in deposits that were largely deposited by cold-based glaciers under cold Amazonian climates (e.g., Kadish et al., 2014;Shean et al., 2005). Hence, we now consider the influence of geologic setting on basal meltwater production by mid-latitude VFFs.

Geologic Setting and Environmental Controls on Melting
Terrestrial graben/rifts are commonly associated with elevated geothermal heat flux, for example, due to heating from intrusive magma bodies (e.g., dikes), volcanic plutons, and/or hydrothermal circulation (e.g., Fernàndez & Banda, 1990). Enhanced geothermal heat flux within terrestrial graben/rifts can influence the production of basal meltwater by glaciers, which occupy them (e.g., Jordan et al., 2010;Schroeder et al., 2014). Identification of mid-latitude VFF-linked candidate eskers within graben/rifts in Phlegra Montes and now Tempe Terra leads us to consider whether a genetic link exists between the setting and the landform. Thus far, a lack of evidence for basal melting of existing VFFs elsewhere in Mars' mid-latitudes implies that broad-scale, cyclic (~120 kyr), obliquity-induced mid-latitude climate warming (e.g., Madeleine et al., 2009) alone was insufficient to melt basal ice during the Amazonian. Hence, mid-latitude Amazonian glaciation was broadly cold-based. A regional warm climate anomaly would have to persist over long timescales in order to induce and maintain temperate conditions at the glacier bed. The lack of morphological evidence for prolonged supraglacial or subglacial melting (e.g., supraglacial channels, fans, and streamlined bedforms) leads us to exclude such a regional climate anomaly as the sole driver of wet-based glaciation in the Tempe Journal of Geophysical Research: Planets 10.1002/2017JE005434 Terra rift. Thus, while we do not exclude climate warming as a possible contributing factor in inducing wetbased glaciation, we argue that perturbation of basal ice temperatures by a geothermal heat flux anomaly (Gallagher & Balme, 2015) is more consistent with the tectonic setting and our geomorphic observations. While primary geothermal anomalies associated with the formation of the ancient rift in Tempe Terra (Hauber et al., 2010) are unlikely to have persisted into the Amazonian, magmato-tectonic processes could have exploited preexisting structural weaknesses in their vicinity to induce geothermal anomalies coincident with recent glaciation. Given uncertainties over the precise ages of features on planetary surfaces, and the current lack of available ground-truth for global geothermal heat flux models, reliable identification of a specific synglacial source of geothermal heat to the rift cannot be achieved with confidence. Therefore, we do not aim to identify a specific geothermal source event. Instead, we highlight evidence that magmato-tectonic activity in Tempe Terra is unlikely to have ceased at the termination of the last major phase of faulting. Thus, we argue that a recent phase of above-average geothermal heat flux in the vicinity of the rift cannot be excluded as a possible heat source for wet-based glaciation.
Multiple generations of Hesperian-to Amazonian-aged graben dissect the region surrounding the Tempe Terra rift (e.g., Scott & Dohm, 1990a, 1990bTanaka, 1990). Giant magmatic dike swarms emanating from the Tharsis volcanic province have been posited as a driver for their initial formation (Mège & Masson, 1996). We identify evidence for recent, localized reactivation of portions of these faults proximal to the rift in Tempe Terra.
A~3,000 km long graben system originating in Tractus Catena (SW) and terminating in the northern plains (NE) passes directly through the center of the rift (Figure 7a). Portions of this system have been associated with middle to late Amazonian (~0.25-0.7 Gyr ago) reactivation of early-Hesperian to mid-Amazonian-aged (~0.7-3.8 Ga) faults by magmatism at Ascraeus Mons (Tanaka, 1990). Small, fresh, crosscutting faults in the rift-marginal fault blocks, and possibly even in the latitude-dependent mantle, support recent localized reactivation of portions of the fault system proximal to and within the rift (Figures 7b and 7c). Additional (up to~40 km-long) fault scarps (e.g., Figure 7d)~10 km south east of the rift, which are subparallel to those within the rift-marginal fault blocks (i.e., possibly related to the same stress field) are comparable in morphological "freshness" to some of the youngest graben scarps on Mars, such as Cerberus Fossae (Figure 7e) (e.g., Taylor et al., 2013). Using standard techniques in NASA Ames Stereo Pipeline (Shean et al., 2016), we generated a 5 m/pixel digital elevation model of the fault scarp in Figure 7d from HiRISE images ESP_052553_2260 and ESP_052065_2260, vertically aligned to High Resolution Stereo Camera image H1528_0000. Average slopes on the fault scarp range from 30°to 55°. Such steep slopes, above the angle of repose, could be attributed to recent formation of the fault scarp, and the resulting short timescales available for scarp degradation and lowering of slopes. Furthermore, the young latitude-dependent mantle, which blankets adjacent surfaces and faults, is absent on the fresh fault scarps, possibly implying faulting since emplacement of mantle materials. A continuous "sheet" of material appears to have been mobilized downslope along an~1.5 km long portion of the fault scarp in Figure 7d (inset, black arrows). We consider it unlikely that local slope instabilities within individual alcoves along the scarp could mobilize a continuous lobe with such a large lateral extent. Therefore, we tentatively propose that seismicity associated with recent fault movements could be a more plausible mechanism for initiating such slope failures.
Such spatially restricted fault reactivation is more likely to have been induced by local tectonism than distal activity in the source regions of the graben (i.e., the Tharsis volcanoes). Impact crater size-frequency distributions for lava flows within the "Tempe Volcanic Province" indicate very recent (100-800 Ma) volcanism 500 km SSW of the rift Plescia, 1981). While the scale of, and distance to, the Tempe Volcanic Province likely excludes it as a direct cause of rift-proximal fault reactivation, or a source of geothermal heat to the rift, this highlights that very recent volcanism, likely derived from shallow melt sources (Plescia, 1981), has occurred within the Tempe Terra region. Noneruptive, intrusive magmatism, which could manifest at the surface as reactivated faults, is likely to have persisted into more recent times than eruptive volcanism. However, identifying a specific mechanism for the origin of the fresh fault scarps would involve significant unknowns, so we do not speculate further.
If our interpretations of the LDA-distal knobbly ridges as terminal moraines and/or our interpretation of the polygonized wedge as older glacial sediments are correct, faults that cut through them (Figure 7f (Table S1) examples of morphologically fresh (i.e., recently active?) fault scarps (white arrows) in and around the Tempe Terra rift, and comparison to young fault scarps in Cerberus Fossae: (a) MOLA elevation on THEMIS daytime infrared context map (Table S1) (Table S1)  wedge appear to exert a structural control upon the morphology of the scarp that passes though it and the overlying knobbly ridge. Thus, it is also plausible that the fault predates the most recent glacial advance within the rift and that its apparently morphologically "fresh" expression within these units can be attributed to modification of scarp degradation processes by the polygonized wedge, and associated disruption of the overlying knobbly ridge. These faults also manifest, with variable scarp morphology, within the LDA-proximal floor unit, which we interpret as a deglaciated foreland. In some locations (predominantly in the northern portion of the unit where surficial sedimentary deposits appear thinner), they have prominent scarps, while in others (predominantly in the southern portion of the unit where LDA-proximal floor materials appear slightly thicker), they are muted. Such muted expressions could result from faulting within unconsolidated (e.g., glacigenic) sediment (synglacial faulting), or deposition over a preexisting scarp (preglacial faulting). Thus, it is difficult to determine whether they predate or postdate emplacement of this unit. No faults manifest within LDA surfaces, but this does not refute small-scale synglacial faulting, particularly given the small vertical expression of the rift floor faults, since viscous relaxation, ice flow, and mantling processes would be expected to erase small fault traces on short timescales.
All three examples of relatively fresh fault scarp populations discussed here were unmapped in Viking-era studies, which concluded that formation of the major graben in this region was largely, yet not entirely, complete by the middle to late Amazonian (Scott & Dohm, 1990a;Tanaka, 1990). It is therefore conceivable that small-scale tectonism might have continued into the late Amazonian. Thus, the existence of positive geothermal anomalies in the vicinity of the rift (for example due to magmatic dike propagation into the near subsurface), in the period during which VFF formed, is plausible. Figure 8a shows the calculated basal ice temperatures assuming that the only source of basal heat is geothermal (i.e., using equation (1) only). Basal temperatures close to 273 K can be produced for some combinations of the range of parameters tested but only for warm surface temperatures (>215 K), thick ice (>1100 m), and high geothermal heat fluxes (>0.08 W m À2 ). These conditions may be unrealistic, although such locally high geothermal heat flux relative to modeled Martian global average (0.023-0.027 W m À2 ; Plesa et al., 2016) may be plausible in the vicinity of a near-surface magmatic dike (e.g., Cassanelli et al., 2015;Wilson & Head, 2002).

Quantitative Constraints on Requisite Geothermal Heat Flux
Our calculations show, however, that internal strain heating can also provide a significant heat source, especially at high assumed driving (shear) stresses, and for warmer mean annual surface temperatures. This occurs due to the highly nonlinear behavior described by equation (2). Shear stress (dependent on the ice thickness and surface slope) is raised to the fourth power in equation (2), and the value of A T increases by 3 orders of magnitude between 220 K and 270 K (e.g., Cuffey & Paterson, 2010). For a driving stress of 100 kPa, for instance, basal temperatures close to 273 K can be produced for surface temperatures of 205 K, or beneath 900 m thick ice, or for geothermal heat fluxes of 0.05 W m À2 , effectively halving the required geothermal anomaly above the modeled global average compared with geothermal heating alone. Internal strain heating can potentially generate over 0.05 W m À2 of additional heating, warming the basal ice by up to 14.5 K for some parameter combinations (Figures 8c and 8d), and making the presence of warm basal ice more likely under a wider range of glaciological, geothermal, and climatic conditions. Considering uncertainty over the effect of impurities upon the rheology of Martian VFF compared to clean terrestrial ice (for which the Cuffey & Paterson, 2010, value of A T is applicable), a factor of 2 increase in A T (i.e., softer ice) increases the maximum warming in the parameter space we investigate to 18.5 K compared to 14.5 K for clean ice (Figure 8d). Halving A T (i.e., harder ice) reduces the maximum warming to 11.1 K.
It is important to note that while our calculations are based on an assumed driving stress, it is not necessary to assume that high driving stresses (which favor increased strain heating) must occur ubiquitously. Karlsson et al. (2015) find that, on average, Martian VFF may comprise softer ice than typical terrestrial glaciers, with average driving stresses of~20-35 kPa. However, even for low-average driving stresses, flow convergence will produce localized areas of ice with high driving stresses due to locally thicker ice and/or steeper surface slopes, and therefore with higher strain heating. Locally, thicker ice in itself also makes warmer basal ice more likely. Convergence could occur where ice flows towards a valley bottom from opposite valley sides or where ice flows out from steeper topography to shallower slopes.
The candidate esker we identify here seems to occur in an area of potential topographic flow convergence. The ridge occurs in an area of strong horizontal flow convergence from the western and eastern flanks of the rift (Figure 3a), and where the along-profile slope decreases markedly (Figure 4). Such a situation would seem very likely to produce a zone of higher shear stresses and locally thicker ice, increasing the strain heating markedly, and facilitating basal melt especially if geothermal heating was also elevated. The fact that the parent LDA would thin toward its margin (leading to colder basal temperatures) does not preclude esker formation, as trapped meltwater has been observed to penetrate the cold margins of terrestrial polythermal glaciers in channelized systems (Bingham et al., 2005(Bingham et al., , 2008Kavanaugh & Clarke, 2001).
The results presented here provide a first-order estimation of the ice thickness, surface temperature, geothermal heat, and driving stress conditions required for basal melting. Higher-order, three-dimensional, thermomechanically-coupled modeling of VFFs under a range of climatic and geothermal forcings, in order to evaluate the likelihood of basal melting and the possible contribution of internal strain heating, is ongoing.

Conclusions
We document a geomorphic assemblage associated with a putative late Amazonian-aged debris-covered glacier within a tectonic rift in Tempe Terra, in Mars' mid-latitudes. We interpret a sinuous ridge emerging from the glacier as an esker formed by deposition of sediment in a pressurized subglacial meltwater Journal of Geophysical Research: Planets 10.1002/2017JE005434 channel during a transient phase of wet-based glaciation. To our knowledge, candidate eskers in Tempe Terra (this study) and Phlegra Montes (Gallagher & Balme, 2015) are the only documented examples that emerge from existing parent VFFs in Mars' mid-latitudes. Their parent VFFs both occupy graben-like rifts within which late Amazonian Tharsis/Tempe Terra and Elysium magmatism could plausibly have enhanced geothermal heat flux. Small, fresh faults in Tempe Terra could be a surface manifestation of such late-stage magmatism. A lack of documented evidence (thus far) for widespread basal melting of mid-latitude VFFs is widely interpreted to indicate that atmospheric warming alone was insufficient for wet-based glaciation. Therefore, enhanced geothermal heat flux, possibly supplemented by internal strain heating, may have been a prerequisite for basal melting. We predict that eskers are most likely to be found in association with VFFs within faulted zones close to recently active volcanic regions. Detailed global-scale surveys of mid-latitude VFFs, focusing on identifying the presence or absence of eskers and/or other morphologies indicative of melting, could help to test this hypothesis, together with more complex three-dimensional, thermo-mechanicallycoupled modeling of the environmental conditions required for basal melting of VFFs. for their thoughtful and constructive reviews and Jim Head for insightful comments on a previous version of the manuscript. We thank Joel Davis for insightful discussions of the inverted channel hypothesis. Funding was provided by STFC grants ST/N50421X/1 (F.E.G.B.) and ST/L000776/1 (M.R.B./A.H./ S.R.L.). S.J.C. is supported by the French space agency CNES. We are grateful to the Mars instrument teams and particularly to the HiRISE team for their fulfillment of our imaging requests. A list of data products used in this study and the full extent of the geomorphic map, crater count plots, and animated gifs of the data in Figures 8a-8d are available as supporting information to this article (Table S1, Figures S1 and S2, and Movies S1-S4, respectively). MOLA and THEMIS data products can be downloaded from the U.S. Geological Survey Planetary GIS Web Server (https://webgis.wr.usgs.gov/). CTX images can be downloaded from the Arizona State University Image Explorer (http://viewer.mars.asu.edu/viewer/ctx). HiRISE images can be downloaded from https://hirise.lpl.arizona.edu/. HRSC data products are can be downloaded from Freie Universität Berlin HRSCview data explorer (http://hrscview.fu-berlin.de/ cgi-bin/ion-p?page=entry2.ion). Map shapefiles, raw impact crater count data, and HiRISE digital elevation models available from corresponding author upon request.