Numerical Bifurcation Analysis of Marine Ice Sheet Models

Numerical Bifurcation Analysis of Marine Ice Sheet Models

T. E. Mulder T. E. Mulder and H. A. Dijkstra Institute for Marine and Atmospheric research Utrecht, Department of Physics, Utrecht University, the Netherlands, 22email: h.a.dijkstra@uu.nlF. W. Wubs Johann Bernoulli Institute for Mathematics and Computer Science, Groningen University, Groningen, The Netherlands 44email:    H. A. Dijkstra and F. W. Wubs T. E. Mulder and H. A. Dijkstra Institute for Marine and Atmospheric research Utrecht, Department of Physics, Utrecht University, the Netherlands, 22email: h.a.dijkstra@uu.nlF. W. Wubs Johann Bernoulli Institute for Mathematics and Computer Science, Groningen University, Groningen, The Netherlands 44email:

The climate variability associated with the Pleistocene Ice Ages is one of the most fascinating puzzles in the Earth Sciences still awaiting a satisfactory explanation. In particular, the explanation of the dominant 100 kyr period of the glacial cycles over the last million years is a long-standing problem. Based on bifurcation analyses of low-order models, many theories have been suggested to explain these cycles and their frequency. The new aspect in this contribution is that, for the first time, numerical bifurcation analysis is applied to a two-dimensional marine ice sheet model with a dynamic grounding line. In this model, we find Hopf bifurcations with an oscillation period of about 100 kyr which may be relevant to glacial cycles.

Keywords: Marine Ice Sheets, Bifurcation Analysis, Multiple Equilibria, Oscillatory Modes

1 Introduction

Very detailed information on past temperatures on Earth has been obtained from marine benthic records, in particular oxygen isotope ratios. Water in ice cores contains two isotopes of oxygen, O and O. The normalized isotope ratio O is calculated as a deviation from a reference sample as


where the reference sample is ‘standard mean’ ocean water. The isotope is lighter than so that water containing is preferentially evaporated and a temperature-dependent fractionation occurs. Changes in reflect the combined effect of changes in global ice volume and temperature at the time of deposition of the sampled material. During very cold conditions, global ice volume is relatively large and hence sea level is low, which enriches water in the ocean with . Also because of the colder temperatures, more remains in the ocean and less becomes locked in the ice. Hence, the ratio in ice cores will decrease (becomes more negative) under colder conditions.

In marine sediment cores the opposite behavior of the isotope ratio is found and increases (becomes more positive) when the temperature decreases (during colder conditions, the concentration of the heavier isotope will increase). In Fig. 1, a time series is shown of a composite O ocean sediment (benthic) record over the last 2 Myr Lisiechi2005 (). A cooling trend is found on which variability in ice cover is superposed. Analysis reveals that this variability is first dominated by a 41 kyr period and after the so-called Mid Pleistocene Transition (MPT) at about 700 kyr, it is dominated by a 100 kyr period.

Figure 1: (a) The LR04 benthic stack over the Pleistocene, constructed by the graphic correlation of 57 globally distributed benthic records (data from Lisiechi2005 ()). The x-axis indicates time BC in kyr (so the time is increasing from right to left).

The European Project for Ice Coring in Antarctica (EPICA) has provided two deep ice cores in East Antarctica from which climate conditions can be reconstructed back to 800 kyr BP Jouzel2007 (). From the reconstructed temperature anomaly time series (Fig. 2), one observes the asymmetry between the slow glaciation and the rapid deglaciations. Glacial-interglacial transitions have affected all components of the climate system and induced relatively large amplitude changes of many variables in these components. One important player in the climate system responsible for the globalization of these transitions is believed to be the atmospheric CO concentration. A composite CO record is also shown in Fig. 2, created from a combination of records from the Dome C and Vostok ice cores. It is observed that the atmosphere CO concentration varies from about 180 ppm to 280 ppm during a glacial-interglacial transition and that an optimal correlation with the O time series occurs near lag zero.

These results lead to many intriguing questions: Why did glacial-interglacial cycles appear in the Pleistocene? Which processes in the climate system caused the glacial-interglacial changes in global mean temperature and ice sheet extent? What caused the transition (the MPT) from the 41 kyr world to the 100 kyr world about 700 kyr ago?

