3.1 Bifurcation diagram

We consider the case of constant accumulation $a=\bar{a}$ , flat sea level, i.e. $S(x_{g})=0$ and a fixed shape of the bedrock $b(x)$ given by (Schoof Reference Schoof2007a ):

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle b(x)=-\!\left(729-2184.8\left({\displaystyle \frac{x}{L_{B}}}\right)^{2}+1031.72\left({\displaystyle \frac{x}{L_{B}}}\right)^{4}-151.72\left({\displaystyle \frac{x}{L_{B}}}\right)^{6}\right), & \displaystyle\end{eqnarray}$$

with $L_{B}=750\times 10^{3}~\text{m}$ . To obtain a starting solution for the continuation it is possible to use the boundary layer theory developed in Schoof (Reference Schoof2007b ). However, we found that an ice sheet surface profile of the form

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle s(x)=G_{1}\sqrt{1-\left({\displaystyle \frac{x}{x_{g}}}\right)^{2}}+G_{2}, & \displaystyle\end{eqnarray}$$

with $x\in [0,x_{g}]$ , is sufficient to give an initial state that will rapidly converge to a stationary (starting) solution using the Newton–Raphson method. The initial ice sheet thickness is then given by $h(x)=s(x)+b(x)$ and the initial velocity follows from assuming steady conditions and integrating the continuity equation:

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle u(x)=\frac{\bar{a}\,x}{h(x)}. & \displaystyle\end{eqnarray}$$

A typical height is chosen with $G_{1}=3\times 10^{3}$ m and a correction $G_{2}$ (in m) is used such that the flotation criterion (2.5) is satisfied: $G_{2}=(\unicode[STIX]{x1D70C}_{w}/\unicode[STIX]{x1D70C}_{i}-1)b(x_{g})$ , where we use that the sea level $S(x_{g})=0$ at the starting solution. Note that the profile (3.2) contains a steep gradient at the grounding line, which is essential for a converging Newton–Raphson iteration to a starting solution. The continuation in $A$ is initialised with a solution at $A_{0}=4.6416\times 10^{-24}$ . All other model parameters are given in table 2.

Figure 3. One-parameter bifurcation diagram (a) and solutions (b) with the parameter $A$ . The bedrock contains an upward slope which admits multiple steady states for a constant parameter (points b, d and f in panel (a)). Eigenvalue analysis shows two saddle-node bifurcations, $L_{1}$ and $L_{2}$ : a single eigenvalue crosses the imaginary axis to the positive right half plane at point c and returns to the left half plane at point e. The number of grid-points is $N=1600$ ; the other parameters are given in table 2.

From the steady starting solution at $A_{0}$ , a pseudo-arclength continuation traces the solution branch in the direction of decreasing $A$ , see figure 3. A decrease in the parameter $A$ in the SSA model corresponds to a decrease in temperature and an increase in ice growth. The bifurcation diagram presented reveals the multiple equilibria regime associated with hysteretic behaviour (Schoof Reference Schoof2007a ). For an interval of values of the parameter $A$ , three equilibria are distinguished for the case $A=8.5014\times 10^{-26}$ , marked with labels b, d and f in figure 3. At point c, an eigenvalue of (2.21) is found to cross the imaginary axis to the right half-plane, coinciding with the annihilation of a stable and an unstable branch in the direction of decreasing $A$ . Hence a saddle-node bifurcation is identified ( $L_{1}$ ). At point e, the same eigenvalue returns to the left half-plane through a second saddle-node bifurcation ( $L_{2}$ ). The values of the parameter $A$ at the labelled points, including the bifurcations, are shown in table 2.

Table 2. Parameter values for the experiment in figure 3, similar to the values chosen in Schoof (Reference Schoof2007a ) and Pattyn et al. (Reference Pattyn, Schoof, Perichon, Hindmarsh, Bueler, de Fleurian, Durand, Gagliardini, Gladstone and Goldberg2012).

Figure 4. Normalised components of the eigenmode that becomes unstable. The corresponding steady states are located at points b–f, described in figure 3 and table 2. We distinguish between a perturbation pattern related to sheet thickness ${\hat{h}}$ and a pattern related to ice velocity $\hat{u}$ . The signs of the eigenvectors are taken such that the perturbation in the grounding line is positive. From the eigenvectors and the equilibrium solution we compute a normalised spatial pattern of the evolution $\unicode[STIX]{x2202}\tilde{h}/\unicode[STIX]{x2202}t$ , together with normalised components $-\unicode[STIX]{x2202}(\hat{u} \bar{h})/\unicode[STIX]{x2202}x$ and $-\unicode[STIX]{x2202}(\bar{u}{\hat{h}})/\unicode[STIX]{x2202}x$ .

3.2 Instability mechanism

The advantage of the approach chosen here is that the spatial patterns of perturbations destabilising the marine ice sheet can be determined from the eigenvectors in (2.21). For the unstable equilibrium (point d) in figure 3, it is of interest to examine the eigenmode with a positive growth factor, showing in detail the characteristics of the instability. The eigenvector is made available using (2.21) and is depicted for the steady states (points b–f) in figure 4. The perturbations in thickness and velocity are taken corresponding to a positive grounding line perturbation. Note that a perturbation of the solution vector has the form $\hat{\boldsymbol{x}}\text{e}^{\unicode[STIX]{x1D70E}t}=[\hat{\boldsymbol{h}},\hat{\boldsymbol{u}},\hat{x}_{g}]^{\text{T}}\text{e}^{\unicode[STIX]{x1D70E}t}$ , with $\hat{x}_{g}$ the scalar grounding line perturbation. An eigenvector with corresponding eigenvalue $\unicode[STIX]{x1D70E}>0$ and $\hat{x}_{g}<0$ gives the destabilising pattern for unstable ice sheet retreat. By adjusting the sign of the eigenvector, such that $\hat{x}_{g}>0$ , we restrict our exposition to destabilising patterns for unstable ice sheet growth.

