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 Reference MonaghanMonaghan, 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 (Reference MonaghanMonaghan, 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 (Reference MonaghanMonaghan, 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 (Reference Monaghan and KajtarMonaghan 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. Reference Monaghan and KajtarMonaghan 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 (Reference MonaghanMonaghan, 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 (Reference MonaghanMonaghan, 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 (Reference Cleary and MonaghanCleary and Monaghan, 1999),

(11)

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

(12)

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

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 (Reference Bueler, Brown and LingleBueler 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 (Reference HuppertHuppert, 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) (Reference HuppertHuppert, 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 (Reference HuppertHuppert, 1982),

(17)

and the solution is then expressed as (Reference HuppertHuppert, 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 Reference HuppertHuppert (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 Reference Robison, Huppert and WorsterRobison and others (2010) and their theoretical prediction.

In Reference Robison, Huppert and WorsterRobison 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, Reference Robison, Huppert and WorsterRobison 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 Reference Robison, Huppert and WorsterRobison 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 (Reference Robison, Huppert and WorsterRobison 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 Reference Robison, Huppert and WorsterRobison 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 Reference Robison, Huppert and WorsterRobison 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.