2.1. Formulation

In this study, we employ a model closely related to that in Hindmarsh (Reference Hindmarsh2004), to describe ice flow over an infinite sloped plane, except that we neglect basal sliding and assume no surface accumulation. The geometry is described by a Cartesian coordinate system, $(\hat x,\hat y, \hat z)$ , where the $\hat z$ -axis is perpendicular to the plane and the $\hat x$ -axis aligns with the flow direction. The bedrock surface is flat, so that the lower boundary of the ice is $\hat z=0$ , and $\hat h(\hat x,\hat y,\hat t)$ represents the ice thickness. The ice layer is shallow, with typical depths ${\mathcal{H}}$ much less than the characteristic length scale ${\mathcal{L}}$ for variations in the plane of the ice; the aspect ratio ${\varepsilon } = {{\mathcal{H}}} / {{\mathcal{L}}} \ll 1$ . We assume further that the bedrock slope is small, of the same order of magnitude as the aspect ratio, and therefore define a scaled slope $s$ such that $\varepsilon {s}$ is the actual slope; see figure 1. Note that we identify dimensional variables with hats, removing this decoration below when we switch to dimensionless variables.

This choice of scaling ensures that both the surface slope and the bedrock slope contribute to the leading-order momentum balance, as shown later, and therefore preserves the analogy with ice-sheet flow in an idealised flow-down-a-slope setting. We use this setting as a simplification of more realistic geophysically relevant flow geometries. In a more realistic model, one would generally include accumulation or ablation, allow the geometry and steady state to vary in the streamwise direction and determine the thickness together with the flow through mass conservation and appropriate upstream and downstream boundary conditions. Our aim here is not to model any particular glacier or ice sheet directly, but to isolate the thermomechanical effects in a simpler setting.

Figure 1.

Sketch of the model geometry. Cartesian coordinates are aligned with an inclined plane corresponding to the ice–rock interface. The shading represents temperature perturbations in both the ice and bedrock.

The ice is modelled as an inertialess, incompressible fluid, so that

(2.1) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \tau } - \boldsymbol{\nabla }\hat p +{\rho }\boldsymbol {g} &=\boldsymbol{0},\end{align}

(2.2) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat u}&=0 ,\end{align}

for $0\lt \hat z\lt \hat h$ , where $\hat p$ is the pressure, $\boldsymbol{\hat \tau }$ is the deviatoric stress tensor, $\rho$ is the ice density and $\boldsymbol{\hat u}=(\hat u, \hat v,\hat w)$ is the velocity field. The force of gravity vector is

(2.3) \begin{align} \boldsymbol {g} = \frac {{g}}{\sqrt {1 + {\varepsilon }^2{s}^2}} \begin{pmatrix} {\varepsilon }{s} \cr -1 \end{pmatrix}, \end{align}

where $g$ is gravitational acceleration.

The heat equation for the ice combines the effects of advection, diffusion and viscous dissipation, and can be written as

(2.4) \begin{align} \hat T_{\hat t}+\hat {\boldsymbol u}\boldsymbol{\cdot }\boldsymbol{\nabla }\hat T={\kappa _{_{\!I}}} {\nabla} ^2\hat T +\frac {1}{{\rho } {c_p}} \boldsymbol{\hat \tau }:\hat {\boldsymbol {D}} \quad (0\lt \hat z\lt \hat h). \end{align}

Here, $\hat T$ is the temperature, $\kappa _{_{\!I}}$ and $c_p$ are the thermal diffusivity and heat capacity of ice, respectively, and $\hat {\boldsymbol {D}}=(\boldsymbol{\nabla }\hat {\boldsymbol u} + ( \boldsymbol{\nabla }\hat {\boldsymbol u})^{\mathsf{T}})/2$ is the strain-rate tensor. We also use subscripts of $\{\hat t,\hat x, \hat y, \hat z\}$ to denote partial derivatives. In the bedrock, the absence of advection simplifies the heat equation to

(2.5) \begin{align} \hat T_{\hat t}={\kappa _{_{\!R}}} {\nabla} ^2\hat T \quad (\hat z\lt 0), \end{align}

where we use a subscript $R$ rather than $I$ to differentiate the material properties of the bedrock from those of the ice.

The ice surface is stress free and held at fixed temperature, assuming an instantaneous relaxation to a constant atmospheric temperature $T_a$ . In combination with the kinematic condition, we impose

(2.6) \begin{align} \left . \begin{aligned} \hat {\boldsymbol{\tau }}\boldsymbol{\cdot }\boldsymbol n - \hat p \; \boldsymbol n &=\boldsymbol{0} \\ \hat T&={T_a} \\ \hat h_{\hat t}+\hat u \hat h_{\hat x}+\hat v \hat h_{\hat y}&=\hat w \end{aligned} \right \} \quad \textrm {at} \ \hat z = \hat h, \end{align}

where $\boldsymbol{n}$ is the outward unit vector normal to the surface. At the ice–bedrock interface, we impose no slip and the continuity of the temperature and heat flux, stated as

(2.7) \begin{align} \left . \begin{aligned} \hat {\boldsymbol{u}}&= \boldsymbol{0} \\ \hat T^-&=\hat T^+ \\ {K_{_{\!R}}} \hat T_{\hat z}^-&= {K_{_{\!I}}} \hat T_{\hat z}^+ \end{aligned} \right \} \quad \textrm {on} \ \hat z = 0, \end{align}

where $K_{_{\!I}}$ and $K_{_{\!R}}$ denote the thermal conductivities of ice and rock, and the $\pm$ superscripts denote the limits approached from above (the ice layer) or below (the bedrock), respectively. The no-slip assumption is unrealistic, as basal sliding would occur when temperatures approach the melting point. However, to isolate the effect of temperature-dependent viscosity, we neglect basal sliding. Deep within the bedrock, we prescribe a geothermal heat flux $G$ , so that

(2.8) \begin{align} -{K_{_{\!R}}}\hat T_{\hat z}\rightarrow {G}, \quad \text{for} \quad \hat z\rightarrow -\infty . \end{align}

To complete the model, we require a constitutive law relating strain rate $\hat {\boldsymbol {D}}$ to stress $\boldsymbol{\hat \tau }$ . We adopt a generalised Newtonian fluid model, in which the effective viscosity is a function of both temperature and stress, so that

(2.9) \begin{align} \hat {\boldsymbol {D}}={\hat {f}}(\hat \tau ,\hat T)\frac {\boldsymbol{\hat \tau }}{\hat \tau }, \quad \hat \tau ^2=\frac {1}{2}\boldsymbol{\hat \tau }:\boldsymbol{\hat \tau }. \end{align}

2.2. Non-dimensionalisation

Following the usual shallow-ice approximation (Fowler & Larson Reference Fowler and Larson1978; Hutter Reference Hutter1983), we define scaled variables as

(2.10) \begin{align} (x,y) = \frac {(\hat x,\hat y)}{ {{\mathcal{L}}}}, \quad (z,h) = \frac {(\hat z, \hat h)}{{{\mathcal{H}}}}, \quad (u,v) = \frac {(\hat u, \hat v)}{ {{\mathcal{U}}}},\quad w = \frac {\hat w}{{\varepsilon }{{\mathcal{U}}}}, \quad t = \frac {{{\mathcal{U}}}\hat t}{{{\mathcal{L}}}}, \end{align}

(2.11) \begin{align} (p, {\varepsilon } \boldsymbol{\tau } ) = \frac {(\hat p, \boldsymbol{\hat \tau })}{{\rho } {g} {{\mathcal{H}}} } , \quad {f}( \tau , T) = \frac {{\hat {f}}(\hat \tau ,\hat T)}{{\hat {f}}({\varepsilon }{\rho }{g}{{\mathcal{H}}},{T_a})} , \quad T = \frac {\hat T - {T_a}}{{[\Delta T]}}, \end{align}

where