Approaches to answers on the Pleistocene Ice Ages problem have a very interesting history which is nicely described in Imbrie1986 (). A connection with the orbital characteristics of the Earth-Sun system was already made in the 19 century, but in the 1930s it was suggested Milankovich1930 () that glaciations occur when the insolation intensity is weak at high northern latitudes during summer. When the 65N insolation is small, ice can persist throughout the year leading to the growth of ice sheets. Favorable conditions for this to happen are when the spin axis is less tilted and the aphelion (the point in the orbit, where the Earth is farthest from the Sun) coincides with summer in the Northern Hemisphere.

The variations in insolation are caused by the changes in orbital characteristics of the Earth, and there are three types of motion relevant for the amount of radiation received at a particular point on Earth. First, the spin axis of the Earth undergoes precession. One full cycle of precession has a period of 27 kyr, but coupled to the movement of the long end of the ellipse around the Sun (in 105 kyr) the net effect is a fluctuation in solar radiation with a period of 23 kyr. In addition, both the obliquity and the eccentricity of the Earth’s orbit undergo periodic variations. The tilt angle changes in 41 kyr between 22 and 24 leading to variations in seasonal contrast, and the eccentricity varies from 0.0 (perfect circle) to about 0.05 with periodicities of 100 kyr and 450 kyr.

Figure 2: Reconstructed temperature (drawn) and atmospheric CO (dotted) concentration from ice cores on Antarctica.

A time series of the insolation at 60N (Fig. 3a) clearly shows variations of about 100 Wm over the last 1 Myr. A comparison of the spectrum of this insolation curve (Fig. 3b) and the spectrum of a record from an ocean sediment core (Fig. 3d) over the last 1Myr (as shown in Fig. 3c) shows that there are clear signatures of the 19 and 23 kyr precession and of the 41 kyr obliquity variations of the Earth’s orbit in the Ocean Drilling Program record. On the other hand, at the 100 kyr time scale, there is hardly any forcing amplitude while the climate signal in the record has the largest amplitude.

(a) (b)
(c) (d)
Figure 3: (a) June insolation at 60 N. (b) Spectrum of the time series in (a). (c) Time series of at ODP677 (83W, 1N) over the last 1 Myr. (d) Spectrum of the time series in (c). Note that , and .

Comparing O records with the insolation time series Raymo2008 (), it is interesting that the amplitude of the ice-cover variations can be very large while the insolation variation is very small. Furthermore, the variations in insolation provide no clue on the transition from the 41 kyr world to the 100 kyr world as they have the same temporal characteristics through the transition.

Although it is clear that the orbital insolation variations must play a role, a simple linear forcing-response relation apparently does not apply. The 100 kyr variations in insolation due to eccentricity are very weak (and it is the only forcing with an annual mean signal) and the 41 kyr and 23 kyr provide only low-frequency variations on the seasonal variations, not on the annual mean insolation. Hence, processes internal to the climate system must play a role in the amplification of the orbitally induced insolation variations.

2 Basic theories of interglacial-glacial cycles

There have been many suggestions on the dominant mechanisms of ice-age variability. The orbital variations in insolation over the globe are at the heart of all these theories. Traditionally, the June insolation at 65N (such as shown in Fig. 3a) has been used as the most important part of this forcing as this determines whether snow will be left at the end of the Northern Hemisphere Summer season. Below, we will refer to this orbital component of the insolation as the M-forcing.

In a first theory, the behavior of the climate system is considered as a transient deviation from a single steady equilibrium due to the M-forcing. One can calculate that a 1% change in solar insolation leads to about a 1C temperature change. The dominant amplitude of eccentricity variations is at a period of about 400 kyr while the next strongest variations occur on about a 100 kyr time scale. The variations in eccentricity do modify the globally and annually-averaged amount of insolation but the amplitudes are very small, in the order of 0.1% of the solar constant. The variations in eccentricity can therefore only account for a climate signal of at most 0.1 K. A similar analysis gives that the total variations in precession and obliquity can only account directly for a signal of at most 0.5C, about an order of magnitude too small. When considering many other processes (e.g., ice sheets, bedrock) in the direct response to the M-forcing, the sensitivity does not increase enough to explain the climate signal Ghil1994b ().

