Meltwater generation in ice stream shear margins: case study in Antarctic ice streams

Liquid water within glacier ice and at the glacier beds exerts a significant control on ice flow and glacier stability through a number of processes, including altering the rheology of the ice and lubricating the bed. Some of this water is generated as melt from regions of rapid deformation, including shear margins, due to heating by viscous dissipation. However, how much meltwater is generated and drained from shear margins remains unclear. Here, we apply a model that describes the evolution of ice temperature, melting, and water transport within deforming ice to estimate the flux of meltwater from shear margins in glaciers. We estimate the flux of meltwater from temperate ice zones in three Antarctic regions: Bindschadler and MacAyeal Ice Streams, Pine Island Glacier, and Byrd Glacier. We show that the flux of meltwater from shear margins in these regions may be as significant as the meltwater produced by frictional heating at the bed, with average fluxes of ∼0.005–0.1 m yr−1. This contribution of shear heating to meltwater flux at the bed may thus affect both the rheology of the ice as well as sliding at the bed, both key controls on fast ice flow.


Introduction
Fast-flowing glaciers and ice streams drain a significant fraction of the Greenland and Antarctic ice sheets [1][2][3][4]. The speed of flow in many ice streams is largely controlled by friction at the bed (the ice-sediment interface) and lateral shear stresses in the margins [5][6][7]. Many of these ice streams flow over complex hydrologic systems that transport water at the glacier bed to the ocean [8][9][10][11][12], and the presence of liquid water in glaciers and in these channels affects both basal friction and lateral shear stresses in the margins. The amount and distribution of water at the bed alters the basal friction of glaciers, as water acts as a lubricant and changes the characteristics of the sediments the ice slides over, thereby affecting the speed of flow [13][14][15][16]. In some regions, water collects into subglacial lakes which may store significant amounts of water. Observations have shown subglacial floods to cause acceleration of glaciers, ultimately affecting the amount of ice discharge [17][18][19][20]. Whether meltwater is held within the ice or drained to the bed also may affect the rheology of ice in the margins by softening the ice [21][22][23][24] and may partially set the width of ice streams that are not topographically controlled [25]. Additionally, meltwater generation in ice shelves may enable flux of freshwater to the ocean, affecting ocean circulation and the ice-ocean interface which may control ice shelf stability [26]. Therefore, meltwater has a significant effect on the flow speed and stability of rapidly deforming glaciers.
NASA's Ice, Cloud, and Land Elevation Satellite (ICESat and ICESat-2) enabled observation of glacier hydrology, including distribution of subglacial lakes and hydrologic networks [9] and estimates of and changes to water volume [10,12,[27][28][29]. While observations have improved our understanding of how the subglacial hydrologic system is changing and affecting glacier flow, a characterization of the physical processes generating the water in these systems remains incomplete and is necessary for projections of subglacial hydrology and its future effects on fast flow given a changing climate.
In Antarctica, much of the water in the subglacial hydrologic system is generated at the glacier bed through melting by geothermal heat and basal friction, while in Greenland, most of the water is generated by surface melt due to warm temperatures and percolates through moulins. However, some meltwater is generated within zones of temperate ice in the ice column [30]. This temperate ice is usually created through shear heating (work done during rapid ice deformation being dissipated as heat) and is generally found in regions of significant shear, such as the base of ice sheets where the ice is frozen to the bed and in the margins of ice streams [31][32][33][34][35]. Previous work has used models of ice temperature to suggest that extensive temperate zones in the margins of Antarctic ice streams may exist due to the rapid deformation in these regions [36,37]. Since temperate ice is porous, meltwater generated in temperate zones can percolate through the temperate zone and eventually drain into the hydrologic channels at the bed [30,34,38]. The contribution of temperate ice zones generated by shear heating to meltwater production, however, is currently poorly characterized.
Schoof & Hewitt [39] derive a model for polythermal glaciers that computes the flux of melt from temperate ice in an ice column, and Hewitt & Schoof [40] present modifications to this analytical model that can be applied to ice sheet models. We apply this coupled model to shear margins in the Antarctic Ice Sheet using remotely sensed observations to estimate the meltwater flux out of shear margins of the Antarctic Ice Sheet. We further compare the contributions of shear margins to basal hydrologic channels as compared to contributions from geothermal heating and frictional sliding at the bed to suggest that shear margins are a significant contributor to melt at the bed. Finally, we use these models to evaluate the effect of vertical velocity structure on estimates of melt at the bed.

