2.1. Ice-shelf model

With $z=0$ corresponding to sea level, the base of the ice shelf is at $z=b(x,y,t)<0$ , and its thickness is $h(x,y,t)>0$ . Assuming hydrostatic equilibrium, the base elevation and the thickness are related by

where ${\it\rho}_{o}$ is the reference ocean density and ${\it\rho}_{i}$ the ice density. Mass balance and horizontal force balance are expressed as

(see MacAyeal Reference MacAyeal1989, for example), where $\boldsymbol{u}(x,y,t)=(u,v)$ is the vertically uniform ice velocity, ${\it\eta}$ the effective viscosity, $g$ the gravitational acceleration and $m$ is the basal melt rate (volume of water melted per unit area, per unit time). The gradient $\boldsymbol{{\rm\nabla}}=(\partial /\partial x,\partial /\partial y)$ is a two-dimensional operator. For the commonly used Glen’s law ice rheology, ${\it\eta}$ is a function of the strain rate, but for this paper we treat the ice as Newtonian so the viscosity is constant.

The ice depth and velocity are assumed known at the grounding line and, as mentioned previously, the ice depth is zero at the front, so boundary conditions are

where $h_{g}(y)$ and $u_{g}(y)$ may in general vary with the transverse coordinate $y$ . Aside from requiring the velocity to be bounded, boundary conditions on the velocity components at $x=X$ are not required due to the degeneracy of the stress components resulting from $h=0$ .

2.2. Melting

Melting is determined through the balance of heat transferred to and from the ice–ocean interface. Following previous authors (e.g. Jenkins, Nicholls & Corr Reference Jenkins, Nicholls and Corr2010), we parameterize the turbulent heat transfer as proportional to the temperature difference and velocity of the plume,

Here $L$ is the latent heat, $c$ is the specific heat capacity of the ocean water, ${\it\gamma}_{T}$ is a constant dimensionless heat transfer coefficient and $\boldsymbol{U}$ and $T$ are the velocity and temperature of the plume (i.e. outside an interfacial viscous and diffusive boundary layer). In (2.4) we have assumed that the ice-shelf temperature is already close to the melting temperature $T_{m}$ , and so ignore conduction into the ice; the effect of colder ice can be effectively included in (2.4) as a small correction to the latent heat (e.g. Jenkins Reference Jenkins2011; Wells & Worster Reference Wells and Worster2011).

If the salinity dependence of the melting point were included one would need, in addition to (2.4), to model the salinity at the interface; a number of parameterizations have been developed for this in the context of both ice shelves and sea ice (McPhee, Morison & Nilsen Reference McPhee, Morison and Nilsen2008; Jenkins et al. Reference Jenkins, Nicholls and Corr2010).

2.3. Plume model

For typical temperature and salinity variations near an ice shelf, the primary control on density is the salinity (with thermal expansion coefficient ${\it\beta}_{T}=3.87\times 10^{-5}~\text{K}^{-1}$ , ${\it\beta}_{S}=7.86\times 10^{-4}~\text{psu}^{-1}$ (Jenkins Reference Jenkins2011), and with typical temperature and salinity differences given in table 1, the percentage changes in density due to temperature and salinity changes are ${\sim}0.01\,\%$ and ${\sim}3\,\%$ , respectively). We therefore neglect thermal expansion and take the density difference between the ambient ocean and the plume as

Here ${\it\rho}_{a}$ and ${\it\rho}$ are the density of the ambient ocean and plume water, respectively, ${\it\rho}_{o}$ is the reference density, ${\it\beta}_{S}$ the haline contraction coefficient and $S_{{\it\Delta}}$ is the salinity deficit of the plume water from the ambient salinity $S_{a}$ .

Table 1.

Parameters, scales and non-dimensional variables used in this paper, roughly based on observation of the high-melt region downstream of the grounding line of Petermann Glacier. [J11] refers to Jenkins (Reference Jenkins2011) and [RS08] refers to Rignot & Steffen (Reference Rignot and Steffen2008). The length scale $x_{0}$ is chosen such that ${\it\gamma}=1$ in (2.9). Alternatively, the smaller length scale $x_{0}=1~\text{km}$ gives rise to the parameter values in parentheses. The thermal transfer coefficient ${\it\gamma}_{T}$ is chosen to be consistent with an observed melt rate of $18~\text{m}~\text{yr}^{-1}$ for the Petermann glacier (Rignot & Steffen Reference Rignot and Steffen2008), given (2.10) and the value $Q_{g0}$ of the subglacial melt flux.

