2.1. Governing equations

2.1.1. Multi-domain formulation

Let us denote by $({h} _{g}, u_{g})$ and $(h_{f}, u_{f})$ the values taken by the functions $(h, u)$ on the grounded portion $\varOmega _{g}$ and the floating portion $\varOmega _{f}$ of the domain $\varOmega$, respectively. With these notations, the governing equations read as follows in the grounded portion:

(2.1a,b) \begin{equation} \left\{ \begin{gathered} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\dfrac{\mathrm{d}}{\mathrm{d}\kern 0.06em x}(u_{g} {h}_{g}) = a, &\quad \text{in } \varOmega_{g},\\ 2 A^{-{1}/{n}} \dfrac{\mathrm{d}}{\mathrm{d}\kern 0.06em x}\left( {h}_{g} \left\vert \dfrac{\mathrm{d} u_{g}}{\mathrm{d}\kern 0.06em x} \right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} u_{g}}{\mathrm{d}\kern 0.06em x}\right) - \tau_{b} \phantom{=0,} &\nonumber\\ \quad- \varLambda {h}_{g}\vert u_{g} \vert^{m-1}u_{g} = \rho g {h}_{g} \dfrac{\mathrm{d}}{\mathrm{d}\kern 0.06em x}(b + {h}_{g}), &\quad \text{in } \varOmega_{g}. \end{gathered}\right.\end{equation}

Equation (2.1a) is a mass-conservation equation, stating that the flux variation of the ice flow must be exactly compensated by the net mass accumulation rate $a$. Equation (2.1b) is a momentum-conservation equation and establishes a balance between the divergence of membrane stress, the friction stress, the lateral-drag stress and the gravitational stress. The factor $A$ and the exponent $n$ are ice viscosity parameters associated with the Glen flow law (usually, $n=3$), $\varLambda$ and m are lateral-drag coefficients, $\rho$ is the ice density, $\rho _{w}$ is the water density, $g$ is the gravitational acceleration and $b=b(x)$ is the prescribed elevation of the underlying bedrock. The models used for the friction stress $\tau _{b}=\tau _{b}(h, u)$ are described in the next section. While we do not explicitly consider lateral drag in the present study, we do include it in the problem formulation, as it allows for an easier comparison with the other results from the literature.

In the floating portion, friction with the ocean and the air is neglected, leading to the following:

(2.2a,b) \begin{equation} \left\{ \begin{gathered} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\dfrac{\mathrm{d}}{\mathrm{d}\kern 0.06em x}(u_{f} h_{f}) = a, &\quad \text{in } \varOmega_{f},\\ 2 A^{-{1}/{n}} \dfrac{\mathrm{d}}{\mathrm{d}\kern 0.06em x}\left( h_{f} \left\vert \dfrac{\mathrm{d} u_{f}}{\mathrm{d}\kern 0.06em x} \right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} u_{f}}{\mathrm{d}\kern 0.06em x}\right) -\varLambda h_{f}\vert u_{f} \vert^{m-1}u_{f} \phantom{=0,} & \nonumber\\ \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad= \rho\left(1 - \dfrac{\rho}{\rho_{w}}\right) g h_{f} \dfrac{\mathrm{d} h_{f}}{\mathrm{d}\kern 0.06em x},&\quad \text{in } \varOmega_{f}. \end{gathered}\right.\end{equation}

Finally, continuity conditions are added at the interface between the regions:

(2.3) \begin{align} {h}_{g} =h_{f}, \quad u_{g}= u_{f}, \quad 2 A^{-{1}/{n}} {h}_{g} \left\vert \dfrac{\mathrm{d} u_{g}}{\mathrm{d}\kern 0.06em x} \right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} u_{g}}{\mathrm{d}\kern 0.06em x}= 2 A^{-{1}/{n}} h_{f} \left\vert \dfrac{\mathrm{d} u_{f}}{\mathrm{d}\kern 0.06em x} \right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} u_{f}}{\mathrm{d}\kern 0.06em x}\quad \text{on } \varSigma. \end{align}

The portions $\varOmega _{g}$ and $\varOmega _{f}$ and their interface $\varSigma$ are defined by a flotation condition:

(2.4a–4) \begin{equation} \left\{ \begin{gathered} &\varOmega_{g} = \{x \in \varOmega: \rho g h > \rho_{w} g \langle -b\rangle\},\\ &\varOmega_{f} = \{x \in \varOmega: \rho g h < \rho_{w} g \langle -b\rangle\},\\ &\varSigma = \bar{\varOmega}_{g} \cap \bar{\varOmega}_{f}.\quad\quad\quad\quad\quad\quad\quad\quad \end{gathered}\right.\end{equation}

The symbol $\langle {\cdot } \rangle = \max ({\cdot }, 0)$ corresponds to the Macaulay brackets. Hence, the grounded portion includes both the parts where the bedrock lies above the sea level (i.e. where $\langle - b\rangle = 0$), as well as the parts where the bedrock lies below the sea level, but where there is too much ice for it to be floating (i.e. where $\rho g h > -\rho _{w} g b$).

In the simplest configuration, such as the one shown in figure 1, the grounded and floating portions can be written as open sets $\varOmega _{g} = (0, x_{gl})$ and $\varOmega _{g} = (x_{gl}, x_{c})$, so that the grounded-line position can be properly defined as the unique element of $\varSigma$: $\varSigma = \{ x_{gl}\}$. We note that, in general, the geometry might be more complex. For example, there could be several isolated points on which the ice sheet switches from a grounded to a floating position and vice versa, leading to multiple grounding lines. A more exotic configuration, not considered here, is the one described by Pegler (Reference Pegler2018a) with the so-called marginal-flotation zones. In that case, the interface $\varSigma$ becomes a set of its own, i.e. the grounding-line width becomes finite.

2.1.2. Boundary conditions

At $x=0$, we assume the ice to be sufficiently slow so that it is virtually motionless (this could also correspond to a symmetry condition):

(2.5)\begin{equation} u = 0, \quad \text{at } x = 0. \end{equation}

At the calving front, equilibrium between the horizontal stress in the ice and the ocean water pressure yields the following Neumann boundary condition:

(2.6)\begin{equation} 2A^{-{1}/{n}} \left\vert \dfrac{\mathrm{d} u}{\mathrm{d}\kern 0.06em x}\right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} u}{\mathrm{d}\kern 0.06em x} = \dfrac{1}{2}\rho \left(1 - \dfrac{\rho}{\rho_{w}}\right)gh, \quad \text{at } x = x_{c}. \end{equation}

Actually, if one considers an ice shelf without lateral drag and restricts the domain to the grounded part $\varOmega _{g}$ only, which we will do in this study, then this boundary condition can still be used, i.e.

(2.7)\begin{equation} 2A^{-{1}/{n}} \left\vert \dfrac{\mathrm{d} u}{\mathrm{d}\kern 0.06em x}\right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} u}{\mathrm{d}\kern 0.06em x} = \dfrac{1}{2}\rho \left(1 - \dfrac{\rho}{\rho_{w}}\right)gh, \quad \text{at } x = x_{gl}. \end{equation}

Indeed, (2.2b) with $\varLambda = 0$ implies that the quantity

(2.8)\begin{equation} \left[2A^{-{1}/{n}} h \left\vert \dfrac{\mathrm{d} u}{\mathrm{d}\kern 0.06em x}\right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} u}{\mathrm{d}\kern 0.06em x} - \dfrac{1}{2}\rho \left(1 - \dfrac{\rho}{\rho_{w}}\right)gh\right] \end{equation}

is conserved through the ice shelf.

2.2. Friction laws

2.2.1. Power-law friction laws

The simplest friction law is the Weertman friction law, for which $\tau _{b}$ is proportional to $\vert u \vert ^p$ with $p > 0$ (Weertman Reference Weertman1957). Usually, $p=1/3$ is chosen. To take into account effective pressure, one can use the so-called Budd friction law (Budd, Keage & Blundy Reference Budd, Keage and Blundy1979) for which

(2.9)\begin{equation} \tau_{b} = C N^{q}\,\vert u \vert^{p-1} u, \end{equation}

with $C$ a friction coefficient, $N$ an effective pressure and $p,q \ge 0$. Two elementary effective-pressure models are presented in § 2.2.5. The Budd friction law can be rewritten as a sliding law, i.e. the velocity can be written as a function of the basal friction stress:

(2.10)\begin{equation} u = C^{-{1}/{p}}N^{-{q}/{p}} \vert \tau_{b} \vert^{{1}/{p}-1} \tau_{b}. \end{equation}

It can also be noted that the law in (2.9) includes as a particular case the Weertman friction law if one sets $p = 1/3$ and $q=0$.

2.2.2. Coulomb friction law

A Coulomb behaviour assumes that there is a yield stress $\tau _{y} = C N$ that must be reached for ice to be sliding:

(2.11a,b) \begin{equation} \left\{ \begin{gathered} \tau_{b} = C N \text{sgn}(u), &\quad \text{if $\vert u \vert > 0$},\\\quad\quad \vert\tau_{b}\vert \le C N, &\quad \text{if $u = 0$}. \end{gathered}\right.\end{equation}

If the ice velocity is non-zero everywhere, then $\tau _{b} = C N \text {sgn}(u)$, which formally corresponds to a Budd friction law with $p=0$ and $q=1$. In the rest of this paper, we will always consider this case.

2.2.3. Hybrid friction laws

Tsai et al. (Reference Tsai, Stewart and Thompson2015) have considered a hybrid law that combines the Weertman and Coulomb friction laws:

(2.12)\begin{equation} \tau_{b} = \min(A_{s}^{{-}p}\vert u \vert^{p}, {C}N)\,\text{sgn}(u), \end{equation}

with $A_{s}$ a friction coefficient that controls the onset of the plastic behaviour. Such a law was originally introduced in Schoof (Reference Schoof2010). Smoothed versions have already been studied analytically and numerically (Schoof Reference Schoof2005, Reference Schoof2010; Gagliardini et al. Reference Gagliardini, Cohen, Råback and Zwinger2007). They take the following form:

(2.13)\begin{equation} \tau_{b} = C \left( \dfrac{\vert u \vert}{\vert u \vert + A_{s}C^{{1}/{p}}N^{{1}/{p}}}\right)^{p}N\,\text{sgn}(u), \end{equation}

or, by introducing $u_0 = A_{s}C^{{1}/{p}}N^{{1}/{p}}$,

(2.14)\begin{equation} \tau_{b} = C \left( \dfrac{\vert u \vert}{\vert u \vert + u_0}\right)^{p}N\,\text{sgn}(u). \end{equation}

This type of law, which exhibits viscoplastic behaviour, is interesting from a modelling perspective because it can be used to include both form and skin drag, even if these are distinct mechanisms (Minchew & Joughin Reference Minchew and Joughin2020). Form drag is associated with friction induced by ice deformation around obstacles and can be modelled with a power-law friction law, while skin drag is associated with friction induced by shear stress at the ice–bedrock interface, and can be modelled with a Coulomb friction law. Recently, Zoet & Iverson (Reference Zoet and Iverson2015, Reference Zoet and Iverson2020) have shown that such laws are in good agreement with experimental results.

2.2.4. Summary

The friction laws that we will consider in this article are shown in figure 2.

Figure 2. The considered friction models. The same notation $C$ is used for the friction coefficient in every friction law although those coefficients are not necessarily comparable to one another.

2.3. Dimensionless formulation

We introduce scales $[x], [h], [u]$, and $[\tau _{b}]$, leading to the dimensionless variables

(2.15)\begin{equation} \hat{x} = \dfrac{x}{[x]}, \quad \hat{h} = \dfrac{h}{[h]}, \quad \hat{b} = \dfrac{b}{[h]}, \quad \hat{u} = \dfrac{u}{[u]}, \quad \hat{\tau}_{b} = \dfrac{\tau_{b}}{[\tau_{b}]}, \end{equation}

and to the following dimensionless ratios:

(2.16a,b) \begin{equation} \left\{ \begin{gathered} \alpha = \dfrac{a}{([u]/[x])[h]}, \quad \beta = \left(\dfrac{\mathrm{d} b}{\mathrm{d}\kern 0.06em x}\right)\dfrac{[x]}{[h]}, \quad \gamma = \dfrac{[\tau_{b}][x]}{\rho g [h]^2},\\ \delta = \dfrac{\rho_{w} - \rho}{\rho_{w}}, \quad \varepsilon = \dfrac{A^{-{1}/{n}}[u]^{{1}/{n}}}{2\rho g [x]^{{1}/{n}}[h]}, \quad \lambda = \dfrac{\varLambda [u]^{m}[x]}{\rho g [h]}. \end{gathered}\right.\end{equation}

These scales and ratios should be characteristic of ice streams. The problem can be further simplified by choosing the scales so that additional constraints on the dimensionless ratios are enforced, e.g. by setting some of them to a unit value. However, we postpone these assumptions to a later stage, where the context will provide justification for them. We also introduce the dimensionless flotation thickness $\hat {h}_{b}$ as $\hat {h}_{b} = (1-\delta )^{-1} \hat {b}$. With these notations, the governing equations become

(2.17a,b) \begin{equation} \left\{ \begin{gathered} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\dfrac{\mathrm{d}}{\mathrm{d} \hat{x}}(\hat{u}_{g}\,\hat{h}_{g}) = \alpha,\\ 4\varepsilon \dfrac{\mathrm{d}}{\mathrm{d} \hat{x}}\left( \hat{h}_{g} \left\vert\dfrac{\mathrm{d} \hat{u}_{g}}{\mathrm{d} \hat{x}} \right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} \hat{u}_{g}}{\mathrm{d} \hat{x}}\right) - \gamma\,\hat{\tau}_{b} - \lambda\,\hat{h}_{g} \vert \hat{u}_{g} \vert^{m-1}\hat{u}_{g} = \hat{h}_{g} \left(\dfrac{\mathrm{d} \hat{h}_{g}}{\mathrm{d} \hat{x}} + \beta\right), \end{gathered}\right.\end{equation}

for $0 < \hat {x} < \hat {x}_{gl}$,

(2.18) \begin{equation} \hat{h}_{g} =\hat{h}_{f}, \quad \hat{u}_{g}= \hat{u}_{f}, \quad \hat{h}_{g} \left\vert \dfrac{\mathrm{d} \hat{u}_{g}}{\mathrm{d}\kern 0.06em x} \right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} \hat{u}_{g}}{\mathrm{d}\kern 0.06em x}= \hat{h}_{f} \left\vert \dfrac{\mathrm{d} \hat{u}_{f}}{\mathrm{d}\kern 0.06em x} \right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} \hat{u}_{f}}{\mathrm{d}\kern 0.06em x}, \end{equation}

at $\hat {x} = \hat {x}_{gl}$,

(2.19a,b) \begin{equation} \left\{ \begin{gathered} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\dfrac{\mathrm{d}}{\mathrm{d} \hat{x}}(\hat{u}_{f}\,\hat{h}_{f}) = \alpha,\\ 4 \varepsilon\, \dfrac{\mathrm{d}}{\mathrm{d} \hat{x}}\left( \hat{h}_{f} \left\vert\dfrac{\mathrm{d} \hat{u}_{f}}{\mathrm{d} \hat{x}} \right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} \hat{u}_{f}}{\mathrm{d} \hat{x}}\right) - \lambda\,\hat{h}_{f} \vert \hat{u}_{f} \vert^{m-1}\hat{u}_{f} = \delta\,\hat{h}_{f} \dfrac{\mathrm{d} \hat{h}_{f}}{\mathrm{d} \hat{x}}, \end{gathered}\right.\end{equation}

for $\hat {x}_{gl} < \hat {x} < \hat {x}_{c}$, with the following boundary conditions:

(2.20a–c) \begin{equation} \left\{ \begin{gathered} \quad\quad\quad\quad\quad\quad\hat{u} = 0, & \quad \text{at } \hat{x} = 0,\\ \left\vert\dfrac{\mathrm{d} \hat{u}}{\mathrm{d} \hat{x}} \right\vert^{({1}/{n})-1}\dfrac{\mathrm{d} \hat{u}}{\mathrm{d} \hat{x}} = \dfrac{\delta\hat{h}}{8 \varepsilon}, & \quad \text{at } \hat{x} = \hat{x}_{gl},\\ \quad\quad\quad\quad\quad\quad\hat{h} = \hat{h}_{b}, & \quad \text{at } \hat{x} = \hat{x}_{gl}. \end{gathered}\right.\end{equation}

2.4. Flux conditions

A flux condition is an expression of the grounding-line flux $q_{gl} \equiv h(x_{gl}) u(x_{gl})$ as a function of the different physical parameters $A, C,\ldots$ that appear in the problem formulation. Such an expression usually takes the form of an approximation that is valid within an asymptotic regime associated with the magnitude of the previously introduced dimensionless ratios. Historically, the first flux condition was derived by Schoof (Reference Schoof2007b) for the Weertman friction law. They considered an unbuttressed ice sheet, i.e. $\lambda = 0$, a scaling and a bed geometry such that $\alpha \sim 1, \gamma \sim 1$ and $\vert \beta \vert \lesssim 1$, and they assumed that $\varepsilon \ll 1$ and $\delta \ll 1$. Tsai et al. (Reference Tsai, Stewart and Thompson2015) derived a flux condition under the same assumptions, but for the Coulomb friction law. They showed that the resulting flux condition was more sensitive compared with the one derived by Schoof (Reference Schoof2007b). The importance of buttressing, i.e. the case $\lambda \neq 0$, was discussed by Pegler (Reference Pegler2016, Reference Pegler2018a,Reference Peglerb), Schoof et al. (Reference Schoof, Davis and Popa2017) and Haseloff & Sergienko (Reference Haseloff and Sergienko2018, Reference Haseloff and Sergienko2022). They showed that taking into account lateral drag could significantly change the dynamics of ice sheets, in particular by modifying the stability criterion that was previously derived for unbuttressed ice sheets (Schoof Reference Schoof2012). The regime of low basal stress, $\gamma \ll 1$, was covered by Sergienko & Wingham (Reference Sergienko and Wingham2019). The same authors also discussed the importance of $\alpha$ and $\beta$, showing that the so-called marine ice-sheet instability hypothesis was not applicable in general (Sergienko Reference Sergienko2022b; Sergienko & Wingham Reference Sergienko and Wingham2022). Schoof et al. (Reference Schoof, Davis and Popa2017) studied the impact that calving laws have on the flux conditions. All these authors, except Tsai et al. (Reference Tsai, Stewart and Thompson2015), have considered the Weertman friction law in their studies. Recently, Sergienko & Haseloff (Reference Sergienko and Haseloff2023) studied the notion of stability of marine ice sheets submitted to a climate forcing for a broad class of friction laws. However, they did not derive flux conditions for the configuration studied in the present paper, which we describe hereafter.

In this paper, we derive flux conditions for the Budd friction law with two elementary effective-pressure models, and show how they can be extended to hybrid friction laws. We will use the same assumptions that were done by Schoof (Reference Schoof2007b) and Tsai et al. (Reference Tsai, Stewart and Thompson2015), namely, we consider an unbuttressed ice sheet ($\lambda = 0$), scales that are such that $\alpha, \beta$, and $\gamma$ are at most of order $O(1)$, and consider the asymptotic regimes $\varepsilon \ll 1$ and $\delta \ll 1$. We will discuss in a later section the validity of these hypotheses, and we will show how the flux conditions can be modified to remain valid in the event that $\alpha, \beta$, and $\gamma$ are not small or moderate.

3.1. Derivation of the flux condition

3.1.1. Equivalent dynamical system for the boundary-layer problem

One can expand the unknown fields as powers of $\varepsilon$ and keep the leading-order terms because $\varepsilon$ is typically very small – approximately $10^{-3}$ for commonly used values of the physical parameters. One then expects an equilibrium between the friction and gravity terms in (3.3b), with the divergence of membrane stress which can be neglected. However, this balance fails in two cases. If $\delta$ is such that $\varepsilon \ll \delta$, then the Neumann boundary condition (3.3d) at the grounding line cannot be fulfilled. This hints at the existence of a boundary layer near the grounding line, in which the membrane-stress divergence becomes relatively important. Furthermore, if the friction stress reaches a zero value at the grounding line (e.g. if $1_{\text {A}} = 1$ and $q \neq 0$), then all the terms appearing in (3.3b) must become very small close to the grounding line, leading again to a boundary layer. In what follows, we place ourselves in one of these two cases so that we expect the presence of a boundary layer close to the grounding line.

To solve a very similar problem, Schoof (Reference Schoof2007b) and Tsai et al. (Reference Tsai, Stewart and Thompson2015) used the method of matched asymptotics: the solution inland, known as the outer solution, was matched with the so-called inner solution associated with the boundary layer. To obtain this inner solution, they introduced a scaling of the form

(3.4)\begin{equation} \hat{x}_{gl} - \hat{x} = \varepsilon^{\kappa_x} X, \quad \hat{h} = \varepsilon^{\kappa_{h}} H, \quad \hat{h}_{b} = \varepsilon^{\kappa_{h}} H_{b}, \quad \hat{b} = \varepsilon^{\kappa_{h}} B, \quad \hat{u} = \varepsilon^{\kappa_{u}} U, \end{equation}

where $\kappa _x, \kappa _h$ and $\kappa _u$ are chosen in a such way that the divergence of membrane stress, the friction stress and the gravity stress are of the same order of magnitude near the grounding line; in other words, they are of all of order $O(\varepsilon ^\kappa )$ for a same exponent $\kappa$. Furthermore, they are chosen such that the flux $Q = H U$ is $O(1)$ at the grounding line. This leads in the current context to the following exponent values:

(3.5)\begin{equation} \kappa_x = \dfrac{n(p-q+2)}{n + (p-q) + 3}, \quad \kappa_u ={-}\dfrac{n}{n + (p-q) + 3}, \quad \kappa_h = \dfrac{n}{n + (p-q) + 3}. \end{equation}

We remark that with the assumed values for $n, p$, and $q$, we have $\kappa _x > 0, \kappa _u < 0$ and ${\kappa _h > 0}$. At leading order, the flux $Q$ is then constant within the boundary layer, and we replace it by the grounding-line flux $Q_{gl}$.

The inner problem can be further transformed. As in Schoof (Reference Schoof2007b) and Tsai et al. (Reference Tsai, Stewart and Thompson2015), the solution to the inner problem is written as a trajectory of a two-dimensional dynamical system of the form $\tilde {X} \mapsto (\tilde {U}, \tilde {W})$, where $\tilde {X}, \tilde {U}$ and $\tilde {W}$ are respectively a scaled spatial coordinate, a scaled velocity and a scaled membrane stress, thus allowing the dynamics of the system to be interpreted in the phase plane $(\tilde {U}, \tilde {W})$. To obtain this dynamical system, the following change of variables is introduced:

(3.6) \begin{equation} \left.\begin{array}{c@{}} X ={H_{gl}}^{({2-q-np})/({p+1})} \tilde{X}, \quad U = {H_{gl}}^{({2-q+n})/({p+1})} \tilde{U},\\ - \vert U_X \vert^{({1}/{n})-1} U_X = H_{gl} \tilde{W},\quad {Q_{gl}} = {H_{gl}}^{({n+(p-q)+3})/({p+1})} {\tilde{Q}{_{gl}}}. \end{array}\right\}\end{equation}