Model for effective pressure, porosity and meltwater flux
We present a model for temperature and meltwater flux in shear margins of ice streams, initially derived by Schoof & Hewitt [39]. The following model was derived and presented in [39] and we describe it here in order to outline the physics behind and assumptions underlying the equations applied in this study.
Heat is introduced into ice primarily through frictional heat at the bed, geothermal heat at the bed, and ice deformation internally. In shear margins of ice streams, internal ice deformation is dominated by lateral shear [41][42][43]. Heat generated in the margins due to lateral shear is then conducted and advected both laterally and vertically. In areas where there is enough heat to bring the ice up to its melting point, a temperate zone can be generated. Temperate ice is defined as a mixture of ice and liquid meltwater that is at the pressure melting temperature. We use a fixed melting temperature of T m = 273 K, the pressure melting point of ice at approximately atmospheric pressure conditions. This model neglects the effects of varying pressure on the melting temperature, a reasonable assumption given that the melting temperature of ice varies by less than 2 K at pressures expected in ice sheets, which is within the uncertainties of the model. Temperate ice is porous, with dimensional porosityφ that is affected by a melting rate m. Ice in temperate zones is approximated as incompressible, although porosity will realistically make temperate ice slightly compressible [40]. The heat and moisture equations are then defined, respectively, as: and ∂φ where tildes imply dimensionality in terms that will later be non-dimensionalized, and with the constraintsT < T m ,φ > 0, and φ(T m −T) = 0. T is ice temperature, u is ice velocity, ρ I is the mass density for ice, c p is the specific heat capacity of ice, K is thermal conductivity, ρ w is the mass density of water, L is the specific latent heat of fusion, J is water flux, and W = σ ijεij is the rate of heating through viscous dissipation (hereafter called shear heating, where σ ij is the deviatoric stress tensor andε ij is the strain-rate tensor). Applying the constitutive relation for ice, W can be written as W = A −1/nε(n+1)/n , where˙ = (1/2)˙ ij˙ ij is the effective strain rate and A is the ice softness parameter which is inversely proportional to ice viscosity η I such that 2η I = A −1/nε(1−n)/n [44]. The parameter n is the viscous stress exponent. In this study, we set the stress exponent n = 3, a value supported by field data [45], and we set the ice softness parameter A to be the value associated with temperate ice [46]. Schoof & Hewitt [39] use these governing equations to estimate the flux of water out of temperate ice zones (J). To illustrate this, we redefine equations (2.1a,b) in terms of enthalpyH, a term that defines the internal energy of the system along with its pressure and volume and which is particularly useful for multiphase fluid flow. We define enthalpy as where T m is the constant melting temperature. This model assumes mass conservation ∂ρ ∂t + ∇ · (ρu) = 0, (2.3) in whichρ = (1 −φ)ρ I and ρ I is assumed to be constant. Applying this, equations (2.1a,b) can be combined as ∂H Schoof & Hewitt [39] define meltwater fluxJ using Darcy's Law, a model for fluid flow through a porous medium due to gravitational body forces and the gradient of effective pressure, as where k 0 is the coefficient of hydraulic permeability, α is the porosity exponent that sets the dependence of hydraulic diffusion on porosity, η w is the water viscosity, g is a scalar for gravity, z is the unit vector in the z-direction (such that gravity is positive vertically), and N is effective pressure, defined as the difference between the pressure in the ice and pore water pressure. This pressure difference can be related to the rate of ice compaction bỹ (2.6) in which η I /φ acts as a bulk ice viscosity and η I is ice viscosity [38,39]. Assuming incompressibility and mass conservation, this can be rewritten asÑ = η Ĩ φ ∇ ·J. Thus, equations (2.4) and (2.5) can be re-expressed as coupled equations: with the boundary conditionJ · n = 0 requiring that there is no water flux at the boundary between cold and temperate ice, where n is the normal vector relative to the cold-temperate boundary. Equations (2.7a,b) represent a general, three-dimensional model for the evolution of temperature, effective pressure and porosity in temperate ice. In this study, we seek to apply an analytical model in order to improve our physical understanding of the processes underlying the generation of melt in shear margins. Applying this simple framework to natural systems, we can explore the dependencies on key parameters as well as the physics that underlies those dependencies.
Therefore, we apply a version of the model derived in [39] with simplifying assumptions. First, we neglect horizontal advective terms and assume that vertical advection is the dominant advective term. Neglecting lateral advection may affect the estimates produced by increasing the extent of temperate ice zones and affecting vertical velocities. The effects of this simplification are discussed further in §5. We further assume that horizontal thermal diffusion is negligible, thermal conductivity K is constant with depth, and there is no melting at the ice surface. Therefore, equations (2.7a,b) can be written in a one-dimensional formulation as where a is rate of ice accumulation. The boundary conditions at the cold-temperate boundary are Further, both porosity and effective pressure are defined as being 0 at and above the coldtemperate boundary. The boundary condition at the surface is T = T s , where T s is the surface temperature. The boundary conditions at the bed are N = N 0 . For simplicity moving forward, Schoof & Hewitt [39] non-dimensionalize equations (2.8a,b). To do so, they define scales for enthalpyH, heightz, ice temperatureT, porosityφ, effective pressureÑ, timet, and meltwater fluxJ based on ice thickness h, specific heat capacity of ice c p , mass density of ice ρ I , a scale for ice viscosity [η I ], thermal conductivity K, the difference between melting temperature and surface temperature T = T m − T s , and latent heat of fusion L. These non-dimensionalizations, defined in [39] arẽ such that H, z, T, φ, N, t, J are dimensionless. The full non-dimensionalization is presented in the electronic supplementary material, and the resulting non-dimensionalized versions of equations defined as a function of the non-dimensional numbers: , the contribution of effective pressure to moisture flux (2.12c) Kη w , the ratio of heat advected in the meltwater to thermal conduction in the ice (2.12d) where w is the vertical velocity, which is defined on a coordinate system such that z is positive upwards. In this work, we apply a steady-state solution found in [39], and in future work will expand to consider time-dependent solutions. The steady-state forms of equations (2.11a,b) are with the respective as described above. Equation (2.13b) is solved analytically with the important assumption that ice viscosity is constant with depth. This assumption is discussed further in §5. The parameter δ in equations (2.13a-c) is small, and therefore at first order, we would like to neglect the term δ∂N/∂z, but this is a singular perturbation [47]. For this reason, solving the  equations directly (what is referred to as the outer solution) provides a reasonable but not entirely accurate solution (figure 1) for both porosity and meltwater flux. To produce fully accurate solutions, the higher-order terms must be evaluated. Schoof & Hewitt [39] derive asymptotic solutions to equations (2.13a,c) in both the outer layer (above the boundary layer) and the inner layer (within the boundary layer; this solution is defined as the inner solution). This model depends on the choices for the vertical velocity w structure. Here, we restate solutions for the case of a constant vertical velocity structure (w = −a) and a linear vertical velocity structure (w = −az) to explore the implications of the structure of w to the resulting meltwater flux at the bed.

