Modelling dynamic ice-sheet boundaries and grounding line migration using the level set method

Abstract

Computing predictions of future sea level that include well-defined uncertainty bounds requires models that are capable of robustly simulating the evolution of ice sheets and glaciers. Ice flow behaviour is known to be sensitive to the location and geometry of dynamic ice boundaries such as the grounding line (GRL), terminus position and ice surface elevation, so that any such model should track these interfaces with a high degree of accuracy. To address this challenge, we implement a numerical approach that uses the level-set method (LSM) that accurately models the evolution of the ice–air and ice–water interface as well as capturing topological changes in ice-sheet geometry. This approach is evaluated by comparing simulations of grounded and marine-terminating ice sheets to various analytical and numerical benchmark solutions. A particular advantage of the LSM is its ability to explicitly track the moving margin and GRL while using a fixed grid finite-difference scheme. Our results demonstrate that the LSM is an accurate and robust approach for tracking the ice surface interface and terminus for advancing and retreating ice sheets, including the transient marine ice-sheet interface and GRL positions.

1. Introduction

The Greenland and Antarctic ice sheets have been losing mass at an accelerated rate (Bevis and others, 2019; Rignot and others, 2019) and ice-sheet margins have recently undergone dramatic changes (Bunce and others, 2018; Konrad and others, 2018). These ice sheets are expected to experience further significant changes in the future (e.g. Edwards and others, 2019). These rapid dynamic changes are occurring not through the slow internal deformation of ice under the force of gravity but rather because of interactions between ice bodies and their boundaries. Whether that be ice–bed, such as hydrologically-accelerated basal sliding (e.g. Schoof, 2010); ice–ocean, such as submarine melt induced by subglacial discharge and/or fjord temperature (e.g. Jenkins, 2011; Straneo and Heimbach, 2013); or ice–air, through mass-balance feedbacks (e.g. Vizcaino and others, 2015). Ice flow dynamics are known to be very sensitive to the interface locations; for example, the stability of marine ice sheets depends fundamentally on the grounding line (GRL) position (Schoof, 2007). Consequently, careful attention must be paid to these interfaces when modelling ice-sheet flow. In particular, any choice of numerical algorithm must be guided by the need to accurately capture dynamically evolving boundaries, and hence to ensure reliable predictions of ice volume and extent and to minimize uncertainty in sea-level rise estimates.

To address these challenges, we aim to demonstrate that the level-set method (LSM) is an effective approach for accurately tracking the evolution of the ice–air and ice–water interfaces, as well as the terminus position for ice sheets and the GRL position of marine ice sheets. The level-set approach is versatile and can be incorporated into any ice-sheet model (shallow ice to full Stokes) regardless of numerical discretization of the governing equations (finite difference or finite element, fixed grid or adaptive mesh). In this paper, we present an LSM implemented in shallow ice models using a finite-difference discretization on a fixed rectangular grid.

The evolution of an ice-sheet free surface is typically modelled by solving the kinematic boundary condition

(1)$$\displaystyle{{\partial \eta } \over {\partial t}} + u\displaystyle{{\partial \eta } \over {\partial x}} + v\displaystyle{{\partial \eta } \over {\partial y}}-w = {\cal M}_j\comma \;$$

where η = h(x, y, t) for the surface elevation of the ice–air interface or η = b(x, y, t) for the elevation of the bottom ice–water interface, (u, v, w) are the ice velocity components, and ${\cal M}_j\lpar {x\comma \;y\comma \;t} \rpar$ denotes the (surface ${\cal M}_{\rm h}$ or bottom ${\cal M}_{\rm b}$) mass-balance function. We note that the geometry of an evolving ice–air or ice–water interface may experience overriding, breaking, merging, separation, discontinuities, vertical fronts and overhangs, and such events cannot be resolved in the standard setup. Furthermore, because of discretization, elevations are only evaluated at the grid points, and as such, the exact location of the terminus and GRL will fall between grid points and is not tracked explicitly. To improve accuracy, a fine-grid resolution is required, which can become prohibitively expensive. This is particularly apparent for marine ice sheets where fixed-grid methods have been shown to be inadequate in capturing GRL migration (Vieli and Payne, 2005; Pattyn and others, 2012; Seroussi and others, 2014). Consequently, ice-sheet modellers have incorporated various approaches to more accurately track the GRL. These include grid refinement near the GRL (e.g. Durand and others, 2009b; Cornford and others, 2013), requiring adaptive remeshing with each displacement of the GRL. The use of sub-element parameterizations together with mesh refinement has been shown to be beneficial (Seroussi and others, 2014), although even here the exact GRL position is not being tracked. Explicit tracking methods include coordinate stretching to transform the moving boundary onto a fixed domain using a prognostic equation for GRL migration (e.g. Hindmarsh, 1996). However, this approach cannot be generalized to 2-D because of complications in handling the complex evolving geometry of the GRL. Other alternatives include the location-based moving mesh approach of Goldberg and others (2009) where the challenge is to define a suitable monitor function to position nodes, or the velocity-based moving-point approach proposed by Bonan and others (2016) to track the ice-sheet margin, but not yet applied to addressing the GRL problem.