At the grounding line $x_{g}$ , the perturbation of the unstable steady state (point d) shows a slight decrease in ice thickness, while at (points b, c, e and f) a slight increase is observed. In the interior of the ice sheet a relatively large increase in ice thickness is visible for the unstable equilibrium (point d), indicating interior ice growth due to an imbalance between global accumulation and ice flux at the grounding line. That is, the reduced grounding line ice discharge is unable to accommodate the increased accumulation. The resulting interior thickness anomaly is central to the unstable growth/retreat mechanism, as we will see next.

The velocity perturbation at the grounding line shows an increase for stable states and a clear decrease for the unstable state (point d). Together with the negative perturbation in thickness this implies that, at (point d), there must be a decrease in flux $uh$ across the grounding line for a positive perturbation $\hat{x}_{g}>0$ . An increase in grounding line position implies a rise in global accumulation, hence the reduction in grounding line flux implies a net ice growth, confirming the marine ice sheet instability hypothesis. Note that a similar result holds if we take the perturbation in $x_{g}$ negative, giving a net ice loss and a retreat from equilibrium.

The continuation approach allows an efficient computation of flux perturbations using (2.21). From a linear stability analysis of (2.1) with perturbation $\tilde{h}=c{\hat{h}},~\tilde{u} =c\hat{u}$ around an equilibrium $\bar{h},\bar{u}$ we obtain an evolution equation for the thickness perturbation $\tilde{h}$ :

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\tilde{h}}{\unicode[STIX]{x2202}t}+c\frac{\unicode[STIX]{x2202}\hat{q}}{\unicode[STIX]{x2202}x}=0,\quad \text{with}\,c>0, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{q}(x)=\hat{u} \bar{h}+\bar{u}{\hat{h}}, & \displaystyle\end{eqnarray}$$

where we neglect higher-order terms. At the unstable steady state (point d) in figure 4, the thickness perturbation (green squares) shows positive growth, whereas the other points show a dampening. These patterns are determined by spatial derivatives of the perturbed advection of the steady thickness $\hat{u} \bar{h}$ (black dash-dotted line) and the advected thickness perturbation $\bar{u}{\hat{h}}$ (red dotted line). The latter clearly dominates the instability in (point d). Note that at the bifurcation points c and e the components of the perturbation flux have a compensating effect.

Figure 5. Perturbations in grounding line and accumulation flux against the parameter $A$ (a) and as a function of $x_{g}$ (b), together with the bedrock slope. The first and second saddle-node bifurcations are marked $L_{1}$ and $L_{2}$ . Dashed lines correspond to growing perturbations of unstable steady states. The perturbations $\hat{q}(\bar{x}_{g})$ and $a|\hat{x}_{g}|$ correspond to a positive grounding line perturbation $\hat{x}_{g}>0$ . At $L_{1}$ and $L_{2}$ , the grounding line bedrock slope is $-195~\text{mm}~\text{km}^{-1}$ and $-65~\text{mm}~\text{km}^{-1}$ respectively. The maximum slope in the unstable regime is $620~\text{mm}~\text{km}^{-1}$ . Note that grounding line flux perturbations depend on the steady grounding line position $\bar{x}_{g}$ , whereas accumulation flux perturbations depend on the grounding line perturbation $\hat{x}_{g}$ .

To investigate how a perturbation changes from stable to unstable through the saddle-node bifurcation $L_{1}$ , we compute the accumulation and grounding line fluxes. The steady state $(\bar{h},\bar{u})$ gives a balance:

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{q}(x)=\bar{u}\bar{h}=ax. & \displaystyle\end{eqnarray}$$

In figure 5 we show perturbations of the balance (3.6) at the grounding line. A perturbation in accumulation flux is given by $a|\hat{x}_{g}|$ , a grounding line flux perturbation by $\hat{q}(\bar{x}_{g})$ in (3.4). The perturbations are plotted against $A/A_{0}$ (figure 5 a) and $\bar{x}_{g}$ (figure 5 b), together with the bifurcation points and the bedrock slope. At the first saddle-node bifurcation $L_{1}$ , the flux $\hat{q}(x_{g})$ becomes smaller than the accumulation $a|\hat{x}_{g}|$ . Beyond this point, a change in accumulation due to $\hat{x}_{g}$ is not balanced by the grounding line flux and, hence, the perturbation $\tilde{\boldsymbol{x}}$ changes from damped to growing. At $L_{2}$ , $\hat{q}(x_{g})$ becomes greater than $a|\hat{x}_{g}|$ and the mode is damped again.

In figure 5(b) we also plot the bedrock slope, taken positive when the bed elevation increases with $x_{g}$ , that is

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle r_{bed}=-b^{\prime }(x_{g}), & \displaystyle\end{eqnarray}$$

with $b(x)$ as in (3.1). Note that the sign switch in $\hat{q}(x_{g})$ coincides with the sign switch in the bedrock slope. The grounding line flux will increase for a positive $\hat{x}_{g}$ between the bifurcation $L_{1}$ and the point of zero bedrock slope, but, since the change is less than the change in accumulation, the ice sheet will grow. Thus, figure 5 confirms that an eigenvalue of the ‘full model’ in Schoof (Reference Schoof2007a ) becomes positive when

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle a|\hat{x}_{g}|-\hat{q}(x_{g})>0. & \displaystyle\end{eqnarray}$$