There have been several suggestions that the existence of multiple steady states gives, together with the M-forcing, rise to glacial cycles. For example, in the model by Paillard1998 () there are three equilibrium states in the climate system: an interglacial state i, a weak glacial state g, and a strong glacial state G. Transitions from i to g occur when the summer insolation at 65N drops below a value . Furthermore, transitions from g to G occur when the ice volume increases above some critical level . Finally, a transition from G to i occurs when the insolation increases above a level . These are the only transitions which are allowed in this model. By forcing the model with the M-forcing, there is a good overall agreement between model and observations (considering the simplicity of the model). By allowing for a slight linear trend in the maximum ice volume and one in the insolation forcing, Paillard1998 () also finds the MPT at around the correct time and the spectra of his model and typical data correspond reasonably well to each other.

When multiple equilibria, a weak periodic forcing, and noise are present there is also the possibility of stochastic resonance. In fact, the discovery of stochastic resonance actually occurred Benzi1982 (); Benzi1983 (); Nicolis1982 () while trying to explain the 100 kyr dominant glacial cycles using an energy balance atmosphere model. The central element of this theory is the amplification of the weak periodic eccentricity component of the M-forcing by noise in the presence of multiple equilibria.

Coupled atmosphere-cryosphere and cryospheric-lithosphere processes can also give rise to internal variability. The question is of course, whether these processes can also give rise to sustained oscillations (through Hopf bifurcations). If so, the time scales and amplitude ranges of these oscillations with ‘realistic’ values of the parameters will be of interest regarding the glacial cycles. Typical results can be found in Ghil1987 () where a sustained oscillation is indeed present under steady forcing. The time scale of the oscillation is about 6-7 kyr and not 100 kyr which shows that the processes captured in this simple model are not able to generate this long time scale. In LeTreut1988 (), the model in Ghil1981 () is forced by an M-type forcing and the 100 kyr period arises due to nonlinear resonances of the external frequency of the forcing and the internal frequency of oscillation Ghil1994b ().

Many other idealized models have been proposed which involve other components of the climate system Saltzman2001 (), for example, those involving the atmospheric concentration of CO and the global ocean state. The model proposed by Paillard2004 () attributes to the Antarctic ice sheet extent a central role in linking climatic and CO glacial-interglacial changes. The model proposed by Omta2013 () investigates the role of marine calcifiers in glacial-interglacial cycles. For many of these models, bifurcation diagrams have been computed showing that oscillations are associated with Hopf bifurcations Crucifix2012 ().

A relatively simple model, where an internal oscillation exists with a time scale of about 100 kyr is that of Gildor2000 (). The model is a box model of the climate where a similar atmospheric-cryospheric model as in the previous subsection is coupled to an ocean model with a sea-ice component. In Gildor2001 (), results of the model forced by constant annual mean insolation (no seasonal or Milankovitch forcing) are presented to assess the degree to which the internal processes (particularly sea-ice) may control glacial cycle variability. A typical result for (near) standard values of the parameters in Gildor2001 () shows that oscillations with a time scale of about 100 kyr are found. The proposed mechanism of the variability is referred to as the sea-ice switch, where rapid sea-ice growth and decay can act as a switch for the precipitation-temperature feedback affecting the growth and/or decay of ice sheets. The 100 kyr time scale is due to the growth and decay of ice sheets which is coupled by relatively rapid sea ice changes.

When this model is forced with Milankovitch insolation changes, nonlinear resonances may occur between the internal oscillation and the orbital forcing leading to time series which qualitatively resemble the observed records Tziperman2006 (), which demonstrates how phase locking to Milankovitch forcing affects glacial cycles in this idealized model. These nonlinear resonances are likely to be present in every model where a strong nonlinear interaction is represented, explaining for example the good agreement between very conceptual models, where only a multiple state switch Paillard1998 () is represented, and observations. In other words, even if the mechanism of the glacial-interglacial variability is incorrect, there may still be a good fit with the isotopic record. Due to the synchronization, a comparison between time series of simple models and isotope records is not mechanistically selective Crucifix2016 ().