Conservation of mass, momentum, salt and heat in the plume are expressed as

Here $e$ is the rate of entrainment of ambient water into the plume, ${\it\kappa}$ is a turbulent eddy diffusivity and eddy viscosity (assumed equal for simplicity), $C_{d}$ is a constant turbulent drag coefficient, $S_{i}$ ( ${\approx}0$ ) is the salinity of the ice and $T_{a}$ is the ambient ocean temperature. The formulation of turbulent stresses follows that of Holland et al. (Reference Holland, Feltham and Jenkins2007), Payne et al. (Reference Payne, Holland, Shepherd, Rutt, Jenkins and Joughin2007) and Gladish et al. (Reference Gladish, Holland, Holland and Price2012). One could potentially argue for alternative formulations of these terms, either with $D$ moved inside the final derivative or with mixed-derivative terms that echo the viscous stretching terms in (2.1c,d ). Under the approximations made below, such alternative formulations would make no difference to our results.

The plume volume flux increases downstream due to entrainment of ambient fluid and addition of meltwater according to (2.6a ). Similar balances describe conservation of salt in (2.6d ) and conservation of heat in (2.6e ) but with the addition of a thickness-weighted eddy diffusion. The turbulent heat flux from the plume to the ice has been eliminated using (2.4) giving rise to the final term in (2.6e ). The three acceleration terms on the right-hand side of (2.6b ) and (2.6c ) are due to the action of buoyancy forces, thickness-weighted eddy viscosity and interfacial drag, respectively. The buoyancy forces are determined by the slope of the lower plume interface $z=b-D$ (figure 1), so include both the slope of the ice base and the gradients of the plume thickness. The latter are small, but are included for their possible role in stabilizing transverse perturbations. Entrainment is parameterized as proportional to the magnitudes of velocity and basal slope, (e.g. Jenkins Reference Jenkins2011)

where $E_{0}$ is a dimensionless entrainment constant.

Turbulent diffusion in the $x$ direction will shortly be neglected and we therefore impose plume boundary conditions only at the grounding line $x=0$ (we make the implicit assumption that $U^{2}>Dg{\it\beta}_{S}S_{{\it\Delta}}$ so that the plume flow is supercritical in the language of shallow water theory, else a condition would be required at the downstream boundary, e.g. Stoker Reference Stoker1957). Conditions for fresh subglacial discharge are

where $Q_{g}(y)$ and $U_{g}(y)$ are the flux (in $\text{m}^{2}~\text{s}^{-1}$ ) and velocity of the outflow across the grounding line, respectively, which may in general vary with the transverse coordinate.

2.4. Non-dimensionalization

To reduce the model to a minimal set of parameters we proceed to non-dimensionalize the variables. Depth and velocity scales for the ice flow, $h_{0}$ and $u_{0}$ , are chosen based upon typical ice-shelf depths and speeds, and given a suitable horizontal length scale $x_{0}$ (discussed below), the time scale is taken as $t_{0}=x_{0}/u_{0}$ . The melt-rate scale $m_{0}$ is determined below from the plume dynamics, and the ice-shelf model then depends upon three dimensionless parameters (see § 2.5)

Here $r$ is the density ratio, ${\it\gamma}$ measures the strength of gravitationally-driven stretching, and ${\it\lambda}$ measures the strength of melting in the ice-shelf mass balance (2.1b ). Two natural choices of length scale $x_{0}$ would be such that ${\it\gamma}=1$ (with $x_{0}$ representing the length scale over which the ice thickness changes), or such that ${\it\lambda}=1$ (with $x_{0}$ representing the length scale of the ice shelf).

To estimate the size of terms in the plume, we take a typical size of subglacial discharge $Q_{g0}$ . This determines a driving buoyancy flux $Q_{g0}g{\it\beta}_{S}S_{a}$ , which gives rise to natural scales for the plume velocity and melt rate