These difficulties addressed above can be readily handled by the use of the LSM which can capture complex evolving geometries without requiring adaptive mesh refinement, and with the further advantage of being relatively straightforward to implement. The LSM is an increasingly popular tool within computational fluid dynamics for tracking the motion of dynamic fluid interfaces and is finding widespread use in many applications (Sethian, 1999b; Gibou and others, 2018). LSMs were first applied in glaciology by Pralong and Funk (2004), who proposed the LSM with the ice–air flow problem as a means of evaluating the steady-state geometry of a glacier. Further consideration of the applicability of LSM to ice flow, compared among other numerical algorithms for free surface flows, can be found in Caboussat and others (2011). More recently, Bondzio and others (2016) used a LSM to simulate the migration of the calving front in a 2-D plan-view modelling study (see also Bondzio and others (2018)). They demonstrate the benefit of the LSM for simulating calving front dynamics in Jakobshavn Isbræ, west Greenland.

In this study, we build on the pioneering exploratory study of Pralong and Funk (2004) by providing an extensive series of simulations testing the accuracy of the LSM for tracking evolving land- and marine-terminating ice-sheet boundaries. We use flow-line and radial ice flow models together with the LSM to simultaneously model ice elevations, continental margins, GRLs and shelf fronts. We highlight implementation details not previously addressed such as the use of the signed distance function, calculation of extended velocities and the fast marching method (FMM) for reinitialization. In Pralong and Funk (2004) and Bondzio and others (2016), the LSM is implemented using the finite-element method on unstructured meshes. Instead, our implementation uses the finite-difference method on a regular fixed grid, which further highlights the strength and versatility of the LSM by demonstrating that moving ice boundaries can be tracked accurately without requiring local mesh refinement.

2. Level set method

In LSMs, the interface is represented implicitly using a level-set function φ(x, t) which is a differentiable function on a space-time domain Ω × ℝ₊, where x ∈ Ω is the spatial domain in 2-D or 3-D. The surface itself is represented as the zero isosurface or level set φ(x, t) = 0, which propagates at a speed directed normal to the surface ∂Ω.

When tracking the ice–air or ice–water interface for an ice sheet or glacier, we can define a level-set function φ with the following properties:

(2)$$\eqalign{& \varphi \lpar {{\bi x}\comma \;t} \rpar \lt 0\quad {\rm for}\quad {\bi x}\in \Omega _{\rm i}\comma \;\cr & \varphi \lpar {{\bi x}\comma \;t} \rpar \gt 0\quad {\rm for}\quad {\bi x}\in \Omega _{\rm c}\comma \;\cr & \varphi \lpar {{\bi x}\comma \;t} \rpar = 0\quad {\rm for}\quad {\bi x}\in \partial \Omega \comma \;} $$

where Ω_i represents the region inside the ice body, Ω_c is the region outside the ice body (consisting of either air or water), and ∂Ω is the ice–air or ice–water interface (see Fig. 1). Any level-set function fitting this description (Eqn 2) can be used as an initial condition for the level-set evolution equation to track an interface (see Section 2.1). However, we will use a special choice corresponding to the signed distance function (cf. Pralong and Funk, 2004) that is numerically advantageous (Vogl, 2016) and improves mass conservation, compared to other choices (Gibou and others, 2018):