(a) Constant vertical velocity
We first assume that the vertical velocity is constant through the ice column and equal to the surface mass balance a. In assuming a constant vertical velocity equivalent to the surface mass balance, the Darcy flux J we estimate here describes both the flux of melt out of temperate ice zones as well as melt that passes the basal boundary from the vertical advection of ice in the ice column. Physically, this would advect more cold ice through the ice column, counteracting shear heating and possibly reducing the amount of melt generated in shear margins. While this assumption is not necessarily physical, it is an end-member case and accounts to some extent for the lack of lateral advection in the model.
Since the boundary layer thickness is defined through δ, to find a solution far from the boundary layer, [39] neglects the terms with the highest order of δ. For the outer solution of porosity, we substitute equation (2.13a) into equation (2.13c), with the condition that porosity must be zero at the boundary between cold ice and temperate ice (φ(z ct ) = 0, where z ct is the depth of the top of the temperate zone): (2.14) Integrating over depth and applying the boundary condition on porosity, they find  [39] rescale dimensionless height z by a factor of δ β . This rescaling allows a focus on the boundary layer, in whichẑ ∼ O(1). The resulting inner solutions are higher-order corrections to outer solutions (2.16) and (2.15). This correction decays exponentially away from the boundary layer at a rate of √ z(z/H)δ −β . We then combine outer and inner solutions for effective pressure and ice porosity into composite solutions valid uniformly throughout the domain. Composite solutions are constructed by adding the outer and inner solutions and subtracting the overlap [47,48]. The full derivation is presented in [39], resulting in the following composite solutions: These are solutions for N and φ that are valid for the full temperate zone. These solutions can be substituted into a one-dimensional form of equation (2.5) to find the resulting meltwater flux out of the temperate zone:

) Linear vertical velocity
We next present a case with a linear vertical velocity profile of w = −az, where z is the nondimensional height. This can be considered more physical, as the vertical velocity goes to w = 0 at the bed. Therefore, mass is conserved and J(z = 0) can be clearly interpreted as being the flux due to water draining from temperate ice zones. We can also expect the estimated fluxes to be between the endmember cases presented in figure 5. Below, we derive the model for a linear vertical velocity case, similar to a case derived in [39], and we apply it to show the effect of this vertical velocity profile. The application of a linear vertical velocity profile here affects equations (2.13a,b), which are now written as and A higher-order correction (the inner solution) is needed in order to produce complete estimates. However, due to the added complexity of a linear vertical velocity profile, the inner solutions cannot be written analytically. A solution has been described in [39]. Since we showed in the previous section that the outer solution is a good approximation of the composite solution, for the case of a linear vertical velocity we only present the outer solutions. However, for future work we propose to present in more detail the effects of vertical velocity profiles incorporating the higher-order terms. The outer solutions for ice porosity and effective pressure with a linear vertical velocity profile have been discussed in [39] and are (2.21) and Equation (2.21) can be solved numerically to find the outer solution for porosity and then the solution can be applied to equation (2.22) to find the outer solution for effective pressure.