(2.12) \begin{align} {[\Delta T]} = \left (\! \left . \frac {1}{{\hat {f}}}\frac {\partial {\hat {f}}}{\partial \hat T} \right )^{-1} \right |_{\{\hat \tau ={\varepsilon }{\rho }{g}{{\mathcal{H}}},\hat T={T_a}\}}\!\!\!\!\!, \quad\! {{\mathcal{H}}} = \frac {{c_p}{[\Delta T]}}{{g}} , \quad\! {{\mathcal{U}}} = 2{{\mathcal{H}}} {\hat {f}}({\varepsilon }{\rho }{g}{{\mathcal{H}}},{T_a}), \quad\!\! {{\mathcal{L}}} = \frac {{{\mathcal{H}}}^2 {{\mathcal{U}}}}{{\kappa _{_{\!I}}}}. \end{align}

Through (2.12), we define scales for $\hat T-{T_a}$ and $\hat {f}$ from the shear rate and its sensitivity to temperature, evaluated at the characteristic shear stress ${\varepsilon }{\rho }{g}{{\mathcal{H}}}$ and atmospheric temperate $T_a$ . The remaining characteristic scales, ${\mathcal{H}}$ , ${\mathcal{L}}$ and ${\mathcal{U}}$ are chosen to ensure unit Brinkmann and Péclet numbers. Thus, $\mathcal H$ and $\mathcal L$ are determined by the material parameters and thermo-mechanical balances in the governing equations, rather than from geometric features of a particular ice sheet. Although this is not the standard scaling in glacier modelling, it is analytically convenient because it retains advection, diffusion and viscous heating at leading order while reducing the number of control parameters.

To leading order in $\varepsilon$ , the pressure is hydrostatic, $p=h-z$ . Substituting into the momentum equations in the $x$ and $y$ -directions, integrating with respect to $z$ and applying the stress-free boundary condition at the surface, we obtain

(2.13) \begin{align} \boldsymbol{\tau } = \begin{pmatrix}\tau _{13} \\ \tau _{23} \end{pmatrix} &= (h-z) \begin{pmatrix} {s} - h_x \\ -h_y \end{pmatrix}. \end{align}

The remaining leading-order system is given by

(2.14) \begin{align} u_x + v_y + w_z &= 0,\\[-6pt] \nonumber \end{align}

(2.15) \begin{align} T_t + uT_x + vT_y + wT_z &= T_{zz} + \tau {f}(\tau , T),\\[-6pt] \nonumber \end{align}

(2.16) \begin{align} \begin{pmatrix}u_z \\ v_z \end{pmatrix} = {f}(\tau , T)\frac {\boldsymbol{\tau }}{\tau } &,\quad \tau = \sqrt {\tau _{13}^2+\tau _{23}^2}, \end{align}

for $0\lt z\lt h$ , and

(2.17) \begin{align} T_t = \kappa T_{zz}, \quad \kappa = \frac {{\kappa _{_{\!R}}}}{{\kappa _{_{\!I}}}}, \end{align}

for $z\lt 0$ . The boundary conditions become

(2.18) \begin{align} h_t + u h_x + v h_y = w , \quad T = 0, \quad \textrm {at} \ z=h, \end{align}

(2.19) \begin{align} T^- = T^+, \quad KT_z^- = T_z^+, \quad \boldsymbol{u} = \boldsymbol 0 , \quad \textrm {at} \ z=0, \end{align}

and

(2.20) \begin{align} T_z \rightarrow -\varGamma /K, \quad \textrm {for} \ z\to -\infty , \end{align}

where

(2.21) \begin{align} \varGamma = \frac {{{\mathcal{H}}} {G}}{{[\Delta T]} {K_{_{\!I}}}} \quad \textrm {and} \quad K = \frac {{K_{_{\!R}}}}{{K_{_{\!I}}}} \end{align}

are the dimensionless geothermal flux and thermal conductivity ratio.

Moreover, the depth-integrated form of mass conservation, which determines the evolution of ice thickness, is

(2.22) \begin{align} h_t + \frac {\partial }{\partial x} \int _0^h u \; \textrm {d} z + \frac {\partial }{\partial y} \int _0^h v \; \textrm {d} z = 0 . \end{align}

In common with other lubrication-type gravity current models, the shallow-ice approximation retains only leading-order terms, discarding contributions such as horizontal diffusion and velocity gradients; validity of this assumption, and its possible limitations at short wavelengths, are discussed in the following sections.

2.3. Constitutive function

Glen’s flow law for non-isothermal ice is usually written in multiplicative form using an Arrhenius-type temperature dependence (Cuffey & Paterson Reference Cuffey and Paterson2010), giving

(2.23) \begin{align} {\hat {f}}(\hat T,\hat \tau )= A \exp \left (\!- \frac {E}{R\hat T}\right )\hat \tau ^{n}\!. \end{align}

Here, $E$ is the activation energy for creep, $R$ is the universal gas constant, $A$ is a constant rate factor and the exponent $n\approx 3$ . In dimensionless form, this becomes

(2.24) \begin{align} {f}(T,\tau )=\exp \left (\frac {T}{1+\delta T}\right )\tau ^{n}, \end{align}

where $\delta = R {T_a}/E$ . The Frank-Kamenetskii approximation takes $\delta \to 0$ , leaving ${f}=\tau ^n e^T$ .

In the following sections, we will use the form (2.24) for $f$ for illustrative numerical computations. We generally employ the full Arrhenius relation, but in some cases, when indicated, we also compare it with the Frank-Kamenetskii approximation. On the other hand, to maintain generality in the analytical results, we leave the form for this function open, and assume only that $f$ is a non-negative function that is strictly increasing in both $T$ and $\tau$ , implying

(2.25) \begin{align} {f}(T,\tau )\geqslant 0,\quad {f}_T(T,\tau )\gt 0,\quad {f}_{\tau }(T,\tau )\gt 0, \quad \text{for all} \quad T,\tau \geqslant 0, \end{align}

and ${f} \to 0$ as $\tau \to 0$ , where subscripts $T$ and $\tau$ once more denote differentiation. In practice, ensuring that ${f}_\tau \gt 0$ for all values of $\tau$ requires regularising the power law in (2.24), replacing $\tau ^n$ for instance by $(\tau ^2+\tau _0^2)^{(n-1)/2}\tau$ with a small regularising parameter $\tau _0 \ll 1$ . This becomes relevant in particular when dealing with the linearised version of the problem later, where we wish to avoid singular behaviour associated with vanishing shear stress at the ice surface.

2.4. Parameters

The independent dimensionless parameters in the model (2.14)–(2.22) are the thermal conductivity and thermal diffusivity ratios, $K$ and $\kappa$ , as well as the geothermal flux $\varGamma$ , the activation energy parameter $\delta$ and power-law index $n$ in the constitutive model (2.23), the slope parameter $s$ and the mean depth of the ice $H$ (cf. § 3). Glaciologically relevant values for these parameters are collated in table 1.

Table 1.

Typical dimensional and dimensionless parameter values.

Table 1. Long description

A table listing typical dimensional and dimensionless parameter values. The table has two columns and 16 rows. Column headers are symbols and units. Row 1: rho, 910 kilogram per cubic meter. Row 2: c_p, 2 times 10 to the power of 3 joule per kilogram kelvin. Row 3: A, 3.62 times 10 to the power of minus 15 per second per pascal to the power of 3. Row 4: E/R, 7.21 times 10 to the power of 3 kelvin. Row 5: T_a, 233 kelvin. Row 6: K_R, 3.1 watt per milli kelvin. Row 7: K_i, 2.18 watt per milli kelvin. Row 8: kappa_r, 1.6 times 10 to the power of minus 6 square meter per second. Row 9: kappa_i, 1.2 times 10 to the power of minus 6 square meter per second. Row 10: G, 0.06 watt per square meter. Row 11: H, 1.53 kilometer. Row 12: U, 2.64 meter per year. Row 13: L, 164 kilometer. Row 14: s, 0 to 2. Row 15: H, 0.1 to 1. Row 16: K, 1.42.

Recalling that $\mathcal H$ and $\mathcal L$ are chosen so that the Brinkmann and Péclet numbers are both unity, the dimensional parameter values in table 1 give $\mathcal H=1.53\,\textrm{km}$ and $\mathcal L=164\,\textrm{km}$ , as reported there. Although these should not be interpreted too literally as the depth and length of a particular glacier, they are of an order consistent with certain regions of the Greenland ice sheet.

