1. Introduction
Much of the ice loss from the Antarctic and Greenland ice sheets occurs through melting at the ice–ocean interface. Melting is enhanced by the development of buoyant plumes, which can entrain warm ocean water and facilitate heat transfer to the ice interface (MacAyeal Reference MacAyeal and Jacobs1985; Jenkins Reference Jenkins1991; Motyka et al. Reference Motyka, Hunter, Echelmeyer and Connor2003). Such plumes may form both at the near-vertical fronts of tidewater glaciers and beneath the gently sloping base of floating ice shelves. In the latter case, melting interacts with ice deformation to mould the shape of the ice shelf itself, and it is this interaction that forms the focus of our study.
Ice shelves form when ice flows into sufficiently deep water that it becomes neutrally buoyant and floats (e.g. Schoof & Hewitt Reference Schoof and Hewitt2013). The line separating grounded ice from floating ice is referred to as the ‘grounding line’, and this is also the place where subglacial meltwater is discharged to the ocean. Even a small quantity of fresh subglacial discharge can be influential in establishing a buoyant plume (Jenkins Reference Jenkins2011).
A number of ice shelves have been observed to have channel-like incisions at their base, oriented roughly in the direction of ice flow (Rignot & Steffen Reference Rignot and Steffen2008; Bindschadler, Vaughan & Vornberger Reference Bindschadler, Vaughan and Vornberger2011; Vaughan et al. Reference Vaughan, Corr, Bindschadler, Dutrieux, Gudmundsson, Jenkins, Newman, Vomberger and Wingham2012; LeBrocq et al. Reference LeBrocq, Ross, Griggs, Bingham, Corr, Ferraccioli, Jenkins, Jordan, Payne, Rippin and Siegert2013). Observationally inferred basal melt rates are enhanced within the channels compared to the neighbouring thicker ice (Rignot & Steffen Reference Rignot and Steffen2008; Dutrieux et al. Reference Dutrieux, Vaughan, Corr, Jenkins, Holland, Joughin and Fleming2013), and plumes containing a mixture of entrained warm deep waters and ice shelf melt have been observed to emerge from under the Pine Island ice shelf in the vicinity of larger channels (Mankoff et al. Reference Mankoff, Jacobs, Tulaczyk and Stammerjohn2012). The existence of channels has an important but unclear effect on the long-term evolution of an ice shelf; bending stresses resulting from uneven melting may lead to crevassing and eventual breakup (Vaughan et al. Reference Vaughan, Corr, Bindschadler, Dutrieux, Gudmundsson, Jenkins, Newman, Vomberger and Wingham2012), while the focussing of buoyant flow may actually lower the average melt rate, thus protecting it from ocean warming to some extent (Gladish et al. Reference Gladish, Holland, Holland and Price2012; Millgate et al. Reference Millgate, Holland, Jenkins and Johnson2013). It is of interest to understand how such channels are created. One hypothesis is that the convective plume on the underside of the ice is unstable to transverse perturbations which promote uneven melting. Another is that localized sources of subglacial outflow at the grounding line generate stronger plumes and promote enhanced melting directly downstream. Alternatively, the channels may reflect topography at the grounding line, whose imprint is advected into the ice shelf, or they may reflect patterns of stress in the ice shelf.
In this paper, we explore the potential for channels to form through the interaction of plume-driven melting and ice-shelf topography. The basic mechanism here is that regions of high plume velocity induce enhanced heat transfer and thus faster melting, which increases the local transverse slope of the ice shelf and leads to flow focussing. There is a partial analogy with hill-slope erosion (Smith & Bretherton Reference Smith and Bretherton1972). We seek to determine how a one-dimensional ice-shelf profile is affected by transverse perturbations in the conditions at the grounding line. Our goal is to gain a basic understanding of the mechanisms rather than to make quantitative predictions; thus we adopt simplified models for the ice shelf and the plume on its underside. The results may also provide insight into related channelization problems where the incision of channels is caused by gravitationally-driven focussing of a wall-bounded current.
A number of previous studies are relevant. Ice-shelf melting has been modelled extensively using depth-averaged ocean layer models (e.g. MacAyeal Reference MacAyeal and Jacobs1985; Jenkins Reference Jenkins1991; Holland, Feltham & Jenkins Reference Holland, Feltham and Jenkins2007), adapted from the classical theory of buoyant plumes (Morton, Taylor & Turner Reference Morton, Taylor and Turner1956; Ellison & Turner Reference Ellison and Turner1959). General circulation models have also been used to study the effect of patterned ice shelf topography or subglacial discharge on the melt rate (e.g. Millgate et al. Reference Millgate, Holland, Jenkins and Johnson2013; Sciascia et al. Reference Sciascia, Straneo, Cenedese and Heimbach2013). Until recently, these models have assumed a fixed ice-shelf geometry. However, two recent numerical studies have examined the genuinely coupled problem in which melting feeds back on the ice shelf (Gladish et al. Reference Gladish, Holland, Holland and Price2012; Sergienko Reference Sergienko2013). These studies involved simultaneous solution of a depth-averaged ice flow model together with the buoyant plume that controls melting. They found that channelization can occur if there are topographic perturbations at the grounding line, or if lateral confinement of the ice shelf gives rise to across-flow variations in ice thickness. In this paper we examine in detail the destabilizing and stabilizing processes that control this channelization, albeit in a necessarily simplified model.
The paper is structured as follows. In § 2 we describe the elements of the model and perform non-dimensionalization so as to obtain a minimal set of parameters and to facilitate the subsequent analysis. In § 3 we consider one-dimensional solutions for a steady ice shelf with plume beneath, finding that under certain simplifications, the melting rate is approximately uniform. In § 4, the model equations are perturbed from the steady state to allow small transverse variations, and in § 5 we solve for such perturbations driven by topographical or discharge perturbations at the grounding line. In § 6 we discuss the results along with possible extensions. Further details and results on temporal stability are included in the appendices A–C.
2. Model
We first outline the basic ingredients of the model, before stating the equations and then non-dimensionalizing them. The problem is posed in two horizontal dimensions, with
$x$
in the direction of ice flow and
$y$
transverse to the flow (see figure 1).
Figure 1. (a) A schematic of the steady one-dimensional ice shelf. Ice thickness
$h$
and velocity
$u$
evolve in the along-shelf coordinate (
$x$
) from their grounding-line values (denoted with subscripts
$g$
) due to the deformation within the ice, and subshelf melting. The melt rate
$m$
is controlled by a buoyant plume layer of thickness
$D$
and velocity
$U$
, initiated by a meltwater source
$Q_{g}$
at the grounding line. The dynamics of the plume layer and melt rate depend on the temperature
$T$
and salinity
$S$
of the plume, the mass of which increases with
$x$
as ocean water (with temperature
$T_{a}$
and salinity
$S_{a}$
) is entrained at rate
$e$
. (b) A perturbation in thickness in the transverse (
$y$
) direction will grow or decay from its grounding-line magnitude
$\tilde{h}_{g}$
to
$\tilde{h}(x)$
, depending on the effects of ice deformation and plume dynamics.
The ice shelf extends from the grounding line at
$x=0$
to the front at
$x=X(y,t)$
, and is modelled as an extensional viscous thin-film flow, characterized by its vertical thickness
$h(x,y,t)$
and depth-averaged horizontal velocity
$\boldsymbol{u}(x,y,t)=(u,v)$
. In the glaciological literature, this model is referred to as the ‘shallow shelf approximation’ (MacAyeal Reference MacAyeal1989). The front position must be determined as part of the problem; in the assumed absence of calving, it is simply the position at which melting causes the ice thickness to reduce to zero. Ice flow from the grounded ice sheet is assumed known, and surface accumulation is ignored.
The plume is assumed to follow the underside of the ice without detaching, and is modelled using a form of the turbulent plume models commonly employed to study ice–ocean interactions (MacAyeal Reference MacAyeal and Jacobs1985; Jenkins Reference Jenkins1991, Reference Jenkins2011; Holland et al.
Reference Holland, Feltham and Jenkins2007; Gladish et al.
Reference Gladish, Holland, Holland and Price2012; Sergienko Reference Sergienko2013). The plume is characterized by its vertical thickness
$D(x,y,t)$
beneath the ice base
$z=b(x,y,t)$
, its thickness-averaged horizontal velocity components
$\boldsymbol{U}(x,y,t)=(U,V)$
, temperature
$T(x,y,t)$
and salinity
$S(x,y,t)$
. The time dependence of these quantities derives from the slow evolution of the ice–ocean interface, with the plume dynamics themselves being treated as quasi-static on this time scale. The plume is assumed to be initiated at the grounding line by the subglacial discharge of buoyant water.
We make a standard Boussinesq approximation, taking the ocean water density as a constant reference value,
${\it\rho}_{o}$
, except when it contributes to the buoyancy of the plume (see (2.5) below). Coriolis forces are also ignored. To simplify the thermodynamics, we take the melting temperature
$T_{m}$
to be constant, ignoring both pressure and salinity dependence, and assume that the ambient ocean’s temperature and salinity are spatially uniform. These are significant simplifications of the real problem, as they ignore changes in effective buoyancy flux that may eventually lead to detrainment further out along the shelf as well as refreezing due to a change in melting temperature. However, we believe they do not detract from the essential dynamics of channelization which is our focus, and they allow for a more transparent analysis. The possibility of including such features in future work is discussed in § 6.
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
$$\begin{eqnarray}b=-({\it\rho}_{i}/{\it\rho}_{o})\,h,\end{eqnarray}$$
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
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial h}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(h\boldsymbol{u})=-({\it\rho}_{o}/{\it\rho}_{i})m,\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial x}\left[2{\it\eta}h\left(2\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right]+\frac{\partial }{\partial y}\left[{\it\eta}h\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right]-(1-{\it\rho}_{i}/{\it\rho}_{o}){\it\rho}_{i}gh\frac{\partial h}{\partial x}=0,\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial x}\left[{\it\eta}h\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right]+\frac{\partial }{\partial y}\left[2{\it\eta}h\left(\frac{\partial u}{\partial x}+2\frac{\partial v}{\partial y}\right)\right]-(1-{\it\rho}_{i}/{\it\rho}_{o}){\it\rho}_{i}gh\frac{\partial h}{\partial y}=0,\qquad & \displaystyle\end{eqnarray}$$
$\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
$$\begin{eqnarray}\displaystyle & \displaystyle h=h_{g},\quad \boldsymbol{u}=u_{g}\boldsymbol{e}_{x}\quad \text{at }x=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \displaystyle h=0\quad \text{at }x=X, & \displaystyle\end{eqnarray}$$
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,
$$\begin{eqnarray}mL=c{\it\gamma}_{T}|\boldsymbol{U}|(T-T_{m}).\end{eqnarray}$$
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
$$\begin{eqnarray}{\it\rho}_{a}-{\it\rho}={\it\rho}_{o}{\it\beta}_{S}S_{{\it\Delta}},\quad S_{{\it\Delta}}=S_{a}-S.\end{eqnarray}$$
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
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(D\boldsymbol{U})=e+m, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(D\boldsymbol{U}U)=Dg{\it\beta}_{S}S_{{\it\Delta}}\left(\frac{\partial b}{\partial x}-\frac{\partial D}{\partial x}\right)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\kappa}D\boldsymbol{{\rm\nabla}}U)-C_{d}|\boldsymbol{U}|U & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(D\boldsymbol{U}V)=Dg{\it\beta}_{S}S_{{\it\Delta}}\left(\frac{\partial b}{\partial y}-\frac{\partial D}{\partial y}\right)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\kappa}D\boldsymbol{{\rm\nabla}}V)-C_{d}|\boldsymbol{U}|V & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(D\boldsymbol{U}S)=eS_{a}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\kappa}D\boldsymbol{{\rm\nabla}}S)+mS_{i}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(D\boldsymbol{U}T)=eT_{a}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\kappa}D\boldsymbol{{\rm\nabla}}T)+mT_{m}-\frac{mL}{c}. & \displaystyle\end{eqnarray}$$
$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)
$$\begin{eqnarray}e=E_{0}|\boldsymbol{U}||\boldsymbol{{\rm\nabla}}b|,\end{eqnarray}$$
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
$$\begin{eqnarray}DU=Q_{g},\quad \boldsymbol{U}=U_{g}\boldsymbol{e}_{x},\quad S=0,\quad T=T_{m},\quad \text{at }x=0,\end{eqnarray}$$
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)
$$\begin{eqnarray}r=\frac{{\it\rho}_{o}}{{\it\rho}_{i}},\quad {\it\gamma}=\frac{(1-{\it\rho}_{i}/{\it\rho}_{o}){\it\rho}_{i}gh_{0}x_{0}}{8{\it\eta}u_{0}},\quad {\it\lambda}=\frac{{\it\rho}_{o}m_{0}x_{0}}{{\it\rho}_{i}h_{0}u_{0}}.\end{eqnarray}$$
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
$$\begin{eqnarray}U_{0}=\left(\frac{Q_{g0}g{\it\beta}_{S}S_{a}}{E_{0}}\right)^{1/3},\quad m_{0}=\frac{c{\it\gamma}_{T}U_{0}(T_{a}-T_{m})}{L}.\end{eqnarray}$$
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
$$\begin{eqnarray}D_{0}=E_{0}h_{0},\quad S_{{\it\Delta}0}=\frac{Q_{g0}S_{a}}{D_{0}U_{0}},\quad T_{{\it\Delta}0}=\frac{{\it\gamma}_{T}x_{0}}{D_{0}}(T_{a}-T_{m}).\end{eqnarray}$$
Upon adopting these scalings, the plume model is found to depend on the following dimensionless parameters (see § 2.5 below):
$$\begin{eqnarray}{\it\nu}=\frac{{\it\kappa}}{U_{0}x_{0}},\quad {\it\delta}=\frac{D_{0}}{h_{0}}(=E_{0}),\quad {\it\epsilon}_{g}=\frac{Q_{g0}}{D_{0}U_{0}},\end{eqnarray}$$
$$\begin{eqnarray}{\it\epsilon}_{m}=\frac{m_{0}x_{0}}{D_{0}U_{0}},\quad {\it\mu}=\frac{C_{d}x_{0}}{D_{0}},\quad {\it\beta}=\frac{c(T_{a}-T_{m})}{L}.\end{eqnarray}$$
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
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial h}{\partial t}+\frac{\partial }{\partial x}(hu)+\frac{\partial }{\partial y}(hv)=-{\it\lambda}m, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial x}\left[2h\left(2\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right]+\frac{\partial }{\partial y}\left[h\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right]-8{\it\gamma}h\frac{\partial h}{\partial x}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial x}\left[h\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right]+\frac{\partial }{\partial y}\left[2h\left(\frac{\partial u}{\partial x}+2\frac{\partial v}{\partial y}\right)\right]-8{\it\gamma}h\frac{\partial h}{\partial y}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}b=-h/r,\quad m=|\boldsymbol{U}|\left(1-\frac{{\it\epsilon}_{m}}{{\it\beta}}T_{{\it\Delta}}\right).\end{eqnarray}$$
The reduced equations for the plume are
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(D\boldsymbol{U})=|\boldsymbol{U}||\boldsymbol{{\rm\nabla}}b|+{\it\epsilon}_{m}m, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(D\boldsymbol{U}U)=DS_{{\it\Delta}}\frac{\partial b}{\partial x}+{\it\nu}\frac{\partial }{\partial y}\left(D\frac{\partial U}{\partial y}\right)-{\it\mu}|\boldsymbol{U}|U, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(D\boldsymbol{U}V)=DS_{{\it\Delta}}\left(\frac{\partial b}{\partial y}-{\it\delta}\frac{\partial D}{\partial y}\right)+{\it\nu}\frac{\partial }{\partial y}\left(D\frac{\partial V}{\partial y}\right)-{\it\mu}|\boldsymbol{U}|V, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(D\boldsymbol{U}S_{{\it\Delta}})=\frac{{\it\epsilon}_{m}}{{\it\epsilon}_{g}}m+{\it\nu}\frac{\partial }{\partial y}\left(D\frac{\partial S_{{\it\Delta}}}{\partial y}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(D\boldsymbol{U}T_{{\it\Delta}})={\it\beta}m+m+{\it\nu}\frac{\partial }{\partial y}\left(D\frac{\partial T_{{\it\Delta}}}{\partial y}\right). & \displaystyle\end{eqnarray}$$
$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
$$\begin{eqnarray}h=h_{g},\quad \boldsymbol{u}=u_{g}\boldsymbol{e}_{x},\quad \text{at }x=0,\end{eqnarray}$$
$$\begin{eqnarray}DU={\it\epsilon}_{g}Q_{g},\quad \boldsymbol{U}=U_{g}\boldsymbol{e}_{x},\quad S_{{\it\Delta}}=1/{\it\epsilon}_{g},\quad T_{{\it\Delta}}={\it\beta}/{\it\epsilon}_{m},\quad \text{at }x=0,\end{eqnarray}$$
and
$$\begin{eqnarray}h=0,\quad \text{at }x=X.\end{eqnarray}$$
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
$$\begin{eqnarray}{\it\epsilon}_{m}\ll {\it\epsilon}_{g},\quad {\it\epsilon}_{m}\ll {\it\beta}\ll 1,\quad {\it\mu}\ll 1.\end{eqnarray}$$
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.
3. One-dimensional steady state
Taking the limits (2.18), and assuming independence of the transverse coordinate
$y$
, the buoyancy forcing and the entrainment parameterization in (2.15) are both proportional to the slope
$\partial b/\partial x$
. This allows the plume equations to be solved as a function of depth
$z=b(x)$
, independent of the actual shape of the interface. Moreover, with
${\it\epsilon}_{m}=0$
, the equation for the temperature deficit (2.15e
) decouples from the other plume equations. To simplify the analysis below, we define the buoyancy
$B=DS_{{\it\Delta}}$
. Changing variables according to
$$\begin{eqnarray}\frac{\partial }{\partial x}=\frac{\partial b}{\partial x}\frac{\partial }{\partial z},\end{eqnarray}$$
the dimensionless plume equations (2.15), with
${\it\epsilon}_{m}={\it\mu}={\it\beta}=0$
, are
$$\begin{eqnarray}\frac{\partial }{\partial z}(DU)=U,\quad \frac{\partial }{\partial z}(DU^{2})=B,\quad \frac{\partial }{\partial z}(BU)=0,\end{eqnarray}$$
representing increase of the mass flux by entrainment, a momentum balance between changing inertia of the plume and the buoyancy force, and conservation of the along-slope buoyancy flux. The boundary conditions are
$$\begin{eqnarray}DU={\it\epsilon}_{g}Q_{g},\quad U=U_{g},\quad B=Q_{g}/U_{g},\quad \text{at }z=b(0).\end{eqnarray}$$
For
${\it\epsilon}_{g}=0$
, the equations (3.2) are satisfied by constant velocity and buoyancy, and increasing plume thickness, that is
$$\begin{eqnarray}D=z-b(0)=(h_{g}-h(x))/r,\quad U=Q_{g}^{1/3},\quad B=Q_{g}^{2/3},\end{eqnarray}$$
recalling that
$z=b(x)=-h(x)/r$
. This does not satisfy the boundary conditions on the velocity or the buoyancy. However, the
${\it\epsilon}_{g}\rightarrow 0$
limit is singular; there is a boundary layer with width of order
${\it\epsilon}_{g}$
in which these conditions are satisfied, and the solution converges to (3.4) outside this boundary layer. In fact, an implicit solution to the system (3.2) with conditions (3.3) can be derived by finding
$U$
in terms of
$Q=DU$
, then substituting into the equation for
$\partial Q/\partial z$
and integrating. One obtains
$B=Q_{g}/U$
, and
$$\begin{eqnarray}U=\frac{Q_{g}^{1/3}({\it\epsilon}_{g}^{3}U_{g}^{3}Q_{g}^{2}-{\it\epsilon}_{g}^{3}Q_{g}^{3}+Q^{3})^{1/3}}{Q},\quad z-b(0)=\int _{{\it\epsilon}_{g}}^{Q/Q_{g}}\frac{Q_{g}^{2/3}q}{({\it\epsilon}_{g}^{3}U_{g}^{3}/Q_{g}-{\it\epsilon}_{g}^{3}+q^{3})^{1/3}}\,\text{d}q.\end{eqnarray}$$
Some examples of this solution are shown in figure 2, for different
$U_{g}$
. As
${\it\epsilon}_{g}\rightarrow 0$
, the solutions converge to that given by (3.4). For the rest of the paper, we take the limit
${\it\epsilon}_{g}\rightarrow 0$
, and use the simple solution (3.4). As a result, we replace the boundary condition (3.3) with matching conditions to the boundary layer as
$x\rightarrow 0$
, and hence
$z\rightarrow b(0)$
. These are
$$\begin{eqnarray}D=0,\quad U=Q_{g}^{1/3},\quad B=Q_{g}^{2/3},\quad \text{at }x=0.\end{eqnarray}$$
Figure 2. The solution (3.5) for plume flux
$Q=DU$
and velocity
$U$
, versus depth
$z=b$
, for
${\it\epsilon}_{g}=0.05$
(larger than in table 1 so as to make the boundary layer visible) and
$Q_{g}=1$
. Regardless of the subglacial discharge velocity
$U_{g}$
, the velocity converges to
$U=1$
outside of an order
${\it\epsilon}_{g}$
boundary layer, and the flux is close to the linear solution in (3.4) (it does not converge exactly to that solution as there is an order
${\it\epsilon}_{g}$
correction which grows with
$z$
; for smaller
${\it\epsilon}_{g}$
the dashed lines are much closer to the solid line).
With these approximate plume dynamics, the dimensionless melt rate
$m=U=Q_{g}^{1/3}$
is constant and scales with the one-third power of the subglacial discharge (cf. Jenkins Reference Jenkins2011). It is interesting to note here that it is the subglacial flux
$Q_{g}$
that is important; the prescribed subglacial discharge velocity
$U_{g}$
is only important within the very narrow boundary layer near the grounding line. Strictly, the melt rate also exhibits this order
${\it\epsilon}_{g}$
wide boundary layer; since
$T_{{\it\Delta}}(0)={\it\beta}/{\it\epsilon}_{m}$
the melt rate (2.14) is zero at
$x=0$
, but (2.15e
) for
$T_{{\it\Delta}}$
implies that it decays at a rate proportional to
$1/Q$
and so
$1/{\it\epsilon}_{g}$
. Outside of the boundary layer,
$T_{{\it\Delta}}$
drops out of (2.14) for melt rate under the limit
${\it\epsilon}_{m}\ll {\it\beta}$
. The physical interpretation here is that the plume temperature is dominated by entrainment and hence lies close to the ambient temperature
$T_{a}$
. The melting is therefore driven by the constant difference
$T_{a}-T_{m}$
from the melting point.
Furnished with this constant melt rate
$m=Q_{g}^{1/3}$
, the equations for the steady one dimensional steady ice shelf (2.13) may also be solved exactly. Taking
$Q_{g}=1$
(as we are permitted in one dimension, given the non-dimensionalization), the mass and momentum equations for the ice flow reduce to
$$\begin{eqnarray}\frac{\partial }{\partial x}(hu)=-{\it\lambda},\quad \frac{\partial }{\partial x}\left(h\frac{\partial u}{\partial x}\right)-2{\it\gamma}h\frac{\partial h}{\partial x}=0.\end{eqnarray}$$
Again, by the choice of scales we may take
$h_{g}=u_{g}=1$
, so that
$h(0)=u(0)=1$
, and
$h(X)=0$
. These are readily integrated to find (e.g. MacAyeal Reference MacAyeal1989)
$$\begin{eqnarray}\displaystyle X=\frac{1}{{\it\lambda}},\quad u=(1+{\it\gamma}X-{\it\gamma}X(1-x/X)^{2})^{1/2},\quad h=\left(\frac{(1-x/X)^{2}}{1+{\it\gamma}X-{\it\gamma}X(1-x/X)^{2}}\right)^{1/2}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
which are plotted in figure 3. The non-dimensional length of the ice shelf is controlled by the parameter
${\it\lambda}$
which measures the importance of basal melting relative to stretching.
Figure 3. Explicit solutions (3.8) for the non-dimensional base-state surface
$s$
and basal elevation
$b=-h/r$
(ice thickness
$h=\overline{h}(x)=s-b$
) and ice velocity
$u=\overline{u}(x)$
, for dimensionless melt rates
${\it\lambda}=0.5,1,2$
, and
${\it\gamma}=1$
.
4. Linearized 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
$$\begin{eqnarray}h(x,y)=\overline{h}(x)+\tilde{h}(x)\text{e}^{\text{i}ky},\end{eqnarray}$$
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
$$\begin{eqnarray}\displaystyle & \displaystyle (\tilde{h}\overline{u}+\overline{h}\tilde{u} )^{\prime }+\text{i}k\overline{h}\tilde{v}=-{\it\lambda}\tilde{U} , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 2[\overline{h}(2\tilde{u} ^{\prime }+\text{i}k\tilde{v})+2\tilde{h}\overline{u}^{\prime }]^{\prime }+\text{i}k\overline{h}(\text{i}k\tilde{u} +\tilde{v}^{\prime })-8{\it\gamma}(\overline{h}\tilde{h})^{\prime }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle [\overline{h}(\text{i}k\tilde{u} +\tilde{v}^{\prime })]^{\prime }+2\text{i}k\overline{h}(\tilde{u} ^{\prime }+2k\text{i}\tilde{v})+2\text{i}k\tilde{h}\overline{u}^{\prime }-8{\it\gamma}\text{i}k\overline{h}\tilde{h}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{U}\tilde{D}^{\prime }+\overline{D}\tilde{U} ^{\prime }+\text{i}k\overline{D}\tilde{V}+r^{-1}\overline{U}\tilde{h}^{\prime }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{D}\overline{U}\tilde{U} ^{\prime }+(2\overline{D}^{\prime }\overline{U}+{\it\nu}k^{2}\overline{D})\tilde{U} +r^{-1}\overline{h}^{\prime }\tilde{B}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{D}\overline{U}\tilde{V}^{\prime }+(\overline{D}^{\prime }\overline{U}+{\it\nu}k^{2}\overline{D})\tilde{V}+\text{i}kr^{-1}\overline{B}\tilde{h}+{\it\delta}\text{i}k\overline{B}\tilde{D}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{B}\tilde{U} ^{\prime }+\overline{U}\tilde{B}^{\prime }+\text{i}k\overline{B}\tilde{V}+{\it\nu}k^{2}\tilde{B}-{\it\nu}k^{2}\overline{D}^{-1}\overline{B}\tilde{D}=0, & \displaystyle\end{eqnarray}$$
$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
$$\begin{eqnarray}\tilde{h}(0)=\tilde{h}_{g},\quad \tilde{u} (0)=0,\quad \tilde{v}(0)=0,\end{eqnarray}$$
$$\begin{eqnarray}\tilde{B}(0)=\tilde{B}_{g}\equiv {\textstyle \frac{2}{3}}\tilde{Q}_{g},\quad \tilde{D}(0)=0,\quad \tilde{U} (0)=U_{g}\equiv {\textstyle \frac{1}{3}}\tilde{Q}_{g}.\end{eqnarray}$$
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
$$\begin{eqnarray}2\tilde{u} ^{\prime }(\overline{X})+\text{i}k\tilde{v}(\overline{X})=2{\it\gamma}\tilde{h}(\overline{X}),\quad \text{i}k\tilde{u} (\overline{X})+\tilde{v}^{\prime }(\overline{X})=0,\quad \tilde{V}(0)=-\frac{\text{i}k\tilde{h}_{g}}{{\it\lambda}+{\it\gamma}},\end{eqnarray}$$
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
$$\begin{eqnarray}\tilde{h}(\overline{X})+\overline{h}^{\prime }(\overline{X})\tilde{X}=0,\end{eqnarray}$$
where
$\tilde{h}(\overline{X})$
is determined from the solution to (4.2). Hence
$$\begin{eqnarray}\tilde{X}=-\frac{\tilde{h}(\overline{X})}{\overline{h}^{\prime }(\overline{X})}.\end{eqnarray}$$
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. Boundary-driven instabilities
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:
-
(i) ice deformation in the absence of plume dynamics ((4.2a–c ) are solved alone, with
$\tilde{U} =0$
); -
(ii) coupled perturbation of ice and plume, with eddy diffusivity neglected (
${\it\nu}=0$
); -
(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
$$\begin{eqnarray}\overline{u}\tilde{h}^{\prime }+(2{\it\gamma}\overline{h}+{\textstyle \frac{1}{2}}\overline{u}^{\prime })\tilde{h}=\overline{u}\tilde{h}^{\prime }+{\textstyle \frac{5}{2}}\overline{u}^{\prime }\tilde{h}=-{\it\lambda}\tilde{U} +O(k^{-1}),\quad k\rightarrow \infty ,\end{eqnarray}$$
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
$$\begin{eqnarray}\tilde{h}(x)\sim \tilde{h}_{g}\overline{u}(x)^{-5/2},\quad k\rightarrow \infty .\end{eqnarray}$$
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
$$\begin{eqnarray}\tilde{h}(x)\sim \tilde{h}_{g}A_{h}k^{-3/4}H(x)\exp (k^{1/2}C(x)),\quad k\rightarrow \infty ,\end{eqnarray}$$
$H$
and
$C$
, and constant
$A_{h}$
are given by 
$$\begin{eqnarray}C(x)=\int _{0}^{x}\left(\frac{-{\it\lambda}\overline{h}^{\prime }}{\overline{u}(1-\overline{h})^{2}}\right)^{1/4}\,\text{d}x,\quad H(x)=\frac{(-\overline{h}^{\prime })^{1/8}}{\overline{u}(1-\overline{h})^{3/4}},\quad A_{h}\approx 0.28{\it\lambda}^{-3/8}({\it\gamma}+{\it\lambda}).\end{eqnarray}$$
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
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{U}\tilde{D}^{\prime }\sim -\text{i}k\overline{D}\tilde{V} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{D}\overline{U}\tilde{U} ^{\prime }\sim -r^{-1}\overline{h}^{\prime }\tilde{B} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{D}\overline{U}\tilde{V}^{\prime }\sim -\text{i}kr^{-1}\overline{B}\tilde{h} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{U}\tilde{B}^{\prime }\sim -\text{i}k\overline{B}\tilde{V}. & \displaystyle\end{eqnarray}$$
$\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
$$\begin{eqnarray}k_{\mathit{max}}\sim {\it\alpha}{\it\nu}^{-2/3},\quad {\it\nu}\rightarrow 0,\end{eqnarray}$$
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
$$\begin{eqnarray}{\it\alpha}\approx 0.91.\end{eqnarray}$$
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
$$\begin{eqnarray}\tilde{h}(x)\sim -\tilde{B}_{g}A_{B}k^{-7/4}H(x)\exp (k^{1/2}C(x)),\quad k\rightarrow \infty ,\end{eqnarray}$$
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}$
.
6. Discussion and conclusion
By employing an idealized model, we have analysed the spatial instability of an under ice-shelf plume that leads to the growth of channels in the ice-shelf base from perturbations in grounding-line conditions. The numerical and asymptotic results derived in this paper help to understand the behaviour of more comprehensive models.
The chief mechanism of the spatial instability is the concentration of fresher, buoyant water due to transverse flow into pre-existing channels, which increases along-flow velocity and thus turbulent heat transfer. A similar instability, though of lower order in transverse wavenumber, is observed due to growth of perturbations in the along-shelf slope itself, which becomes important when there are lateral variations in subglacial discharge at the grounding line, rather than ice thickness.
Transverse ice motion driven by the ice-thickness gradient has a slight smoothing effect to counteract channel growth, but this effect does not preferentially smooth short-wavelength perturbations. Stabilisation of short wavelengths is due to the turbulent eddy diffusivity within the plume layer, and the transverse wavelength of maximum growth has a two-thirds power law dependence on the diffusivity. For physical parameters in table 1 and channelization caused by variations in ice thickness at the grounding line, this leads to a preferred wavelength of
${\sim}1$
–
$6~\text{km}$
for eddy diffusivity
${\it\kappa}=10$
–
$100~\text{m}^{2}~\text{s}^{-1}$
, compared to the
${\sim}5~\text{km}$
spacing observed on Petermann Glacier. Despite the simplified nature of our analysis, the fact that these values are comparable lends some support to this mechanism of wavelength selection. If the diffusivity is too large, then channelization does not occur, as we observed in the case of variable subglacial discharge and
${\it\kappa}=100~\text{m}^{2}~\text{s}^{-1}$
. Of course, even with the simplified physics in this model, the actual pattern that should be seen depends on the relative amplitudes of each wavelength of the initial grounding-line perturbation, and on the nonlinear interaction between this spectrum of wavelengths.
We have concentrated on the development of steady-state channels, seeded by transverse variability in the boundary conditions at the grounding line. Another issue is the possibility of time-dependent channelization due to instability from an initial perturbation, which could potentially occur even if the grounding-line thickness and buoyancy flux are uniform in
$y$
. Sergienko (Reference Sergienko2013) has shown that such ‘self-channelization’ is possible if the base-state thickness varies with
$y$
due to lateral confinement, although numerical simulations found no evidence that a one-dimensional state is unstable. In appendix B we have adapted the linearised perturbations from § 4 to include a temporal growth rate
${\it\sigma}$
, and analyse the temporal stability by determining the growth rates from the corresponding eigenvalue problem. The eigenvalues all appear to have negative real part, suggesting that the one-dimensional state is indeed stable to such perturbations. On the other hand, this stability is likely to be of limited relevance to the real world, where transverse perturbations at the grounding line are presumably rife.
It has been suggested that the shape of the ice–ocean interface affects the average melt rate across the shelf (e.g. Gladish et al. Reference Gladish, Holland, Holland and Price2012). Given that our study deals only with linear perturbations, the melting rate always averages out across the width, so it unfortunately does not allow us to say anything about the net effect of channelization on melting.
We now turn our attention to the limitations of this study, and to possible extensions. Some of the terms that we have neglected (such as thermal expansion, and increase of plume volume flux from melting) can be shown formally to be small. Other simplifications pose more serious restrictions, such as the neglect of turbulent drag and changes in plume buoyancy resulting from melting. We argue that, at least over small spatial length scales from the grounding line, the reduced model includes the dominant dynamics relevant to channel growth and decay. Of course, the assumption of small perturbations that is necessary to allow linearization is not valid far from the grounding line if the growth rate is large, nor if the initial perturbations are large. Thus, our analysis should be seen as examining the tendency for channelizing or smoothing of the ice shelf base rather than predicting its eventual evolution.
Arguably the most serious omissions are the changes in buoyancy along the shelf. These may be due both to the neglected buoyancy source from subshelf melting (the
${\it\epsilon}_{m}/{\it\epsilon}_{g}$
term in (2.15d
)) which increases buoyancy, and due to the ambient salinity stratification, which generally decreases it. The major cost of including these effects in the linearized model such as (4.2) is that the base state no longer has an explicit representation (as in (3.4), (3.8)), and must itself be computed numerically.
Another simplification is our use of Newtonian rheology for the ice flow. The same analysis as performed in this paper could be carried out for a power-law rheology; there is still an exact solution for the one-dimensional base state (MacAyeal Reference MacAyeal1989), and the equivalent linearized system to (4.2) could be derived with a variable effective viscosity. We expect the results to be very similar to those found here, since the perturbations in ice deformation velocity
$(\tilde{u} ,\tilde{v})$
have little effect on short wavelength perturbations, compared to the destabilising effects of the plume dynamics.
The simple parameterizations of turbulent entrainment (2.7) and heat transfer (2.4) are based on those of Jenkins (Reference Jenkins2011). These in turn are approximations of parameterizations with more complicated dependence on the Richardson number (for entrainment) and Reynolds number (for heat transfer) (Holland et al.
Reference Holland, Feltham and Jenkins2007; Payne et al.
Reference Payne, Holland, Shepherd, Rutt, Jenkins and Joughin2007), which were used in the numerical studies of Gladish et al. (Reference Gladish, Holland, Holland and Price2012) and Sergienko (Reference Sergienko2013). With such additional complexities, the base state would again have to be computed numerically. However, the basic form of these parameterizations is by no means certain. With regard to the heat transfer, it is the dependence of melt rate
$m$
on plume velocity
$U$
that is instrumental in producing the channelization outlined in this paper. As long as the melt rate is an increasing function of velocity, the linearization around the base state will produce a melt rate perturbation proportional to the velocity perturbation
$\tilde{U}$
, resulting in a similar expression to (4.2a
) but with modified prefactors. For example, if the salinity dependence of the melting temperature were included, together with a parameterization of the interfacial salinity using the three equation formulation of Jenkins (Reference Jenkins2011), one can show that a modification of the melting rate occurs, but this should not significantly change any of our conclusions.
The neglect of Coriolis forces was necessary to avoid an essentially two-dimensional base state for the ice thickness. Numerical evidence suggests that Coriolis forcing will tend to deflect the path of channels on longer spatial scales, and produce channels with asymmetric cross-sectional profiles, but not quantitatively change the spacing of channels (Sergienko Reference Sergienko2013). For this reason we believe it will not have a serious impact on the basic channel-forming mechanism outlined here, but we acknowledge that Coriolis forces are an inherent aspect of the real system.
Acknowledgements
This publication arises from research funded in part by Award No KUK-C1-013-04, made by King Abdullah University of Science and Technology (KAUST), and in part by the John Fell Oxford University Press (OUP) Research Fund. M.C.D. acknowledges further support from the Engineering and Physical Sciences Research Council (EPSRC) of the UK through grant no. EP/K008595/1. A.J.W. acknowledges financial support through the research program of the European Union FP7 award PCIG13-GA-2013-618610 SEAICE-CFD. I.J.H. was supported by a Marie Curie FP7 Integration Grant within the 7th European Union Framework Programme. The authors are indebted to the three anonymous referees whose comments significantly improved the final paper.
Appendix A. Numerical solution method
The system of equations (4.2) for perturbations to the base thickness and velocity is a linear, but non-constant coefficient, boundary-value problem. Because
$\overline{D}$
vanishes at
$x=0$
and
$\overline{h}$
vanishes at
$x=\overline{X}$
, the order of the system is reduced at these points and the boundary-value problem is therefore singular. Given the form of the base state (3.4), (3.8), an analytic or series solution of even the linearized problem is intractable. In order to compute numerical solutions to (4.2), we adopt a Chebyshev spectral collocation method. The seven unknown variables are represented by their values at
$N$
discrete Chebyshev nodes (mapped to the interval
$[0,\overline{X}]$
), and concatenated into a single solution vector
$$\begin{eqnarray}\boldsymbol{y}=[\tilde{h}~\tilde{u} ~\text{i}\tilde{v}~r\tilde{D}~\tilde{U} ~\text{i}\tilde{V}~\tilde{B}]^{\text{T}}\in \mathbb{R}^{7N-4}.\end{eqnarray}$$
Four variables,
$\tilde{h}(0)$
,
$\tilde{u} (0)$
,
$\tilde{v}(0)$
and
$\tilde{B}(0)$
, are prescribed from the four non-singular left-hand boundary conditions, namely (4.3) and the first of (4.4), so are removed from the solution vector
$\boldsymbol{y}$
; the remaining singular conditions in (4.4) and (4.5) are enforced naturally from solving the remaining components of the linear system at the boundary points. In addition, the problem has been simplified by suitable definition of the dependent variables that appear in (A 1). Firstly, solving for
$r\tilde{D}$
instead of
$\tilde{D}$
, along with the definition
$r\overline{D}=1-\overline{h}$
from (3.4), removes the parameter
$r$
explicitly from (4.2). Secondly, including the imaginary constant
$\text{i}$
with the transverse ice and plume velocities
$\tilde{v}$
and
$\tilde{V}$
results in a linear system (and therefore solution vector
$\boldsymbol{y}$
) with purely real elements (note that all imaginary constants
$\text{i}$
in (4.2) cancel upon making this substitution of dependent variables).
For given values of
${\it\lambda},{\it\gamma},k,{\it\nu}$
and
${\it\delta}$
, the system (4.2) may then be represented by a single matrix equation
$$\begin{eqnarray}\unicode[STIX]{x1D648}\boldsymbol{y}=\boldsymbol{b},\quad \unicode[STIX]{x1D648}\in \mathbb{R}^{(7N-4)\times (7N-4)},\;\boldsymbol{b}\in \mathbb{R}^{7N-4}.\end{eqnarray}$$
The derivatives in (4.2) are represented in
$\unicode[STIX]{x1D648}$
by the (dense) Chebyshev differentiation matrix (Trefethen Reference Trefethen2000), while the rows and columns corresponding to the non-singular boundary conditions are removed and used to construct
$\boldsymbol{b}$
. Typically we use
$N=50$
–
$100$
nodes to compute the results contained in this paper; this runs in a negligible amount of time using MATLAB on a personal computer, which allows us to explore the system over a wide range of parameters.
Appendix B. Stability with respect to initial perturbations
We consider the temporal linear stability of our system to time-dependent perturbations by adding a growth rate
${\it\sigma}$
to the perturbations (4.1), as
$$\begin{eqnarray}h(x,y,t)=\overline{h}(x)+\tilde{h}(x)\text{e}^{{\it\sigma}t+\text{i}ky},\end{eqnarray}$$
and similarly for the other variables. The linearized system of equations is then identical to that in (4.2), except for the first equation (in which the only time derivative appears), which becomes
$$\begin{eqnarray}{\it\sigma}\tilde{h}+(\tilde{h}\overline{u}+\overline{h}\tilde{u} )^{\prime }+k\overline{h}\text{i}\tilde{v}=-{\it\lambda}\tilde{U} .\end{eqnarray}$$
To determine the linear stability of this modified system we compute eigenvalues
${\it\sigma}$
and eigenfunctions
$\tilde{h}$
(and other variables) using a variant of the numerical code outlined in appendix A. We set all boundary conditions to be homogeneous (i.e. no grounding-line perturbations), and impose further the right-hand condition
$\tilde{h}(\overline{X})=1$
to remove the trivial solution. The eigenvalues
${\it\sigma}$
are found by solving the problem without the left-hand condition on
$\tilde{h}(0)$
, and then using a root-finding algorithm to obtain the values of
${\it\sigma}$
such that
$\tilde{h}(0)=0$
. For simplicity we only report results here for the case
${\it\nu}={\it\delta}=0$
, although we found that their inclusion does not qualitatively affect the stability.
Figure 7. (a) The first six eigenvalues
${\it\sigma}$
in the complex plane, for varying
$k$
.
$\text{Re}({\it\sigma})$
appear to tend to
$-\infty$
as
$k\rightarrow 0$
and
$k\rightarrow \infty$
, with maximum around
$k=1$
. (b) An expanded view of the path of the first eigenvalue as
$k$
varies. (c) The first six eigenfunctions
$\text{Re}(\tilde{h})$
for
$k=1$
. Note the increasing frequency of oscillations as
$\text{Im}({\it\sigma})$
increases.
The eigenvalues
${\it\sigma}$
thus found are shown in figure 7. The eigenvalues are complex with negative real part, with corresponding oscillatory eigenfunctions that grow exponentially in
$x$
. This downstream growth leads to a great deal of sensitivity in detecting the eigenvalues; indeed, we found it difficult to resolve the values for
$k\gtrsim 6$
. However, the real parts
$\text{Re}({\it\sigma})$
thus computed are all negative, with a maximum occurring at
$k\approx 1$
for each eigenvalue;
$\text{Re}({\it\sigma})$
tends to
$-\infty$
both as
$k\rightarrow 0$
and as
$k\rightarrow \infty$
. This strongly suggests that all eigenvalues have negative real part for all
$k$
and the system is stable.
The structure of complex eigenvalues and oscillatory eigenfunctions depicted in figure 7 is primarily due to the advection of the underlying ice flow. Indeed, this appears to be the primary reason for the lack of an inherent global instability here; there is no mechanism to cause information about perturbations to propagate backwards and interact with the forwards advection of the interface (cf. Gladish et al. Reference Gladish, Holland, Holland and Price2012).
Appendix C. WKB analysis of channel growth
The large
$k$
behaviour of the system (4.2) is determined by using a WKB-type ansatz, that is, expansions of the form
$$\begin{eqnarray}\tilde{h}(x)\sim Ak^{a}(H_{0}(x)+k^{-1/2}H_{1}(x))\exp (k^{1/2}C(x)),\quad k\rightarrow \infty ,\end{eqnarray}$$
and similar for other variables (an examination of the dominant balance indicates that
$\tilde{U}$
and
$\tilde{V}$
are of order
$k^{1/2}\tilde{h}$
, while
$\tilde{B}$
and
$\tilde{D}$
are of order
$k\tilde{h}$
). As with the standard application of WKB analysis (e.g. Holmes Reference Holmes1995), the exponent function
$C$
is determined from the leading order problem in large
$k$
, while the first coefficient function
$H_{0}$
(and those of the other variables) comes from the equations at next order. After considerable algebra, we find that
$$\begin{eqnarray}C^{\prime }(x)=\left(\frac{-{\it\lambda}\overline{h}^{\prime }}{\overline{u}(1-\overline{h})^{2}}\right)^{1/4},\quad H_{0}(x)=\frac{(-\overline{h}^{\prime })^{1/8}}{\overline{u}(1-\overline{h})^{3/4}}.\end{eqnarray}$$
The prefactor constant
$A$
and exponent
$a$
must be determined from the boundary conditions. Note for large
$k$
, (5.1) implies
$$\begin{eqnarray}\tilde{h}(0)=\tilde{h}_{g},\quad \tilde{h}^{\prime }(0)=-\frac{{\it\lambda}}{2}\tilde{B}_{g}-\frac{5}{2}{\it\gamma}\tilde{h}_{g}+O(k^{-1})\end{eqnarray}$$
(recall that
$\overline{u}(0)=1$
,
$\overline{u}^{\prime }(0)={\it\gamma}$
, and that
$\tilde{U} (0)=\tilde{B}_{g}/2$
from (4.4)). However, these conditions cannot be applied directly to (C 1) as the system (4.2) possesses an order
$k^{-1}$
boundary layer at
$x=0$
due to its singular nature there. From the base state solution (3.8) we have
$\overline{u}(0)=1$
,
$\overline{h}^{\prime }(0)=-({\it\gamma}+{\it\lambda})$
,
$\overline{D}(0)=0$
and
$\overline{D}^{\prime }(0)=({\it\gamma}+{\it\lambda})/r$
. Defining the inner variable
${\it\xi}$
by
$$\begin{eqnarray}x=\sqrt{\frac{{\it\lambda}+{\it\gamma}}{{\it\lambda}}}k^{-1}{\it\xi},\end{eqnarray}$$
the leading order plume equations (4.2e–g
) and large
$k$
ice flow equation (5.1) can be rearranged into a single fourth-order differential equation
$$\begin{eqnarray}{\it\xi}^{2}\tilde{h}_{{\it\xi}{\it\xi}{\it\xi}{\it\xi}}+6{\it\xi}\tilde{h}_{{\it\xi}}+4\tilde{h}_{{\it\xi}{\it\xi}}-\tilde{h}=0,\end{eqnarray}$$
with boundary conditions from (C 3)
$$\begin{eqnarray}\tilde{h}(0)=\tilde{h}_{g},\quad \tilde{h}_{{\it\xi}}(0)=-k^{-1}\sqrt{\frac{{\it\lambda}+{\it\gamma}}{{\it\lambda}}}\left(\frac{{\it\lambda}}{2}\tilde{B}_{g}+\frac{5}{2}{\it\gamma}\tilde{h}_{g}\right).\end{eqnarray}$$
The perturbation in buoyancy
$\tilde{B}_{g}$
appears only in the
$O(k^{-1})$
derivative condition; thus, the channel growth is dominated by the perturbations in ice thickness (assuming
$\tilde{h}_{g}\neq 0$
). For
$\tilde{h}(0)=1$
and
$\tilde{h}_{{\it\xi}}(0)=0$
, the far field behaviour of (C 5) is
$$\begin{eqnarray}\tilde{h}\sim 0.28{\it\xi}^{-3/4}\text{e}^{2\sqrt{{\it\xi}}},\quad {\it\xi}\rightarrow \infty .\end{eqnarray}$$
On the other hand, if
$\tilde{h}(0)=0$
and
$\tilde{h}_{{\it\xi}}(0)=1$
, the far field behaviour of (C 5) is
$$\begin{eqnarray}\tilde{h}\sim 0.42{\it\xi}^{-3/4}\text{e}^{2\sqrt{{\it\xi}}},\quad {\it\xi}\rightarrow \infty .\end{eqnarray}$$
Here the factors of
$0.28$
and
$0.42$
have been determined numerically by comparing the asymptotic behaviour to the exact series solutions of (C 5).
Transforming
${\it\xi}\mapsto x$
and matching to the behaviour of the outer problem as
$x\rightarrow 0$
gives the coefficients
$A$
and exponent of
$k$
for each case. For
$\tilde{h}_{g}\neq 0$
, the leading-order behaviour in large
$k$
is given by (C 6), thus
$$\begin{eqnarray}A\approx 0.28\tilde{h}_{g}{\it\lambda}^{-3/8}({\it\gamma}+{\it\lambda}),\quad a=-3/4.\end{eqnarray}$$
For
$\tilde{h}_{g}=0$
, the behaviour follows from (C 7), thus
$$\begin{eqnarray}A\approx -0.21\tilde{B}_{g}{\it\lambda}^{1/8}({\it\gamma}+{\it\lambda})^{3/2},\quad a=-7/4.\end{eqnarray}$$
These are the two results given in (5.3) and (5.6), respectively.