Effect of vertical velocity structure in idealized setting
To compare the accuracy of the outer and composite solutions, we compare these solutions with results from a numerical model ( figure 1). The numerical model solves equations (2.8a,b) in a finite volume implementation with a mesh of dz = h/256 and where steady state is defined such that the sum of squares of differences between iterations is less than 10 −8 . The full numerical model is described in [38]. In this idealized simulation, we let h = 200 m, surface slope be 4 • , surface temperature be 272.15 K, and surface accumulation be 0.2 m yr −1 . Further, we set α = 2.33. We then calculate the relevant non-dimensional parameters by equation (2.12). This rate of shear heating and ice thickness results in a temperate zone forming at approximately 68% of the ice thickness (figure 1a). For effective pressure and porosity, the analytical expressions follow the numerics closely. Effective pressure is large at the cold-temperate boundary (z ct ) and decreases approximately linearly as it gets closer to the boundary layer. In the boundary layer, effective pressure begins to decrease exponentially (figure 1b). Porosity, on the other hand, is zero at the cold-temperate boundary and increases approximately linearly until the boundary layer, at which point it increases with a linear and an exponential term (figure 1c). The composite solutions for pressure and porosity follow the outer solution above the boundary layer and the inner solution within the boundary layer. The flux of meltwater out of the temperate zone increases down the ice column, with q composite (0) = −7.28 (figure 1d). Since the flux at the bed computed by numerics is q numerics = −7.32, our analytical estimate carries a small (less than 1%) error. Note that, as described above, the inner solution is a higher-order correction, and the outer solution provides a good approximation to the numerical estimate of flux (q outer (0) = −7.50). Comparison to numerical estimates for another rate of shear heating is presented in the electronic supplementary material to demonstrate the generality of the agreement between analytics and numerics.
In figure 2, we compare the solutions for a linear vertical velocity case and a constant vertical velocity case. Here, ice temperature is computed analytically for the constant vertical velocity case and numerically for the linear vertical velocity case. The linear vertical velocity changes the  temperature profile only very slightly, producing a temperate zone thickness of 0.6191, where a constant vertical velocity produces a nearly identical temperate zone thickness of 0.6186. As expected, a linear vertical velocity profile produces a larger (though negligibly so) temperate zone, due to less cold ice being advected down the ice column. While the ice temperature profiles are largely the same, the porosity profiles (and therefore pressure and flux profiles) differ significantly. The linear vertical velocity case produces larger porosities at the bed, likely due to the reduction in the advection of less porous (or non-porous) ice down the ice column. The difference in fluxes is significant, with the constant vertical velocity producing a non-dimensional flux of approximately −4 and the linear vertical velocity producing a non-dimensional flux of approximately −8. In this case, the difference is likely due to the fact that, in the linear vertical velocity, w → 0 at the bed, meaning that the Darcy flux is primarily responsible for advecting enthalpy through the bottom of the ice column (whereas in the case of a constant vertical velocity, the vertical velocity partially accounted for this). This will result in the linear vertical velocity case producing a larger flux of meltwater at the bed. This explanation holds in the cases where ice temperature is not affected by the vertical velocity profile; we will see later in this paper that this is not generally true.