The conservation equations further suggest appropriate scalings for the plume thickness $D$ , salinity deficit $S_{{\it\Delta}}$ and temperature deficit, which we define as $T_{{\it\Delta}}=T_{a}-T$ . These scalings are

Upon adopting these scalings, the plume model is found to depend on the following dimensionless parameters (see § 2.5 below):

Here ${\it\nu}$ is the dimensionless eddy diffusivity and ${\it\delta}$ represents the size of the plume-thickness correction to the buoyancy term; both are typically small and will henceforth be retained only when multiplying $y$ derivatives, to account for their possible smoothing effect on transverse perturbations. The parameters ${\it\epsilon}_{g}$ and ${\it\epsilon}_{m}$ represent the size of the subglacial discharge and subshelf meltwater flux, respectively, relative to the entrained flux. Both of these ratios tend to be small, but these fluxes make important contributions to the buoyancy and temperature. Finally, ${\it\mu}$ is the dimensionless drag coefficient and ${\it\beta}$ is the inverse Stefan number.

Representative values of these parameters are estimated in table 1 for two different choices of length scale. We base our calculations on the values obtained by taking $x_{0}$ to be the natural stretching length scale of the ice shelf, such that ${\it\gamma}=1$ in (2.9) and the one-dimensional ice shelf has a dimensionless slope of order unity near the grounding line. However, we anticipate that our analysis is most appropriate over shorter length scales (we are primarily interested in the development of channels just downstream of the grounding line) so we use parameter estimates on a shorter length scale $x_{0}=1~\text{km}$ to further simplify the model in § 2.6. The neglect of ambient temperature stratification and pressure dependence of the freezing point are likely to be good approximations at this shorter length scale.

2.5. Non-dimensional model

The non-dimensional model is summarized here, with all subsequent variables henceforth corresponding to dimensionless versions of their earlier dimensional counterparts. Equations for the ice shelf are

with

The reduced equations for the plume are

The last two equations for salinity and temperature deficits, $S_{{\it\Delta}}$ and $T_{{\it\Delta}}$ , respectively, are obtained by subtracting $S_{a}$ and $T_{a}$ times the mass equation (2.6a ) from (2.6d ) and (2.6e ) for $S$ and $T$ , and then non-dimensionalizing.

The non-dimensional boundary conditions are

and

The boundary values $h_{g}$ , $u_{g}$ , $Q_{g}$ and $U_{g}$ , which may depend on the transverse coordinate $y$ , have been scaled with $h_{0}$ , $u_{0}$ , $Q_{g0}$ and $U_{0}$ , respectively.

2.6. Simplification

The problem stated in § 2.5 could be solved numerically to determine the coupled evolution of ice shelf and plume. Indeed, the studies of Gladish et al. (Reference Gladish, Holland, Holland and Price2012) and Sergienko (Reference Sergienko2013) solved similar sets of equations with a number of additional terms. Since our purpose is to gain a more detailed understanding of the processes and parameter dependencies, we take a different approach that forces us to make some further simplifications at this point.

Specifically, we wish to study the channelization of an otherwise $y$ -independent state due to perturbations at the grounding line, and for this purpose it is highly advantageous to have an analytical one-dimensional solution. Such a solution can be found if one assumes that the buoyancy source is dominated by subglacial discharge rather than by melting from the shelf, that the entrained ambient water provides the plume with a large thermal inertia compared to heat transfer at the ice–ocean interface and if one ignores turbulent drag. Mathematically, this corresponds to taking the asymptotic limits

Although not necessary to obtain an analytic solution, we also take the limit ${\it\epsilon}_{g}\ll 1$ . Under these assumptions, the plume mass is dominated by the entrained flux, and the buoyancy flux is dominated by the initial grounding-line source ( ${\it\epsilon}_{m}\ll {\it\epsilon}_{g}\ll 1$ ); the subshelf meltwater does not have an appreciable effect on plume temperature ( ${\it\epsilon}_{m}\ll {\it\beta}$ ); and the drag does not affect the velocity ( ${\it\mu}\ll 1$ ).