In the remainder of the paper, we take $s$ and $H$ to be the main two control parameters of the ice-flow problem (because they are not material parameters but depend on the geometrical setting), treating all of $\{\kappa ,K,\varGamma ,\delta ,n\}$ as secondary parameters that we do not vary. As we will show in the following sections, the appearance of different regimes in the flux–slope relation for the base flow depends on the values of some of these parameters and, in extreme cases, certain regimes may disappear altogether. However, the qualitative stability properties associated with each regime remain robust.

3.1. Multiple branches of the base-flow solution

Sample solutions to the base-flow problem in (3.2)–(3.4) are presented in figure 2, which plots flux against slope for ${H}=(1/2)$ . For ${s}\ll 1$ and using the Arrhenius relation ( $\delta \neq 0$ ), the base-flow solution is

(3.6) \begin{align} {\varTheta } \sim \varGamma ({H}-z) \quad \textrm {and} \quad U \sim {s}^n {H}^{n+1} \int _{1-z/{H}}^1 \zeta ^n \exp \left ({\frac {\varGamma\!{H}\zeta }{1+\delta \varGamma\!{H}\zeta }}\right ) \textrm {d} \zeta , \end{align}

with a flux and basal temperature of

(3.7) \begin{align} Q &\sim {s}^{n} {H}^{n+2} \int _0^1 \zeta ^{n+1} \exp \left ({\frac {\varGamma\!{H}\zeta }{1+\delta \varGamma\!{H}\zeta }}\right ) \textrm {d} \zeta ,\\[-6pt] \nonumber \end{align}

(3.8) \begin{align} {\varTheta }(0) &\sim \varGamma\!{H} + {s}^{n+1} {H}^{n+3} \int _0^1 \zeta ^{n+2} \exp \left ({\frac {\varGamma\!{H}\zeta }{1+\delta \varGamma\!{H}\zeta }}\right ) \textrm {d} \zeta.\end{align}

As shown in figure 2, the corresponding solution branch only extends to ${s}\approx 2.1$ , at which point, the solution folds back on itself, generating a second branch of solutions with higher fluxes that decrease with increasing slope. At a second critical slope, the solution folds a second time, to generate a third branch with extremely high flux. On this branch, $\varTheta$ becomes so large that it is acceptable to take the limit ${\varTheta }\gg 1$ in (3.4)–(3.5), giving

(3.9) \begin{align} Q \sim \frac {s^n H^{n+2} e^{1/\delta }}{n+2} \quad \textrm {and} \quad {\varTheta }(0) \sim \varGamma\!{H} + \frac {s^{n+1} {H}^{n+3} e^{1/\delta }}{n+3}. \end{align}

The potential for the existence of multiple branches in the flux–slope relation was shown in previous studies (Clarke et al. Reference Clarke, Nitsan and Paterson1977; Yuen & Schubert Reference Yuen and Schubert1979; Fowler Reference Fowler1980; Fowler & Larson Reference Fowler and Larson1980; Hindmarsh Reference Hindmarsh2004). The dot-dashed lines in panels (a) and (b) indicate the melting-point threshold for the parameters in table 1. Above this line, the solution exceeds the melting point somewhere in the ice. In non-dimensional terms, this corresponds to $T \approx 5.31$ for the chosen parameters. Because the highest temperature occurs at the base, the threshold corresponds to a solution with basal temperature exactly at the melting point. Our model therefore breaks down above this threshold. For the parameter values considered here, branch 1 lies entirely below the threshold, and a portion of branch 2 does as well. Branch 3, however, lies entirely above it, so its analysis is purely mathematical rather than physically relevant.

Figure 2 also shows the temperature and velocity profiles of the three branches at ${s}=2$ . The unrealistic nature of branch 3 is evident in the extreme scale used in panel (d), whose magnitude is not physically plausible. The temperature field also illustrates the differing roles played by viscous dissipation, with the linear temperature profile in branch 1 (panel (c)) corresponding to negligible viscous heating and differing markedly from the concave profiles caused by significant dissipation on the other two branches. An asymptotic solution for the temperature profile for large basal temperatures such as those in branch 2 (panel (c)) can be found in Appendix C, showing that the concave part of the profile is in fact a boundary layer near the base of the ice, underlying a conduction-dominated region that extends across the remainder of the ice column.

Figure 3 shows how the basal temperature varies with slope, and compares the results using the Arrhenius relation with those for the Frank-Kamenetskii approximation. Notably, with the Frank-Kamenetskii approximation, the third branch vanishes, and the second branch continues monotonically to smaller slopes and higher basal temperatures. The figure further includes the asymptotes at small values of ${\varTheta }(0)-\varGamma\!{H}$ , from (3.8), and large values of ${\varTheta }(0)-\varGamma\!{H}$ , from (3.9).

Figure 3.

Relative base temperature, $\varTheta (0) - \varGamma H$ , plotted against the slope $s$ for ${H} = ({1}/{2})$ and the remaining parameter settings provided in table 1. The blue solid line represents the results obtained using the Arrhenius relation, while the black dashed line corresponds to those with the Frank-Kamenetskii approximation. The dot-dashed line shows Joseph’s bound (Joseph Reference Joseph1966) and the thicker red dotted lines indicate the asymptotes for small and large $\varTheta (0) - \varGamma H$ .

Analysis by Joseph (Reference Joseph1966) can be used to place bounds on the regions where there are multiple solutions. We translate Joseph’s analysis to the current problem in Appendix B, which shows that there can be no solutions to the right of the dot-dashed line added to figure 3. Note that a key element of Joseph’s bound is that a certain function of the base temperature must remain positive (see Appendix B). This is violated for large ${\varTheta }(0)$ with the Arrhenius law. The bound therefore disappears at sufficiently large base temperatures, allowing the upper branch to pass to the right of the dot-dashed line.

Figure 4(a) presents diagrams of $Q$ against $s$ for various values of the base height $H$ . As the base height increases, the fold points move closer together until they collide at a codimension-2 cusp bifurcation, and vanish for larger base heights: beyond that threshold, the second branch disappears, and the first and third branches merge. Figure 4(c) is similar to figure 4(a), except that contours of $H$ are shown in the $(s,Q)$ -plane.

Figure 4.

(a) Flux–slope relations for various values of the mean height $H$ , with the remaining parameter values given in table 1. The dotted red line tracks the fold points. The blue dot-dashed line shows a trajectory on the $(s,H,Q)$ surface corresponding to the non-uniform steady state computed from the two-dimensional problem with a mass source term (3.10) with $a(x)=7(1-x^4)$ ; details are given in the Appendix D. (b) Height (blue dot-dashed line) and flux (solid black line) as functions of $x$ for this non-uniform steady state. (c) The trajectory associated with the non-uniform steady state, plotted over contours of $H$ in the $(s,Q)$ -plane. The blue square marks the ice divide, the blue star marks the point where the flux attains its maximum and the blue triangle marks the margin. These markers are included in every panel for reference.

3.2. Quasi-static flux

As we have seen, the steady-state base solution predicts a flux through (3.4) that depends on the slope $s$ and ice depth $H$ : $Q=\mathcal{Q}({s},{H})$ . If this flux law were to extend to the full, time-dependent, two-dimensional problem with a mass source term $a$ , then we could write

(3.10) \begin{align} h_t + Q_x = a, \quad Q=\mathcal{Q}({S},h), \quad {S}={s} - h_x. \end{align}

The simplified system (3.10) permits a straightforward exploration of linear stability. For simplicity omit the source term $a$ once more and perturb around a uniform steady state by setting $h = {H} + {\eta }(x,t)$ , we obtain

(3.11) \begin{align} {\eta }_t - \mathcal{Q}_s{\eta }_{xx} + \mathcal{Q}_{H} {\eta }_x = 0, \end{align}

where

(3.12) \begin{align} \mathcal{Q}_s = \frac {\partial \mathcal{Q}}{\partial s} \quad \text{and} \quad \mathcal{Q}_{H} = \frac {\partial \mathcal{Q}}{\partial {H}} \end{align}

