1. Introduction
To better evaluate sea level rise (Vaughan & Arthern Reference Vaughan and Arthern2007) in a climate change regime, it is crucial to develop ice-flow models that are mechanically complete while being computationally tractable. Ice is known to behave like a slow non-Newtonian fluid, governed by Glen’s flow law (Glen Reference Glen1953). Thus, the velocity and the pressure of the ice satisfy the nonlinear Stokes equations. In practice, solving these equations requires considerable computational resources and complex meshing procedures to be implemented at the large scale (in time and space), which is necessary when modelling ice sheets. As a consequence, a number of simplified models based on different types of scaling have been proposed in the last few decades to improve simulations of real glaciers. Simplifications in the past have often been based on the small aspect ratio (between the characteristic height and length) of glaciers. More precisely, after writing the solution of the Stokes equations as an asymptotic expansion in the aspect ratio
${\it\epsilon}$
, high-order terms in
${\it\epsilon}$
are neglected. More recently, a complementary scaling based on the ratio
${\it\lambda}$
of vertical to horizontal stress proved to be further relevant to derive efficient models (Schoof & Hindmarsh Reference Schoof and Hindmarsh2010).
Two categories of models derive from the scaling based on the aspect ratio
${\it\epsilon}$
. In the first category, each order that was not neglected gives rise to a system of equations: the shallow ice approximation (SIA) for
${\it\epsilon}^{0}$
(Fowler & Larson Reference Fowler and Larson1978; Hutter Reference Hutter1983), the first-order SIA for
${\it\epsilon}^{1}$
and the second-order SIA for
${\it\epsilon}^{2}$
(Baral, Hutter & Greve Reference Baral, Hutter and Greve2001; Egholm et al.
Reference Egholm, Knudsen, Clark and Lesemann2011; Ahlkrona, Kirchner & Lötstedt Reference Ahlkrona, Kirchner and Lötstedt2013), which can be solved iteratively at much lower costs than the original Stokes problem. In contrast, the models of the second category exploit the loss of higher-order terms by eliminating unknown variables and by reducing the dimension of the mathematical model. In this category, the first level of simplification assumes hydrostatic vertical normal stresses (Greve & Blatter Reference Greve and Blatter2009), with the simplification that the pressure variable is eliminated from the Stokes equations. However, this hydrostatic approximation is almost never used in practice. Instead, the first-order approximation (FOA) (Blatter Reference Blatter1995; Pattyn Reference Pattyn2003) model, which further assumes negligible horizontal derivatives of the vertical velocity compared with vertical derivatives of the horizontal velocity, is usually preferred. Compared with the hydrostatic approximation, the third component of the velocity has disappeared from the FOA system. Despite the fact that the unknowns of the FOA are reduced to the horizontal components of the velocity, the FOA model is still mathematically three-dimensional (3D). As a consequence, solving the FOA still requires meshing complex and shallow geometries which change in time. To remove such complexity, the dimension of the mathematical model can be reduced by further mechanical simplifications. For instance, the shallow shelf approximation (SSA) (Morland Reference Morland, Veen and Oerlemans1987; MacAyeal Reference MacAyeal1989; Weis, Greve & Hutter Reference Weis, Greve and Hutter1999), which accounts only for longitudinal (or membrane) stresses, is depth-integrated and thus two-dimensional (2D). In contrast, the
${\it\epsilon}^{0}$
SIA, which accounts only for vertical shear stresses, reduces to a one-dimensional (1D; i.e. vertical) mathematical model, independently in each column of ice. In the literature (Schoof & Hindmarsh Reference Schoof and Hindmarsh2010), the SSA is called a ‘membrane’ model while the SIA is called a ‘lubrication’ model. Those two models (SIA and SSA) are popular for describing the dynamics of ice sheets and ice shelves since the size and the complexity of the system to solve are definitely reduced compared with any 3D models.
In practice, the vertical shear components of the stress tensor are significant where ice is frozen to the ground, while the longitudinal components are dominant in fast sliding and floating areas. As a consequence, both components must be combined if the entire domain is to be modelled. This has motivated the construction of ‘hybrid’ models, which account for both kind of stresses, while being mathematically 2D. The simplest hybrid model consists of the linear combination
$\text{SIA}+\text{SSA}$
, which is arrived at by adding together the velocities of each model (Bueler & Brown Reference Bueler and Brown2009). Unfortunately, this model does not include the simultaneous coupling between the vertical shear and the longitudinal stresses. As a result, this model cannot capture the 3D ice flows that occur in deep and narrow valleys or in the vicinity of grounding lines (Jouvet & Gräser Reference Jouvet and Gräser2013; Pattyn et al.
Reference Pattyn, Perichon, Durand, Favier, Gagliardini, Hindmarsh, Zwinger, Albrecht, Cornford, Docquier, Fürst, Goldberg, Gudmundsson, Humbert, Hütten, Huybrechts, Jouvet, Kleiner, Larour, Martin, Morlighem, Payne, Pollard, Rückamp, Rybak, Seroussi, Thoma and Wilkens2013). In contrast, the L1L2 (Hindmarsh Reference Hindmarsh2004) or some variants such as those proposed by Pollard & Deconto (Reference Pollard and Deconto2009), Schoof & Hindmarsh (Reference Schoof and Hindmarsh2010) or Goldberg (Reference Goldberg2011) include the vertical shear stress in the computation of the effective viscosity of the SSA. All these hybrid models have in common that they solve a single nonlinear elliptic 2D problem, and that the velocity profile is reconstructed a posteriori (Schoof & Hindmarsh Reference Schoof and Hindmarsh2010; Winkelmann et al.
Reference Winkelmann, Martin, Haseloff, Albrecht, Bueler, Khroulev and Levermann2011; Cornford et al.
Reference Cornford, Martin, Graves, Ranken, Brocq Le, Gladstone, Payne, Ng and Lipscomb2013). The hierarchy of the previously mentioned models is illustrated in figure 1.
Overview of the hierarchy of ice flow models: Stokes, FOA, SIA, multilayer model, ‘L1L2-like’ and SSA. The dimension of the mathematical model is indicated in the exponent.
In this paper, a new hybrid model generalising the SSA is introduced. The SSA model assumes a vertically constant velocity profile, such that it only accounts for the longitudinal components of the stress while neglecting the vertical components. To recover these components, the velocity profile of the new model is partitioned and assumed to be vertically piecewise constant. Inspired by an ocean model (Audusse et al. Reference Audusse, Bristeau, Perthame and Sainte-Marie2011), this approach consists of seeing the ice thickness as a pile of thin layers which can spread out, contract and slide over each other. Similar to the SSA, the model is obtained by integrating the FOA model vertically, and each of the layers locally. The boundary terms appear when integration gives rise to interlayer tractions. To account for sliding between layers, only the layer-orthogonal shear stress components are retained. These components redefine the interlayer tractions. The final model consists of a tridiagonal system of 2D nonlinear elliptic equations (defined later by (2.73)–(2.75)), whose number corresponds to the number of layers. By construction, this multilayer model naturally generalises the SSA, which corresponds to the one-layer case of the model. However, unlike the SSA, the multilayer is hybrid since it combines the longitudinal and the vertical shear stresses. (Note that in the literature (Hindmarsh Reference Hindmarsh2004; Egholm et al. Reference Egholm, Knudsen, Clark and Lesemann2011), the terminology ‘multilayer’ is sometimes used as a synonym of ‘hybrid’, as defined in this paper.) In addition, unlike ‘L1L2-like’ models (Hindmarsh Reference Hindmarsh2004), the two types of stresses are additively decoupled within each layer of the multilayer system.
The main difficulty in deriving the model consists of redefining the tractions describing the interlayer sliding. A zeroth-order traction (in the surface slope) can be redefined by simply retaining only the vertical shear stress components from the stress tensor like in the SIA approach. However, a more complex higher-order traction can be redefined. To do so, only the normal shear stress is kept in the local frame which is tangent to the interlayer interface. If the layers are chosen such that there are aligned with the streamlines, then, this redefinition of the traction is second-order in the interlayer surface slope. Interestingly, the second-order terms, which add no further computational complexity to the system, lead the multilayer model to surpass the FOA in the ‘infinite parallel-sided slab’ configuration. Indeed, in this simplified set-up the multilayer solution converges to the Stokes solution when increasing the number of layers. By contrast, increasing the vertical resolution in the FOA solution does not lead to convergence to Stokes solution, even in this configuration.
This paper is organised as follows: the model is derived in § 2 and then an analytical solution is devised for the ‘infinite parallel-sided slab’ setting in § 3. In § 4, a new insight into the multilayer approach is given by deriving the continuous equations from which the multilayer model with zeroth-order interlayer traction is a vertical semi-discretisation. Lastly, mechanical performance of the multilayer model for the flowline ISMIP-HOM experiments is tested against traditional higher-order models in § 5.
2. Model derivation
In this section, a generic 3D system of ice sheet and ice shelf is considered. The most complete viscous ice-flow model, nonlinear Stokes, and the FOA, which is a simplification, are described in §§ 2.1 and 2.2, respectively. Then an integration procedure derives the multilayer model in § 2.3. The redefinition of the interlayer tractions is reported in § 2.4 while the boundary conditions are rewritten in the multilayer setting in § 2.5. Lastly, the multilayer model is summarised in § 2.6.
Let
$V$
be a 3D domain of ice defined by
$$\begin{eqnarray}V=\{(x,y,z),\underline{s}(x,y)\leqslant z\leqslant \overline{s}(x,y)\},\end{eqnarray}$$
where
$(x,y)$
denote the horizontal coordinates,
$z$
denotes the vertical coordinate (positive upward), and
$\underline{s}(x,y),\overline{s}(x,y)$
are the elevations of the lower and upper ice surfaces, respectively. Call
$b(x,y)$
the elevation of the bedrock. Note that
$\underline{s}=b$
holds where ice is grounded and
$\underline{s}>b$
where ice is floating. The flotation of ice satisfies the Archimedes principle,
$$\begin{eqnarray}\underline{s}=\max \left\{b,-\frac{{\it\rho}}{{\it\rho}_{w}}h\right\},\end{eqnarray}$$
where
$h:=\overline{s}-\underline{s}$
is the ice thickness and
${\it\rho}$
and
${\it\rho}_{w}$
are the densities of ice and water, respectively; see figure 2. Relation (2.2) says that if the buoyancy
$-{\it\rho}_{w}gb$
is less than the force exerted by ice
${\it\rho}gh$
, then ice is grounded, otherwise ice is floating and
${\it\rho}/{\it\rho}_{w}$
of the ice thickness is below sea level.
Cross-section of an ice sheet and an ice shelf, with notation.
The boundary of
$V$
is divided into the upper interface
$$\begin{eqnarray}{\it\Gamma}_{s}=\{(x,y,z),z=\overline{s}(x,y),\underline{s}(x,y)<\overline{s}(x,y)\},\end{eqnarray}$$
the lower interface
$$\begin{eqnarray}{\it\Gamma}_{0}\cup {\it\Gamma}_{m}\cup {\it\Gamma}_{f}=\{(x,y,z),z=\underline{s}(x,y),\underline{s}(x,y)<\overline{s}(x,y)\}\end{eqnarray}$$
and the lateral boundary
${\it\Gamma}_{l}$
, which can include a possible vertical ice cliff at the calving front, see figure 2. At the lower interface, ice might be frozen to the ground, sliding on the ground or floating on the water:
${\it\Gamma}_{0}$
denotes the non-sliding part,
${\it\Gamma}_{m}$
the sliding grounded part and
${\it\Gamma}_{f}$
the floating part. In addition, for
$k\in \{s,0,m,f,l\}$
,
${\it\Omega}$
and
${\it\Omega}_{k}$
denote the projections of
$V$
and its boundaries
${\it\Gamma}_{k}$
on the horizontal plane
$(Oxy)$
, respectively.
In what follows, the velocity field of the ice fluid in
$V$
is denoted by
$\boldsymbol{u}=(u_{x},u_{y},u_{z})$
, the pressure field by
$p$
, the derivative for variable
$i$
by
$\partial _{i}~(i\in \{x,y,z\})$
, and Einstein’s summation convention is adopted.
2.1. Stokes approximation
The Stokes model comes from the momentum conservation equation when acceleration terms are ignored. If also the incompressibility condition is imposed, mass and momentum balance yield
$$\begin{eqnarray}\displaystyle & -\partial _{j}{\it\sigma}_{ij}={\it\rho}g_{i},\quad \text{in }V, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \partial _{i}u_{i}=0,\quad \text{in }V, & \displaystyle\end{eqnarray}$$
${\it\sigma}_{ij}$
are the components of the Cauchy stress tensor and
$(g_{x},g_{y},g_{z})=(0,0,-g)$
, where
$g$
is the gravitational constant. Let
${\it\tau}_{ij}$
be the components of the deviatoric stress tensor, and define
${\it\tau}_{ij}$
and the pressure
$p$
by
$-3p={\it\sigma}_{ii}$
and 
$$\begin{eqnarray}{\it\sigma}_{ij}={\it\tau}_{ij}-p{\it\delta}_{ij},\end{eqnarray}$$
where
${\it\delta}_{ij}$
is the Kronecker symbol. The mechanical behaviour of ice satisfies a viscosity relation:
$$\begin{eqnarray}{\it\tau}_{ij}=2{\it\mu}\unicode[STIX]{x1D60B}_{ij},\end{eqnarray}$$
where
$\unicode[STIX]{x1D60B}_{ij}$
denotes the components of the strain rate tensor defined by
$$\begin{eqnarray}\unicode[STIX]{x1D60B}_{ij}={\textstyle \frac{1}{2}}(\partial _{j}u_{i}+\partial _{i}u_{j}).\end{eqnarray}$$
The nonlinear viscosity
${\it\mu}$
is determined by the shear-thinning power law (Glen Reference Glen1953):
$$\begin{eqnarray}{\it\mu}={\textstyle \frac{1}{2}}A^{-1/n}\left|{\textstyle \frac{1}{2}}\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ji}\right|^{(1/n-1)/2},\end{eqnarray}$$
where
$A>0$
and
$n\geqslant 1$
are two parameters called the rate factor and Glen’s exponent, respectively. In reality,
$A$
is not constant since it depends on ice temperature (Glen Reference Glen1953; Fowler & Larson Reference Fowler and Larson1978; Paterson Reference Paterson1994). However, for the sake of simplicity, it is assumed in this paper that the ice is isothermal.
The boundary conditions that supplement (2.5)–(2.8) and (2.10) are the following. No force applies on the ice–air interface,
$$\begin{eqnarray}{\it\sigma}_{ij}n_{j}=0,\quad \text{on }{\it\Gamma}_{s},\end{eqnarray}$$
where
$$\begin{eqnarray}\boldsymbol{n}=(n_{x},n_{y},n_{z})^{\text{T}}=(-\partial _{x}\overline{s},-\partial _{y}\overline{s},1)^{\text{T}}\end{eqnarray}$$
is an outer normal vector along
${\it\Gamma}_{s}$
. Along the lower surface interface let
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{n}=(n_{x},n_{y},n_{z})^{\text{T}}=\frac{(\partial _{x}\underline{s},\partial _{y}\underline{s},-1)^{\text{T}}}{\sqrt{1+(\partial _{x}\underline{s})^{2}+(\partial _{y}\underline{s})^{2}}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{t}^{x}=(t_{x}^{x},t_{y}^{x},t_{z}^{x})^{\text{T}}=(1,0,\partial _{x}\underline{s})^{\text{T}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{t}^{y}=(t_{x}^{y},t_{y}^{y},t_{z}^{y})^{\text{T}}=(0,1,\partial _{y}\underline{s})^{\text{T}}, & \displaystyle\end{eqnarray}$$
$\boldsymbol{n}$
is normalised since this is needed later (in (2.17)) unlike
$\boldsymbol{t}^{x}$
and
$\boldsymbol{t}^{y}$
. The no-slip condition on the frozen base is 
$$\begin{eqnarray}u_{i}=0,\quad \text{on }{\it\Gamma}_{0}.\end{eqnarray}$$
A nonlinear friction condition applies where the ice is sliding (Hutter Reference Hutter1983):
$$\begin{eqnarray}\displaystyle & u_{i}n_{i}=0,\quad \text{on }{\it\Gamma}_{m}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\tau}_{ij}n_{j}t_{i}^{k}=-C|u_{j}u_{j}|^{(1/m-1)/2}u_{i}t_{i}^{k},\quad \text{on }{\it\Gamma}_{m}, & \displaystyle\end{eqnarray}$$
$k\in \{x,y\}$
, where
$m>0$
, and where
$C=C(x,y)>0$
is a predetermined distribution of basal strength in this paper. The condition on the floating interface reads (Weis et al.
Reference Weis, Greve and Hutter1999; Greve & Blatter Reference Greve and Blatter2009): 
$$\begin{eqnarray}{\it\sigma}_{ij}n_{i}={\it\rho}_{w}gzn_{j},\quad \text{on }{\it\Gamma}_{f}.\end{eqnarray}$$
Finally, essentially the same condition as (2.19) applies at the calving front below sea level, while the normal stress is zero at the grounded margins such that
$$\begin{eqnarray}{\it\sigma}_{ij}n_{i}={\it\rho}_{w}g\min (z,0)n_{j},\quad \text{on }{\it\Gamma}_{l}\end{eqnarray}$$
where
$(n_{x},n_{y},n_{z})=(n_{x},n_{y},0)$
is an outer normal vector to
${\it\Gamma}_{l}$
. Condition (2.20) also applies as stated to grounded margins where
$z>0$
so
${\it\sigma}_{ij}n_{i}=0$
.
2.2. FOA
Call
${\it\epsilon}=[h]/[x]$
the aspect ratio of
$V$
, where
$[h]$
and
$[x]$
denote its typical height and length of
$V$
. A dimensionless scaling (Blatter Reference Blatter1995; Schoof & Hindmarsh Reference Schoof and Hindmarsh2010) shows that
$$\begin{eqnarray}\displaystyle & \partial _{j}{\it\sigma}_{zj}=\partial _{z}{\it\sigma}_{zz}+O({\it\epsilon}^{2})\quad \text{in }V, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\sigma}_{zj}n_{j}={\it\sigma}_{zz}n_{z}+O({\it\epsilon}^{2}),\quad \text{on }{\it\Gamma}_{s}\cup {\it\Gamma}_{0}\cup {\it\Gamma}_{m}\cup {\it\Gamma}_{f}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\tau}_{zj}n_{j}={\it\tau}_{zz}n_{z}+O({\it\epsilon}^{2}),\quad \text{on }{\it\Gamma}_{s}\cup {\it\Gamma}_{0}\cup {\it\Gamma}_{m}\cup {\it\Gamma}_{f} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}D_{iz}={\textstyle \frac{1}{2}}\partial _{z}u_{i}+O({\it\epsilon}^{2}),\quad i\in \{x,y\}~\text{in }V.\end{eqnarray}$$
From now on (in §§ 2.2 and 2.3) the remainders
$O({\it\epsilon}^{2})$
in (2.21)–(2.24) are neglected such that the Stokes problem and its boundary conditions simplify. Indeed, using (2.21), the third equation of (2.5) simplifies into
$$\begin{eqnarray}\partial _{z}{\it\sigma}_{zz}={\it\rho}g,\quad \text{in }V,\end{eqnarray}$$
which is a hydrostatic approximation, while, using (2.22), the third equation of (2.11) becomes
$$\begin{eqnarray}{\it\sigma}_{zz}=0,\quad \text{on }{\it\Gamma}_{s}.\end{eqnarray}$$
Integrating vertically (2.25) with (2.26) yields
$$\begin{eqnarray}{\it\sigma}_{zz}={\it\tau}_{zz}-p=-{\it\rho}g(\overline{s}-z),\quad \text{in }V.\end{eqnarray}$$
By (2.6) and (2.8), it follows that
$$\begin{eqnarray}p={\it\rho}g(\overline{s}-z)-{\it\tau}_{xx}-{\it\tau}_{yy},\quad \text{in }V.\end{eqnarray}$$
Thus,
$p$
can be eliminated from the two first equations of (2.5):
$$\begin{eqnarray}\displaystyle & \partial _{x}(2{\it\tau}_{xx}+{\it\tau}_{yy})+\partial _{y}{\it\sigma}_{xy}+\partial _{z}{\it\sigma}_{xz}={\it\rho}g\partial _{x}\overline{s}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \partial _{x}{\it\sigma}_{xy}+\partial _{y}(2{\it\tau}_{yy}+{\it\tau}_{xx})+\partial _{z}{\it\sigma}_{yz}={\it\rho}g\partial _{y}\overline{s}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & (2{\it\tau}_{xx}+{\it\tau}_{yy})n_{x}+({\it\sigma}_{xy})n_{y}+{\it\sigma}_{xz}n_{z}=0,\quad \text{on }{\it\Gamma}_{s}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & ({\it\sigma}_{xy})n_{x}+(2{\it\tau}_{yy}+{\it\tau}_{xx})n_{y}+{\it\sigma}_{yz}n_{z}=0,\quad \text{on }{\it\Gamma}_{s}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & (2{\it\tau}_{xx}+{\it\tau}_{yy})n_{x}+({\it\sigma}_{xy})n_{y}+{\it\sigma}_{xz}n_{z}=-C|u_{j}u_{j}|^{(1/m-1)/2}u_{i}t_{i}^{x},\quad \text{on }{\it\Gamma}_{m}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & ({\it\sigma}_{yx})n_{x}+(2{\it\tau}_{yy}+{\it\tau}_{xx})n_{y}+{\it\sigma}_{yz}n_{z}=-C|u_{j}u_{j}|^{(1/m-1)/2}u_{i}t_{i}^{y},\quad \text{on }{\it\Gamma}_{m}. & \displaystyle\end{eqnarray}$$
$\underline{s}=-({\it\rho}/{\it\rho}_{w})h$
, which derives from the floating condition (2.2), the third equation of (2.19) becomes 
$$\begin{eqnarray}{\it\sigma}_{zz}={\it\tau}_{zz}-p={\it\rho}_{w}g\underline{s}=-{\it\rho}gh,\quad \text{on }{\it\Gamma}_{f}.\end{eqnarray}$$
Again using (2.28) to eliminate the pressure from the two first equations of (2.19), using (2.2) and (2.35) imply
$$\begin{eqnarray}\displaystyle & (2{\it\tau}_{xx}+{\it\tau}_{yy})n_{x}+({\it\sigma}_{xy})n_{y}+{\it\tau}_{xz}n_{z}=0,\quad \text{on }{\it\Gamma}_{f}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & ({\it\sigma}_{yx})n_{x}+(2{\it\tau}_{yy}+{\it\tau}_{xx})n_{y}+{\it\tau}_{yz}n_{z}=0,\quad \text{on }{\it\Gamma}_{f}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & (2{\it\tau}_{xx}+{\it\tau}_{yy})n_{x}+({\it\sigma}_{xy})n_{y}=({\it\rho}_{w}g\,\min (z,0)+{\it\rho}g\,(\overline{s}-z))n_{x},\quad \text{on }{\it\Gamma}_{l}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & ({\it\sigma}_{yx})n_{x}+(2{\it\tau}_{yy}+{\it\tau}_{xx})n_{y}=({\it\rho}_{w}g\min (z,0)+{\it\rho}g(\overline{s}-z))n_{y},\quad \text{on }{\it\Gamma}_{l}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\it\mu} & = & \displaystyle {\textstyle \frac{1}{2}}A^{-1/n}\Big[{\textstyle \frac{1}{2}}(\partial _{x}u_{x})^{2}+{\textstyle \frac{1}{2}}(\partial _{y}u_{y})^{2}+{\textstyle \frac{1}{2}}(\partial _{x}u_{x}+\partial _{y}u_{y})^{2}\nonumber\\ \displaystyle & & \displaystyle +\,{\textstyle \frac{1}{4}}(\partial _{y}u_{x}+\partial _{x}u_{y})^{2}+{\textstyle \frac{1}{4}}(\partial _{z}u_{x})^{2}+{\textstyle \frac{1}{4}}(\partial _{z}u_{y})^{2}\Big]^{(1/n-1)/2}.\end{eqnarray}$$
In summary, in the FOA model a small-aspect-ratio approximation yields a hydrostatic approximation for the normal stress (2.27), equivalently relation (2.28) for the pressure, and this in turn allows the pressure to be eliminated from all equations. Furthermore the model reduces to two scalar equations for the horizontal velocity, but in three variables
$(x,y,z)$
, so that numerical schemes for the FOA must discretise a 3D domain. Instead of doing so, in the next section, the FOA is first vertically integrated over a pile of layers which cover the ice thickness. In contrast with the FOA, the resulting model, called multilayer, only horizontal discretisation is needed.
2.3. Depth integration over layers
The domain of ice is now divided in the vertical direction into
$L$
layers of thickness
$h^{1},\dots ,h^{L}$
such that
$$\begin{eqnarray}\mathop{\sum }_{l=1,\dots ,L}h^{l}=h,\end{eqnarray}$$
see figure 3. Call
$s^{l}=\underline{s}+\sum _{j=0}^{l}h^{\,j}$
the elevation of the upper surface of layer
$l$
for
$l=0,\dots ,L$
, with the convention
$h^{0}=0$
. Several modes of division of the ice thickness are possible, including uniform depth defined by
$h^{l}=h/L$
. However, a vertical division (2.41) chosen such that layers are aligned with the streamlines will yield better model accuracy, see § 2.4.
Multilayer splitting of the ice thickness.
The derivation of the SSA model is based on the assumption of a vertically constant velocity profile (Morland Reference Morland, Veen and Oerlemans1987; MacAyeal Reference MacAyeal1989; Weis et al.
Reference Weis, Greve and Hutter1999). Instead, here
$(u_{x},u_{y})$
is assumed to be vertically layer-wise constant, equal to
$(u_{x}^{l},u_{y}^{l})$
on layer
$l$
:
$$\begin{eqnarray}u_{k}(x,y,z)=\mathop{\sum }_{l=1,\dots ,L}u_{k}^{l}(x,y)\mathbf{1}_{(s^{l-1},s^{l}]}(z),\end{eqnarray}$$
for
$k\in \{x,y\}$
, where
$\mathbf{1}_{I}(z)$
equals 1 if
$z\in I$
and 0 otherwise, see figure 3. The discontinuities of the velocity lead to undefined tractions between the layers. The redefinition of such tractions is addressed in § 2.4.
Consider an arbitrary layer indexed by
$l\in \{1,\dots ,L\}$
. Using Leibniz’s rule, the integration of (2.29) vertically over the layer
$l$
yields
$$\begin{eqnarray}\displaystyle & & \displaystyle 2\partial _{x}\left(\int _{s^{l-1}}^{s^{l}}{\it\tau}_{xx}\,\text{d}z\right)-2[{\it\tau}_{xx}]_{z=s^{l}}\partial _{x}s^{l}+2[{\it\tau}_{xx}]_{z=s^{l-1}}\partial _{x}s^{l-1}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\partial _{x}\left(\int _{s^{l-1}}^{s^{l}}{\it\tau}_{yy}\,\text{d}z\right)-[{\it\tau}_{yy}]_{z=s^{l}}\partial _{x}s^{l}+[{\it\tau}_{yy}]_{z=s^{l-1}}\partial _{x}s^{l-1}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\partial _{y}\left(\int _{s^{l-1}}^{s^{l}}{\it\sigma}_{xy}\,\text{d}z\right)-[{\it\sigma}_{xy}]_{z=s^{l}}\partial _{y}s^{l}+[{\it\sigma}_{xy}]_{z=s^{l-1}}\partial _{y}s^{l-1}\nonumber\\ \displaystyle & & \displaystyle \quad +\,[{\it\sigma}_{xz}]_{z=s^{l}}-[{\it\sigma}_{xz}]_{z=s^{l-1}}={\it\rho}gh^{l}\partial _{x}\overline{s},\end{eqnarray}$$
where
$[\cdot ]_{z=s^{l}}$
(respectively
$[\cdot ]_{z=s^{l-1}}$
) stands for the limit
$z\mapsto s^{l}$
(respectively
$z\mapsto s^{l-1}$
) with
$z<s^{l}$
(respectively
$z>s^{l-1}$
). Using the fact that
$(u_{x}^{l},u_{y}^{l})$
is
$z$
-independent in the layer
$l$
, equation (2.43) takes the form:
$$\begin{eqnarray}\partial _{x}\left(h^{l}\left(2{\it\tau}_{xx}+{\it\tau}_{yy}\right)\right)+\partial _{y}\left(h^{l}\left({\it\sigma}_{xy}\right)\right)+{\it\Sigma}_{x}^{l,0}+{\it\Sigma}_{x}^{l,-1}={\it\rho}gh^{l}\partial _{x}\overline{s},\end{eqnarray}$$
and, for
$q=0,-1$
,
$$\begin{eqnarray}{\it\Sigma}_{x}^{l,q}=(-1)^{-q}{\it\alpha}^{l+q}\left[(2{\it\tau}_{xx}+{\it\tau}_{yy})n_{x}^{l+q}+{\it\sigma}_{xy}n_{y}^{l+q}+{\it\sigma}_{xz}n_{z}^{l+q}\right]_{z=s^{l+q}},\end{eqnarray}$$
where
$$\begin{eqnarray}\boldsymbol{n}^{l}=(n_{x}^{l},n_{y}^{l},n_{z}^{l})^{\text{T}}=\frac{(\partial _{x}s^{l},\partial _{y}s^{l},-1)^{\text{T}}}{{\it\alpha}^{l}},\end{eqnarray}$$
is the outer normal unit vector to the upper boundary of layer
$l$
, and
$$\begin{eqnarray}{\it\alpha}^{l}=\sqrt{1+(\partial _{x}s^{l})^{2}+(\partial _{y}s^{l})^{2}}.\end{eqnarray}$$
Similarly, integrating (2.30) vertically over the layer
$l\in \{1,\dots ,L\}$
leads to
$$\begin{eqnarray}\partial _{x}\left(h^{l}\left({\it\sigma}_{xy}\right)\right)+\partial _{y}\left(h^{l}\left({\it\tau}_{xx}+2{\it\tau}_{yy}\right)\right)+{\it\Sigma}_{y}^{l,0}+{\it\Sigma}_{y}^{l,-1}={\it\rho}gh^{l}\partial _{y}\overline{s},\end{eqnarray}$$
where, for
$q=0,-1$
,
$$\begin{eqnarray}{\it\Sigma}_{y}^{l,q}=(-1)^{-q}{\it\alpha}^{l+q}\left[{\it\sigma}_{yx}n_{x}^{l+q}+(2{\it\tau}_{yy}+{\it\tau}_{xx})n_{y}^{l+q}+{\it\sigma}_{yz}n_{z}^{l+q}\right]_{z=s^{l+q}}.\end{eqnarray}$$
Using the simplification due to (2.23),
${\it\Sigma}_{k}^{l,0}$
and
${\it\Sigma}_{k}^{l,-1}$
correspond to the tractions at the top and bottom of layer
$l$
:
$$\begin{eqnarray}{\it\Sigma}_{k}^{l,q}=(-1)^{-q}{\it\alpha}^{l+q}\left[{\it\tau}_{ij}n_{i}^{l+q}t_{j}^{k,l+q}\right]_{z=s^{l+q}},\end{eqnarray}$$
where
$\boldsymbol{n}^{l}$
is defined by (2.46) and
$$\begin{eqnarray}\displaystyle & \boldsymbol{t}^{x,l}=(t_{x}^{x,l},t_{y}^{x,l},t_{z}^{x,l})^{\text{T}}=(1,0,\partial _{x}s^{l})^{\text{T}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{t}^{y,l}=(t_{x}^{y,l},t_{y}^{y,l},t_{z}^{y,l})^{\text{T}}=(0,1,\partial _{y}s^{l})^{\text{T}} & \displaystyle\end{eqnarray}$$
$l$
.
2.4. Interlayer tractions
The continuity of the stress across the layers implies
$$\begin{eqnarray}{\it\Sigma}_{k}^{l,0}=-{\it\Sigma}_{k}^{l+1,-1},\quad \text{for all }l=1,\dots ,L-1,~\text{for all }k=x,y,\end{eqnarray}$$
such that the layers are coupled. From (2.50), (2.8) and (2.10),
${\it\Sigma}_{k}^{l,0}$
and
$-{\it\Sigma}_{k}^{l+1,-1}$
should be equal to
$$\begin{eqnarray}A^{-1/n}{\it\alpha}^{l}\left|{\textstyle \frac{1}{2}}\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ji}\right|^{(1/n-1)/2}\unicode[STIX]{x1D60B}_{ij}n_{j}^{l}t_{i}^{k,l}.\end{eqnarray}$$
However,
$\unicode[STIX]{x1D60B}_{ij}$
is not defined between the layers because of the discontinuity of the velocity field (2.42). The goal of this section is to redefine
${\it\Sigma}_{k}^{l,0}$
and
$-{\it\Sigma}_{k}^{l+1,-1}$
into a meaningful quantity called
$S_{k}^{l}$
. In what follows, two possible redefinitions are described: a simple one which is zeroth order in the surface slope defined by (2.57) and a more complex second-order one defined by (2.59).
Both redefinitions are based on the following key hypothesis: the layer-parallel components of the stress are set at zero:
$$\begin{eqnarray}\unicode[STIX]{x1D60B}_{ij}T_{ix}T_{jx}=0,\quad \unicode[STIX]{x1D60B}_{ij}T_{ix}T_{jy}=0,\quad \unicode[STIX]{x1D60B}_{ij}T_{iy}T_{jy}=0,\end{eqnarray}$$
where
$\{T_{ix}\}_{i}$
and
$\{T_{jy}\}_{j}$
are two independent vectors orthogonal to
$\boldsymbol{n}^{l}$
. We remark that, in (2.55),
$\unicode[STIX]{x1D60B}_{ij}T_{ix}T_{jx}$
,
$\unicode[STIX]{x1D60B}_{ij}T_{ix}T_{jy}$
and
$\unicode[STIX]{x1D60B}_{ij}T_{iy}T_{jy}$
equal the
$(x,x)$
,
$(x,y)$
and
$(y,y)$
components of the matrix
$\unicode[STIX]{x1D60B}_{ij}$
, however, expressed in the local frame generated by
$\{T_{ix}\}_{i}^{l}$
,
$\{T_{jy}\}_{j}^{l}$
, and
$\boldsymbol{n}^{l}$
, and then aligned with the upper surface of layer
$l$
. As a consequence, (2.55) says that the matrix
$\unicode[STIX]{x1D60B}_{ij}$
(or equivalently the stress tensor) expressed in the local frame contains only four non-zero entries that are
$(x,z)$
,
$(y,z)$
,
$(z,x)$
and
$(z,y)$
, see appendix A. This is justified by the fact that the layer-orthogonal shear components dominate in the stress expressed in the local frame tangential to the layer boundary since the layers can slide on each other.
2.4.1. Zeroth-order interlayer tractions to define the multilayer
$^{\ast }$
model
To define
$S_{k}^{l}$
at zeroth order, one neglects
$O({\it\delta})$
terms where
${\it\delta}$
denotes the scale for the slopes in the surface elevation of layers, such that
$\boldsymbol{n}^{l}=(0,0,1)$
,
$\boldsymbol{t}^{x,l}=(1,0,0)$
and
$\boldsymbol{t}^{y,l}=(0,1,0)$
. By using (2.55) and approximating the derivative with respect to
$z$
by the finite difference:
$$\begin{eqnarray}\frac{\partial _{z}u_{i}}{2}=\frac{u_{i}^{l+1}-u_{i}^{l}}{h^{l+1}+h^{l}},\quad i\in \{x,y\},\end{eqnarray}$$
${\it\Sigma}_{k}^{l,0}$
and
$-{\it\Sigma}_{k}^{l+1,-1}$
can be redefined by
$$\begin{eqnarray}S_{k}^{l}=A^{-1/n}\left|\left(\frac{u_{x}^{l+1}-u_{x}^{l}}{h^{l+1}+h^{l}}\right)^{2}+\left(\frac{u_{y}^{l+1}-u_{y}^{l}}{h^{l+1}+h^{l}}\right)^{2}\right|^{(1/n-1)/2}\left(\frac{u_{k}^{l+1}-u_{k}^{l}}{h^{l+1}+h^{l}}\right).\end{eqnarray}$$
Later, the multilayer model based on the simplified zeroth-order interlayer traction (2.57) is labelled ‘multilayer
$^{\ast }$
’. This approximation is similar to that involved when deriving the SIA (Hutter Reference Hutter1983; Greve & Blatter Reference Greve and Blatter2009); see also § 3.
2.4.2. Second-order interlayer tractions to define the general multilayer model
In order to redefine the interlayer traction
$S_{k}^{l}$
at a higher order in
${\it\delta}$
, two further hypothesis are formulated. First, the multilayer vertical splitting (2.41) is chosen such that ‘the layers are aligned with the three-dimensional direction of the flow’:
$$\begin{eqnarray}u_{j}n_{j}^{l}=0.\end{eqnarray}$$
Second, ‘the layers are sufficiently thin’ such that the typical aspect ratio of layers
${\it\epsilon}/L$
is small compared with
${\it\delta}$
. As a consequence, terms
$O({\it\delta}{\it\varepsilon}/L)$
are neglected, while terms
$O({\it\delta}^{2})$
are retained. Using those two further assumptions, it is shown in appendix A that the second-order interlayer tractions
${\it\Sigma}_{k}^{l,0}$
and
$-{\it\Sigma}_{k}^{l+1,-1}$
can be redefined by
$$\begin{eqnarray}S_{k}^{l}=2({\it\alpha}^{l})^{2}{\it\nu}^{l}\left[M_{ik}^{l}\left(\frac{u_{i}^{l+1}-u_{i}^{l}}{h^{l}+h^{l+1}}\right)\right],\quad i\in \{x,y\}\end{eqnarray}$$
where
${\it\alpha}^{l}$
is defined by (2.47),
$$\begin{eqnarray}{\it\nu}^{l}=\frac{1}{2}A^{-1/n}\left({\it\alpha}^{l}\right)^{1/n-1}\left|M_{ij}^{l}\left(\frac{u_{i}^{l+1}-u_{i}^{l}}{h^{l}+h^{l+1}}\right)\left(\frac{u_{j}^{l+1}-u_{j}^{l}}{h^{l}+h^{l+1}}\right)\right|^{(1/n-1)/2},\quad i,j\in \{x,y\},\end{eqnarray}$$
and
$$\begin{eqnarray}M_{ij}^{l}={\it\delta}_{ij}+(\partial _{i}s^{l})(\partial _{j}s^{l}),\quad i,j\in \{x,y\},\end{eqnarray}$$
where
${\it\delta}_{ij}$
is Kronecker’s delta. Obviously, if one neglects
$O({\it\delta}^{2})$
terms in (2.59) and (2.60), then
$$\begin{eqnarray}{\it\alpha}^{l}=1,\quad M_{ij}^{l}={\it\delta}_{ij},\end{eqnarray}$$
such that (2.59) reduces to (2.57), which we have labelled as ‘multilayer
$^{\ast }$
’. The second-order interlayer traction (2.59) defines the most general model, and is labelled ‘multilayer’ (without star) in what follows. In fact,
$S_{k}^{l}$
defined in (2.57) corresponds to the vertical shear stresses in the global frame, assuming the other components to be zero. If the surface gradients are small (i.e. if
${\it\delta}^{2}$
is small), then the global frame is close enough to the local frame tangential to the layer boundary so that (2.57) can be used to redefine
$S_{k}^{l}$
. In contrast, if these gradients are not negligible, the slope at the layer interface must be accounted for, and so
$S_{k}^{l}$
is better used in the form (2.59); see § 5.
2.5. Boundary conditions
The boundary condition at the top of the highest layer, expressible in terms of
${\it\Sigma}_{k}^{L,0}$
, is considered first. Returning to (2.45) and the free-stress condition (2.31), (2.32), it follows that
$$\begin{eqnarray}{\it\Sigma}_{k}^{L,0}=0,\quad \text{for }k=x,y.\end{eqnarray}$$
The boundary conditions expressing no-slip or sliding at the bottom of the lowest layer via
${\it\Sigma}_{k}^{1,-1}$
are now considered. On
${\it\Omega}_{f}$
, no sliding between the lowest layer and the base occurs, and the approach in § 2.4 can be applied considering the bedrock as a fixed layer 0, i.e. with
$u_{x}^{0}=u_{y}^{0}=0$
. Thus,
${\it\Sigma}_{k}^{1,-1}$
takes the form
$$\begin{eqnarray}S_{k}^{0}=2({\it\alpha}^{0})^{2}{\it\nu}^{0}\left[M_{ik}^{0}\left(\frac{u_{i}^{1}-u_{i}^{0}}{h^{1}+h^{0}}\right)\right]=2({\it\alpha}^{0})^{2}{\it\nu}^{0}\left[M_{ik}^{0}\left(\frac{u_{i}^{1}}{h^{1}}\right)\right],\quad i\in \{x,y\}\end{eqnarray}$$
where
$$\begin{eqnarray}{\it\nu}^{0}=\frac{1}{2}A^{-1/n}\left({\it\alpha}^{0}\right)^{1/n-1}\left|M_{ij}^{l}\left(\frac{u_{i}^{1}}{h^{1}}\right)\left(\frac{u_{j}^{1}}{h^{1}}\right)\right|^{(1/n-1)/2},\quad i,j\in \{x,y\}\end{eqnarray}$$
since
$h^{0}=0$
and
$u_{i}^{0}=0$
by convention. On the sliding part
${\it\Omega}_{m}$
, conditions (2.33) and (2.34) with (2.45) lead to
$$\begin{eqnarray}{\it\Sigma}_{k}^{1,-1}=C\;{\it\alpha}^{0}\;|u_{j}u_{j}|^{(1/m-1)/2}\left(u_{i}^{1}t_{i}^{0,k}\right),\quad i,j\in \{x,y\}.\end{eqnarray}$$
Now (2.66) can be rewritten into the horizontal velocity components. Since the velocity field is tangential to
$(\boldsymbol{t}^{0,x},\boldsymbol{t}^{0,y})$
(condition (2.17)), it follows,
$$\begin{eqnarray}{\it\Sigma}_{k}^{1,-1}=\left({\it\alpha}^{0}\right)\bar{{\it\nu}}^{0}\left[M_{ik}^{0}\left(u_{i}^{1}\right)\right],\quad i\in \{x,y\},\end{eqnarray}$$
where
$$\begin{eqnarray}\bar{{\it\nu}}^{0}=C\left|M_{ij}^{0}\left(u_{i}^{1}\right)\left(u_{j}^{1}\right)\right|^{(1/n-1)/2},\quad i,j\in \{x,y\}.\end{eqnarray}$$
Finally, the conditions (2.36) and (2.37) on the floating part
${\it\Omega}_{f}$
become
$$\begin{eqnarray}{\it\Sigma}_{k}^{1,-1}=0,\quad \text{for }k=x,y,~\text{on }{\it\Omega}_{l}.\end{eqnarray}$$
At the calving front or at the glacier margins
${\it\Omega}_{l}$
, equations (2.38) and (2.39) are integrated from
$s^{l-1}$
to
$s^{l}$
, and the following boundary condition is obtained:
$$\begin{eqnarray}\displaystyle & h^{l}\left(2{\it\tau}_{xx}+{\it\tau}_{yy}\right)n_{x}+h^{l}\left({\it\sigma}_{xy}\right)n_{y}=F^{l}n_{x}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & h^{l}\left({\it\sigma}_{xy}\right)n_{x}+h^{l}\left({\it\tau}_{xx}+2{\it\tau}_{yy}\right)n_{y}=F^{l}n_{y}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}F^{l}={\textstyle \frac{1}{2}}{\it\rho}g\left[(\overline{s}-s^{l})^{2}-(\overline{s}-s^{l+1})^{2}\right]+{\textstyle \frac{1}{2}}{\it\rho}_{w}g\left[(\min (s^{l+1},0))^{2}-(\min (s^{l},0))^{2}\right].\end{eqnarray}$$
2.6. Summary
From now on,
${\it\Sigma}_{k}^{l,0}$
and
$-{\it\Sigma}_{k}^{l+1,-1}$
defined by (2.50) are replaced by
$S_{k}^{l}$
defined by (2.59). Finally, using (2.7), (2.8), (2.10), the multilayer solution
$u_{k}^{l}$
solves the following
$2\times 2$
-block tridiagonal system of equations:
$$\begin{eqnarray}\partial _{j}\left(2{\it\mu}^{L}h^{L}\left(\frac{\partial _{j}u_{k}^{L}+\partial _{k}u_{j}^{L}}{2}+(\partial _{i}u_{i}^{l}){\it\delta}_{jk}\right)\right)-2({\it\alpha}^{L-1})^{2}{\it\nu}^{L-1}\left[M_{ik}^{L-1}\left(\frac{u_{i}^{L}-u_{i}^{L-1}}{h^{L}+h^{L-1}}\right)\right]={\it\rho}gh^{L}\partial _{k}\overline{s},\end{eqnarray}$$
for all
$l\in \{2,\dots ,L-1\}$
:
$$\begin{eqnarray}\displaystyle & & \displaystyle \partial _{j}\left(2{\it\mu}^{l}h^{l}\left(\frac{\partial _{j}u_{k}^{l}+\partial _{k}u_{j}^{l}}{2}+(\partial _{i}u_{i}^{l}){\it\delta}_{jk}\right)\right)-2({\it\alpha}^{l})^{2}{\it\nu}^{l}\left[M_{ik}^{l}\left(\frac{u_{i}^{l}-u_{i}^{l+1}}{h^{l}+h^{l+1}}\right)\right]\nonumber\\ \displaystyle & & \displaystyle \quad -\,2({\it\alpha}^{l-1})^{2}{\it\nu}^{l-1}\left[M_{ik}^{l-1}\left(\frac{u_{i}^{l}-u_{i}^{l-1}}{h^{l}+h^{l-1}}\right)\right]={\it\rho}gh^{l}\partial _{k}\overline{s},\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle & & \displaystyle \partial _{j}\left(2{\it\mu}^{1}h^{1}\left(\frac{\partial _{j}u_{k}^{1}+\partial _{k}u_{j}^{1}}{2}+(\partial _{i}u_{i}^{1}){\it\delta}_{jk}\right)\right)-2({\it\alpha}^{1})^{2}{\it\nu}^{l}\left[M_{ik}^{l}\left(\frac{u_{i}^{1}-u_{i}^{2}}{h^{1}+h^{2}}\right)\right]\nonumber\\ \displaystyle & & \displaystyle \quad -\,2({\it\alpha}^{0})^{2}{\it\nu}^{0}\left[M_{ik}^{0}\left(\frac{u_{i}^{1}}{h^{1}}\right)\right]\times \mathbf{1}_{{\it\Omega}_{f}}-({\it\alpha}^{0})\bar{{\it\nu}}^{0}\left[M_{ik}^{0}\left(u_{i}^{1}\right)\right]\times \mathbf{1}_{{\it\Omega}_{m}}={\it\rho}gh^{l}\partial _{k}\overline{s},\end{eqnarray}$$
where
$i,j,k\in \{x,y\}$
, with summation convention over
$i,j$
,
$$\begin{eqnarray}{\it\mu}^{l}={\textstyle \frac{1}{2}}A^{-1/n}\left[{\textstyle \frac{1}{2}}\left(\partial _{x}u_{x}^{l}\right)^{2}+{\textstyle \frac{1}{2}}\left(\partial _{y}u_{y}^{l}\right)^{2}+{\textstyle \frac{1}{2}}\left(\partial _{x}u_{x}^{l}+\partial _{y}u_{y}^{l}\right)^{2}+{\textstyle \frac{1}{4}}\left(\partial _{y}u_{x}^{l}+\partial _{x}u_{y}^{l}\right)^{2}\right]^{(1/n-1)/2},\end{eqnarray}$$
and
$\mathbf{1}_{R}(x,y)$
equals 1 if
$(x,y)\in R$
and 0 otherwise. At the calving front or at the glacier margins, equations (2.70) and (2.71) can be reformulated by
$$\begin{eqnarray}\left(2{\it\mu}^{l}h^{l}\left(\frac{\partial _{j}u_{k}^{l}+\partial _{k}u_{j}^{l}}{2}+(\partial _{i}u_{i}^{l}){\it\delta}_{jk}\right)\right)n_{j}=F^{l}n_{k},\quad \text{on }{\it\Gamma}_{l},\end{eqnarray}$$
where
$i,j,k\in \{x,y\}$
, with summation convention over
$i,j$
, and
$F^{l}$
is defined by (2.72).
The system (2.73)–(2.75) is similar to the equation of the SSA (MacAyeal Reference MacAyeal1989; Schoof Reference Schoof2006), which corresponds to the one-layer model (i.e. when
$L=1$
). Indeed, the SSA model consists of a single elliptic nonlinear equation while the multilayer model consists of a system of elliptic nonlinear equations. Unlike the SSA, the system (2.73)–(2.75) has additional terms, which couple the layers, and which represent the normal shear stresses. In contrast with other hybrid models like the L1L2 (Hindmarsh Reference Hindmarsh2004) or the variants proposed in Pollard & Deconto (Reference Pollard and Deconto2009) or Schoof & Hindmarsh (Reference Schoof and Hindmarsh2010), the terms for the longitudinal and the vertical shear stresses are decoupled in an additional way. Note that (2.75), which applies on the lowest layer, includes no-slip and sliding conditions. When neglecting the
$O({\it\delta}^{2})$
components, the friction term of (2.75) reduces to the common expression
$C|(u_{x}^{1})^{2}+(u_{y}^{1})^{2}|^{(1/m-1)/2}u_{k}^{1}\times \mathbf{1}_{{\it\Omega}_{m}}$
, see Cornford et al. (Reference Cornford, Martin, Graves, Ranken, Brocq Le, Gladstone, Payne, Ng and Lipscomb2013), for example.
Finally, it must be stressed that the vertical discretisation of the multilayer system is momentum-conservative by construction (because of (2.53)).
3. Infinite parallel-sided slab
In this section, an analytic solution of the multilayer system is found for the ‘infinite parallel-sided slab’ and compared with the solution of the Stokes equations. The simplified setting (Hutter Reference Hutter1983; Greve & Blatter Reference Greve and Blatter2009) relies on the following assumptions.
-
(i) The flow is 2D in the vertical
$x{-}z$
plane (no
$y$
dependency) and the horizontal domain is infinite
${\it\Omega}=\, ]-\infty ,+\infty [$
. -
(ii) The bedrock slope
$\partial _{x}b$
and the ice thickness
$h$
are constants, such that
$\partial _{x}b=\partial _{x}s$
. -
(iii) The no-slip condition (2.16) applies everywhere on the bedrock.
-
(iv) The bedrock is everywhere above sea level such that there is no floating ice shelf.
Here, the vertical splitting is chosen uniform, defined by
$h^{l}=h/L$
. Since the geometry shows no
$x$
variation, the solution along the layers
$(u_{x}^{1},\dots ,u_{x}^{L})$
is independent of
$x$
. All derivatives can be removed in the multilayer system (2.73)–(2.75), which reads: find
$(u^{1},\dots ,u^{L})$
so that
$$\begin{eqnarray}\left.\begin{array}{@{}rcl@{}}\displaystyle -{\it\alpha}^{2}\left({\it\alpha}^{2}\frac{u_{x}^{L}-u_{x}^{L-1}}{2A(h/L)}\right)^{1/n} & = & {\it\rho}g(h/L)(\partial _{x}b)\\ \cdots \, & & \cdots \\ \displaystyle -{\it\alpha}^{2}\left({\it\alpha}^{2}\frac{u_{x}^{l}-u_{x}^{l-1}}{2A(h/L)}\right)^{1/n}+{\it\alpha}^{2}\left({\it\alpha}^{2}\frac{u_{x}^{l+1}-u_{x}^{l}}{2A(h/L)}\right)^{1/n} & = & {\it\rho}g(h/L)(\partial _{x}b)\\ \cdots \, & & \cdots \\ \displaystyle -{\it\alpha}^{2}\left({\it\alpha}^{2}\frac{u_{x}^{1}}{A(h/L)}\right)^{1/n}+{\it\alpha}^{2}\left({\it\alpha}^{2}\frac{u_{x}^{2}-u_{x}^{1}}{2A(h/L)}\right)^{1/n} & = & {\it\rho}g(h/L)(\partial _{x}b),\end{array}\right\}\end{eqnarray}$$
where
${\it\alpha}=\sqrt{1+(\partial _{x}b)^{2}}$
and
${\it\alpha}=1$
for the multilayer and the simplified multilayer
$^{\ast }$
models, respectively.
It is easy to verify that the solution of the system above is
$$\begin{eqnarray}u_{x}^{l}=-\frac{2A({\it\rho}g(\partial _{x}b))^{n}}{{\it\alpha}^{2(n+1)}}\left(\frac{h}{L}\right)\left(\frac{1}{2}(h)^{n}+\left(\frac{h(L-1)}{L}\right)^{n}+\cdots +\left(\frac{h(L-l+1)}{L}\right)^{n}\right).\end{eqnarray}$$
Formula (3.2) also arises from the quadrature of
$$\begin{eqnarray}\partial _{z}u_{x}=-\frac{2A[{\it\rho}g(\partial _{x}b)(s-z)]^{n}}{{\it\alpha}^{2(n+1)}},\end{eqnarray}$$
with a rectangle rule on each layer. As a matter of fact, the Stokes solution also satisfies (3.3), see Greve & Blatter (Reference Greve and Blatter2009, p. 146). As a consequence, the
$L\rightarrow \infty$
limit of the multilayer solution (3.2) equals the exact solution of the Stokes system. In contrast the limit of
$\text{multilayer}^{\ast }~({\it\alpha}=1)$
solution equals the SIA solution (Hutter Reference Hutter1983; Greve & Blatter Reference Greve and Blatter2009)
$$\begin{eqnarray}\partial _{z}u_{x}=-2A[{\it\rho}g(\partial _{x}b)(s-z)]^{n}.\end{eqnarray}$$
Thus, although the multilayer solution (3.3) converges to the exact solution of the Stokes equations when refining the vertical multilayer splitting (2.41), this convergence does not hold with the FOA since the remainders (2.21)–(2.24) are non-zero. (In fact, one can show that the FOA solution equals the Stokes solution divided by a factor
$(1+4(\partial _{x}s^{l})^{2}+4(\partial _{y}s^{l})^{2})^{2}$
.) This shows that in this special case the reconstruction of the interlayer traction (appendix A) recovers second-order terms that are neglected in the FOA.
4. Continuous formulation of the multilayer
$^{\ast }$
model
Since the multilayer model relies on finite differences to approximate vertical derivatives and to reconstruct the interlayer tractions, a question naturally arises from such an approach: can the multilayer model be derived from a vertical semi-discretisation of a 3D model with an extruded mesh? The model derivation of § 2 and the analysis of § 3 show that if this is the case, then, this 3D model is neither the Stokes nor the FOA model. On the one hand, the multilayer model ‘is not included’ in the FOA because of the redefinition of the tractions at second-order. Indeed, the ‘infinite parallel-sided slab’ of § 3 provides an example of where the multilayer is exact for Stokes (up to quadrature errors) while the FOA is not. On the other hand, the multilayer models are clearly not a semi-discretisation of the Stokes model.
In contrast, the multilayer
$^{\ast }$
model does not involve any second-order terms. We can now show it is a vertical semi-discretisation of a 3D model deriving from the FOA. Indeed, one rewrites (in the flowline setting for simplicity) the
$l$
th equation of the multilayer
$^{\ast }$
system:
$$\begin{eqnarray}\displaystyle & & \displaystyle 2h^{l}A^{-1/n}\partial _{x}\left(\left|\partial _{x}u_{x}^{l}\right|^{1/n-1}\partial _{x}u_{x}^{l}\right)-A^{-1/n}\left|\frac{u_{x}^{l}-u_{x}^{l+1}}{h^{l}+h^{l+1}}\right|^{1/n-1}\left(\frac{u_{x}^{l}-u_{x}^{l+1}}{h^{l}+h^{l+1}}\right)\nonumber\\ \displaystyle & & \displaystyle \quad -\,A^{-1/n}\left|\frac{u_{x}^{l}-u_{x}^{l-1}}{h^{l}+h^{l-1}}\right|^{1/n-1}\left(\frac{u_{x}^{l}-u_{x}^{l-1}}{h^{l}+h^{l-1}}\right)={\it\rho}gh^{l}\partial _{x}\overline{s},\end{eqnarray}$$
where
$h^{l}$
and
$\partial _{x}$
can be interchanged in the first term of (4.1) since
$O({\it\delta})$
terms are neglected in the multilayer
$^{\ast }$
model. It appears that (4.1) can be seen as a vertical semi-discretisation by finite differences of the vertically integrated version of
$$\begin{eqnarray}2A^{-1/n}\partial _{x}\left(|\partial _{x}u_{x}|^{1/n-1}\partial _{x}u_{x}\right)+A^{-1/n}\partial _{z}\left(|\partial _{z}u_{x}/2|^{1/n-1}(\partial _{z}u_{x}/2)\right)={\it\rho}g\partial _{x}\overline{s},\end{eqnarray}$$
where
$\partial _{x}u_{x}^{l}$
approximates
$\partial _{x}u_{x}$
. One can compare (4.2) with the FOA equation (2.29), i.e.
$$\begin{eqnarray}\displaystyle & & \displaystyle 2A^{-1/n}\partial _{x}\left(\left[\left(\partial _{x}u_{x}\right)^{2}+\left(\frac{\partial _{z}u_{x}}{2}\right)^{2}\right]^{(1/n-1)/2}\left(\partial _{x}u_{x}\right)\right)\nonumber\\ \displaystyle & & \displaystyle \quad +\,A^{-1/n}\partial _{z}\left(\left[\left(\partial _{x}u_{x}\right)^{2}+\left(\frac{\partial _{z}u_{x}}{2}\right)^{2}\right]^{(1/n-1)/2}\left(\frac{\partial _{z}u_{x}}{2}\right)\right)={\it\rho}g\partial _{x}\overline{s}.\end{eqnarray}$$
As a matter of fact, the viscosity which contains both vertical and horizontal derivatives in the FOA is simplified by dropping the crossing terms which multiply vertical and horizontal derivatives in the multilayer
$^{\ast }$
model. As a consequence, the multilayer
$^{\ast }$
model is expected to be the most accurate when one of the two types of stresses (horizontal stress or vertical shear) dominates (i.e. for very small
${\it\lambda}$
or very high
${\it\lambda}$
, see Schoof & Hindmarsh (Reference Schoof and Hindmarsh2010)). In contrast, the model is expected to less accurate when those two types of stress are in the same order of magnitude, see § 5. Interestingly, the FOA equation (4.3) and the multilayer
$^{\ast }$
(4.2) coincide in the case of a Newtonian fluid (
$n=1$
) since the viscosity is constant.
5. ISMIP-HOM experiments results
ISMIP-HOM (Pattyn et al. Reference Pattyn, Perichon, Aschwanden, Breuer, de Smedt, Gagliardini, Gudmundsson, Hindmarsh, Hubbard, Johnson, Kleiner, Konovalov, Martin, Payne, Pollard, Price, Rückamp, Saito, Souček, Sugiyama and Zwinger2008) experiments consist of modelling exercises based on various idealised ice geometries and boundary conditions in order to generate different types of ice flows, which can be met in real glacier modelling. These benchmark experiments have become popular for conducting comparative studies of the performances of ice flow models, see e.g. Gagliardini & Zwinger (Reference Gagliardini and Zwinger2008) and Goldberg (Reference Goldberg2011). In this section, numerical multilayer solutions computed for the diagnostic flowline ISMIP-HOM experiments (B, D and E) are compared with those obtained using the FOA and the Stokes models.
In all ISMIP-HOM experiments, the multilayer solutions were computed using a column-wise extension of the Newton multigrid solver presented in Jouvet & Gräser (Reference Jouvet and Gräser2013) for the SSA. According to the experiment, three types of multilayer vertical splitting (2.41) were used: (a) the ‘uniform’ splitting defined by
$h^{l}=h/L$
, (b) the ‘exact’ splitting defined by the streamlines of the Stokes solution (Gagliardini & Zwinger Reference Gagliardini and Zwinger2008) and (c) the ‘bed-aligned’ splitting defined by a pile of layers made of constant-thickness lower layers and degenerating, partly zero-thickness upper layers. For the sake of convenience, the ‘
$L$
-
$\text{layer}^{\,s}$
’ denotes the multilayer with
$L$
layers and the splitting of type
$s\in \{u,e,b\}$
, where ‘
$u$
’, ‘
$e$
’ and ‘
$b$
’ denote ‘uniform’, ‘exact’ and ‘bed-aligned’, respectively. In addition, a star is added to the exponent, e.g. the 16-
$\text{layer}^{\,u,\ast }$
, when neglecting
$O({\it\delta}^{2})$
terms in the interlayer terms (2.59), i.e. setting
${\it\alpha}^{l}=1$
and
$M_{ij}^{l}={\it\delta}_{ij}$
instead of (2.47) and (2.61).
For each experiment, the horizontal segment
${\it\Omega}$
was divided into 256 equal sized segments to generate a 1D mesh. On the one hand, the resulting mesh was used to compute the multilayer solutions. On the other hand, a triangular 2D mesh was built by extruding uniformly vertical layers of the 1D mesh between the lower and the upper surfaces, and splitting each rectangle into two triangles. This 2D mesh was used to compute the FOA solution using a nonlinear Gauß–Seidel solver. Lastly, the Stokes solutions published by Gagliardini & Zwinger (Reference Gagliardini and Zwinger2008) (or the model ‘oga1’ from Pattyn et al. (Reference Pattyn, Perichon, Aschwanden, Breuer, de Smedt, Gagliardini, Gudmundsson, Hindmarsh, Hubbard, Johnson, Kleiner, Konovalov, Martin, Payne, Pollard, Price, Rückamp, Saito, Souček, Sugiyama and Zwinger2008)) are used for comparison purposes. The numerical convergence of the multilayer and the FOA solutions was assessed by examining the discrepancy between this solution and that obtained by doubling the horizontal resolution and the number of layers. For further confidence, the gradient and the viscosity fields of the multilayer solutions were also checked. As result, all fields were always found convergent when refining the horizontal mesh, the multilayer vertical splitting or both simultaneously. For all experiments, the following physical parameters were used:
${\it\rho}=910~\text{kg}~\text{m}^{-3}$
,
$n=3$
,
$A=3.17\times 10^{-24}~\text{Pa}^{-3}~\text{s}^{-1}$
and
$g=9.81~\text{m}~\text{s}^{-2}$
.
In experiment B, the geometry is defined by
$$\begin{eqnarray}\displaystyle & \overline{s}(x)=-x\tan (0.5^{\circ }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \underline{s}(x)=\overline{s}(x)-1000+500\sin (2{\rm\pi}x/\bar{L}), & \displaystyle\end{eqnarray}$$
$x\in {\it\Omega}=[0,\bar{L}]$
, the no-slip condition is prescribed on the bedrock, i.e.
${\it\Omega}_{f}={\it\Omega}$
, periodic boundary conditions connect the left- and right-hand sides of
${\it\Omega}$
, and the stress-free condition is prescribed on the top surface, see Pattyn et al. (Reference Pattyn, Perichon, Aschwanden, Breuer, de Smedt, Gagliardini, Gudmundsson, Hindmarsh, Hubbard, Johnson, Kleiner, Konovalov, Martin, Payne, Pollard, Price, Rückamp, Saito, Souček, Sugiyama and Zwinger2008) for further details. Figure 4(a,c,e) displays the surface horizontal velocities of the 16-
$\text{layer}^{\,u,\ast }$
, the 16-
$\text{layer}^{u}$
, the 16-
$\text{layer}^{e}$
, the 16-
$\text{layer}^{b}$
, the FOA and the Stokes models for
$\bar{L}=10$
, 40 and 160 km. To sum up, figure 5 displays the 16-
$\text{layer}^{u}$
solution over the entire domain for
$\bar{L}=10$
and 40 km.
ISMIP-HOM experiment B (a,c,e) and D (b,d, f) surface horizontal velocities of the 16-
$\text{layer}^{\,u,\ast }$
, the 16-
$\text{layer}^{u}$
, the 16-
$\text{layer}^{b}$
(only for experiment B), the 16-
$\text{layer}^{e}$
(only for experiment B), the FOA and the Stokes models for
$\bar{L}=10~\text{km}$
(a,b), 40 km (c,d) and 160 km (e,f), respectively.
ISMIP-HOM experiment B results. (a) Streamlines of the Stokes velocity field for
$\bar{L}$
= 10; (b,c) horizontal velocity field of the 16-
$\text{layer}^{u}$
model for
$\bar{L}=10$
and 40 km, respectively. For the sake of convenience, the domain was stretched vertically.
Figure 4(a,c,e) indicates that the simplified 16-
$\text{layer}^{\,u,\ast }$
model leads to greater disagreements with the FOA and the Stokes solutions than all 16-layer models at small wavelengths (i.e. when the interlayer slopes are high) while no differences between the 16-
$\text{layer}^{u}$
and the 16-
$\text{layer}^{\,u,\ast }$
solutions can be observed for larger wavelengths (i.e. when the interlayer slopes are negligible). Recall that, in contrast with the other 16-layer models, the 16-
$\text{layer}^{\,u,\ast }$
neglects
$O({\it\delta}^{2})$
terms in the interlayer traction. This proves that these terms are relevant and must be kept, in particular, where the geometry is steep. As seen in § 4, the discrepancy between the 16-
$\text{layer}^{\,u,\ast }$
and the FOA solutions is the direct consequence of the decoupling between (4.2) and (4.3) which differentiate both models. In contrast with the 16-
$\text{layer}^{\,u,\ast }$
, the other 16-layer solutions rely on the assumption that the layers are aligned with the streamlines. To get rid of the effects of this assumption, one looks first at the 16-
$\text{layer}^{e}$
solution which uses the ‘exact alignment’, i.e. the one drawn by streamlines of the Stokes model, see figure 5(a). Interestingly, when
$\bar{L}=10~\text{km}$
, the 16-
$\text{layer}^{e}$
solution matches the Stokes one well in the first half of the domain, but differs by approximately 20 % in the second part, where a depression occurs in the bedrock. However, this disagreement decreases (without vanishing) for the higher wavelengths (
$\bar{L}=40$
and
$\bar{L}=160~\text{km}$
). The agreement of the 16-
$\text{layer}^{e}$
and the Stokes solutions where the FOA deviates in the first half of the domain cannot be interpreted as a superiority of the multilayer model over the FOA, as observed with the ‘infinite parallel-sided slab’ in § 3. Indeed, the multilayer and Stokes models differ by more than one feature such that several sources of error might cancel and lead accidentally to a good agreement between both solutions. Also, this agreement is no longer true in the second half of the domain, where the depression area occurs. Figure 5(a) shows that the streamlines of the Stokes velocity field used to build the ‘exact alignment’ are significantly curved at the top of the depression. In contrast, the ‘uniform alignment’ which is used to compute the 16-
$\text{layer}^{\,u}$
solution leads to nearly flat layers close to the top. As a consequence, using the ‘uniform alignment’ instead of the ‘exact one’ is expected to introduce some artificial resistance and to slow down the flow. Figure 4 confirms that the magnitude of the 16-
$\text{layer}^{\,u}$
solution is damped compared to the reference 16-
$\text{layer}^{\,e}$
one in the depression area. However, the impact of this incorrect alignment of layers is negligible for the higher wavelengths
$L$
. Lastly, the bed-aligned multilayer splitting, which builds a pile of layers made of constant-thickness lower layers and zero-thickness upper layers, is considered in order to compute the last multilayer solution. This empirical alignment is motivated by the fact that the streamlines seem to be approximatively aligned to bed, as shown by figure 5(a). As a result, figure 4 indicates that the 16-
$\text{layer}^{b}$
solution obtained with bed-aligned division nearly fits the one obtained with the exact alignment. As a conclusion, the bed-aligned splitting proves to be a better empirical choice than the uniform splitting in the case of experiment B.
In experiment D, the geometry is defined by
$$\begin{eqnarray}\displaystyle & \overline{s}(x)=-x\tan (0.1^{\circ }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \underline{s}(x)=\overline{s}(x)-1000, & \displaystyle\end{eqnarray}$$
$x\in {\it\Omega}=[0,\bar{L}]$
, the slip condition is prescribed everywhere on the bedrock, i.e.
${\it\Omega}_{m}={\it\Omega}$
, with
$m=1$
and 
$$\begin{eqnarray}C(x)=1000(1+\sin (2{\rm\pi}x/\bar{L})),\end{eqnarray}$$
periodic boundary conditions connect the left- and right-hand sides of
${\it\Omega}$
, and the stress-free condition is prescribed on the top surface, see Pattyn et al. (Reference Pattyn, Perichon, Aschwanden, Breuer, de Smedt, Gagliardini, Gudmundsson, Hindmarsh, Hubbard, Johnson, Kleiner, Konovalov, Martin, Payne, Pollard, Price, Rückamp, Saito, Souček, Sugiyama and Zwinger2008) for further details. Again, a uniform multilayer splitting, which here matches the bed-aligned one, was chosen in order to calculate the multilayer solutions. Figure 4(b,d, f) displays the surface horizontal velocities of the
$L$
-
$\text{layer}^{\,u,\ast }$
, the
$L$
-
$\text{layer}^{u}$
, the FOA and the Stokes models for the largest
$L$
and
$\bar{L}=10$
, 40 and 160 km. After that, figure 6 displays the multilayer solution over the entire domain for
$\bar{L}=10$
and 40 km.
ISMIP-HOM experiment D results for (a)
$\bar{L}=10$
km and (b)
$\bar{L}=40$
km. Horizontal velocity field of the
$\text{multilayer}^{u}$
model. For the sake of convenience, the domain was stretched vertically.
In contrast to experiment B, the convergence with respect to the number of layers
$L$
was found (not shown) slower for low-wavelength
$\bar{L}$
. Indeed, figure 6 shows that the solution strongly varies with
$z$
for
$\bar{L}=10~\text{km}$
, but is nearly constant with respect to
$z$
for higher values of
$\bar{L}$
. As a consequence, the case
$\bar{L}=10$
naturally requires further layers for a given convergence threshold. Moreover, unlike experiment B, the simplified 16-
$\text{layer}^{\,u,\ast }$
and the 16-
$\text{layer}^{u}$
show indistinguishable solutions. This is due to the very slight slope, which renders the
$O({\it\delta}^{2})$
terms inactive in the interlayer tractions. For all wavelengths, the multilayer solutions match the FOA and the Stokes solutions well, and it can be verified (not shown) that opting for the ‘exact splitting’ instead of the uniform splitting (like previously in experiment B) hardly improves the solution at all.
Experiment E is conducted along the central flowline of a 5-km-long temperate glacier in Switzerland (Haut Arolla glacier). Model inputs (longitudinal surface and bedrock profiles) are given (Pattyn et al.
Reference Pattyn, Perichon, Aschwanden, Breuer, de Smedt, Gagliardini, Gudmundsson, Hindmarsh, Hubbard, Johnson, Kleiner, Konovalov, Martin, Payne, Pollard, Price, Rückamp, Saito, Souček, Sugiyama and Zwinger2008). In the first experiment (denoted E1), the no-slip condition is prescribed everywhere on the bedrock, i.e.
${\it\Omega}_{f}={\it\Omega}$
, while in the second experiment (denoted E2), the no-slip condition applies everywhere except in the zone defined by
${\it\Omega}_{m}=\{x,\;2200~\text{m}\leqslant x\leqslant 2500~\text{m}\}$
, where the perfect slip condition with
$C=0$
is prescribed. The stress-free condition is used on the top surface. Again, a uniform vertical splitting was chosen first in order to calculate the multilayer solutions, see figure 7. Figure 8 displays the surface horizontal velocities of the 16-
$\text{layer}^{\,u,\ast }$
, the 16-
$\text{layer}^{u}$
, the 16-
$\text{layer}^{e}$
, the 16-
$\text{layer}^{b}$
, the FOA and the Stokes models for E1 and E2. At the end, figure 9 displays the 16-
$\text{layer}^{u}$
and Stokes solutions over the entire domain for experiments E1 and E2.
Multilayer splitting (2.41): uniform (‘
$u$
’), bed-aligned (‘
$b$
’), ‘exact’ for E1 and for E2 in ISMIP-HOM experiment E. For the sake of convenience, the domain was stretched vertically. (a) Empirical alignment. (b) Exact alignment.
ISMIP-HOM experiment E1 (a) and E2 (b) surface horizontal velocities of the 16-
$\text{layer}^{\,u,\ast }$
, the 16-
$\text{layer}^{\,u}$
, the 16-
$\text{layer}^{\,b}$
, the 16-
$\text{layer}^{\,e}$
, the FOA and the Stokes models.
ISMIP-HOM Experiments E1 (a) and E2 (b) Horizontal velocity fields obtained using the 16-
$\text{layer}^{u}$
model (top) and the Stokes model (bottom). For the sake of convenience, the domain was stretched vertically.
Figure 8 indicates that the results of 16-
$\text{layer}^{\,u,\ast }$
model differs from those of the other 16-layer models, in particular, it overestimates the solution in the first area
$\{x,\;0\leqslant x\leqslant 1~\text{km}\}$
for both experiments E1 and E2. As in experiment B, this confirms the relevance of the
$O({\it\delta}^{2})$
terms in the interlayer traction. In order to find out the causes of the misfit in the central part, the 16-
$\text{layer}^{e}$
solution is computed with the ‘exact alignment’ of the layers, which is given by the streamlines of the Stokes solutions, see figures 7(b) and 8. Interestingly, this solution matches the Stokes solution very well in magnitude, but fails to reproduce the short-wavelength variations of the central area. Regarding the conclusions of experiment B, the slight mismatch of the multilayer in experiment E1 likely results from the small depression in the bedrock of the same area, see figure 9. Concerning experiment E2, it is interesting to notice that only the Stokes solution shows the influence of the abrupt changes in the basal conditions on the surface. In contrast, all of the other solutions (including the FOA) show smooth surface velocities. In addition, the multilayer solution obtained using the ‘exact alignment’ fits the Stokes solution better than the FOA one for experiment E2. As in experiment B, using a ‘uniform alignment’ instead of the ‘exact one’ damps and, then, deteriorates the 16-layer solution for both experiments, as shown in figure 8. In contrast, the bed-aligned multilayer splitting proves to be a better empirical choice for experiment E1 (as in experiment B) even if a gap between the 16-
$\text{layer}^{e}$
and the 16-
$\text{layer}^{b}$
solutions is still visible. This is due to the fact that the streamlines are approximatively aligned to the bed in the central area, see figure 7. However, using such an alignment improves much less the solution in experiment E2 since the basal boundary conditions disturb the streamlines such that they are no longer aligned to the bed topography.
6. Conclusions and perspectives
As is often the case in geophysics, ice flow models are derived by following the two-step procedure: the mechanical modelling strictly precedes the numerical modelling. It is sometimes useful to go back to the physics after the discretisation, e.g. when preconditioning the linear systems (Brown, Smith & Ahmadia Reference Brown, Smith and Ahmadia2013). As in Audusse et al. (Reference Audusse, Bristeau, Perthame and Sainte-Marie2011), the common order of modelling was partly inverted here since the vertical discretisation comes first followed by the mechanical modelling. Using a multilayer splitting of the ice thickness, a new hybrid ice flow model generalising the SSA was derived by depth-integrating the 3D FOA model locally. The resulting model consists of a system of 2D equations similar to the SSA ones. Advantageously, the vertical shear stress components, which are ignored in the SSA, are recovered in the multilayer model through the reconstructed tractions. Keeping only the normal shear stress components to describe the interlayer sliding, the tractions can be simply redefined at zeroth order for the slope of layers. The resulting model (called multilayer
$^{\ast }$
) can be seen as the vertical semi-discretisation of a decoupled version of the FOA (§ 4). However, one can redefine the interlayer tractions at second order if the layer boundaries are in alignment with the streamlines. Based on this more general redefinition, the multilayer solution even equals the Stokes one in the ‘infinite parallel-sided slab’ setting of § 3. When running the model for prognostic flowline ISMIP-HOM benchmark experiments, the multilayer solutions based on vertical splittings which are empirically aligned to the streamlines show good agreement with the higher-order solutions if no severe depression occurs in the bedrock.
As a follow up of this paper, several aspects of the multilayer approach introduced in this paper are developed in ongoing work. First, methods for upgrading any SSA solver to the multilayer system will be developed and their numerical performances will be assessed. In particular, the computational performance gain when using the multilayer approach instead of any 3D model will be quantified. In the present paper, the multilayer solutions were restricted to simple flowline diagnostic experiments. To complete this work, model comparisons for 3D prognostic real glaciers including a shelf part will be conducted. Another crucial aspect of the multilayer approach is the choice of the vertical subdivision and its impact on the results. This splitting must be done in such a way that the ice flows and the layers are approximately aligned. Although the ‘bed-aligned splitting’ proved to be a good empirical choice in the no-sliding case of ISMIP-HOM experiments, a gap was still visible with the ‘exact alignment’. Moreover, the quality of this approximation was called into question when prescribing basal sliding conditions. Thus, the question of the optimal choice of the vertical splitting and its impact on the results will be investigated in a future study.
Acknowledgements
The author wishes to acknowledge Ed Bueler, H. Blatter, K. Hutter and three anonymous referees for their detailed and helpful comments on the manuscript, C. Gräser for computational assistance, S. Braun-Clarke for proofreading the English, and R. Kornhuber for support. The Elmer/ice code was used to compute the streamlines in ISMIP-HOM experiments. Elmer’s team is acknowledged for its support during an Elmer course. The author was supported by the Deutsche Forschungsgemeinschaft (project KL 1806 5-1).
Appendix A. Second-order interlayer tractions
In § 2.4, the interlayer tractions (2.57) were derived at zeroth order for the slope of layers. Here the second-order tractions (2.59) are derived from the three following hypotheses. First, the layer-parallel stresses in the local frame are set at zero:
$$\begin{eqnarray}\unicode[STIX]{x1D60B}_{ij}T_{ix}T_{jx}=0,\quad \unicode[STIX]{x1D60B}_{ij}T_{ix}T_{jy}=0,\quad \unicode[STIX]{x1D60B}_{ij}T_{iy}T_{jy}=0,\end{eqnarray}$$
where
$\{T_{ix}\}_{i}$
and
$\{T_{jy}\}_{j}$
are two unit independent vectors orthogonal to
$\boldsymbol{n}^{l}$
such that they all shape an orthogonal basis. Second, the multilayer vertical splitting (2.41) is chosen such that ‘the layers are aligned with the three-dimensional direction of the flow’:
$$\begin{eqnarray}u_{j}n_{j}^{l}=0.\end{eqnarray}$$
Third, ‘the layers are sufficiently thin’ such that
${\it\epsilon}/L$
, which corresponds to the typical aspect ratio of layers, is small compared to
${\it\delta}$
, which corresponds to the typical slope of layers. As a consequence, high-order terms
$O({\it\delta}{\it\varepsilon}/L)$
(called later
${\mathcal{R}}^{l}$
,
$\bar{{\mathcal{R}}}^{l}$
,
$\bar{\bar{{\mathcal{R}}}}^{l}$
and
$\bar{\bar{\bar{{\mathcal{R}}}}}^{l}$
) are neglected, but second-order terms
$O({\it\delta}^{2})$
are retained.
Since
$|\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ji}/2|$
is invariant under changes of basis, then
$$\begin{eqnarray}\left|\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ji}/2\right|=\left|(\unicode[STIX]{x1D60B}_{mn}\unicode[STIX]{x1D61B}_{mi}\unicode[STIX]{x1D61B}_{nj})(\unicode[STIX]{x1D60B}_{mn}\unicode[STIX]{x1D61B}_{mj}\unicode[STIX]{x1D61B}_{ni})/2\right|\!,\end{eqnarray}$$
where the orthogonal matrix
$\unicode[STIX]{x1D61B}_{ij}$
is built with the components of the orthogonal basis composed by
$\{\unicode[STIX]{x1D61B}_{ix}\}_{i}$
,
$\{\unicode[STIX]{x1D61B}_{iy}\}_{i}$
and
$\{\unicode[STIX]{x1D61B}_{iz}:=n_{i}^{l}\}_{i}$
. In addition to (A 1),
$\unicode[STIX]{x1D60B}_{mn}\unicode[STIX]{x1D61B}_{mz}\unicode[STIX]{x1D61B}_{nz}=0$
since the trace operator is also invariant under changes of basis. Due to those simplifications, the symmetric matrix
$D_{mn}T_{mi}T_{nj}$
contains only four non-zero entries such that (A 3) can be rewritten as
$$\begin{eqnarray}\left|\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ji}/2\right|=(\unicode[STIX]{x1D60B}_{mn}\unicode[STIX]{x1D61B}_{mx}\unicode[STIX]{x1D61B}_{nz})^{2}+(\unicode[STIX]{x1D60B}_{mn}\unicode[STIX]{x1D61B}_{my}\unicode[STIX]{x1D61B}_{nz})^{2}=(\unicode[STIX]{x1D60B}_{mn}\unicode[STIX]{x1D61B}_{mx}n_{n}^{l})^{2}+(\unicode[STIX]{x1D60B}_{mn}\unicode[STIX]{x1D61B}_{my}n_{n}^{l})^{2}.\end{eqnarray}$$
From the invariance of the Euclidean distance under an orthogonal transformation, it follows that
$$\begin{eqnarray}\left|\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ji}/2\right|=(\unicode[STIX]{x1D60B}_{ij}n_{j}^{l})(\unicode[STIX]{x1D60B}_{ij}n_{j}^{l}).\end{eqnarray}$$
As a consequence, equation (2.54) becomes
$$\begin{eqnarray}S_{k}^{l}=A^{-1/n}{\it\alpha}^{l}\left|(\unicode[STIX]{x1D60B}_{ij}n_{j}^{l})(\unicode[STIX]{x1D60B}_{ij}n_{j}^{l})\right|^{(1/n-1)/2}\unicode[STIX]{x1D60B}_{ij}n_{j}^{l}t_{i}^{k,l}.\end{eqnarray}$$
Using definition (2.9) of
$\unicode[STIX]{x1D60B}_{ij}$
and a differentiation rule leads to
$$\begin{eqnarray}\unicode[STIX]{x1D60B}_{ij}n_{j}^{l}={\textstyle \frac{1}{2}}(\partial _{j}u_{i})n_{j}^{l}+{\textstyle \frac{1}{2}}(\partial _{i}u_{j})n_{j}^{l}={\textstyle \frac{1}{2}}(\partial _{j}u_{i})n_{j}^{l}+{\textstyle \frac{1}{2}}(\partial _{i}(u_{j}n_{j}^{l}))-{\textstyle \frac{1}{2}}(\partial _{i}n_{j}^{l})u_{j}.\end{eqnarray}$$
By (A 2), the middle term
$(\partial _{i}(u_{j}n_{j}^{l}))/2$
is zero. The other terms
$(\partial _{j}u_{i})n_{j}^{l}/2$
and
$(\partial _{i}n_{j}^{l})u_{j}/2$
are now treated in turn.
First,
$(\partial _{j}u_{i})n_{j}^{l}$
is the derivative of
$\boldsymbol{u}$
in the direction
$\boldsymbol{n}^{l}$
:
$$\begin{eqnarray}{\textstyle \frac{1}{2}}(\partial _{j}u_{i})n_{j}^{l}={\textstyle \frac{1}{2}}\partial _{\tilde{z}}u_{i},\end{eqnarray}$$
where
$\tilde{z}$
is the local variable defined by
$\tilde{z}=n_{x}^{l}x+n_{y}^{l}y+n_{z}^{l}z$
. In (A 8), the derivative with respect to
$\tilde{z}$
(according to the direction orthogonal to the interface) is approximated by the finite difference:
$$\begin{eqnarray}\frac{\partial _{\tilde{z}}u_{i}}{2}=\frac{u_{i}^{l+1}-u_{i}^{l}}{\tilde{h}^{l+1}+\tilde{h}^{l}},\quad i\in \{x,y,z\}.\end{eqnarray}$$
It should be stressed that (A 9) depends on the local coordinate
$\tilde{z}$
and not on
$z$
. However, the Taylor expansion of
$u_{k}^{l}$
at
$\tilde{z}$
can be written
$$\begin{eqnarray}u_{i}^{l}(\tilde{z})=u_{i}^{l}(x,y,z)+{\mathcal{R}}^{l},\quad i\in \{x,y,z\},\end{eqnarray}$$
where
$$\begin{eqnarray}{\mathcal{R}}^{l}=O\left(\frac{(\partial _{j}s^{l}\partial _{j}u_{i}^{l})h^{l}}{{\it\alpha}^{l}}\right)=O\left([u_{i}^{l}]\frac{{\it\delta}{\it\epsilon}}{L}\right)\!,\end{eqnarray}$$
which can be neglected. As a consequence,
$u_{i}^{l}$
and
$u_{i}^{l+1}$
are considered as locally constant in the plane orthogonal to
$\boldsymbol{n}^{l}$
, so (A 9) holds in the primary variables
$(x,y,z)$
too. In addition, calling
$\tilde{h}^{l}$
the thickness function in the local frame, the Taylor expansion of
$\tilde{h}^{l}$
and
$\tilde{h}^{l+1}$
at
$\tilde{z}$
yields
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{h}^{l}=\frac{h^{l}}{{\it\alpha}^{l}}+\bar{{\mathcal{R}}}^{l},\quad \text{where}~\bar{{\mathcal{R}}}^{l}=O\left(\frac{(\partial _{j}s^{l}\partial _{j}h^{l})h^{l}}{{\it\alpha}^{l}}\right)=O\left([h^{l}]\frac{{\it\delta}{\it\epsilon}}{L}\right)\!, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{h}^{l+1}=\frac{h^{l+1}}{{\it\alpha}^{l}}+\bar{\bar{{\mathcal{R}}}}^{l},\quad \text{where}~\bar{\bar{{\mathcal{R}}}}^{l}=O\left(-\frac{(\partial _{j}s^{l}\partial _{j}h^{l+1})h^{l+1}}{{\it\alpha}^{l}}\right)=O\left([h^{l+1}]\frac{{\it\delta}{\it\epsilon}}{L}\right)\!, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\frac{1}{2}(\partial _{j}u_{i})n_{j}^{l}={\it\alpha}^{l}\left(\frac{u_{i}^{l+1}-u_{i}^{l}}{h^{l}+h^{l+1}}\right).\end{eqnarray}$$
Representation of the global frame (left) and the local frame (right).
Second, computing the term
$(\partial _{i}n_{j}^{l})u_{j}/2$
with the definition of
$\boldsymbol{n}^{l}$
given by (2.46), and using (A 2), leads to
$$\begin{eqnarray}\frac{1}{2}(\partial _{i}n_{j}^{l})u_{j}=\frac{1}{2{\it\alpha}^{l}}(\partial _{ij}s^{l})u_{j},\quad i\in \{x,y\},\quad \frac{1}{2}(\partial _{z}n_{j}^{l})u_{j}=0,\end{eqnarray}$$
where
$(\partial _{ij}s^{l})$
is the Hessian matrix of
$s^{l}$
in the horizontal plane. From (A 14) and (A 15), it follows that
$$\begin{eqnarray}\frac{1}{2}(\partial _{j}u_{i})n_{j}^{l}-\frac{1}{2}(\partial _{i}n_{j}^{l})u_{j}=\left(\frac{{\it\alpha}^{l}}{h^{l}+h^{l+1}}\right)\left((u_{i}^{l+1}-u_{i}^{l})+\bar{\bar{\bar{{\mathcal{R}}}}}^{l}\right)\end{eqnarray}$$
where
$$\begin{eqnarray}\bar{\bar{\bar{{\mathcal{R}}}}}^{l}=O\left(-\left(\frac{h^{l}+h^{l+1}}{2({\it\alpha}^{l})^{2}}\right)(\partial _{ij}s^{l})u_{j}\right)=O\left([u_{i}]\frac{{\it\delta}{\it\epsilon}}{L}\right)\!,\end{eqnarray}$$
which can be neglected.
Finally, using (A 14) and (A 16), equation (A 6) becomes
$$\begin{eqnarray}S_{k}^{l}=A^{-1/n}{\it\alpha}^{l}\left|\left({\it\alpha}^{l}\frac{u_{i}^{l+1}-u_{i}^{l}}{h^{l}+h^{l+1}}\right)\left({\it\alpha}^{l}\frac{u_{i}^{l+1}-u_{i}^{l}}{h^{l}+h^{l+1}}\right)\right|^{(1/n-1)/2}\left({\it\alpha}^{l}\frac{u_{i}^{l+1}-u_{i}^{l}}{h^{l}+h^{l+1}}\right)t_{i}^{l,k}\end{eqnarray}$$
for
$i\in \{x,y,z\}$
and
$k\in \{x,y\}$
. In addition, equation (A 2) implies that the third component
$u_{z}^{l+1}-u_{z}^{l}$
can be eliminated such that (A 18) can be rewritten as
$$\begin{eqnarray}S_{k}^{l}=A^{-1/n}{\it\alpha}^{l}\left|M_{ij}^{l}\left({\it\alpha}^{l}\frac{u_{i}^{l+1}-u_{i}^{l}}{h^{l}+h^{l+1}}\right)\left({\it\alpha}^{l}\frac{u_{j}^{l+1}-u_{j}^{l}}{h^{l}+h^{l+1}}\right)\right|^{(1/n-1)/2}M_{ik}^{l}\left({\it\alpha}^{l}\frac{u_{i}^{l+1}-u_{i}^{l}}{h^{l}+h^{l+1}}\right)\!,\end{eqnarray}$$
for
$i,j,k\in \{x,y\}$
, where
$M_{ij}^{l}={\it\delta}_{ij}+(\partial _{i}s^{l})(\partial _{j}s^{l})$
. Finally, the interlayer tractions (A 19) correspond to (2.59).