As seen in the next section, the solution for the plume in this limit becomes remarkably simple. However, the neglect of some of these terms is hard to reconcile with the estimates in table 1, where ${\it\mu}$ is not small and ${\it\epsilon}_{m}$ is barely less than ${\it\epsilon}_{g}$ . Nevertheless, we justify studying this limit on three counts: firstly, the need to make some simplification in order to make headway; secondly, because the essential dynamics driving channelization are still encapsulated; and thirdly and most significantly, because this limit gives a good approximation to the full behaviour relatively close to the grounding line, where channelization is most likely to occur. This can be seen by comparing with solutions to the full system, or by estimating the dimensionless parameters based on the smaller length scale $x_{0}$ , shown in parentheses in table 1. The buoyancy source from subshelf melting and the effect of turbulent drag have an increasingly important effect along the length of the shelf, but in reality, other factors that we have already neglected also come into play on that larger scale, such as the ambient stratification, Coriolis forces and iceberg calving.

We adopt the limits (2.18) for the remainder of this study. Essentially, this means neglecting all terms in ${\it\epsilon}_{m}$ , ${\it\epsilon}_{g}$ , ${\it\beta}$ and ${\it\mu}$ in (2.14), (2.15), and (2.16), although some care is required to correctly account for the boundary conditions. Terms in ${\it\nu}$ and ${\it\delta}$ are retained for their possible stabilizing influence on perturbations.

4.1. Perturbed equations

We now reintroduce the transverse coordinate $y$ and study two-dimensional perturbations from the one-dimensional state found above. Assuming small amplitude perturbations and steady state in time, we derive a linear system of equations and boundary conditions for their amplitude as a function of along-shelf distance $x$ , and use this in § 5 to analyse the downstream development of transverse variability at the grounding line. This is not a temporal linear stability analysis; rather, we are interested in the growth or decay of perturbations as $x$ increases. In appendix B, we analyse an equivalent time-dependent problem, finding no evidence of temporal instability arising from global perturbations in the initial condition in the absence of boundary perturbations at the grounding line.

We look for steady solutions to (2.13) and (2.15) in the form

with a similar expansion for each of the variables. Here, $\overline{h}(x)$ denotes the steady solution from (3.8), $\tilde{h}(x)$ the $x$ dependence of the perturbation (to be solved for) and $k$ is the transverse wavenumber. Similar notation is used for the other variables. A perturbation with more general $y$ dependence could be written as a superposition of such Fourier mode solutions.

The equations are again simplified by treating buoyancy $B=DS_{{\it\Delta}}$ as a single variable. We also subtract the linearized version of the equation for mass (2.15a ) from that for momentum (2.15b ). From (2.13) and (2.15) the linearized corrections to the base state satisfy

and respectively (where primes denote a derivative with respect to $x$ ). In (4.2) we have used $\overline{U}=\overline{B}=1$ are constant, but have retained them for easier interpretation of the physical meaning of each term. We have also used $b=-r^{-1}h$ , changing basal elevation to total ice thickness. There is no need to include an equation for the temperature deficit correction as its effect on the melt rate (2.14) is small; the cause of variable heat transfer is due to variable plume velocity, rather than temperature.

4.2. Boundary conditions for perturbations

Boundary conditions for the linearized equations follow from the original conditions (2.16a – f ). Having made the approximation ${\it\epsilon}_{g}\rightarrow 0$ , however, the latter conditions on the plume variables are now replaced with the matching conditions (3.6a−c ). (Strictly speaking, the use of these matching conditions assumes a certain asymptotic ordering between the smallness of ${\it\epsilon}_{g}$ and the size of the linear perturbations; we do not concern ourselves with this detail.)

We allow for perturbations in the grounding-line ice thickness and subglacial discharge in the form $h_{g}=1+\tilde{h}_{g}\text{e}^{\text{i}ky}$ , $Q_{g}=1+\tilde{Q}_{g}\text{e}^{\text{i}ky}$ , while supposing there are no perturbations to the ice velocity there. Thus the necessary conditions are

It might be of concern that there are only six conditions here, whilst the linearized system (4.2) is of ninth order (i.e. (4.2b , c ) are both second order, while the other five equations are first order). This is again due to the singular nature of the problem: in (4.2), the highest derivatives of the ice velocity components are multiplied by the base-state ice thickness $\overline{h}$ , which is zero at the right-hand boundary $x=\overline{X}$ ; and similarly, the derivatives of the plume velocity components are multiplied by plume thickness $\overline{D}$ , which is zero at the left-hand boundary $x=0$ . As is common for such singular boundary value problems, effective conditions are provided by demanding that the system be satisfied and variables remain bounded at the boundary. In this case, these consistency conditions are