(3)$$\varphi \lpar {{\bi x}\comma \;t} \rpar = \left\{{\matrix{ {-{\rm d}\lpar {{\bi x}\comma \;\partial \Omega } \rpar \comma \;} & {{\rm if}\;\;{\bi x}\in \Omega_i\comma \;} \cr {{\rm d}\lpar {{\bi x}\comma \;\partial \Omega } \rpar \comma \;} & {{\rm if}\;\;{\bi x}\in \Omega_c.} \cr } } \right.$$

The value of the level-set function corresponds to the Euclidean distance d(x, ∂Ω) between any given spatial location x and the corresponding closest point on the interface ∂Ω, with the sign chosen as negative for inside the ice and positive for outside.

Fig. 1. Basic geometry and definition of the level-set function φ(x, t) for a generic ice sheet.

2.1. Level-set evolution

The value of the level-set function at any point on an interface with location x(t) ∈ Ω must satisfy

$$\varphi \lpar {{\bi x}\lpar t \rpar \comma \;t} \rpar = 0.$$

Differentiating this equation in time and applying the chain rule, we obtain

$$\displaystyle{{\partial \varphi } \over {\partial t}} + \nabla \varphi \lpar {{\bi x}\lpar t \rpar \comma \;t} \rpar \cdot {{\bi x}}^{\prime}\lpar t \rpar = 0.$$

Supposing that F is the speed in the outward normal direction, then

$${\bi x}^{\prime}\lpar t\rpar \cdot {\bi n} = F\quad {\rm where}\quad {\bi n} = \nabla \varphi /{\rm \Vert }\nabla \varphi {\rm \Vert }.$$

Therefore, the evolution equation for the level-set function φ can be written as

(4)$$\displaystyle{{\partial \varphi } \over {\partial t}} + F\parallel \nabla \varphi \parallel\, = 0.$$

The normal speed F should depend on the ice velocity field as well as any accumulation and ablation (Pralong and Funk, 2004). To this end, let M(x, t) denote the mass-balance function and u(x, t) the ice velocity field, so that the speed function can be written as

(5)$$F = \lpar {{\bi u}\lpar {{\bi x}\comma \;t} \rpar + {\bi M}\lpar {{\bi x}\comma \;t} \rpar \hat{z}} \rpar \cdot \displaystyle{{\nabla \varphi } \over {\parallel \nabla \varphi \parallel }}\comma \;$$

where $\hat{z}$ is the unit vector in the vertical direction and we have assumed vertical accumulation and ablation (Pralong and Funk, 2004). M(x, t) is determined by the surface and basal mass-balance functions, ${\cal M}_{\rm h}$ and ${\cal M}_{\rm b}$, respectively. M(x, t) and u(x, t) can be derived from the specific ice-sheet model or fit to observational data. The level-set equation is solved throughout the computational domain which requires speed function values, F, at all grid points both inside and outside the ice. This necessitates the computing of so-called extended speed, the full details of which are provided in Appendix A1. The numerical implementation of the level-set equation is described in Appendix A2 and we include an explanation of a re-initialization procedure that uses the FMM in Appendix A3.

2.2. Coupling between ice-sheet model and LSM

For the computational experiments performed in this study, the LSM described above requires a 2-D ice velocity field u = (u, w) that comes from an ice flow model, which we obtain using either the shallow ice approximation (SIA) or the shallow shelf approximation (SSA), see Appendix B1 and B2, respectively. The ice velocities u and w, the speed function F and the level-set function φ are defined on a 2-D grid, with the ice thickness H and the height of the upper surface h defined on a 1-D grid. The overall solution procedure and details regarding the coupling between ice-sheet model and the LSM are as follows:

1. (1) Computing the velocities, u ^k and w ^k: At any time t ^k, the ice thickness H ^k and the height of the upper surface h ^k are known. H ^k, h ^k and the horizontal derivatives of h ^k (determined using the central difference approximation) are used to calculate the horizontal and vertical velocities at t ^k, u ^k and w ^k, respectively. In the radially symmetric SIA experiments, u ^k is computed using Eqn (B3) and w ^k using Eqn (B4). Whereas for the marine ice sheet and ice-shelf experiments, we apply the Picard iteration of the SSA, Eqn (B11), to compute u ^k. The derivatives of u ^k and the bedrock are then determined using the central difference approximation and w ^k is computed using Eqn (B9), for ice stream, and Eqn (B10), for ice shelf.

2. (2) Computing the speed, F ^k: The velocities, u ^k and w ^k, together with the given surface mass balance function, are used to determine F ^k for grid points interior the ice sheet using Eqn (5). For grid points exterior to the ice sheet F ^k is determined using the procedure described in Appendix A1.

3. (3) Computing φ^k+1: The discretized level-set evolution Eqn (A5) and the numerical procedure described in Appendix A2 are used to determine the level set at time t ^k+1, φ^k+1.

4. (4) Extract H ^k+1: The zero level set is determined with subgrid scale precision from the zero contour line of φ^k+1. The zero level set identifies the ice surface and is used to determine H ^k+1 and h ^k+1. The terminus or GRL position is identified as the location where the zero level set intersects the bottom domain boundary. Although we note that our ice velocity solver does not use any subgrid method; for example, the exact GRL position is not used in the SSA computation.

5. (5) Re-initialization: The FMM algorithm (Appendix A3) can be used to rebuild the level-set function φ. This step is only necessary every 50–100 time steps depending on the complexity of the problem and the required accuracy.

3. Numerical results

Our aim in this section is to validate the LSM as an effective method for capturing ice-sheet evolution, including GRL migration, by comparing results from numerical simulations against various analytical and benchmark solutions. We first compare the LSM result for an idealized glacier test case with a prescribed velocity field and mass balance, after which we examine the behaviour and performance of the LSM for tracking both grounded and marine ice-sheet boundaries.

3.1. An idealized test case

Following Pralong and Funk (2004), we first consider an idealized glacier test case in order to focus in on the coupling between the level-set calculation and the ice flow problem. This test case fixes the ice flow solution with the given velocity field u(x, z) = x ² + z ² and w(x, z) = 0 and glacier surface height h(x, t) = x − x ² + xt. The corresponding extended mass-balance function is chosen as ${\cal M}\lpar {x\comma \;z\comma \;t} \rpar = x + \lpar {x^2 + z^2} \rpar \lpar {1-2x + t} \rpar$, so that Eqn (1) is satisfied identically in 2-D (Picasso and others, 2004). The evolution of this ‘glacier surface’ is then simulated using the level-set Eqn (4) with the imposed flow field u, w and mass-balance function ${\cal M}\lpar {x\comma \;z\comma \;t} \rpar$ so that the evolved surface height may be compared with the analytical solution h(x, t).

The numerical results, from the initial position (at t = 0) to the final position (at t = 2), show excellent agreement with the analytical solution (Fig. 2a). We calculate the discrete ℓ₁-norm of the absolute error for h at t = 2 on uniform spaced grids of 60 × 60, 75 × 75, 90 × 90 and 105 × 105 using a time step Δt = 0.005 without re-initialization (Fig. 2b). From these results, we estimate the rate of convergence of the ℓ₁ error norm to be ${\cal O}\lpar {\Delta x^{1.3}} \rpar$ (and ${\cal O}\lpar {\Delta x^{1.9}} \rpar$ for the ℓ₂ error).

Fig. 2. (a) Evolution of the interface due to an imposed velocity field and surface mass balance. The interface position is shown at equally-spaced times between t = 0 and 2. The points (‘°’) represent the analytical solution and the solid lines correspond to the numerical approximation of the LSM. The simulation uses a spatial grid 60 × 60 and time step Δt = 0.005. (b) Verification of the rate of convergence ${\cal O}\lpar {\Delta x^{1.3}} \rpar$ using the ℓ₁-norm error.

3.2. Halfar similarity solution

We next apply the LSM to a well-known time-dependent solution of the SIA model called the Halfar similarity solution. To find this exact solution, we first write the no slip SIA equation (Schoof and Hewitt, 2013)

(6)$$\displaystyle{{\partial H} \over {\partial t}} = \nabla .\lpar {\Gamma H^{n + 2}\vert \nabla h\vert^{n-1}\nabla h} \rpar + {\cal M}\lpar {{\bi x}\comma \;t} \rpar \comma \;$$

where Γ = (2A/(n + 2))(ρg)ⁿ and H denotes the ice thickness (refer to Appendix B1 for details). Halfar (1981) derived a similarity solution for this problem in the case of a flat bed (b(x) = 0) and no surface mass balance $\lpar {\cal M}\lpar {\bi x}\comma \;t\rpar = 0\rpar$. Supposing that H(0, t ₀) = H ₀ for t ₀ > 0 and the distance from the origin r = (x ² + y ²)^1/2, the 2-D Halfar solution to the SIA is

(7)$$H\lpar {r\comma \;t} \rpar = H_0\left({\displaystyle{{t_0} \over t}} \right)^\alpha \left[{1-{\left({{\left({\displaystyle{{t_0} \over t}} \right)}^\beta \displaystyle{r \over {R_0}}} \right)}^{(n + 1)/n}} \right]^{n/(2n + 1)}\comma \;$$

where $t_0 = \lpar {\beta /\Gamma } \rpar \lpar {7/4} \rpar ^3R_0^4 \;H_0^{{-}7}$ is a characteristic time (Bueler and others, 2005). Note that the values of the parameters α = 1/9 and β = 1/18 are such that the factors t ^−α and t ^−β change very slowly for large times t. Other parameters used in this computation are the Glen's flow law exponent n = 3 and ice softness A = 10⁻¹⁶ Pa⁻³ a⁻¹, ice density ρ = 910 kg m⁻³, gravitational acceleration g = 9.81 m s⁻², H ₀ = 3600 m and R ₀ = 750 km.

We choose initial time t = 100 in Eqn (7) so that H(r, t = 100) is the initial ice thickness. To compute the surface elevation of the Halfar similarity solution using the LSM, we evolve the level-set Eqn (4) with the SIA horizontal and vertical velocities (Eqns (B3) and (B4), respectively) to compute the new surface at time t, identified with φ(r, z, t) = 0. The level-set function is computed on the domain [0,1000] km × [0,5000] m with 200 × 100 grid size (Δr = 5 km along the radial axis and Δz = 50 m along the vertical axis).

We compute an extended velocity field outside the ice, which is then used to generate a signed distance function (Eqn (3)) depicted in Figure 3a based on the initial level-set function φ(r, z, t = 100). This initial surface is then evolved in time and Figure 3b compares the surface elevation of the ice sheet at t = 1000 and 10000 years (where φ(r, z, 1000) = 0 and φ(r, z, 10000) = 0) with the Halfar solution (Eqn (7)). The absolute error between the elevations from the exact Halfar solution and the LSM approximation at t = 10000 years is small, limited to at most 30 m at the ice divide (Fig. 3c). The error increases towards the margin due to steeper surface gradients near the terminus (Fig. 3c). In spite of this, the position of the margin is obtained to within a relative error of 0.29% (comparing 894.14 km (Halfar) with 896.71 km (LSM)). These results are in close agreement with the exact solution. We repeated the experiment with uniformly spaced grids of 160 × 80, 200 × 100, 250 × 125 and 300 × 150 and found the rate of convergence of the absolute error in the position of the margin to be ${\cal O}\lpar {\Delta r^{1.1}} \rpar$.

Fig. 3. (a) A contour plot of the initial level-set function φ at t = 100 with the red line denoting the zero level set, φ = 0, of the ice surface; (b) surface elevations of the evolving ice sheet; (c) the absolute error between the Halfar exact solution and the computed LSM solution at t = 10000 years.

3.3. Radially symmetric ice-sheet experiments

In this section, we perform a moving-margin experiment described in the European Ice Sheet Modelling INiTiative (EISMINT) intercomparison project (Huybrechts and Payne, 1996). The aim of this experiment is to find a steady-state ice-sheet surface solution for a given mass-balance function.

Before addressing the EISMINT experiment, we first introduce the steady-state solution for an ice sheet with flat bedrock which occurs when the net mass of ice remains constant over some period of time. In other words, the rate of change of the ice thickness (∂H/∂t) = 0, so that the SIA mass balance (Eqn (6)) reduces at steady state to

(8)$${\cal M}\lpar r \rpar -\displaystyle{1 \over r}\displaystyle{{\partial \lpar {rH\bar{u}} \rpar } \over {\partial r}} = 0\comma \;$$

where r is the radial coordinate and $\bar{u}$ is the vertically averaged ice velocity. Assuming that the surface mass balance ${\cal M}\lpar r \rpar$ is independent of time and the bedrock is flat, the solution for the steady ice thickness profile is

(9)$$\eqalign{H\lpar r \rpar = \,& \left({{\left({\displaystyle{{2\lpar {n + 1} \rpar } \over {n\rho g}}} \right)}^n\displaystyle{{n + 2} \over {2A}}} \right)^{1/2\lpar {n + 1} \rpar } \cr & \times \left({\mathop \int \limits_r^R {\left({\displaystyle{1 \over \xi }\int_0^\xi {{\cal M}\lpar \eta \rpar \eta \;d\eta } } \right)}^{1/n}\;d\xi } \right)^{n/2\lpar {n + 1} \rpar }\comma \;}$$

where R denotes the margin position at steady state (Bonan and others, 2016). The EISMINT intercomparison project imposes the surface mass-balance function

(10)$${\cal M}\lpar r\rpar = \min \lpar {0.5, {10}^{{-}2}\cdot \lpar 450-r\rpar } \rpar \;{\rm m}{\mkern 1mu} {\rm a}^{ - 1}\comma \;$$

for which the bisection method (a simple numerical root-finding algorithm) can be applied to $\int _0^R {\cal M}\lpar \eta \rpar \eta \;d\eta = 0$ to find the steady-state margin position of R = 579.81 km. We note that the integral $\int _0^\xi {\cal M}\lpar \eta \rpar \eta \;d\eta$ in Eqn (9) can be evaluated numerically using (Eqn (10)), and so the steady H profile can be estimated accurately using Simpsons 1/3 rule (see Bonan and others, 2016). Parameters used in the EISMINT benchmarks are the Glen's flow law exponent n = 3 and ice softness A = 10⁻¹⁶ Pa⁻³ a⁻¹, ice density ρ = 910 kg m⁻³ and gravitational acceleration g = 9.81 m s⁻².

We perform two experiments on a flat bedrock with no sliding. The first is the EISMINT moving margin experiment, designed with no initial ice mass (Huybrechts and Payne, 1996). The second is similar but initialized with the following ice mass profile (Bonan and others, 2016):

(11)$$H_0\lpar r\rpar = 1000\left({1-{\left({\displaystyle{r \over {450}}} \right)}^2} \right)\;{\rm m}.$$

Figures 4a, b depict the evolution of the ice-sheet geometry for these two initial conditions and both are run for 20 000 years to reach the steady state. The steady-state ice divide thickness was found to be 2986.91 and 2987.81 m for the two experiments, which both lie within the range 2982.3 ± 26.4 m given by the EISMINT intercomparison, and are extremely close to the numerically integrated value obtained from Eqn (9) (Table 1). Similarly, the margin position of both experiments is very close to the numerically integrated reference value in Table 1. The absolute error between the simulated result and the reference ice thickness across the profile is mostly <1 m and never rises above 4.1 m (Fig. 4c). The relative ℓ₁-norm errors of the surface elevation with and without initial ice mass experiments are 0.038 and 0.036%, respectively. After repeating the experiment, without initial ice mass, for grid sizes 192 × 60, 216 × 60, 240 × 60 and 264 × 60, we determined the rate of convergence of the absolute error in the position of the margin to be ${\cal O}\lpar {\Delta r^{1.5}} \rpar$. Both Table 1 and Figure 4 show that our LSM method is able to achieve an excellent equivalent estimation of the EISMINT intercomparison result without using coordinate stretching or grid refinement near the terminus.

Fig. 4. Ice surface solutions for the EISMINT moving-margin experiment with (a) zero initial ice mass and (b) initial ice mass given by Eqn (11). The LSM simulated profiles are shown every 1000 years (blue lines) until steady state at t = 20000 years (red line) and the steady-state reference solution is represented by circles. (c) The absolute error at the steady state between the LSM (without and with an initial ice mass) and the numerical reference value.

Table 1. Steady-state results from the EISMINT moving-margin experiment. A comparison between the benchmark solutions (see Table 5 in Huybrechts and Payne, 1996), the reference solution from numerical integration using Eqn (9), and the LSM solutions with grid size 240 × 60 (Δr = 2.7 km and Δz = 60 m) obtained without an initial ice mass and with an initial ice mass from Eqn (11)

We now investigate two further experiments, following Bonan and others (2016), which use the EISMINT surface mass balance with no initial ice mass and a non-flat (fixed) bedrock elevation

(12)$$b\lpar r\rpar = 2000-2000\left({\displaystyle{r \over {300}}} \right)^2 + 1000\left({\displaystyle{r \over {300}}} \right)^4-150\left({\displaystyle{r \over {300}}} \right)^6\;{\rm m}.$$

The first experiment considers no sliding, and the second includes sliding with a bed friction parameter C = 7.624 × 10⁶ Pa s^1/3 m^−1/3. Ice-sheet profiles for these non-flat bedrock experiments reach a steady state by roughly t = 10000 years (see Fig. 5). The ice divide thickness at steady state was found to be 4026.25 m for the non-sliding case and 3801.72 m for the sliding case. The margin position of the non-sliding and sliding experiments is 571.81 and 578.47 km, respectively. We have found that the LSM produces smooth changes along the ice interface in contrast to the moving-point approach of Bonan and others (2016) which show linear gradients of the ice-sheet surface near the margin as a result of the mesh size (refer to Figs 2, 4 from Bonan and others (2016)). This is because our ice interface is determined with sub-gridscale accuracy as we interpolate between the 2-D level-set values to determine the zero level set, whereas the ice interface in Bonan and others (2016) is linearly interpolated between the 1-D grid points of the ice thickness evolution equation.

Fig. 5. Ice surface solutions for the EISMINT moving-margin experiment with non-flat bedrock (a) without basal sliding and (b) with basal sliding. The LSM simulated profiles are shown every 1000 years (blue lines) until steady state at t = 10000 years (red line). The grid size is 240 × 60 (Δr = 2.5 km and Δz = 50 m).

3.4. Marine ice-sheet experiments

For the remainder, we shift our focus to simulating marine ice sheets using the governing equations for the SSA described in Appendix B2.

3.4.1. Steady ice shelf

The velocity and thickness of a steady 1-D marine ice shelf can be computed analytically due to the relative simplicity of the SSA Eqn (B7) (Van der Veen, 1983). The exact solution depends on the velocity and ice thickness at the GRL, which we take to be u _g = 50 m a⁻¹ and H _g = 500 m, respectively. We propose two experiments, one where the surface mass balance ${\cal M} = 0$ and a second with ${\cal M} = 0.3\;{\rm m}\;{\rm a}^{{-}1}$. In both cases, the initial condition is an ice-shelf block having dimensions [0, 50] km × [(ρ/ρ _w)H _g, ((ρ/ρ _w) − 1)H _g] m, which satisfies the flotation criteria. The bottom boundary is considered a free surface in the water without any basal melting. The horizontal velocity is computed using the Picard iteration (B11) and then the vertical velocity is determined using Eqn (B10). The grid spacings for the LSM solver are Δx = 0.5 km and Δz = 5 m with the same grid spacing Δx _v = 0.5 km used for the velocity solver. The ice surface is evolved in time using the LSM (Eqn (4)) and Figures 6a, b depict the surface profiles at various times for the two experiments with and without accumulation. The steady-state ice thickness at x = 50 km for the case without accumulation is found to be 223.14 m and with accumulation $\lpar {\cal M} = 0.3\;{\rm m}\;{\rm a}^{{-}1}\rpar$ is 290.73 m. The absolute error is slightly higher near the GRL as expected due to the steeper ice thickness gradient at this location (Fig. 6c). The relative ℓ₁-norm errors in the ice–air shelf interface with and without accumulation are 0.24 and 0.34%, respectively, and the rates of convergence using the ℓ₁-norm error are ${\cal O}\lpar {\Delta x^{1.1}} \rpar$ and ${\cal O}\lpar {\Delta x^{1.3}} \rpar$, respectively. These numerical results show very good agreement with the exact steady-state solution for these ice-shelf test cases.

Fig. 6. Evolution of the ice shelf interface using the shallow shelf approximation for cases (a) zero accumulation and (b) accumulation${\cal M} = 0.3\;{\rm m}\;{\rm a}^{{-}1}$. The initial shelf is a rectangular block of ice and the interface is displayed every 50 years, with the steady state highlighted in red. The points (‘°’) show the exact ice shelf solution for comparison. (c) Absolute error of the steady state (t = 1000 years) for both experiments.

3.4.2. Marine ice-sheet benchmark experiment

Our final experiment is a study of the full marine ice sheet that includes a grounded ice stream attached to a floating ice shelf. The goal here is to examine the ability of the LSM to accurately track GRL migration. We will study the hysteresis effect for a 2-D symmetrical marine ice sheet and compare results with the benchmark Experiment 3 (EXP 3) from the Marine Ice Sheet Intercomparison Project (MISMIP) (Pattyn and others, 2012). The setup uses an overdeepened bed with polynomial shape

(13)$$\eqalign{b\lpar x\rpar = 729-2184.8\left({\displaystyle{x \over {750}}} \right)^2 + 1031.72\left({\displaystyle{x \over {750}}} \right)^4-151.72\left({\displaystyle{x \over {750}}} \right)^6{\rm m\comma \;}}$$

and model parameter values given in Table 2. The experiment consists of a sequence of 13 steps (or time intervals) of a given length, in each of which the ice sheet has a different value of the Glen's flow rate constant A (the data are summarized in Table 3). Horizontal and vertical grid spacings for the level-set discretization are Δx = 7.5 km and Δz = 60 m. To ensure a sufficiently accurate horizontal velocity, for this experiment, we use a finer grid for the velocity solver with spacing Δx _v = 1.875 km, although only velocities coincident with the level-set grid are used by the level-set solver. We use a time step of Δt = 5 years and every 500 years the level-set function is re-computed using the FMM (see Appendix A3). The solution is initialized with a 50 m thick grounded ice layer that extends to the location where it becomes afloat at position x = 479.1 km.

Table 2. Parameter values of the marine ice-sheet experiments

Table 3. Values of the Glen's flow law rate constant A and time intervals used for each step of the MISMIP EXP 3 benchmark, corresponding to the simulations displayed in Figure 7 (Pattyn and others, 2012)

The model simulation proceeds as follows. Starting from the initial values, the solution is computed using the value of A listed for step 1 in Table 3 and over the corresponding time interval. The code is then restarted with the new value of A listed in step 2 and using the result from the first step as the initial state. This procedure continues until the end of the final step and the results are shown in Figure 7.

Fig. 7. Simulated steady-state profiles of the MISMIP EXP 3 results. Steps 1–13 correspond to the parameter changes listed in Table 3.

The corresponding GRL position is plotted as a function of 1/A in Figure 8 in order to demonstrate the hysteresis phenomenon and compare with the MISMIP results in Figure 5 of Pattyn and others (2012). In Figure 8, the black S-shaped curve represents the path according to the boundary layer theory of Schoof (2007) with our modelled steady-state GRL positions correctly located on the upper and lower branches of this approximate analytic solution. The SSA fixed grid models used in MISMIP EXP 3 (EBU1 (Δx = 6 km), DPO4 (Δx = 0.1 km) and FPA5 (Δx = 0.3 km)) are either unable to reproduce hysteresis or reproduce it qualitatively with solutions lying several tens of kilometres from the boundary layer theory (see red plots in Fig. 8). In contrast, our lower resolution fixed-grid (Δx = 7.5 km) level-set approach produces modelled positions that closely match those only achieved with the highly resolved (Δx ≤ 1.2 km) moving grid methods or the finest adaptive grid (Δx = 0.15 km) SSA models in the MISMIP (see Fig. 5 and Table 2 in Pattyn and others (2012)). Our results in Figure 8 are also plotted with SSA-H FPA4 (Δx = 1.2 km), a MISMIP participating SSA model with the Pollard and DeConto heuristic (a GRL parameterization that uses the matched asymptotics of Schoof (2007)), and show excellent agreement (Pattyn and others, 2012).

Fig. 8. Hysteresis in the grounding line position as a function of forcing viscosity (A ⁻¹) for MISMIP EXP 3. The black line is from the boundary layer theory of Schoof (2007); ‘○’ points represent results from our LSM simulations; the red points ‘●’, ‘×’ and ‘◇’ depict results from the fixed grid MISMIP participating models SSA FPA5, SSA EBU1 and SSA DPO4, respectively; and ‘□’ points are from the MISMIP participating model SSH-H FPA4 which uses the Pollard and DeConto heuristic (see Fig. 5 in Pattyn and others (2012)).

Another test of the numerical results is to compare the difference in GRL position between steps 2 and 12, where any differences are expected to be a result of numerical approximation (Durand and others, 2009a). For our LSM configuration, we calculate this gap to be 3.17 km (2.13 km between steps 1 and 13) using our relatively coarse uniform grid (Δx = 7.5 km for the LSM and Δx _v = 1.875 km for the velocity solver). We note that in Figure 6 of Durand and others (2009a), the grid size has to be as low as 40 m to achieve a similar degree of accuracy and a mesh size of 25 − 200 m close to the GRL was used by Gagliardini and others (2016); however, these are full-Stokes models and so use a different physical approximation needing higher resolution.

To investigate further, we examine the sensitivity of the GRL position (x _g) of the steady-state ice-sheet profiles between steps 2 and 12 for different mesh resolution. For this reason, we compare results using a grid size of (Δx, Δz) = (10 km, 80 m), (7.5 km, 60 m), (5 km, 40 m) and (3.75 km, 30 m) for the LSM with a fixed grid size Δx _v = 15/8 km for the velocity solver. Similarly, we also consider different velocity grid sizes with Δx _v = 15/8, 15/16, 15/32 and 15/64 km and a fixed grid size (Δx, Δz) = (3.75 km, 30 m) for the LSM. Results are presented in Figures 9a, b where steady x _g are plotted as a function of the horizontal grid size Δx of the LSM and the grid size Δx _v of the velocity solver, respectively. In Figure 9a, the GRL position gap between steps 2 and 12 reduces as the grid size is reduced and reaches 1.89 km at the LSM grid Δx = 3.75 km. Numerical results also depend on the accuracy of velocities that we use for the LSM. When steady x _g are presented as a function of the grid size of the velocity solver (for a fixed LSM grid size Δx = 1.875 km) the gap reduces to 150 m and we find that the GRL converges to a value near 730 km at the lowest resolution (Fig. 9b).

Fig. 9. Evolution of the steady x _g as a function of (a) the horizontal mesh size Δx of the LSM for fixed mesh size Δx _v = 1.875 km of the velocity solver and (b) the mesh size Δx _v of the velocity solver for fixed mesh size Δx = 3.75 km of the LSM. Blue circles (red squares) represent results obtained for simulations starting from the steady state obtained at step 2 (step 12). The dashed line depicts results obtained using Schoof's boundary layer (BL) theory reported in Durand and others (2009a).

We have shown that by using the LSM on a relatively coarse fixed grid, we can determine the evolving GRL position with fine-scale accuracy. The zero level set determines the ice interface and is found by interpolating level-set values computed on the fixed grid. The GRL position is identified as the location where the zero level set meets the bottom of the domain. In spite of our coarse grid, this method allows us to determine the GRL position with subgrid scale precision. The LSM requires an additional dimension and therefore greater computational time is needed compared to solving the kinematic boundary condition (Eqn (1)). This additional computational time cost is especially apparent when solving the relatively fast SIA and SSA, but would not be so significant if the LSM were coupled to the more computationally demanding full-Stokes equations. Regardless, simulating GRL migration has required ice-sheet models using irregular and adaptive mesh refinements, which come with considerable computational cost and complexity, whereas we have shown that a regular fixed-grid model, using the LSM, can accurately track advancing and retreating GRL positions. Furthermore, the coarser grid also results in computational savings from the longer time steps allowed by the CFL condition. The other alternative is to use a subgrid parameterization, such as using a heuristic rule based on boundary layer theory valid for steady states. The level-set approach used here does not rely on a parameterization or employ any other special treatment at the GRL. The method described has not been tried for two horizontal dimensions, but the 3D LSM could be implemented by extending the array structures and gradient operators for tracking the propagating surfaces in two horizontal dimensions (Sethian, 1999b; Bondzio and others, 2016).

4. Conclusions

We have devised a new level-set algorithm for tracking an evolving ice-sheet surface and GRL position, based on an underlying fixed-grid finite-difference discretization. Other fixed-grid methods tend to be less competitive relative to moving grid methods for dynamic interface problems like ice-sheet models, and they can only obtain comparable accuracy with moving-grid methods by using highly resolved grids near the GRL. Our level-set approach is able to track the ice-sheet margin and GRL location dynamically for both grounded and marine ice sheets without the need for grid refinement and any subgrid parameterization or heuristic rule. The method is tested by comparing numerical simulations with analytical and benchmark solutions. In particular, we compared model solutions for grounded ice sheets with the EISMINT benchmark (Huybrechts and Payne, 1996) and for marine ice sheets with the MISMIP intercomparison benchmark (Pattyn and others, 2012). These experiments demonstrate that the LSM is an accurate approach for capturing the ice-sheet marginal position, while exploiting the efficiency in using an underlying fixed grid that is coarse relative to other methods that employ a uniformly or locally adapted fine mesh. We have therefore shown that the LSM is an accurate and efficient approach for tracking the ice surface interface, terminus positions and GRLs for grounded and marine ice sheets.