Estimates of basal flux in Antarctic ice streams
We apply the model for effective pressure, porosity and meltwater flux to find estimates of meltwater production in temperate zones of shear margins in Antarctic glaciers. We estimate meltwater flux J| z=0 at the bed along three Antarctic shear margins: Bindschadler Ice Stream (figure 5a), Byrd Glacier (figure 5b) and Pine Island Glacier (figure 5c). The model takes as an input strain rate, which here is computed from Landsat-8 velocity fields [49] and has a spatial resolution of 240 × 240 m, ice thickness, which here is found from BedMachine v01 [50], and surface mass balance, which here is found from RACMO estimates [51]. We assume that κ = 0.52, α = 2.33, and δ = 0.001. For the basal boundary condition for effective pressure, N 0 , we assume that N 0 = 1, as   3c, figure 3d).
In the linear vertical velocity case, the temperate ice zones extend over a wider area, producing larger fluxes of greater than 0.02 m yr −1 . This suggests, unlike the example in figure 2, the vertical velocity profile here affects the temperature profile. Thus, the difference in meltwater flux between the constant vertical velocity case and the linear vertical velocity case is a combination of the Darcy flux accounting for the advection of all enthalpy in the column in the case of w → 0, as discussed in the previous section, along with the change in ice temperature due to the advection of cold ice through the ice column counteracting shear heating in the case of w == a.
To further examine the structures of these temperate ice zones, we compute full depth profiles of effective pressure N, porosity φ and meltwater flux along three Antarctic shear margins: Pine Island Glacier (figure 4a), Bindschadler Ice Stream (figure 4b) and Byrd Glacier (figure 4c). The shapes of the depth profiles are strongly dependent on the ill-constrained parameter κ, and in the electronic supplementary material we show estimates of depth profiles in Pine Island Glacier for a wide range of κ. While varying κ does not have a strong influence on estimates of flux at the bed, it does affect estimates through the column, particularly of effective pressure.
All three shear margins contain temperate zones. In Bindschadler Ice Stream, this temperate zone is concentrated near the grounding line due to an increase in strain rate driving shear heating [37,49]. In Byrd Glacier, there is a region of zero temperate zone thickness (for a constant vertical velocity) or smaller temperate zone thickness (for a linear vertical velocity) in the middle of the transect. Pine Island Glacier has a temperate zone for most of the transect, out to approximately 80 km upstream of the grounding line, due to its fast flow.
In all three shear margins, effective pressure is higher in the ice column (approx. 1 − 4 kPa) and decreases sharply near the bed, due to the prescribed boundary condition of N 0 = 1. Porosity is high in both Pine Island Glacier and Byrd Glacier where there are large strain rates and thus    rates of shear heating may be largest in Pine Island Glacier, due to larger strain rates, the flux of meltwater is also dependent upon ice thickness. Since Byrd Glacier has, on average, thinner ice than Pine Island Glacier, this translates into a similar magnitude of meltwater flux as in Pine Island Glacier. There is a region of increased porosity near the bed in Bindschadler Ice Stream approximately 5-10 km from the grounding line, which correlates with a region of acceleration seen in Landsat-8 velocity fields. Previous work has suggested that this region of acceleration may be due to changes in effective pressure at the bed due to subglacial hydrologic channels [38], and therefore the estimates provided here may provide a link to estimating the effect of meltwater flux from shear margins on glacial acceleration. On Pine Island Glacier, there is a reduced meltwater flux near the grounding line due to a region of reduced porosity, which when combined with the large effective pressure at the bed translates to a near-zero meltwater flux.
As we have illustrated in this study, the magnitudes of basal meltwater flux are sensitive to the choice of the vertical velocity structure. Here, we have evaluated two structures: a constant vertical velocity w = −a and a linear vertical velocity w = −az. These represent only two potential cases for the vertical velocity profile, and thus it is valuable to look at what the full range of possible meltwater fluxes may be for varying vertical velocity structures. The end member cases are w = 0, in which there is no advection of cold ice through the ice column, and w = −a, in which all of the ice accumulated at the surface is advected through the ice column. In figure 5, we examine the range of potential fluxes in three Antarctic shear margins between these two endmember cases. Because our model becomes undefined at w = 0, we let the w ≈ 0 be defined as w being some very small fraction of a.
The case of w ≈ 0 produces larger fluxes than w = −a, likely due to a combination of two factors. First, w ≈ 0 decreases the advection of cold ice down the ice column, allowing for faster rates of heating. Therefore, we would expect larger temperate zones in this case (illustrated further in the electronic supplementary material, Section D). Second, w ≈ 0 requires the Darcy flux to account for the advection of enthalpy through the bottom of the ice column, resulting in a larger flux. In the electronic supplementary material, Section D, we show that most of the variation in fluxes comes from the reduction in advection of nonporous ice, rather than a change in the temperate zone thickness.
In many cases, the range is large, with Pine Island Glacier producing a range of fluxes from 0.05 to -0.3 m yr −1 and Bindschadler Ice Stream producing a range from 0.01 to 0.04 m yr −1 . In Byrd Glacier, there is a less than or equal to 0.1 m yr −1 difference between the two estimates. Without clear constraints on the profile of vertical velocity, we take these ranges as the estimates of uncertainties of meltwater flux given uncertainties produced from our parametrization of vertical velocity.