where we have used $\overline{u}^{\prime }={\it\gamma}\overline{h}$ , and $\overline{h}^{\prime }(0)=-({\it\lambda}+{\it\gamma})$ from the base solution (3.8).

Finally, the end condition (2.17) is linearized to find the perturbation $\tilde{X}$ in the length of the shelf. Inserting the expansion $X\sim \overline{X}+\tilde{X}\text{e}^{\text{i}ky}$ into the perturbed form for $h$ in (4.1), and expanding as a Taylor series, the condition $h=0$ at $x=X$ becomes

where $\tilde{h}(\overline{X})$ is determined from the solution to (4.2). Hence

Physically, this simply expresses the fact that where the ice thickness has been locally decreased (in a channel), the shelf would be shorter, and vice versa.

5.1. Perturbations in grounding line thickness: numerical results

In this section we look for steady-state perturbations of the form (4.1) forced by an ice-thickness perturbation $\tilde{h}_{g}$ at the grounding line, and we establish the effects of ice deformation, plume dynamics and transverse diffusion, for different transverse wavenumbers $k$ .

Numerical solution of the linearized system (4.2) is described in appendix A. Figure 4 shows an example solution for the perturbations in ice thickness $\tilde{h}$ , and plume velocity $\tilde{U}$ for transverse wavenumber $k=10$ . Also plotted in panel (e) are the thickness and plume velocity perturbations in the $(x,y)$ plane, over two wavelengths. This example solution exhibits the typical behaviour of the system (4.2). From the grounding line, perturbations in thickness $\tilde{h}$ initially decay due to ice deformation, but as transverse flow $\tilde{V}$ concentrates fresher, more buoyant water into the channels, the increased along-channel velocity $\tilde{U}$ increases heat transfer into the channels, thereby leading to an increase in the perturbation of ice thickness (i.e. enhanced ice thinning). This positive feedback loop is moderated by the effects of diffusion in the transverse direction.

Figure 4.

A typical numerical solution to the linearized problem (4.2) for the magnitude of perturbations to (a) the ice thickness $\tilde{h}$ , (b) ice velocity components $(\tilde{u} ,\tilde{v})$ , (c) plume velocity components $(\tilde{U} ,\tilde{V})$ and (d) buoyancy $\tilde{B}$ for transverse wavenumber $k=10$ , given an initial perturbation $\tilde{h}_{g}=1$ , $\tilde{B}_{g}=0$ . (e) Ice thickness perturbation in the $(x,y)$ plane, with the plume velocity perturbation overlaid as directed streamlines. The flow is concentrated into the channels (darker regions). Parameters are ${\it\lambda}=0.37$ , ${\it\gamma}=1$ , ${\it\nu}=0.02$ , ${\it\delta}=0.036$ .

To further understand the importance of each effect, we compute the solution of (4.2) over a range of wavenumbers $k$ and parameter values, comprising the following three cases:

1. (i) ice deformation in the absence of plume dynamics ((4.2a–c ) are solved alone, with $\tilde{U} =0$ );

2. (ii) coupled perturbation of ice and plume, with eddy diffusivity neglected ( ${\it\nu}=0$ );

3. (iii) coupled perturbation of ice and plume, with eddy diffusivity included.

In each case, we assume a unit perturbation in the ice thickness at the grounding line and use the value ${\it\lambda}=0.37$ . For the latter two cases, we also explore the effect of including or excluding the buoyancy correction term ${\it\delta}$ .

The results of the numerical calculations for each of these cases are displayed in figure 5. The dependence of the thickness perturbation $\tilde{h}$ on both downstream distance $x$ and transverse wavenumber $k$ is important, so we plot the $x$ dependence for a given wavelength ( $k=10$ ), and the $k$ dependence (that is, the power spectrum) for a given position $x$ . The midway point $x=\overline{X}/2\approx 1.4$ is chosen for this purpose. Using the values in table 1, this corresponds to a distance approximately $15~\text{km}$ downstream of the grounding line, comparable in scale to the high-melt region under the Petermann ice tongue (Rignot & Steffen Reference Rignot and Steffen2008).