Using a continuation of steady states with the original SSA equations, we find that the flux perturbations and their relative magnitude fully describe the instability mechanism, confirming the analysis in Schoof (Reference Schoof2012), but without the need for asymptotic expansions. In addition, the eigenvectors reveal destabilising interior patterns with, most interestingly, interior thickness anomalies and their advection. These turn out to play a major role in the growth and retreat of the ice sheet.

3.3 Gravitational sea level effect

Next, instead of a fixed $S(x)=0$ , we use the gravitational sea level effect (GSLE) formulation given by (2.8). With the added sea level we repeat the continuation in the parameter $A$ (§ 3.1). The new bifurcation diagram is shown in figure 6(a) (red triangles). The width of the multiple equilibria regime has decreased significantly and the bifurcations have shifted in the direction of decreasing $A$ . The same grounding line advance now requires a greater decrease in $A$ compared to the fixed sea level case. The parameter difference between the saddle-node bifurcations with fixed sea level is $|A_{L_{2}}-A_{L_{1}}|/A_{0}=3.34\times 10^{-2}$ ; for the variable sea level case we find $|A_{S_{2}}-A_{S_{1}}|/A_{0}=3.18\times 10^{-3}$ . Here we denote the parameter values at the bifurcations $L_{1},L_{2},S_{1}$ and $S_{2}$ with $A_{L_{1}},A_{L_{2}},A_{S_{1}}$ and $A_{S_{2}}$ respectively. The unstable regime has a horizontal extent of approximately 326 km in the fixed sea level case; for the variable sea level case we find 248 km.

Figure 6. (a) Bifurcation diagrams in the parameter $A$ . Fixed sea level bifurcation diagram (black) and bifurcation diagram with variable sea level (red). (b) Accumulation and grounding line flux perturbations as a function of $x_{g}$ for variable sea level, together with the bedrock slope and the resulting effective topographical slope. The first and second saddle-node bifurcations are marked $L_{1}$ ( $S_{1}$ ) and $L_{2}$ ( $S_{2}$ ) for fixed (variable) sea level. Dashed lines correspond to (perturbations of) unstable steady states.

Changes in the instability mechanism due to the presence of a varying sea level are investigated using the eigenvectors from the linear stability analysis (2.21). In figure 6(b), the perturbation in flux intersects the perturbation in accumulation at $S_{1}$ and the ice sheet state becomes unstable. At the second saddle-node bifurcation $S_{2}$ the reverse occurs. Both sign switches in bedrock slope are now located in the stable regime. This confirms that, for our choice of sea level response, a reverse bed slope is reduced from a sufficient to a necessary condition for the instability mechanism.

In figure 6(b) we also plot the effective topographic slope, taken positive when the bed elevation increases with $x_{g}$ , that is

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle r_{eff}=-(b^{\prime }(x_{g})+S^{\prime }(x_{g})). & \displaystyle\end{eqnarray}$$

Similar to the fixed sea level case, there is a direct correspondence between the grounding line flux perturbation and the effective topography. Note also that the magnitudes of the fluxes in the unstable regime are strongly affected by the magnitude of the slope. Hence, the dampening effect of the added variable sea level acts through the response of the unstable flux perturbations to the significantly decreased effective slope, in agreement with the analysis in Gomez et al. (Reference Gomez, Mitrovica, Huybers and Clark2010). Whether the dampening extends to the stochastic case is explored in the next section.

4.1 Additive noise in the accumulation rate

The linearised system of SDAEs (cf. (2.22)) that arises when adding noise to the accumulation $a$ is given by

(4.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D648}\,\text{d}\tilde{\boldsymbol{x}}=\unicode[STIX]{x1D645}\tilde{\boldsymbol{x}}~\text{d}t+\left[\begin{array}{@{}c@{}}B_{a}\\ 0\\ 0\\ \end{array}\right]\text{d}\boldsymbol{W},\quad \text{with}~\tilde{\boldsymbol{x}}=\left[\begin{array}{@{}c@{}}\tilde{\boldsymbol{h}}\\ \tilde{\boldsymbol{u}}\\ \tilde{x}_{g}\end{array}\right]\text{ and}~\unicode[STIX]{x1D648}=\left[\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D644} & 0 & M_{str}\\ 0 & 0 & 0\\ 0 & 0 & 0\end{array}\right], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D645}$ denotes the Jacobian matrix. As the accumulation is taken uniform over the sheet, the additive noise is integrated with respect to a one-dimensional Wiener process, that is, $N_{w}=1$ and $B_{a}\in \mathbb{R}^{N}$ with elements $(B_{a})_{i}=\unicode[STIX]{x1D702}_{ac}$ , where $\unicode[STIX]{x1D702}_{ac}$ ( $\text{my}^{-1}$ ) determines the noise amplitude and is taken at 10 % of the accumulation constant.

The matrix $\unicode[STIX]{x1D648}$ in the left-hand side of (4.1) is not invertible. In order to achieve a problem of the form (2.23), we need to bring $\unicode[STIX]{x1D648}$ in block-diagonal form and eliminate the system of algebraic equations. The transformation $\tilde{\boldsymbol{x}}=\unicode[STIX]{x1D64D}\boldsymbol{z}$ with

(4.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64D}=\left[\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D644} & 0 & -M_{str}\\ 0 & \unicode[STIX]{x1D644} & 0\\ 0 & 0 & \unicode[STIX]{x1D644}\end{array}\right], & & \displaystyle\end{eqnarray}$$

gives a system of SDAEs in a more convenient form:

(4.3) $$\begin{eqnarray}\displaystyle \left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D644} & 0\\ 0 & 0\end{array}\right]\left[\begin{array}{@{}c@{}}\text{d}\,\boldsymbol{z}_{1}\\ \text{d}\,\boldsymbol{z}_{2}\end{array}\right]=\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D645}_{11} & \unicode[STIX]{x1D645}_{12}\\ \unicode[STIX]{x1D645}_{21} & \unicode[STIX]{x1D645}_{22}\end{array}\right]\unicode[STIX]{x1D64D}\left[\begin{array}{@{}c@{}}\boldsymbol{z}_{1}\\ \boldsymbol{ z}_{2}\end{array}\right]\text{d}t+\left[\begin{array}{@{}c@{}}B_{a}\\ 0\end{array}\right]\text{d}\boldsymbol{W}, & & \displaystyle\end{eqnarray}$$

where the second row combines the algebraic constraints and provides an expression for $\boldsymbol{z}_{2}$ in terms of $\boldsymbol{z}_{1}$ . Substituting this expression in the first row eliminates $\boldsymbol{z}_{2}$ from (4.3) (Baars et al. Reference Baars, Viebahn, Mulder, Kuehn, Wubs and Dijkstra2017) and gives an OU-process in $\boldsymbol{z}_{1}$ :

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \text{d}\boldsymbol{z}_{1}=\unicode[STIX]{x1D64E}\boldsymbol{z}_{1}\,\text{d}t+B_{a}\,\text{d}\boldsymbol{W},\quad \text{with}\,\unicode[STIX]{x1D64E}=\tilde{\unicode[STIX]{x1D645}}_{11}-\tilde{\unicode[STIX]{x1D645}}_{12}\tilde{\unicode[STIX]{x1D645}}_{22}^{-1}\tilde{\unicode[STIX]{x1D645}}_{21},\tilde{\unicode[STIX]{x1D645}}=\unicode[STIX]{x1D645}\unicode[STIX]{x1D64D}. & \displaystyle\end{eqnarray}$$

The associated Lyapunov problem with $\unicode[STIX]{x1D655}_{11}=E[\boldsymbol{z}_{1}\boldsymbol{z}_{\mathbf{1}}^{\text{T}}]$ is given by

(4.5) $$\begin{eqnarray}\displaystyle S\unicode[STIX]{x1D655}_{11}+\unicode[STIX]{x1D655}_{11}S^{\text{T}}+B_{a}B_{a}^{\text{T}}=0. & & \displaystyle\end{eqnarray}$$

A solution to (4.5) is obtained using the RAILS solver (Baars et al. Reference Baars, Viebahn, Mulder, Kuehn, Wubs and Dijkstra2017). To obtain the blocks $\unicode[STIX]{x1D655}_{12},\unicode[STIX]{x1D655}_{21}$ and $\unicode[STIX]{x1D655}_{22}$ , the algebraic constraints in (4.3) are used. The full covariance matrix of the perturbation $\tilde{\boldsymbol{x}}$ is then given by

(4.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63E}=E[\tilde{\boldsymbol{x}}\tilde{\boldsymbol{x}}^{\text{T}}]=E[\unicode[STIX]{x1D64D}\boldsymbol{z}\,\boldsymbol{z}^{\text{T}}\unicode[STIX]{x1D64D}^{\text{T}}]=\unicode[STIX]{x1D64D}\unicode[STIX]{x1D655}\unicode[STIX]{x1D64D}^{\text{T}}. & & \displaystyle\end{eqnarray}$$

The first empirical orthogonal function (EOF), i.e. the dominant eigenvector corresponding to the rightmost eigenvalue of the covariance matrix, is shown in figure 7 (closed markers) for three stable equilibria with their grounding line at points $P_{1}=40.96~\text{km}$ , $P_{2}=16.384~\text{km}$ and $P_{3}=2048~\text{m}$ to the left of the grounding line position at $S_{1}$ . Noise is added to the accumulation with standard deviation $\unicode[STIX]{x1D702}_{ac}=\bar{a}/10~\text{my}^{-1}$ . The position of the equilibria and the explained variance of the EOFs are given in table 3. In the case of noise in accumulation, the EOFs correspond very well to the eigenvectors of the linear stability problem (dotted in figure 7).

Figure 7. First empirical orthogonal functions (EOF) for grounding lines at points $P_{1}$ (left), $P_{2}$ (centre) and $P_{3}$ (right) to the left of the first saddle-node bifurcation $S_{1}$ , together with the dominant eigenvector of the Jacobian matrix corresponding to the rightmost eigenvalue (dotted). Noise is added to the uniform accumulation with standard deviation $\unicode[STIX]{x1D702}_{ac}=\bar{a}/10~\text{my}^{-1}$ (closed markers) or to the sea level at the grounding line with $\unicode[STIX]{x1D702}_{sl}=1~\text{m}$ (open markers). The separate components corresponding to thickness (black squares) and velocity (red triangles) are normalised and their sign is taken according to a positive grounding line perturbation $\tilde{x}_{g}>0$ . The explained variances of the EOFs at $P_{1}$ , $P_{2}$ and $P_{3}$ are given in table 3.

Table 3. Horizontal position of the points $P_{1}$ , $P_{2}$ and $P_{3}$ . The grounding line at the bifurcation $S_{1}$ is positioned at 991.261 km from the origin. The explained variance (scaled rightmost eigenvalue) corresponding to each EOF, $\unicode[STIX]{x1D70E}_{expl}^{2}$ , is given for noise in accumulation and noise in sea level (see also figure 8).

