A smoothed-particle hydrodynamics model for ice-sheet and ice-shelf dynamics

Abstract

Mathematical modeling of ice sheets is complicated by the nonlinearity of the governing equations and boundary conditions. Standard grid-based methods require complex front-tracking techniques and have a limited capability to handle large material deformations and abrupt changes in bottom topography. Consequently, numerical methods are usually restricted to shallow ice-sheet and ice-shelf approximations. We propose a new smoothed-particle hydrodynamics (SPH) model for coupled ice-sheet and ice-shelf dynamics. SPH, a fully Lagrangian particle method, is highly scalable and its Lagrangian nature and meshless discretization are well suited to the simulation of free surface flows, large material deformation and material fragmentation. In this paper, we use the SPH model to study ice-sheet/ice-shelf behavior, and the dynamics of the grounding line. The steady-state position of the grounding line obtained from SPH simulations is in good agreement with laboratory observations for a wide range of simulated bedrock slopes and density ratios, similar to those of ice and sea water. The numerical accuracy of the SPH algorithm is verified by simulating the plane-shear flow of two immiscible fluids and the propagation of a highly viscous blob of fluid along a horizontal surface. In the experiment, the ice was represented with a viscous Newtonian fluid. For consistency, in the described SPH model the ice is also modeled as a viscous Newtonian fluid. Typically, ice sheets are modeled as a non-Newtonian fluid, accounting for the changes in the mechanical properties of the ice. Implementation of a non-Newtonian rheology in the SPH model is the subject of our ongoing research.

Introduction

Most existing ice-sheet models use grid-based Eulerian discretizations and are based on various approximations to a steady-state Stokes equation (Marshall, 2005). The approximations are associated with the shallowness, namely the small depth-to-width ratio, of ice sheets to simplify the equations and reduce computational cost (Bueler and Brown, 2009). The simplest models use the first-order shallow-ice approximation (Morland and Johnson, 1980) and the shallow-shelf approximation (Weis and others, 1999). Under certain conditions, shallow-ice and shallow-shelf approximations may lead to significant errors. Examples of this include: large ice-sheet aspect ratio and/or large bedrock slope; tidewater glaciers; ice streams; surge dynamics; the dynamics of flow across the grounding line; and the dynamics in the vicinity of ice-sheet divides (Marshall, 2005). Recently, several models based on the Stokes equation or its high-order approximations have been proposed (Bueler and Brown, 2009; Durand and others, 2009; Pollard and DeConto, 2009; Price and others, 2011; Seddik and others, 2011; Seroussi and others, 2011), but the presence of the free surface (ice/air interface) still remains a major computational challenge for the Eulerian grid-based methods.

We propose a new smoothed-particle hydrodynamics (SPH) model for coupled ice-sheet and ice-shelf dynamics. SPH is a fully Lagrangian method that uses meshless (particle) discretization of the computational domain (Monaghan, 2005). SPH particles are used as discretization points for solving conservation equations. In the proposed SPH model, we solve the full momentum conservation equation, coupled with the continuity equation subject to the free surface boundary condition at the ice/water interface and the no-slip boundary condition at the ice/rock interface. Continuity of velocity and the viscous stress are imposed at the ice/water interface. Lagrangian particle methods do not require interface-tracking algorithms for modeling free- surface (e.g. Monaghan, 1994; Tartakovsky and Meakin, 2005) and moving boundary (multiphase) problems (e.g. Colagrossi and Landrini, 2003; Tartakovsky and others, 2007, 2009). In SPH, each fluid has its own set of particles and there is no need to define the interface. This makes SPH very efficient for ice-sheet modeling.

The marine ice sheets terminate in the ocean and exhibit a balance between snowfall inland and calving of icebergs into the ocean. Ice sheets thin as they flow toward the ocean and eventually detach from the bedrock to form freely floating ice shelves. The location of the detachment is called the grounding line. Recently, a laboratory experiment was performed by Robison and others (2010) to explore the dynamics of the grounding line between marine ice sheets and the freely floating ice shelves into which they develop. In this experiment, the ice sheets and shelves were represented by a viscous Newtonian fluid flowing down a ramp into a tank containing a denser and much less viscous fluid. The viscous fluid floated off and formed a shelf at the grounding line, which proceeded steadily (without retreating) and eventually reached a steady position. The experiment captures the fundamental dynamical properties of flow across a grounding line between a horizontal shear- dominated ice sheet and an extensional deviatoric shear- dominated ice shelf. We use the SPH model to simulate processes studied in the experiment, and we find excellent agreement between the steady positions of the grounding line obtained numerically and experimentally for a wide range of ramp (bedrock) slopes, density ratios and ice-sheet flow rates. Numerical results also agree with theoretical prediction for the location of the grounding line obtained by Robison and others (2010).