Figure 5.

Numerical computation of the amplitude of a $k$ th mode perturbation, given a unit amplitude at the grounding line and fixed parameters ${\it\lambda}=0.37$ , ${\it\gamma}=1$ : (a) the evolution along the ice shelf, for a given wavenumber $k=10$ , for various values of buoyancy correction term ${\it\delta}$ and eddy diffusivity ${\it\nu}$ , as indicated by solid-dashed lines and symbols respectively; (b) the magnitude of perturbations at a given position $x=\overline{X}/2$ , showing unbounded growth as $k\rightarrow \infty$ for ${\it\nu}=0$ , while the inclusion of diffusion smooths out high wavenumber perturbations, leading to a selected wavenumber $k_{\mathit{max}}$ where the amplitude is largest (values of ${\it\delta}$ and ${\it\nu}$ are the same and indicated in the same manner as in (a)); (c) the numerical and WKB approximation (5.3), with ${\it\nu}={\it\delta}=0$ , grow exponentially in $k^{1/2}$ as $k\rightarrow \infty$ ; (d) numerically computed values of the wavenumber $k_{\mathit{max}}$ for positive ${\it\nu}$ , asymptotically fit a $-2/3$ power-law as ${\it\nu}\rightarrow 0$ , as shown on a logarithmic scale. The coefficient ${\it\alpha}=0.91$ is determined by fitting a power law to the two smallest values of ${\it\nu}$ .

The effect of ice deformation in isolation, depicted in figure 5, is to cause perturbations to decay downstream from their grounding-line value. This decay rate, however, is largely independent of wavenumber $k$ , and does not provide sufficient smoothing to prevent channelization when the plume dynamics are included. In that case, but without including eddy diffusion ( ${\it\nu}=0$ ), perturbations grow rapidly downstream, and the amplitude increases with increasing $k$ (we show below that the growth is exponential in $k^{1/2}$ as $k\rightarrow \infty$ ).

The inclusion of the eddy diffusivity ${\it\nu}=0.02$ has the effect of smearing short wavelength features, and thus limits the development of high wavenumber perturbations. The amplitude is maximised for a certain value $k_{\mathit{max}}$ , and decreases for larger $k$ (smaller wavelengths). As $k\rightarrow \infty$ , the mode amplitude tends toward that given purely by the ice deformation. For a smaller value ${\it\nu}=0.002$ , the preferentially selected wavenumber $k_{\mathit{max}}$ is larger, and the peak amplitude is an order of magnitude larger. Note that the value of $k_{\mathit{max}}$ depends slightly on the value of $x$ at which the amplitude is measured; we always measure it at $x=\overline{X}/2$ for consistency.

Our final numerical experiment regards the thus far neglected buoyancy correction term proportional to ${\it\delta}$ . While the inclusion of this term moderates the amplitude of perturbations, it does not prevent the unbounded growth for large $k$ . The numerical results for ${\it\delta}=0.036$ are included in figure 5 for each of the previous values of ${\it\nu}$ ; it does not qualitatively affect the behaviour.

5.2. Asymptotic results for large $k$

These numerical results may be understood in light of the asymptotic behaviour of the system (4.2) in the limit that wavenumber $k$ is large. From the two stress equations for the ice flow (4.2b , c ), the ice velocity components $\tilde{u}$ and $\tilde{v}$ must be of order $k^{-2}\tilde{h}$ and $k^{-1}\tilde{h}$ , respectively. This is consistent with the magnitudes seen in the numerical solution plotted in figure 4 for $k=10$ . Combining (4.2a ) and (4.2c ) results in

where we have used $\overline{u}^{\prime }={\it\gamma}\overline{h}$ from the base state solution (3.8). The first term on the left-hand side represents the advection of ice thickness by the background flow, while the second term represents decay of perturbations due to both transverse flow (the $\overline{h}$ term) and stretching (the $\overline{u}^{\prime }$ term). At leading order, these contributions balance plume-driven melting (the $-{\it\lambda}\tilde{U}$ term). For large $k$ , (5.1) may therefore be used in place of (4.2a–c ).