4.2 Additive noise in sea level

Another source of variability is the sea level at the grounding line $S(x_{g})$ . In the model, the sea level enters through the flotation condition (2.5). Adding noise to the sea level gives a system of SDAEs, where the stochastic forcing appears in the algebraic constraints. Using the transformation $\tilde{\boldsymbol{x}}=R\boldsymbol{z}$ and again combining the algebraic constraints, we find:

(4.7) $$\begin{eqnarray}\displaystyle \left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D644} & 0\\ 0 & 0\end{array}\right]\left[\begin{array}{@{}c@{}}\text{d}\boldsymbol{z}_{1}\\ \text{d}\boldsymbol{z}_{2}\end{array}\right]=\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D645}_{11} & \unicode[STIX]{x1D645}_{12}\\ \unicode[STIX]{x1D645}_{21} & \unicode[STIX]{x1D645}_{22}\end{array}\right]\unicode[STIX]{x1D64D}\left[\begin{array}{@{}c@{}}\boldsymbol{z}_{1}\\ \boldsymbol{ z}_{2}\end{array}\right]\text{d}t+\left[\begin{array}{@{}c@{}}0\\ B_{s}\end{array}\right]\text{d}\boldsymbol{W}. & & \displaystyle\end{eqnarray}$$

The noise forcing now enters through the flotation condition and is determined via $B_{s}\in \mathbb{R}^{(N+1)}$ , with elements

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle (B_{s})_{i}=0,\quad i=1,\ldots ,N-1, & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle (B_{s})_{N}=-\frac{1}{2}\left(1-{\displaystyle \frac{\unicode[STIX]{x1D70C}_{i}}{\unicode[STIX]{x1D70C}_{w}}}\right)\unicode[STIX]{x1D70C}_{w}g\unicode[STIX]{x1D702}_{sl}, & \displaystyle\end{eqnarray}$$

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle (B_{s})_{N+1}=-2\frac{\unicode[STIX]{x1D70C}_{w}}{\unicode[STIX]{x1D70C}_{i}}\unicode[STIX]{x1D702}_{sl}. & \displaystyle\end{eqnarray}$$

The final two entries of $B_{s}$ follow directly from incorporating sea level noise in the two discretised right boundary conditions in appendix A, specifically (A 23) and (A 24), where the stochastic terms are separated to give the entries in $B_{s}$ . Eliminating $\boldsymbol{z}_{2}$ gives an OU-process in $\boldsymbol{z}_{1}$ :

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle \text{d}\boldsymbol{z}_{1}=\unicode[STIX]{x1D651}\boldsymbol{z}_{1}\,\text{d}t+\tilde{B}_{s}\,\text{d}\boldsymbol{W},\quad \text{with}~\unicode[STIX]{x1D651}=\tilde{\unicode[STIX]{x1D645}}_{11}-\tilde{\unicode[STIX]{x1D645}}_{12}\tilde{\unicode[STIX]{x1D645}}_{22}^{-1}\tilde{\unicode[STIX]{x1D645}}_{21},\tilde{B}_{s}=-\tilde{\unicode[STIX]{x1D645}}_{12}\tilde{\unicode[STIX]{x1D645}}_{22}^{-1}B_{s}. & \displaystyle\end{eqnarray}$$

The covariance matrix of the perturbation $\tilde{\boldsymbol{x}}$ is subsequently found by solving

(4.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D651}\unicode[STIX]{x1D655}_{11}+\unicode[STIX]{x1D655}_{11}V^{\text{T}}+\tilde{B}_{s}\tilde{B}_{s}^{\text{T}}=0, & \displaystyle\end{eqnarray}$$

computing the remaining blocks in $\unicode[STIX]{x1D655}$ using the algebraic constraints in (4.7) and computing $\unicode[STIX]{x1D63E}=\unicode[STIX]{x1D64D}\unicode[STIX]{x1D655}\unicode[STIX]{x1D64D}^{\text{T}}$ . The first EOF is plotted in figure 7 (open markers), again for stable equilibria at $P_{1}$ , $P_{2}$ and $P_{3}$ . The standard deviation for noise in sea level is taken $\unicode[STIX]{x1D702}_{sl}=1~\text{m}$ . For each EOF in figure 7, the corresponding explained variance is given by the solid lines in figure 8. An additional (dashed) eigenvalue curve in figure 8 shows the explained variance of the second dominant EOF for sea level noise.

Figure 8. Explained variance of EOFs due to noise in accumulation and sea level, i.e. rightmost eigenvalues of the covariance matrices, scaled with the total sum of the eigenvalues (trace). The open circles mark the points $P_{1}$ , $P_{2}$ , $P_{3}$ from left to right. Two eigenvalues are plotted for noise in sea level, as the variability is influenced by a second mode to the left of the first bifurcation $S_{1}$ .

Figure 9. (a) Ratio of the grounding line standard deviation due to noise in sea level $\unicode[STIX]{x1D70E}_{x_{g}}^{sl}$ and noise in accumulation rate $\unicode[STIX]{x1D70E}_{x_{g}}^{ac}$ , for several points along the branch before the first saddle-node bifurcation $S_{1}$ . The noise amplitudes are $\unicode[STIX]{x1D702}_{ac}=\bar{a}/10~\text{my}^{-1}$ and $\unicode[STIX]{x1D702}_{sl}=1~\text{m}$ . Ratios of the total variability and the ratios for velocity variability at the grounding line are plotted as well. (b) Flux perturbations at the steady grounding line, computed from the first EOF for both noise in accumulation rate and noise in sea level at several points along the branch before the first bifurcation.