denote the partial derivatives of the flux law, evaluated for the base flow. The sign of $\mathcal{Q}_s$ therefore dictates stability, implying that the first and third branches are expected to be stable. The intermediate second branch is expected to be unstable, and, indeed, the quasi-static approximation suggests an ill-posed backward-diffusion problem in which the growth rate diverges quadratically with wavenumber.

As shown in Appendix D, after including a prescribed mass source term $a(x)$ , the system (3.10) also allows a non-uniform steady state with compact support to be computed on a flat bed (with $s = 0$ ), starting at a central ice divide ( $x=0$ ) and imposing vanishing ice thickness and flux at a free boundary representing the lateral edge or margin of the ice. A sample solution is shown in figures 4(b) and 4(c). At the ice divide, the flux vanishes by symmetry and surface slope vanishes, so the trajectory begins on branch 1. As $x$ increases away from the divide, marked by the blue square, the flux initially increases as $a \gt 0$ . The local surface slope $S = -h_x$ increases as well, so the trajectory moves upward along branch 1 until it crosses into branch 2. At the point where the flux attains its maximum, marked by the blue star, the trajectory turns and then proceeds toward smaller values of $Q$ as the profile thins toward the margin, eventually returning to branch 1. Finally, at the margin, both $H$ and $Q$ vanish, and the trajectory terminates at the blue triangle. In this way, the trajectory provides a convenient geometric description of how the non-uniform steady state samples the local slab solutions and whether it remains on branch 1 or crosses the fold and enters branch 2. Clearly, steady-state solutions that cross onto branch 2 are possible (figure 4), making the properties of both branches 1 and 2 relevant to less idealised settings than the simple parallel-sided slab of (3.1)–(3.3).

The same stability problem (3.11) with the last term $\mathcal{Q}_{H} {\eta }_x$ omitted, arises in the short-wave perturbation limit of the non-trivial steady state constructed in Appendix D and displayed in figure 4. The stability properties remain the same, with the non-uniform steady-state solution being violently unstable to short-wavelength perturbations where the solution follows branch 2 of the flux–thickness–slope relationship. However, the underlying assumption of a quasi-static, low-Péclet-number temperature field ignores all effects of thermal relaxation, prompting us to turn to a more detailed analysis below that shows why the simple intuition we may gain from linearising a steady-state flux law as above is inaccurate.

4.1. Zero-wavenumber modes; $k=l=0$

The results from the quasi-static stability arguments of § 3.2 can be placed on a firmer basis by looking for solutions to the linear stability problem that depend only on $z$ . That is, we set $k=l=0$ , reducing (4.2)–(4.10) to

(4.11) \begin{align} \sigma {\theta }- {\theta }_{zz}- {F}_T{\bar \tau } {\theta }=0,& \quad 0\lt z\lt {H},\\[-6pt] \nonumber \end{align}

(4.12) \begin{align} {\theta }_{z}-{K}{\sqrt {\frac {\sigma }{\kappa }}}{\theta }=0,& \quad \text{at} \quad z=0, \\[-6pt] \nonumber\end{align}

(4.13) \begin{align} {\theta }=0,& \quad \text{at} \quad z={H}. \end{align}

We will make repeated use later of the fact that a closely related problem is obtained by taking the equations (3.2)–(3.3) for the base flow and differentiating with respect to the parameter $Q$ , in which case

(4.14) \begin{align} -{\varTheta }_{{Q} zz} - {F}_T {\bar \tau } {\varTheta }_{{Q}} - {s}_{{Q}} ({F}_{\tau }{\bar \tau }+{F}) ({H}-z) &=0, \quad 0 \lt z \lt {H},\end{align}

(4.15) \begin{align} {\varTheta }_{{Q} z} &= 0, \quad \text{at} \quad z = 0,\end{align}

(4.16) \begin{align} {\varTheta }_{{Q}} &= 0, \quad \text{at} \quad z = {H},\end{align}

where ${\varTheta }_{Q}$ is the derivative of the base-flow solution with respect to flux $Q$ . Moreover, substitution of (3.2) into (3.4) yields

(4.17) \begin{align} Q=-\frac {1}{s}\int _0^H\varTheta _{zz}\text{d}z=-\frac {1}{s}\left (\varTheta _z(H)+\varGamma \right ). \end{align}

Multiplying this expression by $s$ and differentiating with respect to $Q$ yields the additional condition

(4.18) \begin{align} \varTheta _{Qz}=-s_Q Q-s \quad \text{at} \quad z = H. \end{align}

In the eigenvalue problem (4.11)–(4.13), the lower boundary condition ensures that this system is not of Sturm–Liouville type (Zettl Reference Zettl2005). However, the problem does take such a form if we consider the limit $K\to 0$ , which is equivalent to imposing a constant heat flux at the base of the ice rather than solving a diffusion problem. We consider this special case first. In doing so, we also reconsider the parametric dependence of the base-flow solutions. In particular, because these solutions vary non-monotonically with $s$ , we recast the base-state problem as one that is parametrised (implicitly) by $Q$ , rather than $s$ , which renders the parametric dependence single valued. The slope then becomes an output parameter of the base-flow problem, easily determined by reflecting plots like figure 2(a) about the line $Q = s$ .

4.1.1. The Sturm–Liouville problem

The $k=l=K=0$ Sturm–Liouville problem is

(4.19) \begin{align} \sigma {\theta } - {\theta }_{zz} - {F}_T {\bar \tau } {\theta } &= 0, \quad 0 \lt z \lt {H}, \\[-6pt] \nonumber\end{align}

(4.20) \begin{align} {\theta }_z &= 0, \quad \text{at} \quad z = 0,\\[-6pt] \nonumber \end{align}

(4.21) \begin{align} {\theta } &= 0, \quad \text{at} \quad z = {H}. \end{align}

Fold points occur when $s_Q = 0$ . Comparing (4.19)–(4.21) with (4.14)–(4.16) at a fold point, we see that ${\varTheta }_{Q}$ satisfies (4.19)–(4.21). Moreover, (4.18) shows that $\varTheta _Q$ is not identically zero. It follows that there must be a zero eigenvalue at a fold point.

By multiplying (4.19) by ${\varTheta }_{Q}$ , (4.14) by $\theta$ and integrating, we also find that

(4.22) \begin{align} \sigma \int _0^{{H}}{\theta }{\varTheta }_{{Q}}\text{d}z= - {s}_{Q} \int _0^{{H}}{\theta }({F}_{\tau }{\bar \tau }+{F}) ({H}-z) \text{d}z. \end{align}

Sturm–Liouville theory guarantees that eigenfunctions associated with the largest eigenvalue maintain a single sign. Our choice of constitutive function also implies that ${F}_{\tau }{\bar \tau }+{F}\gt 0$ . Therefore, (4.22) confirms that a fold point must occur in the bifurcation diagram whenever the largest eigenvalue vanishes. For low values of the base flux $({Q}\ll 1)$ , the third term in (4.19) can be omitted, resulting in a stable system with largest eigenvalue $\sigma =-(\pi /{H})^2$ . Consequently, the solution along the first branch must be stable for modes with $k=l=K=0$ .

A further derivative with respect to $Q$ leads to the relation

(4.23) \begin{align} \sigma _{{Q}}\int _0^{{H}}{\varTheta }_{{Q}}{\theta }\text{d}z = - {s}_{{Q}{Q}} \int _0^{{H}}({F}_{\tau }{\bar \tau }+{F}) ({H}-z) {\theta }\text{d}z , \end{align}

at fold points (since $s_Q = 0$ there). From (4.16) and (4.18), we conclude that, at the fold point, $\varTheta _Q$ is positive in a neighbourhood of $z=H$ . Moreover, at the fold point, $\varTheta _Q$ is the eigenfunction associated with the largest eigenvalue and therefore does not change sign. As a result, $\varTheta _Q\geqslant 0$ throughout the domain. Hence $\sigma _{{Q}}\gt 0$ , implying that $\sigma$ becomes positive along the second branch. This stability pattern is illustrated in figure 2, where the solid line represents stable solutions, while the dashed line indicates unstable solutions.

Figure 5.

