Numerical model (2-D)

We consider a 2-D flowline numerical model. We use the open-source finite element full-Stokes solver Elmer/Ice (Gagliardini and others, Reference Gagliardini2013), with which we conducted a suite of numerical experiments. In particular, we used Elmer/Ice to solve the incompressible Stokes system:

(1)\begin{equation} \nabla \cdot \boldsymbol{\sigma} + \rho_\mathrm{i} \mathbf{g} = \mathbf{0}, \qquad \nabla \cdot \mathbf{u} = 0\,, \end{equation}

where $\boldsymbol{\sigma}, \rho_\mathrm{i}, \mathbf{g}\ \mathrm{and \,} \mathbf{u}$ are the Cauchy stress tensor, ice density, gravity and velocity field, respectively, and $\nabla\cdot$ is the divergence operator. The domain is of length $L_x$, and the system is subjected to periodic boundary conditions on both the horizontal ( $x$) and vertical ( $z$) components of the velocity field, $u(x,z)$ and $w(x,z)$, respectively. The simulation mesh is structured with spatially varying horizontal resolution across five zones: fine resolution (10 m) across the center $50$ km of the domain, intermediate resolution (100 m) in 50 km zones on both sides adjacent to the central region and coarse resolution (500 m) in the far-field beyond. The mesh is extruded vertically to the ice thickness $h$ with 40 vertical levels.

We ran a suite of experiments designed to investigate the relationship between the geometrical parameters of the flow and the resulting stress response. We later compare these empirical relationships to those that arise from the analytical model. To this end, we ran simulations over a range of parameter values for the ice thickness $h$, surface slope $\alpha$, patch length $\ell$ and basal sliding states $\gamma, C_\mathrm{b}$ for a total of 192 simulations, 96 for each rheology (Table 1). Basal sliding is characterized by a Weertman-type (e.g., Weertman, Reference Weertman1957) sliding law, which relates the basal shear stress $\tau_\mathrm{b}$ to velocity,

(2)\begin{equation} \tau_\mathrm{b} = C_\mathrm{b}|u_\mathrm{b}|^{1/n-1}u_\mathrm{b}\,, \end{equation}

where the subscript $\mathrm{b}$ indicates quantities along the basal interface and we follow a common simplification of the Weertman sliding law, taking the exponent to be $1/n$ (e.g., Weertman, Reference Weertman1974; Schoof, Reference Schoof2007b; Tsai and others, Reference Tsai, Stewart and Thompson2015; Sergienko, Reference Sergienko2022), where $n$ is the flow law exponent (Glen, Reference Glen1955). Here, we define $\gamma$ as the ratio between the contributions of deformation and sliding to velocity, which facilitates later comparison to the analytical model strictly valid when $\gamma \ll 1$,

(3)\begin{equation} \gamma \equiv \frac{[u_{\text{deform}}]}{[u_{\text{slide}}]} = \frac{B^{-n}\tau_\mathrm{b}^n h}{(\tau_\mathrm{b}/C_\mathrm{b})^n} = \frac{hC_\mathrm{b}^n}{B^n}\,, \end{equation}

where the brackets $[]$ indicate characteristic scale of the variable inside the brackets. All experiments were conducted with an ice sheet of constant thickness and surface slope (simply by tilting the gravity vector acting upon the flow according to $\alpha$, see, e.g., Fig. 1) and were subjected to periodic boundary conditions at the inflow (left) and outflow (right) of the ice body. We repeated this experiment for two ice rheologies: Newtonian ( $n=1$) and Glen ( $n=3$). The flow, in both cases, results in a strong perturbation in the velocity and stress fields at the center of the domain where the slippery patch resides, but returns to the background flow state, consistent with sliding-accommodated shallow-ice approximation (SIA), toward the edges of the domain (Fig. 2).

Table 1.

Range of parameter values used in the numerical experiments.

Table 1 Long description

Parameters are listed with their symbols, value sets, and units. Ice thickness, symbol h, has values 0.3, 1.0, 2.0, and 3.0 kilometers. Surface or bed slope, symbol alpha, has values 0.1, 0.5, and 1.0 degrees. Patch length, symbol ell, has values 2.5, 5.0, 10.0, and 25.0 kilometers. Flow law exponent, symbol n, is 1 or 3 with no units. Flow law parameter, symbol B to the minus n, is 5.0 times 10 to the minus 15 for n equals 1 and 2.4 times 10 to the minus 24 for n equals 3, with units pascals to the minus n seconds to the minus 1. Sliding ratio, symbol gamma, has values 0.1 and 0.5 with no units. Basal sliding coefficient, symbol C sub b, is defined in terms of B, gamma, h, and n, with units pascal seconds to the one over n meters to the minus one over n.

Figure 2.

Example simulation domain for $h=1000$ m, $\alpha=0.5^\circ$, $\ell=5$ km, $\gamma=0.1$. (a, d) Horizontal velocity field $u(x,z)$ response to the slippery patch for Newtonian $n=1$ and Glen $n=3$ rheologies, respectively. (b, e) Deviatoric normal stress $\tau_{xx}(x,z)$ response to the slippery patch for Newtonian and Glen rheologies, respectively. (c, f) Vertical profiles of horizontal velocity at the (periodic) inflow/outflow (dashed) boundaries compared to analytical basal sliding-accommodated SIA solution (red).

Figure 2 Long description

Example Domain (h equals 1000 m, alpha equals 0.5 degree, l equals 5 km, gamma equals 0.1). Three plots are labeled n equals 1 and three plots are labeled n equals 3. For n equals 1, horizontal velocity is shown as a two dimensional field plot. The horizontal axis label is x over l over 2 in kilometer, with ticks at minus 200, minus 100, 0, 100, 200. The vertical axis label is Elevation in meter, with ticks at 0, 250, 500, 750, 1000. The velocity scale is labeled u in m yr superscript minus 1, with labeled values 0, 25, 50, 75, 100. For n equals 1, normal stress is shown as a two dimensional field plot. The horizontal axis label is x over l over 2 in kilometer, with ticks at minus 200, minus 100, 0, 100, 200. The vertical axis label is Elevation in meter, with ticks at 0, 250, 500, 750, 1000. The stress scale is labeled sigma xx in kPa, with labeled values minus 100, 0, 100, 200. For n equals 1, a line plot is labeled Velocity at plus 200 m. The horizontal axis label is u in m yr superscript minus 1, with ticks at 125, 130, 135. The vertical axis label is x in kilometer, with ticks at minus 200, minus 100, 0, 100, 200. The legend lists SIA, Elmer left and Elmer right. The three curves lie close together and form an S shaped profile, with u near 125 at x near minus 200 and u near 135 at x near 200. For n equals 3, horizontal velocity is shown as a two dimensional field plot. The horizontal axis label is x over l over 2 in kilometer, with ticks at minus 200, minus 100, 0, 100, 200. The vertical axis label is Elevation in meter, with ticks at 0, 250, 500, 750, 1000. The velocity scale is labeled u in m yr superscript minus 1, with labeled values 400, 500, 600, 700, 800. For n equals 3, normal stress is shown as a two dimensional field plot. The horizontal axis label is x over l over 2 in kilometer, with ticks at minus 200, minus 100, 0, 100, 200. The vertical axis label is Elevation in meter, with ticks at 0, 250, 500, 750, 1000. The stress scale is labeled sigma xx in kPa, with labeled values minus 100, 0, 100, 200. For n equals 3, a line plot is labeled Velocity at plus 200 m. The horizontal axis label is u in m yr superscript minus 1, with ticks at 370, 380. The vertical axis label is x in kilometer, with ticks at minus 200, minus 100, 0, 100, 200. The legend lists SIA, Elmer left and Elmer right. The three curves lie close together and form an S shaped profile, with u near 370 at x near minus 200 and u near 380 at x near 200.

