Introduction
Determining the location and the rate of migration of the grounding line – a transition between the grounded part of a marine ice sheet and its floating ice shelf – is a fundamental question in understanding ice-sheet interactions with the rest of the climate system. In circumstances in which the characteristic scales of ice-sheet internal processes are short in comparison with those of wider climate variability, it is possible to consider ice sheets as an isolated system and, through investigation of the stability of their steady-state configurations, learn about their intrinsic processes, specific to ice dynamics, and unrelated to the external climate variability. Weertman (Reference Weertman1974) was the first to provide a theoretical description of stability of steady-state marine ice sheets. He concluded that ‘A stable ice sheet can occur if the bed slopes away from the center of the ice sheet. The generalization of our results to other bed shapes is rather obvious.’ The potential practical importance of this statement has resulted in a number of more detailed investigations (e.g., Chugunov and Wilchinsky, Reference Chugunov and Wilchinsky1996; Schoof, Reference Schoof2007b; Nowicki and Wingham, Reference Nowicki and Wingham2008; Tsai and others, Reference Tsai, Stewart and Thompson2015; Favier and others, Reference Favier, Pattyn, Berger and Drews2016; Haseloff and Sergienko, Reference Haseloff and Sergienko2018). In particular, using matched asymptotic expansions, Schoof (Reference Schoof2007b) has developed a boundary layer theory for an ice stream undergoing rapid sliding at its grounding line. In circumstances where the grounded ice sheet and floating ice shelf are laterally unconfined, the steady-state ice flux at the grounding line is a single-valued monotonically increasing function of ice thickness at the grounding line (Schoof, Reference Schoof2007b; Tsai and others, Reference Tsai, Stewart and Thompson2015). Using the results of this theory Schoof (Reference Schoof2012) performed a linear stability analysis, whose result confirmed Weertman's heuristic argument that the grounding lines on retrograde beds (the bed elevation increases in the direction of ice flow) are unstable. This outcome is widely known as the ‘marine ice-sheet instability hypothesis’.
While these studies have considered two distinct types of flow of the grounded ice, one dominated by vertical shear (e.g., Weertman, Reference Weertman1974; Nowicki and Wingham, Reference Nowicki and Wingham2008; Wilchinsky, Reference Wilchinsky2009), and the other one dominated by basal sliding (e.g., Schoof, Reference Schoof2007b; Tsai and others, Reference Tsai, Stewart and Thompson2015), their common assumptions are that the driving stress is of the same order as the basal stress, and that both are larger than the divergence of the longitudinal stress. At the grounding line the longitudinal stress must balance the hydrostatic pressure. Satisfying these two balances led Schoof (Reference Schoof2007b) to the conclusion that the ice thickness at the grounding line must be small in comparison with the ice thickness in the interior, with the consequence that the ice thickness gradient is large in the vicinity of the grounding line. However, this is in marked contrast to the observed low surface gradients (Rose, Reference Rose1979; MacAyeal and others, Reference MacAyeal, Bindschadler and Scambos1995) of many present-day West Antarctic ice streams, particularly those which flow to the Ross Ice Shelf (Fig. 1). Tsai and others (Reference Tsai, Stewart and Thompson2015) have proposed that the absence of the large surface gradient is a consequence of the Coulomb friction regime in the well-lubricated environment near the grounding line, but these gradients, which are of the order of 1–5·10−3 (Fig. 1a), extend over hundreds of kilometers and do not exhibit significant changes in the vicinity of the grounding line (Fig. 1c). This suggests that these ice streams are characterized by extremely weak beds over extended distances (MacAyeal, Reference MacAyeal1989), and this leads us to the view that analyzing the stability of these ice streams requires a reconsideration of the balance of the terms in the momentum balance.
Fig. 1. Characteristics of Siple Coast ice streams. (a) Magnitude of the surface gradient 
$\vert \vec \nabla S\vert $; (b) magnitude of the driving stress (kPa) computed as 
$\tau_d =\rho gH \vert \vec\nabla S \vert $ derived from BEDMAP 2 data set (Fretwell and others, Reference Fretwell2013); (c) surface (blue line) and bed (red line) elevation along a flowline on Willans Ice Plain (red line in panel (b)). Siple Coast ice streams: McIS – MacAyeal Ice Stream; BIS – Bindschadler Ice Stream; WIS – Whillans Ice Stream. Their driving stress is in ~ 0–20 kPa limit. Inset shows the map of Antarctica (Haran and others, Reference Haran, Bohlander, Scambos and Fahnestock2005), black rectangle shows the location of the Siple Coast region.
Using analytical and numerical models, we examine the same configuration of an unconfined ice stream floating into an unconfined ice shelf as considered by Schoof (Reference Schoof2007b) and Tsai and others (Reference Tsai, Stewart and Thompson2015), but extend the analysis to a regime in which the basal shear is on the same order of magnitude as the divergence of extensional stress through the whole length of the ice stream. Our results show that steady-state configurations of marine ice streams are substantially different from those in the regime of strong basal and driving stress. In the regime of low or absent basal shear, ice surface slopes, and consequently the driving stress, are small through the length of the ice stream, and not only in the vicinity of the grounding line. Their flow is controlled by the effects of the bed geometry – elevation, slopes and curvature – that here we collectively term ‘form drag’, basal drag, the accumulation rate, its gradients and higher derivatives. We find that at the grounding line, the ice flux is a multivalued implicit function of ice thickness. Our analysis show that the stability of distinct steady-state configurations is determined by the sign of a parameter denoted Λ that depends on form drag, basal shear, accumulation rate and their first and higher derivatives. In contrast to the marine ice-sheet instability hypothesis that suggests that, in the case of net accumulation, ice sheets are stable on prograde beds only (the bed elevation decreases in the direction of ice flow), ice streams experiencing low basal shear can be both stable (Λ <0) and unstable (Λ >0) on prograde and retrograde beds. We also find that the reciprocal of the stability parameter Λ can be interpreted in terms of the e-folding time of perturbed steady-state configurations, and thus |Λ|−1 indicates how fast or slow this perturbed configuration returns to or evolves from, their steady states. It provides a measure as to whether the evolution of the perturbed configurations can be viewed as insensitive to climate variability.
Model description
We consider an unconfined ice stream that flows into an unconfined ice shelf (Fig. 2) and use the model formulation of Schoof (Reference Schoof2007b), which we briefly describe below. Under assumptions of negligible vertical shear appropriate for ice stream and ice shelf flow (MacAyeal, Reference MacAyeal1989), the vertically integrated momentum balance is
$$ \eqalign{& 2\left( {A^{-(1/n)}h{\left| {u_x} \right|}^{(1/n)-1}u_x} \right)_x-\tau _{\rm b}-\rho gh\left( {h + b} \right)_x = 0, \cr & \quad x_d \le x \le x_j,} $ $$
$$ 2\left(A^{-({1}/{n})}h\left\vert u_x\right\vert ^{({1}/{n})-1}u_x\right)_x -\rho g^\prime hh_x=0, \quad x_j\leq x\leq x_c, $$where u(x) is the depth-averaged ice velocity, h(x) ice thickness, b(x) is the bed elevation (negative below sea level and positive above sea level), A −1/n is the ice stiffness parameter (assumed to be constant), n is an exponent of Glen's flow law, g is the acceleration due to gravity, τb is basal shear, g′ is the reduced gravity defined as
$$ g^\prime=\delta g $$where
$$ \delta = {\rho_{\rm w}-\rho \over \rho_{\rm w}} $$is the buoyancy parameter, ρ and ρ w are the densities of ice and water, respectively, x d is the location of the ice divide (assumed to be x d = 0), x c is the location of the calving front and x j is the location of the junction between the ice stream and ice shelf. In the absence of a universal ‘sliding law’ – a relationship between shear at the ice/bed interface and characteristics of ice and bed – several formulations are commonly used. Throughout the derivations presented below we adopt a formulation used by Schoof (Reference Schoof2007b)
$$ \tau_{\rm b}=C\left\vert u\right\vert ^{m-1}u, $$where C and m are constant parameters. The first term in (1a) is the divergence of the extensional, or longitudinal stress, τx, the second term is the basal shear, τb, and the final terms in the both equations are the driving stress, τd.
Fig. 2. Model geometry: b – bed elevation (b <0), h – ice thickness, x d – the ice divide location, x g – the grounding line location; x c – the calving front location.
The mass balance is
$$ h_{\rm t}+\left(uh\right)_x =\left\{ \matrix{ \dot a & 0\leq x\leq x_{\rm j} \cr \dot m & x_{\rm j}< x\leq x_{\rm c},} \right. $$where 
$\dot a$ and 
$\dot m$ are the net accumulation/ablation (positive for accumulation) rates of the ice stream and ice shelf, respectively. (The discontinuity at the junction x j is due to the fact that basal ablation by the ocean is possibly discontinuous across the junction.)
The boundary conditions are
$$ (h+b)_x =0,\quad u=0, \quad x=0 $$
$$ 2 A^{-({1}/{n})}h\left\vert u_x\right\vert ^{({1}/{n})-1}u_x={1 \over 2}\rho g^\prime h^2, \quad x=x_{\rm c}. $$The condition at the calving front, x c, requires the longitudinal stress in the ice to balance the pressure deficit at the front. The conditions at the divide, x d, are statements that there is no driving stress there, and that no flow enters or leaves the domain from the left. These are the natural conditions to describe a divide. A detailed treatment near the divide in fact requires the solution of a full-Stokes problem (Fowler, Reference Fowler1992). However, the primary concern of this paper is not the detailed behavior near the divide, and we will use (1)–(6), as a sufficient description.
At the junction, x j continuity conditions have to be prescribed. The number of conditions is determined by the number of the boundary conditions for the ice stream and the ice shelf.
$$ u_{{\rm stream}}(x_{\rm j})=u_{{\rm shelf}}(x_{\rm j}) $$
$$ h_{{\rm stream}}(x_{\rm j})=h_{{\rm shelf}}(x_{\rm j}) $$
$$ \tau_{{\rm stream}}(x_{\rm j})=\tau_{{\rm shelf}}(x_{\rm j}). $$The problems (1)–(7) provide solutions for any given location of the junction x j. However, the solution h(x j) must in addition satisfy the flotation condition
$$ h(x_{\rm j}) = - {\rho_{\rm w} \over \rho} b(x_{\rm j}). $$The grounding line location x g is then determined as the root of this equation in x j. We will use the symbols h g, u g and τg to mean ice thickness, velocity and the longitudinal stress at the grounding line.
Additionally, the following inequalities have to be satisfied
$$ h(x) \gt -{\rho_{\rm w} \over \rho}b(x), \quad 0 \lt x \lt x_{\rm g}, $$
$$ h(x)\leq - {\rho_{\rm w} \over \rho}b(x), \quad x_{\rm g}\leq x \lt x_{\rm c}. $$These conditions reflect the fact that the ice stream is grounded and the ice shelf is floating. Similar to all previous studies of unconfined ice shelves (e.g., MacAyeal and Barcilon, Reference MacAyeal and Barcilon1988; Schoof, Reference Schoof2007b) we integrate the momentum balance of the ice shelf (1b), and using the boundary condition at the calving front, x c, (6b) and the continuity conditions at the grounding line (7), and reduce the problem to the ice stream part only, with the boundary conditions at the grounding line – the flotation condition (8) and the stress condition
$$ 2 A^{-({1}/{n})}h\left\vert u_x\right\vert ^{({1}/{n})-1}u_x={1 \over 2}\rho g^\prime h^2, \quad x=x_{\rm g}. $$The model (1)–(10) is identical to one used by Schoof (Reference Schoof2007b).
In what follows we use precise numerical solutions of (1)–(10) to compare with approximate asymptotic expressions, and describe them as ‘exact’, acknowledging their numerical origin. To obtain the ‘exact’ solutions, we use the finite-element solver Comsol™, (COMSOL, 2018), whose accuracy is well-established (e.g. Schäfer and Turek, Reference Schäfer, Turek and Hirschel1996).
We now proceed to non-dimensionalize (1)–(10) using as characteristic scales the bathymetry of a continental shelf [b], the length of an ice steam [x] and its velocity [u]. In terms of these
$$ [h]=[b], \quad [\dot a] ={[b][u] \over [x]}, \quad [t] = {[x] \over [u]}. $$With the two non-dimensional parameters
$$ \varepsilon= {2A^{-({1}/{n})}[u]^{1/n} \over \rho g[h][x]^{1/n}} $$
$$ \gamma = {[C][u]^m[x] \over \rho g[h]^2}, $$relating to deformation and sliding, respectively, the non-dimensional momentum balance and mass balance become
$$ \varepsilon\left(h\left\vert u_x\right\vert ^{({1}/{n})-1}u_x\right)_x - \gamma\left\vert u\right\vert ^{m-1}u-h(h+b)_x=0, \quad 0\leq x\leq x_{\rm g}, \quad h\leq- {b \over 1-\delta} $$and
$$ h_{\rm t}+\left(uh\right)_x = \dot a, \quad 0\leq x\leq x_{\rm g}. $$The non-dimensional boundary conditions at x = 0 and x = x g are
$$ (h+b)_x =0,\quad u=0, \quad x=0 $$
$$ \varepsilon h\left\vert u_x\right\vert ^{({1}/{n})-1}u_x={1 \over 2}\delta h^2, \quad x=x_{\rm g}, $$
$$ h=-{b(x) \over 1-\delta} \quad x=x_{\rm g}. $$The inequality constraint (9) is
$$ h(x) \gt -{b(x) \over 1-\delta}, \quad 0 \lt x \lt x_{\rm g}. $$The non-dimensional problem (13)–(16) has three non-dimensional parameters: ε, expression (12a), γ, expression (12b) and δ, expression (3). Their magnitudes are determined by the characteristic scales of the problem. The buoyancy parameter δ, defined by (3), is determined by the values of densities of ice, ρ = 917 kg m−3, and sea-water, ρ w ≈1020 kg m−3, which results in δ ≈ 0.1. To determine the value of ε, we take as the characteristic scales for the present-day Antarctic ice streams to be D = 1000 m, L = 500 km and V = 500 m a−1 (Fretwell and others, Reference Fretwell2013) together with the ice rheological parameters n = 3 and A = 1.35 · 10−25 Pa−3 s−1 (that corresponds to the ice temperature T = −20°C). These then provide the value for parameter ε is ~1.4·10−2 ≪1.
The choice of γ is not determined by direct observations, because the basal shear, τb, is not directly observed, and must be inferred from the surface observations, e.g., by means of inversion techniques from surface observations (e.g. MacAyeal, Reference MacAyeal1992; Sergienko and others, Reference Sergienko, Bindschadler, Vornberger and MacAyeal2008; Sergienko and Hindmarsh, Reference Sergienko and Hindmarsh2013). We obtain a value of C based on the results of limited direct observations of the basal conditions (e.g., Engelhardt and others, Reference Engelhardt, Humphrey, Kamb and Fahnestock1990) and results of inverse modeling studies (e.g, MacAyeal, Reference MacAyeal1992; MacAyeal and others, Reference MacAyeal, Bindschadler and Scambos1995; Joughin and others, Reference Joughin, MacAyeal and Tulaczyk2004; Sergienko and others, Reference Sergienko, Creyts and Hindmarsh2014) indicating that basal shear under ice streams has small magnitudes. For instance, using Stokes' equations as a forward model, Sergienko and others (Reference Sergienko, Creyts and Hindmarsh2014) find that the inferred values of the basal shear in the trunks of MacAyeal and Bindschadler ice streams is ~0.5–5 kPa. Consequently, characteristic scales of the sliding parameter C are 2·104–2·105 Pa m−1/3 s1/3, and according to (12b), γ ~2.8·10−2 – 2.8 ·10−1.
We seek to provide an approximate solution of the problem (13)–(16), and this requires us to assign relative orders or the three parameters, ε, γ and δ.
As shown above ε and γ are naturally regarded as small, and we take ε ~ γ ~ o(1); but the choice of the order of δ is open to us. One choice is to follow Schoof (Reference Schoof2007b) and take δ as O(1). As we noted in section ‘Introduction’, this leads, as a consequence of the resulting imbalance of terms in the boundary condition (15b), to suppose the existence of a boundary layer close to x g in which the longitudinal and basal stress terms in the momentum balance (13) are of the same order as the driving stress (Schoof, Reference Schoof2007b; Tsai and others, Reference Tsai, Stewart and Thompson2015). The alternative is to take δ as a small parameter, and, indeed, with our choice of scales δ is closer to ε than to unity. We accordingly take δ ~ ε. This will remove the need for a boundary layer at the junction. The effectiveness of this assignment may be judged from the numerical results we describe later.
The problem can then be reduced to one in which δ is the sole parameter through rescaling in the following fashion using the parameter ε/δ:
$$ h=\left({\varepsilon \over \delta}\right)^\alpha H,\quad b=\left({\varepsilon \over \delta}\right)^\alpha B,\quad u=\left({\varepsilon \over \delta}\right)^\beta U,\quad x=\left({\varepsilon \over \delta}\right)^\kappa X,\quad t=\left({\varepsilon \over \delta}\right)^\sigma T. $$ Substituting these expressions into the momentum balance (13), mass balance (14) and the boundary conditions (15), balancing the first and second terms 
$\displaystyle {(1+\alpha -\kappa +(({\beta -\kappa })/{n})=1+m\beta )}$ in the momentum balance, all terms in the mass balance (α − σ = α + β − κ = 0), and the both sides in the stress (15b) condition at x g (1 + α + ((β − κ)/n) = 2α) we get
$$ \alpha=\sigma={n \over n+1},\quad \beta=-{1 \over (m+1)(n+1)}\left(=-{n \over n+1}\right),\quad \kappa={n(m+1)-1 \over (m+1)(n+1)}\left(=0\right), $$where expressions in parentheses are for the case of γ = 0, so that, in this case, x does not need rescaling. The rescaled problem becomes
$$ \delta\left(H\left\vert U_X\right\vert ^{({1}/{n})-1}U_X\right)_X -\delta \Gamma\left\vert U\right\vert^{m-1}U-H(H+B)_X=0, $$
$$ H_T+(UH)_X=\dot a $$
$$ (H+B)_X =0,\quad U=0, \quad X=0 $$
$$ H\left\vert U_X\right\vert ^{({1}/{n})-1}U_X={1 \over 2} H^2,\quad X=X_{\rm g} $$
$$ H=-{B \over 1-\delta} \quad X=X_{\rm g} $$
$$ H \gt -{B \over 1-\delta},\quad 0<X<X_{\rm g} $$where 
$\displaystyle {\Gamma ={\gamma }/{\varepsilon }}$ and the subscripts X and T are derivatives with respect to the coordinate X and time T. In contrast to the original problem (13)–(16) a small parameter δ appears only in the momentum balance (19a). This reduction, which was introduced by Schoof (Reference Schoof2007b, Appendix A) in connection with the corresponding problem in which the basal stress is of the same order as the driving stress, is exact: the solution for any pair of values of ε and δ may be obtained from the solution of the problem (19) through the scale assignment (17). However if, as we do here, one seeks an approximate solution of the problem (19) as a perturbation series expansion in the small parameter δ, one must formally assign δ ~ ε if the other terms in (19) are to be O(1).
To illustrate how different scales of basal shear affect the steady-state configurations we solve the problem (1)–(9) numerically (the details of numerical solver and used parameters are described in Appendix A) with two different values of the sliding parameter C. In the first experiment we use the same value as used by Schoof (Reference Schoof2007b) C = 7.6 · 106 Pa m−1/3 s1/3; in the second experiment, we use the value of C three orders of magnitude smaller C = 7.6 · 103 Pa m−1/3 s1/3. The computed ice-stream geometry and the momentum balance components for the two cases are markedly different (Fig. 3). In the case of a large value of the sliding coefficient, ice thickness is significantly larger in the interior than the magnitude of the bed elevation throughout the ice stream, and also than the ice thickness in the immediate vicinity of the grounding line. Consequently, the surface gradients are large in the vicinity of the grounding line (Fig. 3a). In contrast, in the case of low value of C, surface elevation neither varies substantially through the length of the ice stream nor in the immediate vicinity of the grounding line (Fig. 3d). Also, the ice thickness has the same order of magnitude as bed elevation, and for the chosen bed geometry, it is larger at the grounding line than in the interior of the ice stream.
Fig. 3. Surface elevation and components of the momentum balance for a high C = 7.6 · 106 Pa m−1/3 s1/3 (panels (a)–(c)) and low C = 7.6 · 103 Pa m−1/3 s1/3 (panels (d)–(e)) values of the sliding coefficient. Panels (a) and (d) surface (s) and bed (b) elevation (m), arrows show the ice-flow direction; (b), (c) and (e) components of the momentum balance, τx, τb and τd (kPa); insets in panel (b) show locations of a zoom shown in panel (c) τx for C = 7.6 · 106 Pa m−1/3 s1/3. The simulation parameters are the following: bed elevation b(x) = b 0 + b acos (πx/L), with b 0 = −500 m, b a = 250, and L=500 km; ice stiffness parameter A = 1.35 · 10−25 Pa−3 s−1 (corresponds to T ice ≈ −20°C); and accumulation rate 
$\dot a=0.94$ m a−1.
The momentum balance in the case of large C (Fig. 3b) is predominantly between the driving and basal stresses. These are of the order of 100–200 kPa throughout the length of the ice stream, and achieve their maximum at the grounding line. Although the divergence of the longitudinal stress increases near the grounding line (Fig. 3c), its magnitude remains significantly (two orders) smaller than the magnitudes of the driving and basal stresses. In the case of low C (Fig. 3e), the driving and basal stresses are much smaller than in the case of large C (less than 1 kPa), and the divergence of the longitudinal stress is of the same order as the basal shear (Fig. 3e), in line with the conclusions of MacAyeal (Reference MacAyeal1989) who found that ‘…basal drag and horizontal deviatoric stress gradients can contribute equally in balancing gravitational driving stress.’ These results suggest that the two regimes of high and low basal and driving stresses are very different; they also suggest that the power-law sliding law is capable of producing low surface gradients and low driving stress not only in the vicinity of the grounding line, but through the length of the ice stream, if the sliding coefficient C is small.
Steady states
In this section, we investigate steady solutions to (19) by setting H T in (19b) to zero. We treat δ as a small parameter and seek a solution as series expansion in δ, i.e.
$$ H\sim H_0+\delta H_1+\delta^2 H_2 +\cdots $$
$$ U\sim U_0+\delta U_1+\delta^2 U_2 +\cdots $$where the coefficients of powers of δ are supposed of O(1). Substituting H and U into the momentum-balance Eqn (19a), the mass-balance equation
$$ (UH)_X=\dot a $$and the boundary conditions (19c)–(19e), and equating the coefficients of each power of δ leads to various orders approximations of the problem (19a), (21) and (19c)–(19e). Detailed descriptions of the zeroth- and first-order solutions are provided in Appendix A. In this section we provide a brief summary.
The zeroth-order problem is
$$ H_0(H_0+B)_X=0 $$
$$ (U_0H_0)_X=\dot a $$
$$ U_0=0, \quad X=0 $$
$$ (H_0+B)_X=0, \quad X=0 $$
$$ H_0\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{0X}={1 \over 2}H_0^2,\quad X=X_{\rm g} $$
$$ H_0=-B, \quad X=X_{\rm g}. $$Its solutions, H 0 and U 0, are determined from (22a) together with (22f), and (22d)
$$ H_0(X)=-B(X) $$
$$ U_0(X)={\int_0^X\dot adx \over H_0(X)}={Q \over H_0}, $$where 
$Q=\int _0^X \dot a dx$ is ice flux.
Expression (23a) describes a solution in which the stream is neutrally buoyant: the thickness equals the bed depth and the surface elevation is zero. It provides no description of the roles of the divergence of the longitudinal stress and the basal shear stress in determining the position of the grounding line, for which we need to determine the first-order solutions. The zeroth-order solution does, however, provide an important insight. Using (23), the stress condition (22e) can be written as
$$ QB_{X}-\dot aB=\left({1 \over 2}\right)^n\left\vert B\right\vert ^{n+2},\ X=X_g, $$which provides an expression for the grounding flux. Q(X g) is not a function of B(X g) alone, as is the case with a high basal shear stress (Schoof, Reference Schoof2007b), but depends in addition on the basal gradient at the grounding line. This arises because, at this order, the surface is flat, and it is the basal gradient, rather than the surface gradient, that accounts for the ice thickness gradient. We will return to the implications of this observation having first obtained the first-order solutions.
The first-order problem is
$$ H_{1X}={1 \over H_0}\left[\left(H_0\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{0X}\right)_X -\Gamma\left\vert U_0\right\vert ^{m-1}U_0\right] $$
$$ (U_1H_0+U_0H_1)_X=0 $$
$$ U_1=0, \quad X=0 $$
$$ H_{1X}=0, \quad X=0 $$
$$ H_1\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{0X}+{1 \over n}H_0\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{1X}=H_1H_0, \quad X=X_{\rm g} $$
$$ H_1=-B, \quad X=X_{\rm g}. $$All the components of the first-order expressions are determined from the zeroth-order solutions U 0 and H 0 (23). Substitution of these solutions and their derivatives into (25a) gives
$$ \eqalign{ H_{1X} = &{\vert B\vert ^{-(({n+2})/{n})} \over n}\left\vert QB_{X}-\dot aB\right\vert ^{({1}/{n})-1}\big[(2-n)B_{X}(QB_{X}-\dot aB) -\cr &B(QB_{XX}-B\dot a_X)\big]-\Gamma{Q^m \over \vert B\vert ^{m+1}}.} $$This expression has two terms that are determined by bed properties – its geometry (the first term) and basal shear (the second term). The first term, which includes bed elevation B (note that B is a negative quantity), its slope B X, and curvature B XX, can be regarded as a form drag, which is the resistance to ice flow that arises due to the shape of the bed (e.g. Schoof, Reference Schoof2002a). The last term of (26) is the basal shear expressed in terms of ice flux 
$\displaystyle {Q=U_0H_0}$. The fact that the basal shear is an additive term suggests that in the low basal and driving stress regimes, ice flow is less sensitive to a specific form of basal shear (‘sliding law’) as long as its magnitude is low.
In general (26) is not compatible with (25d), suggesting the presence of a boundary layer at the divide at first order. However, as we have already noted, our interest is not in the details at the divide, and we therefore impose a further constraint that 
$(n-2)B_X\dot a+ B\dot a_X=0$ at X=0. This is satisfied if, e.g., the bed elevation and accumulation rate have no gradients at the divide, i.e., B X(X = 0) = 0 and 
$\dot a_X(X=0)=0$.
Taking into account the boundary condition (25c), a solution of (25b) is
$$ U_{1}=-U_0{H_1 \over H_0}. $$H 1 is determined by integrating (26) using the flotation condition (25f).
Except for the case B = const, which is of limited interest, (26) cannot be integrated analytically. H 1 is ice surface elevation to leading order. It is determined by the integral of expression (26) that itself depends on only the local bed geometry and basal shear. The surface elevation is a non-local function of these bed properties, and depends on their integral effects through the length of the ice stream.
At this point the zeroth- and first-order solutions are determined for any given location of the junction between an ice stream and ice shelf, but we have yet to determine the position of the grounding line X g. For this purpose we use the remaining stress boundary condition (19d), and determine the grounding line location through the simultaneous solution of this condition and the momentum and mass-balance solutions in the interior of ice stream.
Combining the zeroth- and first-order stress conditions (22e) and (25e), the flotation conditions, (22f) and (25f), the zeroth-order solutions (23) and the first-order solution for U 1 (27), the stress condition (19d) becomes
$$ \eqalign{ \vert QB_{X}&-\dot aB\vert ^{({1}/{n})-1}\left[\left(1-{\delta \over n}\right)\left(QB_{X}-\dot aB\right)-{\delta \over n}Q\left(H_{1X}+B_{X}\right)\right] =\cr & {1+\delta \over 2}\left\vert B\right\vert ^{(({n+2})/{n})}, \quad X=X_{\rm g}.} $$where H 1X is given by (26) evaluated at X = X g. Values of X g are determined as the roots of (28). It is possible that the choices of the functions B and 
$\dot a$ can produce zero, one or more steady-state locations of the grounding line.
Comparing (28) and (24) one can see that acknowledging the role of the basal shear stress and the form drag further increases the parameters on which the flux at the ground line depends to include the basal curvature and the basal sliding coefficient Γ. The dimensional form of (28) is
$$ \eqalign {&\left\vert qb_{x}-\dot ab\right\vert ^{(1-n)/{n}} \left[ \left( 1-{\delta \over n} \right) \left(qb_{x}-\dot ab\right) - {\delta \over n} q \left( h_{1x}+b_{x}\right) \right] = \cr & \quad \rho g^\prime A^{({1}/{n})} \vert b\vert ^{({n+2})/{n}} {1+\delta \over 4}, \quad x=x_{\rm g},} $$where
$$\eqalign{ h_{1x}&={1 \over \rho g^\prime \vert b\vert }\left\{{2A^{-({1}/{n})} \over n\vert b\vert ^{2 \over n}}\left\vert qb_{x}-\dot ab\right\vert ^{({1}/{n})-1}\big[(2-n)b_{x}(qb_{x}-\dot ab)-\right. \cr & \left. \quad b(qb_{xx}-b\dot a_x)\big]-C\left({q \over \vert b\vert }\right)^m\right\},} $$where b(x) <0.
Expression (29) may be compared with the relationship derived for unconfined configurations in a regime with a large magnitude basal shear (Schoof, Reference Schoof2007b)
$$ q=\left({A\left(\rho g\right)^{n+1}\delta^n \over 4^nC}\right)^{{1}/({m+1})}\left({\vert b\vert \over 1-\delta}\right)^{({m+n+3})/{m+1}}, \quad x=x_{\rm g}, $$where we have applied the flotation condition on the RHS of (31) so to make it comparable with (29). Equation (29) is considerably more complex than (31), and in particular it depends upon the bed shape, and the accumulation rate and its gradient at the grounding line, in addition to the bed elevation itself. According to these, ice streams with the same bed depth at the grounding line can have different fluxes and vice versa; and the relationship can be a multivalued function. Expression (28) results in an ice flux with a finite magnitude even when the basal shear is zero, in contrast to the infinite flux which is the consequence of (31) as C →0.
Comparison of approximate and ‘exact’ solutions
To assess the behaviors of the zeroth- and first-order solutions, and to compare these with the ‘exact’ solutions, we compute the grounding line position, and the relationship between flux and thickness at the grounding line, for the sinusoidal bed illustrated in Figure 4a, as the accumulation rate is varied. The grounding line positions for the case of the sliding parameter C=0 are shown in Figure 4a, the relationships between ice flux and thickness for various values of the sliding parameter C in Figures 4b and 4c. As Figure 4b shows, the relationship between flux and thickness is not a simple, monotonic curve. In this case, for a given flux, steady states can be found at two thicknesses, or for a given thickness, steady states can be found for two values of flux. As C increases, the relative importance of form drag and basal shear changes, and the difference between multiple values decreases. Figures 4a and 4b show solutions obtained from the zeroth- and first-order solutions. The multivalued behavior is qualitatively captured in the zeroth-order solution (black asterisks), although the addition of the role of form drag, which is accounted for when the first-order solution is included and C = 0, and the basal shear stress as C increases are quantitatively different. The quantitative difference is apparent too in the grounding line locations. Figure 4a shows three examples of the grounding line positions computed from the zeroth-order expression (24) (black asterisks labeled ‘1–3’). These positions are substantially different from the ‘exact’ positions and approximate ones computed using the first-order solutions (filled blue circles and diamonds labeled ‘1–3’).
Fig. 4. Conditions at the grounding line. (a) The grounding line position for sliding parameter C=0 and discreet values of the accumulation rate 
$\dot a$; (b) and (c) the relationship between ice flux and ice thickness computed for various values of the accumulation rate 
$\dot a$ and sliding parameter C. The values of 
$\dot a$ are discretely chosen, ranging from 0.18 to 39.8 m a−1. The size of symbols in panels (a) and (b) correspond to the values of 
$\dot a$. In panels (a) and (b) circles are ‘exact’ computations, diamonds are the roots of expression (29), asterisks are the roots of the zero-order expression (24). In panel (a) black asterisks with black labels ‘1–3’ are three examples of the roots of the zeroth-order expression (24); filled symbols with blue labels ‘1–3’ are the roots of expression (29) with the same values of 
$\dot a$ as those denoted by black asterisks. The labels ‘1–3’ in panel (b) correspond to the labeled grounding line locations in panel (a). Circles in panel (c) are ‘exact’ solutions and solid lines are the boundary layer relationship (31) derived by Schoof (Reference Schoof2007b). Notice logarithmic scale on the vertical axis on panel (c). The bed elevation b(x) has the same analytic form and stiffness parameter A has the same value as those in Figure 3. In all experiments the ice flow is from left to right.
Figures 4a and 4b also include a comparison with the ‘exact’ solutions (denoted by circles). These show a good agreement between the ‘exact’ results and expression (29). This close coincidence indicates that the scale assignment of δ provides a description that closely resembles the actual dynamics. That accurate solutions require a solution to O(δ) reflects the fact that the value of δ is 0.1. This numerical value is material, in that accurate solutions then need to account for the form drag, and the basal shear if it is present. In detail, there are values of accumulation rate for which the zeroth- and first-order solutions do not provide a solution for a steady-state grounding line, although the ‘exact’ solutions do provide one. These are indicated by fewer diamond and asterisk symbols indicating the approximate and zero-order solutions compared to the circle symbols indicating the ‘exact’ solutions in Figures 4a and 4b.
Figure 4c extends the range of 4b to allow a comparison for values of C of which the basal shear stress is much larger than the contribution of form drag. It shows that the effects of form drag are completely suppressed and the flux–thickness relationship follows expression (31) (purple line in Fig. 4c). On the other hand, for smaller values of C expression (31) (solid lines in Fig. 4c) consistently overestimates ice flux (and, as we have noted, results in infinity if the basal shear is zero). (Because (29) is not applicable for high values of C, in Figure 4c we have used ‘exact’ solutions for consistency.)
In circumstances where the basal shear is small, the surface slopes and driving stress are determined by the gradient of the longitudinal stress τx (the first term on the RHS of (25a) and (30)), which in turn is determined by the ice stiffness parameter A and spatial variability of the bed elevation. Gudmundsson (Reference Gudmundsson2003) and Sergienko (Reference Sergienko2012) show that the effects of the variability in bed elevation on ice flow increase as the basal shear decreases. Figure 5 illustrates the effects of form drag on ice-stream flow. It shows the analytical (blue curves in panels 5a–5b and 5d–5e) and ‘exact’ solutions (red curves in panels 5a–5b and 5d–5e) with the sliding parameter C=0 for two bed elevations that have the same slope, but one of which has low-amplitude long-period (more than 10 ice thicknesses) undulations 
$\displaystyle {b(x)=b_0+ b_xx+b_a \sin ({\pi x}/{L})}$ where b 0 = −100 m, b x = −0.0005, b a = 5 m, and L=5 km (Figs 5a–5c) and whereas the second one does not 
$\displaystyle {b(x)=b_0+ b_xx}$, where b 0 = −200 m, b x = −0.0005 (Figs 5d–5f). The configurations of the two ice streams are significantly different: the grounding line position in the case of the bed with undulations is ~116 km upstream of the grounding line position in the case of the bed without undulations; the maximum ice flow speed is 26% slower in the former case than in the latter case; the maximum magnitude of the divergence of the longitudinal stress (and consequently driving stress) is more than two orders of magnitude larger in the former case than in the latter, although it is noticeable that the surface elevation in both cases is very similar. When the basal shear stress is low, the detailed shape of the bed has a strong influence on the location of the grounding line.
Fig. 5. Effects of form drag. Steady-state configurations computed with C=0 for two bed elevation profiles: b(x) = b 0 + b x x + b asin (πx/L) panels (a)–(c) and b(x) = b 0 + b x x panels (d)–(f). (a) and (d) profiles of surface and bed elevation, arrows show the ice-flow direction; (b) and (e) ice velocity; (c) and (f) the momentum-balance components. The terms of the momentum-balance (panels (c) and (f)) computed numerically. Note that panels (c) and (f) show − τd.
Linear stability analysis of steady states
The conclusion that a grounding line is unconditionally unstable if the bed slope is retrograde, due to Weertman (Reference Weertman1974) and given a more precise form by Schoof (2012), depends on the result (31) that the grounding line flux is a monotonically increasing function of the depth of the bed. However, as our results show, this result holds only when the basal shear and driving stress have large magnitudes. When the driving and basal shear stresses are small, the flux–thickness relationship (29) explicitly depends on the basal slope b x and curvature b xx, the accumulation rate 
$\dot a$ and its gradient 
$\dot a_x$, the basal shear stress, and the exponent n of the ice rheology (n >2 or n ≤2). Their various combinations could equally result in an increase or decrease in the ice flux with increase in ice thickness. One might anticipate that the stability of the grounding line position may exhibit a similar complexity.
To investigate the stability of the steady-state grounding line positions, we consider a time-variant problem with small deviations from a steady-state configuration, i.e.,
$$ H=H_{{\rm ss}}(X)+\sigma\tilde H(X,T), $$
$$ U=U_{{\rm ss}}(X)+\sigma\tilde U(X,T), $$
$$ X_{\rm g}=X_{{\rm gss}}+\sigma\tilde X_{\rm g}(T) $$where σ is a small perturbation parameter and H ss, U ss and X gss are the solutions of (19a), (21) and (19c)–(19e). Substituting these expressions into the momentum, mass-balance equations and the boundary conditions and collecting the terms of O(σ) yields the following perturbation problem
$$ \eqalign{\delta\left[\tilde H\left\vert U_{{\rm ss}X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}X}+{1 \over n}H_{{\rm ss}}\left\vert U_{{\rm ss}X}\right\vert ^{({1}/{n})-1}\tilde U_X\right]_X- \cr \tilde H_XH_{{\rm ss}}-\tilde H\left(H_{{\rm ss}}+B\right)_X-m\delta\Gamma \left\vert U_{{\rm ss}}\right\vert ^{m-1}\tilde U=0 }$$
$$ \tilde H_T+\left(\tilde UH_{{\rm ss}}+\tilde HU_{{\rm ss}}\right)_X\;=0 $$
$$ \tilde H_X=\tilde U=0,\quad X=0 $$
$$ {1 \over n} H_{{\rm ss}}\left\vert U_{{\rm ss}X}\right\vert ^{({1}/{n})-1}\tilde U_X+\left (H_{{\rm ss}}\left\vert U_{{\rm ss}X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}X}\right)_X\ \tilde X_{\rm g}=$$
$$ H_{{\rm ss}}\left(H_{{\rm ss}X}\tilde X_{\rm g}+{1 \over 2}\tilde H\right),\quad X=X_{{\rm gss}} $$
$$ \tilde H+H_{{\rm ss}X}\tilde X_{\rm g}=-{B_X \over 1-\delta}\tilde X_{\rm g},\quad X=X_{{\rm gss}}. $$ As in the previous section, we use the method of perturbation series expansion to construct to O(δ) solutions of (33). This is algebraically quite complex, and we give the details in Appendix B. However, the form of the solution for 
$\tilde X_{\rm g}$ is quite simple. It is,
$$ \tilde X_{{\rm g}}\simeq-{\tilde H(0) \over H_{{\rm ss}1X}+B_X}e^{\Lambda T} $$where 
$\tilde H(0)$ is a small initial perturbation, and the denominator is evaluated at X gss. As is apparent from this form, the stability of the steady state is entirely determined by the sign of a parameter Λ in (34). For Λ <0, they are stable and for Λ >0 they are unstable. This parameter is given by
$$ \eqalign{ \Lambda &=\left(B_XX_{{\rm gss}}-B\right)^{-1}\left\{\dot a-\left({2 \over n}-1\right){Q_{{\rm gss}} \over B}B_X \ +\right. \cr & \left. \left(3n-{2 \over n}+3\right) {B_X \left(Q_{{\rm gss}}B_X-B\dot a\right) \over B\left(H_{{\rm ss}1X}+B_X\right)} + \right. \cr & \left({1 \over n}-1\right) Q_{{\rm gss}} {Q_{{\rm gss}}B_{XX}-B\dot a_X \over Q_{{\rm gss}}B_{X}-B\dot a} + {1 \over H_{{\rm ss}1X} +B_{X}} \times\cr & \left. \left[Q_{{\rm gss}}\left({1 \over n}B_{XX}+H_{{\rm ss}1XX}\right)-\left({1 \over n}-1\right)B\dot a_X\right]\right\}} $$and is a function of the steady-state parameters taken at X = X gss. In dimensional terms Λ is
$$ \eqalign{ \Lambda&=\big(b_{x}x_{\rm g}-b\big)^{-1}\Bigg\{\dot {a}-\Bigg({2 \over n}-1\Bigg) {q \over b}b_{x}+\Bigg(3n-{2 \over n}+3\Bigg)\times \cr &{b_{x}\big(qb_{x}-b\dot {a}\big) \over b\big(\delta h_{1x}+ b_{x}\big)}+ \Bigg({1 \over n}-1\Bigg)q{qb_{xx}-b\dot {a}_x \over qb_{x}-b\dot {a}}+\cr & {1 \over n(\delta h_{1x}+ b_{x})}\big[q\big( \delta n h_{1xx}+b_{xx}\big)+\big(n-1\big)b\dot {a}_x\big]\Bigg\}} $$where h 1x is defined by (30).
As we surmised, the stability condition in the low stress regime depends in a complex way on the accumulation and basal geometry. In the low-driving and basal stress regime, stability is determined by a combination of the bed geometry – its slope b x and curvature b xx – the accumulation rate 
$\dot a$ and its gradient 
$\dot a_x$; it also implicitly depends on the basal shear τb and its gradient through h 1x and h 1xx terms. Figure 6a illustrates the stability condition Λ <0 for the case of bed elevation shown in Figure 4a and the sliding parameter C = 1.8 · 104 Pa m−1/3 s1/3. The grounding line can be both stable (open circles) and unstable (crossed circles) at locations with positive and negative bed slopes.
Fig. 6. Stability condition. (a) Bed elevation and steady-state grounding line positions for C = 1.8 · 104 Pa m−1/3 s1/3 and different discrete values of the accumulation rate 
$\dot a$; stable positions are open circles, unstable positions are crossed circles. (b) Time evolution of the grounding line marked by a red triangle in panel (a) perturbed from its steady state; (c) same as (b) for a position marked by blue triangle. In panels (b) and (c) red lines are numerical solutions; blue lines are fitting curves ~ e Λt, where Λ is computed using expression (36); dashed lines indicate steady-state positions marked by triangles in panel (a).
To validate the stability condition, which is derived based on approximations of the perturbation theory, we obtain the ‘exact’ solutions by solving the time-dependent problem (1)–(10) numerically for an ice stream whose grounding line is perturbed from its steady state. These simulations are performed for the configurations with the steady-state grounding line positions marked by triangles in Figure 6a. In both simulations, at t = 0 the grounding line position is displaced 1 km upstream of its steady-state position and the steady-state ice thickness and velocity are used as initial conditions. Figure 6b illustrates the grounding line behavior perturbed from an unstable position. After initial advance toward the steady-state position, the grounding line retreats away from it with the distance increasing exponentially in time (Λ >0). In contrast, the grounding line perturbed from the stable position (Fig. 6c) advances beyond its steady-state position and then retreats to it with the distance decreasing exponentially in time (Λ <0) and stabilizes at its steady-state location. In performing these ‘exact’ time-dependent simulations, we have found consistently that they contain the initial transients illustrated in Figures 6b–6c that are not accounted for in the linearized theory that leads to (34). However, we have also found, consistently, that the sign of Λ determined from (36) correctly determines the stability, and that, following the transient, the value of Λ determined from (36) correctly describes the time-dependent solutions. Due to the presence of the transients in the ‘exact’ solutions, the exponential functions illustrated in Figures 6b–6c are fit to the form 
$\tilde X_{{\rm g}0} e^{\Lambda (t-t_0)}$ (where X g0 and t 0 are fitting constants), with Λ determined from (36).
This behavior is quite distinct from that associated with high basal and driving stresses. The expression (36) can be compared with the stability criterion for the large basal and driving stress regimes, which is 
$q^\prime (h(x_{\rm g}\dot a(x_{\rm g})))\gt \dot a$ (Schoof, Reference Schoof2012), which cannot be satisfied on retrograde beds. In contrast, the stability of the two configurations illustrated in Figures 6b–6c is the reverse. The grounding line position on the prograde bed (6b) is unstable; that on the retrograde bed (6c) is stable.
Equation (34) shows that the reciprocal of the parameter Λ, 
$\displaystyle {t_e=\left \vert \Lambda \right \vert ^{-1}}$ is the e-folding time of adjustment of a marine ice sheet perturbed from its steady-state configuration. In contrast to the case of high basal stress, in which only the sign of Λ can be inferred (Schoof, Reference Schoof2012), (36) provides an expression for the magnitude of this constant, and it is of some interest to consider the magnitude of these. Figure 7 shows t e as the sliding parameter C and the accumulation rate 
$\dot a$ are varied. The values of t e vary from months to hundreds of years. These values provide, in any given case, some insight into the extent to which the stability analysis has value in considering the evolution of the sheet independently of variations in external forcing. At the lower end, which is short in comparison with variations of climate (parameterized here through the accumulation rate 
$\dot a$ and implicitly through the temperature-dependent ice-stiffness parameter A) or sea level, one may expect the sheet to exhibit stable or unstable behavior as identified by the stability analysis; at the upper end, an observed variation in the position of the grounding line may reflect the variations in climate or sea level and have little to do with any internal instability associated with its configuration.
Fig. 7. e-folding time (years) of all grounding line positions shown in Figure 4b.
Discussion and conclusions
Motivated by the present-day Siple Coast ice streams characterized by low surface slopes and consequently driving stress, we focused this study on the grounding line dynamics of unconfined marine ice sheets with grounded ice streams in a regime of low driving and basal stresses. Using the perturbation series expansion method, we have constructed analytic expressions for various approximation orders of the ice-stream thickness, velocity and flux. Our analysis suggests that the dynamics of such marine ice sheets is markedly different from those with ice streams experiencing high driving and basal stresses, configurations extensively studied previously (Weertman, Reference Weertman1974; Schoof, Reference Schoof2007b, Reference Schoof2012; Tsai and others, Reference Tsai, Stewart and Thompson2015). Our results demonstrate that in circumstances where the basal shear is low or absent through the length of the ice stream, there is no boundary layer in the vicinity of the grounding line as in the regime of high driving and basal stresses. The approximate expressions for ice thickness show that low surface slopes and the driving stress are a consequence of low or absent basal shear (at the zeroth-order the former are zero), and they increase in response to increase of the magnitude of the basal shear. Therefore, observed low surface slopes on the present-day ice stream are indicative of low basal shear.
In contrast to the high driving-stress regime in which the effects of the longitudinal stress divergence are small, in the low driving-stress regime, its magnitude can be of the same order as the other components of the ice-stream momentum balance. In circumstances where the basal shear is zero, the longitudinal stress divergence balances the driving stress at higher orders. Its magnitude increases with increasing spatial variability of the underlying bed geometry. The longitudinal stress divergence can be viewed as form drag that together with basal shear (when it is non-zero) inhibits ice flow (e.g. Schoof, Reference Schoof2002a; Gudmundsson, Reference Gudmundsson2003; Sergienko, Reference Sergienko2012). Low-amplitude bed undulations (on the order of a few meters) can result in form drag on the order of a few kPa (Figs 5a–c). This has the direct consequence that bed elevation measurements, which are presently accurate to 30–100 m, (Fretwell and others, Reference Fretwell2013), are insufficient to accurately determine the form drag, and also suggests that the accuracy of the surface observations used to constrain bed elevation based on mass-conservation principles (e.g. Morlighem and others, Reference Morlighem2017) may not be sufficient to resolve small amplitude bed topography.
The importance of the shape of the bed and its slope manifests itself in the functional dependence of ice flux on ice thickness at the grounding line (Eqns (24), (28) and (29)): it is no longer a single-valued monotonic function as in the high driving stress regime (31), and depending on the bed geometry and the accumulation rate, ice streams can have the same flux through the grounding line with significantly different ice thicknesses (Fig. 4). Our results show that in the low driving stress regime, the ice flux at the grounding line is consistently lower than predicted by expression (31), and is finite in the case of the zero basal shear, in contrast to the infinite flux according to (31). These results suggest that ice-sheet models that employ (31) relationship to parameterize the grounding line dynamics (e.g., Pollard and DeConto, Reference Pollard and DeConto2009) consistently overestimate the ice flux through the grounding line for ice streams with low driving and basal stress. Such models are also prone to unstoppable grounding line migration, which is an explicit consequence of the use of (31).
Linear stability analysis of steady-state configurations (section ‘Linear stability analysis of steady states’) suggests that in the low driving and basal stress regime, unconfined marine ice sheets do not conform exactly to the marine ice-sheet instability hypothesis, according to which marine ice sheets with grounding lines on beds with retrograde slopes (the bed elevation increases in the direction of ice flow) are always unconditionally unstable (Weertman, Reference Weertman1974; Schoof, Reference Schoof2007b,, Reference Schoof2012). Instead, in the regime studied here, the grounding line stability is determined by a combination of form drag (i.e., the geometric bed properties – elevation, slopes and curvature), basal shear and the accumulation rate and their gradients (Eqns (35) and(36)). These results demonstrate that the bed slope alone is not a sufficient indicator of the marine ice-sheet stability (Fig. 6).
The results of our linear-stability analysis have direct application to present-day ice sheets, despite being derived for highly simplified and idealized configurations. In the low basal stress regime, we are able to calculate the magnitude, in addition to the sign, of the e-folding time of the ice sheet in response to a small perturbation. Our analysis shows that in the low driving and basal stress regime the e-folding time can range from months to hundreds of years, depending on the bed properties (its geometry and basal shear), accumulation, its gradients and its higher spatial derivatives (Fig. 7). At the low end of this range, sub-monthly timescales are significantly shorter than those of climate or sea level variability. Conversely, with e-folding times of some hundreds of years or longer, an observed grounding line variation may result entirely from externally-forced changes (e.g., Fyke and others, Reference Fyke, Sergienko, Löfveström, Price and Lenaerts2018). This leads us to the view that care is needed, in considering the cause of grounding line retreat such as that presently observed in West Antarctica, to distinguish the effects of external forcing from those of a configurational instability.
Acknowledgments
We would like to thank Scientific Editor Jonathan Kingslake and two anonymous referees for thoughtful comments and useful suggestions that greatly improved the manuscript. OS is supported by awards NA08OAR4320752 and NA13OAR439 from the National Oceanic and Atmospheric Administration, US Department of Commerce. The statements, findings, conclusions and recommendations are those of the authors and do not necessarily reflect the views of the National Oceanic and Atmospheric Administration, or the US Department of Commerce.
Appendix A. Steady-state problem
Here we describe the zeroth and first orders of the steady-state problem (19a), (21) and (19c)–(19e) and their solutions. First, we construct solutions of the zeroth-order problem and express them through the external parameters of the problem (e.g., the accumulation rate, bed elevation). Next, we use the zeroth-order solutions to construct the first-order solutions. Substituting the zeroth- and first-order solutions to expression (19d) we derive an implicit condition for the location of the grounding line X g. The same expression represents a functional relationship between the ice flux and thickness at the grounding line. All expressions derived here are non-dimensional. The dimensional forms are listed in section ‘Steady states’ of the main text.
The zeroth-order problem is
$$ H_0(H_0+B)_X=0 $$
$$ (U_0H_0)_X=\dot a $$
$$ U_0=0, \quad X=0 $$
$$ (H_0+B)_X=0, \quad X=0 $$
$$ H_0\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{0X}={1 \over 2}H_0^2,\quad X=X_{\rm g} $$
$$ H_0=-B, \quad X=X_{\rm g}. $$Its solutions, H 0 and U 0, are determined from (A1a) and (A1b)
$$ H_0(X)=-B(X) $$
$$ U_0(X)={\int_0^X\dot adx \over H_0(X)}={Q \over H_0}, $$where 
$Q=\int _0^X \dot a dx$. Using
$$ U_{0X}= {Q_X \over H_0}-Q{H_{0X} \over H_0^2}={\dot a \over H_0}-Q{H_{0X} \over H_0^2}. $$The stress condition at the grounding line (A1e) can be written as
$$ \dot aH_0-QH_{0X}=\left({1 \over 2}\right)^nH_{0}^{n+2},\ X=X_g $$This expression represents the zeroth-order relationship between the ice flux and ice thickness at the grounding line. It does not depend on Γ, and consequently the grounding line positions determined as roots of this expression are the same for any values of sliding. The first-order problem is
$$ H_{1X}={1 \over H_0}\left[\left(H_0\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{0X}\right)_X -\Gamma\left\vert U_0\right\vert ^{m-1}U_0\right] $$
$$ U_1=-U_0{H_1 \over H_0} $$
$$ U_1=0, \quad X=0 $$
$$ H_1\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{0X}+{1 \over n}H_0\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{1X}=H_1H_0, \quad X=X_{\rm g} $$
$$ H_1=-B, \quad X=X_{\rm g}. $$Expression for U 1X is determined by differentiating (A5b)
$$ U_{1X}=-U_{0X}{H_1 \over H_0}+U_0{H_{0X}H_1 \over H_0^2}-U_0{H_{1X} \over H_0} $$and accounting for the flotation condition, (22f) and (A5e), this expression at X = X g is
$$ U_{1X}(X_{\rm g})=-U_{0X}+U_0{H_{0X} \over H_0}-U_0{H_{1X} \over H_0}. $$We note that all components of the above expression are determined from the zeroth-order solutions U 0 and H 0 (A2). Thus,
$$ U_{0X}={Q_X \over H_0}-Q{H_{0X} \over H_0^2}={\dot a \over H_0}-Q{H_{0X} \over H_0^2} $$where 
$Q_X=\dot a$ was taken into account.
The expression for H 0|U 0X|(1/n)−1U 0X can be written as
$$ H_0\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{0X}=\left\vert \dot aH_0^{n-1}-H_0^{n-2}QH_{0X}\right\vert ^{({1}/{n})-1}\left(\dot aH_0^{n-1}-H_0^{n-2}QH_{0X}\right), $$and the expression for (H 0|U 0X|(1/n)−1U 0X)X is
$$ \eqalign{ &\left(H_0\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{0X}\right)_X={H_0^{n-3} \over n}\left\vert \dot aH_0^{n-1} -H_0^{n-2}QH_{0X}\right\vert ^{({1}/{n})-1}\times\ \cr & \qquad\qquad\qquad \left[(n-2)H_{0X}(\dot aH_0-QH_{0X})+H_0(H_0\dot a_X-QH_{0XX})\right].} $$Substituting these expressions into (A5a) and rearranging terms yield
$$ \eqalign{ H_{1X}= & {H_0^{-(({n+2})/{n})} \over n}\left\vert \dot aH_0-QH_{0X}\right\vert ^{({1}/{n})-1}\bigg[(n-2)H_{0X}(\dot aH_0-QH_{0X})\,+ \cr & H_0(H_0\dot a_X-QH_{0XX})\bigg]-\Gamma{Q^m \over H_0^{m+1}}.} $$To O(δ), the boundary condition (19d) is
$$ H_0\left\vert U_{0X}\right\vert ^{({1}/{n})-1}U_{0X}\left(1+{\delta \over n}{U_{1X} \over U_{0X}}\right)= {H_0^2 \over 2}(1+\delta). $$Substituting (A9) and (A7) into (A12) and rearranging terms leads to
$$ \eqalign{&\left\vert \dot aH_0-QH_{0X}\right\vert ^{({1}/{n})-1}\left[\left(1-{\delta \over n}\right)\left(\dot aH_0-QH_{0X}\right)-{\delta \over n}Q\left(H_{1X}-H_{0X}\right)\right]\ =\cr&\quad {1+\delta \over 2}H_0^{(({n+2})/{n})}, \quad X=X_{\rm g}}, $$where H 1X is defined by (A11) and H 0 is defined by (A2a). This equation implicitly defines the location of the grounding line X g with Q expressed as 
$\displaystyle {Q(X_{\rm g})=\int _0^{X_{\rm g}}\dot adx}$.
The zeroth-order terms H 0 and U 0 are determined by (A2). The first-order term H 1(X) is recovered by integrating expression (A11) with the boundary condition (A5e) at X = X g, and the first-order term U 1(X) is determined by expression (A5b).
In order to asses the validity and accuracy of the approximate analytical expressions constructed above, we compare them with the ‘exact’, numerical solutions that are obtained using a finite-element solver COMSOL™(COMSOL, 2018). In all simulations, the grid resolution is spatially variable: it is 200 m through 95% of the length of the domain, and 1 m in the 5% closest to the grounding line position x g.
For a given bed elevation profile, b(x), accumulation rate 
$\dot a$, (assumed to be constant), the sliding law parameters C and m (chosen to be m = 1/n), and ice stiffness parameter A, the grounding line position is determined from expression (29) (we use the Matlab™ routine fzero to find roots of (29)) and the h k and u k, (k = 0, 1) are determined as described above. Unless otherwise noted, we use the bed elevation profile 
$\displaystyle {b(x)=b_0+b_a \cos ({\pi x}/{L})}$, with b 0 = −500 m, b a = 250 and L = 500 km, and the ice stiffness parameter A = 1.35 · 10−25 Pa−3 s−1 (which corresponds to T ice ≈ −20°C).
Figure 8 shows numerical and analytical ice thickness and velocity profiles computed for various values of the sliding parameter C and the accumulation rate 
$\dot a=1.89$ m a−1. The difference in the grounding line position computed with two methods is 738 m for C=0 and increases to 3.4 km for C = 7.6 · 104 Pa m−1/3 s1/3 (corresponds to a value of the non-dimensional parameter Γ = 10, Eqn (19a)) that are 0.16 and 0.47% respectively. The maximum difference in the ice thickness (Fig. 8c) and ice velocity (Fig. 8d) increases from 0.92 m and 2.76 m a−1 for C=0 to 10 m and 36.4 m a−1 for C = 7.6 · 104 Pa m−1/3 s1/3. This comparison suggests that the analytic expressions adequately approximate solutions of the problem (1), (5), (6)–(9) with sufficient accuracy.
Fig. 8. Steady-state profiles of (a) ice thickness (m) and (b) ice velocity (m a−1) computed for the same bed elevation and ice stiffness parameter as in Figure 3 and accumulation rate 
$\dot a=1.9\, {\rm m}\, {\rm a}^{-1}$ analytically (solid lines) and numerically (dashed lines) for various values of the sliding coefficient C (Pa m−1/3 s1/3). (c) and (d) differences between numerical and analytical results: (c) ice thickness Δh = h numerical − h analytical (m); (D) ice velocity Δu = u numerical − u analytical (m a−1).
Appendix B. Stability analysis
This section describes a stability analysis of steady-state solutions obtained in sections ‘Steady states’ and Appendix A ‘Steady-state problem’. Small deviations from a steady-state configuration can be written as
$$ H=H_{{\rm ss}}(X)+\sigma\tilde H(X,T),\quad U=U_{{\rm ss}}(X)+\sigma\tilde U(X,T),\quad X_{\rm g}=X_{{\rm gss}}+\sigma\tilde X_{\rm g}(T) $$where H ss, U ss and X gss are solutions of (19a), (21) and (19c)–(19e). Their substitution to the momentum and mass balances and boundary conditions (19) result in the following perturbation problem
$$ \eqalign{\delta\Big[\tilde H\big\vert U_{{\rm ss}X}\big\vert ^{({1}/{n})-1}U_{{\rm ss}X}+{1 \over n}H_{{\rm ss}}\big\vert U_{{\rm ss}X}\big\vert ^{({1}/{n})-1}\tilde U_X\Big]_X- \cr \tilde H_XH_{{\rm ss}}-\tilde H\left(H_{{\rm ss}}+B\right)_X-m\delta\Gamma \left\vert U_{{\rm ss}}\right\vert ^{m-1}\tilde U=0} $$
$$ \tilde H_T+\left(\tilde UH_{{\rm ss}}+\tilde HU_{{\rm ss}}\right)_X=0 $$
$$ \tilde H_X=\tilde U=0,\quad X=0 $$
$$\eqalign{& {1 \over n} H_{{\rm ss}} \left\vert U_{{\rm ss}X}\right\vert ^{({1}/{n})-1} \tilde U_X + \left( H_{{\rm ss}} \left\vert U_{{\rm ss}X} \right\vert ^{({1}/{n})-1} U_{{\rm ss}X} \right)_X \tilde X_{\rm g} = \cr & \quad H_{{\rm ss}} \left( H_{{\rm ss}X} \tilde X_{\rm g} + {1 \over 2} \tilde H \right), \quad X = X_{{\rm gss}}} $$
$$ \tilde H+H_{{\rm ss}X}\tilde X_{\rm g}=-{B_X \over 1-\delta}\tilde X_{\rm g},\quad X=X_{{\rm gss}} $$where (B2d) is obtained in the following way. To O(σ)
$$\eqalign{&H\left\vert U_X\right\vert ^{({1}/{n})-1} U_X \approx H_{{\rm ss}} \left\vert U_{{\rm ss}X} \right\vert ^{({1}/{n})-1} U_{{\rm ss}X} + \cr &\quad \sigma \left( \tilde H \left\vert U_{{\rm ss}X} \right\vert ^{({1}/{n})-1} U_{{\rm ss}X} + {1 \over n} \left\vert U_{{\rm ss}X} \right\vert ^{({1}/{n})-1} \tilde U_X \right),\quad X=X_{\rm g}} $$
$$ {H^2 \over 2}\approx{1 \over 2}(H^2_{{\rm ss}}+2\sigma\tilde HH_{{\rm ss}}),\quad X=X_{\rm g} $$Because 
${H_{{\rm ss}}\left \vert U_{{\rm ss}X}\right \vert ^{({1}/{n})-1}U_{{\rm ss}X}=({H_{{\rm ss}}^2}/{2})}$, 
${\tilde H\left \vert U_{{\rm ss}X}\right \vert ^{({1}/{n})-1}U_{{\rm ss}X} = {1 \over 2}\tilde HH_{{\rm ss}}}$.
To O(σ), 
$F(X_{\rm g})=F(X_{{\rm gss}})+\sigma \tilde X_{\rm g}F_X(X_{{\rm gss}})$, and the above expressions become
$$ \eqalign{ & H\left\vert U_X\right\vert ^{({1}/{n})-1}U_X\approx H_{{\rm ss}}\left\vert U_{{\rm ss}X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}X}+ \cr &\sigma\left(\tilde X_{\rm g}\left(H_{{\rm ss}}\left\vert U_{{\rm ss}X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}X}\right)_X+\tilde H\left\vert U_{{\rm ss}X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}X}+{1 \over n}H_{{\rm ss}}\left\vert U_{{\rm ss}X}\right\vert ^{({1}/{n})-1}\tilde U_X \right),\cr &\quad X=X_{{\rm gss}} } $$
$$ {H^2 \over 2}\approx{1 \over 2}\left[H^2_{{\rm ss}}+2\sigma\left(\tilde X_{\rm g}H_{{\rm ss}}H_{{\rm ss}X}+\tilde HH_{{\rm ss}}\right)\right],\quad X=X_{{\rm gss}} $$Using the steady-state condition at X gss, 
$\displaystyle {H_{{\rm ss}}\left \vert U_{{\rm ss}X}\right \vert ^{({1}/{n})-1}U_{{\rm ss}X}= {1 \over 2} H_{{\rm ss}}^2}$, the above expressions become (B2d).
As before, we write the steady-state solutions as
$$ H_{{\rm ss}}=H_{{\rm ss}0}+\delta H_{{\rm ss}1}+\cdots $$
$$ U_{{\rm ss}}=U_{{\rm ss}0}+\delta U_{{\rm ss}1}+\cdots $$Substituting these expressions into the flotation condition (B2e) leads to
$$ \tilde X_{\rm g}\simeq-{1 \over \delta}{\tilde H \over B_X+H_{{\rm ss}1X}}. $$Expression (B6) shows that 
$\tilde X_{\rm g}$ and 
$\tilde H$ cannot be of the same order in δ. For 
$\tilde X_{\rm g}$ to be O(1), 
$\tilde H$ must be O(δ) to leading order and because of the mass balance (B2b) so does 
$\tilde U$. Hence, the series expansion of the perturbed variables should take a form
$$ \tilde H=\delta\tilde H_{1}+\delta^2\tilde H_2\ldots $$
$$ \tilde U=\delta\tilde U_{1}+\delta^2\tilde U_2\ldots $$
$$ \tilde X_{\rm g}=\tilde X_{{\rm g}0}+\delta\tilde X_{{\rm g}1}+\cdots. $$Substitution of these and (B5) into the perturbation problem (B2) leads to the first-order problem
$$ -\tilde H_{1X}H_{{\rm ss}0}=0 $$
$$ \tilde H_{1T}+(\tilde H_1U_{{\rm ss}0}+\tilde U_1H_{{\rm ss}0})_X=0 $$
$$ \tilde H_{1X}=\tilde U_1=0,\quad X=0 $$
$$ \!\!\!\eqalign{ &{1 \over n}H_{{\rm ss}0}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}\tilde U_{1X}+\tilde X_{{\rm g}0}\left[H_{{\rm ss}1}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}0X}+ {1 \over n}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}1X} \right]_X\,=\cr & {1 \over 2}H_{{\rm ss}0}\tilde H_1+\tilde X_{{\rm g}0}\left(H_{{\rm ss}0}H_{{\rm ss}1X}+H_{{\rm ss}1}H_{{\rm ss}0X}\right),\quad X=X_{{\rm gss}}} $$
$$ \tilde X_{{\rm g}0}=-{\tilde H_1 \over H_{{\rm ss}1X}+B_X} $$The system (B8) is not as over-determined as it appears. Because H ss0 ≠ 0, the momentum balance (B8a) reduces to 
$\tilde H_{1X}=0$, and the first of the boundary conditions (B8c) is automatically satisfied. Solutions of (B8a) and (B8b) result in the following expressions for 
$\tilde H_1$ and 
$\tilde U_1$
$$ \tilde H_1=F(T) $$
$$ \tilde U_1= -{1 \over H_{{\rm ss}0}}\left[\int_0^Xdx\tilde H_{1T}+\tilde H_1U_{{\rm ss}0}\right] $$where F(T) is an arbitrary function of T determined from the boundary conditions and 
$\tilde H_{1T}=F^\prime (T)$.
Expression for 
$\tilde U_{1X}$ is determined from (B9b)
$$ \tilde U_{1X}= {H_{{\rm ss}0X} \over H_{{\rm ss}0}^2}\left[\int_0^Xdx\tilde H_{1T}+\tilde H_1U_{{\rm ss}0}\right]-{1 \over H_{{\rm ss}0}}\left[\tilde H_{1T}+\tilde H_1U_{{\rm ss}0X}\right] $$Using the zeroth-order steady-state condition (A1e),
$$ {1 \over n}H_{{\rm ss}0}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}={1 \over 2n}{H_{{\rm ss}0}^2 \over U_{{\rm ss}0X}} $$and substituting (B9a), (B10) and (B8e) to (B8d) leads to the following equation for F(T)
$$ \eqalign{ {H_{{\rm ss}0}^2 \over 2nU_{{\rm ss}0X}}&\left[{F^\prime \over H_{{\rm ss}0}^2}\left(H_{{\rm ss}0X}X_{{\rm gss}}-H_{{\rm ss}0}\right)+{F \over H_{{\rm ss}0}^2}\left(H_{{\rm ss}0X}U_{{\rm ss}0}-H_{{\rm ss}0}U_{{\rm ss}0X}\right)\right]\ =\cr & {F \over H_{{\rm ss}1X}+B_X}\bigg\{\left[H_{{\rm ss}1}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}0X}+{1 \over n}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}1X} \right]_X -\cr &H_{{\rm ss}0}H_{{\rm ss}1X}+H_{{\rm ss}1}H_{{\rm ss}0X}\bigg\}+{F \over 2}H_{{\rm ss}0}} $$After rearranging terms Eqn (B12) becomes
$$ \eqalign{ F^\prime&\bigg(H_{{\rm ss}0X}X_{{\rm gss}}-H_{{\rm ss}0}\bigg)=-F\bigg\{U_{{\rm ss}0}H_{{\rm ss}0X}-(n+1)H_{{\rm ss}0}U_{{\rm ss}0X} -\cr &{2nU_{{\rm ss}0X} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}}\bigg[\bigg(H_{{\rm ss}1}\bigg\vert U_{{\rm ss}0X}\bigg\vert ^{({1}/{n})-1}U_{{\rm ss}0X}+{1 \over n}H_{{\rm ss}0}\bigg\vert U_{{\rm ss}0X}\bigg\vert ^{({1}/{n})-1}U_{{\rm ss}1X}\bigg)_X- \cr &H_{{\rm ss}0}\bigg(H_{{\rm ss}0X}+H_{{\rm ss}1X}\bigg)\bigg]\bigg\}.} $$A solution of this equation is
$$ F(T)=F(0)e^{\Lambda T} $$where
$$ \eqalign{ \Lambda=&\bigg(H_{{\rm ss}0}-H_{{\rm ss}0X}X_{{\rm gss}}\bigg)^{-1}\bigg\{U_{{\rm ss}0}H_{{\rm ss}0X}-(n+1)H_{{\rm ss}0}U_{{\rm ss}0X}- \cr &{2nU_{{\rm ss}0X} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}}\bigg[\bigg(H_{{\rm ss}1}\bigg\vert U_{{\rm ss}0X}\bigg\vert ^{({1}/{n})-1}U_{{\rm ss}0X}+{1 \over n}H_{{\rm ss}0}\bigg\vert U_{{\rm ss}0X}\bigg\vert ^{({1}/{n})-1}U_{{\rm ss}1X}\bigg)_X -\cr & H_{{\rm ss}0}\bigg(H_{{\rm ss}0X}+H_{{\rm ss}1X}\bigg)\bigg]\bigg\}} $$The sign of the parameter Λ determines stability of the grounding line position: if perturbations from a steady-state configuration decay with time this configuration is stable, if they grow, it is unstable, i.e., for stable configurations Λ <0 and for unstable configurations Λ >0. Equation (B15) includes both the zeroth- and the first-order solutions of the steady-state problem. In order to express it through the parameters of the problem (e.g., bed elevation, accumulation rate, etc.) we substitute expressions for the steady-state terms obtained in Appendix A, and derive an expression for Λ using following intermediate steps. Using
$$ \eqalign{ \bigg(H_{{\rm ss}1}&\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}0X}+{1 \over n}H_{{\rm ss}0}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}1X}\bigg)_X\,= \cr &H_{{\rm ss}1X}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}0X}+{1 \over n}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}0XX}\,+ \cr & \quad{1 \over n}H_{{\rm ss}0X}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}1X}+{1 \over n}\left({1 \over n}-1\right)H_{{\rm ss}0}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-2}U_{{\rm ss}1X}U_{{\rm ss}0XX}\,+ \cr &\quad\quad{1 \over n}H_{{\rm ss}0}\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}1XX}} $$and rearranging terms
$$ \eqalign{ \Lambda=&\left(H_{{\rm ss}0}-H_{{\rm ss}0X}X_{{\rm gss}}\right)^{-1}\bigg\{U_{{\rm ss}0}H_{{\rm ss}0X}-H_{{\rm ss}0}U_{{\rm ss}0X}+{3nH_{{\rm ss}0}H_{{\rm ss}0X}U_{{\rm ss}0X} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}}\,- \cr & {H_{{\rm ss}0} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}}\bigg[H_{{\rm ss}0}\left(U_{{\rm ss}0XX}+U_{{\rm ss}1XX}\right)+H_{{\rm ss}0X}U_{{\rm ss}1X}\,+ \cr &2\left({1 \over n}-1\right)\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}0XX}U_{{\rm ss}1X}\bigg]\bigg\}.} $$With
$$ \eqalign{ U_{{\rm ss}1XX}=&-U_{{\rm ss}0XX}{H_{{\rm ss}1} \over H_{{\rm ss}0}}-2{U_{{\rm ss}0X} \over H_{{\rm ss}0}}\left(H_{{\rm ss}1X}-H_{{\rm ss}0X}{H_{{\rm ss}1} \over H_{{\rm ss}0}}\right)\ +\cr &2{U_{{\rm ss}0}H_{{\rm ss}0X} \over H_{{\rm ss}0}^2}\left(H_{{\rm ss}1X}-H_{{\rm ss}0X}{H_{{\rm ss}1} \over H_{{\rm ss}0}}\right) -{U_{{\rm ss}0} \over H_{{\rm ss}0}}\left(H_{{\rm ss}1XX}-H_{{\rm ss}0XX}{H_{{\rm ss}1} \over H_{{\rm ss}0}}\right)} $$
$$ \eqalign{ &U_{{\rm ss}1XX}(X_{{\rm gss}})=-U_{{\rm ss}0XX}-2{U_{{\rm ss}0X} \over H_{{\rm ss}0}}\left(H_{{\rm ss}1X}-H_{{\rm ss}0X}\right)\,+\cr &\quad2{U_{{\rm ss}0}H_{{\rm ss}0X} \over H_{{\rm ss}0}^2}\left(H_{{\rm ss}1X}-H_{{\rm ss}0X}\right) -{U_{{\rm ss}0} \over H_{{\rm ss}0}}\left(H_{{\rm ss}1XX}-H_{{\rm ss}0XX}\right)} $$the above expression becomes
$$ \eqalign{ \Lambda=&\left(H_{{\rm ss}0}-H_{{\rm ss}0X}X_{{\rm gss}}\right)^{-1}\bigg\{H_{{\rm ss}0}U_{{\rm ss}0X}-U_{{\rm ss}0}H_{{\rm ss}0X}+{3nH_{{\rm ss}0}H_{{\rm ss}0X}U_{{\rm ss}0X} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}} -\cr & {H_{{\rm ss}0} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}}\bigg[U_{{\rm ss}0}\left(H_{{\rm ss}0XX}-H_{{\rm ss}1XX}\right)+H_{{\rm ss}0X}U_{{\rm ss}1X}\, +\cr &2\left({1 \over n}-1\right)\left\vert U_{{\rm ss}0X}\right\vert ^{({1}/{n})-1}U_{{\rm ss}0XX}U_{{\rm ss}1X}\bigg]\bigg\}.} $$Using Eqn (A7) for U ss1X(X gss) and the steady-state boundary condition (A1e) again
$$ \eqalign{ \Lambda=&\left(H_{{\rm ss}0}-H_{{\rm ss}0X}X_{{\rm gss}}\right)^{-1}\bigg\{H_{{\rm ss}0}U_{{\rm ss}0X}+{(3n+1)H_{{\rm ss}0}H_{{\rm ss}0X}U_{{\rm ss}0X} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}} \,+\cr &\bigg({1 \over n}-1\bigg)U_{{\rm ss}0}H_{{\rm ss}0}{U_{{\rm ss}0XX} \over U_{{\rm ss}0X}}-{H_{{\rm ss}0} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}} \times \cr &\bigg[U_{{\rm ss}0}\bigg(H_{{\rm ss}0XX}-H_{{\rm ss}1XX}\bigg)-\bigg({1 \over n}-1\bigg)H_{ss0}U_{{\rm ss}0XX}\bigg]\bigg\}.} $$With
$$ \eqalign{ &U_{{\rm ss}0XX} = \left( {\dot a \over H_{{\rm ss}0}} - U_{{\rm ss}0} {H_{{\rm ss}0X} \over H_{{\rm ss}0}} \right)_X = \cr & \quad {1 \over H_{{\rm ss}0}} \left( \dot a_X - 2U_{{\rm ss}0X} H_{{\rm ss}0X} - U_{{\rm ss}0} H_{{\rm ss}0XX} \right)} $$Λ can be expressed as
$$ \eqalign{ \Lambda=&\left(H_{{\rm ss}0}-H_{{\rm ss}0X}X_{{\rm gss}}\right)^{-1}\bigg\{H_{{\rm ss}0}U_{{\rm ss}0X}-2\left({1 \over n}-1\right)U_{{\rm ss}0}H_{{\rm ss}0X}\,+\cr &{(3n-{2 \over n}+3)H_{{\rm ss}0}H_{{\rm ss}0X}U_{{\rm ss}0X} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}} +\left({1 \over n}-1\right){U_{{\rm ss}0} \over U_{{\rm ss}0X}}\left(\dot a_X-U_{{\rm ss}0}H_{{\rm ss}0XX}\right)\,-\cr & {H_{{\rm ss}0} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}} \left[U_{{\rm ss}0}\left({1 \over n}H_{{\rm ss}0XX}-H_{{\rm ss}1XX}\right)-\left({1 \over n}-1\right)\dot a_X\right]\bigg\}.} $$Using expression (A8) for U ss0X(X g) and rearranging terms
$$ \eqalign{ \Lambda=&\left(H_{{\rm ss}0}-H_{{\rm ss}0X}X_{{\rm gss}}\right)^{-1}\bigg\{\dot a-\left({2 \over n}-1\right)U_{{\rm ss}0}H_{{\rm ss}0X}\,+\cr &\left(3n-{2 \over n}+3\right){H_{{\rm ss}0X}\left(\dot a-U_{{\rm ss}0}H_{{\rm ss}0X}\right) \over H_{{\rm ss}1X}-H_{{\rm ss}0X}} \,+ \cr & \left({1 \over n}-1\right)H_{{\rm ss}0}U_{{\rm ss}0}{\dot a_X-U_{{\rm ss}0}H_{{\rm ss}0XX} \over \dot a-U_{{\rm ss}0}H_{{\rm ss}0X}}- {H_{{\rm ss}0} \over H_{{\rm ss}1X}-H_{{\rm ss}0X}} \, \times\cr & \left[U_{{\rm ss}0}\left({1 \over n}H_{ss0XX}-H_{{\rm ss}1XX}\right)-\left({1 \over n}-1\right)\dot a_X\right]\bigg\}} $$where
$$ \eqalign{ H_{{\rm ss}1XX}=&{H_{{\rm ss}0}^{-{n+2 \over n}} \over n}\bigg\vert \dot aH_{{\rm ss}0}-Q_{{\rm gss}}H_{{\rm ss}0X}\bigg\vert ^{{1 \over n}-1}\bigg[(n-2)\bigg(\dot aH_{{\rm ss}0}-Q_{{\rm gss}}H_{{\rm ss}0X}\bigg)\, \times\cr & \bigg(H_{{\rm ss}0XX}-{n+2 \over n}{H_{{\rm ss}0X} \over H_{{\rm ss}0}}\bigg) +(H_{{\rm ss}0}\dot a_X-Q_{{\rm gss}}H_{{\rm ss}0XX})\bigg\{\bigg(1+n-{2 \over n}\bigg)H_{{\rm ss}0X}\,+ \cr &\bigg({1 \over n}-1\bigg)H_{{\rm ss}0}{H_{{\rm ss}0}\dot a_X-Q_{{\rm gss}}H_{{\rm ss}0XX} \over \dot aH_{{\rm ss}0}-Q_{{\rm gss}}H_{{\rm ss}0X}}-{n+2 \over n}\bigg\} \,+\cr & H_{{\rm ss}0}\bigg(\dot a_XH_{{\rm ss}0X}+\dot a_{XX}H_{{\rm ss}0}-\dot aH_{{\rm ss}0XX}-Q_{{\rm gss}}H_{{\rm ss}0XXX}\bigg)\bigg] \,- \cr & m\Gamma\bigg\vert U_{{\rm ss}0}\bigg\vert ^{m-1}{U_{{\rm ss}0X} \over H_{{\rm ss}0}}+\Gamma\bigg\vert U_{{\rm ss}0}\bigg\vert ^{m-1}U_{{\rm ss}0}{H_{{\rm ss}0X} \over H_{{\rm ss}0}^2}.} $$In terms of 
$\displaystyle {Q_{{\rm gss}}=\int _0^{X_{{\rm gss}}}\dot adx}$ this expression is
$$ \eqalign{ \Lambda=&\bigg(H_{{\rm ss}0}-H_{{\rm ss}0X}X_{{\rm gss}}\bigg)^{-1}\bigg\{\dot a-\bigg({2 \over n}-1\bigg){Q_{{\rm gss}} \over H_{{\rm ss}0}}H_{{\rm ss}0X}\,+\cr &\bigg(3n-{2 \over n}+3\bigg){H_{{\rm ss}0X}\bigg(H_{{\rm ss}0}\dot a-Q_{{\rm gss}}H_{{\rm ss}0X}\bigg) \over H_{{\rm ss}0}\bigg(H_{{\rm ss}1X}-H_{{\rm ss}0X}\bigg)}\,+\quad\quad\quad\quad \cr &\bigg({1 \over n}-1\bigg)Q_{{\rm gss}}{H_{{\rm ss}0}\dot a_X-Q_{{\rm gss}}H_{{\rm ss}0XX} \over H_{{\rm ss}0}\dot a-Q_{{\rm gss}}H_{{\rm ss}0X}} \,- \cr &{1 \over H_{{\rm ss}1X}-H_{{\rm ss}0X}}\bigg[Q_{{\rm gss}}\bigg({1 \over n}H_{{\rm ss}0XX}-H_{{\rm ss}1XX}\bigg)-\bigg({1 \over n}-1\bigg)H_{{\rm ss}0}\dot a_X\bigg]\bigg\}\quad\quad\quad.} $$The dimensional form of Λ is expression (36).
Open access