At leading order, the following leading-order system is then obtained:

(3.7a–d) \begin{equation} \left\{ \begin{gathered} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\dfrac{\mathrm{d} \tilde{U}}{\mathrm{d} \tilde{X}} ={-} \vert \tilde{W} \vert^{n-1} \tilde{W}, &\quad \text{for } \tilde{X} >0,\\ \dfrac{\mathrm{d} \tilde{W}}{\mathrm{d} \tilde{X}} ={-}\dfrac{\vert \tilde{W}\vert^{n+1}}{\tilde{U}} - \dfrac{1}{4}\dfrac{\tilde{U}}{{\tilde{Q}{_{gl}}}} \left(\dfrac{{\tilde{Q}{_{gl}}}}{\tilde{U}} - 1_{\text{A}}\right)^{q}\vert \tilde{U}\vert^{p-1}\tilde{U}\phantom{=0,} & \nonumber\\\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad +\dfrac{{\tilde{Q}{_{gl}}}\vert \tilde{W}\vert^{n-1}\tilde{W}}{4 \tilde{U}^2}, &\quad \text{for } \tilde{X} >0,\\\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (\tilde{U},\tilde{W})=({\tilde{Q}{_{gl}}},{\delta}/{8}), & \quad \text{at } \tilde{X} = 0,\\\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (\tilde{U},\tilde{W}) \rightarrow (0,0), & \quad \text{as } \tilde{X} \rightarrow + \infty. \end{gathered}\right.\end{equation}

Equation (3.7d) is a matching condition and follows from the fact that the inner and outer solutions must be of the same order in an intermediate region. Because $U = \hat {u} \varepsilon ^{-\kappa _u}$ and the outer solution is such that $\hat {u} \sim 1$, we must enforce $U \rightarrow 0$, and therefore $\tilde {U} \rightarrow 0$, outside of the boundary layer. Similarly, $U_X = O(\varepsilon ^{\kappa _x-\kappa _u})$, and thus $\tilde {W} \rightarrow 0$ outside of it.

3.1.2. Flux condition

The rescaled flux at the grounding line, $\tilde {Q}_{gl}$, appears as a free parameter in (3.7). In the following section we will provide a justification for the existence of a trajectory that follows the flow defined by (3.7a) and (3.7b) and satisfies the boundary condition (3.7c) for a unique value of $\tilde {Q}_{gl}$ dependent on the effective-pressure model and the parameters $n, p, q$ and $\delta$. Then

(3.8)\begin{equation} \tilde{Q}_{gl} \equiv \tilde{Q}_{gl}(1_{\text{A}}, n, p, q, \delta). \end{equation}

This numerical value can be computed using the numerical method described in the Appendix B. Using (3.6), it is possible to switch back to the original variables. The flux at the grounding line is then given by the following expressions, for the N$_\text {A}$ and N$_\text {B}$ effective-pressure models, respectively:

(3.9a)$$\begin{gather} q_{gl} = {\tilde{Q}{_{gl}}} (\rho g)^{-({q-1})/({p+1})}(2\rho g)^{{n}/({p+1})}C^{-{1}/({p+1})}A^{{1}/({p+1})} h_{gl}^{({n+(p-q)+3})/({p+1})}, \end{gather}$$

(3.9b)$$\begin{gather}q_{gl} = {\tilde{Q}{_{gl}}} (\rho g)^{-({q-1})/({p+1})}(2\rho g)^{{n}/({p+1})}[C(1-c)^{q}]^{-{1}/({p+1})}A^{{1}/({p+1})} h_{gl}^{({n+(p-q)+3})/({p+1})}. \end{gather}$$

3.1.3. Impact of the relative ice–water density difference

Tsai et al. (Reference Tsai, Stewart and Thompson2015) also showed a way to derive the approximate dependence of $\tilde {Q}_{gl}$ on $\delta$. The idea is to remark that if $\delta$ is treated as a small parameter in (3.7), then $\tilde {U}_{\tilde {X}} \approx 0$ within the boundary layer. This observation supports the introduction of a new scaling so that this term becomes $O(1)$ at the grounding line. With

(3.10)\begin{align} \tilde{X} = (\delta/8)^{r_1} \check{X}, \quad \tilde{Q}_{gl} = (\delta/8)^{r_2} \check{Q}_{gl}, \quad \tilde{U} = (\delta/8)^{r_2} [\check{Q} - (\delta/8) \check{U}], \quad \tilde{W} = (\delta/8) \check{W}, \end{align}

a distinguished limit can be obtained, in which the dominant powers of $\delta$ balance each other. For the N$_\text {A}$ model, a distinguished limit is achieved for

(3.11)\begin{equation} r_1 = [{(p-q+1)-np}]/({p+1}) \quad \text{and} \quad r_2 = ({n - q})/({p+1}). \end{equation}

For the N$_\text {B}$ model, a distinguished limit is obtained for

(3.12)\begin{equation} r_1 = [{(p+1)-np}]/({p+1}) \quad \text{and} \quad r_2 = {n}/({p+1}). \end{equation}

Finally, the following flux at the grounding line is obtained for the N$_\text {A}$ and N$_\text {B}$ effective-pressure models:

(3.13a) \begin{align} q_{gl} &= {\check{Q}{_{gl}}} \left(\dfrac{1-\rho/\rho_{w}}{8}\right)^{({n-q})/({p+1})} (\rho g)^{-({q-1})/({p+1})}(2\rho g)^{{n}/({p+1})}\nonumber\\ &\quad \times C^{-{1}/({p+1})}A^{{1}/({p+1})} h_{gl}^{({n+(p-q)+3})/({p+1})}, \end{align}

(3.13b) \begin{align} q_{gl} &= {\check{Q}{_{gl}}} \left(\dfrac{1-\rho/\rho_{w}}{8}\right)^{{n}/({p+1})} (\rho g)^{-({q-1})/({p+1})}(2\rho g)^{{n}/({p+1})}\nonumber\\ &\quad \times[C(1-c)^{q}]^{-{1}/({p+1})}A^{{1}/({p+1})} h_{gl}^{({n+(p-q)+3})/({p+1})}. \end{align}

This scaling ${\tilde {Q}{_{gl}}} = (\delta /8)^{r_2} {\check {Q}{_{gl}}}$, that is, the way in which ${\tilde {Q}{_{gl}}}$ depends on $\delta$, is verified numerically in figure 3.

Figure 3. Comparison between values of ${\tilde {Q}{_{gl}}}$ obtained numerically (circles) and the scaling ${\tilde {Q}{_{gl}}} \propto (\delta /8)^{r_2}$ (lines) for several friction laws and effective-pressure models. The lines obey the equation ${\tilde {Q}{_{gl}}} = {\tilde {Q}{_{gl}}}\vert _{\delta =0.1} (\delta /0.1)^{r_2}$ with $r_2=(n-1_{\text {A}}q)/(p+1)$. In (b), the Weertman and the Budd results coincide, as expected.

3.2. Analysis of the leading-order dynamical system

We now consider the analysis of the dynamical system governed by the system of (3.7). More precisely, we motivate the existence of a solution for a unique value of the grounding-line flux $\tilde {Q}_{gl}$ by considering separately the case where the friction stress vanishes, or not, at the grounding line.

3.2.1. Strategy

To study the leading-order dynamical system, we first rewrite the system of equations in a way that allows the dynamics close to the origin in the $(\tilde {U}, \tilde {W})$ phase plane, i.e. for $\tilde {X}\rightarrow + \infty$, to be studied. To this end, we rewrite this system in terms of new variables $\mathcal {X}, \xi, \varPsi$ and $\mathcal {Q}_{gl}$. The interpretation of these variables is the following: $\mathcal {X}$ plays the role of a spatial coordinate, $\xi$ is a rescaled velocity, $\varPsi$ is a measure of the ratio of friction stress over gravity stress and $\mathcal {Q}_{gl}$ is a rescaled grounding-line flux. The specific form that these variables take will be described separately for the case in which friction vanishes at the grounding line, and the case in which it does not. A problem of the following form is then obtained:

(3.14a–d) \begin{equation} \left\{ \begin{gathered} \dfrac{\mathrm{d} \xi}{\mathrm{d} \mathcal{X}} = {F}_\xi(\xi, \varPsi, \mathcal{X}; \mathcal{Q}_{gl}),& \quad \text{for } \mathcal{X} > 0,\\ \dfrac{\mathrm{d} \varPsi}{\mathrm{d} \mathcal{X}} = {F}_\varPsi(\xi, \varPsi, \mathcal{X}; \mathcal{Q}_{gl}),& \quad \text{for } \mathcal{X} > 0,\\ (\xi,\varPsi) = ({G}_\xi(\mathcal{Q}_{gl}),{G}_\varPsi(\mathcal{Q}_{gl})), & \quad \text{at } \mathcal{X} = 0,\\\quad\quad\quad\quad\quad\quad\quad\quad (\xi,\varPsi) \rightarrow (0,1), & \quad \text{as } \mathcal{X} \rightarrow +\infty. \end{gathered}\right.\end{equation}

We then identify the point $(\xi,\varPsi )=(0,1)$ as a fixed point, and study the dynamics of the flow defined by (3.14a) and (3.14b) close to that point. It turns out that the only way to reach the fixed point is through a centre manifold that is unique. Therefore, if a solution to the problem defined by (3.14) exists, it necessarily goes through this centre manifold. The question then amounts to finding whether an orbit that reaches this centre manifold, i.e. that obeys (3.14a), (3.14b) and (3.14d), can satisfy the boundary condition (3.14c). This last condition is in fact satisfied for exactly one value of the grounding-line flux $\mathcal {Q}_{gl}$. To show this, we introduce a mapping $D$ as follows:

(3.15)\begin{equation} D: (0, +\infty) \rightarrow \mathbb{R}: \mathcal{Q}_{gl} \mapsto D(\mathcal{Q}_{gl}) \equiv f(\mathcal{Q}_{gl})[{\varPsi}^{c}(G_\xi(\mathcal{Q}_{gl}); \mathcal{Q}_{gl}) - {G}_\varPsi(\mathcal{Q}_{gl})], \end{equation}

in which $f$ is a strictly positive or a strictly negative function and ${\varPsi }^{c}(\xi, \mathcal {Q}_{gl})$ is the $\varPsi$ coordinate of the centre manifold at position $\xi$. To satisfy (3.14c), it is then necessary and sufficient that $D(\mathcal {Q}_{gl})=0$ for some $\mathcal {Q}_{gl}$. If, in addition, $D$ is a strictly monotonic function, then this root is unique. Overall, this means that there is exactly one value of $\mathcal {Q}_{gl}$ that leads to a solution of (3.14), and the solution to the leading-order dynamical system exists and is unique.

To simplify the notations in what follows, we define $c_1$ and $c_2$ by

(3.16)\begin{equation} c_1 = 1+({p-q+3})/{n} \quad \text{and} \quad c_2 = 1-({p-q+3})/{n}. \end{equation}

We note that, for the assumed ranges of values of $n, p$ and $q$, the following inequalities hold:

(3.17)\begin{equation} c_1 > 1 \quad \text{and} \quad -1 < c_2 < 1. \end{equation}

3.2.2. Non-vanishing friction at the grounding line

We first consider the case of a non-vanishing friction stress at the grounding line, that is, a friction model with either an exponent $q = 0$, so that there is no dependence with respect to the effective pressure, or with the N$_\text {B}$ effective-pressure model. We note that this case shares similarities with the study considered in Schoof et al. (Reference Schoof, Davis and Popa2017), where the authors have included a lateral-drag term in their momentum balance. This term is of the form $\varLambda h \vert u \vert ^{m-1} u$, which is analogous to a Budd friction law with $p = m, q=1$ and the N$_\text {B}$ effective-pressure model. In fact, it can be noted that the Budd friction law taken with the N$_\text {B}$ effective-pressure model is effectively equivalent to considering a friction term dominated by lateral drag.

We introduce $\xi, \varPsi$ and $\mathcal {Q}_{gl}$ as

(3.18)\begin{equation} \xi = \tilde{Q}_{gl}^{({q-2})/{n}-1}\tilde{U}^{c_1}, \quad \varPsi = \tilde{Q}_{gl}^{-({q-2})/{n}} {\tilde{W}}\,{\tilde{U}^{1-c_1}}, \quad \mathcal{Q}_{gl} = \tilde{Q}_{gl}^{({p+1})/{n}}, \end{equation}

and $\mathcal {X}$ as

(3.19a,b) \begin{equation} \left\{ \begin{gathered} \mathcal{X} = \int_0^{\tilde{X}} s(\xi({X}), \varPsi({X}))\,\text{d}{X},\quad\quad\quad\quad\quad\quad\quad\quad\quad\\ s(\xi, \varPsi)=\tilde{Q}_{gl}^{({q-2+np})/{n c_1}} \xi^{({n(p-q)-(p-q+3)+n})/({(p-q+3)+n})}\vert\varPsi\vert^{n-1}\varPsi. \end{gathered}\right.\end{equation}

The system (3.7) then becomes

(3.20a–d) \begin{equation} \left\{ \begin{gathered} \dfrac{\mathrm{d} \xi}{\mathrm{d} \mathcal{X}} ={-}c_1 \xi^2,& \quad \text{for } \mathcal{X} > 0,\\ \dfrac{\mathrm{d} \varPsi}{\mathrm{d} \mathcal{X}} ={-} c_2 \xi\varPsi - \dfrac{1}{4}\vert \varPsi \vert^{{-}n-1}\varPsi + \dfrac{1}{4},& \quad \text{for } \mathcal{X} > 0,\\ (\xi,\varPsi) = (\mathcal{Q}_{gl},\mathcal{Q}_{gl}^{{-}1}\delta/8), & \quad \text{at } \mathcal{X} = 0,\\ (\xi,\varPsi) \rightarrow (0,1), & \quad \text{as } \mathcal{X} \rightarrow +\infty. \end{gathered}\right.\end{equation}

It can be remarked that $\mathcal {Q}_{gl}$ completely disappears from the differential equations and is only present in the boundary conditions. This system is similar to the system considered by Schoof (Reference Schoof2011), where they considered the Weertman friction law. The only differences are the values of the parameters $c_1$ and $c_2$ which, in our case, could depend on $q$ if we consider the N$_\text {B}$ effective-pressure model. The method used in Schoof (Reference Schoof2011) to show the existence and uniqueness of a solution can still be applied. We briefly describe it, the calculations being analogous.

The idea of Schoof (Reference Schoof2011) to show existence and uniqueness properties of a similar system is to consider the characterisation of $(\xi, \varPsi )=(0, 1)$ as a fixed point that can only be reached through a centre manifold that is unique, as well as the evolution of the product $\varPsi \xi$ along that manifold. They showed that this product was equal to zero at the fixed point, and increasing without bound for increasing values of $\xi$ along that orbit. It then follows that there is exactly one value of $\mathcal {Q}_{gl}$ that satisfies (3.20c), which shows the existence and uniqueness of a solution. These ideas can still be applied to the more general case that is considered here.

The reasoning can also be made with respect to the mapping $D$ defined in (3.15) by choosing $f(\mathcal {Q}_{gl}) = \mathcal {Q}_{gl}$. Indeed, the centre manifold is independent of $\mathcal {Q}_{gl}$, so $\varPsi ^{c}(\xi ; \mathcal {Q}_{gl}) \equiv \varPsi ^{c}(\xi )$. Furthermore, the mapping $\xi \mapsto \xi \varPsi ^{c}(\xi )$ increases without bound with $\xi$. Therefore, the mapping

(3.21)\begin{equation} \mathcal{Q}_{gl} \quad \mapsto \quad D(\mathcal{Q}_{gl}) = \mathcal{Q}_{gl}\varPsi^{c}(\mathcal{Q}_{gl}) - ({\delta}/{8})\end{equation}

also increases without bound with $\mathcal {Q}_{gl}$. Because $\xi \varPsi ^{c}(\xi )=0$ for $\xi = 0$, we also have $D(0)=-\delta /8 < 0$. Hence, $D$ has exactly one root, which concludes the discussion.

3.2.3. Vanishing friction at the grounding line

We now consider friction laws that vanish at the grounding line, namely friction laws that involve the N$_\text {A}$ effective-pressure model (in particular, we consider that $q \neq 0$). In that case, it cannot be shown that the product $\varPsi \xi$ increases monotonically with $\xi$ along an orbit that reaches the centre manifold. Geometrically, the hyperbola $\varPsi = (\delta /8)/\xi$ will not necessarily intersect the solution trajectory at a single location.

We propose another strategy. Specifically, we consider another change of variables for $\xi$, namely, $\xi = ({\tilde {U}}/{\tilde {Q}_{gl}})^{{1}/{2}}$, and we take $f(\mathcal {Q}_{gl}) = 1$ in (3.15). This change of variables is similar to the one described in the supplementary material of Schoof et al. (Reference Schoof, Davis and Popa2017). We will also limit ourselves to the Budd friction law with a linear dependence with respect to the effective pressure, that is, $q=1$. For that value, we note that $c_2 > 0$. The system (3.7) then becomes

(3.22a–d) \begin{equation} \left\{ \begin{gathered} \quad\quad\quad\quad\quad\quad\quad\dfrac{\mathrm{d} \xi}{\mathrm{d} \mathcal{X}} ={-}\dfrac{1}{2} \mathcal{Q}_{gl} \xi^{2c_1+1},& \quad \text{for } \mathcal{X} > 0,\\ \dfrac{\mathrm{d} \varPsi}{\mathrm{d} \mathcal{X}} ={-} c_2\,\mathcal{Q}_{gl}\xi^{2c_1}\varPsi - \dfrac{1}{4}\vert \varPsi \vert^{{-}n-1}\varPsi (1 - \xi^2) + \dfrac{1}{4},& \quad \text{for } \mathcal{X} > 0,\\\quad\quad\quad\quad\quad\quad\quad (\xi,\varPsi) = (1,\mathcal{Q}_{gl}^{{-}1}\delta/8),& \quad \text{at } \mathcal{X} = 0,\\\quad\quad\quad\quad\quad\quad\quad (\xi,\varPsi) \rightarrow (0,1),& \quad \text{as } \mathcal{X} \rightarrow +\infty, \end{gathered}\right.\end{equation}

and the mapping $D$ becomes

(3.23)\begin{equation} \mathcal{Q}_{gl} \quad \mapsto \quad D(\mathcal{Q}_{gl}) = {\varPsi}^{c}(1; \mathcal{Q}_{gl}) - ({\delta}/{8})\mathcal{Q}_{gl}^{{-}1}. \end{equation}

As before, the only fixed point in the system is the point $(\xi, \varPsi )=(0, 1)$, which corresponds to the boundary condition (3.22d). We can again determine that the only way to reach this point is through a centre manifold. In contrast to the previous case, $\mathcal {Q}_{gl}$ appears in the definition of the flow defined by (3.22a) and (3.22b), so the centre manifold depends on $\mathcal {Q}_{gl}$. To demonstrate that $D$ possesses exactly one root, we then show, based on asymptotic expansions, that the following properties hold: (i) $D$ is a continuous mapping, (ii) $\text {d}D/\text {d}\mathcal {Q}_{gl} > 0$ for all $\mathcal {Q}_{gl} > 0$, (iii) $\lim _{\mathcal {Q}_{gl} \rightarrow + \infty } D(\mathcal {Q}_{gl}) > 0$ and (iv) $D(\delta /8) < 0$. The details of this analysis can be found in the Appendix A.

3.3. Existence of a boundary layer

It can be remarked that, for some configurations, we obtain $\check {Q}_{gl} \approx 1$. In fact, these configurations are those that are such that friction at the grounding line does not vanish, i.e. they correspond to a friction law with $q = 0$, or with $q > 0$ but with the effective-pressure model N$_\text {B}$. In that case, no boundary layer is needed close to the grounding line, and the membrane-stress divergence can be neglected. Indeed, the flux condition can be obtained by simply combining a balance between the friction and the gravity stresses and the boundary conditions at the grounding line. For the Budd friction law, that approach yields

(3.24)\begin{equation} C N^{q} u^{p} \approx{-} \rho g h \dfrac{\mathrm{d}}{\mathrm{d}\kern 0.06em x}(b + h). \end{equation}

With the assumption that the bedrock slope $\text {d}b/\text {d}x$ is negligible ($\vert \beta \vert \lesssim 1$) and that the flux divergence $\text {d} q_{adv}/\text {d}x$ is not too large ($\alpha = 1$), and in particular much smaller than $q_{adv} (\text {d}u/\text {d}x)/u$, we have

(3.25)\begin{equation} \dfrac{\mathrm{d}}{\mathrm{d}\kern 0.06em x}(b + h) \approx \dfrac{\mathrm{d} h}{\mathrm{d}\kern 0.06em x} \approx q_{adv}\dfrac{\mathrm{d}}{\mathrm{d}\kern 0.06em x} \left(\dfrac{1}{u}\right). \end{equation}

Using this relation in (3.24) and combining it with the grounding-line boundary condition (2.7) leads to the following relation at the grounding line:

(3.26)\begin{equation} C N^{q} \left(\dfrac{q_{gl}}{h_{gl}}\right)^{p} \approx \rho g \dfrac{h_{gl}^3}{q_{gl}} \left(\dfrac{1}{4} \rho \left(1 - \dfrac{\rho}{\rho_{w}}\right)g \right)^{n} h_{gl}^{n} A, \end{equation}

that is,

(3.27)\begin{equation} q_{gl} \approx \left(\dfrac{\rho g}{C}\right)^{{1}/({p+1})} N^{-{q}/({p+1})}\left(\dfrac{1}{4} \rho\left(1 - \dfrac{\rho}{\rho_{w}}\right)g \right)^{{n}/({p+1})}A^{{1}/({p+1})}h_{gl}^{({n+p+3})/({p+1})}. \end{equation}

This relation corresponds to our flux condition (3.13b) with $\check {Q}_{gl}\approx 1$, as announced. In fact, the observation that the membrane-stress divergence can be neglected to derive the flux condition has been remarked by Schoof (Reference Schoof2007b, Reference Schoof2011) in their derivation for the Weertman friction law, and later by Sergienko & Wingham (Reference Sergienko and Wingham2022) who revisited their boundary-layer analysis. In particular, Sergienko & Wingham (Reference Sergienko and Wingham2022) have shown that, for $\varepsilon \ll \delta$, which is what is assumed here, the boundary layer is very weak, and this observation can be explained by the nonlinearity of the governing equations. Furthermore, the boundary layer will become increasingly weak as $\delta$ becomes smaller.

However, this analysis is not valid for a combination of friction law and effective-pressure model that is such that friction stress vanishes at the grounding line. For these configurations, there is another stress regime in the vicinity of the grounding line. In our analysis, this takes the form of a boundary layer, which is necessary to obtain the correct flux condition. If it were not the case, then one would obtain $\check {Q}_{gl}=1$. It follows that the fact that $\check {Q}_{gl}$ takes a value distinct from unity reflects the importance of membrane-stress divergence in the boundary layer. This key observation was already made by Tsai et al. (Reference Tsai, Stewart and Thompson2015) for the Coulomb law, and is here confirmed for the more general Budd friction law.

The distinction between these two distinct behaviours can be observed in solutions to the different formulations of the problem that arise in the derivation of the flux conditions. First, let us consider the solutions to the problem written in its dimensionless form, namely to the system (3.3). We consider the Weertman law with $p=1/3$, the Coulomb law and the Budd law with $p=1/3$ and $q=1$, with both the N$_\text {A}$ and N$_\text {B}$ hydrological models. We take $\beta = -10^{-1}, \varepsilon = 6\times 10^{-4}$ and $\delta = 10^{-1}$. The solutions of the problem are shown in figure 4. The most striking difference concerns the ratio of the membrane-stress divergence and the gravity stress: this ratio is almost equal to zero in the entire grounded domain for the Weertman friction law, as well as for the Coulomb and Budd friction laws with the N$_\text {B}$ model. On the other hand, it becomes significant close to the grounding line for the Coulomb and Budd friction laws when they are coupled with the N$_\text {A}$ model, i.e. when the friction stress vanishes at the grounding line.

Figure 4. Solutions to the dimensionless problem for various friction laws, with the N$_\text {A}$ (continuous lines) and the N$_\text {B}$ effective-pressure model (dashed lines). Panel (c) shows the ratio of the membrane-stress divergence and the gravity stress.

A similar observation can be made if the problem is formulated in terms of $(\tilde {U}, \tilde {W})$, i.e. by considering the system (3.7). The solutions are shown in figure 5. Qualitatively, the solutions associated with vanishing grounding-line friction exhibit a stronger curvature in their trajectories. Importantly, the far-field solutions, shown in dotted lines and corresponding to a simple friction–gravity balance, do not represent well the dynamics close to the grounding line located at $\tilde {U}=\tilde {Q}_{gl}$.

Figure 5. Solutions to the problem formulated in terms of $(\tilde {U}, \tilde {W})$, for various friction laws, with the N$_\text {A}$ (continuous lines) and the N$_\text {B}$ effective-pressure model (dashed lines). The dotted lines correspond to the far-field solutions associated with the coupling of the Coulomb and Budd friction laws with the N$_\text {A}$ model.

Finally, this observation is also present in the version of the problem used in the analysis presented in the previous subsection, namely (3.20). Indeed, the solution is then obtained as a portion of an orbit that reaches the fixed point located at $(\xi, \varPsi ) = (0,1)$ through its centre manifold. The solution trajectory can be parametrised by $\xi \in [0, \mathcal {Q}_{gl}]$, where $\mathcal {Q}_{gl}$ is the $\xi$ coordinate of the intersection of the centre manifold with the hyperbola whose equation is $\xi \varPsi =\delta /8$. An asymptotic analysis reveals that the centre manifold is such that $\varPsi \sim 1$ for small $\xi$. It thus follows that, for small values of $\delta /8$, the solution trajectory is included in the region which is such that $\varPsi \sim 1$. Because $\varPsi$ is a scaled version of the ratio between friction and gravity stresses, the divergence of membrane stress can be neglected over the whole domain, even close to the grounding line. This argument is similar to the one developed in Schoof (Reference Schoof2011) for the Weertman friction law.

5.1. Non-vanishing friction law with $\gamma \sim 1$: negligible membrane-stress divergence

Let us consider a general Budd friction law $\tau _{b}={C}N ^{q}\vert u \vert ^{p-1}u$ for which the divergence of membrane stress can be neglected in the momentum-balance equation, so that no boundary layer is needed close to the grounding line to obtain the flux condition (i.e. for which $\check {Q}_{gl} \approx 1$ in the flux conditions that have been derived). To fulfil this condition, we consider a case where the effective pressure at the grounding line, $N_{gl}$, is non-zero, so that the friction stress does not vanish, and where we have $\gamma \sim 1$, so that the friction stress indeed balances the gravity stress. The N$_\text {B}$ effective-pressure model falls into the category of effective-pressure models that are such that $N_{gl} \neq 0$. In that case, the combination of the mass-balance equation, the momentum-balance equation and the grounding-line boundary condition yields the following algebraic equation at the grounding line:

(5.1)\begin{equation} C N_{gl}^{q} q_{gl}^{p+1} + \rho g q_{gl}h_{gl}^{p+1} \left(\dfrac{\text{d}b}{\text{d}x}\right)_{gl} = \rho g \left(\dfrac{1}{4}\rho \left(1 - \dfrac{\rho}{\rho_{w}}\right)g\right)^{n}A h_{gl}^{n+p+3} - \rho g h_{gl}^{p+2} a. \end{equation}

A similar expression can be found in Schoof (Reference Schoof2007a,Reference Schoofb) and in Sergienko & Wingham (Reference Sergienko and Wingham2022). If $a$ and $(\text {d}b/\text {d}x)_{gl}$ cannot be neglected, then no expression relating the grounding-line flux $q_{gl}$ to the grounding-line thickness $h_{gl}$ that is both exact and explicit can be obtained. We note that (5.1) can be written as

(5.2)\begin{equation} \dfrac{q_{gl}}{q_{{gl,c}}} \left(\left(\dfrac{q_{gl}}{q_{{gl,c}}}\right)^{p} + \dfrac{\beta}{\beta_{c}}\right) = 1 - \dfrac{\alpha}{\alpha_{c}}, \end{equation}

in which $q_{{gl,c}}$ is a reference flux, given by

(5.3)\begin{equation} q_{{gl,c}} = \left(\dfrac{\rho g}{C}\right)^{{1}/({p+1})} N_{gl}^{-{q}/({p+1})}\left(\dfrac{1}{4} \rho\left(1 - \dfrac{\rho}{\rho_{w}}\right)g \right)^{{n}/({p+1})}A^{{1}/({p+1})}h_{gl}^{({n+p+3})/({p+1})}. \end{equation}

In the case where the friction stress is non-zero at the grounding line, this reference flux is equal to the flux that would be obtained if $a$ and $(\text {d}b/\text {d}x)_{gl}$ could be neglected in (5.1), i.e. this is the expression of the flux that we have derived in § 3. The ratios ${\alpha }/{\alpha _{c}}$ and ${\beta }/{\beta _{c}}$ are defined as

(5.4)\begin{equation} \dfrac{\alpha}{\alpha_{c}} = \dfrac{a}{\left(\dfrac{1}{4}\rho \left(1 - \dfrac{\rho}{\rho_{w}}\right)g\right)^{n}A h_{gl}^{n+1}} \quad \text{and} \quad \dfrac{\beta}{\beta_{c}} = \dfrac{({\mathrm{d} b}/{\mathrm{d}\kern 0.06em x})_{gl}\,q_{{gl,c}}h_{gl}^{{-}1}}{\left(\dfrac{1}{4}\rho \left(1 - \dfrac{\rho}{\rho_{w}}\right)g\right)^{n}A h_{gl}^{n+1}}. \end{equation}

These ratios provide a way to quantify the impact of the hypotheses made to derive the flux conditions on these flux expressions, more precisely the discrepancy with respect to the reference flux value $q_{{gl,c}}$. This difference will be small if $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ are both small. We note that the denominators in (5.4) are proportional to the strain rate at the grounding line, so $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ can respectively be interpreted as a normalised measure of the variation of ice velocity associated with the net mass accumulation rate and the bedrock slope. These ratios can also be written with respect to the dimensionless numbers introduced in § 2: we have

(5.5)\begin{align} \dfrac{\alpha}{\alpha_{c}} = \alpha\left(\dfrac{\varepsilon}{\delta/8}\right)^{n} \left(\dfrac{h_{gl}}{[h]}\right)^{{-}n} \quad \text{and} \quad \dfrac{\beta}{\beta_{c}} = {\beta} {\gamma^{-{1}/({p+1})}} \left(\dfrac{\varepsilon}{\delta/8}\right)^{{np}/({p+1})} \left(\dfrac{h_{gl}}{[h]}\right)^{-{(p+q-np+1)}/({p+1})}, \end{align}

so in the limit of $\varepsilon \rightarrow 0$, the ratios $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ must tend towards zero. However, that limit is never reached in practice. Equation (5.4) then allows us to compute, quantitatively, the importance of $a, ({\mathrm {d} b}/{\mathrm {d} x})_{gl}$ and $C$ on the derivation of flux conditions, as a function of the original dimensional variables. Analogously, (5.4) allows us to assess the importance of $\alpha, \beta$ and $\gamma$ on the validity of our flux conditions.

The dimensionless ratios $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ can also be used to derive new flux conditions which are approximately valid in the case where those ratios are not small. Indeed, we can formally write that

(5.6)\begin{align} q_{gl} = \check{Q}_{gl} \left(\dfrac{\rho g}{C}\right)^{{1}/({p+1})} N_{gl}^{-{q}/({p+1})}\left(\dfrac{1}{4} \rho\left(1 - \dfrac{\rho}{\rho_{w}}\right)g \right)^{{n}/({p+1})}A^{{1}/({p+1})}h_{gl}^{({n+p+3})/({p+1})}, \end{align}

where $\check {Q}_{gl}=\check {Q}_{gl}({\alpha }/{\alpha _{c}},{\beta }/{\beta _{c}})$ is a correction factor. Note that, by construction, $\check {Q}_{gl}=q_{gl} / q_{{gl,c}}$. The expression of $\check {Q}_{gl}$ can be approximated by considering the algebraic equation (5.2). On the one hand, this equation can be solved numerically for several fixed values of ${\alpha }/{\alpha _{c}}$ and ${\beta }/{\beta _{c}}$ (figure 8a). The number of acceptable solutions of (5.2), i.e. of real and strictly positive solutions values for $\check {Q}_{gl}$, depends on the value of $\alpha /\alpha _{c}$. In fact, $\alpha /\alpha _{c} = 1$ plays the role of a bifurcation point. Indeed, for $\alpha /\alpha _{c} < 1$, there is exactly one acceptable solution $\check {Q}_{gl}$. It is also found that in that case, $\check {Q}_{gl}$ decreases with both ${\alpha }/{\alpha _{c}}$ and ${\beta }/{\beta _{c}}$. For $\alpha /\alpha _{c}=1$, there is exactly one acceptable solution, provided that $\beta /\beta _{c} < 0$; otherwise, there is no solution. For $\alpha /\alpha _{c} > 1$, we observe a folding of the solution branch $\beta /\beta _{c} \mapsto \tilde {Q}_{gl}({\alpha }/{\alpha _{c}},{\beta }/{\beta _{c}})$, which becomes multi-valued for $\beta /\beta _{c} < (\beta /\beta _{c})_{\ast }$, and for which there is no solution for $\beta /\beta _{c} > (\beta /\beta _{c})_{\ast }$. The critical value $(\beta /\beta _{c})_{\ast }$ is given by

(5.7)\begin{equation} \left(\dfrac{\beta}{\beta_{c}}\right)_{{\ast}} ={-}(p+1){p^{-{p}/({p+1})}} \left(\dfrac{\alpha}{\alpha_{c}}-1\right)^{{p}/({p+1})}. \end{equation}

On the other hand, (5.2) can be solved approximately based on asymptotic analysis. Specifically, the following asymptotic expressions hold. For $\alpha /\alpha _{c} < 1$,

(5.8)$$\begin{gather} \check{Q}_{gl} \sim \left(1 - \dfrac{\alpha}{\alpha_{c}}\right)^{{1}/({p+1})} - \dfrac{1}{p+1}\left(1 - \dfrac{\alpha}{\alpha_{c}}\right)^{({p-1})/({p+1})} \dfrac{\beta}{\beta_{c}}, \quad \text{for } \left\vert\dfrac{\beta}{\beta_{c}}\right\vert \ll 1, \end{gather}$$

(5.9)$$\begin{gather}\check{Q}_{gl} \sim \left(-\dfrac{\beta}{\beta_{c}}\right)^{{1}/{p}}, \quad \text{for }\dfrac{\beta}{\beta_{c}} <0 \text{ and }\left\vert\dfrac{\beta}{\beta_{c}}\right\vert \gg 1, \end{gather}$$

(5.10)$$\begin{gather}\check{Q}_{gl} \sim \left(1-\dfrac{\alpha}{\alpha_{c}}\right) \left(\dfrac{\beta}{\beta_{c}}\right)^{{-}1}, \quad \text{for }\dfrac{\beta}{\beta_{c}} >0 \text{ and } \left\vert\dfrac{\beta}{\beta_{c}}\right\vert \gg 1. \end{gather}$$

For $\alpha /\alpha _{c} = 1, \check {Q}_{gl} = (-{\beta }/{\beta _{c}})^{{1}/{p}}$, provided that ${\beta }/{\beta _{c}} <0$. For $\alpha /\alpha _{c} >1$, the upper and lower solution branches obey the following relations:

(5.11)$$\begin{gather} \check{Q}_{gl} \sim \left(-\dfrac{\beta}{\beta_{c}}\right)^{{1}/{p}}, \quad \text{for }\dfrac{\beta}{\beta_{c}} <0 \text{ and } \left\vert\dfrac{\beta}{\beta_{c}}\right\vert \gg 1, \end{gather}$$

(5.12)$$\begin{gather}\check{Q}_{gl} \sim \left(1-\dfrac{\alpha}{\alpha_{c}}\right) \left(\dfrac{\beta}{\beta_{c}}\right)^{{-}1}, \quad \text{for } \dfrac{\beta}{\beta_{c}} <0 \text{ and } \left\vert\dfrac{\beta}{\beta_{c}}\right\vert \gg 1. \end{gather}$$

It can be expected that the lower solution branch is seldom reached in practice as it corresponds to relatively large values of $a$ (as $\alpha /\alpha _{c} > 1$) but to relatively small values of the flux $q_{gl}$ (as $\check {Q}_{gl} \lesssim 1$). Combining these expressions together, a closed-form formula can be obtained to approximate the value of $\check {Q}_{gl}$. Assuming that we are in the case where there is a least one solution for $\check {Q}_{gl}$, i.e. considering the case $\alpha /\alpha _{c} < 1$ or $\beta /\beta _{c} < (\beta /\beta _{c})_{\ast }$, we suggest the following expression (figure 8(b), dashed line):

(5.13)\begin{align} \check{Q}_{gl} \approx \left\{ \begin{array}{@{}ll} (1 - {\alpha}/{\alpha_{c}})^{{1}/({p+1})} - \dfrac{1}{p+1}(1 - {\alpha}/{\alpha_{c}})^{({p-1})/({p+1})} {\beta}/{\beta_{c}}+(-{\beta}/{\beta_{c}})^{{1}/{p}}, & \text{for ${\beta}/{\beta_{c}}<0$}, \\ {(1 -{\alpha}/{\alpha_{c}})^{{1}/({p+1})}}[{1 + (1 - {\alpha}/{\alpha_{c}})^{-{p}/({p+1})}{\beta}/{\beta_{c}}}]^{{-}1}, & \text{for ${\beta}/{\beta_{c}}\ge 0$}. \end{array} \right. \end{align}

This expression can then be used to obtain the new flux condition (5.6), which is still approximately valid for values of $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ which are not small.

Figure 8. Effect of $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ on $\check {Q}_{gl}$, for a non-vanishing Budd friction law with $p = 1/3$. The coloured continuous lines are obtained by solving numerically (5.2). The dashed black line is obtained using (5.13). (a) Various values of $\alpha / \alpha _c$ and (b) zoom on the case $\alpha / \alpha _c = 0.25$.

5.2. Vanishing friction law: non-negligible membrane-stress divergence

We now consider the case where the divergence of membrane stress cannot be neglected in the momentum-balance equation. Specifically, we consider the Budd friction law combined with the N$_\text {A}$ effective-pressure model. In that case $N_{gl}=0$, so it does not make sense to use the reference flux $q_{{gl,c}}$ defined in (5.3). Instead, we define it as

(5.14) \begin{align} q_{{gl,c}} = \left(\dfrac{1-\rho/\rho_{w}}{8}\right)^{{n}/({p+1})} (\rho g)^{-({q-1})/({p+1})}(2\rho g)^{{n}/({p+1})}C^{-{1}/({p+1})}A^{{1}/({p+1})} h_{gl}^{({n+(p-q)+3})/({p+1})}. \end{align}

Note that, in contrast to the previous subsection, this is not the expression of the flux that was derived in § 3. It is rather a reference flux that is used to define $\beta /\beta _{c}$, without any specific physical interpretation. Again, we include the effect of the assumptions into a prefactor $\check {Q}_{gl}$ such that

(5.15) \begin{align} q_{gl} &= \check{Q}_{gl}\left(\dfrac{1-\rho/\rho_{w}}{8}\right)^{{({n-q})/({p+1})}} (\rho g)^{-({q-1})/({p+1})}(2\rho g)^{{n}/({p+1})}C^{-{1}/({p+1})} \nonumber\\ & \quad \times A^{{1}/({p+1})} h_{gl}^{({n+(p-q)+3})/({p+1})}, \end{align}

with $\check {Q}_{gl}=\check {Q}_{gl}({\alpha }/{\alpha _{c}},{\beta }/{\beta _{c}})$. The previous discussion holds if the divergence of membrane stress can be neglected in the momentum-balance equation. In general, and in particular for the Budd friction law with the N$_\text {A}$ effective-pressure model, that is not the case. Still, we can follow a strategy similar to the one used in § 4 to derive the flux conditions of hybrid friction laws to take into account the effect of ${\alpha }/{\alpha _{c}}$ and ${\beta }/{\beta _{c}}$: we can treat these ratios as parameters of the problem, and consider a mapping of the form $(\tilde {\alpha }, \tilde {\beta }) \mapsto \tilde {Q}_{gl}(\tilde {\alpha }, \tilde {\beta })$. More precisely, if we keep the terms associated with the net mass accumulation rate and the bedrock slope in the derivation of the flux condition described in § 3, we obtain the following system of equations, in place of (3.7):

(5.16a–d) \begin{equation} \left\{ \begin{gathered} &\hphantom{00000000000000\;}\dfrac{\mathrm{d} \tilde{U}}{\mathrm{d} \tilde{X}} ={-} \vert \tilde{W} \vert^{n-1} \tilde{W}, \quad\text{for } 0 < \tilde{X} < \tilde{Q}_{gl}/\tilde{\alpha},\\ &\dfrac{\mathrm{d} \tilde{W}}{\mathrm{d} \tilde{X}} ={-}\dfrac{1}{4}\dfrac{\tilde{U}}{{\tilde{Q}{_{gl}}-\tilde{\alpha}\tilde{X}}} \left(\dfrac{{\tilde{Q}{_{gl}}-\tilde{\alpha}\tilde{X}}}{\tilde{U}} - 1_{\text{A}} \left\langle 1 + \dfrac{\tilde{\beta}}{1-\delta}\tilde{X}\right\rangle\right)^{q}\nonumber\\ & \hphantom{00000\,}\times\, \vert\tilde{U}\vert^{p-1}\tilde{U} -\dfrac{\vert \tilde{W}\vert^{n+1}}{\tilde{U}}+\tilde{\alpha}\dfrac{\tilde{W}}{\tilde{Q}_{gl}-\tilde{\alpha} \tilde{X}}-\dfrac{1}{4}\dfrac{\tilde{\alpha}}{\tilde{U}} \nonumber\\ &\hphantom{00000\;} +\dfrac{({\tilde{Q}{_{gl}}-\tilde{\alpha}\tilde{X})}\vert \tilde{W}\vert^{n-1}\tilde{W}}{4 \tilde{U}^2} - \dfrac{\tilde{\beta}}{4}, \quad\text{for } 0 < \tilde{X} < {\tilde{Q}_{gl}}/{\tilde{\alpha}}, \\ &\hphantom{000000000000\,\;}(\tilde{U},\tilde{W})=({\tilde{Q}{_{gl}}},{\delta}/{8}), \quad\text{at } \tilde{X} = 0,\\ &\hphantom{00000000000000000000000\;\;}\tilde{U}=0, \quad\text{at } \tilde{X} = \tilde{Q}_{gl}/\tilde{\alpha}, \end{gathered}\right.\end{equation}

with

(5.17)\begin{equation} \tilde{\alpha} = \left(\dfrac{\delta}{8}\right)^{n}\dfrac{\alpha}{\alpha_{c}} \quad \text{and} \quad \tilde{\beta} = \left(\dfrac{\delta}{8}\right)^{{np}/({p+1})} \dfrac{\beta}{\beta_{c}}. \end{equation}

This system of equations is fundamentally different from (3.7). Indeed, it is formally equivalent to the initial system of equations presented in § 2, for unbuttressed ice sheets, since no additional assumption has been made. By contrast, the system of (3.7) used in § 3 to obtain the flux conditions was only valid within the boundary layer near the grounding line and in the limit of $\varepsilon \rightarrow 0$. The system (5.16) is also more complex in two respects. On the one hand, the dynamical system defined by (5.16a) and (5.16b) is non-autonomous, since $\tilde {X}$ appears in the definition of $\text {d}\tilde {W}/\text {d}\tilde {X}$. On the other hand, this system depends on the additional parameters $\tilde {\alpha }$ and $\tilde {\beta }$. Because $\tilde {\beta }$ is proportional to the bedrock slope $\text {d}b/\text {d}x$ which depends on the $x$ coordinate, in general, $\tilde {\beta } = \tilde {\beta }(\tilde {X})$.

However, the analysis can be simplified by considering linear bed geometries, so that $\tilde {\beta }$ is constant. Let us fix the values of both $\tilde {\alpha }$ and $\tilde {\beta }$. The system of (5.16) is then a parametric system which only possesses solutions for specific values of $\tilde {Q}_{gl}$. Despite the differences that have been mentioned, we have found that the shooting method introduced in § 3 and described in Appendix B was still applicable to the system (5.16). We can thus obtain these particular values $\tilde {Q}_{gl}$. Then, we convert the mapping $(\tilde {\alpha }, \tilde {\beta })\mapsto \tilde {Q}_{gl}(\tilde {\alpha }, \tilde {\beta })$ back the mapping $(\alpha /\alpha _{c}, \beta /\beta _{c})\mapsto \check {Q}_{gl}(\alpha /\alpha _{c}, \beta /\beta _{c})$ by using (5.17) and $\tilde {Q}_{gl} = (\delta /8)^{(n-q)/(p+1)}\check {Q}_{gl}$, which was derived in § 3. We have represented the effect of $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ on $\check {Q}_{gl}$ in figure 9(a) using the aforementioned numerical method.

Figure 9. Effect of $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ on $\check {Q}_{gl}$, for the Budd friction law with the effective-pressure model N$_\text {A}$, $n=1/3, p = 1/3, q=1$, and $\delta = 0.1$. The coloured continuous lines are obtained by finding the values of $\tilde {Q}_{gl}$ that yield a solution to (5.16). The dashed black line is obtained using (5.18). (a) Various values of $\alpha / \alpha _c$ and (b) zoom on the case $\alpha / \alpha _c = 0.25$.

In contrast to the case of non-vanishing friction laws, it is not easy to derive asymptotic expressions for $\check {Q}_{gl}$ for large or small values of $\beta /\beta _{c}$, as one has to solve (5.16), which is significantly more complex than an algebraic equation. Instead, we parametrise $\check {Q}_{gl}$ using a curve-fitting approach with simple expressions. We suggest the following expression (figure 9(b), dashed line):

(5.18)\begin{equation} \check{Q}_{gl} \approx \left\{ \begin{array}{@{}ll} \check{Q}_{gl}^0(1 - 3.72\,{\beta}/{\beta_{c}}), & \quad \text{for ${\beta}/{\beta_{c}}<0$}, \\ \check{Q}_{gl}^0[1 + 17.76(1-{\alpha}/{\alpha_{c}})^{{-}1} {\beta}/{\beta_{c}}]^{{-}1}, & \quad \text{for ${\beta}/{\beta_{c}}\ge 0$}, \end{array} \right. \end{equation}

with $\check {Q}_{gl}^0\equiv \check {Q}_{gl}\vert _{(\alpha /\alpha _{c},\beta /\beta _{c})=(0,0)} = 0.71$.

5.3. Non-vanishing friction law with $\gamma \ll 1$

Sergienko & Wingham (Reference Sergienko and Wingham2019) have considered flux conditions for the Weertman friction law in a regime of low basal and gravity stress. Specifically, they considered $\varepsilon \sim \delta \sim \gamma \ll 1$, leading to the divergence of membrane stress being of the same order as the friction stress, but much smaller than the gravity stress. This is a different regime from ours: in § 3 we have assumed that $\gamma \sim 1$ and considered a scaling that is such that the divergence of membrane stress, the friction stress and the gravity stress have the same order of magnitude.

They have obtained, as a zeroth-order solution, the following expression:

(5.19)\begin{equation} q_{gl}\left(\dfrac{\text{d}b}{\text{d}x}\right)_{gl} + a (1-\delta)h_{gl} = \left(\dfrac{1}{4}\rho g \left(1 - \dfrac{\rho}{\rho_{w}}\right)\right)^{n} A[(1-\delta)h_{gl}]^{n+2}. \end{equation}

In the limit $\delta \ll 1$, this equation becomes

(5.20)\begin{equation} q_{gl} = \left[1-\dfrac{a}{\left(\dfrac{1}{4}\rho \left(1 - \dfrac{\rho}{\rho_{w}}\right)g\right)^{n}A h_{gl}^{n+1}} \right] \left[\dfrac{\left({\text{d}b}/{\text{d}x}\right)_{gl}}{\left(\dfrac{1}{4}\rho \left(1 - \dfrac{\rho}{\rho_{w}}\right)g\right)^{n}A h_{gl}^{n+2}}\right]^{{-}1}. \end{equation}

This is exactly our equations (5.10) and (5.12), i.e. this flux condition can be associated with the regime $\vert \beta /\beta _{c}\vert \gg 1$ of a Budd friction law which does not vanish at the grounding line and in which the membrane-stress divergence is negligible. This scaling can be motivated by (5.5): $\vert \beta /\beta _{c}\vert \propto \gamma ^{-1/(p+1)}$.

6.1. Set-up

The values chosen for the physical parameters are typical for numerical experiments with marine ice sheets, and are similar to the ones presented in Pattyn et al. (Reference Pattyn2012). We take $n=3, \rho =900\ {\rm kg}\ {\rm m}^{-3}$, $\rho _{w}=1000\ {\rm kg}\ {\rm m}^{-3}$ and $g=9.8\ {\rm m}\ {\rm s}^{-2}$. Glen's viscosity parameter is set to $A = 4.9\times 10^{-25}\ {\rm Pa}^{-3}\ {\rm s}^{-1}$, and the net mass accumulation rate is set to $a=9.51\times 10^{-9}\ {\rm m}\ {\rm s}^{-1}$. In terms of the friction laws, we consider the (W), (C), (B), (T), (RC1) and (RC2) friction laws with $p=1/3$ and $q=1$, and with both the N$_\text {A}$ and the N$_\text {B}$ effective-pressure models. The friction coefficient for the (W) friction law is set to $C = 7.624\times 10^6\ {\rm Pa}\ {\rm m}^{-1/3}\ {\rm s}^{1/3}$. For the other friction laws, the friction coefficient will be specified for each specific numerical experiment. The hydrology parameter $c$ is set to $0.96$. Three bed elevation profiles are considered (figure 10). The first one is a polynomial bed that will be used to compare the flux conditions in an idealised configuration. The second one depends linearly on $x$ and will be used to check the effect of the bed slope (and thus of $\beta$) on the flux conditions. The third one is similar to the linear one, but an oscillatory signal has been added on top of it. It will be used to investigate the effect of local variability in the bedrock profile.

Figure 10. Bed profiles considered in the numerical experiments. (a) Polynomial bed, (b) linear bed and (c) linear bed with oscillations.

Results are obtained either from the flux conditions themselves, or from the numerical solution of the initial problem ((2.1)–(2.6)). For the spatial discretisation, we use an in-house finite-element code. The mesh is uniform with a constant element size of 180 m.

6.2. Flux conditions for the Budd and hybrid friction laws

The first experiment compares the flux conditions obtained in §§ 3 and 4 with results of numerical simulations. It mimics the experiment 3 of the Marine Ice Sheet Model Intercomparison Project (Pattyn et al. Reference Pattyn2012), which is a benchmark for the comparison of marine ice-sheet flowline models. We considered the polynomial bed profile (figure 10a), fixed all the parameters to their reference values, except for the ice rheology parameter $A$ which is varied. For each particular value of $A$, a steady-state ice-sheet solution was obtained and the grounding-line position was retrieved. On the one hand, this position was obtained numerically, thanks to the finite-element solution. On the other hand, we computed the grounding-line position from the flux conditions: from the mass-conservation equation, we have the global balance

(6.1)\begin{equation} q_{gl}(h_{gl}) = a\,x_{gl}, \end{equation}

where we have written $q_{gl}=q_{gl}(h_{gl})$ to emphasise the dependency on the grounding-line ice thickness. The flotation condition $h_{gl} = - (\rho _{w}/\rho ) b(x_{gl})$ then allowed us to obtain an algebraic equation for $x_{gl}$:

(6.2)\begin{equation} q_{gl}(-({\rho_{w}}/{\rho})b(x_{gl})) = a\,x_{gl}. \end{equation}

We solved this nonlinear equation using a Newton–Raphson procedure.

It remains to choose the values of the friction coefficients for all the friction laws except for the Weertman one. This is quite delicate, because the friction coefficients associated with different friction laws are not necessarily comparable to one another; in particular, they do not have the same dimensions. For the Weertman friction law, (6.2) has a solution $x_{gl}\approx 800$ km for $A=10^{-25}\ {\rm Pa}^{{-3}}\ {\rm s}^{{-1}}$. We then chose the friction coefficients $C$ for the Coulomb friction law and the Budd friction law so as to obtain this solution as well. The obtained friction parameters are shown in table 2. For the hybrid friction laws, we considered the same friction coefficient $C$ as the one obtained for the Coulomb friction law because the Coulomb friction law is a limit case of the hybrid friction laws. The coefficient $A_s$ was chosen such that $A_s^{-p}$ had the same value as the Weertman friction coefficient, again by identification of the hybrid friction law as a Weertman friction law. Finally, we considered $u_0 = 10^{-5}\ {\rm m}\ {\rm s}^{-1}$, which is a typical value for the velocity in marine ice sheets. All these values are summarised in tables 2 and 3.

Table 2. Numerical values of the friction coefficients used for the verification of the flux conditions.

Table 3. Numerical values of the additional friction parameters $A_s$ and $u_0$ used for the verification of the flux conditions.

The results are shown in figure 11. The grounding-line positions obtained using the flux conditions match the results from the numerical simulations. We note that the physical parameters and the bed profile considered in this numerical experiment are consistent with the assumptions made during the derivation of the flux conditions, namely, the net mass accumulation rate and the bedrock slopes are not too large, and the friction coefficient is not too small. With respect to the discussion of § 5, the experiments have been conducted in a regime for which $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ are small.

Figure 11. Comparison of the evolution of the grounding-line position with respect to $A$ for the different friction laws and effective-pressure models, using the flux condition (lines) and results of a finite-element discretisation of the original problem (circles, crosses). The results for the N$_\text {A}$ and N$_\text {B}$ effective-pressure models are respectively shown in blue and in green.

As a side note, it can be observed that the curves all have the same shape, which could suggest that the choice of friction laws actually has little impact on the mechanical equilibrium of marine ice sheets, and in particular on flux conditions. However, this similarity is not the result of the impact of friction laws but rather stems from the methodology used. The flux conditions associated with different friction laws differ in two aspects: the exponent on the grounding-line thickness, and the dependence of the factor that multiplies this thickness with respect to the physical parameters ($A, C,\ldots$). The considered bedrock does not show a strong variability, so that the exponent on top of the grounding-line thickness has a limited effect. Moreover, by construction, the friction coefficients were chosen uniformly and in such a way that the curves pass through the same point, which effectively leads to a similar factor in front of the grounding-line thickness. This explains the similarity between the curves shown in figure 11. In practice, however, the friction coefficients are not uniform, but, rather, are tuned spatially so as to obtain a similar thickness and velocity profile compared with some observations. This results in very different dynamics. We refer interested readers to Brondex et al. (Reference Brondex, Gagliardini, Gillet-Chaulet and Durand2017) where these differences are discussed.

6.3. Effect of $\alpha, \beta$ and $\gamma$

We now conduct a series of numerical experiments to determine numerically the situations in which the assumptions made to derive the flux conditions in § 3 are not valid, and to confirm that the new expressions, namely (5.6) and (5.15) combined respectively with the corrections factors defined in (5.13) and (5.18), can be applied to correct these flux conditions. Practically, we check that they lead to the same grounding-line flux value as the numerical results. We call the flux conditions derived in § 5 ‘enriched’ flux conditions. First, we consider the linear bed profile (figure 10b), whose elevation is given by $b(x) = b_0 + b_1(x/L)$ with $b_0 = 720$ m, $b_1=-900$ m and $L=750$ km. We vary three physical parameters: the net mass accumulation rate $a$, the bedrock slope $\text {d}b/\text {d}x$ and the friction coefficient $C$. The goal is to reach a regime in which $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ are not small so that the flux conditions derived in § 3 are not valid anymore. Then, we consider the more realistic ‘rough’ bedrock profile, as well as different values for the friction coefficient. We always consider the Budd friction law with both the N$_\text {A}$ and N$_\text {B}$ effective-pressure models. We choose a reference friction coefficient of $C_{0}^{\text {A}}=1.73$ m$^{-1/3}$ s$^{1/3}$ in the first case and $C_{0}^{\text {B}} = 43.22$ m$^{-1/3}$s$^{1/3}$ in the second case. These values where chosen such that $C_{0}^{\text {A}}\rho g h_{gl} \approx C_{0}^{\text {B}}(1-c)\rho g h_{gl} \approx 7.624\times 10^6\,$Pa m$^{-1/3}$s$^{1/3}$ for $h_{gl}=500$ m.

First, we consider the reference physical parameters previously introduced, and we modify the values of $a, \text {d}b/\text {d}x$ and $C$ in the following way. We first consider $a$, and vary its value within the interval $a_0 \le a \le 10 a_0$, where $a_0$ is the reference value introduced in the set-up subsection. For each fixed value of $a$, we let the ice sheet evolve until it reaches a steady state. This leads to a collection of grounding-line fluxes, which are compared with the grounding-line fluxes that would have been obtained thanks to our flux conditions. For the N$_\text {A}$ effective-pressure model, we use (5.15) combined with (5.18), while for the N$_\text {B}$ effective-pressure model, we use (5.6) combined with (5.13). We then perform a similar procedure for $\text {d}b/\text {d}x$ and $C$, which are respectively varied in the ranges $10 (\text {d}b/\text {d}x)_0\le \text {d}b/\text {d}x \le (\text {d}b/\text {d}x)_0$ and $0.5\times 10^{-1}C_0 \le C \le C_0$, with $(\text {d}b/\text {d}x)_0 = b_1/L$. In this former case, only the slope of the linear bed is varied; the value $b(0)=b_0$ is left unchanged. By increasing the value of $a$, of $\vert \text {d}b/\text {d}x \vert$ and reducing the value of $C$, we attempt to reach a regime in which $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ cannot be neglected. The results are shown in figure 12. It can be observed that, for the parameters considered, the ratio $q_{{gl,fc}}/q_{gl}$ stays close to one when the N$_\text {B}$ effective-pressure model is used, even when we use the flux condition derived in § 3. By contrast, this ratio departs significantly from one when the slope or the friction parameter are varied in a simulation in which the N$_\text {A}$ effective-pressure model is considered. That is not the case if we use the enriched flux conditions, as those lead to a ratio that is always close to one.

Figure 12. Comparison between the fluxes $q_{{gl,fc}}$ obtained thanks to the flux conditions derived in § 3 (circles) and thanks to the enriched flux conditions (crosses), and the grounding-line fluxes $q_{gl}$ obtained numerically, when $a/a_0, {(\text {d}b/\text {d}x)}/{(\text {d}b/\text {d}x)_0}$ and $C/C_0$ are varied (first line). Ratios $\alpha /\alpha _{c}$ and $\beta /\beta _{c}$ corresponding to each numerical solution (second line). We have considered the Budd friction law with the N$_\text {A}$ (blue) and the N$_\text {B}$ (green) effective-pressure models.

In practice, we expect a relatively variable bedrock elevation; hence, a linear configuration might not be appropriate. To investigate the impact of this bedrock variability, we consider the bedrock profile shown in figure 10(c). Its elevation is given by $b(x) = b_0 + b_1(x/L) + b_2\sin (2{\rm \pi} x/L_{o})$, where $b_0, b_1$ and $L$ have the same values as before, $b_2=300$ m, and where $L_{o}$ is varied between $100$ and $300$ km. The physical parameters are the same as the ones used previously when varying the net mass accumulation rate $a$. We observe in figure 13 similar findings compared with the previous numerical experiment. Firstly, the ratio $q_{{gl,fc}}/q_{gl}$ calculated using the flux conditions derived in § 3 deviates further from a unit value as the bedrock has a larger slope variation. Secondly, the effect is much more pronounced in the case of the N$_\text {A}$ effective-pressure model. Lastly, the use of corrective factors in flux conditions enables satisfactory results, namely a $q_{{gl,fc}}/q_{gl}$ ratio that remains close to unity.

Figure 13. Effect of local variability in the geometry profile on the ratio between the fluxes $q_{{gl,fc}}$ obtained based on the flux conditions derived in § 3 (circles) and on the enriched flux conditions (crosses), and the flux $q_{gl}$ obtained numerically. We have considered the Budd friction law with the N$_\text {A}$ (blue) and the N$_\text {B}$ (green) effective-pressure models; (a) $L_o=100$ km, (b) $L_o=200$ km and (c) $L_o=300$ km.