Analytical model (1-D)

To derive the one-dimensional (1-D) model, we begin this section by considering the depth-averaged shallow-shelf approximation (SSA) (MacAyeal, Reference MacAyeal1989)

(4)\begin{equation} \partial_x(2h\tau_{xx}) - \tau_b - \rho g h \partial_x s = 0\,, \end{equation}

where $g$ is gravity, $h$, $\rho$ and $s$ are ice thickness, density and surface, respectively, and $\tau_{xx}$ is longitudinal deviatoric stress. This approximation of Stokes requires a dominant contribution of sliding to the overall velocity, namely $\gamma\ll 1$ (Eqn (3)). As mentioned above (and in Appendix A), this depth-averaged model uses the following parameterization for the basal shear term,

(5)\begin{equation} \tau_\mathrm{b} = C|u|^{1/n-1}u\,, \end{equation}

where $C, u$ are depth-averaged quantities with $C=C_\mathrm{b}(u_\mathrm{b}/u)^{1/n}$ (Eqn (B.3)) relating the depth integrated quantity $C$ in our 1-D model with the basal quantities $C_\mathrm{b}, u_\mathrm{b}$ from the Elmer/Ice simulations. The flow law exponent $n$ comes from Glen’s flow law (Glen, Reference Glen1955), which relates the strain rate to the applied stress as

(6)\begin{equation} \tau_{xx} = B|\dot{\epsilon}_{xx}|^{1/n - 1}\dot{\epsilon}_{xx}\,, \end{equation}

with $B$ the ice hardness parameter and $\dot{\epsilon}_{xx}=\partial_xu$ the deviatoric strain rate, where $u$ is horizontal velocity. We take the simplifying case of constant ice thickness and surface slope, $s_x=-\alpha$. With these relationships, we rewrite the SSA equation in terms of velocity as

(7)\begin{equation} 2Bh(|u_x|^{1/n-1}u_x)_x - C|u|^{1/n-1}u + \rho g h \alpha = 0\,, \end{equation}

where derivatives are here and henceforth denoted with subscripts. The basal traction condition is piecewise constant across the domain,

(8)\begin{equation} C(x)=\begin{cases} 0, & x\in[0,\ell], \\[6pt] C \gt 0, & x\in(-\infty,0]\cup[\ell,\infty) \,. \end{cases} \end{equation}

We are predominantly interested in the spatial structure of the perturbed stress that arises from the presence of a spatially finite patch of free-slip. Our model considers the stress profile across the patch as setting the boundary condition for the surrounding BLs upstream and downstream from the patch. Inside the patch $C=0$ and the SSA reduces to

(9)\begin{equation} (2h\tau_{xx})_x = -\tau_\mathrm{d}, \qquad x \in [0, \ell] \end{equation}

where the $\tau_\mathrm{d}=\rho g h\alpha$ is the driving stress. Independent of $n$, this yields a linear membrane stress profile across the patch (see also Fig. 4b and f, middle panel) of

(10)\begin{equation} \tau_{xx}(x) = \tau_{xx}(0) - \frac{\tau_\mathrm{d}}{2h}x\,. \end{equation}

The antisymmetry of the membrane stress about the patch center, $\tau_{xx}(\ell/2) = 0$, gives the boundary condition at the upstream patch edge,

(11)\begin{equation} \tau_{xx}(0) = -\tau_{xx}(\ell) = \frac{\tau_\mathrm{d}\ell}{4h} =\frac{\rho g \alpha \ell}{4}\,, \end{equation}

such that Eqn (10) can be rewritten as

(12)\begin{equation} \tau_{xx} = \frac{\rho g \alpha \ell}{4}\left(1-2\frac{x}{\ell}\right)\,, \end{equation}

yielding a characteristic scale for the stress perturbation as

(13)\begin{equation} [\tau_{xx}] \equiv \frac{\rho g \alpha \ell}{4}\,.\\ \end{equation}

Outside the patch where our basal traction is given by $C(x)=C$, the governing equation is

(14)\begin{equation} (2h\tau_{xx})_x - C|u|^{1/n-1}u = -\tau_d\,, \qquad x\in (-\infty, 0] \cup [\ell, \infty) \end{equation}

subjected to the boundary condition above (Eqn (11)). Using the constitutive relation $\tau_{xx} = B|u_x|^{1/n-1}u_x$, the relevant length and velocity scales of this region can be found from the balance between deviatoric stress gradient term, $2Bh(|u_x|^{1/n-1}u_x)_x$, and basal drag, $C|u|^{1/n-1}u$, term in Eqn (7),

(15)\begin{equation} \frac{2Bh[u]^{1/n}}{[x]^{1/n}\cdot[x]} = C[u]^{1/n}\,, \end{equation}

from which $[u]^{1/n}$ cancels, giving the characteristic lengthscale of this region as

(16)\begin{equation} [x] = \left( \frac{2Bh}{C} \right)^{\frac{n}{n+1}}\,. \end{equation}

This region is host to a BL at the edge (upstream and downstream) of the patch, where the flow transitions from the far field, which is characterized by a balance between basal drag and driving stress to the patch interior, which is characterized by a balance between longitudinal stresses and driving stress. Thus, across this BL, there is a balance between longitudinal stresses and basal drag. The velocity scale is then determined by the membrane stress boundary condition at $x=0$ (Eqn (11)),

(17)\begin{equation} 2Bh\frac{[u]^{1/n}}{[x]^{1/n}} = \frac{\tau_\mathrm{d}\ell}{2}\,, \end{equation}

which gives

(18)\begin{equation} [u] = [x]\left( \frac{\tau_\mathrm{d} \ell}{4Bh} \right)^{n}= [x]\left( \frac{\rho g \alpha \ell}{4B} \right)^{n}\,. \end{equation}

Nondimensionalizing with

(19)\begin{equation} u = [u]\tilde{u} \end{equation}

(20)\begin{equation} x = [x]\tilde{x} \end{equation}

(21)\begin{equation} \tau_{xx} = [\tau_{xx}]\tilde{\tau}_{xx} \end{equation}

with $[u]$, $[x]$ and $[\tau_{xx}]$ defined in Eqns (18), (16) and (13) respectively, gives the dimensionless governing equation outside the patch,

(22)\begin{equation} (|\tilde{u}_{\tilde{x}}|^{1/n-1}\tilde{u}_{\tilde{x}})_{\tilde{x}} - |\tilde{u}|^{1/n-1}\tilde{u} = -2\delta\,, \end{equation}

where

(23)\begin{equation} \delta\equiv\frac{[x]}{\ell} = \frac{(2Bh/C)^{n/(n+1)}}{\ell} \end{equation}

is the ratio of the BL extent to the patch length (not to be confused with the Dirac delta function). Here, we decompose the velocity into a sum of two components, namely the background far-field velocity $u_\infty$ and the perturbed BL velocity $u'$. We do this in order to facilitate later comparison between this analytical solution and the depth-averaged numerical simulations with the background velocity subtracted away. In particular, we take

(24)\begin{equation} \tilde{u}(\tilde{x}) = \varepsilon + \tilde{u}'(\tilde{x})\,, \end{equation}

where

(25)\begin{equation} \varepsilon \equiv \frac{u_\infty}{[u]}\,, \end{equation}

relates the relative scales of the far-field and perturbation velocity components. Substitution into the dimensionless governing equation, noting that $\tilde{u}_{\tilde{x}} = \tilde{u}'_{\tilde{x}}$, gives