In order to obtain a consistent comparison with the experimental data, in the SPH model the ice sheet and ice shelf are modeled as a viscous Newtonian fluid. In reality, ice sheets have a complex rheology and are traditionally modeled as a non-Newtonian fluid with rheology represented by Glen’s power law (Huybrechts and others, 1996; Goldsby and Kohlsted, 2001). However, the rheology of the ice sheet and ice shelf primarily affects the vertical velocity profiles within them, and the effect of the rheology on the depth-averaged quantities and position of the grounding line is unknown. Implementation of the non-Newtonian rheology in the SPH method is relatively straightforward (Shao and Lo, 2003; Hosseini and others, 2007). The effect of the non-Newtonian rheological properties of ice on the ice-sheet dynamics is a subject of our ongoing research.

This paper is organized as follows. First we describe the SPH model. Next, the model is verified against finite- difference and analytical solutions for the plane-shear flow of two immiscible fluids and the propagation of gravity current along a rigid solid surface. Finally, we compare the dynamics of the grounding line obtained from SPH simulations with the laboratory experiment of Robison and others (2010).

SPH Discretization of the Navier-Stokes and Continuity Equations

SPH equations of motion

In the proposed model, the SPH method is used to discretize the governing Navier-Stokes and continuity equations. SPH uses a meshless (particle) discretization of the computational domain (see Monaghan, 2005, for a review). The SPH interpolation scheme, allows approximation of a continuous field, A(r), using the values of A at a set of discretization points. Here r _i, is the position of discretization point i, A_i = A(r _i ), and ρ_i and m_i are the density and the mass of the phase associated with point i. Because each point possesses a mass and volume, m_i /ρ_i it is natural to think of discretization points as physical particles. W is the bell-shaped SPH weighting function with compact support kh (W(|r - r _i | > kh, h) = 0), where |r| is the magnitude of a vector r, the value of k depends on a functional form of W and h depends on position r. The interpolation scheme assumes summation over all the SPH particles within a distance kh of r. We choose a cubic spline kernel, which is defined as

(1)

where x = |r|/h, D is the spatial dimension and a_D is a constant that assures proper normalization of the smoothing function (a_D = 2/3, 10/7π and 1/π for D = 1, 2 and 3, respectively). Here only particles within 2h of the central particle contribute to the smoothing kernel, which is with continuous second derivatives for all r. The gradient of the kernel is

(2)

The SPH approximation of continuous fields and their spatial derivatives allow the continuity and Navier-Stokes equations to be written in the form of a system of ordinary differential equations (Monaghan, 2005):

(3)

(4)

Where

(5)

fluid pressure and μ_i is the fluid viscosity at position r_i; ρ ¹ ₀ and ρ ² ₀ represent the reference densities of fluids 1 and 2, respectively. In Eqn (4), the first term inside the summation is a repulsive force acting between particles i and j that is obtained from the discretization of the pressure gradient term in the Navier-Stokes equation. The second term is a repulsive force that is needed only if the density difference between two fluids is large (e.g. a density ratio of 5 to 10 or more (Monaghan, 1994)). The third term is the viscous force defined in Eqn (5). The attractive feature of Eqn (4) is that the forces have antisymmetric form with respect to indexes i and j (forces acting between the SPH particles satisfy Newton’s third law). Because of this, Eqn (4) satisfies the linear and angular momentum conservation laws exactly.

Ice, water, bedrock and the ocean bottom are represented by different sets of particles with different masses, viscosities and densities. In Eqns (3) and (4), Σ _j indicates summation over all neighboring particles of particle i. Particles representing solid walls are fixed in space. They enter into the calculation of the densities of fluid particles (Eqn (3)) and forces acting on the fluid particles (Eqn (4)). v is the fluid velocity vector, P is the pressure, g is the gravitational acceleration vector and π is the artificial viscosity.

For computational efficiency, we approximate an incompressible fluid as a slightly compressible fluid. The equation of state, was used to close the system of Eqns (3) and (4). In the equations of state, ρ ₀ is the reference density, and P ₀ is the magnitude of pressure given by Here c_i , is the speed of sound of particle i. This equation of state limits the relative density fluctuation to Under such conditions, the SPH fluids behave like incompressible fluids.