The structure of the sea level noise induced EOFs is noticeably different from the eigenvector patterns of the Jacobian matrix, in particular at points $P_{1}$ and $P_{2}$ (left and centre panel in figure 7). Overall, the explained variance of the dominant EOF at the three positions is considerably less than in the accumulation case. The difference between the EOFs of both noise types is due to the spatial noise distribution in the Ornstein–Uhlenbeck process. The noise in accumulation is distributed homogeneously, whereas the noise in sea level is applied with a spatial pattern given by $-\tilde{\unicode[STIX]{x1D645}}_{12}\tilde{\unicode[STIX]{x1D645}}_{22}^{-1}B_{s}$ . With this spatial pattern a second EOF is excited that competes with the dominant EOF before the bifurcation, see figure 8.

4.3 Unique regime: stochastic variability

With the availability of covariance matrices at any stable steady state, we can investigate the effect of the two distinct noise sources on the grounding line variability. In figure 9(a), we plot the ratio of the grounding line standard deviations $\unicode[STIX]{x1D70E}_{x_{g}}$ due to noise in sea level and accumulation rate:

(4.13) $$\begin{eqnarray}\displaystyle r_{x_{g}}=\unicode[STIX]{x1D70E}_{x_{g}}^{sl}/\unicode[STIX]{x1D70E}_{x_{g}}^{ac}. & & \displaystyle\end{eqnarray}$$

Next to the grounding line noise response, we are interested in the total induced variability and hence depict the ratio of the covariance matrix traces

(4.14) $$\begin{eqnarray}\displaystyle r_{tr}=\mathop{\sum }_{i}(\unicode[STIX]{x1D70E}_{i}^{sl})^{2}/\mathop{\sum }_{i}(\unicode[STIX]{x1D70E}_{i}^{ac})^{2},\quad i=1,\ldots ,2N+1, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{i}^{2}$ is the dimensional variance of every state vector component. The ratio of the velocity standard deviations at the grounding line are given as well:

(4.15) $$\begin{eqnarray}\displaystyle r_{u}=\unicode[STIX]{x1D70E}_{u_{N}}^{sl}/\unicode[STIX]{x1D70E}_{u_{N}}^{ac}. & & \displaystyle\end{eqnarray}$$

The quantities are plotted along the bifurcation diagram up to the saddle-node bifurcation $S_{1}$ . For our amplitude choices $\unicode[STIX]{x1D702}_{ac}$ , $\unicode[STIX]{x1D702}_{sl}$ , the effect of noise in the accumulation rate on the variability in the grounding line position is considerably greater than that due to sea level noise, in agreement with the analysis in Hogg et al. (Reference Hogg, Shepherd, Gourmelen and Engdahl2016). In the vicinity of the saddle-node bifurcation, the ratio $r_{u}$ shows a local maximum due to sea level variability. As the total dimensional variability is mainly determined by the grounding line variance, the effect is not present in the ratio $r_{tr}$ .

The noise induced perturbations in grounding line flux and accumulation are computed from the relevant EOFs and plotted in figure 9(b). The sea level induced maximum in $r_{u}$ is related to the difference in noise induced flux perturbations $\hat{q}(x_{g})$ , which directly depend on the ice thickness at the grounding line. The flux perturbation under sea level variability $\hat{q}^{sl}(x_{g})$ remains large, close to the bifurcation, compared to its accumulation counterpart $\hat{q}^{ac}(x_{g})$ .

Note that the fluxes calculated from the EOF due to noise in accumulation mimic the perturbations in figure 6(b), since this mode follows the eigenvectors of the Jacobian matrix. Furthermore, analogous to the results in figures 5 and 6(b) and Schoof (Reference Schoof2007a ), we find that, at the bifurcation, the noise induced accumulation and grounding line flux perturbations behave according to the instability mechanism (Schoof Reference Schoof2012), with equality at the bifurcation.

Recall that figure 6 shows a significant change in the multiple equilibria regime when a GSLE is represented in the model. Using the Lyapunov techniques described in the previous section we can now investigate the stochastic properties of the different model configurations. For every computed stable point in the bifurcation diagram we solve the Lyapunov equation and determine the variances for noise in the accumulation rate and noise in the sea level variability. The results are shown in figures 10 and 11.

Figure 10. Comparison of the model statistics without (black) and with a gravitational sea level effect GSLE (red), with additive noise in the accumulation. (a) Grounding line standard deviation $\unicode[STIX]{x1D70E}_{x_{g}}$ at every stable point in the bifurcation diagram, where $L_{1},S_{1}$ and $L_{2},S_{2}$ mark the first and second bifurcation encountered when decreasing $A$ . (b) Grounding line flux standard deviation $\unicode[STIX]{x1D70E}_{uh}$ . (c) Overlay of the data in (a), such that the bifurcation points lie on the dashed line. (d) Overlay of the data in (b).

The introduction of a sea level equation causes a slight change in the stochastic properties under noise in accumulation rate. In general, as the multiple equilibria regime is shifted and decreased, an ice sheet will have a smaller grounding line variance and grounding line position for the same parameter value $A$ (cf. figure 6). Hence, the adaptation of the effective topographic slope causes the sheet to become more resilient under noise.

The standard deviations of the grounding line flux shown in figures 10(b) and 10(d) are computed from

(4.16) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D70E}_{uh}^{2}=\text{Var}(uh) & = & \displaystyle \text{Var}((\bar{u}+\tilde{u} )(\bar{h}+\tilde{h}))\end{eqnarray}$$