In the absence of plume dynamics, that is, the first case discussed in the previous subsection, $\tilde{U} =0$ and (5.1) may be integrated to give

Although not shown, this formula agrees with the numerical solution (marked ‘ice’) in figure 5 when $k$ is large, being within 10 % of the numerical solution when $k>6$ , and within 1 % when $k\gtrsim 60$ . As observed in the numerical solution, the perturbation decays as the base-state ice velocity $\overline{u}$ increases (see figure 3), but this decay rate is asymptotically independent of $k$ and so does not preferentially smooth out high wavenumber perturbations. Thus, at leading order, the viscous stretching damps all perturbations equally.

Next we include the plume dynamics with no eddy diffusivity or buoyancy correction term ( ${\it\nu}={\it\delta}=0$ ). The unbounded growth depicted in figure 5 is explained by again taking the large $k$ limit. We carry this out by assuming a WKB (Wentzel–Kramers–Brillouin) ansatz, in which each variable is supposed to have exponential dependence on some function of $x$ and a power of $k$ , so that differentiation increases a term’s order in $k$ . The determination of dominant balances in (4.2d–g ) shows that $\tilde{U}$ and $\tilde{V}$ must be of order $k^{1/2}\tilde{h}$ , while $\tilde{B}$ and $\tilde{D}$ are of order $k\tilde{h}$ . The exponent in the WKB expansion, and thus differentiation, is of order $k^{1/2}$ . Details of the full WKB analysis are included in appendix C. The process is similar to that routinely performed on second-order problems, although the application of boundary conditions is hampered by the the existence of an order $k^{-1}$ boundary layer near $x=0$ due to the singular nature of the system (4.2d–g ). The final result is that

where the functions $H$ and $C$ , and constant $A_{h}$ are given by

This result is only valid outside the boundary layer near $x=0$ , as described in appendix C. Most important to note is that the exponent is proportional to $k^{1/2}$ and, roughly, to $x^{1/2}$ , indicating that the perturbation grows exponentially as $x$ increases downstream (note that the integrand in $C(x)$ scales as $x^{-1/2}$ as $x\rightarrow 0$ , so the integral, while improper, is defined). In figure 5(c) we plot the numerical and asymptotic solutions for fixed $x=\overline{X}/2$ against $k^{1/2}$ , with the exponential dependence manifest as a straight line when the $\tilde{h}$ axis is scaled logarithmically. While the numerical and asymptotic solutions are distinguishable (the first neglected term in the expansion of $H(x)$ in (5.3b−d ) is order $k^{-1/2}$ ), it is clear that the exponential dependence on $k^{1/2}$ has been faithfully captured.

5.3. Channelizing mechanism

The asymptotic analysis presented above also provides a way to understand the positive feedback loop that causes the channelizing instability. By carrying out the WKB analysis described in appendix C, the dominant terms in each equation of the plume model (4.2d–g ) are found to be

In the third balance (5.4c ), the buoyancy force resulting from transverse thickness variability ( $\text{i}kr^{-1}\overline{B}\tilde{h}$ ) accelerates fluid in the transverse direction ( $\overline{D}\overline{U}\tilde{V}^{\prime }$ ), driving flow up into the channels. The fourth balance (5.4d ) then shows that the buoyancy carried by the increased transverse flow ( $\text{i}k\overline{B}\tilde{V}$ ) has the effect of focusing buoyancy as $x$ increases ( $\overline{U}\tilde{B}^{\prime }$ ). In the second balance (5.4b ), the increased buoyancy in channels ( $-r^{-1}\overline{h}^{\prime }\tilde{B}$ ) accelerates the fluid in the along-channel direction ( $\overline{D}\overline{U}\tilde{U} ^{\prime }$ ). The increased velocity then increases the melt rate, which enlarges the perturbations in thickness according to (5.1), and completes the positive feedback loop. The first balance (5.4a ) simply describes the impact of transverse flow on increasing or decreasing local plume thickness, which does not play a dominant role in the channelizing instability.

5.4. Influence of eddy diffusivity