Although many bifurcation studies have been done on low-order models, there appear to have been no studies where numerical bifurcation analysis has been applied to spatially extended models of ice sheets. In this contribution , we make a first step in this direction, by looking for multiple equilibria and Hopf bifurcations in a two-dimensional model of a marine ice sheet.

3 Methodology

Consider in Fig. 4 a two-dimensional marine ice sheet situated on a bedrock topography in a Cartesian coordinate system. The ice sheet and bedrock are taken symmetrical, with a symmetry axis at . The grounding line is indicated by , the ice thickness by , and the bedrock by .

Figure 4: A two-dimensional marine ice sheet. The ice thickness is given by and the horizontal ice velocity by . Up until the grounding line , the ice sheet rests on the bedrock . At the ice sheet extends into a floating ice shelf.

3.1 Model

An introduction to ice-sheet and glacier modeling can be found in Greve_2009 () and Veen_2013 (). The dynamics of a marine ice sheet is modeled using the shallow-shelf approximation (SSA), which is obtained by simplifying the full Stokes problem for gravity driven ice flow Greve_2009 (). The two-dimensional SSA, as implemented in Schoof_2007 (), is used as a benchmark problem for the marine ice sheet intercomparison project (MISMIP, MISMIP_2012 ()). Since a thorough comparison with other results is available, this will be our model of choice for the bifurcation analysis. Conservation of mass gives:


where is the ice thickness, the ice velocity and the accumulation rate. On a downward sloping bed (Fig. 4), the accumulation and the ice flux at the grounding line are in equilibrium: a positive perturbation of the grounding line increases both the accumulation and the flux, leading to a zero net ice growth.

Conservation of momentum gives:


where and are coefficients of Glen’s flow law, a constitutive relation describing the rheology of ice (Veen_2013 (), Greve_2009 (), typically ). The ice density is given by , is the gravitational acceleration and the bedrock taken positive in the downward direction. The consecutive terms in the momentum balance (3) represent longitudinal stresses, vertical shear stresses, and the driving stress respectively. The parameters and determine the sliding of the ice. Together, and effectively describe the size of the transition zone, i.e., the region in which the grounded sheet becomes afloat and transforms into an ice shelf.

The left boundary of the domain is located at an ice divide, a location in the ice with zero horizontal flow. Hence, we take at . As we assume symmetry around , we also require that


We will denote the grounding line by . At the ice sheet becomes afloat and the following flotation condition holds:


From an integration of the shelf flow Schoof_2007 (), an extra condition at the grounding line is obtained:


Note that the profile of the shelf is not given by the SSA model, but obtained using an equilibrium analysis Veen_2013 ().

3.2 Non-dimensional equations

The first difficulty one encounters is the unknown right boundary of the problem given by the grounding line position . As discussed in Schoof_2007 () and Vieli_2005 (), a moving grid approach can be used to track . Using a transformation , the original domain is mapped onto the fixed domain . As a result, the problem now has three unknowns: , and . The differential operators are transformed using the chain rule:

where and denote the independent variables in the transformed domain. Since we only transform in space we have that . The transformations of (2)-(6) are then given by:


To improve numerical accuracy the equations are non-dimensionalized. Let


with typical thickness m, horizontal extent m and typical timescale y. Substituting these expressions into (7) gives


where we let . Similarly, the non-dimensionalized version of (8) is given by


where we introduce the new constants


Finally, at the boundaries we obtain




3.3 Numerical implementation

From here on we omit the tildes and assume the unknowns are non-dimensional. As in Schoof_2007 (), the domain is discretized using a staggered grid with a fixed mesh-width: . The left boundary is taken at the vertex and the right boundary at the cell center . Thus, for , we have vertices at and cell centers at . The discretized solution values for ice thickness are located at the vertices , while the values for ice velocity are positioned at the cell centers .

The transformed and non-dimensionalized continuity equation (13) is discretized using a central difference for the stretching and an upwind discretization for the flux:


At the left boundary symmetry requires

Using these expressions we can resolve the dependence on nonexistent grid-points. For and , mass conservation is therefore given by


Note that at the rightmost vertex (), the right hand side of the discretized mass conservation (20) does not contain any dependencies on nonexistent grid-points. In the left hand side we will need to use a one-sided difference for the stretching term.

Define . The momentum conservation (14) is discretized using central differences:


At the left boundary we let . At the right boundary we impose the following discretization of (18) with a substituted flotation condition (cf. (17)):


The discretization contains unknown , unknown and an unknown grounding line position . To achieve a closed system of equations, the flotation criterion at the cell center is prescribed using an extrapolation of the thickness, which gives the closing requirement:


Finally we obtain a problem of the form


The unknown functions are discretized: . The real-valued matrix and nonlinear operator are given by


Here, is the discretization of the stretching in the left-hand side of Equations (20)-(22). is given by the right-hand side of discretizations (20)-(22). Similarly, and are given by the right hand sides of the discretized momentum equation and flotation criterion (23)-(24).

3.4 Pseudo-arclength continuation

The discretized equations give a problem of the form


where and are linear and non-linear operators respectively. We explicitly introduce the parameter dependence since we are interested in solution branches satisfying . For example, our first parameter of interest will be the temperature, which is present in the coefficient in Glen’s flow law.

Various continuation techniques exist to trace a stationary solution branch while varying a parameter. A successful approach is to parameterize a solution branch with a pseudo-arclength parameter : and impose an approximate normalization condition on the tangent, to close the system of equations: , where is an initial known stationary solution, the tangents w.r.t. the arclength parameter at and a specified step size Dijkstra_2005 (); Kuznetsov_2004 ().

To find a new point on the solution branch a predictor-corrector method is used. A suitable tangent predictor is given by

Note that the prediction is denoted by , whereas an actual new solution will be denoted by . The correction onto the solution branch is made through the solution of the nonlinear system given by


A Newton-Raphson root finding procedure, initialized with the prediction , gives the following iteration:


where , and is the Jacobian matrix of . If this iteration converges a new stationary solution has been found. At a fold bifurcation the Jacobian matrix will have a zero eigenvalue, yet the system in (31) remains non-singular, and the continuation is able to trace the solution branch into its unstable domain.

When a stationary solution , satisfying


has been found, its stability can be investigated with a perturbation and a Taylor expansion around the stationary solution:


Solutions of (33) are of the form . Substitution gives a generalized eigenvalue problem:


The stability of a stationary solution depends on the sign of the real part of the eigenvalues. If we find an eigenvalue with positive real part the solution contains a growing mode and is thus unstable.

4 Results

4.1 Bifurcation diagram

First we consider the case of constant accumulation, represented by a constant , together with a fixed bedrock . The bedrock profile is chosen as in Schoof_2007 ():


with . To find a stationary solution, an ice sheet surface profile of the form


is used as initial guess for the Newton-Raphson iteration solving . Typically, and is chosen such that the flotation criterion is satisfied: , where is the bedrock at the grounding line. The analytic profile contains a steep gradient at the grounding line, which is essential for a good convergence of the root finding process.

Figure 5: One-parameter bifurcation diagram (left) and solutions (right), decreasing the parameter in the SSA model which corresponds to a decrease in temperature and an increase in ice growth. The bedrock contains an upward slope which admits multiple steady states ((b),(d) and (f)) for a constant parameter. Eigenvalue analysis shows two saddle-node bifurcations: a single eigenvalue crosses the imaginary axis to the positive right half plane at (c) and returns to the left half plane at (e). The number of grid-points is ; the other parameters are given in Table 1.

After an initial equilibrium profile is obtained for , a pseudo-arclength continuation traces the solution branch in the direction of decreasing , see Figure 5. The presented bifurcation diagram confirms the multiple equilibria regime associated with the hysteretic behavior shown in Schoof_2007 (). For a fixed value of the parameter , three equilibria are distinguished and marked as (b), (d) and (f). At the point (c) a saddle-node bifurcation occurs and an eigenvalue is observed to cross the imaginary axis to the right half plane. At (e), the same eigenvalue returns to the right half plane through a second saddle-node bifurcation. The values of the parameters are summarized in Table 1.

Table 1: Parameter values for the experiment in Figure 5, similar to the values chosen in Schoof_2007 () and MISMIP_2012 ().

4.2 Numerical accuracy

In order to investigate the numerical accuracy of the discretization and continuation methodology we perform a convergence experiment, see Table 2. Let the error in the approximated bifurcation point be proportional to a power of the mesh-width:


where denotes the actual bifurcation point and and are constants. We define a difference between subsequent mesh-halvings and let . Then, the ratio between consecutive differences only depends on the power :

Table 2: Convergence of the first bifurcation (point (c) in Figure 5).

From Table 2 we suspect as becomes large, corresponding to . The scheme must therefore be of first-order accuracy, which is undoubtedly due to the first-order upwind discretization in the continuity equation (20). Unfortunately, the upwind discretization of the ice flux is essential to the stability of the scheme. We conclude that a higher order upwind scheme should be considered in (20).

4.3 Mechanism

The advantage of the approach chosen here is that the spatial patterns of perturbations destabilizing the marine ice sheet can be determined from the eigenvectors in (34).

Figure 6: Normalized components of the eigenmode that becomes unstable. The corresponding steady states are located at the points (b)-(f), described in Figure 5 and Table 1. We distinguish between a perturbation pattern related to sheet thickness and a pattern related to ice velocity . 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 normalized spatial pattern of the evolution , together with normalized components and .

For the unstable equilibrium (d) in Figure 5 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 (34) and is depicted for the steady states (b),(c),(d),(e) and (f) in Figure 6. 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 , with the scalar grounding line perturbation. An eigenvector with corresponding eigenvalue and gives the destabilizing pattern for unstable ice sheet retreat. By adjusting the sign of the eigenvector, such that , we restrict our exposition to destabilizing patterns for unstable ice sheet growth.