(4.17) $$\begin{eqnarray}\displaystyle & {\approx} & \displaystyle \bar{h}^{2}~\text{Var}(\tilde{u} )+\bar{u}^{2}~\text{Var}(\tilde{h})+2\bar{h}\bar{u}~\text{Cov}(\tilde{u} ,\tilde{h}),\end{eqnarray}$$

where we distinguish between steady ( $\bar{u},\bar{h}$ ) and stochastic components ( $\tilde{u} ,\tilde{h}$ ), and ignore higher-order terms. The variability in flux shows a decrease just before a bifurcation, followed by a drastic increase at the bifurcations. This in contrast to the grounding line variability, which shows a steady increase when approaching a bifurcation. The overlay of the bifurcations in figure 10(c,d) shows a lack of qualitative differences in the stochastic properties. However, for both model configurations we find a major asymmetry in the flux variability around the bifurcations, a feature that is not present in the grounding line variability.

Figure 11. Comparison of the model statistics without (black) and with a gravitational sea level effect GSLE (red), with additive noise at the sea level. (a) Standard deviation $\unicode[STIX]{x1D70E}_{x_{g}}$ at every stable point in the bifurcation diagram, where $L_{1},S_{1}$ and $L_{2},S_{2}$ denote the first and second bifurcation encountered when decreasing $A$ . (b) Standard deviation of the grounding line flux $\unicode[STIX]{x1D70E}_{uh}$ . (c) Overlay of the data in (a), the bifurcations lie on the dashed line. (d) Overlay of the data in (b).

The difference in variances between the models with the two different sea level representations in the case of noise in sea level are more pronounced. Both grounding line and flux variance show a decreasing trend as the sheet becomes larger under decreasing $A$ , see figures 11(a) and 11(b). An ice sheet will vary less under sea level noise as its size increases. The overlays in figures 11(c) and 11(d) show that the adaptation of the effective topographic slope causes a significant additional damping, which cannot be solely due to a difference in sheet size, since at $A/A_{0}\sim 1$ the sheets are equal in size for both model configurations. As noise in sea level is applied at the grounding line, the model response is likely to be more sensitive to changes in slope than the response to noise in accumulation. Note that, again, an asymmetry is present in the flux variability.

In general, the results in figures 10 and 11 show typical grounding line variabilities ranging from $\mathit{O}(100~\text{m})$ for noise in sea level, to $\mathit{O}(1000~\text{m})$ for noise in accumulation, depending on the dynamical regime. A continuation of steady states accompanied by the solution of the Lyapunov equations enables a straightforward separation of variability sources. This may be useful to explain observed grounding line migration (Hogg et al. Reference Hogg, Shepherd, Gourmelen and Engdahl2016), as will be elaborated in the discussion below.

4.4 Multiple equilibria: transition probabilities

For noise in accumulation we find roughly the same variability with a GSLE equation as in the case without such a representation. In contrast, variability due to noise in sea level is clearly dampened by the added GSLE equation. Here, the adaptation of the topographic slope appears to not only dampen the growth of perturbations, but also reduce the spread of the probability distributions along the branch. In figure 6, however, the distance between the stable branches is reduced in the GSLE configuration. This could have an effect on noise induced transitions in the multiple equilibrium regime, which we will investigate in this section. Although there is less sea level induced variability in the GSLE configuration, a smaller spatial distance between high ( $2N+1$ ) dimensional ellipsoidal confidence regions may indicate higher transition probabilities (Kuehn Reference Kuehn2012).

These confidence ellipsoids are determined by the probability density functions (PDFs) and indicate, with a certain confidence, where a state may be during a transient computation. The distance between ellipsoids associated with equilibria at different branches, but with the same parameter value, provides information about the relative frequency of noise induced transitions (Kuehn Reference Kuehn2012). Hence, calculating the distance between ellipsoids can be worthwhile when transient computations are expensive. Note, however, that such confidence regions are limited to representing the PDF under the assumption of linearised dynamics according to (2.22). For more information about these ellipsoids, see for instance Cowan (Reference Cowan1998).

Assuming we have computed a stable deterministic steady state $(\bar{\boldsymbol{x}},A)$ and the associated covariance matrix $\unicode[STIX]{x1D63E}$ , then the confidence ellipsoid can be represented as

(4.18) $$\begin{eqnarray}\displaystyle & \displaystyle E=\big\{\boldsymbol{x}\in \mathbb{R}^{2N+1}:\boldsymbol{v}^{\text{T}}\boldsymbol{x}\leqslant \boldsymbol{v}^{\text{T}}\bar{\boldsymbol{x}}+\sqrt{\boldsymbol{v}^{\text{T}}\tilde{\unicode[STIX]{x1D63E}}\boldsymbol{ v}}~\forall ~\boldsymbol{v}\in \mathbb{R}^{2N+1}\big\}, & \displaystyle\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D63E}}=Q_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D63E}$ is the positive semi-definite shape matrix associated with the confidence ellipsoid and $Q_{\unicode[STIX]{x1D6FC}}=3239$ is a scalar that can be obtained from the $\unicode[STIX]{x1D712}^{2}$ -distribution for a confidence level $1-\unicode[STIX]{x1D6FC}=0.683$ (Cowan Reference Cowan1998).

To determine the distance $d$ between ellipsoids, a convex minimisation problem of the form

(4.19) $$\begin{eqnarray}\displaystyle & \displaystyle -d=\min _{\Vert \boldsymbol{v}\Vert _{2}=1}\left(-\boldsymbol{v}^{\text{T}}\bar{\boldsymbol{x}}_{1}+\sqrt{\boldsymbol{v}^{\text{T}}\tilde{\unicode[STIX]{x1D63E}_{1}}\boldsymbol{ v}}+\boldsymbol{v}^{\text{T}}\bar{\boldsymbol{x}}_{2}+\sqrt{\boldsymbol{v}^{\text{T}}\tilde{\unicode[STIX]{x1D63E}_{2}}\boldsymbol{ v}}\right) & \displaystyle\end{eqnarray}$$