Growth rates of the largest eigenvalues as a function of $Q$ . The solid line represents the growth rate from the Sturm–Liouville system ( $k=l=K=0$ ) described by (4.19)–(4.21), and the dashed lines correspond to the system described by (4.11)–(4.13) for $K=6.4,3.2,1.6,0.8,0.4,0.2,0.1$ . These computations were performed for $H=(1/2)$ and the remaining parameter settings listed in table 1. The vertical dot-dashed lines indicate the locations of the fold points, and the horizontal dot-dashed line indicates zero growth rate. Also shown is the slope of the base-flow solution.

Figure 5 shows the growth rate of numerical solutions to (4.19)–(4.21) across the three branches, illustrating the transition of stability at each fold point. At the second fold point, we observe numerically that the leading eigenvalue again reverses sign, implying a return to stability along the third branch. This feature cannot be established directly from the preceding arguments because the second fold point could in principle correspond to the vanishing of either the largest or the second-largest eigenvalue.

Fowler & Larson (Reference Fowler and Larson1980) previously analysed the Sturm–Liouville problem in (4.19)–(4.21). Without the benefit of numerical solutions to confirm their results, Fowler & Larson conclude that the base flow is always stable, even along the second branch. The root of their disagreement with our results lies in the properties of the eigenfunctions of the Sturm–Liouville problem (our $\theta$ , denoted by $\chi$ in their notation). Fowler & Larson (Reference Fowler and Larson1980) claim that $\chi$ has self-similar structure under changes in a control parameter $\alpha$ (see their (29)); $\alpha$ is equivalent to $s^{n+1}$ here, where $s$ is slope and $n$ the exponent in (2.24). The function $\chi$ , however, satisfies a particular version of our (4.19) ((24) in Fowler & Larson (Reference Fowler and Larson1980)) in which the base profile ${\varTheta }(z)$ (denoted by $\theta$ in their work) appears as a coefficient. The self-similar representation in their (29) is built on a rescaling of the vertical coordinate, $\xi \mapsto (\beta /\alpha )^{1/(n+3)}\xi$ , and holds only if the base profile is invariant under this rescaling (only if $\theta (\xi )$ and $\theta \! ((\alpha /\beta )^{1/(n+3)}\xi )$ coincide), which is not the case for the non-constant profiles at hand. This invalidates the remainder of their stability argument.

Indeed, our results above establish that the second branch is unstable, as already indicated by the numerical work in Hindmarsh (Reference Hindmarsh2004). Moreover, the linear stability analysis of the Couette, Poiseuille and Hagen–Poiseuille configurations, conducted by Joseph (Reference Joseph1964), yields results consistent with ours. Lastly, our findings are in agreement with those of Joseph (Reference Joseph1965a ), who investigated the linear stability of temperature distributions in conducting solids with nonlinear heat generation. Although his model involves a heat generation term dependent only on temperature, his extension to non-planar geometries introduces a spatial dependence to the heating term, making it identical to our formulation.

4.1.2. Finite bedrock diffusion

We return now to the more general, non-self-adjoint problem (4.11)–(4.13), retaining diffusion in the bedrock. Multiply (4.11) by the complex conjugate of $\theta$ and integrate, to find on application of the boundary conditions that

(4.24) \begin{align} \sigma \int _0^{{H}}|{\theta }|^2\text{d}z+\frac {K}{\sqrt {\kappa }}\sqrt {\sigma }|{\theta }(0)|^2 +\int _0^{{H}}|{\theta }_{z}|^2\text{d}z=\int _0^{{H}}{F}_{T }{\bar \tau }|{\theta }|^2\text{d}z. \end{align}

Assuming $\sigma$ has zero real part, so that $\sigma =i\sigma _{i}$ with $\sigma _{i}\in \mathbb{R}$ , and taking the imaginary part of (4.24), we obtain

(4.25) \begin{align} |\sigma _{i}|\int _0^{{H}}|{\theta }|^2\text{d}z +\frac {K}{\sqrt {2\kappa }}\sqrt {|\sigma _{i}|}|{\theta }(0)|^2=0, \end{align}

thus $\sigma _{i}=0$ . In other words, if $\sigma$ crosses the imaginary axis, it must do so through the origin. At a zero eigenvalue, the eigenvalue problems (4.11)–(4.13) and (4.19)–(4.21) are, however, identical, and we can therefore conclude that we can still only have $\sigma = 0$ in the non-self-adjoint problem (4.11)–(4.13) when $s_Q = 0$ as per (4.22). A change of stability at zero wavenumber can only happen at a fold point even when bedrock diffusion is retained in the problem.

Similarly to § 4.1.1, we can show that

(4.26) \begin{align} \int _0^{{H}}[\lambda ^2{\varTheta }_{{Q}}+ {s}_{Q} ({F}_{\tau }{\bar \tau }+{F}) (H-z)]{\theta }\text{d}z =-\frac {K}{\sqrt {\kappa }}\lambda {\theta }(0){\varTheta }_{{Q}}(0), \end{align}

where $\lambda =\sqrt {\sigma }$ has positive real part. Differentiating with respect to $Q$ at the fold ( $\lambda ={s}_{{Q}}=0$ ) results in

(4.27) \begin{align} {s}_{{Q}{Q}} \int _0^{{H}}({F}_{\tau }{\bar \tau }+{F}) (H-z){\theta } \text{d}z = -\frac {K}{\sqrt {\kappa }}{\theta }(0){\varTheta }_{{Q}}(0)\lambda _{{Q}}. \end{align}

At fold points, the eigenfunction $\theta$ corresponding to $\sigma = 0$ satisfies the same problem as ${\varTheta }_{Q}$ , and hence ${\theta }\propto {\varTheta }_{Q}\gt 0$ . Thus, ${\theta }(0)\neq 0$ (since otherwise the problem is homogeneous with two homogeneous boundary conditions, admitting only the trivial solution ${\theta } \equiv 0$ ), $\lambda _{{Q}}$ is real and

(4.28) \begin{align} \text{sgn}\left ( \lambda _{{Q}} \right ) = \text{sgn} (-{s}_{{Q}{Q}}).\end{align}

Therefore, on passing the first fold point (where $-{s}_{{Q}{Q}} \gt 0$ ), the real part of $\lambda$ (and $\sigma$ ) becomes positive and remains so along the second branch. On the first branch below the fold, however, this argument implies that $\textrm{Re} (\lambda )$ must be negative, which contradicts the requirement in the boundary conditions that $\lambda$ has a positive real part. Therefore, the problem has no solutions that decay for $z\to -\infty$ . The numerical solutions shown in figure 5 suggest that the unstable mode intersects and disappears into the continuous spectrum of non-spatially decaying modes at the fold points, consistent with the observation that we cannot find a solution for $\lambda$ with positive real part.

4.2. Modes without transverse structure $(l=0)$

For perturbations with streamwise spatial structure ( $k\gt 0$ ), the range of unstable modes is much larger. In particular, when $K\gt 0$ , all three branches become unstable over a finite range of wavenumbers $k$ . The instability along the lowest and intermediate branches is illustrated in figures 6 and 7(a), where unstable modes extend from $k=0$ up to an $O(1)$ critical wavenumber. Along the lower branch, $\sigma _r := \textrm{Re} (\sigma ) \propto k^{3/2}$ for $k\to 0$ . In contrast, $\sigma _r$ converges to the $O(1)$ growth rate of the zero-wavenumber problem considered in § 4.1.2 for the intermediate branch (figure 6 b).

Figure 6.

(a) Growth rate of the most unstable mode as a function of streamwise wavenumber $k$ , for three values of the flux along the first $({Q}=10^{-0.1},10^{0},10^{0.1}$ ; solid) and second $({Q}=10^{0.2},10^{0.3},10^{0.4}$ ; dot-dashed) branches, with ${H}=( {1}/{2})$ and the parameters specified in table 1. The grey dashed line marks the zero axis. The growth rates are replotted using logarithmic axes in (b).