We do not attempt a similar WKB analysis of the system (4.2) with ${\it\nu}\neq 0$ , which becomes rather involved. However, given the large $k$ scalings of the terms in (4.2e–g ) obtained above, it can be seen that the diffusion terms become of the same order as the dominant terms in each equation when $k^{2}{\it\nu}=O(k^{1/2})$ . Thus diffusion will dominate if $k^{3/2}{\it\nu}\gg 1$ , and we expect the maximum amplitude wavenumber to be related to the diffusivity by

for some constant ${\it\alpha}$ . We verify (5.5a ) using the numerical solutions and an optimization algorithm to find $k_{\mathit{max}}$ (see figure 5 d). The comparison also gives us a numerical approximation of the factor ${\it\alpha}$ (for the given ${\it\lambda}=0.37$ , ${\it\gamma}=1$ , and $x=\overline{X}/2$ ) of

For ${\it\nu}=0.02$ (corresponding to a dimensional eddy diffusivity ${\it\kappa}=100~\text{m}^{2}~\text{s}^{-1}$ ), this results in a selected wavenumber $k_{max}\approx 12$ . Returning to the dimensions used in table 1, this corresponds to a wavelength of approximately $6~\text{km}$ , close to the ${\sim}5~\text{km}$ spacing observed on Petermann Glacier. An order of magnitude decrease in the diffusivity ${\it\nu}=0.002$ ( ${\it\kappa}=10~\text{m}^{2}~\text{s}^{-1}$ ) leads to a selected wavenumber $k_{\mathit{max}}\approx 54$ , or a wavelength of approximately $1~\text{km}$ .

Figure 6.

Numerical computation of the amplitude of a $k$ th mode perturbation in ice thickness $\tilde{h}$ given a unit amplitude perturbation in the buoyancy $\tilde{B}_{g}=-1$ and no thickness perturbation $\tilde{h}_{g}=0$ at the grounding line: (a) amplitude versus distance $x$ for $k=10$ for various values of buoyancy correction term ${\it\delta}$ and eddy diffusivity ${\it\nu}$ , as indicated by solid-dashed lines and symbols respectively; (b) amplitude versus wavenumber $k$ for midpoint $x=\overline{X}/2$ (values of ${\it\delta}$ and ${\it\nu}$ are the same and indicated in the same manner as in (a)). Eddy diffusion with ${\it\nu}=0.02$ is sufficient to dominate even at relatively small wavenumbers, so no selected wavenumber $k_{\mathit{max}}$ arises in that case.

5.5. Perturbations in the buoyancy flux

One can also consider grounding-line perturbations in the subglacial discharge $\tilde{Q}_{g}$ or equivalently buoyancy $\tilde{B}_{g}$ . We perform the same numerical experiments as the previous section calculating the amplitude of a $k$ th mode perturbation at $x=\overline{X}/2$ , taking $\tilde{B}_{g}=-1$ (the sign is reversed so that the ice thickness perturbation has the same sign as the previous subsection for comparison). As before, we consider diffusion-free perturbations before adding in diffusion. The results are plotted in figure 6.

The perturbation in $\tilde{h}$ still increases exponentially in $k^{1/2}$ in the absence of diffusion, although the prefactor is an order of $k$ smaller than when the grounding-line thickness is perturbed directly. The addition of diffusivity ${\it\nu}=0.02$ as in table 1 is sufficient to suppress the increase with $k$ completely, while for intermediate values of ${\it\nu}$ a preferentially selected wavelength $k_{\mathit{max}}$ does appear.

The smaller size of the developing ice-thickness perturbation compared to the case of imposed thickness changes at the grounding line is revealed in the large $k$ WKB approximation. As seen in appendix C the asymptotic behaviour in this case is

where $H$ and $C$ are the same as in (5.3b−d ), and $A_{B}\approx 0.21{\it\lambda}^{1/8}({\it\lambda}+{\it\gamma})^{3/2}$ . The power of $k$ on the prefactor in (5.6) is one less than that in (5.3a ).

Again, the inclusion of the buoyancy correction term ${\it\delta}$ has only a minor effect on the overall picture, failing to remove the exponential growth of large wavenumber perturbations in the absence of diffusion, and having an appreciable effect on the diffusive system only when ${\it\nu}$ is much smaller than ${\it\delta}$ .