One layer of solid (boundary) particles is placed along the fluid/solid boundary, and forces acting between the boundary particles and fluid particles near the boundary are added to Eqn (4) to ensure that the normal velocity is zero (Monaghan and Kajtar, 2009). To satisfy the no-slip boundary condition, the boundary particles are included in the calculation of viscous forces acting on the fluid particles. The boundary force is calculated as

(6)

where particle j is the neighbor boundary particle of fluid particle is proportional to a Wendland one-dimensional (1-D) cubic kernel. The boundary particle spacing, Δp _b, is less than the fluid particle spacing by a factor of 2 or more. Monaghan and Kajtar (2009) demonstrated that for Δp _b satisfying this condition, the total boundary force on a fluid SPH particle acts in the direction perpendicular to the boundary within a negligible error.

Finally, the smoothing length, h, determines the resolution and the number of neighbors that contribute to the properties at a point. The efficiency and the accuracy would therefore be greater if h was chosen so that it depended on the local particle number density. Typically, h_i is calculated from

(7)

where D is the number of dimensions (Monaghan, 1992). A symmetric kernel is obtained by replacing h with the arithmetic mean of the smoothing lengths for the two particles as h = (h_i + h_j )/2 (Monaghan, 1992).

Integration and time-stepping

The second-order leapfrog integrator is used to integrate the SPH equations of motion. For known position, r, velocity, v, and acceleration, f, at time t, the leapfrog integrator gives the positions and velocities at time t + Δt as

(8)

The smoothing length, h_i (t), and density, ρ_i (t), are integrated as

(9)

where R_i (t) = dρ _i (t)/dt.

With this integrator, it is crucial that the time-step size is chosen correctly, both to ensure the accuracy of the evolution and to ensure numerical stability. The time-step is chosen as the minimum of the Courant-Friedrichs-Lewy condition,

(10)

the viscosity condition (Cleary and Monaghan, 1999),

(11)

And the boundary force condition (Monaghan and kajtar, 2009),

(12)

Where V_b = max(|v _ij |, V _max).

Numerical Results

To verify our SPH algorithm we applied two tests: (1) the plane-shear flow of two fluids with different densities and viscosities, a fixed, no-slip, horizontal bottom boundary and an upper free surface; and (2) the evolution of a blob of viscous fluid spreading along a horizontal boundary.

Plane-shear flow of two fluids

We study a two-dimensional (2-D) plane-shear flow of two immiscible fluids with a free surface. A layer of fluid 1 is bounded by an impermeable plane boundary at the bottom moving with velocity u ₀ and a layer of fluid 2 at the top. The shear flow is initiated by a sudden stop of the bottom plane. For this problem, the (2-D) Navier-Stokes and continuity equations can be reduced to 1-D equations:

(13)

Subject to initial and boundary conditions:

(14)

The parameters of the two fluids are given in Table 1. We use the SPH model to solve the 2-D continuity and Navier-Stokes equations and use finite-difference (FD) solutions of Eqns (13) and (14) to verify the SPH solution. The FD scheme has first- order accuracy in time and second-order accuracy in space and has the grid size Δy = 0.005. In the SPH simulations, pressure and velocity are assumed to be periodic in the x- direction. The computational domain in the x-direction is set to 10h (to avoid the effect of the periodic boundary condition), and SPH particles exiting domain at x = 10h with velocity vi are returned in the domain at x = 0 with the same velocity. Figure 1 compares an SPH solution of the continuity and Navier-Stokes equations, obtained with different resolution, Δp, with a FD solution of the equivalent 1-D Eqns (13) and (14) for three different times. In the SPH solution, Δp is the initial spacing between the SPH particles. The figure shows that SPH solutions approach the FD solution with decreasing Δp (increasing resolution of the SPH solution). Figure 2 shows the L ₂ relative error versus Δp at three times. At t = 40.7s, the error is found to decrease as Δp ^1.0, while at t = 309.3 s it decreases as Δp ^1.5, and at t = 773.2 s the error decreases as Δp ^0.6. We find that the shear free surface flow is simulated accurately by the SPH model and the accuracy is improved with increasing spatial resolution.

Table 1. Parameters of the two fluids

Fig. 1. The velocity profiles of the flow at different times.

Fig. 2. The L ₂ relative error of the velocity.

Gravity current along a rigid surface