(26)\begin{equation} (|\tilde{u}'_{\tilde{x}}|^{1/n-1}\tilde{u}'_{\tilde{x}})_{\tilde{x}} - |\varepsilon + \tilde{u}'|^{1/n-1}(\varepsilon + \tilde{u}') = -2\delta\,. \end{equation}

We note that the dimensionless background balance, obtained by evaluating the governing equation far from the patch where $\tilde{u}' \to 0$, requires

(27)\begin{equation} |\varepsilon|^{1/n-1}\varepsilon = 2\delta\,, \end{equation}

which allows us to rewrite the governing equation as

(28)\begin{equation} (|\tilde{u}'_{\tilde{x}}|^{1/n-1}\tilde{u}'_{\tilde{x}})_{\tilde{x}} = |\varepsilon + \tilde{u}'|^{1/n-1}(\varepsilon + \tilde{u}') - |\varepsilon|^{1/n-1}\varepsilon\,. \end{equation}

For $n=1$, the right-hand side simplifies exactly: $(\varepsilon + \tilde{u}') - \varepsilon = \tilde{u}'$, independent of $\varepsilon$. For $n \gt 1$, since $\varepsilon = (2\delta)^n \ll 1$ when $\delta \ll 1$ (Eqn (27)), the right-hand side of Eqn (28) reduces to $|\tilde{u}'|^{1/n-1}\tilde{u}'$, giving the leading-order BL equation (dropping tildes henceforth) to be

(29)\begin{equation} (|u'_x|^{1/n-1}u'_x)_x = |u'|^{1/n-1}u'\,, \end{equation}

which is subjected to the following boundary conditions:

(30)\begin{equation} u' = 0 \quad (x\to-\infty) \end{equation}

(31)\begin{equation} |u'_x|^{1/n-1}u'_x = 1 \quad (x=0)\,, \end{equation}

where Eqn (31) is the dimensionless version of Eqn (11).

Taking typical scales of $h\sim10^3$ m, $\ell\sim10^4$ m, $B\sim7.5\times10^7$ Pa s $^{1/3}$, $C\sim5\times10^6$ Pa s $^{1/3}$ m $^{-{1/3}}$ gives $[x]\sim2.3\times10^3$ m (Eqn (16)). In order to strictly use SSA to describe this system, we must have $h \ll [x] \ll \ell$. With the scales estimated from above, we have $h\sim1\times10^3$ m, $[x]\sim2.3\times10^3$ m and $\ell\sim10^4$ m, meaning $\delta\sim0.23$. While the parameter combinations may not place us deeply within the asymptotic regime of $\delta\ll 1$, the resulting reduced-order model is reasonable, illustrative and interpretable, but not asymptotically strict. Thus, the remainder of this derivation shall be concerned with evaluating analytically the system under these assumptions, but note that the BL scaling (Eqn (16)) sufficiently captures relationships outside this parameter regime (see, e.g., Fig. 4c, d, g and h). For convenience, we define

(32)\begin{equation} y\equiv|u'_x|^{1/n-1}u'_x \end{equation}

(33)\begin{equation} \Rightarrow u'_x = |y|^{n-1}y\,, \end{equation}

which allows us to write our leading-order governing equation as

(34)\begin{equation} u'_x \frac{d}{du'}(|u'_x|^{1/n-1}u'_x) = |y|^{n-1}y\frac{dy}{du'} = |u'|^{1/n-1}u' \end{equation}

(35)\begin{equation} \Rightarrow \frac{1}{n+1}\frac{d}{du'}(|y|^{n+1}) = |u'|^{1/n-1}u' \,. \end{equation}

Integration of this equation (and applying the boundary condition from Eqn (30)) yields

(36)\begin{equation} |y|^{n+1} = n|u'|^{1/n+1} \end{equation}

(37)\begin{equation} \Rightarrow u'_x = |y|^{n-1}y = n^{\frac{n}{n+1}}|u'|^{1-1/n}u'^{1/n} \,, \end{equation}

where we can drop the absolute values given the solution (subjected to the boundary conditions) will be positive everywhere, and we write, for simplicity

(38)\begin{equation} u'_x = n^{\frac{n}{n+1}}u' \,. \end{equation}

Integration of the ordinary differential equation (ODE) gives

(39)\begin{equation} u'(x) = u'_0\exp{(n^{\frac{n}{n+1}}x)} \,, \end{equation}

where

(40)\begin{equation} u'_0\equiv u'(0)=\left(\frac{1}{n}\right)^{\frac{n}{n+1}} \end{equation}

Figure 3.

Surface longitudinal stress $\tau_{xx}^\mathrm{surf}$ from the 2-D Elmer/Ice numerical simulation ensemble for Newtonian ( $n=1$, left) and Glen ( $n=3$, right) rheologies. (a, d) Dimensional surface stress profiles across the slippery patch, colored by ice thickness $h$ with marker shape depicting surface slope $\alpha$ and linestyle depicting sliding ratio $\gamma$. Profiles with $\delta=[x]/\ell \lt 0.3$ are shown. (b, e) Maximum surface stress $\max|\tau_{xx}^{\,\mathrm{surf}}|$ versus the analytical prediction $\rho g \alpha \ell/4$ with proportionality constants $\kappa^s_{1,\tau}$ and $\kappa^s_{3,\tau}$ for Newtonian and Glen rheologies, respectively. (c, f) Decay lengthscale (DL; e-folding decay distance from the stress peak) versus the theoretical decay lengthscale $[x] = (2Bh/C)^{n/(n+1)}$ with $\kappa^s_{1,L}$ and $\kappa^s_{3,L}$ for Newtonian and Glen rheologies, respectively. Asymmetric error bounds (gray shade) denote the 16th–84th percentile of $\kappa$ prefactor estimates, i.e., the central 68% ( $\pm1\sigma$-equivalent range) of the distribution. We have restricted all plots, including the log–log plots to display simulation results satisfying $h \lt [x] \lt \ell$.

Figure 3 Long description

The image consists of six graphs comparing surface stress and decay lengthscale for Newtonian and Glen′s law rheologies. The top row shows surface stress profiles: (a) for Newtonian rheology with n equals 1 and (d) for Glen′s law with n equals 3. The x-axis is labeled x over l divided by 2 in kilometers and the y-axis is labeled tau subscript xx superscript surf in kilopascals, ranging from negative 1000 to 1000. The profiles are colored by ice thickness and show variations across the slippery patch. The bottom row includes scatter plots: (b) and (e) show stress versus analytical predictions, with the x-axis labeled rho g alpha l over 4 in kilopascals (analytical) and the y-axis labeled stress in kilopascals. (c) and (f) depict decay lengthscale versus theoretical decay lengthscale, with the x-axis labeled [2Bh over C] superscript n over n plus 1 in meters (analytical) and the y-axis labeled decay lengthscale in meters. Data points are marked with different shapes and colors, indicating various conditions and parameters. The plots illustrate the relationship between simulated and analytical values for both rheologies.

Figure 4.

Depth-averaged longitudinal stress $\bar{\tau}_{xx}$ from the 2-D Elmer/Ice simulation ensemble for (left) Newtonian $n = 1$ and (right) Glen $n = 3$ rheologies. Top row (a, e): Dimensional profiles across the patch, colored by ice thickness $h$ with marker shape denoting surface slope $\alpha$ and linestyle denoting $\gamma$. Middle row (b, f): Nondimensionalized stresses in the three regions of upstream boundary layer, patch interior and downstream boundary layer, each with analytical solution (solid black curve; Eqn (43)). Profiles with $\delta = [x]/\ell \lt 0.3$ are shown and collapse well onto the 1-D model. Bottom row (c, d, g, h): Log–log plots showing $\max|\bar{\tau}_{xx}| = \kappa \, \rho g \alpha \ell / 4$ and decay lengthscale $\mathrm{DL} = \kappa \, [x]$, with proportionality constants $\kappa^{da}_{1,\tau}$, $\kappa^{da}_{1,L}$ (Newtonian) and $\kappa^{da}_{3,\tau}$, $\kappa^{da}_{3,L}$ (Glen). Asymmetric error bounds (gray shade) denote the 16th–84th percentile of prefactor estimates. To comply with the modeling assumptions, we have restricted all plots, including the log–log plots, to display simulation results satisfying $h \lt [x] \lt \ell$.

Figure 4 Long description

Two line graphs titled, Depth Averaged Stress. One graph is labeled, Newtonian open parenthesis n equals 1 close parenthesis. One graph is labeled, Glen′s law open parenthesis n equals 3 close parenthesis. In both graphs, the horizontal axis is labeled, x minus ell over 2 open parenthesis kilometer close parenthesis, ranging from minus 40 to 40 with tick labels at minus 40, minus 30, minus 20, minus 10, 0, 10, 20, 30, 40. The vertical axis is labeled, tau xx bar open parenthesis kilopascal close parenthesis, ranging from minus 1000 to 1000 with tick labels at minus 1000, minus 500, 0, 500, 1000. Each graph contains multiple curves with circular markers and a horizontal reference line at 0. In each graph, curves are near 0 around minus 40, rise to a positive peak near x minus ell over 2 around minus 10, cross the 0 line near x minus ell over 2 around 0, drop to a negative trough near x minus ell over 2 around plus 10, then approach 0 again by plus 40. The image B showing three adjacent line plots labeled, Upstream BL, Patch interior and Downstream BL. The vertical axis is labeled, tau xx bar over open parenthesis rho g h bar alpha bar ell close parenthesis, ranging from minus 1.0 to 1.0 with tick labels at minus 1.0, minus 0.5, 0, 0.5, 1.0. The Upstream BL horizontal axis is labeled, x over ell, ranging from minus 6 to minus 2 with tick labels at minus 6, minus 4, minus 2. The Patch interior horizontal axis is labeled, x over ell, ranging from minus 2 to 2 with tick labels at minus 2, 0, 2. The Downstream BL horizontal axis is labeled, open parenthesis 2 minus x over ell close parenthesis, ranging from 0 to 6 with tick labels at 0, 2, 4, 6. Multiple curves are plotted together with a black curve. In Upstream BL, curves increase toward about 1.0 near x over ell equals minus 2. In Patch interior, curves form an approximately straight descending line from about 1.0 at x over ell equals minus 2 to about minus 1.0 at x over ell equals 2. In Downstream BL, curves increase from about minus 1.0 near open parenthesis 2 minus x over ell close parenthesis equals 0 toward 0 as the horizontal value approaches 6. The image C showing a log log scatter plot labeled, Stress. The horizontal axis is labeled, rho g h bar alpha bar ell open parenthesis kilopascal close parenthesis, with tick labels at 10 to the power 1, 10 to the power 2, 10 to the power 3. The vertical axis is labeled, max open brace tau xx bar close brace open parenthesis numerical close parenthesis, with tick labels at 10 to the power 1, 10 to the power 2, 10 to the power 3. Multiple point sets lie close to a diagonal black line across the plotted range. The image D showing a log log scatter plot labeled, Decay lengthscale. The horizontal axis is labeled, ell open parenthesis 2 B h bar over C close parenthesis to the power open parenthesis 1 over n plus 1 close parenthesis open parenthesis analytical close parenthesis, with tick labels at 10 to the power 0, 10 to the power 1, 10 to the power 2. The vertical axis is labeled, lambda open parenthesis kilometer close parenthesis, with tick labels at 10 to the power minus 1, 10 to the power 0, 10 to the power 1. Multiple point sets lie close to a diagonal black line. A shaded vertical band appears near the lowest horizontal values. The image E showing three adjacent line plots labeled, Upstream BL, Patch interior and Downstream BL, with the same axis labels and ranges as image B: vertical axis, tau xx bar over open parenthesis rho g h bar alpha bar ell close parenthesis, from minus 1.0 to 1.0; Upstream BL horizontal axis, x over ell, from minus 6 to minus 2; Patch interior horizontal axis, x over ell, from minus 2 to 2; Downstream BL horizontal axis, open parenthesis 2 minus x over ell close parenthesis, from 0 to 6. Curves follow the same overall shapes as in image B, including values near 1.0 at x over ell equals minus 2 in Upstream BL, an approximately straight decline from about 1.0 to about minus 1.0 across Patch interior and recovery from about minus 1.0 toward 0 in Downstream BL. The image F showing a log log scatter plot labeled, Stress, with the same axis labels and tick labels as image C: horizontal axis, rho g h bar alpha bar ell open parenthesis kilopascal close parenthesis, with 10 to the power 1, 10 to the power 2, 10 to the power 3; vertical axis, max open brace tau xx bar close brace open parenthesis numerical close parenthesis, with 10 to the power 1, 10 to the power 2, 10 to the power 3. Multiple point sets lie close to a diagonal black line. The image G showing a log log scatter plot labeled, Decay lengthscale, with the same axis labels and tick labels as image D: horizontal axis, ell open parenthesis 2 B h bar over C close parenthesis to the power open parenthesis 1 over n plus 1 close parenthesis open parenthesis analytical close parenthesis, with 10 to the power 0, 10 to the power 1, 10 to the power 2; vertical axis, lambda open parenthesis kilometer close parenthesis, with 10 to the power minus 1, 10 to the power 0, 10 to the power 1. Multiple point sets lie close to a diagonal black line. A shaded vertical band appears near the lowest horizontal values.

is determined by the boundary condition at $x=0$, namely $|y(0)|^{n+1} = 1^{n+1} = n|u'_0|^{\frac{n+1}{n}}$. This gives the velocity structure in the BL (as discussed further in Appendix B). For the deviatoric stress, we use the constitutive relation in dimensionless form (tildes dropped)

(41)\begin{equation} \tau_{xx} = |u'_x|^{1/n-1}u'_x = y(x)\,. \end{equation}

Substituting the velocity solution gives the dimensionless stress decay profile

(42)\begin{equation} \tau_{xx}(x) = n^{\frac{1}{(n+1)}}(u'(x))^{\frac{1}{n}} = \exp{\left( n^{-\frac{1}{(n+1)}}x \right)}\,, \end{equation}

where we drop the absolute values signs because $u' \gt 0$ and $u'_x \gt 0$ in the upstream BL. Downstream from the trailing edge of the patch ( $x=\ell$), the stress profile is solved for in the same way, with the stress decaying at the same rate into $x \gt \ell$.

Summary

In summary, the relevant spatial scale for the stress perturbation differs across the domain. Namely, $[x]=\left( \frac{2Bh}{C} \right)^{\frac{n}{n+1}}$ (Eqn (16)) dictates the stress decay lengthscale upstream and downstream from the patch (Eqn (42)), and $\ell$ characterizes the linear stress variation lengthscale across the patch interior (Eqn (12)). Returning to dimensional variables $\tau_{xx}$ and $x$, we write the summary for the deviatoric stress profile across the full domain:

(43)\begin{equation} \boxed{ \frac{\tau_{xx}(x)}{\rho g \alpha \ell/4} = \begin{cases} \exp\left(n^{-\frac{1}{n+1}}\frac{x}{[x]}\right), & x \leq 0 \\[10pt] 1 - 2\frac{x}{\ell}, & 0 \leq x \leq \ell \\[10pt] -\exp\left(-n^{-\frac{1}{n+1}}\frac{(x-\ell)}{[x]}\right), & x \geq \ell. \end{cases} } \end{equation}

Inside the BL regions, the spatial coordinate $x$ is measured from the nearest patch edge and scaled by $[x]=\left( \frac{2Bh}{C} \right)^{\frac{n}{n+1}}$ (Eqn (16)), while inside the patch the coordinate is measured from the upstream edge and scaled by the patch length $\ell$. This analytical solution is rendered as the solid black curve in Fig. 4b and f.

Surface stress

For the numerical experiments, we used the flow parameters that were presented in Table 1. The stress at the ice surface exhibits regions of tension and compression on the upstream and downstream sides of the slippery patch, respectively (Fig. 3a and b). This is expected, as the ice accelerates across the first half of the patch due to the basal traction loss and decelerates across the second half of the patch as it begins to feel the basal shear stress effects on the other side of the patch. Though our 1-D analytical model does not apply to the ice surface stresses, we observed from the simulations that the analytical scaling relationships between the ice geometry (e.g., $h, \alpha, \ell$) and the stress response hold. In particular, the magnitude of the surface stress scales linearly with $\rho g \alpha \ell/4$ and the decay length scales with $[x]$ (Eqn (16)) for experiments with Newtonian and Glen rheologies, albeit each with different prefactors. To obtain these prefactors, for each simulation, we compute the ratio of the simulated quantity to the theoretical scale (e.g., $\max|\tau_{xx}^{\mathrm{surf}}| / (\rho g \alpha \ell / 4)$). The associated prefactor $\kappa$ is then the median of those ratios across the ensemble of simulations, with gray bands spanning the 16th–84th percentiles, i.e., the central 68% ( $\pm1\sigma$-equivalent range) of the distribution. To quantify the decay lengthscales from the numerical experiments, we take the length upstream from the patch onset over which the perturbed stress is above $\max{|\tau_{xx}|}/e$ because the 1-D model implies a decaying exponential function, so this is the natural measure of the associated decay lengthscale. The stress magnitude relates to the theoretical stress scale with prefactors of $\kappa^s_{1,\tau}=0.78$ and $\kappa^s_{3,\tau}=0.95$ for the Newtonian and Glen experiments, respectively, while the decay lengths relate to the theoretical decay lengthscale with prefactors of $\kappa^s_{1,L}=1.40$ and $\kappa^s_{3,L}=3.22$ for the Newtonian and Glen experiments, respectively. These prefactors are all $\mathcal{O}(1)$, indicating that, to the first order, the characteristic scales $[\tau_{xx}]\equiv \rho g \alpha\ell/4$ and $[x] \equiv\left( \frac{2Bh}{C} \right)^{\frac{n}{n+1}}$ provide good estimates of surface stress magnitude and decay lengthscale, respectively. The superscript on the prefactors indicates whether it is fit for surface ( $s$) or depth-averaged ( $da$) quantities, while the subscript indicates the rheology ( $n=1,3$) and if the prefactor is for the stress magnitude ( $\tau$) or decay length scale ( $L$).

We note that we restrict the analysis to simulations satisfying the minimum requirement of $h \lt [x] \lt \ell$, which ensures the BL (decay) length scale is sufficiently broad such that we can reasonably use a depth-averaged model (i.e., SSA), but not broader than the slippery patch length $\ell$. For the profile panels (Fig. 3a and d), we further require that $\delta\equiv \frac{[x]}{\ell} \lt 0.3$ to display cases where BLs are well separated from the patch interior. The scaling panels (Fig. 3b, c, e and f) include all simulations, which satisfy the first filter ( $h \lt [x] \lt \ell$). These two restrictions are likewise applied to our analysis of the depth-averaged stresses, below.

Depth-averaged stress

For the depth-averaged stresses, we once again compare the scalings from the 1-D model to results from the simulations, which we now have depth-averaged. The depth-averaged stress magnitude scales linearly with $\rho g \alpha \ell/4$ and the decay length scales linearly with $[x]$ again for both rheologies with prefactors $\kappa^{da}_{1,\tau}, \kappa^{da}_{1,L}$ and $\kappa^{da}_{3,\tau}, \kappa^{da}_{3,L}$ for Newtonian and Glen rheologies, respectively (Figure 4c, d, g and h). The prefactor values for these depth-averaged results echo the patterns observed for the surface stresses, where the stress magnitudes are similar between rheologies ( $\kappa^{da}_{1,\tau}=0.93, \kappa^{da}_{3,\tau}=0.94$ for Newtonian and Glen, respectively). The decay length prefactors show more discrepancy ( $\kappa^{da}_{1,L}=1.04, \kappa^{da}_{3,L}=2.86$ for Newtonian and Glen, respectively).

The nondimensionalized stress profiles also agree well with the 1-D analytical model (Fig. 4b and f), confirming the existence of up (down)stream BLs with exponential stress decay (recovery) and linear stress decline across the patch. The collapsed numerical profiles are tighter to the 1-D analytical model for the Newtonian experiments than for the Glen experiments, consistent with and indeed expected from the model derivation (see, e.g., Eqns (25)–(29)). In particular, for $n=1$, the BL equation is exact and the background velocity terms cancel identically, independent of $\delta$. However, for $n \gt 1$, this cancellation is approximate, strictly requiring $\varepsilon = (2\delta)^n\ll1$, leading to increasing error with increasing $\delta$ (see e.g., Fig. C1). In Fig. 4a and e, we only included simulations with $\delta \lt 0.3$. The simulations with Glen rheology are thus subject to two compounding effects: the approximation is inexact for $\delta \gt 0$ and $\delta$ is systematically larger due to the broader extent of the BL for Glen rheologies, which degrades the agreement between the 1-D model and the depth-averaged simulations with Glen rheology, though the overall agreement remains good across both rheologies. A summary of the prefactors for both surface and depth-averaged stresses is provided in Table 2.

The stress magnitude prefactors are generally less than 1 while the prefactors on the decay lengths are generally greater than 1. This is to be expected. The 1-D model stress magnitude is an overestimate due to the fact that we do not fully resolve the dynamics near the transition point, which would require a full-Stokes treatment that smooths the stress between the upstream and patch-interior parts, leading to this overestimate of the 1-D model. The decay lengthscale, by definition, is the distance over which the perturbed stress remains above $1/e$ times the maximum stress. The $e$-folding lengthscale from the 1-D model (Eqn (42)) is $n^{1/(n+1)} [x]\geq[x]$.

Table 2.

Numerically determined prefactors $\kappa$ for scaling relationships of maximum stress perturbation $\max(\tau_{xx})$ and decay lengthscale DL.

Table 2 Long description

Four cases list numerically determined prefactors kappa for two quantities: maximum tau xx and decay lengthscale, each reported as a median with 16th to 84th percentile bounds. For n 1, surface, the maximum tau xx prefactor is 0.78 with plus 0.12 and minus 0.13, and the decay lengthscale prefactor is 1.40 with plus 0.50 and minus 0.31. For n 3, surface, the maximum tau xx prefactor is 0.95 with plus 0.07 and minus 0.03, and the decay lengthscale prefactor is 3.22 with plus 2.64 and minus 0.67. For n 1, depth averaged, the maximum tau xx prefactor is 0.93 with plus 0.04 and minus 0.07, and the decay lengthscale prefactor is 1.04 with plus 0.56 and minus 0.06. For n 3, depth averaged, the maximum tau xx prefactor is 0.94 with plus 0.03 and minus 0.06, and the decay lengthscale prefactor is 2.86 with plus 0.64 and minus 0.78.

Notes: DA and S denote depth-averaged and surface quantities, respectively. Values are medians across the simulation ensemble; superscripts and subscripts give the 16th–84th percentile bounds.

Implications

From the modeling work above, it is clear that the stress perturbation responding to a basal slippery patch is critically dependent on the ice thickness, surface slope, slippery patch length and surrounding basal sliding conditions. The guiding question of ‘what is the structure of stress perturbations from slippery patches?’ stems from the need to understand how the ice sheet might respond to continued Arctic warming and meltwater forcings. In particular, our work was initially motivated by the lake drainage cascade hypothesis, whereby the stress perturbation resulting from the drainage of one lake might trigger the drainage of another (e.g., Christoffersen and others, Reference Christoffersen, Bougamont, Hubbard, Doyle, Grigsby and Pettersson2018). The communication of stresses between lakes in response to their drainage is still yet to be fully understood. It remains unclear how and how far these stresses may be transmitted in response to HLD-induced meltwater injections (which we treat here as slippery patches). Modeling of membrane stress response to basal traction loss by Christoffersen and others (Reference Christoffersen, Bougamont, Hubbard, Doyle, Grigsby and Pettersson2018) suggests there may be a mechanism by which lakes located several tens of kilometers apart may transmit stresses and provide a sufficient triggering mechanism. However, in situ observations of HLD-induced stresses suggest coupling lengthscale on the order of only several (2–3) ice thicknesses (Stevens and others, Reference Stevens2024; Reference Stevens2026). In this study, we have aimed to contribute to this ongoing conversation in the literature, using our simple model to reveal the structure of an idealized slippery patch-induced stress perturbation.

As discussed above, we find that the magnitude of the stress perturbation is set by $\rho g \alpha \ell/4$, and decays into the upstream ice over a lengthscale of $(2Bh/C)^{n/(n+1)}$. Though the stress decay lengthscale increases with the local ice thickness, the exponential nature of the stress response suggests that we should be particular in our definition of the term ‘coupling length’. Say, for example, we are interested in the absolute distance over which the perturbed upstream stress is above a critical value $\tau_\mathrm{val}$, such that hydrofracture-triggered lake drainage can occur. This distance, referred to here as the absolute coupling length (ACL), is then a function not only of the stress decay lengthscale but also the maximum magnitude of the stress perturbation. The maximum stress magnitude and therefore the ACL are therefore unique functions of location on the ice sheet given a specific $(h, \alpha, \ell)$,

(44)\begin{equation} \mathrm{ACL} = \frac{[x]}{n^{-1/(n+1)}} \ln\!\left(\frac{\rho g \alpha \ell / 4}{\tau_\mathrm{val}}\right) = \frac{(2Bh/C)^{n/(n+1)}}{n^{-1/(n+1)}} \ln\!\left(\frac{\rho g \alpha \ell / 4}{\tau_\mathrm{val}}\right)\!. \end{equation}

As an illustrative example, we use $\tau_\mathrm{val}=30$ kPa due to the onset of fracture initiation in linear elastic fracture modeling of ice (e.g., Lai and others, Reference Lai2020). Note that more sophisticated mechanistic understanding of the critical stress required for hydrofracture-triggered lake drainage is needed to determine $\tau_\mathrm{val}$. ACL would decrease with $\tau_\mathrm{val}$ so if one were to use $\tau_\mathrm{val}=100$ kPa the resulting ACL would be smaller. Using Eqn (44), we generate phase spaces that demonstrate the sensitivity of ACL to different parameter combinations in our model (Fig. 5). In particular, we assume a fixed patch length of $\ell=5$ km and display ACL as a function of thickness $h$ and slope $\alpha$ (Fig. 5a and b) and then assume a fixed surface slope and display ACL as a function of thickness $h$ and patch length $\ell$ (Fig. 5c and d). For all phase spaces, we assume a sliding ratio of $\gamma=0.1$.

Regions of particular interest in the phase spaces are the intersection of regions bounded by: (1) typical driving stresses; (2) observed GrIS $h, \alpha$ values; and (3) model validity. We computed thickness and surface slope values consistent with typical GrIS driving stresses in the range $\tau_\mathrm{d}=\rho g h \alpha=50$–150 kPa (e.g., Cuffey and Paterson, Reference Cuffey and Paterson2010). For observed values, we derived the 90th-percentile $(h,\alpha)$ envelope of the GrIS from BedMachine v6 surface and bed topography (Morlighem and others, Reference Morlighem2017) and ITS_LIVE (Version 2) surface velocities (Gardner and others, Reference Gardner2025). In particular, we computed the joint distribution of thicknesses $h$ and along-flow surface slope $\alpha$ from BedMachine v6 and ITS_LIVE, respectively. The ice surface field was coarsened and smoothed to a resolution of 1.5 km and surface slopes were computed from this smoothed DEM. The slopes were then projected onto the local ice-flow direction $\hat{u}=(v_x,v_y)/|\mathbf{v}|$ which was obtained from ITS_LIVE, giving the along-flow slope $\alpha=\arctan(|\nabla s\cdot\hat{u}|)$. Locations with speeds $|\mathbf{v}| \lt 5$ m yr $^{-1}$ were excluded. The resulting $(h,\alpha)$ pairs over grounded ice were binned (120 $\times$ 120 bins spanning $h\in[200,3000]$ m and $\alpha\in[0.05^\circ, 1.0^\circ]$). The 90th-percentile envelope was then defined as the smallest set of histogram bins whose cumulative count comprises 90% of all grounded ice pixels. For model validity, we use the metric $\delta=[x]/\ell \lt 1$, as a reference for the location on the phase space where the BL length exceeds the patch length, above which the 1-D model deteriorates.

In the bottom row (Fig. 5c and d), the vertical dashed line marks $\ell=10$ km, an approximate upper bound on slippery patch lengths, informed by in situ observations of slippery patch extent (Stevens and others, Reference Stevens2015; Reference Stevens2024). The magenta outline traces the region where these observational and model-validity constraints are simultaneously satisfied. Within this region, ACLs are up to 14 km for $n=1$ and up to 12 km for $n=3$.

Figure 5.

Phase space maps of the absolute coupling length (ACL, Eqn (44)) for the analytical 1-D boundary layer model with $\gamma=0.1$ and stress threshold value $\tau_{val}=30$ kPa. (a, b) ACL as a function of ice thickness $h$ and surface slope $\alpha$ for a patch length of $\ell=5$ km for Newtonian ( $n=1$) and Glen ( $n=3$) rheologies, respectively. (c, d) ACL as a function of $h$ and patch length $\ell$ for a fixed surface slope $\alpha=0.5^\circ$. In all panels, the black dash-dot curve marks the boundary where $\delta=1$, outside of which the boundary layer lengthscale exceeds the patch length ( $\delta \gt 1$) and the 1-D model breaks down significantly for the Glen rheology (Eqns (22) and (23)). In (a, b), the black dashed curves bound typical GrIS driving stress range $\tau_d = \rho g h\alpha$ between $50$ and $150$ kPa; the solid black contour delineates the 90th percentile envelope of $(h,\alpha)$ values derived from BedMachine v6 (Morlighem and others, Reference Morlighem2017) and ITS_LIVE (Gardner and others, Reference Gardner2025). In (c, d), the vertical dashed line marks $\ell=10$ km, an approximate upper bound for slippery patches, and the solid horizontal lines bound the 5th–95th percentile thickness range of GrIS grounded ice pixels with along flow slope $\alpha=0.5^\circ\pm0.05^\circ$. The magenta outline traces the region where all observational and model validity constraints are simultaneously satisfied, and thus is the region of interest in the phase space.

Figure 5 Long description

Contour plot titled, Absolute Coupling Length, Stress Above 30 kPa, gamma equals 0.1, with text, Top fixed ell equals 5 km. The horizontal axis label reads, Surface Slope alpha degrees, ranging from 0.2 to 1.0. The vertical axis label reads, Thickness h m, ranging from 500 to 3000. A vertical scale bar shows values from 4 to 36. Contour labels visible include 4, 8, 12, 16, 20, 24, 28, 32, 36. The contour values increase toward higher thickness and higher surface slope, with lower labeled contours near thickness around 500 and surface slope around 0.2 and higher labeled contours near thickness around 3000 and surface slope around 1.0. A magenta outlined shape labeled, region of interest, spans surface slope roughly 0.2 to 0.9 and thickness roughly 500 to 1200. Within this outlined region, visible contour labels include 4, 8, 12 and 16. A legend lists line entries including, GrIS tau d equals 10 kPa, GrIS tau d equals 100 kPa, GrIS BedMachine and ITS LIVE and, region of interest. Contour plot titled, Absolute Coupling Length, Stress Above 30 kPa, gamma equals 0.1, with text, Bottom fixed alpha equals 0.5 degrees. The horizontal axis label reads, Surface Slope alpha degrees, ranging from 0.2 to 1.0. The vertical axis label reads, Thickness h m, ranging from 500 to 3000. A vertical scale bar shows values from 4 to 36. Contour labels visible include 4, 8, 12, 16, 20, 24, 28, 32, 36. The contour values increase toward higher thickness and higher surface slope, with lower labeled contours near thickness around 500 and surface slope around 0.2 and higher labeled contours near thickness around 3000 and surface slope around 1.0. A magenta outlined shape labeled, region of interest, spans surface slope roughly 0.2 to 0.9 and thickness roughly 500 to 1200. Within this outlined region, visible contour labels include 4, 8, 12 and 16. Contour plot titled, Absolute Coupling Length, Stress Above 30 kPa, gamma equals 0.1. The horizontal axis label reads, Patch length ell km, ranging from 2.5 to 20.0. The vertical axis label reads, Thickness h m, ranging from 500 to 3000. A vertical scale bar shows values from 4 to 36. Contour labels visible include 4, 8, 12, 16, 20, 24, 28, 32, 36. A vertical dashed line is drawn at patch length 5.0. The contour values increase toward higher thickness and toward larger patch length, with lower labeled contours near thickness around 500 and patch length around 2.5 and higher labeled contours near thickness around 3000 and patch length around 20.0. A magenta outlined shape labeled, region of interest, spans patch length roughly 2.5 to 10.0 and thickness roughly 500 to 1500. Within this outlined region, visible contour labels include 4, 8, 12 and 16. Contour plot titled, Absolute Coupling Length, Stress Above 30 kPa, gamma equals 0.1. The horizontal axis label reads, Patch length ell km, ranging from 2.5 to 20.0. The vertical axis label reads, Absolute coupling length L m, ranging from 500 to 3000. A vertical scale bar shows values from 4 to 36. Contour labels visible include 4, 8, 12, 16, 20, 24, 28, 32, 36. The contour values increase toward higher absolute coupling length and toward larger patch length, with lower labeled contours near absolute coupling length around 500 and patch length around 2.5 and higher labeled contours near absolute coupling length around 3000 and patch length around 20.0. A magenta outlined shape labeled, region of interest, spans patch length roughly 2.5 to 10.0 and absolute coupling length roughly 500 to 1500. Within this outlined region, visible contour labels include 4, 8, 12 and 16. Across images A and B, the plotted variables are surface slope alpha degrees on the horizontal axis and thickness h m on the vertical axis and both show increasing contour values toward higher thickness and higher surface slope. Across images C and D, the horizontal axis is patch length ell km and both show increasing contour values toward larger patch length and higher values on the vertical axis. The magenta region of interest appears in all four images and contains contour labels from 4 to 16.

To illustrate how the ACL dependencies play out, and indeed compete, in a spatial context, we evaluate the analytical model along an idealized convex ice-sheet profile (e.g., Weertman, Reference Weertman1976). We assume a perfectly plastic ice sheet with a divide-to-terminus distance of $L=300$ km resting on a flat bed with yield stress $\tau_y = 50$ kPa. The thickness and surface slope are thus related by $h(x) = \sqrt{2\tau_y x / (\rho g)}$ and $\alpha(x) = \tau_y / (\rho g h)$, where $x$ is the distance from the terminus (Fig. 6). This profile traces a 1-D path through the phase space (Fig. 5), with a steep thin margin and gentle thick interior.

Figure 6a and b show the ice-sheet cross section colored by the ACL and maximum stress magnitude, respectively, assuming a 5 km slippery patch with $\gamma=0.1$ and $n=3$. Tracing from the terminus, the ACL first increases toward the interior until it reaches a maximum value of $\sim 3$ km around 75 km from the terminus, at which point the ACL decreases monotonically toward the interior (Fig. 6a, c and d). Despite the fact that the decay lengthscale, $[x]$, increases with thickness $h$, we see here the competing effect of the decreasing stress magnitude. Stress magnitude is a strong function of surface slope $\alpha$ which decreases monotonically from terminus to interior (Fig. 6b, c and d). As a result, the ACL reaches a peak value somewhere between the interior and the terminus, corresponding to conditions typical of the GrIS ablation zone which is host to many supraglacial lakes.

Figure 6.

Absolute coupling length (ACL) and maximum longitudinal deviatoric stress $\tau_{xx}$ evaluated along an idealized convex ice-sheet profile (Weertman, Reference Weertman1976) with divide-to-terminus distance $L=300$ km, yield stress $\tau_y=50$ kPa, Glen ( $n=3$) rheology, sliding ratio $\gamma=0.1$, patch length $\ell=5$ km and stress threshold $\tau_{val}=30$ kPa (Eqn (44)). (a) Surface topography colored by ACL, showing a peak of $\sim 3$ km in the ablation zone $\sim 75$ km from the terminus before declining toward the interior. (b) Same topography colored by the maximum stress magnitude $\max{\tau_{xx}}=\rho g \alpha\ell/4$, which decreases monotonically from the margin toward the flattening interior. (c) ACL (blue, left axis) and maximum stress (black, right axis) as functions of ice thickness $h$, demonstrating competing effects of increasing decay lengthscale and decreasing stress magnitude with growing thickness. (d) ACL and stress maximum as functions of slope.

Figure 6 Long description

Absolute Coupling Length (tau subscript d equals 30 kPa) and Maximum Stress Across Convex Ice Sheet Profile (Weertman 1976) The image A showing a curve labeled, Topography with Absolute Coupling Length. The horizontal axis is labeled, Distance From Terminus [km], with tick labels 300, 250, 200, 150, 100, 50, 0. The vertical axis is labeled, Topography Profile [m], with tick labels 0, 500, 1000, 1500, 2000. A single smooth curve slopes downward from near 2000 at 300 to near 0 at 0. A vertical color scale at right is labeled, ACL [km], with tick labels 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0. The image B showing a curve labeled, Topography with Max Stress. The horizontal axis is labeled, Distance From Terminus [km], with tick labels 300, 250, 200, 150, 100, 50, 0. The vertical axis is labeled, Topography Profile [m], with tick labels 0, 500, 1000, 1500, 2000. A single smooth curve slopes downward from near 2000 at 300 to near 0 at 0. A vertical color scale at right is labeled, max(tau subscript x) [kPa], with tick labels 25, 50, 75, 100, 125, 150. The image C showing a two line graph labeled, ACL and max(tau subscript x) vs. h. The horizontal axis is labeled, h [m], with tick labels 0, 500, 1000, 1500, 2000. The left vertical axis is labeled, Absolute Coupling Length [km], with tick labels 1.0, 1.5, 2.0, 2.5, 3.0. The right vertical axis is labeled, Max Stress [kPa], with tick labels 0, 100, 200, 300, 400, 500. One line rises from about 1.0 at h equals 0 to a peak near 3.0 around h between 1000 and 1500, then declines to about 2.0 by h equals 2000. A second line increases gradually from near 0 at h equals 0 to about 100 at h equals 2000. The image D showing a two line graph labeled, ACL and max(tau subscript x) vs. alpha. The horizontal axis is labeled, alpha [degrees], with tick labels 1, 2, 3, 4. The left vertical axis is labeled, Absolute Coupling Length [km], with tick labels 1.0, 1.5, 2.0, 2.5, 3.0. The right vertical axis is labeled, Max Stress [kPa], with tick labels 0, 100, 200, 300, 400, 500. One line decreases from about 3.0 at alpha equals 1 to about 1.0 at alpha equals 4. A second line rises from about 1.0 at alpha equals 1 to near 3.0 around alpha near 3, then declines to about 2.0 at alpha equals 4.

Some observations suggest that a higher stress exponent, $n=4$, may better describe ice deformation in Greenland (Gillet-Chaulet and others, Reference Gillet-Chaulet, Hindmarsh, Corr, King and Jenkins2011; Bons and others, Reference Bons2018). Our framework permits a direct analytical estimate of its effect on the ACL. At fixed sliding ratio $\gamma$, the basal and depth-averaged friction coefficients are both proportional to the ice hardness, $C_b=B(\gamma/h)^{1/n}$ and $C=C_b(u_b/u)^{1/n}$ (Eqn (B.3)), so $B$ cancels and the decay length (Eqn (16)) reduces to $[x] = (2Bh/C)^{n/(n+1)} = h\left[2^n(n+2+2\gamma)/((n+2)\gamma)\right]^{1/(n+1)}$, depending only on $h$, $n$ and $\gamma$. The ACL (Eqn (44)) then depends on $n$ only through the slowly varying prefactor $\left[2^n\,n\,(n+2+2\gamma)/((n+2)\gamma)\right]^{1/(n+1)}$, which approaches $2$ as $n\to\infty$ so that the coupling length tends to $2h\ln(\rho g\alpha\ell/4\tau_\mathrm{val})$ (Fig. 7). At representative ablation-zone conditions ( $h=1$ km, $\alpha=0.5^\circ$, $\ell=5$ km, $\gamma=0.1$), this gives an ACL of ${\approx}\,4.3$ km for $n=4$, which is ${\approx}\,8\%$ shorter than the ${\approx}\,4.7$ km for $n=3$ (Fig. 7).

Figure 7.

Absolute coupling length (ACL) versus flow-law exponent $n$ for sliding ratios $\gamma = 0.1,\,0.3,\,0.5$, at fixed ice thickness $h=1$ km, surface slope $\alpha=0.5^\circ$, patch length $\ell=5$ km and stress threshold $\tau_\mathrm{val}=30$ kPa. At fixed $\gamma$, the ACL depends only on $(h,\alpha,\ell,\gamma,n)$. Filled circles mark the integer exponents, with the $n=1$, $n=3$ and $n=4$ values labeled on the $\gamma=0.1$ curve; the dashed line is the $n\to\infty$ limit $2h\ln(\rho g\alpha\ell/4\tau_\mathrm{val})$. The coupling length is only weakly sensitive to $n$, and non-monotonic at larger $\gamma$.

Figure 7 Long description

H equals 1 kilometer, alpha equals 0.5 degree, ell equals 5 kilometer, tau val equals 30 kilopascal. The x-axis is labeled, Flow law exponent n, ranging from 1.0 to 5.0 with tick marks every 0.5. The y-axis is labeled, Absolute coupling length, kilometer, ranging from 2.5 to 5.5 with tick marks every 0.5. A legend in the upper right lists three line series: gamma equals 0.1, gamma equals 0.3 and gamma equals 0.5. For gamma equals 0.1, the curve starts near 5.5 at n equals 1.0, passes through a marked point at n equals 2.0 with label n equals 2: 5.47 kilometer, passes through a marked point at n equals 3.0 with label n equals 3: 4.71 kilometer, passes through a marked point at n equals 4.0 with label n equals 4: 4.34 kilometer and ends near 4.1 at n equals 5.0. For gamma equals 0.3, the curve starts near 3.35 at n equals 1.0, has a marked point near 3.7 at n equals 2.0, a marked point near 3.65 at n equals 3.0, a marked point near 3.55 at n equals 4.0 and ends near 3.45 at n equals 5.0. For gamma equals 0.5, the curve starts near 2.75 at n equals 1.0, has a marked point near 3.2 at n equals 2.0, a marked point near 3.25 at n equals 3.0, a marked point near 3.2 at n equals 4.0 and ends near 3.15 at n equals 5.0. A dashed horizontal line near the bottom is labeled, n to infinity: 2h ln rho g alpha ell over 4 tau val equals 2.37 kilometer.

Our model assumes that a lake drainage forms a spatially finite region of free slip among an otherwise homogeneous basal traction condition. While this is a reasonable way to represent the total loss of basal traction over an HLD injection site due to an increase in local effective pressure sufficient to counter the full ice overburden (e.g., Lai and others, Reference Lai2021; Stevens and others, Reference Stevens2024), it ignores treatment of the subglacial hydrological system as either ‘channelized’ or ‘cavity’ based. The subglacial system has been modeled (Schoof, Reference Schoof2010; Hewitt, Reference Hewitt2013) and observed (Hoffman and others, Reference Hoffman, Catania, Neumann, Andrews and Rumrill2011; Chandler and others, Reference Chandler2013; Andrews and others, Reference Andrews2014; Reference Andrews2018) to evolve between these end member states over the course of a melt season. However the ACL metric may still be useful for estimating the slippery patch-driven stress communication lengthscales, at least qualitatively, in the context of the GrIS subglacial hydrological system. In the early-to-mid melt season, when many HLD events occur (Das and others, Reference Das2008; Doyle and others, Reference Doyle2013; Stevens and others, Reference Stevens2015; Chudley and others, Reference Chudley2019), the subglacial system is characterized by strong ice-bed coupling, parameterized in our model by the sliding law coefficient, $C$. A larger value of $C$ predicts a smaller ACL, implying that lake drainages earlier in the melt season may induce shorter-lengthscale stress perturbations than those in the later melt season, when the ice-bed frictional coupling is reduced.

Several important physical processes are not captured in our present model and represent natural extensions of this work. For one, ours is not a time-dependent model. The transient evolution of the stress perturbation following a lake drainage event and the timescale over which the perturbation arises and decays remain to be characterized. Likewise, the viscoelasticity of ice and the motion of the recently injected subglacial water both vertically and laterally are not considered. Related, we use a simplified representation of the basal traction and ignore variation in basal topography, which in our view warrants further analysis. The current work is limited also by being 1-D along-flow (2-D numerically, but still along-flow) and thus ignores lateral stresses. Addressing these extensions, some of which are the subject of current ongoing and future planned work, will be necessary for a more complete characterization of the impacts of lake drainage-induced injections of meltwater to the GrIS bed.