Discussion (a) Contributions of shear margins to meltwater at the bed
These results quantify the contribution of shear heating to meltwater in basal channels and cavities. Other potential sources of meltwater at the bed include melting by geothermal heat flux and melting by friction at the ice-bed interface that is generated due to ice sliding over the basal surface [54,55]. In some regions, surface melt may also be an important source, though we neglect it given cold surface temperatures in Antarctica. The contributions from geothermal heat flux and frictional melting at the bed were modelled by Robel et al. [56], incorporating melting by geothermal heat flux, sliding at the bed and conduction into the ice. We apply the same model to estimate the comparative contributions of geothermal heat and sliding at the bed to our new estimates of melting by shear heating in margins. We neglect the effect of thermal conduction into the ice as there is no conduction into a temperate ice zone due to the constant temperature. Following [56], we then define the meltwater supply at the bed s as (2.18). To find melting by sliding at the bed, we assume that vertical shearing is negligible such that u s ≈ u b , where u s is surface velocity. This is a reasonable assumption in most West Antarctic ice streams, where the ice slips over weak till and sediments [6,7,57]. This assumption provides an upper bound on melt rate from frictional heating. To find basal stress τ b , we assume that the till that the ice slides over can be approximated as a perfectly plastic material, an assumption supported by laboratory experiments [13,54,58,59] and inferences from observations [60,61], such that the yield stress τ * = c 0 + μN, where c 0 is the cohesion of the material and μ is the internal friction coefficient. Laboratory studies have found c 0 ≈ 0 and μ ≈ 1/2 [54,62], meaning the yield stress τ * = 1 2 N. Assuming that the basal sediments are yielding, a reasonable assumption for the ice flowing at the speeds that West Antarctic ice streams do [14,58], then τ * = τ b = (1/2)N, where τ b is the basal shear stress. In this study, we impose that effective pressure at the bed is N 0 = 1, an assumption further discussed later in this section.
To estimate melting by geothermal heat flux, we take values of geothermal heat flux inferred from observational and field studies [55,63]. [63] estimated geothermal heat flux underneath Thwaites Glacier, another glacier in West Antarctica, using radar and found G ≈ 114 mW m −2 . [55] presented estimates of geothermal heat flux from geomagnetic data and borehole measurements and found that flux ranges from ≈60-100 mW m −2 in West Antarctica. However, there is still significant uncertainty in estimates of geothermal heat flux [64]. Therefore, we compute a range of meltwater flux from geothermal heating, assuming that geothermal heat flux falls between 20 and 120 mW m −2 . These values encompass most of the estimates of geothermal heat flux in non-volcanic regions of the Earth [65].
The contributions of geothermal heating, frictional heating from sliding at the bed and shear heating (estimated in this study with both a constant and linear vertical velocity) to meltwater flux are presented in figure 6. Contributions from all three sources are approximately of a similar magnitude (approx. 10 −4 -10 −2 m yr −1 ). In most cases, contributions from shear margins are larger than contributions from frictional heating at the bed (in the case of Bindschadler Ice Stream and Pine Island Glacier, shear margin contribution is approximately 1 − 2 orders of magnitude larger than contributions from frictional heating at the bed). Further, we are again not accounting for feedbacks between these processes and water content. Lubrication at the bed by liquid water may reduce frictional heating and thus reduce the contribution of frictional heating to meltwater production at the bed.
In Bindschadler Ice Stream, the contribution of shear margins decreases with distance from the grounding line, as strain-rates decrease. A similar trend occurs in Byrd Glacier, though far upstream (approx. 80 km), there is a region of high meltwater flux contribution from shear heating. In Pine Island Glacier, the contribution of shear heating is approximately two orders of magnitude larger than the contribution of frictional heating at the bed. Near the grounding line, meltwater flux decreases sharply. This decrease can be seen in estimates of porosity ( figure 4c), where there is a narrow band of reduced porosity near the grounding line. This is due to an anomalous, sudden decrease in ice thickness. The resolution of the surface velocity data (250 m × 250 m) is not fine enough to determine whether this reduction has a clear physical basis. In general, the linear vertical velocity produces a larger meltwater flux by approximately half to one order of magnitude than the constant vertical velocity case.
It is important to note that the contributions from frictional heating at the bed are strongly dependent upon the choice for the basal boundary condition for effective pressure N 0 . By choosing N 0 = 1, we are suggesting significantly less friction at the bed than there likely is in natural systems. Therefore, we would expect, in real ice streams, for the contribution from frictional heating to be at most an order of magnitude larger, though still of the same or lower magnitude than the contributions from shear margins estimated here. In regions where surface melting is significant (such as in the Greenland Ice Sheet), the contribution from surface melt would likely be dominant over the three sources presented here, though on the Antarctic Ice Sheet that contribution is likely small.
These results suggest that shear heating is a non-negligible contribution of meltwater to the bed in shear margins. Furthermore, these results support previous works that propose feedbacks between rapidly deforming ice and subglacial hydrology. Melting from shear margins drains to the bed and can provide a control on the width and position of ice streams, particularly in ice streams that are not topographically controlled [31,34,36,38]. Further, water content in subglacial till and in subglacial channels can affect the rate of ice flow over the bed by affecting the strength of subglacial till and thus the rate of basal sliding [25,34,38,[66][67][68]. These results also suggest that shear heating must be accounted for as a contribution to subglacial lakes and may be a way to model the volume of water in subglacial lakes. Improving the modelling of the contribution of water to subglacial lakes may be a step towards predicting changes to flow speed due to floods from subglacial lakes [17,20,69]. Finally, the localization of meltwater delivery in shear margins has important implications for the spatial variability of basal shear stress and the patterns of subglacial channels underneath Antarctic ice streams. The estimates provided here may provide a step towards illuminating the processes affecting basal properties, shear margin migration and freshwater input into the ocean. Further, these results can provide an input into subglacial routing models which may estimate where water drained from temperate ice flows and how much ends up in subglacial lakes. We reserve for future work an exploration of where the meltwater goes.

(b) Summary of key assumptions and future work
This study makes a few assumptions that impact these results, though we expect the main conclusions (that shear margins may provide a significant contribution to meltwater at the bed) to hold. Here, we discuss these model assumptions and how this impacts the analysis of the results presented in this study.

(i) Lateral advection
One of the most significant is the assumption that lateral advection of cold ice into the shear margins along the length of the shear margins is negligible. While this assumption allows us to explore meltwater flux estimates in natural systems in a simple enough framework to understand the dependencies of estimates on physical processes, it also affects the estimates presented in this study. In particular, a number of studies have suggested that the lateral advection of cold ice into the margin may reduce ice temperature, counteracting some of the deformational heating and thus reducing thickness of temperate zones and the flux of meltwater [25,31,33,[70][71][72][73]. Further, lateral advection of nonporous ice may also affect estimates of ice porosity. Lateral advection of cold ice is thus an important physical process to capture in order to generate fully reliable estimates of temperature profiles and meltwater fluxes. This is particularly relevant in ice streams that are not topographically confined, such as Bindschadler and MacAyeal Ice Streams and Pine Island Glacier. This is distinct from Byrd Glacier, which sits in the Trans-Antarctic Mountains and thus has minimal rates of lateral advection. While we leave for future work the incorporation of lateral advection into this model formulation, there are a number of studies that provide steps towards incorporating horizontal advective terms. As a few examples: Meyer & Minchew [37] provide a parametrization in the one-dimensional set-up, in which they define lateral advection as a correction to the Brinkman number. Suckale et al. [33] present a model that resolves a twodimensional section across an ice stream margin and parametrizes the effect of lateral advection on ice temperature, though not on mass conservation. Haseloff et al. [71] explicitly incorporate the lateral advection of cold ice into a thermomechanical model that accounts for its effects on ice temperature as well as mass conservation. Based on these studies, it is plausible that the temperate ice zones are overestimated in this study (and may, in fact, not be present in the shear margins of interest). Incorporating a framework for lateral advection into the model presented here is a necessary next step to fully realizing the presence of temperate zones and ice temperature profiles in Antarctica.
(ii) Basal boundary condition for effective pressure, N 0 The boundary layer solution for porosity applied in this study and initially derived and presented by Schoof & Hewitt [39] is valid for N 0 δ −(1/2) . In this study, we use a value of N 0 = 1, which translates to dimensional effective pressure at the bed of approximately 50-1000 Pa, depending on the ice thickness in the ice stream.
Since N 0 = 1 implies a dimensional effective pressure at the bed of approximately 50-1000 Pa, this implies that the basal shear stress τ b ∼ 25-500 Pa, based on the assumption of a perfectly plastic bed as discussed previously. This does align with some in situ measurements of basal effective pressure [74]. However, this is significantly lower than basal shear stress estimates found through inverse methods. Basal shear stress in ice streams (particularly West Antarctic ice streams) is estimated to be τ b ∼ 10 3 -10 4 Pa [5][6][7]75], meaning that our assumption here is lower than these estimates from observations by approximately 1 − 2 orders of magnitude. Larger effective pressure at the bed would likely cause ice compaction at the bed, impeding the draining of water to the bed and reducing the flux estimates presented here. Accommodating basal boundary conditions for effective pressure greater than that used in this study would require the derivation of a different boundary layer solution, and this would be a fruitful direction for future work to build on this model and improve the accuracy of the estimates presented here.

(iii) Porosity and effective pressure model parameters
The estimates presented here are affected by uncertain parameters, in particular the porosity exponent α, the non-dimensional ratio of heat advected in meltwater to heat conducted in ice κ, the non-dimensional number representing the thickness of the basal boundary layer δ, and the non-dimensional numbers that affect rates of shear heating and its effect on ice temperature Br, Pe.
The parameters κ and α represent the permeability of glacier ice. Their exact values are uncertain. The porosity exponent α likely lies between α = 2 and α = 3, consistent with a Kozeny-Carman model of permeability [40,76]. We define α = 2.33, in line with [39]. We compute κ from equation (2.12), and the specific values are presented in the electronic supplementary material. The effect of varying these parameters was explored in [39]. Larger values of κ and α result in somewhat larger meltwater flux estimates, though the flux at the bed is not very sensitive to either parameter and is more sensitive to κ than to α. The depth structure of flux is affected by the choice of κ.
The non-dimensional number δ physically represents the contribution of compaction to meltwater flux, and it directly relates to the thickness of the basal boundary layer by δ (1/2)z being the non-dimensional height of the boundary layer. This connection arises due to the fact that, within the boundary layer, pressure gradients are the primary control on meltwater flux, while outside of the boundary layer, meltwater flux is primarily driven by gravity. The parameter δ controls the height at which pressure gradients control meltwater flux. In the (iv) Flow-rate parameter, A We assume a constant ice softness parameter A and take a simple, isotropic form of Glen's Law that computes ice viscosity based solely on stress and strain rate. However, factors such as anisotropy (fabric), ice temperature and ice porosity likely affect ice viscosity. Incorporating the effect of anisotropy, ice temperature and ice porosity in ice softness will likely increase the incidence of temperate zones and increase meltwater flux from temperate ice in shear margins. Currently, anisotropy is parametrized through fabric enhancement factors in the flow law, though the values of these enhancement factors are uncertain. Therefore, we keep enhancement factors implicit within the value of A in this study in order to reduce uncertain parameters in the model and we reserve for future work a more complete study of the effect of anisotropy on rates of deformation and, therefore, on the generation of temperate zones. Given that ice in shear margins is approximately in simple shear, many studies suggest that a scalar enhancement to the ice softness parameter A is appropriate for these conditions [77][78][79]. Any enhancement to strain rates not captured in our constant value for A would likely increase ice temperatures from a given strain rate and therefore increase the thickness of existing temperate zones and create temperate zones where we do not estimate any in this study.
Including ice temperature in the ice softness parameter A through an Arrhenius relation has been done in previous studies, and these studies have shown that the generation of temperate ice from shear margins is likely [25,33,38]. Incorporating temperature into A may also produce larger and more extensive temperate zones than we present here, because shear heating would soften the ice and enable faster deformation, thereby increasing the Brinkman number Br and directly increasing the thickness of temperate zones. In the electronic supplementary material, Section E, we show that regions and extent of temperate ice zones estimated are larger when we allow for coupling between A and ice temperature, suggesting larger meltwater fluxes than those presented here.
A similar effect occurs when incorporating ice porosity into ice softness A for regions of temperate ice. The effect of porosity on ice softness is generally parameterized by A = A 0 (1 + cφ) where c ≈ 200 [46,80], and therefore increasing porosity in temperate zones, as seen in this study, would result in an increase in ice softness [25]. Softer ice would deform faster, increasing the rate of shear heating and thus increasing the presence of temperate zones. Given that, in this study, we use a constant ice softness A, the estimates presented in this study can be thought of as a lower bound on meltwater flux from shear margins in the absence of horizontal advection, and we reserve for future work a full consideration of the feedbacks involved with porosity, ice temperature and ice softness.
(v) Stress exponent, n In this study, we take an isotropic form of Glen's flow Law with the stress exponent n = 3. While this value is commonly used in glaciological literature and is supported by laboratory experiments [45,46], recent studies have suggested that in ice sheet conditions, the value of the stress exponent n may be closer to n = 4 [81][82][83]. Ice rheology affects these results through estimates of ice temperature and therefore estimates of the thickness of temperate zones. Larger values of n allow for higher rates of shear heating and thus would produce more significant temperate zones. This would increase our estimates of meltwater flux to the bed and extend the regions in which considering polythermal glacier structures is necessary. Further, there has been some consideration for the effect of ice temperature and ice deformation rate on the stress exponent n. In particular, studies have considered how the stress exponent n changes with liquid water content in the ice and have found significant softening and acceleration due to the presence of liquid water [21,80,84] and potential changes in the deformation mechanism, resulting in estimates of n as low as n = 1.1 for large water contents [22,24].

Conclusion
In this study, we quantify the contribution of shear heating in the margins of ice streams to the flux of meltwater to the bed, which often flows in channels that transport water to the ocean and affect the rates of ice flow and the geometry of ice streams. This work has two aims. First, we demonstrate that the contribution of shear heating to melt at the bed of ice sheets is comparable to other mechanisms that are currently accounted for in models (such as geothermal heating and frictional sliding at the bed). Further, the meltwater supply from shear margins is highly localized, and previous studies have suggested that the localized increase in water content at the bed may influence yielding of the bed and the location of ice stream margins, which has implications for large-scale ice stream flow [25]. Second, we show that the parameterization of vertical velocity partially controls the estimates of meltwater flux at the bed, suggesting that a better understanding of the structure of vertical advection is necessary to better constrain estimates of meltwater flux.
In applying this analytical model, we made a number of assumptions, the most significant of which are assuming a constant ice softness parameter A, neglecting lateral advection, and applying a vertical velocity that is constant with depth. A necessary direction of future work is to incorporate the effects of ice temperature, porosity, and anisotropy into ice rheology. This would likely increase the estimates of meltwater flux due to increasing the rate of shear heating in softer ice and increasing the thickness of temperate ice zones, as shown in the electronic supplementary material. Another important direction is the incorporation of lateral advection, which may reduce estimates of meltwater flux due to the advection of cold ice into the shear margin. We have shown some preliminary results, based on models from Schoof & Hewitt [39], for the application of a linear vertical velocity profile, and future work will evaluate the effects of vertical and longitudinal advective profiles as well. Finally, a useful direction for future work is to consider time-dependence and the evolution of meltwater flux.
The results presented here suggest a need to incorporate models for polythermal ice flow, as done in [39,40], and to constrain ice flow properties in temperate ice. We note that the application of asymptotics provides a higher-order correction on meltwater estimates, which likely has less of an effect on meltwater estimates than data uncertainties and the simplification of the model itself. While we still show full composite solutions to present estimates accurate to a higher-order, we note that in further applications, the outer solutions may provide good approximations to the full solution. These results may help quantify the feedbacks between rapid, localized deformation and basal properties and the effects of these feedbacks on ice flow. Finally, the estimates presented here may shed light on the amount of meltwater in subglacial channels and subglacial lakes, which is currently uncertain and has the potential impact ice flow significantly.
Data accessibility. The source code for the model presented in this study, along with the data needed to generate figure 4, are openly available at https://github.com/megr090/MeltwaterFlux. No new data were produced for this study, and data used in this study are publicly available through their respective publications.
The data are provided in electronic supplementary material [85].