Next, we study a viscous gravity current propagating over a rigid horizontal surface, a problem that is often used to verify ice-sheet codes (Bueler and others, 2007). This test imitates a volume of ice that spreads as a sheet over a horizontal surface under gravity. In the absence of surface tension, this problem allows a 2-D analytical solution (Huppert, 1982). We simulate the gravity current propagation with the SPH model in two dimensions and use the analytical solution to verify the SPH model.

The analytical solution is obtained using the lubrication- theory approximation as follows. A current of density ρ is considered to be propagating over a rigid horizontal surface with depth H(x, t). Here the vertical velocities are negligibly small and the pressure is hydrostatic. The balance between the pressure gradient and the viscous forces then gives the governing nonlinear partial differential equation for H(x, t) (Huppert, 1982):

(15)

where ν is the kinematic viscosity of the viscous fluid. For an incompressible fluid the continuity equation yields

(16)

where q is the constant volume of the current. A similarity solution of H(x, t) can be found in terms of the similarity variable (Huppert, 1982),

(17)

and the solution is then expressed as (Huppert, 1982)

(18)

where ηL is the value of η at x = L(t), and L(t) satisfies H(L, t) ≡ 0. The function ϕ(y) with y = η/ηL satisfies

(19)

and

(20)

The analytical solution of Eqns (19) and (20) is obtained as:

(21)

Figure 3 compares the ice-sheet thickness over time, H(x, t), obtained from the SPH simulation and the analytical solution at different times. Figure 4 shows the margin length, L(t), as function of time obtained analytically, and from the SPH simulations with different resolutions. These figures demonstrate that the SPH solution converges to the analytical solution with decreasing Δp (increasing resolution).

Fig. 3. The surface profile of the ice sheet over time as the sheet spreads under gravity. The dashed curves represent the solution of Huppert (1982). Here, L ₀ is the initial margin length and Δp is the initial particle spacing.

Fig. 4. Length of the ice sheet as a function of time obtained analytically and from SPH simulations with different initial particle spacing, Δp.

Dynamics of the grounding line

To demonstrate the capabilities of the SPH ice-sheet model, we use the model to investigate the dynamics of the grounding line, and compare it with the laboratory observations of Robison and others (2010) and their theoretical prediction.

In Robison and others’ (2010) laboratory experiment, a Newtonian viscous fluid (golden syrup) was supplied from a reservoir onto a ramp, down which it flowed into a tank filled with aqueous potassium carbonate solution. The syrup is more viscous but less dense than potassium carbonate solution and is able to float on the top of the ‘ocean’. A constant head of golden syrup was maintained in the supply chamber. The flux rate was adjusted by altering the head height. At a certain point (the grounding line) along the ramp, the syrup floated off and formed a shelf. A side-view camera was used to record and measure the positions of the grounding line. The grounding line was found to move down the ramp but eventually reach a steady position.

By approximating the flow of an ice sheet as shear- dominated and the flow of ice shelves as extensional with zero shear to leading order, Robison and others (2010) obtained an approximate analytical solution for the steady position of the grounding line,

(22)

where the dimensionless grounding line

(23)

and

(24)

Here α is the bed slope, g is the gravity acceleration, g′ = g(ρ _w - ρ)/ρ _w is the reduced gravity and ρ and ν are the density and kinematic viscosity of the viscous fluid, respectively; ρ _w is the density of the dense fluid and q ₀ is the flux of viscous fluid.

The set-up of a typical simulation is illustrated in Figure 5, which has the same dimensions as the laboratory set-up of Robison and others (2010). Viscous light (gray) fluid flows down a rigid ramp with slope α into the denser but less viscous ‘ocean’ (blue). The gray fluid is supplied through a funnel from a reservoir where a constant level of fluid is maintained. An overflow tank is attached on the right to keep the constant ‘sea level’.

Fig. 5. The set-up used in our simulation. A viscous light fluid (gray) flows down onto a ramp and into a dense fluid (blue).

To initialize a simulation, the particles were placed on the vertices of a grid of squares with lattice size Δp. Cray particles were placed inside the reservoir and the funnel, and blue particles were placed inside the basin. The overflow tank was initially empty. One layer of ‘boundary’ particles was placed along the solid boundary (the dark gray line) with spacing Δp _b = Δp/3. At time t = 0, gravity drove the fluid particles (gray) to flow down the ramp and into the basin of blue particles. The flux of gray particles was controlled with the friction force, -kV_i, added into the momentum conservation, Eqn (4), for all particles within the funnel. The constant level of viscousfluid in the reservoir was maintained by inserting gray particles whenever the level drop exceeded Δp. The number of particles in the simulations ranged from 15 000 to 20 000, and the CPU time was ~4 hours on a standard single-processor desktop computer. The typical time-step used in the simulations was Δt = 0.0001 s.