is solved. Here $\bar{\boldsymbol{x}}_{1}$ is an equilibrium at one branch and $\bar{\boldsymbol{x}}_{2}$ the equilibrium at the other branch, for the same parameter value $A$ . $\tilde{\unicode[STIX]{x1D63E}_{1}}=Q_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D63E}_{1}$ and $\tilde{\unicode[STIX]{x1D63E}_{2}}=Q_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D63E}_{2}$ are the respective shape matrices.

Figure 12. Distance patterns between the ellipsoids, calculated for the case without GSLE (black) and the case with GSLE (red), for noise at the accumulation (solid line) and at the sea level (dashed line). The noise amplitudes for the different sources ( $\unicode[STIX]{x1D702}_{ac}$ , $\unicode[STIX]{x1D702}_{sl}$ ) are amplified according to the amplifications in the transient results (see figure 13).

The results in figure 12 show that the distance between the ellipsoids is clearly reduced when approaching a bifurcation point, leading to larger transition probabilities (Kuehn Reference Kuehn2012). The shapes of the distances reveal asymmetries between the bifurcations. The case without GSLE shows a slight minimum at $L_{2}$ . The case with GSLE has a similar, but more pronounced asymmetry with a minimum distance at $S_{2}$ .

The main drawback of the ellipsoid distance computation is that it is based on local dynamical and stochastic behaviour of the model, combined with the global distance between branches. Possible nonlinear processes that affect the noise induced perturbations are, at best, only partially represented. Moreover, to compare the different model configurations, we restrict ourselves to normalised distance patterns in figure 12. Actual distance values have little meaning and serve merely as a test function for transition probabilities (Kuehn Reference Kuehn2012).

To verify the relative transition frequencies that the ellipsoid distances in figure 12 suggest, we perform transient computations using a backward Euler–Maruyama scheme, a well-established time integration technique for systems of SDAEs. Here, the deterministic part of (2.12) is treated implicitly, whereas the stochastic part is integrated explicitly with respect to a Brownian path. For a transition count we define a transition as a trajectory $\unicode[STIX]{x1D703}(t)$ from the upper branch to the lower branch and vice versa, similar to the definition in Kuehn (Reference Kuehn2012). Let $\boldsymbol{p}^{u}$ and $\boldsymbol{p}^{l}$ denote steady states in the upper and lower branch respectively and define some small tolerance $\unicode[STIX]{x1D70C}$ . A trajectory is counted as a transition from the upper to the lower branch $T^{u\rightarrow l}$ when it reaches a neighbourhood of $\boldsymbol{p}^{l}$ , $||\unicode[STIX]{x1D703}(t)-\boldsymbol{p}^{l}||_{2}<\unicode[STIX]{x1D70C}$ , after visiting a neighbourhood of $\boldsymbol{p}^{u}$ , $||\unicode[STIX]{x1D703}(t)-\boldsymbol{p}^{u}||_{2}<\unicode[STIX]{x1D70C}$ , within some time $T_{e}$ . $T^{l\rightarrow u}$ is defined similarly. In order to get a reliable number of transitions, we need to multiply the standard deviations $\unicode[STIX]{x1D702}_{ac},\unicode[STIX]{x1D702}_{sl}$ by amplification factors $f_{ac},f_{sl}$ , respectively. In the case of noise in the accumulation we choose $f_{ac}=5$ , for the sea level noise we choose $f_{sl}=20$ . The results are shown in figures 13(a) and 13(b).

Figure 13. Sum of the transitions from the upper to the lower branch and vice versa, relative to the total number of transients $N_{trans}$ , i.e.  $(\#T^{u\rightarrow l}+\#T^{l\rightarrow u})/N_{trans}$ . (a) Transitions under noise in accumulation with GSLE (red) and without GSLE (black), using a backward Euler–Maruyama time integration with time step $\unicode[STIX]{x0394}t=1000~y$ , end time $T_{e}=100\,000~y$ and amplified standard deviation $f_{ac}\unicode[STIX]{x1D702}_{ac}$ , $f_{ac}=5$ . (b) Transitions under noise in sea level with GSLE (red) and without GSLE (black), with amplified noise $f_{sl}\unicode[STIX]{x1D702}_{sl}$ , $f_{sl}=20$ . The total number of transient experiments for each point on the branch is $N_{trans}=960$ . A neighbourhood is defined in terms of grounding line positions, where the tolerance is chosen to be $\unicode[STIX]{x1D70C}=10~\text{km}$ .

Although the distance between the probability distributions shows a decrease, this is not sufficient to compensate for the decrease in variability. Hence, the number of transitions is significantly reduced in the model with GSLE. This holds for both noise in accumulation rate (figure 13 a) and noise in sea level (figure 13 b), although the reduction of variability is more drastic in the sea level noise case.

From the results in figure 13 we may conclude that it is more likely to travel from the upper branch (large sheet) to the lower branch (small sheet) than vice versa, which holds for both noise types. The asymmetry in figure 12 for both model configurations is similar. As described above, the same asymmetry can be found in the flux variability. Figures 10(d) and 11(d) show that the flux variability near $L_{2}(S_{2})$ at the upper branch is generally higher than at the lower branch near $L_{1}(S_{1})$ , which is likely due to the weaker effective topographical slope at $L_{2}(S_{2})$ .