The stability properties of the lowest branch and the low-wavenumber behaviour are different when diffusion in the bedrock is ignored and a constant heat flux is imposed at the base of the ice ( $K=0$ ). In this case, the lowest branch is stable for sufficiently small flux (or slope), and a transition to instability arises at a critical value of $Q$ that lies below the flux of the first fold point; see figure 7(b). Between that critical flux and the fold point, a window of unstable wavenumbers appears, with $\sigma _r\propto k^2$ for $k\to 0$ .

The two variants of long-wave instability for $K\gt 0$ and $K=0$ can be confirmed using asymptotic expansion for $k\ll 1$ . We provide the analysis for the $K \gt 0$ case below.

Figure 7.

Contours of constant growth rate on the $(Q,k)$ -plane, for (a) $K=1.42$ and (b) $K=0$ , with ${H}=1/2$ , while the remaining parameters are specified in table 1. The (black) dot-dashed line indicates the first fold point. To the left of this line, the red dotted contours correspond to the long-wave asymptotic predictions, while to the right they represent the zero-wavenumber modes described in § 4.1.2. In each case, the growth rate along the dotted curves is matched to the corresponding contour value.

4.2.1. Long-wave instability of the first branch

To explore linear stability in the long-wave limit, we use the asymptotic expansions

(4.29) \begin{align} \sigma &= k\sigma _0+k^{3/2}\sigma _{1/2}+O(k^2), \quad {\theta }={\theta }_0+k^{1/2}{\theta }_{1/2}+O(k), \quad u'=u'_0+k^{1/2}u'_{1/2}+O(k),\nonumber\\ w' &= kw'_0+k^{3/2}w'_{1/2}+O(k^2), \quad Q'=Q'_0+k^{1/2}Q'_{1/2}+O(k). \end{align}

In view of our choice of normalisation of the eigenvalue problem ( ${\eta }=1$ ), the perturbation to the ice thickness $\eta$ does not need to be expanded. At leading order we obtain

(4.30) \begin{align} - {\theta }_{0,zz}&={F}_T{\bar \tau } {\theta }_0+({F}_{\tau }{\bar \tau }+{F})\tau '_0, \\[-6pt] \nonumber\end{align}

(4.31) \begin{align} u'_{0,z}&={F}_{\tau }\tau '_0+{F}_T {\theta }_0, \\[-6pt] \nonumber\end{align}

(4.32) \begin{align} \tau '_0&={\eta }{s}, \end{align}

for $0\lt z\lt {H}$ , with

(4.33) \begin{align} {\theta }_{0,z}=u'_0&=0, \quad \text{at} \quad z=0, \\[-6pt] \nonumber\end{align}

(4.34) \begin{align} {\theta }_{0}+{\varTheta }_z{\eta }&=0, \quad \text{at} \quad z={H}, \end{align}

and

(4.35) \begin{align} \sigma _0 {\eta }+iQ'_0=0,\quad Q'_0=\int _0^{{H}}u'_0\text{d}z+{U}({H}){\eta }. \end{align}

By differentiating the base-flow equations with respect to the parameter $H$ , while keeping all other parameters fixed, we derive the following relationship:

(4.36) \begin{align} \frac {1}{{\eta }}\left [{\theta }_0, u'_0, Q'_0\right ] = \left [{\varTheta }_{{H}}, {U}_{\!{H}}, {Q}_{\!{H}}\right ]. \end{align}

This result implies that the leading-order eigenvalue is $\sigma _0=-i{Q}_{\!{H}}$ , which is purely imaginary ( $Q_{{H}}$ is real). Hence, we must proceed to next order to assess the stability of the system.

At $O(k^{1/2})$ , we find

(4.37) \begin{align} - {\theta }_{1/2,zz}={F}_T{\bar \tau } {\theta }_{1/2},\quad u'_{1/2,z}={F}_T {\theta }_{1/2}, \end{align}

with

(4.38) \begin{align} {\theta }_{1/2,z}(0)=K\sqrt {\sigma _0/\kappa }\; {\theta }_0(0), \quad u'_{1/2}(0)=0, \end{align}

(4.39) \begin{align} {\theta }_{1/2}({H})=0, \end{align}

and

(4.40) \begin{align} \sigma _{1/2}{\eta }+i\int _0^{{H}}u'_{1/2}\text{d}z=0. \end{align}

Integrating by parts

(4.41) \begin{align} \int _0^{{H}}u'_{1/2}\text{d}z = \int _0^{{H}}\int _0^z{F}_T(\zeta ) {\theta }_{1/2}(\zeta )\text{d}\zeta \text{d}z =\int _0^{{H}}({H}-z){F}_T(z) {\theta }_{1/2}(z)\text{d}z. \end{align}

Thus,

(4.42) \begin{align} \sigma _{1/2}=-\frac {i}{{s} {\eta }}\int _0^{{H}}{F}_T {\bar \tau } {\theta }_{1/2}\text{d}z , \end{align}

since ${\bar \tau }/{s}={H} -z$ . The key is therefore to compute the integral on the right, or more specifically, to find the sign of its imaginary part. This requires us to determine the sign of the imaginary part of ${\theta }_{1/2}$ , given that ${F}_T{\bar \tau }\gt 0$ and ${\eta }=1$ . It turns out that $\textrm{Im} ({\theta }_{1/2})\gt 0$ along the first branch, a result that is established in Appendix A. Numerically we find similar behaviour along the third branch.

Figure 8.

Numerically computed growth rates plotted against $Q$ (solid blue lines) and compared with the predictions from long-wave asymptotics (green lines and circles) and the zero-wavenumber modes of § 4.1.2 (red lines and diamonds), for three values of $k$ (as indicated), $H=(1/2)$ and the other parameter settings listed in table 1.

Figure 8 compares the predictions of the long-wave analysis with numerical computations of $\sigma _r$ . Also shown are the results from the zero-wavenumber analysis of § 4.1.2. One might assume that the instabilities of the first and second branches arise from different modes owing to the different behaviour in the limit $k\to 0$ on the two branches. Figure 8 shows that this is not the case, and that the most unstable modes are actually connected across the first fold point. This connection can be established, and the two disparate scalings with $k$ reconciled, because the long-wave growth rate from (4.42) diverges as we approach the fold from the left, whereas the zero-wavenumber growth rate of the second branch vanishes on approaching from the right (see figure 8).

Note that the comparison between the asymptotic predictions and numerical results worsens for $k=0.01$ and ${Q}=O(k)$ in figure 8, because, at these fluxes, there is a cross-over to a different asymptotic limit that applies for $Q \ll k \ll 1$ . Also, the result in (4.42) highlights how the long-wave instability is insensitive to the detailed rheology of the ice. Indeed, the fluid could actually be Newtonian ( $n=1$ ) with a temperature-dependent viscosity, and the flux–slope relation of the base states can be monotonically increasing and single valued (as is the case at small mean depths in figure 4), and the wave-like instability mechanism still persists.

4.2.2. Physical mechanism

The long-wave analysis above also helps to clarify the physical mechanism behind the instability: as illustrated in figure 9, flux increases with ice depth, causing regions near the crest of a wave to flow faster than regions near a trough. Consequently, areas where ice depth decreases downstream tend to accumulate material; areas of increasing ice depth tend to thin. The ice surface grows or shrinks to accommodate these variations, causing a wave to travel downstream. This wave is captured by the leading-order relations in the long-wave expansion in (4.30)–(4.35).

The presence of the wave affects the heat exchange between the ice and the bedrock, with excess heat generated within the ice underneath the wave crests, and diffusing into the underlying bedrock. This is the physical content of the $O(k^{1/2})$ corrections in (4.37)–(4.40). But the heat diffusion into the bed lags behind the wave crests (see figure 9), leading to additional heat retained in the ice on the upstream side of the wave, and cooler temperatures persisting below the downstream side of the wave. The combination reduces the viscosity where the layer is thinning, and enhances it where the layer thickens. This positive feedback, represented by (4.42) and facilitated by the diffusion of heat through the bedrock, causes the instability.

Figure 9.

Colour map of the temperature perturbation $\theta$ on the $(kx,z)$ -plane, for ${H}=(1/2)$ , $k=0.1$ , $Q=1$ , $t=0$ and the parameters in table 1. Black contours show the horizontal velocity perturbation $u'$ , and we have included an exaggerated surface perturbation to illustrate the relevant phase shifts.