Figure 6 compares the steady position of the grounding line obtained from the SPH simulations and measured in the experiment (Robison and others, 2010). The good overall comparison between numerical and experimental results indicates that the SPH model accurately captures the position of the contact line for a wide range of slopes, fluxes and density ratios.

Fig. 6. The steady positions of the grounding line calculated in our simulations, compared with those measured in the experiment of Robison and others (2010). Different densities, ρ, and bed slopes, α, were chosen according to the values used in the experiment (Table 2) for varied values of A + ∊. x _G is dimensionless, as defined in Eqn (23).

Table 2. Values used for the input parameters in the grounding line simulations

The convergence of the simulation was also examined. The position of the grounding line as a function of time was calculated at different spatial resolutions, Δp. In all the simulations, the position of the grounding line increased with time to an asymptotic value. Figure 7 shows that the asymptotic values (the steady-state position of the grounding line) converge with decreasing Δp. We find that a sufficiently accurate value of the steady position of the grounding line can be obtained with Δp = 0.026, and we chose this resolution for all the simulations.

Fig. 7. The position of the grounding line as a function of time calculated at different spatial resolutions, Δp. All the variables are made dimensionless by Eqn (23).

The black solid curve in Figure 6 is a linear fit to the experimental data, has a slope -β = -1.05 and follows a scaling law:

(25)

The position of the grounding line, x _G, obtained from the SPH simulations has a similar scaling behavior, with β = 1.1. The slight discrepancy in the values of β for experimental and numerical data and the significant scatter of the experimental data around the line given by Eqn (25) can be explained by several factors. The SPH simulations are 2-D, while the experiments are 3-D. In the experiment, the side-walls of the tank exert viscous shear stresses on the fluid, and, in some experiments, propagation of the syrup on the water surface was unstable (fingering of the syrup was observed). Both the experimental and numerical data show scaling behavior similar to the power law, Eqn (22), obtained theoretically by Robison and others (2010). However, the theory results in γ = 1, while both the experimental and numerical results predict γ ≈ 3. This discrepancy can be explained by the shallow-ice approximation used in the theoretical model.

Conclusions

A SPH model for ice-sheet dynamics was developed. The model was verified using analytical and numerical solutions for plane-shear flow of two immiscible fluids and the propagation of a gravity current along a rigid solid surface. The SPH simulations converge to exact solutions with sufficientspatial resolution. We have shown that the essential dynamics of marine ice sheets can be captured by SPH simulations. The SPH model was validated using data from an experiment of Robison and others (2010), mimicking the flow of an ice sheet and ice shelf. The experimental and numerical results show the migration of the grounding line between a grounded sheet and a freely floating shelf toward a steady position. The steady-state positions of the grounding line, obtained from the SPH simulations and the laboratory experiments, agree well for a wide range of bedrock slopes, density ratios and ice-sheet fluxes.

Typically, ice sheets are modeled as a non-Newtonian fluid (Huybrechts and others, 1996; Coldsby and Kohlsted, 2001), accounting for nonlinear dependence of the deviatoric stress on the strain rate. While in contact with bedrock, the flow of an ice sheet is resisted by shear stresses, both internal and at its base. However, beyond the grounding line, the ice shelf is characterized by a weak internal straining flow, with negligible tangential stress exerted upon it by the air above it or the ocean beneath. The effect of the non-Newtonian rheological properties of ice on the ice-sheet dynamics is the subject of our ongoing research.

Several important processes, such as basal melt- ing/freezing, ice accumulation and basal sliding were absent in the experiments of Robison and others (2010). Including these processes in the SPH ice-sheet model is the subject of our future research. An SPH framework for modeling melting/freezing processes was proposed by Monaghan and others (2005). The SPH deposition model for mineral growth of Tartakovsky and others (2007, 2008) can be used to model ice accumulation due to precipitation. Basal sliding in the SPH model can be modeled by prescribing the basal sliding velocity to the boundary particles. Implementation of a Robin boundary condition, which is commonly used in ice-sheet models to describe basal sliding (Cagliardini and others, 2007), in SPH is described by Ryan and others (2010).