1. Introduction
The Greenland and Antarctic ice sheets have been losing mass at an accelerated rate (Bevis and others, Reference Bevis2019; Rignot and others, Reference Rignot2019) and ice-sheet margins have recently undergone dramatic changes (Bunce and others, Reference Bunce, Carr, Nienow and Ross2018; Konrad and others, Reference Konrad2018). These ice sheets are expected to experience further significant changes in the future (e.g. Edwards and others, Reference Edwards2019). 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, Reference Schoof2010); ice–ocean, such as submarine melt induced by subglacial discharge and/or fjord temperature (e.g. Jenkins, Reference Jenkins2011; Straneo and Heimbach, Reference Straneo and Heimbach2013); or ice–air, through mass-balance feedbacks (e.g. Vizcaino and others, Reference Vizcaino2015). 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, Reference Schoof2007). 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
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, Reference Vieli and Payne2005; Pattyn and others, Reference Pattyn2012; Seroussi and others, Reference Seroussi, Morlighem, Larour, Rignot and Khazendar2014). 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, Reference Durand, Gagliardini, Zwinger, Le Meur and Hindmarsh2009b; Cornford and others, Reference Cornford2013), 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, Reference Seroussi, Morlighem, Larour, Rignot and Khazendar2014), 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, Reference Hindmarsh1996). 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 (Reference Goldberg, Holland and Schoof2009) 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 (Reference Bonan, Baines, Nichols and Partridge2016) 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, Reference Sethian1999b; Gibou and others, Reference Gibou, Fedkiw and Osher2018). LSMs were first applied in glaciology by Pralong and Funk (Reference Pralong and Funk2004), 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 (Reference Caboussat, Jouvet, Picasso, Rappaz and Magoules2011). More recently, Bondzio and others (Reference Bondzio2016) used a LSM to simulate the migration of the calving front in a 2-D plan-view modelling study (see also Bondzio and others (Reference Bondzio, Morlighem, Seroussi, Wood and Mouginot2018)). 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 (Reference Pralong and Funk2004) 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 (Reference Pralong and Funk2004) and Bondzio and others (Reference Bondzio2016), 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:
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, Reference Pralong and Funk2004) that is numerically advantageous (Vogl, Reference Vogl2016) and improves mass conservation, compared to other choices (Gibou and others, Reference Gibou, Fedkiw and Osher2018):
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
Differentiating this equation in time and applying the chain rule, we obtain
Supposing that F is the speed in the outward normal direction, then
Therefore, the evolution equation for the level-set function φ can be written as
The normal speed F should depend on the ice velocity field as well as any accumulation and ablation (Pralong and Funk, Reference Pralong and Funk2004). 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
where $\hat{z}$
is the unit vector in the vertical direction and we have assumed vertical accumulation and ablation (Pralong and Funk, Reference Pralong and Funk2004). 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) 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) 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) 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) 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) 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 (Reference Pralong and Funk2004), 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 2 + z 2 and w(x, z) = 0 and glacier surface height h(x, t) = x − x 2 + 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, Reference Picasso, Rappaz, Reist, Funk and Blatter2004). 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 ℓ1-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 ℓ1 error norm to be ${\cal O}\lpar {\Delta x^{1.3}} \rpar$
(and ${\cal O}\lpar {\Delta x^{1.9}} \rpar$
for the ℓ2 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 ℓ1-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, Reference Schoof and Hewitt2013)
where Γ = (2A/(n + 2))(ρg)n and H denotes the ice thickness (refer to Appendix B1 for details). Halfar (Reference Halfar1981) 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 0) = H 0 for t 0 > 0 and the distance from the origin r = (x 2 + y 2)1/2, the 2-D Halfar solution to the SIA is
where $t_0 = \lpar {\beta /\Gamma } \rpar \lpar {7/4} \rpar ^3R_0^4 \;H_0^{{-}7}$
is a characteristic time (Bueler and others, Reference Bueler, Lingle, Kallen-Brown, Covey and Bowman2005). 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−16 Pa−3 a−1, ice density ρ = 910 kg m−3, gravitational acceleration g = 9.81 m s−2, H 0 = 3600 m and R 0 = 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, Reference Huybrechts and Payne1996). 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
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
where R denotes the margin position at steady state (Bonan and others, Reference Bonan, Baines, Nichols and Partridge2016). The EISMINT intercomparison project imposes the surface mass-balance function
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, Reference Bonan, Baines, Nichols and Partridge2016). Parameters used in the EISMINT benchmarks are the Glen's flow law exponent n = 3 and ice softness A = 10−16 Pa−3 a−1, ice density ρ = 910 kg m−3 and gravitational acceleration g = 9.81 m s−2.
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, Reference Huybrechts and Payne1996). The second is similar but initialized with the following ice mass profile (Bonan and others, Reference Bonan, Baines, Nichols and Partridge2016):
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 ℓ1-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, Reference Huybrechts and Payne1996), 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 (Reference Bonan, Baines, Nichols and Partridge2016), which use the EISMINT surface mass balance with no initial ice mass and a non-flat (fixed) bedrock elevation
The first experiment considers no sliding, and the second includes sliding with a bed friction parameter C = 7.624 × 106 Pa s1/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 (Reference Bonan, Baines, Nichols and Partridge2016) 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 (Reference Bonan, Baines, Nichols and Partridge2016)). 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 (Reference Bonan, Baines, Nichols and Partridge2016) 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, Reference Van der Veen1983). The exact solution depends on the velocity and ice thickness at the GRL, which we take to be u g = 50 m a−1 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 ℓ1-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 ℓ1-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, Reference Pattyn2012). The setup uses an overdeepened bed with polynomial shape
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, Reference Pattyn2012)
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 (Reference Pattyn2012). In Figure 8, the black S-shaped curve represents the path according to the boundary layer theory of Schoof (Reference Schoof2007) 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 (Reference Pattyn2012)). 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 (Reference Schoof2007)), and show excellent agreement (Pattyn and others, Reference Pattyn2012).
Fig. 8. Hysteresis in the grounding line position as a function of forcing viscosity (A −1) for MISMIP EXP 3. The black line is from the boundary layer theory of Schoof (Reference Schoof2007); ‘○’ 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 (Reference Pattyn2012)).
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, Reference Durand, Gagliardini, de Fleurian, Zwinger and Le Meur2009a). 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 (Reference Durand, Gagliardini, de Fleurian, Zwinger and Le Meur2009a), 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 (Reference Gagliardini2016); 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 (Reference Durand, Gagliardini, de Fleurian, Zwinger and Le Meur2009a).
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, Reference Sethian1999b; Bondzio and others, Reference Bondzio2016).
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, Reference Huybrechts and Payne1996) and for marine ice sheets with the MISMIP intercomparison benchmark (Pattyn and others, Reference Pattyn2012). 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.
Acknowledgements
We gratefully appreciate the clear and detailed feedback from the anonymous reviewers and the Scientific Editor Ian Hewitt. This work has been supported by a Natural Sciences and Engineering Research Council (NSERC) Discovery Grant.
Appendix A. Numerical implementation of the LSM
A.1. Extended speed
The level-set Eqn (4) requires that the speed F (Eqn (5)) is defined for all level sets throughout the computational domain Ω, not just the zero level set or on one side of the interface. Firstly, the mass-balance source term is prescribed on the interface only and must be smoothly extended within the ice for numerical stability. As such M is determined by linearly interpolating vertically between ${\cal M}_{\rm h}$
on the ice–air interface and ${\cal M}_{\rm b}$
on the ice–water (or ice–bedrock) interface. Note that this is an artificial measure that helps to smooth the derivatives in the level-set method solver and not a real change in ice mass. In our case, the ice velocity components are obtained from an ice-sheet model and thus defined inside the ice region (Ωi). Hence, these velocities, added to M, must now be extended outside the ice domain (Ωc) (Adalsteinsson and Sethian, Reference Adalsteinsson and Sethian1999). Given a level-set function φ, our goal is then to construct the extended speed F ext such that
where we require that F ext matches the given speed F on the zero level set,
This new speed field F ext is known as the ‘extended speed’ (see Fig. 10).
Fig. 10. Constructing extended speeds. The solid line inside the domain represents the ice–air interface or zero level set. Suppose F is known at ‘○’ points inside the ice then F ext must be extended to ‘*’ points outside the ice.
A desirable feature of F ext is that it should move the neighbouring level sets in such a way that the signed distance function is preserved. Following Zhao and others (Reference Zhao, Chan, Merriman and Osher1996), φ(x(t), t) remains a signed distance function if and only if
which in 2-D becomes
Since the interface may experience topological changes, a few different cases must be considered when determining F ext (Sethian, Reference Sethian1999b). As an example, suppose (i − 1, j) and (i, j − 1) are the points where F is known (inside the ice), then using finite differences we can approximate the extended speed F ext at the position (i, j) (outside the ice) by
where Δx and Δz are the horizontal and vertical grid spacings, respectively. This approach results in Eqn (A2) being satisfied for all points outside the ice.
A.2. Numerical scheme
The LSM is a versatile numerical technique that can be implemented in concert with a variety of discretizations including finite differences, finite elements, moving meshes, etc. For the sake of simplicity, we have chosen to use a fixed, rectangular, Euclidean mesh in which all grid cells are of equal size although the grid spacing in each direction may be different. After defining discrete values of φ and F at every gridpoint in the computational domain, we use a discrete form of the governing equations to evolve φ forward in time, and hence transport the interface across the underlying grid.
We use an explicit Runge–Kutta (RK) type scheme to determine φ(x, t + Δt) based on known previous values of φ(x, t), the speed in the outward normal direction F, and the gradient $\nabla \varphi$
. For a given time, t k, let φk = φ(t k) and after some time increment Δt, we denote new values as φk+1 = φ(t k + Δt). We implement the Total Variation Diminishing Runge–Kutta (TVD-RK) scheme of second order, also known as Heun's method or the modified Euler method (Osher and Fedkiw, Reference Osher and Fedkiw2006)
where
and ${\parallel} \nabla \varphi _{ij}\parallel{ = } \sqrt {{\lpar {{\lpar {\varphi_x} \rpar }_{ij}} \rpar }^2 + {\lpar {{\lpar {\varphi_z} \rpar }_{ij}} \rpar }^2}$
. Here (φx)ij and (φz)ij denote the spatial derivatives of φ at the position (x i, z j). As is usual for explicit time-stepping schemes, the allowable time step Δt is restricted in practice by a Courant–Friedrichs–Lewy (CFL) condition that depends on the spatial grid size Δx and the flow speed.
Moving on to the spatial discretization, traditional finite-difference methods based on fixed stencil interpolations work well for globally smooth problems, but at second or higher order spatial accuracy, these schemes are necessarily oscillatory near a discontinuity. We therefore approximate spatial derivatives (φx)ij and (φz)ij in Eqn (A5) using the Essentially Non-Oscillatory (ENO) scheme (Osher and Fedkiw, Reference Osher and Fedkiw2006). In this approach, a higher order non-oscillatory interpolant for piecewise smooth functions is used to approximate φ, and is then differentiated piecewise to obtain a corresponding discrete approximation for $\nabla \varphi$
. In essence, the ENO approach extends first-order accurate upwind differencing to second-order spatial accuracy in a way that suppresses oscillations.
We next discuss suitable choices for initial and boundary values of φ. The level-set function can be initialized simply using Eqn (3) with ∂Ω as the initial surface for the ice sheet. Then, using the extended velocity discussed in Section A.1, the level set in every subsequent time step is guaranteed to remain a signed distance function.
Every node at the edge of the computational domain must be assigned a suitable boundary condition. We choose to use a special form of linear extrapolation described by Mitchell (Reference Mitchell2004) that adds an appropriate number of ‘ghost nodes’ beyond the edge of the grid when working on nodes near the edge. The values of φ at ghost nodes are determined by linear extrapolation from the computational boundary with a slope direction that matches the sign of the level set at the boundary node. Suppose (x i, z j) for i = 1, 2, …, p and j = 1, 2, …, q are nodes in the domain and (x 0, z j), (x p+1, z j), (x i, z 0) and (x i, z q+1) are ghost nodes, then the values at the ghost nodes left of the domain are given by
for j = 1, 2, …, q, and we similarly define the values at the other ghost nodes (x p+1, z j), (x i, z 0) and (x i, z q+1). This is not a traditional PDE boundary condition, however, it is quite useful in level-set computations for domains with inflow boundaries that have no physically appropriate boundary conditions, as it remains stable whereas regular linear extrapolation may cause stability issues (Mitchell, Reference Mitchell2004).
A.3. Re-initialization using the Fast Marching Method
The level sets that are located near the zero level set move with speeds that can considerably distort and stretch the level-set function φ. Under such circumstances, φ can develop noisy features and steep gradients that are detrimental to finite-difference approximations and so can fail to preserve the signed distance function. As a result, it may be necessary to periodically re-initialize the level-set function, which involves stopping the calculation at some point in time and rebuilding the level-set function according to the signed distance function (Sethian, Reference Sethian1999b), thereby ensuring that φ remains smooth enough to allow its spatial derivatives to be computed with sufficient accuracy.
Although there are several ways this re-initialization could be carried out in practice, we implement the Fast Marching Method (FMM), which is known to be very effective for this purpose. The FMM offers a fast approach for rebuilding φ having computational cost of ${\cal O}\lpar {N\log N} \rpar$
, where N is the total number of grid points (Adalsteinsson and Sethian, Reference Adalsteinsson and Sethian1999). To re-initialize the signed distance function φ, the FMM solves the eikonal equation ${\parallel} \nabla \varphi \parallel{ = } 1$
on either side of the interface ∂Ω (Sethian, Reference Sethian1999a; Vogl, Reference Vogl2016). The FMM algorithm considers three categories of grid points: CLOSE, FAR and ACCEPTED. The ACCEPTED points are initially assigned to grid nodes that immediately surround the zero level set, CLOSE points are then one gridpoint further away, and the remaining nodes in the domain are labelled FAR (see Fig. 11). The shortest path from each ACCEPTED grid node to the contour of the zero level set (ice–air or ice–water interface) is determined using a non-linear optimization solver and used to assign the signed-distance value to each ACCEPTED point. The procedure continues with the following steps that efficiently and systematically marches CLOSE and FAR points to ACCEPTED in order to assign the signed-distance value at all grid points in the domain.
(1) The signed-distance value of CLOSE points is calculated based on the known signed-distance at neighbouring ACCEPTED points and the grid size.
(2) Let TRIAL be the CLOSE point with the smallest value of φ.
(3) Any FAR points that directly neighbour TRIAL are relabelled CLOSE.
(4) Relabel TRIAL points to ACCEPTED.
(5) Repeat steps 1–4 until all points become ACCEPTED.
Fig. 11. Initialization of the Fast Marching Method, where ‘○’ denote the initial ACCEPTED points, ‘x’ the CLOSE points and ‘*’ indicating the FAR points.
The accuracy of this approach means re-initialization is required less often. For our simulations, we used the FMM to rebuild the level-set function every 50–100 time steps.
Appendix B. Ice-sheet equations
B.1. Shallow ice approximation (SIA)
The SIA treats the ice sheet as a shallow film that flows and spreads under its own weight (Hutter, Reference Hutter1983). We denote the sheet thickness at position (x, y) and time t by H(x, y, t). Then, in terms of z measured vertically upward from sea level, the height of the upper surface in contact with the atmosphere is represented by z = h(x, y, t) and the lower bedrock surface by z = b(x, y) (which has zero normal velocity, assuming negligible melt there). Referring to Figure 12, these three height variables are related by
The ice deformation is determined by the incompressible Stokes equations, coupled with Glen's flow law (Glen, Reference Glen1958) under the shallow ice assumption. In the isothermal case, the horizontal velocity components U = (u, v) as
in the case where there is no sliding relative to the underlying bedrock. The other variables and parameters in the equations are the gravitational acceleration g = (0, 0, − g), ice density ρ, creep parameter A and Glen's law exponent n ≈ 3.
Fig. 12. Geometry of the shallow ice-sheet flow problem.
For a grounded ice sheet that is radially symmetric about the ice divide, denoted r = 0. The radial symmetry implies that the sheet geometry depends only on r so that h = h(r, t), H = H(r, t) and b = b(r). For the case of non-sliding ice, the radial velocity is ${\bi U} = u \hat{\bi r}$
where $\hat{\bi r}$
denotes the unit vector in the radial direction. At the ice divide (r = 0), a symmetry condition is imposed
To handle the case when there may be some slip at the ice-sheet base, we consider a friction law that relates basal stress τ b to the sliding velocity u b at the bed by means of the relationship $\tau _{\rm b} = f\lpar {u_{\rm b}} \rpar = Cu_{\rm b}^m$
, where the bed friction parameter C depends on the local bed roughness and a bed friction exponent, m = (1/n) (Schoof and Hewitt, Reference Schoof and Hewitt2013). Using the stress balance and Glen's flow law, the sliding velocity is u b = −((ρgH/C)(∂h/∂r))1/m, so that the radial velocity can be written in a more general form that captures both the sliding and non-sliding cases:
The vertical velocity w(r, z, t) may then be obtained from u(r, z, t) using the incompressibility condition, which in cylindrical coordinates is
Integrating this equation in z and applying the vertical no-flow boundary condition at the bed z = b(r) yields the corresponding expression for
B.2. Shallow shelf approximation (SSA)
The SSA is used to consider a 2-D symmetrical marine ice sheet. Denoting the horizontal coordinate in the flow direction by x, symmetry implies that
The ice stream or grounded portion of the marine ice sheet occupies the region 0 ≤ x < x g, where x g denotes the grounding line position. The momentum conservation equation of SSA for the grounded ice sheet (0 ≤ x < x g) was derived by MacAyeal (Reference MacAyeal1989) as
where h = b + H as in the SIA model, C is the bed friction parameter and H ≥ −(ρ w/ρ)b where ρ w is the sea water density. The SSA Eqn (B6) represents a balance between longitudinal strain rates that are determined by the integrated ice hardness (the coefficient A −1/nH), the slipperiness of the bed (the coefficient C and exponent m), and the geometry of the ice sheet (the thickness H and surface slope ∂h/∂x).
For the unbuttressed freely floating ice shelf that occupies the region x g < x < x c, where x = x c denotes the calving front, we have H < −(ρ w/ρ)b and h = (1 − ρ/ρ w)H. There is no basal friction and so the term C|u|m−1u vanishes and the sole driving stress for ice shelves is ρ(1 − ρ/ρ w)gH(∂H/∂x), giving the momentum conservation equation for the shelf (x g < x < x c) as
At the calving front, there is an imbalance between hydrostatic pressures in ice and water due to the buoyancy of ice, hence
At the grounding line, where ice stream couples to ice shelf, we have H = −(ρ w/ρ)b, as well as the boundary condition Eqn (B8) applied at x = x g.
Assuming ice to be an incompressible material, the vertical velocity for the ice stream with rigid bedrock ((∂b/∂t) = 0) and no melting at the bottom surface is determined by
For the ice shelf, the vertical velocity is given by
where w sl is the vertical velocity at sea level that can be determined from the known surface and basal mass balances, ${\cal M}_{\rm h}$
and ${\cal M}_{\rm b}$
, respectively (Greve and Blatter, Reference Greve and Blatter2009).
For a given ice thickness H(x) (determined by the LSM in our case) the velocity u(x) is determined by solving the non-linear partial differential Eqns (B6) and (B7). Our solution approach implements an iterative numerical method, often called a Picard iteration. Denote the current velocity iterate as u (k) and the previous iterate as u (k−1), then the Picard iteration for Eqn (B6) (and similarly for Eqn (B7)) is:
where W (k−1) = 2A −1/nH|∂u/∂x (k−1)|1/n−1. For grounded ice, 0 < x < x g, we assume the ice is held by basal resistance only to obtain the initial velocity estimate, u (0)(x) = ( − C −1ρgH(∂h/∂x))1/m and the boundary conditions from Eqns (B5) and (B8) are u(0) = 0, and (∂u/∂x)(x g) = A((1/4)ρ(1 − ρ/ρ w)gH)n. For floating ice, x g < x < x c, an initial guess for velocity comes from assuming a uniform strain rate provided by the calving front condition:
where u g denotes the ice velocity at the grounding line, and the boundary conditions are u(x g) = u g and (∂u/∂x)(x c) = A((1/4)ρ(1 − ρ/ρ w)gH)n.
Open access