4.2.3. Short-wave perturbations

As observed in figures 6 and 7, modes with zero transverse wavenumber $l$ in general become unstable only if the streamwise wavenumber lies within a finite range (whose size depends on model parameters), and large wavenumbers $k$ are invariably stable. This contrasts markedly with the case of temperature-dependent sliding at the base of the ice (Mantelli & Schoof Reference Mantelli and Schoof2019), in which shallow-flow models remain unstable even at arbitrarily short streamwise wavelengths. To clarify the origin of this difference, we give an abbreviated analysis of the relevant asymptotic limit $k \gg 1$ here.

The relevant leading-order behaviour is $\sigma =O(k^2)$ , ${\theta }=O(1)$ , $u'= O(k)$ , $w'= O(k^2)$ , $\tau ' = O(k)$ , $Q'= O(k)$ with $\eta = 1$ , yielding, at leading order,

(4.43) \begin{align} \sigma {\eta }+i k Q'=0,\quad Q'=\int _0^{{H}}u'\text{d}z,\quad u'_z={F}_{\tau }\tau ',\quad \tau '=-ik ({H}-z){\eta }, \end{align}

with $u'=0$ at $z=0$ . This leads to

(4.44) \begin{align} \sigma =-k^2 \int _0^{{H}}({H} -z)^2{F}_{\tau }({\bar \tau }(z),{\varTheta }(z)) \text{d}z\lt 0, \end{align}

which decouples entirely from the thermal problem: temperature appears only parametrically through the base solution $\varTheta$ in (4.44). The leading-order behaviour here is simply that of a diffusive highly viscous gravity current model, with a negative growth rate whose magnitude grows quadratically with wavenumber.

The fact that the flow problem can be decoupled from the thermal problem hinges on the leading-order versions of (4.2), (4.5) and (4.9) furnishing a suitable solution for $w'$ and $\theta$ . The heat equation (4.2) demands that $\theta$ be determined by advection only at leading order, so that