At the grounding line , the perturbation of the unstable steady state (d) shows a slight decrease in ice thickness, while at (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 (d), indicating interior ice growth due to an imbalance between global accumulation and ice flux at the grounding line.

The velocity perturbation at the grounding line shows an increase for stable states and a clear decrease for the unstable state (d). Together with the negative perturbation in thickness this implies that, at (d), there must be a decrease in flux across the grounding line for a positive perturbation . 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 negative, giving a net ice loss and a retreat from equilibrium.

The continuation approach allows an efficient computation of flux perturbations using (34). From a linear stability analysis of (2) with perturbation around an equilibrium we obtain an evolution equation for the thickness perturbation :


where we neglect higher order terms. At the unstable steady state (d) in Figure 6, 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 (black dash-dotted line)and the advected thickness perturbation (red dotted line). The latter clearly dominates the instability in (d). Note that at the bifurcation points (c) and (e) the components of the perturbation flux have a compensating effect.

To investigate how a perturbation changes from stable to unstable through the saddle-node bifurcation , we compute the accumulation and grounding line fluxes. The steady state gives a balance:


(a) (b)

Figure 7: Perturbations in grounding line and accumulation flux against the parameter (a) and as a function of (b), together with the bedrock slope. The first and second saddle-node bifurcations are marked and . Dashed lines correspond to growing perturbations of unstable steady states. The perturbations and correspond to a positive grounding line perturbation . At and , the grounding line bedrock slope is mm/km and mm/km respectively. The maximum slope in the unstable regime is 620 mm/km. Note that grounding line flux perturbations depend on the steady grounding line position , whereas accumulation flux perturbations depend on the grounding line perturbation .

In Figure 7 we show perturbations of the balance (41) at the grounding line. A perturbation in accumulation flux is given by , a grounding line flux perturbation by in (40). The perturbations are plotted against (Figure 7a) and (Figure 7b), together with the bifurcation points and the bedrock slope. At the first saddle-node bifurcation , the flux becomes smaller than the accumulation . Beyond this point, a change in accumulation due to is not balanced by the grounding line flux and, hence, the perturbation changes from damped to growing. At , becomes greater than and the mode is damped again.

In Figure 7b we also plot the bedrock slope, taken positive when the bed elevation increases with , that is


with as in (35). Note that the sign switch in coincides with the sign switch in the bedrock slope. The grounding line flux will increase for a positive between the bifurcation 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 7 confirms that an eigenvalue of the ‘full model’ in Schoof_2007 () becomes positive when .

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_2012 (). In addition, the eigenvectors reveal destabilizing interior patterns with, most interestingly, interior thickness anomalies and their advection. These turn out to play a major role in the unstable growth and retreat of the ice sheet.

4.4 Glacial isostatic adjustment

The simplest model describing the interaction between an ice sheet and the underlying bedrock is a local lithosphere, relaxing asthenosphere (LLRA) model. An equilibrium argument Greve_2009 () gives:


with the density of the asthenosphere, the equilibrium bedrock and the initial, load-free bedrock. Due to the highly viscous asthenosphere the equilibrium bedrock is reached after a significant response time . The evolution of the bedrock can be modeled using


with the actual bedrock profile.

Adding this elastic bedrock to the SSA model means introducing a new unknown and a new discretized differential equation. To remain consistent with the implementation in Section 3.3, we need to perform the transformation , giving


Non-dimensionalizing (45) is straightforward. Using a central difference for the stretching we obtain the following discretization:


Symmetry at the left boundary gives a vanishing spatial derivative. At the grounding line the pressure exerted by a column of ice equals that of the column of water: . Substituting the flotation criterion gives the following discretization for the right boundary of the bedrock equation:


Equation (47) acts as a Dirichlet boundary condition in the stationary case. Note that, with this boundary, we assume the bedrock is adjusted regardless of the type of material that is on top of it. Hence, the initial load-free bedrock cannot be subject to a water load. At the grounding line, ice and water exert the same pressure such that the bedrock extends continuously into its submerged part.

Again we let the problem have the form


but now with , , and , given by


The functions and give the discretizations of the stretchings in the left hand side of the continuity and bedrock equation, gives the discretization of the elastic bedrock equation (46). The other functions are the same as in Section 3.3, except with their dependence on made explicit.

Instead of having a constant accumulation , as in vdBerg_2006 () we will give it a linear dependence on the surface height and a ceiling :


where is an accumulation gradient and the equilibrium height. Below we have mass loss, i.e. ablation, the opposite of accumulation. The balance between accumulation and ablation is governed by air temperature, which we assume to decrease with increasing surface height.

The implementation of the height-dependent accumulation is straightforward, mainly requiring a few extra dependencies on and in the Jacobian matrix. To obtain a smooth transition from the linear function to the ceiling , we use an approximation to the Heavyside function:


with moderately small.

Recall the bedrock as given in (35). Now that we have implemented the isostatic adjustment, we need to define a smooth original bedrock that is not subject to a water load. Eventually, the bedrock will only be partially submerged, but adjusting only within the submerged sub-interval introduces unfavorable discontinuities. For that reason we will obtain using a global adjustment.

Consider the stationary case of (44). Replacing the ice load with a water load gives the required adjustment:


The added components open up a multitude of possible continuation parameters. Here we will restrict ourselves to a continuation in the equilibrium height , to pursue oscillatory solutions due to the load accumulation feedback Ghil_1994 (): as ice thickness increases due to accumulation, the surface height may decrease as a result of isostatic adjustment (with a delay ), leading to a decrease in accumulation. This feedback suggests the possibility of oscillatory solutions within a certain parameter regime.

A Hopf bifurcation occurs when a steady periodical solution emerges from a fixed point. In the spectrum given by the eigenvalue analysis, a Hopf bifurcation corresponds to a complex conjugate pair crossing the imaginary axis from the left to the right half plane. The result of a continuation in is shown in Figure 8. An additional list of parameters is given in Table 3.

Figure 8: Bifurcation diagram (left) and solution (right) of a continuation in equilibrium height , starting at . A Hopf bifurcation is detected at (a): and the corresponding solution is plotted on the right. Instances of the oscillatory perturbation corresponding to the complex conjugate eigenpair are plotted in Figure 9. The number of grid-points is , the other parameters are given in Table 3.
asthenosphere relaxation timescale
   asthenosphere density
   maximum accumulation
accumulation gradient
Glen’s flow law rheology parameter
Table 3: Parameter values for the experiment in Figure 8. The accumulation ceiling and gradient are chosen similar to vdBerg_2006 (). Parameters not mentioned here remain equal to the ones given in Table 1.
Figure 9: Oscillating perturbation in ice thickness (left) and velocity (right), given by the eigenvector corresponding to the complex conjugate eigenpair at the Hopf bifurcation . Several instances of the perturbation are obtained using (53). The angular frequency is given by .

A Hopf bifurcation is detected at , with complex conjugate eigenpair . The corresponding complex eigenmode describes the perturbation destabilizing the solution. The real part of the perturbation gives a disturbance profile Dijkstra_2005 ():


At the bifurcation we have , then and give two instances of the oscillatory perturbation: at and at , see Figure 9. The non-dimensional period is given by . As the timescale in the experiment is taken , the dimensional period is .

The obtained oscillating perturbation demonstrates the dynamics of the load accumulation feedback. For the ice thickness, the perturbation in the interior seems to be constantly ahead of the perturbation at the grounding line. Note that the region between the ice divide and the point at which the surface attains the equilibrium height () is subject to a net accumulation. From the perturbation in thickness we see that the feedback is clearly driven by the accumulation.

Beyond there is a net ablation. In this region the perturbations in thickness and velocity show a drastic change in behavior. A steep increase in velocity occurs when the ice thickness decreases (), which corresponds to a large increase in flux to facilitate the mass loss. Similarly, when the sheet grows (), the flux needs to reduce to facilitate growth.

Figure 10: Oscillating perturbation in the bedrock, given by the eigenvector corresponding to the complex conjugate eigenpair at the Hopf bifurcation.

In the ablation region the perturbation in ice thickness shows a few peculiar oscillations. At the grounding line the sheet thickness seems to be slightly less than at the ice divide, perhaps to facilitate the appropriate flux. However, just before the grounding line and after the effect of a reverse accumulation feedback seems to be visible. Due to the negative accumulation the bedrock is adjusted differently, see Figure 10. In the oscillatory bedrock perturbation there seems to be some irregular influence from the free boundary. Nevertheless, a change in adjustment around is visible in at least two phases of the oscillation at and .

5 Summary and Discussion

From the theories of ice-age cycles, it is clear that ice-sheet dynamics plays a central role in the explanation on how the variations in insolation lead to multi-millennial variability of the climate system, in particular on the 100 kyr time scale. A problem with the current theories is that there are many different conceptual models which can give a dominant 100 kyr variability, but it is difficult to falsify them based on the proxy data record Crucifix2016 ().

A step forward is to determine spatial patterns of variability associated with the glacial cycles, similar to what has been done for other problems of climate variability such as El Niño Dijkstra2002_ARFM () and the Atlantic Multidecadal Variability Frankcombe2010b (). This obviously requires a next level of models in the hierarchy (see Chapter 6 of Dijkstra2013B ()), at least formulated in terms of partial differential equations (spatially extended models), also (often) referred to as intermediate complexity models (ICMs).

In this contribution, we applied techniques of numerical bifurcation theory to study the bifurcation behavior of solutions of such an ICM two-dimensional ice sheet model Schoof_2007 (). The complication arises here from the grounding line dynamics at , which is in principle a free boundary problem. Here, it is handled by using a transformation , where the original domain is mapped onto the fixed domain .

In the version of the model where the bottom topography is fixed, the results provide insight into the mechanism of marine ice sheet instability. There are robust intervals in parameter space, where multiple equilibria occur, corresponding to a large ice sheet and a small ice sheet. Hence, for these parameter values, transitions can occur where the ice-sheet decreases substantially in shape and length. Tracing unstable equilibria and performing a stability analysis enables the investigation of the actual structure of the growing perturbation. We have shown numerically that a positive eigenvalue is associated with the instability criterion Schoof_2012 () for the full problem, and that the advected thickness perturbation (the term , where is the steady state velocity and the ice thickness perturbation) dominates the instability process.

Transition behavior under the effects of noise in the accumulation parameter on the grounding line motion under stable conditions have been studied in Mulder2018 (). The magnitude (for typical noise levels) of these motions is in the order of 1000 m, which are similar amplitudes as observed for example in Hogg_2016 (). It appears to be more likely to jump from a large ice sheet state to a small ice sheet state than vice versa. Grounding line flux variability shows a related asymmetry, likely due to differences in local bedrock conditions and/or global ice sheet extent.

When the load-accumulation feedback is included by extending the ice-sheet with a dynamical bottom topography, oscillatory instabilities occur through a Hopf bifurcation. Here, the eigenvectors associated with the instability provide an interesting spatial pattern of variations of the ice sheet, with largest amplitudes near the grounding line. The time scale here is connected to the relaxation time scale