1. Introduction
Numerical simulations of the flow of shear-thinning ice over the complex topography of ice-sheet beds using 3-dimensional full-Stokes models are computationally expensive, but in a wide variety of scenarios the geometry of the situation allows for simplified calculations to capture the dynamics (Hindmarsh, Reference Hindmarsh2004). Ice is often much shallower than the horizontal scales on which driving stresses and average basal conditions vary, motivating shallow approximations where pressure gradients are hydrostatic (Cuffey and Paterson, Reference Cuffey and Paterson2010; Fowler, Reference Fowler2011). If the ice is frozen to the bed and vertical shear dominates the deformation, the shallow ice approximation (SIA) accurately captures the dynamics (Raymond, Reference Raymond1996). In fast-flowing regions of ice sheets where basal resistance is low, the dominant stresses are a combination of extensional stress and lateral shear, so the shallow shelf approximation (SSA) or depth-integrated membrane models provide lower-dimensional formulations (Morland, Reference Morland1987; MacAyeal, Reference MacAyeal1989; Blatter, Reference Blatter1995; Pattyn, Reference Pattyn2003; Hindmarsh, Reference Hindmarsh2004; Schoof and Hindmarsh, Reference Schoof and Hindmarsh2010).
In simple geometries, such as wide, uniform ice streams, a further advantage of the SIA and SSA models is that they can be solved exactly to give closed-form expressions for the ice velocity and flux (Raymond, Reference Raymond1996). These expressions can then be used as a simple, essentially 1D representation of ice dynamics in models of related dynamical processes, such as the response of tidewater glaciers to calving (Benn and others, Reference Benn, Hulton and Mottram2007), the impact of basal melt during Heinrich events (Mann and others, Reference Mann, Robel and Meyer2021) or interactions of the ice with its bed at the start of surges (Minchew and Meyer, Reference Minchew and Meyer2020). However, the simplified expressions for the velocity and flux of ice have often been used beyond the regimes in which they are asymptotically valid, in particular in the transition from SSA to SIA as basal drag increases. Simply summing the two expressions (e.g. Mann and others, Reference Mann, Robel and Meyer2021) as an ad-hoc transition does not correctly capture the particular dynamics occurring in this rapidly sliding, high basal shear regime, as we show numerically and mathematically in this work. Furthermore, by understanding the transition between these flow regimes, the detailed flow profiles may be used to infer basal properties.
Small amounts of basal drag, such as would be experienced by an ice stream sliding over a bed of yielding sediment, can be viewed as a perturbation to the SSA expression for freely sliding flow. A natural question is at what point the flow becomes significantly modified, and whether this happens before reaching the fully non-sliding regime. Schoof and Hindmarsh (Reference Schoof and Hindmarsh2010) mathematically examined the effect of basal drag as a higher-order correction to purely extensional flow. Here, motivated by the geometry of ice streams and their importance in wider models of ice sheet and climate dynamics, we look instead at the case where the dominant resistance to flow comes from lateral shear at the sidewalls, and the rheology of the shear margins is the primary control on the velocity of the stream (Meyer and Minchew, Reference Meyer and Minchew2018; Hunter and others, Reference Hunter, Meyer, Minchew, Haseloff and Rempel2021).
Since ice is shear-thinning, basal drag impacts not just the overall force balance but also the depth-averaged rheology. This motivates the ‘L1L2’ modelling framework (Hindmarsh, Reference Hindmarsh2004) in which a linear approximation to the vertical shear acts to modify the ice viscosity. This framework has been implemented in large-scale ice-sheet models, including BISICLES, to capture ice flow dynamics lying between the SSA and SIA regimes (e.g. Goldberg and Sergienko, Reference Goldberg and Sergienko2011; Cornford and others, Reference Cornford2013). The hybrid model captures the transition between SSA and SIA in a less computationally intensive way than full-Stokes simulations, but remains a predominantly numerical technique. Here we mathematically analyse the ‘L1L2’ method to derive simple, yet more accurate expressions for the ice velocity and flux in this transitional regime.
In this paper, we look at the flow of shear-thinning ice over a bed of uniform shear strength, through a wide, shallow rectangular stream. We use asymptotic analysis to derive closed-form expressions for the velocity field and flux through the stream, both when basal resistance is small and as it increases towards the driving stress. We show that these expressions recover both the SSA and SIA limits as expected, but we also obtain a new, intermediate regime in which basal shear-softening is significant. We further solve for the flow-field numerically, confirming the existence of this regime, and show that our new expression significantly reduces the error in predicted velocities.
2. Governing equations and approach
We consider the steady flow of shear-thinning ice in a shallow, rectangular channel of depth H and half-width W which is uniform in the along-flow direction, with equal bed and surface slopes (Fig. 1). This means the pressure is hydrostatic, and the velocity field is only in the along-stream direction, u = (u(y, z), 0, 0). As the flow is uniform in this direction, we have no extensional strains in this geometry, a simplification that allows us to retain and evaluate the impact of vertical strain instead. While longitudinal and vertical strains can become comparable if basal resistance is very low and the along-stream variation is large, along-stream changes in velocity often occur over length scales much greater than the width or depth of ice streams, with extensional stresses effectively acting to smooth variations in driving stress (Kamb and Echelmeyer, Reference Kamb and Echelmeyer1986).
Diagram of the ice stream geometry, showing the u = 0 boundary condition on the sidewalls y = ±W, the surface velocity profile u s(y) at z = H, the profile of sliding speed u b(y) at z = 0, and the point at y = y u marking the changeover in basal boundary condition from τ b = μN in the central region to u b = 0 close to the sidewalls. An example vertical velocity profile is drawn to show the concentration of vertical shear close to the bed.
Given these assumptions, the stress field τ within the ice satisfies
where the driving stress,
is due to gradients in surface height S, density of ice ρ i and gravitational acceleration g. We take a power-law shear-thinning rheology for the ice
where A is the viscosity prefactor, which in general may be a function of temperature and grain size, but here we take as constant. The shear-thinning exponent n is usually taken to be 3 for ice (Glen, Reference Glen1955), although using values from 1.8 to 4 represents the different modes of deformation that make up ice flow (Goldsby and Kohlstedt, Reference Goldsby and Kohlstedt2001; Millstein and others, Reference Millstein, Minchew and Pegler2022). Even higher effective shear exponents could represent the effect of damage or shear heating in the margins of the ice stream, where high shear rates warm and soften the ice further (Minchew and others, Reference Minchew, Simons, Riel and Milillo2017). We will show numerical results for n = 3, but give expressions for general n as far as practical.
2.1. Boundary conditions
On the two sidewalls at y = ±W we apply a no-slip boundary condition, motivated by the fact that ice streams and glaciers are usually either topographically constrained between high bedrock or surrounded by nearly stagnant ice. At the upper surface of the ice, z = H, we impose a no-stress condition as the ice is in contact with the atmosphere. At the base of the stream z = 0, we link the basal shear stress, τ = τ b, to the basal sliding speed, u = u b, through a sliding law representing the rheology and/or geometry of the bed in contact with the base of the ice stream.
Over beds of deformable sediment, an appropriate basal boundary condition is a plastic traction law,
where μ is a friction coefficient and N = p i − p w is the effective pressure, which is the difference between the pressure exerted by the weight of the ice p i and the water pressure at the sliding surface p w. This will form the basis for our analysis here, since deformable sediments are found below many fast-flowing streams (Kamb, Reference Kamb, Alley and Bindschadler2001; Iverson, Reference Iverson2010).
In large-scale models, the basal shear stress and sliding speed are often related via a power-law,
originally derived from ice sliding over a hard bed with small-scale topography (Weertman, Reference Weertman1957) but presently extended to other scenarios by modifying the value of the exponent m. Large values of m ≥ 8 have been used to approximate plastic deformation (Rosier and others, Reference Rosier, Gudmundsson and Green2015; Minchew and others, Reference Minchew2016; Joughin and others, Reference Joughin, Smith and Schoof2019). A numerical study of the impact of this power-law sliding in a similar channel geometry is given in Adhikari and Marshall (Reference Adhikari and Marshall2012), and we will discuss the value of our plastic bed calculations for reproducing these results.
If the basal shear strength is large compared to the driving stress τ d, (4) reduces to a bed that is unyielded everywhere, and we obtain a no-slip boundary condition, u b = 0, as would be the case in the SIA. When τ b → 0, we recover a no-stress basal boundary condition, where the SSA is strictly valid. In this paper, we recover both of these well-known limits but focus on the intermediate regime where both u b/u s (u s the surface speed) and τ b/τ d are large. We seek to find closed-form approximations for the ice velocity field, in particular the surface speed u s (which is most readily observable), basal sliding u b and total ice flux Q, as a function of the width W and height H of the stream, the shear thinning exponent n and the ratio of basal shear stress to driving stress τ b/τ d.
Our reference point will be the expressions for centreline surface velocity u mid and flux Q obtained by simply summing the SIA and SSA results for flow through a uniform ice stream (Raymond, Reference Raymond1996; Cuffey and Paterson, Reference Cuffey and Paterson2010),
which neglect the effect of basal shear softening. Our goal is to improve on the accuracy of these results while still retaining a readily useable closed-form expression.
3. Numerical method
We numerically solve for the flow of ice in a rectangular channel with variable basal traction to assess the accuracy of our asymptotic expressions. We use COMSOL to solve (1) in the half-channel 0 < y < W, 0 < z < H, invoking the symmetry of the flow, as a generalised Poisson equation for u. The effective viscosity is regularised by the addition of a small constant $\epsilon _1 = 10^{-9}$
to prevent divergence at z = H, y = 0,
The boundary conditions are no-stress on y = 0 and z = H, no-slip on y = W and a traction law on z = 0 that represents a regularised form of (4),
We used an automatically generated triangular mesh on COMSOL's extra-fine setting (5 grid points across the depth of the stream), and with relative error for convergence of 10−5. The regularisation parameters $\epsilon _1 = 10^{-9}$
and $\epsilon _2 = 10^{-4}$
were also chosen such that decreasing their value by an order of magnitude did not affect the calculated velocities by more than a factor of 10−5. Using the finest available mesh setting also did not affect calculated velocities by more than a factor of 2 × 10−4 as compared to the extra-fine setting.
The numerical results are not the focus of this work but serve as a benchmark for our derived expressions for velocity and flux, as described in the following sections.
4. Mathematical method
To begin our mathematical analysis, we start by depth-integrating (1) to get
where τ b appears as the integration constant found by evaluating this expression at z = 0. Evaluating (10) at z = H, where τ xz = 0, we find that
Where the bed is yielded, so that τ b = μN uniformly, we can integrate across the stream to arrive at
and so the problem reduces to evaluating the depth-integrated lateral shear stress, and in particular evaluating the relative effects of vertical and lateral shear on the viscosity.
4.1. Approximation of vertical shear stress and viscosity
If basal sliding is significant across the width of the channel, we assume that gradients in the surface velocity profile u s(y) provide an accurate reflection of the lateral shear rate throughout the full depth of the ice.
Motivated by numerically calculated stress fields (Fig. 2), the much larger horizontal than vertical lengthscales, and the approach of Hindmarsh (Reference Hindmarsh2004), we further assume that the vertical shear stress that feeds into the rheology is approximately linear with depth. Hence we may approximate the lateral and vertical shear stresses, and hence the ice rheology, as
where the tilde denotes that this is now based on the approximate vertical shear stress, so e.g. $\tilde {\eta }$
is the approximate ice viscosity. This is similar to the L1L2 model of Hindmarsh (Reference Hindmarsh2004) but here vertical shear stresses are scaled with μN rather than τ d.
Left, numerically calculated τ xz, and right, difference between τ xz and the approximation $\tilde {\tau }_{xz}$
(13), for increasing values of 1 − μN/τ d = 10−2.5, 10−2, 10−1.5, 10−1, 10−0.5, 1. All values are scaled so that τ d = 1. The agreement is generally excellent, but τ xz is noticeably less than $\tilde {\tau }_{xz}$
in a region O(H) near the sidewall, and over unyielded regions of the bed where τ b < μN. Near the centre of the stream the divergence of viscosity leads to error in the numerically calculated field (particularly visible when μN = 0 and the exact solution is τ xz = 0 everywhere). All for W/H = 10 and n = 3.
4.2. Horizontal shear rate
Writing $\tilde {\eta } = \tau _{xy}/{\partial u_s\over \partial y},\;$
we can approximate Eqns. (13)–(14) in two limits depending on the relative sizes of τ xy and $\tilde {\tau }_{xz}$
. If vertical shear stresses are much greater than horizontal, $\tilde {\tau }_{xz}\gg \tau _{xy},\;$
which exhibits strong depth-dependence. In contrast, if lateral shear stresses are much greater than vertical, $\tau _{xy}\gg \tilde {\tau }_{xz},\;$
which is uniform in depth.
We define a transition depth z t at which a transition in the dominant stress occurs, and below which basal shear dominates, given by
Note that in regions where the lateral stress is large, basal shear never dominates and z t would be negative. Thus, we have two regimes depending on the sign of z t. Where lateral stresses are smaller than μN, the leading order expression for the depth-integrated stress is given by integrating the dominant term on either side of z t, so starting from (12) we have
while if z t ≤ 0 then lateral shear dominates the viscosity everywhere, so
If μN ≪ τ d, then (23) applies almost everywhere and the leading order expression for the surface velocity profile is given by
which is the SSA approximation for an ice stream with constant basal traction μN (Raymond, Reference Raymond1996). While the modification to τ d is negligible in the limit for which (24) is asymptotically valid, this expression is used even when μN ~ τ d (e.g. Minchew and others, Reference Minchew, Meyer, Robel, Gudmundsson and Simons2018).
However, as basal drag increases, the vertical shear stress starts to dominate in the central region of the stream where lateral stresses are lowest. This shear-softening decreases the depth-integrated viscosity, and leads to a significant increase in velocity compared to the prediction of the SSA approach when (τ d − μN)W ~ τ dH (Fig. 3).
Including shear softening improves the match to numerically calculated values of (a) the centreline surface velocity u s(0) and (b) the total flux Q, particularly when τ d − μN ~ 0.1τ d. Here W = 10H and n = 3. SSA and SIA limits of the surface velocity are shown as dashed and dotted lines.
The SSA result becomes an even poorer approximation as basal drag increases further, and the bed becomes unyielded over a significant portion of the stream. Using (12), which assumes that τ b = μN everywhere, overestimates the resistance to flow and thus additionally underestimates the ice velocity; to correct this we must look further at the vertical structure of the flow.
4.3. Vertical shear rate and basal boundary condition
Now that we have a measure of both τ xy and τ xz, we can calculate the effective viscosity and thus evaluate the vertical deformation of the ice by finding ∂u/∂z. This allows us to estimate the basal velocity u b and to better characterise whether the bed should be yielded (τ b = μN) or unyielded (u b = 0).
At the sides of the stream, where lateral shear dominates the viscosity, and still approximating the vertical shear stress as linear in depth, we can use (13) solve for the approximate vertical shear,
and then integrating this vertical shear rate down from the surface we obtain a velocity field
resulting in a basal sliding speed
Note that this becomes negative very close to y = W; in effect this determines the point beyond which the bed is unyielded and u b = 0 should apply instead of τ b = μN. We can find the approximate edge of the unyielded region y = y u by solving u b(y u) = 0, so that
To leading order in W − y u we can solve for y u approximately and get
which, for its simplicity, we will take through the rest of the paper in all calculations and plots using y u. Figure 4 shows that (29) accurately captures the behaviour as y u → W, as well as agreeing that y u → 0 as μN → τ d, but determining the exact boundary of the yielded bed is a free boundary problem (Schoof, Reference Schoof2006) and in general can only be solved for numerically, and hence expressions left in terms of y u could be improved with higher-order calculations of this yield boundary's location.
Match between the approximate expression for y u (29), and the numerically calculated edge of the yielded bed, for W = 10H, n = 3. The disagreement at larger μN is expected since when y u moves away from the sidewalls the approximations that lead to (28) no longer hold, but (29) does capture that y u → 0 there, in agreement with numerical results.
To find the internal deformation in the centre of the stream, where basal shear begins to affect the ice viscosity, there are two regimes to consider: below z t the basal shear dominates, while close to the surface the lateral shear sets the viscosity. Considering only these dominant contributions to the viscosity, we can obtain the leading order approximation to ∂u/∂z above z t by
so that the velocity field, found by integrating down from the surface profile, is given by
in particular at z = z t where
To continue down below z t, where the leading order form of the viscosity is dominated by basal drag, we similarly calculate the vertical shear as
then integrate down from z = z t to find
Thus, the basal velocity is
Note that the deformation velocity, u d = u s − u b is
The prefactor is the SIA surface speed, and the second term inside the brackets is much less than 1 until the very edges of this central, basal drag dominated region. Similarly the deformation flux q d, the integral of u(z) − u b through the depth of the stream, is dominated by the SIA term so
5. Results
5.1. Surface velocity and sliding speed
Surface velocity is the primary observation available to constrain ice dynamics and basal processes, while sliding speed is a key determinant of the total flux through the ice stream. We now combine our expressions for stress and velocity to obtain the surface velocity at the centreline of the ice stream. We have calculated the horizontal shear rate over the yielded regions where τ b = μN, the extent y u of the unyielded region of the bed, and the vertical shear at the edge of this region, where both u b = 0 and τ b = μN. Thus, we integrate up from the bed at y u to find u s(y u), then integrate inwards to find the surface velocity profile over the yielding region, giving
There are two regimes to consider depending on the dominant contribution to the depth-integrated viscosity at y u. When the basal drag is small, (23) dominates and we have a small modification to the SSA solution,
Since y u is very close to W, and τ d − μN ≫ μN, the deformation velocity is small compared to the horizontal shear and so that the leading order expression for the centreline velocity is exactly the SSA result of
and this is also the leading order term in u b.
When the basal drag is large, (21) dominates the horizontal shear, and so
and thus
Using the expression for y u from (29), we find that (42) provides a good match to the numerical results when τ d − μN ≪ τ d (Fig. 5).
Fractional error in expressions for u mid as compared to the numerically calculated values (1 – expression/numerical value, so overestimates are negative errors), for a range of aspect ratios and basal strengths, with n = 3. (a) Sum of SIA and SSA, Eqn. (6), (b) including basal shear softening, but with y u set to W (45), and (c) the full expression (43) using (29) for y u. Both (b) and (c) represent improvements over (a), with (c) the closest match to the numerical results. When the aspect ratio is not large, all the approximations start to break down. The maximum error in (b) is −0.40 at W/H = 4 but reduces to −0.17 at W/H = 11. The error in (c) ranges from −0.24 to 0.085 but remains within ±0.1 for W/H > 5.75.
We suggest that an appropriate expression for the centreline velocity that transitions smoothly between the two regimes can be found by adding together (40) and (42), giving
We see that the first term, which dominates when μN ≪ τ d, is the SSA expression, while the final term, which remains when τ d = μN, is the SIA expression. In effect, we have introduced a new term to represent the intermediate regime, in which basal drag is significant enough to affect the viscosity, while basal sliding remains high enough to affect the surface velocity. Using (29) as the approximation for y u, this becomes a closed-form expression for the surface velocity as a function of basal strength which could be readily used as a lumped description of ice stream dynamics in a model of wider climate processes. The corresponding expression for the sliding speed at the centre of the stream is found by subtracting off the deformation velocity,
We also consider a further simplification to (43) by approximating y u as W,
This will systematically overestimate the surface velocity, but in wide channels, this can still provide an improved match to the numerical results compared to neglecting the basal shear-thinning entirely (Fig. 5).
Figure 5 shows the difference between the numerically calculated values of u mid and the three approximations: SIA plus SSA, (45) and (43). SIA plus SSA significantly underestimates the centreline velocity for wide streams with intermediate sliding rates – clearly showing that SSA is valid when μN/τ d ≪ 1, while SIA is valid when μN/τ d ≫ H/W, missing the intermediate regime if W ≫ H, where the errors can be upwards of 40%. As expected, (45) is an overestimate but increasingly close to the numerical results as W/H increases. (43) provides the closest match to the numerical results, remaining within 10% of the calculated values and only significantly deviating where W/H < 5 and the shallow approximations break down.
5.2. Total flux
While surface velocity is most easily measurable, the evolution of ice stream thickness depends on the total flux through the stream. We now calculate the ice flux through the stream, again considering the two regimes depending on the size of μN compared to τ d.
When μN is small, vertical deformation is negligible, and so the total flux Q is found by integrating u sH across the width of the stream, with u s given by (39). Thus, we recover the SSA result,
However, if the bed is strong, we divide the flux between the region over the yielded bed, and the region over the unyielded margins. Over the central, yielded region, we have contributions to the flux from both the basal sliding (41), and the internal deformation (37), such that
In the margins, we have u b = 0 and large vertical shear stress. While lateral shear will start to dominate in a region of order H close to the sidewalls, until then the SIA holds and the velocity is approximately y-independent, with
Thus, the flux in this outer region is approximately given by
where the −αH represents the reduced flux due to matching onto the no-slip sidewalls. On dimensional grounds, we anticipate that the sidewall correction should tend to a constant multiple of H as W increases. Numerically, it appears that α ≈ 1.4 (Fig. 6). While in the asymptotic limit of W ≫ H, this correction is negligible, for moderately wide channels the inclusion of this term provides a simple increase in the accuracy of our expression (Fig. 7). Further analytic progress in this regime is beyond the scope of this paper, and further numerical calculations of the non-sliding regime (taking instead a shape factor-based approach to the effect of the sidewalls) are found in Nye (Reference Nye1965). For wide channels, this is a small correction, and we are primarily interested in the more rapidly sliding regime.
With a no-slip base, the effect of the sidewalls on the flux as summarised by the value of α in Eqn. (49) tends towards a constant value as W/H increases.
Fractional error in expressions for Q as compared to the numerically calculated values, for a range of aspect ratios and basal strengths, with n = 3. (a) Sum of SIA and SSA, Eqn. (7) – note this does not include the sidewall modification to the SIA flux, (b) including basal shear softening, but with y u set to W (52), and (c) the full expression (51) using (29) for y u. Both (b) and (c) represent improvements over (a), with (c) the closest match to the numerical results. To guide the eye, the fractional error in (b) ranges from −0.29 (at the smallest values of W/H) to 0.017 and in (c) from −0.036 to 0.098.
The total flux through the ice stream in the high basal drag regime therefore is approximately given by
Again we suggest that summing the fluxes derived from each limit gives a single expression that transitions accurately between regimes, so
As shown in Figure 7, this expression is within 10% of the numerically calculated results, far better than just summing the SSA and SIA values.
As before, a simpler version using y u ≈ W is given by
and can also provide a good approximation to the numerical results, particularly for wide streams (Fig. 7b). However, it is always an overestimate. By contrast, neglecting basal shear-thinning by setting y u = 0 and simply summing the SSA and SIA expressions is an underestimate that captures the intermediate regime increasingly poorly for wider streams (Fig. 7a).
5.3. Cross-stream surface velocity profiles
The surface velocity profile reveals more information about the bed and ice rheology than the centre-line speed alone. In our numerical results, we see appreciable changes in surface velocity profile as the bed strength increases (Fig. 8), with the width of the central region of quasi-uniform velocity first narrowing in the intermediate regime and then increasing to almost the full width of the stream in the SIA limit.
Numerically calculated profiles of surface velocity for increasing values of 1 − μN/τ d, compared to the SSA profile (39), the numerically calculated surface velocity with u b = 0, and the model we propose in (53), taking into account basal shear softening. The velocities are scaled so that the SIA limit is u mid = 1.
As previously noted, a calculation of the velocity field close to the sidewalls, over the unyielded sections of the bed, is beyond the scope of our present analysis (see also the discussion section). Thus, to explain the full surface velocity profiles including the shear margins, we rely on the numerically calculated profiles for the u b = 0 limit.
Initially, as the basal drag increases, basal shear-softening reduces the effective viscosity in the centre of the stream, leading to larger shear rates there, and more rounded velocity profiles (u mid − u s ~ y (n+2)/2) compared to the freely sliding case (u mid − u s ~ y n+1). As bed strength increases further, the unyielded regions of the bed encroach towards the centre of the stream. The shear margins are concentrated over the non-sliding regions (c.f. Truffer and others, Reference Truffer, Echelmeyer and Harrison2001), and numerically we see that the shape of these margins closely resembles surface velocity profile for a completely static bed. We therefore suggest that an approximate expression for the surface velocity profile over the entire stream can be found, similarly to the method by which the centre-line velocity was found in (43), by summing the SSA solution, the numerically calculated profile for non-sliding bed, and a term representing effect of basal shear softening where y < y u (using (29) for y u),
The predictions of (53) give good agreement with the numerical velocity profiles across the full width of the stream and throughout the full range of bed strengths (Fig. 8).
Simple expressions for surface velocity have been used to estimate ice depth or basal properties (e.g. Li and others, Reference Li, Ng, Li, Qin and Cheng2012) before recourse to full-Stokes modelling, and our improved estimate of surface velocity could be used to improve the accuracy of such work. However, with uncertainties in ice rheology, bed geometry and basal conditions, it is not possible to simultaneously invert for all unknowns from only a single value, namely maximum surface velocity. In large-scale inversions, the spatial heterogeneity of surface velocity is used as a constraint. In a similar manner, the width of the shear margins and the degree of uniformity across the centre of the stream could potentially provide sufficient constraints to simultaneously constrain ice rheology and basal strength. This could be used, for example, to extend the work of Millstein and others (Reference Millstein, Minchew and Pegler2022), approximating ice rheology, to grounded ice streams, or to re-examine Minchew and others (Reference Minchew, Meyer, Robel, Gudmundsson and Simons2018) beyond the use of the SSA approximation. However, the heterogeneity of real ice stream beds will complicate this kind of analysis.
6. Discussion
We have seen that for ice flowing in a stream over a plastic bed, there are three regimes of behaviour for the ice flow depending on the yield strength of the bed as compared to the driving stress. If the bed exerts very little resistance on the ice, the flow is primarily uniform with depth and resisted by lateral shear from the sidewalls (the SSA regime). If the bed strength is high and the ice cannot slide, vertical shear alone sets the flow field throughout the majority of the ice stream (the SIA regime). Between these limits is a regime where sliding is rapid and lateral stresses are still the primary control on ice speed, but basal drag causes shear of significant enough magnitude to soften the ice, thereby altering the response of the ice stream away from the margnis. We can summarise the effect of this softening through a new term in the expressions for surface velocity (43) and flux (51).
6.1. Ice-stream geometry and ice rheology
We have focused on a rectangular channel geometry which permits several simplifications in our approach. Firstly, we take the depth of the ice to be uniform, which considerably simplifies the integration of the shear profile across the width of the channel. While a cross-stream variation in ice depth would not alter the approach significantly, it would result in more involved analytical calculations. It would also be more complex to apply the traction condition on a sloping boundary, although if this slope is small then to leading order this correction will be linear in the slope angle. Secondly, in our model setup there is a clear separation between the bed of the stream and the sidewalls, allowing us to cleanly impose different boundary conditions on each surface. While subglacial sediments may be concentrated in the deepest parts of the glacier bed, with sidewalls sloping away from the bed (e.g. Truffer and others, Reference Truffer, Echelmeyer and Harrison2001), many ice streams in Antarctica are not confined by topography but found between regions of slower-flowing ice. In fact, one might consider these streams, generated by variations in bed strength, as an extreme example of a very wide channel in which only a narrow central section of the bed has yielded, and any actual sidewalls play a negligible role. Interestingly, this suggests that basal stress-controlled streams could show a different dependence of speed on their width and bed strength than topographically controlled streams – although the controls on the width of the yielded bed are very different from those we consider here. The flow of ice over a flat bed with variable bed strength is considered in Schoof (Reference Schoof2004) – the self-consistent determination of the yielded regions of the bed are key to solving for the stress field in this geometry.
We have taken our flow to be uniform in the along-flow direction, which reduces the dimensionality of the problem and removes any extensional stresses. This simplifies the expression for the effective viscosity, as we have assumed either basal shear or lateral shear must dominate in any region of the stream. In fact, at the centre of an ice stream, extensional stresses will dominate if along-flow velocity variations over a distance comparable to the stream width are large relative to maximum flow speeds. Variation on this magnitude is uncommon inland but occurs as ice streams accelerate towards their grounding lines (Mouginot and others, Reference Mouginot, Rignot, Scheuchl and Millan2017). Including these extensional stresses also introduces longitudinal terms in the force balance. We suggest that a comparison of extensional to lateral stresses in a freely sliding regime is better left to two-dimensional membrane models. Once the dominant horizontal stress relevant to the geometry is identified, basal drag can be introduced either using this work, or the analysis of Schoof and Hindmarsh (Reference Schoof and Hindmarsh2010), who looked at basal drag as a perturbation to extensional flow.
In taking a simple power-law rheology for the ice, we have ignored the complex dependence of ice viscosity on a host of factors, but most pertinently temperature, water content, damage and fabric (grain size and orientation). All of these factors are particularly likely to differ between the bulk of the ice stream and the regions of highest shear, namely in the shear margins and close to the bed (Harrison and others, Reference Harrison, Echelmeyer and Larsen1998; Minchew and others, Reference Minchew, Meyer, Robel, Gudmundsson and Simons2018). By considering these as additional mechanisms by which the ice viscosity depends on shear rate, one could attempt to parameterise their effects by using an altered power n in the ice rheology. As such, we have left all the expressions in our results as holding for general n. However, the flow of ice advects heat, damage and water downstream, so the dependence of viscosity on shear rate can be non-local in a way which cannot be fully captured by this steady, uniform model.
We note that all these assumptions in our model setup are common when working with simplified models of ice stream flow. However, it is worth considering, and exploring with more detailed numerical models, whether including basal shear-softening and more accurately capturing flow in this geometry is significant compared to the error introduced by simplifying the ice stream geometry and rheology in the first place.
6.2. Approximations in the analysis
Our mathematical analysis is not of the full-Stokes Eqn. (1) directly, but of an approximation using a linear shapefactor for the vertical shear stress. This approximation can be motivated by the separation of scales between the horizontal and vertical velocity gradients for wide streams, and the numerical stress fields (Fig. 2) provide reassurance that the approximation is reasonable. Further, the L1L2 model on which it is based is well-tested compared to full-Stokes simulations (Hindmarsh, Reference Hindmarsh2004; Cornford and others, Reference Cornford2013; Robinson and others, Reference Robinson, Goldberg and Lipscomb2022). However, this approximation breaks down for narrow streams over strong beds, as can be seen in Figure 5.
Indeed, a question remains as to what happens close to the sidewalls when the bed is not sliding. As discussed in the results, close to the sidewalls the surface velocity and basal stress go to zero across a region of width ~H, which leads to errors in the estimated flux through the ice stream. Numerical analysis of this problem dates back to Nye (Reference Nye1965) and Kozicki and others (Reference Kozicki, Chou and Tiu1966), but to our knowledge no analytic results exist for the flow of shear-thinning fluids near a corner. The problem is inherently two-dimensional, with significant vertical gradients in horizontal stress that the L1L2 approach is not designed to handle. Further study of these dynamics are certainly interesting from a fluid-dynamical perspective, but beyond the scope of the present paper.
We have calculated the centreline surface velocity in a manner that avoids considering the velocity field above the unyielded region, but our estimate of y u becomes less accurate as the region widens. Since our estimate of y u does tend to 0 as μN → τ d, we have the correct leading order expressions for u mid and Q in this regime. However, the estimate of sliding speed depends directly on y u and could therefore be improved by more detailed calculation of the dynamics above the unyielded bed (c.f. Haseloff and others, Reference Haseloff, Hewitt and Katz2019).
6.3. Applicability to other sliding laws
While previous work employing simple expressions for sliding speed and flux through ice streams has mainly centred on flow over a bed of uniform strength, larger scale ice-sheet models more frequently use regularised sliding laws. Here we consider to what extent our results for a constant-strength bed are applicable to sliding laws of the form u b = f(τ b).
If we view our results as a calculation for the basal sliding speed u b(τ b), flux Q(τ b), and surface velocity u s(τ b) for a given basal traction, we can look for self-consistent solutions to
effectively approximating the basal traction as its value in the centre of the bed. This allows us to predict the surface velocity and flux as functions of the sliding law and driving stress, which we can again compare to numerically calculated values.
Based on our expression for maximum sliding speed (44), (54) becomes
an implicit equation for u b or τ b which we cannot directly solve in general, but remains a computationally cheap problem compared to numerical solution of the Stokes equations.
Considering a power-law sliding law of the form
we see generally reasonable agreement between the numerical and self-consistently calculated values of Q and u mid (Fig. 9). However, the clearest difference between the numerical velocity fields for the plastic bed and power-law sliding is in the surface velocity profiles for very strong beds (large C). This can be attributed to the difference in the shape of the basal velocity profiles as u b → 0. Over a plastic bed, the margins stop yielding while the centre of the bed continues to slide, the ice above which is resisted primarily by lateral shear – this leads to fairly rounded surface velocity profiles. Over a power-law bed, once τ b ~ τ d, one can rearrange (56) to find a constant sliding speed u b = (τ d/C)m over the majority of the bed, and without lateral variations in sliding speed, the vertical shear stress dominates the force balance, leading to a correspondingly uniform surface velocity,
with shear margins of order H close to the sidewalls.
We can also look at the impact of mixed sliding laws of the form
which would correspond to a visco-plastic material at the bed (Iverson, Reference Iverson2010; Warburton and others, Reference Warburton, Hewitt and Neufeld2023). As might be expected, this gives qualitatively similar results to the plastic bed when C is small, and similar results to the power-law bed when C is large. In both regimes the self-consistent method for calculating τ b leads to reasonable agreement with numerically calculated values of Q and u mid (Fig. 10).
The comparisons suggest that our simple expressions for flux and velocity can produce useful results for more general sliding laws, by seeking self-consistent values of basal drag and basal sliding. It again suggests that surface velocity profiles are a key indicator of basal boundary conditions.
Alternatively, we can view the expression for sliding speed as a function of basal shear stress (55) as a mixed boundary condition that an ice stream exerts on its bed. It could be used as a more ice-stream-representative boundary condition in studies of till dynamics, beyond forcing with a constant shear stress or speed.
7. Conclusions
By considering the impact of basal drag on depth-integrated viscosity, we have produced improved expressions for the velocity field in an ice stream flowing over a bed of constant yield strength which capture the transition from the sidewall-resisted regime to the basal shear dominated regime. These more accurate expressions for the surface velocity (43), basal sliding speed (44) and total flux (51) through the stream are simple enough to include in wider models of the interaction of ice streams with their environment, and could also be used to rapidly invert for bed properties given surface observations. While we have focused on a simple geometry and simple ice rheology, these results still provide an improved estimate of ice behaviour under these oft-assumed conditions.
Acknowledgements
K. L. P. W. is supported by a postdoctoral fellowship from the Dartmouth Society of Fellows. Some of this work formed part of K. L. P. W.'s doctorate, supported by the Natural Environment Research Council (grant no. NE/L002507/1), on which Grae Worster and Andrew Hogg provided helpful comments. We also thank Logan Mann and Brent Minchew for useful conversations.
Open access