(4.45) \begin{align} {\theta } = -\frac {w' {\varTheta }_z}{\sigma }, \end{align}

while (4.9) demands at leading order that ${\theta } = - {\varTheta }_z \eta$ at $z = H$ and ${\theta } = 0$ at $z = 0$ . Ordinarily, we would expect boundary layers to appear at $z = 0$ and $z = H$ , having reduced the second-order heat equation to an algebraic expression in (4.45). Indeed, thermal boundary layers are the cause of the short-wave instabilities in § 4 of Mantelli & Schoof (Reference Mantelli and Schoof2019) for the case of temperature-dependent basal sliding. Here, in the absence of sliding and of heat dissipation at the bed itself, the purely advective solution (4.45) does not run into such problems: the vertical velocity, given by

(4.46) \begin{align} w'(z) = - ik \int _0^z u'(z') \textrm {d} z', \end{align}

satisfies $w' = 0$ at $z = 0$ , and hence the boundary condition ${\theta } = 0$ holds there. Moreover, at $z = H$ , $w' = -ikQ' = \sigma \eta$ and the surface boundary condition ${\theta } - w' {\varTheta }_z/\sigma = - {\varTheta }_z \eta$ also holds automatically.

4.3. Three-dimensional perturbations ( $k,l\gt 0$ )

The diffusive behaviour of the problem in one horizontal dimension deduced in § 4.2.3 is in fact deceptive, as we see when considering fully three-dimensional perturbations. Sample growth rates for such perturbations are plotted in the $(k,l)\text{-}$ plane in figures 10(a) and 10(d). The remaining panels present the corresponding results for the two components of the phase velocity, $-\textrm{Im} (\sigma )\boldsymbol{k}/|\boldsymbol{k}|^2$ (where $\boldsymbol{k} = (k,l)$ is the usual wave vector), each scaled by the mean speed of the base flow, $Q/H$ . Panels (a,b,c) shows the results for an example along the first branch; panels (d,e,f) illustrates a case along the second branch. For the first branch, we recover the window of unstable wavenumbers in $k$ for low $l$ ; increasing the transverse wavenumber $l$ eventually suppresses the instability, rendering stable all the modes with higher $l$ . By contrast, although an increase in transverse wavenumber initially still suppresses the unstable window in $k$ along the second branch, a different regime of instability appears at higher $l$ . This suppression can be traced to the transverse flux in the depth-integrated mass conservation (4.10): a transverse thickness perturbation drives a transverse flux that acts diffusively on $\eta$ , contributing a stabilising term $\propto -l^2$ to the growth rate, even though no explicit transverse diffusion appears in the model. Stability boundaries, illustrating the two regimes for several values of flux $Q$ , are shown in figure 11.

Figure 10.

Colour maps over the $(k,l)$ -plane for (a,d) growth rate $\textrm{Re} (\sigma )$ , and (b,e) streamwise and (c,f) cross-stream phase speed, of the most unstable mode at (a,b,c) $Q=1$ and (d,e,f) $Q=2$ ; $H=(1/2)$ and the other parameter settings listed in table 1. The phase speeds, $-\textrm{Im} (\sigma )(k,l)/(k^2+l^2)$ , are scaled by the mean speed of the base flow, $Q/H$ . Note that in (d) the growth rate plateaus at a finite value ( $\approx 0.1$ , dark blue) as $l$ increases. Note also the logarithmic colour scale; regions in white correspond to negative growth rates.

Figure 11.

Stability boundaries on the $(k,l)$ -plane for different values of the base flux along the first $({Q}=10^{-0.1},10^{-0.05},10^{0},10^{0.05},10^{0.1},10^{0.15}$ ; red, solid) and second branches $({Q}=10^{0.2},10^{0.25},10^{0.3}, 10^{0.35}, 10^{0.4},10^{0.45}$ ; blue, dashed) for $H=( 1/2)$ and the remaining parameter settings listed in table 1.

The two different unstable regimes were identified previously by Hindmarsh (Reference Hindmarsh2004). As in Hindmarsh’s work, we will continue to refer to them as regimes A and B. Regime A appears for all branches and is characterised by small to order-one values of both wavenumbers. Regime B arises only for the intermediate branch, and is characterised by a flat plateau in growth rate at higher transverse wavenumbers and small to order-one values of $k$ . As illustrated in figure 10, the unstable modes of regime A have significant phase speeds (up to six times the average base-flow speed). By contrast the phase speeds of the unstable modes in regime B are much lower, and approach zero for large enough transverse wavenumber $l$ .

Figure 10 also shows that regime A typically features higher growth rates than regime B. This property is illustrated further in in figure 12(a), which plots the ratio of the maximum growth rate in regime B to that in regime A as a density over the $(H,Q)\text{-}$ plane. The parameter space shown in this figure corresponds to base flows along the second branch for which the base temperature lies below the melting point of ice. Figure 12(b) also plots the streamwise wavenumber $k$ at which maximum growth occurs (in regime A). The maximum growth rate typically occurs for $l=0$ and $k=\text{ord}(1)$ , except at the smallest mean depths, where $k$ for the most unstable mode becomes larger.

Figure 13 shows how the instability varies across the $(H,Q)$ -plane along the first and second branches, with the red dashed line marking the fold points. Figure 13(a) displays the maximum growth rate over all wavenumbers. The fold-point curve is clearly associated with the contours of largest growth, indicating that the strength of the instability is strongly linked to the branch structure. The instability is weaker on the first branch, but strengthens as one approaches the fold point and crosses onto the second branch. There, the maximum growth rate continues to increase, reaches a peak and then decreases. On the other hand, the maximum growth rate decreases monotonically as $H$ increases. In all cases, this maximum occurs at $l=0$ , and figure 13(b) shows the corresponding phase speed scaled by the mean speed of the base flow, $-H\textrm{Im} (\sigma _{max })/(kQ)$ . Over the parameter range explored, this scaled phase speed varies only moderately, taking values between approximately 6 and 10, and increasing as $Q$ decreases and $H$ increases. The more significant variation in scaled phase speed occurs between regimes A and B; however, since the maximum overall growth rate always belongs to regime A, this difference remains moderate in the present figure.

Figure 12.

Contours on the $(H,Q)\text{-}$ plane of (a) the ratio of maximum growth rate in regime B to that in regime A, and (b) the streamwise wavenumber $k$ maximising growth. The white region below the contours corresponds to the first branch (where there is no regime B); that above locates where base temperatures surpass the melting point.

In view of the resolution-dependent thermal spokes in the numerical solutions of Payne & Baldwin (Reference Payne and Baldwin2000), the most pertinent feature of the three-dimensional dispersion relation is, however, not the fastest growth rate, but the fact that growth rates remain positive and bounded away from zero for arbitrarily large $l$ in regime B. This motivates us to delve further into the relevant short-wave asymptotics for transverse perturbations.

For $l \gg 1$ and $k = O(1)$ in regime B, numerical experiments suggest that the solution to (4.2)–(4.10) behaves as $\sigma = O(1)$ , $\theta = O(1)$ , $u' = O(1)$ , $v' = O(l^{-1})$ , $w = O(1)$ , $Q' = O(1)$ , $\tau ' = O(l^{-2})$ with $\eta = O(l^{-2})$ , resulting in a leading-order problem of the form

(4.47) \begin{align} \sigma {\theta } + i k {U} {\theta } + w' {\varTheta }_z & = {\theta }_{zz} + {F}_T {\bar \tau } {\theta }, \\[-6pt] \nonumber\end{align}

(4.48) \begin{align} u_z' & = {F}_T {\theta }, \\[-6pt] \nonumber\end{align}

(4.49) \begin{align} v_z' & = - i l s^{-1}{F} \eta , \\[-6pt] \nonumber\end{align}

(4.50) \begin{align} i k u'+ilv' + w'_z & = 0, \\[-6pt] \nonumber\end{align}

(4.51) \begin{align} \tau ' & = -i k({H} - z) {\eta } + {\eta } {s} ,\end{align}

subject to $u' = v' = w' = 0$ at $z= 0$ , as well as

(4.52) \begin{align} {\theta }_z = K{\theta } \sqrt {\sigma /\kappa } \quad \text{at }z = 0, \quad {\theta } = 0 \quad \text{at } z = H, \end{align}

while the mass conservation (4.10) becomes

(4.53) \begin{align} ik \int _0^H u' \textrm {d} z + il \int _0^H v' \textrm {d} z = 0. \end{align}

From (4.53) and (4.50), we can then conclude that

(4.54) \begin{align} \eta & = -\frac {sik}{l^2}\frac {\int _0^H (H-z) F_T(z) {\theta }(z) \textrm {d} z}{\int _0^H (H-z)F(z)\textrm {d} z},\end{align}

(4.55) \begin{align} w' = & - ik \int _0^z (z-z')F_T(z'){\theta }(z')\textrm {d} z' + i k \frac {\int _0^H (H-z') F_T(z') {\theta }(z) \textrm {d} z'}{\int _0^H (H-z')F(z')\textrm {d} z'} \int _0^z (z-z')F(z')\textrm {d} z'.\end{align}

As can be seen, substituting (4.55) into (4.47) together with the boundary conditions (4.52), yields an integro-differential eigenvalue problem for $\theta$ and $\sigma$ contains no dependence on $l$ .

This problem rationalises why the eigenvalue $\sigma$ becomes independent of $l$ for $l\gg 1$ in figure 10(d), emphasising how the growth rate does not vanish at large wavenumbers. Moreover, in the simultaneous limit of small $k$ , $ik U$ and $w'$ both vanish. In that case, the integro-differential eigenvalue problem above becomes identical to the simpler one in (4.11)–(4.13). In other words, the behaviour of regime B for short transverse wavelengths and long streamwise wavelengths ( $l\gg 1$ and $k\ll 1$ ) is actually given by the zero-wavenumber ( $k=l=0$ ) problem of § 4.1.2. In fact, setting $k = 0$ and looking for modes with $\eta = 0$ in the full stability problem (4.2)–(4.10) again leads to the same problem (4.11)–(4.13) regardless of the value of $l$ .

As shown in § 4.1.2, (4.11)–(4.13) admit an eigenvalue with positive real part on the second, intermediate branch but not the first and third branches when $k = l = 0$ . In fact, there is no solution to (4.11)–(4.13) on the first and third branches. As confirmed by our numerical results as well as those in Hindmarsh (Reference Hindmarsh2004), these results now carry over to the limit $k\ll 1$ regardless of the value of $l$ : regime B exists only on the second, intermediate branch, and corresponds to instability without any short-wavelength stabilisation in the transverse direction.

Figure 13.

Colour plots on the $(H,Q)$ -plane of (a) the maximum growth rate over all wavenumbers, and (b) the corresponding streamwise phase speed scaled by the mean speed of the base flow. The red dashed line marks the first fold point; the region above it corresponds to branch 2, and the region below it to branch 1.

4.4. Nonlinear evolution of transverse perturbations

Physically, the $k = \eta = 0$ case corresponds to adjacent columns of ice (at different values of $y$ ) decoupling from each other. We can exploit that insight to go beyond the linearised model we have looked at so far, and consider a simple version of the nonlinear problem (2.14)–(2.22) in order to better understand the implications of the missing short-wavelength cutoff. Suppose we restrict ourselves to initial conditions in which ice thickness $h$ is a constant and temperature is independent of the streamwise coordinate $x$ , as represented by the $k = 0$ mode with $\eta = 0$ in the linearised model. Then, $h$ is constant in time and space, with $T = T(y,z,t)$ and $v = w = 0$ at all times. Consequently, $T$ satisfies

(4.56) \begin{align} \begin{aligned} T_t - T_{zz} =\, & \bar {\tau }f(\bar {\tau },T) & \text{for } 0 \lt z \lt h, \\ T_t -\kappa T_{zz} =\, & 0 & \text{for } z \lt 0,\quad\kern7.5pt \end{aligned} \end{align}

subject to $T = 0$ at $z = h$ , $T_z \rightarrow -\varGamma /K$ as $z \rightarrow -\infty$ , and continuity of temperature and flux at $z = 0$ . Here $\bar {\tau }(z) = s(h-z)$ is the base solution for shear stress, which also remains constant in time.

The transverse coordinate $y$ does not appear explicitly anywhere in this geometrically reduced problem, either by way of parameters or through derivatives, and $T$ still evolves independently in each ice column. As a result, if the initial temperature varies with $y$ , then at each transverse location the solution may relax toward either the lower or the upper branch in figure 2(a), leading to a piecewise continuous steady state given continuous initial data. The solution presumably still depends continuously on the initial data for any finite time $t$ , and the problem does not appear to be ill posed as suggested by Hindmarsh (Reference Hindmarsh2004, Reference Hindmarsh2011).

The discontinuities involved in a piecewise continuous solution can be regularised by restoring the transverse diffusion term $T_{yy}$ . The resulting system

(4.57) \begin{align} \begin{aligned} T_t - \varepsilon ^2 T_{yy} &= T_{zz}+\bar {\tau }f(\bar {\tau },T) && \text{for } 0\lt z\lt h, \\ T_t - \kappa \varepsilon ^2 T_{yy} &= \kappa T_{zz} && \text{for } z\lt 0, \end{aligned} \end{align}

subject to the same boundary and interface conditions, is analogous to an Allen–Cahn-type system (Bartels Reference Bartels2015). Here, the local dynamics is governed not by a scalar reaction term, but by a differential operator in $z$ . Nevertheless, the analogy is appropriate because the local problem admits two stable steady states in the form of distinct vertical temperature profiles, corresponding precisely to the lower and upper branch base states. One may therefore expect the solution to separate into transverse regions that relax toward different branches, with sharp jumps replaced by narrow transition layers between them. In this sense, the dynamics is analogous to phase separation.

In addition to transverse diffusion of heat, lateral discontinuities in $T$ and therefore in $u$ will also be regularised by transverse shear stresses $\tau _{12}$ , which have also been omitted in our leading-order model, but affect the calculation of $\bar {\tau }(y,z) = \tau _{12}$ in (4.57) when retained.