Articles | Volume 18, issue 3
Research article
04 Mar 2024
Research article |  | 04 Mar 2024

Smoothed particle hydrodynamics implementation of the standard viscous–plastic sea-ice model and validation in simple idealized experiments

Oreste Marquis, Bruno Tremblay, Jean-François Lemieux, and Mohammed Islam

The viscous–plastic (VP) rheology with an elliptical yield curve and normal flow rule is implemented in a Lagrangian modelling framework using the smoothed particle hydrodynamics (SPH) meshfree method. Results show, from a perturbation analysis of SPH sea-ice dynamic equations, that the classical SPH particle density formulation expressed as a function of sea-ice concentration and mean ice thickness leads to incorrect plastic wave speed. We propose a new formulation for particle density that gives a plastic wave speed in line with theory. In all cases, the plastic wave in the SPH framework is dispersive and depends on the smoothing length (i.e., the spatial resolution) and on the SPH kernel employed in contrast to its finite-difference method (FDM) implementation counterpart. The steady-state solution for the simple 1D ridging experiment is in agreement with the analytical solution within an error of 1 %. SPH is also able to simulate a stable upstream ice arch in an idealized domain representing the Nares Strait in a low-wind regime (5.3 m s−1) with an ellipse aspect ratio of 2, an average thickness of 1 m and free-slip boundary conditions in opposition to the FDM implementation that requires higher shear strength to simulate it. In higher-wind regimes (7.5 m s−1) no stable ice arches are simulated – unless the thickness is increased – and the ice arch formation showed no dependence on the size of particles, in contrast to what is observed in the discrete-element framework. Finally, the SPH framework is explicit, can take full advantage of parallel processing capabilities and shows potential for pan-Arctic climate simulations.

1 Introduction

Sea ice is an important component of the Earth’s system to consider for accurate climate projection. Generally, numerical models used for geophysical sea-ice simulations and projections are based on a system of differential equations assuming a continuum. The equations that predict the sea-ice dynamics are a combination of the momentum equations, which describe the drift of sea ice under external and internal forces, and the continuity equations which ensure mass conservation. The external forces (per unit area) generally include surface air stress, water drag, sea surface tilt and the Coriolis effect, and the internal forces are related to the ice stress term. This internal stress term is based on various constitutive relations which can differ between models. The more commonly used constitutive laws are the standard viscous–plastic model (Hibler1979) or modifications thereof (e.g., elastic–viscous–plastic or EVP and elastic–plastic–anisotropic or EPA; Hunke and Dukowicz1997; Tsamados et al.2013). They are typically discretized on an Eulerian mesh using the finite-difference method (FDM). FDM is the simplest method to discretize and solve partial differential equations numerically. However, it is based on a local Taylor series expansion to approximate the continuum equations and construct a topologically rectangular network of relations between nodes (e.g., Arakawa grids).

Even though the VP (and EVP) rheologies are commonly used to describe sea-ice dynamics and are able to capture important large-scale deformation features (Bouchat et al.2022; Hutter et al.2022), they still have difficulties representing smaller-scale properties (Schulson2004; Weiss et al.2007; Coon et al.2007) such as linear kinematic features (LKFs) unless run at very high resolution (≈2kmRingeisen et al.2019; Hutter et al.2022). To improve the simulation of small-scale ice features and to alleviate the problem of FDM with complex geometries (Peiró and Sherwin2005), the community also considered new sea-ice rheologies (Schreyer et al.2006; Girard et al.2011; Dansereau et al.2016; Ringeisen et al.2019) and explored different space discretization frameworks like the finite-element method (FEM) (Rampal et al.2016; Mehlmann et al.2021), the finite-volume method (FVM) (Hutchings et al.2004; Losch et al.2010; Adcroft et al.2019) or the discrete-element method (DEM) (Hopkins and Thorndike2006; Herman2016; Damsgaard et al.2018).

In recent decades, the spatial resolution of sea-ice models has become comparable to the characteristic length of the ice floes. This makes the continuum assumption of current FDM, FVM and FEM models questionable. Also, Eulerian models are known to have difficulties determining the precise locations of inhomogeneity, free surfaces, deformable boundaries and moving interfaces (Liu and Liu2010). These shortcomings have led to an increased interest in the DEM approach. Another advantage of using DEMs is that the granularity of the material (Overland et al.1998) is directly represented using discrete rigid bodies from which the physical interactions are calculated explicitly in the hope that large-scale properties naturally emerge. In practice, the emergent properties of a granular medium still depend on the assumed floe size and the nature of collisions in contrast to the continuous numerical methods which can indirectly account for floe interactions through the formulation of a constitutive law. Nevertheless, DEMs easily capture the formation of cracks, leads and large deformations, making them a consistent framework for the numerical simulation of granular material like sea ice (Fleissner et al.2007).

Despite the shortcomings of the continuum approaches, FDM, FVM and FEM are still the most commonly used frameworks in the community because they have been developed and tested for a longer period and are well understood, computationally more efficient and easily coupled for large-scale simulations. In an attempt to take advantage of both continuum and discrete formulations, blends between the two approaches have been proposed – e.g., the finite- and discrete-element methods (Lilja et al.2021) or the material-point method (Sulsky et al.2007). Those framework, however, still use a mesh to solve the dynamic equations in addition to considering sea ice as discrete elements, making them even more computationally expensive. Finally, a fairly new approach in the context of sea-ice modelling – also taking from both continuum and discrete frameworks – uses a Lagrangian meshfree continuous method called smoothed particle hydrodynamics (SPH) developed by Lucy (1977) and Gingold and Monaghan (1977). This meshfree method is known to facilitates the numerical treatment and description of free surfaces (Liu and Liu2010) which are common in sea-ice dynamics with polynyas, LKFs, free-drifting ice floes and unbounded ice extent. As in DEM, the physical quantities are carried out by particles in space (an analogy for ice floes in the real world) but evolve according to the same dynamic equations used in the continuum approach. Furthermore, the method has the advantage of treating the system of equations in a Lagrangian framework discretized explicitly, making it well-suited for parallelization.

SPH has been applied successfully for modelling other granular materials such as sand, gravels and soils (Salehizadeh and Shafiei2019; Yang et al.2020; Sheikh et al.2020). In the context of mesoscale and larger sea-ice modelling, Gutfraind and Savage (1997) initiated the SPH study of sea-ice dynamics using a VP rheology based on a Mohr–Coulomb failure criterion. The ice concentration and thickness were fixed at 100 % and 1 m with a continuity equation expressed in terms of a particle density. The internal ice strength between particles was derived diagnostically from ice density assuming ice was a nearly incompressible material. Later, Wang et al. (1998) developed a sea-ice model of the Bohai Sea (east coast of China) using an SPH viscous–plastic rheology (Hibler1979) with continuity equations for ice concentration and mean thickness, as well as ice strength calculated from static ice jam theory (Shen et al.1990). Following Wang et al. (1998), Ji et al. (2005) implemented a new viscoelastic–plastic rheology that was in better agreement with observations from the Bohai Sea. Recently, Staroszczyk (2017) proposed a sea-ice model considering ice behaving as a compressible non-linear viscous material with a (dimensionless) contact-length-dependent parameterization for floe collisions and rafting (Gray and Morland1994; Morland and Staroszczyk1998). In all of the above, except for Gutfraind and Savage (1997), the same ice particle density definition is used.

In this work, we use the standard VP sea-ice model with an elliptical yield curve and normal flow rule (Hibler1979) as a proof of concept. Further development of the SPH model should consider a broader range of rheologies. We propose a reformulation of the ice particle density that is internally consistent with the model physics. One goal of the study is to investigate differences coming from the numerical framework. To this end, we theoretically investigate the plastic wave propagation, a fundamental property of a sea-ice physical model, using a 1D perturbation analysis and we test the model in a ridging and ice arch experiment following earlier works by Lipscomb et al. (2007), Dumont et al. (2009), Rabatel et al. (2015), Dansereau et al. (2017), Williams et al. (2017), Damsgaard et al. (2018), Ranta et al. (2018), Plante et al. (2020), and West et al. (2022). We chose to investigate the SPH method performance on a ridging experiment since it has an analytical steady-state solution that can be used to establish the model accuracy and it is possible to evaluate whether the coupling with the mass equations is coherent. We also test SPH performance on ice arch simulations since this classic problem is an example of large-scale features resulting from small-scale interactions involving fractures of the material. The two experiments allow a direct comparison with previous work and identify advantages and disadvantages with the continuum and discrete sea-ice dynamics. The long-term goal is to lay the foundation for an SPH sea-ice formulation that can be used in large-scale models.

The paper is organized as follows. In Sect. 2, a description of the sea-ice VP rheology, momentum and continuity equation implementations in the SPH framework is presented. Results of a plastic wave propagation analysis, ridging experiments and ice-arching simulations are presented in Sect. 3. Finally, Sect. 4 discusses the SPH advantages and limitations of the SPH framework, future model development, and the main conclusions from the work.

2 Model

2.1 Momentum and continuity equations

Following Plante et al. (2020), we consider sea ice to behave as a 2D granular material described by the 2D momentum equation (neglecting the Coriolis and sea surface tilt terms):

(1) ρ i h D u D t = σ + τ ,

where ρi is the ice density, h is the mean ice thickness (ice volume over an area), u=ux^+vy^ is the ice velocity vector, σ is the vertically integrated internal stress tensor acting in the y^ direction on a face with a unit-outward-normal pointing in the x^ direction, τ is the sum of water stress and surface air stress, and DDt=t+u is the Lagrangian derivative operator. The Coriolis and sea surface tilt terms are neglected from the momentum equation to ease the comparison with the analytical solution and simple 1D problem. Note that using the Lagrangian derivative operator naturally incorporates the advection of momentum in the ice dynamics – a term that is typically neglected for most continuum-based Eulerian sea-ice models. The surface air stress and the water stress can be written using the bulk formulation as (McPhee1979)


where ρa and ρw are air and water densities, ua and uw are air and water velocity vectors, Ca and Cw are the air and water drag coefficients, and u is neglected in the formulation of the wind stress since uua. The continuity equations for the mean ice thickness h and the ice concentration A can be written as


where the thermodynamic source terms are omitted. Note that the thickness and concentration are independent prognostic variables in a two-category model (Hibler1979), resulting in a singularity when thickness reaches zero. To avoid this singularity and for a more mathematically correct treatment of the mass equation, Gray and Morland (1994) introduced a continuous solution where the concentration asymptotes to zero and one. In the following, we ignore melting and consider cases where only convergent motion is present and the use of a two-category model does not have an impact on the simulated results.

2.2 Constitutive laws

The constitutive relations for the viscous–plastic ice model with an elliptical yield curve, a normal flow rule and tensile strength can be written as (Beatty and Holland2010)


where ϵ˙ij is the symmetric part of the strain-rate tensor, ζ and η are the non-linear bulk and shear viscosities, Pr is the replacement pressure, kt is the tensile strength factor, and δij is the Kronecker delta. Following Bouchat and Tremblay (2017) we write


where P=P*hexp(-C(1-A)) is the ice strength (Hibler1979), P* and C are respectively the ice compressive strength and ice concentration parameters, S is the ice shear strength, and e is the ellipse aspect ratio. In the limit where Δ goes to zero, ζ and η tend to infinity. To avoid this situation, the deformation Δ is capped to Δmin=2×10-9s-1. Using the Δ* formulation, the replacement pressure Pr can be written as

(12) P r = P Δ Δ * ,

which ensures that the stresses are zero when the strain rates are zero.

2.3 Governing differential equations: SPH framework

To solve the system of equations in the SPH framework, equations involving spatial derivatives (Eqs. 1, 4, 5 and 7) are reformulated (see Appendix A for more details on the SPH theory) using Eqs. (A5A6A7) with the particle subscripts p and q (see Fig. A1), and a temporal evolution for the ice particle position is defined:


It is important to make the distinction between the intrinsic ice density ρi and the particle densities ρp. For consistency reasons with the standard VP rheology, we consider the following definition of density independent of ice concentration in contrast to previous work (Wang et al.1998; Ji et al.2005; Staroszczyk2017) (see Results section for discussion):

(18) ρ p = ρ i h p .

By formulating density as in Eq. (18), the continuity Eq. (15) has the same form as the more commonly used continuity density equation (Monaghan2012):

(19) D ρ p D t = - ρ p u p = q = 1 N m q ( u p - u q ) p W pq ,

except for the fact that the divergence of the velocity field is scaled by the ice material density ρi DρpDt=ρiDhpDt. Overall, a particle can be seen as an unresolved collection of floes scattered within the support domain 𝒜 that can converge, ridge over one another, break and drift apart. Note that since the particle density ρp definition is independent of Ap, the concentration can be interpreted as a quantity that measures the compactness of the sea ice at the particle location. It describes the probability of ice floes carried by a particle, which is a point in space, to come in “contact” with ice floes of another particle and get repulsed according to the ice strength.

2.4 Numerical approach

Following Hosseini et al. (2019), we use a second-order predictor–corrector scheme to evolve in time the SPH ice system of equations (see Algorithm 1 below). This integration scheme takes a given function f (here f can be x, u, A and h) and uses a predictor step to calculate its value fn+1/2 at time t=(n+12)Δt (where Δt is the time step) followed by a correction step to calculate the solution fn+1 at time t=(n+1)Δt from fn+1/2:

(20)fpn+1/2=fpn+Δt2DfpnDt+O(Δt2),(21)fp correctedn+1/2=fpn+Δt2Dfpn+1/2Dt,(22)fpn+1=2fp correctedn+1/2-fpn+O(Δt3).

In the above equations, Ot2) and Ot3) represent higher-order terms, which are ignored in the proposed scheme. Following Lemieux and Tremblay (2009), a simple 1D model taking into account only the viscous term – the most restrictive condition – leads to the following stability criterion:

(23) Δ t ρ i h l min 2 η max = e 2 ρ i l min 2 Δ min P * ( 1 + k t ) ,

where lmin is the minimum smoothing length across all the particles. The stability criterion imposes a strict limitation on the time step (10-4 to 10−2 s for particles with a radius of 1 to 10 km); this cannot be avoided using a pseudo-time step because particles in an SPH framework are irregularly placed and move within the domain at each time step. This makes the parallelization of the particle interactions algorithm mandatory for any practical applications. On the positive side, the explicit time stepping also eliminates possible convergence issues of the numerical solver. A pseudo-code for the proposed algorithm is shown below (Algorithm 1).

Algorithm 1Sea-ice SPH.

Domain shape and boundaries, Spatial resolution, Total integration time
initialize particle and boundary according to input
for i=0 to IntegrationTime do
for j=0 to nInteraction do
end for
for all particles do
physicalQuantities (externalForce,internalForce)
end for
timeStep smoothingLength
monitor particle interaction statistics
end for

2.5 Particle interactions

Following Rhoades (1992), we use the bucket search algorithm parallelized using shared memory multiprocessing (OpenMP) to find all the neighbours of each particle in favour of the tested tree algorithm (Cavelan et al.2019), which involve pointers and complex memory structure that are not easy to manipulate in OpenMP. The proposed OpenMP parallelization is rudimentary, and one time step in a domain with 40 000 particles takes ≈0.1 s. For this reason, the model requires more computational resources for the effective resolution when compared with a continuum approach. This could be greatly improved by taking advantage of CPU clusters (Yang et al.2020) or GPUs (Xia and Liang2016).

After the neighbour search, the interactions between pairs of particles are computed using the Wendland C6 kernel – Wendland kernels have the best stability properties for wavelengths smaller than the smoothing kernel (Dehnen and Aly2012) – which is written as


where αd is a normalization factor depending on the dimension of the problem. Note that R (=κ|rp-rq|/lp) is the normalized distance between particles in the referential rprq. Consequently, we always integrate from 0 to lp (the smoothing length) independently of the kernel instead of 0 to κlp as shown by Liu and Liu (2010). The constant αd becomes 78κ27πl2 in 2D, with a factor of κ2 difference from the usual definition. Note that the scaling factor κ has a value of 1 for the Wendland C6 kernel. The choice of kernel was validated using stability tests with six different kernels including the original Gaussian kernel (Gingold and Monaghan1977); a quartic-spline Gaussian approximation (Liu and Liu2010); a quintic-spline Gaussian approximation (Morris et al.1997); a quadratic kernel (Johnson and Beissel1996); and the Wendland C2, C4 and C6 kernels (Wendland1995).

2.6 Smoothing length

The smoothing or correlation length is a key element of SPH and has a direct influence on the accuracy of the solution and the efficiency of the computation. For instance, if lp is too small, there may not be enough particles in the support domain violating the kernel moment requirements. If the smoothing length lp is too large, all the local properties of particles would be smoothed out over a large number of neighbours and the computation time would increase with the number of interactions. In 2D, the optimal number of neighbours interacting with any particle p should be about 20 to balance the precision and the computational cost (Liu and Liu2003). We therefore implement a variable smoothing length that evolves in time and space to maintain this approximate number of neighbours. To this end, we keep the mass of particles constant in time and evaluate the smoothing length from the particle density. Note that keeping the mass of a particle constant has the advantage of ensuring mass conservation. This assumption is justified in our case since we are only interested in sea-ice dynamics, and ridging changes the area covered by ice floes but not their mass. However, fixing the ice mass is only valid when neglecting the thermodynamics and needs to be modified for synoptic-scale simulations.

The initial mass of a particle is defined from the ice area it represents within its support domain (Δ𝒜p in Fig. 1). To avoid creating porosity in the medium, we divide the space in equal square area (=Lp2) that covers the whole domain. Since we want approximately 20 neighbours for every particle, we introduce α (=3 in all simulations), a parameter that stands for the approximate number of particles desired in any direction within the support domain. The parameter α can also be interpreted as the proportionality constant between the particle spacing Lp and the smoothing length lp. Note that to increase the accuracy of the particle approximation, α can be increased by any desired factor (see Fig. 1). The mass carried by a particle is therefore written as

(26) m p = Δ A p ρ i h 0 p = L p 2 ρ i h 0 p ,

where h0p is the initial mean thickness of the particle. The smoothing length is then updated at each time step diagnostically from

(27) l p = α L p = α m p ρ p .

Figure 1Graphical representation of the initial position of the particles and the relevant parameter for the smoothing length evolution: the ice area carried by the particle Δ𝒜p (solid orange square), the parameter α (=2 in this schematic for visibility), the support domain 𝒜 (dashed orange line), the smoothing length lp (red arrow) and the initial distance between particles Lp. Black circles are neighbouring particles q, and the orange circle is the current particle p. Note that, as for Fig. A1, the particle sizes in this schematic are also arbitrary.


The smoothing length lp is capped to 10 times its initial value when the particle density tends to zero. This capping prevents conservation of mass for a density lower than 1 % of its initial value (see Eq. 26). We justify this capping because such small densities do not affect the ice dynamics.

2.7 Boundary treatment

We implemented the boundary treatment of Monaghan and Kajtar (2009) because of its simplicity, versatility and low computational cost. The boundaries are set up by placing stationary particles with fixed smoothing length lb and a mass mb equal to the average ice particle mass mp. The boundary smoothing length lb is chosen in a way that only one layer of ice particles initially interact with the boundary (this makes lb resolution-dependent). The boundary particles are (equally) spaced apart by a factor of 1/4 of their smoothing length (lb/4). In this manner, all ice particles p within a support domain lb will interact with approximately four boundary particles (denoted by the subscript b) at a time, resulting in a net normal repulsive force FNp,

(28) F N p = b = 1 N b μ ( r p - r b ) | r p - r b | 2 W p b 2 m b m p + m b ,

which is added to their momentum equation. In Eq. (28), μ is a constant with units of kg m4 s−2 used to adjust the repulsion strength and is also simulation-dependent because it needs to counterbalance the particle acceleration and prevent them from escaping the domain. This free parameter is not suited for complex pan-arctic simulations but is sufficient in our idealized experiments. For all the simulations, a free-slip boundary condition, i.e., no tangential friction force between boundary particle and ice particle, is applied.

Table 1Physical parameters used in ridging and arch simulations.

Values of the parameters used for the simulations are the same as the ones presented in Williams et al. (2017) to facilitate comparison in the Results section.

Download Print Version | Download XLSX

3 Results

3.1 Plastic wave propagation

We first compare the plastic wave speed for the VP dynamic equations with and without the SPH approximations. To this end, we do a perturbation analysis for a 1D case with a fixed sea-ice concentration (A=1). In this case, the 1D SPH sea-ice dynamic equations (Eqs. 1316) form a system of three equations and three unknowns (x, u and h):


where xpq is a short form for xpxq and

(32) Γ = P * 2 [ ± ( e - 2 + 1 ) 1 / 2 - 1 ] .

In the above, we made use of the following 1D normal stress for convergent plastic motion (see Gray1999; Williams et al.2017, for 1D normal stress derivation):

(33) σ = σ x x = P * 2 [ ± ( e - 2 + 1 ) 1 / 2 - 1 ] h = Γ h .

Linearizing around a mean state (u=0 and h=h0), considering small perturbations (δx, δu and δh) and ignoring the second-order term, we obtain


where ΔAq=mqρih0=mqρq (Eq. A4) and where we have used the binomial expansion 1h=1h0+δh1h0(1-δhh0). Following Williams et al. (2017), we do a perturbation analysis on the system of Eqs. (34)–(36) and assume a wave solution of the form δf=f^exp(i(kx-ωt)), where i is the imaginary number; k is the wavenumber; ω is the angular velocity; and f is a dummy variable standing for u, x and h. Substituting δf in Eqs. (34)–(36), the resulting set of equations in the reference frame following the ice motion reduces to


Note that since the ice is initially at rest, the Lagrangian and the Eulerian frameworks are equivalent. For a large enough wavelength (so that the perturbation can be resolved across multiple particles with high accuracy, i.e., λlp and N→∞), the summations can be approximated by integrals over the space; i.e., q=1NAqxpq|xpq| becomes -dxpq. Taking advantage of the kernel properties – i.e., all moments higher than 0 vanish – we can write Eqs. (38)–(39) as


where the integrals have been converted to Fourier transform using F(Wxpq)=-(Wxpq)exp(-ikxpq)dxpq=ikF(W)=ikW̃. Finally, Eqs. (37), (40) and (41) represent a system of three equations for three unknowns (x^, u^, h^) that we solve by substitution. This leads to the following form for the phase speed of the plastic wave (ωk):

(42) c SPH = ω k = ± W ̃ - Γ ρ i ( 2 W ̃ - 1 ) .

For wavelengths much larger than the smoothing length (λ1klp), the Fourier transform of the kernel tends to 1 (W̃1) and the SPH formulation reduces to the viscous–plastic theory without SPH approximations (see for instance Williams et al.2017), i.e.,

(43) c VP = ± - Γ ρ i ,

with a plastic wave propagation speed cVP≈5.7m s−1 for typical sea-ice parameters (see Table 1). Consequently, a major difference of SPH with the FDM framework is that the plastic wave speed is dispersive with a phase velocity cSPH that is dependent on the wavelength and the smoothing length. In general, only the plastic waves with a wavelength between approximately 1 and 11 times the smoothing length will have their travelling speed modified by more than 1 %. More specifically, in the limit where the wavelength λ approaches the smoothing length lp, the plastic wave speed increases in the SPH framework to a maximum value of ≈6.7m s−1 (see Fig. 2a). Note that for wavelengths smaller than the smoothing length, the summations in Eqs. (40)–(41) cannot be written as integrals, but the particles still respond partially to the perturbations. This sometimes leads to the tensile and the zero-energy mode instabilities (Swegle et al.1995). As mentioned above, Dehnen and Aly (2012) showed that Wendland kernels can diminish the tensile instability and the pairing of particles. A deeper analysis of unresolved waves (λ<lp) in the context of sea-ice SPH dynamic equations is beyond the scope of the current study.

For the more general case when the base state allows for a variable concentration (linearized around a mean state A=A0) and considering the classical – denoted by a superscript C – particle density definition (ρpC=ρihpAp) used by Wang et al. (1998), Ji et al. (2005), and Staroszczyk (2017), the plastic wave speed becomes

(44) c A,SPH C = ± W ̃ - Γ * ρ i ( C A 0 - 3 + 2 W ̃ ) ,

where Γ*=Γexp(-C(1-A0)). We argue that the plastic wave speed cA,SPHC obtained with the classical density definition does not converge (see Fig. 2b) to the viscous–plastic theory, cA,VP, derived from FDM (see Williams et al.2017, for a derivation),

(45) c A,VP = ± - Γ * ρ i ( C A 0 + 1 ) ,

because the ice concentration is taken into account in both the definition of ρpC and implicitly in the definition of the average thickness hp. When we consider the new formulation of particle density independent of concentration as proposed above (Eq. 18), the wave speed equation becomes

(46) c A,SPH = ± W ̃ - Γ * ρ i ( C A 0 - 1 + 2 W ̃ ) ,

which reduces to the FDM VP theory (Eq. 45) when the wavelength is large compared to the smoothing length (see Fig. 2c). Note that the perturbation analysis presented above is not valid for the classical density definition proposed by Wang et al. (1998), Ji et al. (2005), and Staroszczyk (2017) since they use a different set of momentum, continuity and constitutive equations to describe sea ice. In a similar manner to the plastic wave speed with a fixed concentration (Eq. 42), the wave speed cA,SPH (Eq. 46) is dispersive and the wavelengths between 1 and 11 times the smoothing length are those that are mostly affected (more than 1 %). However, in this case, the plastic wave speed is damped for wavelengths close to the smoothing length for a mean concentration state higher than 0.1. Note that while the plastic wave speed is defined for all A, it does not have a physical meaning for A<0.85 since there are negligible ice–ice interactions.

Figure 2SPH plastic wave speed as a function of the normalized wavelength (λ/lp) for the Wendland C6 kernel. Panel (a) shows the classical VP rheology with fixed concentration (Eq. 42) normalized by the FDM plastic wave speed with fixed concentration (Eq. 43), (b) shows the classical VP rheology with a variable concentration and the density definition ρpC=ρihpAp (Eq. 44) normalized by the FDM plastic wave speed with a variable concentration (Eq. 45), and (c) shows the classical VP rheology with a variable concentration and the density definition ρp=ρihp (Eq. 46) normalized by the FDM plastic wave speed with a variable concentration (Eq. 45). Different homogeneous base states of concentration A0 are shown varying from 0 to 1.


Figure 3Idealized domain of the ridging experiment. The blue circles represent the ice particles, and the black ones are the boundary particles. The grey arrow shows the wind forcing. More particles than shown in this schematic were used during the simulation.


Figure 4Temporal evolution of simulated sea-ice thickness along the central horizontal line of the domain for (a) the ridge experiment initialized with a concentration A=1 and average thickness h=1m and (b) the ridge experiment initialized with a concentration A=0.5 and average thickness h=0.5m. The wall is located at x=0, and the wind speed is -5x^m s−1. The theory follows Eq. (47).


3.2 Ridging experiments

We validate our implementation of the SPH model (with the new definition of particle density ρp) in a 1D ridging experiment for which we can validate against the simulated field from a viscous–plastic sea-ice model based on the FDM – the 1D version of McGill-SIM model used in the SIREx studies (Bouchat et al.2022; Hutter et al.2022) – and against the analytical solution (see Williams and Tremblay2018, for a derivation):

(47) - d σ d x = ρ a C a | u a | u a d h d x = 2 ρ a C a | u a | u a P * ( e - 2 + 1 + 1 ) ,

i.e., a linear profile in thickness with a slope proportional to the square of the wind velocity and inversely proportional to the ice strength. We consider a rectangular domain of 1000 by 2000 km including the boundary (the ice field is 1900 km to ensure that no particles escape on the open side) with 37 240 particles; an initial homogeneous smoothing length lp of 21.429 km (spacing lp/α=7.14km); and a smaller – to limit boundary effect – boundary particle smoothing length lb of 4 km (spacing lb/4=1.0km) to represent the wall (see Fig. 3). Particles are initialized with an average thickness h=1m and a concentration A=1. They are forced against the wall by a constant unidirectional wind of 5 m s−1. Note that the water stress is removed in the simulation for a faster convergence to the steady state, which enables higher resolution – a water current of 0 m s−1 would slow down the ice and the ridge formation since it is driven by the advection speed. The Coriolis force should normally also have to be considered with this domain size at classical polar latitude – the Rossby number is 𝒪(10−2) – but is neglected in this idealized experiment to conserve the symmetry of the solution and compare it to the theoretical 1D equation (Eq. 47). In results presented below (Figs. 45), the particle properties are averaged over a grid of approximately 10 by 5 km cells for plotting purposes.

Results show that the simulated thickness field converges to the analytical solution (within an error of ≈1 %) after ≈5 d with a slope of 1.33×10-3m km−1, compared to 1.34×10-3m km−1 for the theory – lower-resolution simulations were run for a longer time and also converged to this stable state (results not shown). This is comparable to the precision obtained by the 1D McGill-SIM FDM model which reaches a slope of 1.35×10-3m km−1. Artifacts are observed close to the boundary where the repulsive force prevents the particles from reaching the “wall”. Additionally, when a particle comes into contact with the boundary with a certain inertia (due to the 1/r dependence of the boundary force), we observe oscillations in the motion of particles which can propagate far in the domain (e.g., Fig. 4a, at x[50,300]km and t=[30,45]h). The oscillations are damped, and the energy is dissipated by the rheology term with time until an equilibrium is reached. Note that reintroducing the water drag diminishes the oscillation coming from the boundary but does not remove them completely. A more physical boundary treatment is beyond the scope of this study.

Figure 5Evolution in time of (a) the thickness normalized by concentration rate of change in time d(h/A)dt, (b) the average thickness h and (c) the concentration A at x=300km. The rate of change in time is computed from dfdt(x,t)=f(x,t+Δt)-f(x,t-Δt)2Δt.


We also repeated the ridge experiment with the same forcing and total sea-ice volume but letting the sea-ice concentration evolve with time. Specifically, the initial average thickness and concentration were set to h=0.5m and A=0.5. This ensures that both h and A covary in time such that hA remains constant – note that A and h follow the same continuity Eqs. (15) and (16) or (4) and (5) when omitting the SPH approximations, and therefore they should vary identically in time until A reaches 1 – in the marginal ice zone (MIZ), which we define as the area where the sea-ice concentration ranges between 0.15 and 0.85 and where low ridging by ice collision occurs (see Fig. 4b). To accomplish this, the domain was extended to 4000 km (3800 km, excluding the boundaries) and the initial particles spacing changed from 7.14 km to 10.0 km for a corresponding initial smoothing length lp of 30.0 km and total number of particles of 38 000. In this configuration, the model converges to a steady-state solution in ≈10 d with a slope of 1.36×10-3m km−1, in agreement with theory within an error of ≈1 % (see Fig. 4b). Results at x=300km, away from boundary effects, show that (as desired) thickness and concentration evolve coherently – d(h/A)dt0 – before ice concentration reaches ≈85 % (see Fig. 5a). At that point (t≈22h), ice–ice interactions emerge and the ridging process starts d(h/A)dt>0. One key difference with the simulation initialized at A=1 is a thickness build-up (above 1 m) at the edge of the ridge in the MIZ. At this location, the continuity equation for sea-ice concentration is capped, while that of the mean ice thickness remains continuous. This results in a local increase in ice thickness to ≈1.1m. This process is akin to the wave radiation drag in the MIZ (Sutherland and Dumont2018). A detailed analysis of simulations in simple convergent ice flow in the MIZ with ice concentration close to 100 % will be considered in future work.

In the ridge building phase, the speed of advance of the ridge front increases until a maximum concentration is reached after ≈70h (see Fig. 5c). Subsequently, the ice drift speed reduces, and the rate of advance of the ridge slows down. When the ice thickness gradient is in balance with the surface wind stress (after ≈200h), d(h/A)dt reaches a steady state. Overall, we can observe three stages in the ridge formation (see Fig. 5): first, a rapid compaction stage when ice particles are drifting close to their free-drift speed since the ice strength is weak; second, a transition stage between A≈0.85 and 1.00 when ridging occurs in the MIZ analogous to the wave radiation drag mentioned above; and third, a ridging stage with changes in ice thickness that are about 1 order of magnitude higher than during the transition stage. Note that the amplitude of oscillations between particles within the domain or at the boundaries in ridging experiments diminishes when incorporating the water drag (a damping term). The water drag also increases the time needed to reach steady state because the ice drift speed is slower.

3.3 Arch experiments

We next compare the SPH approach with the FDM and DEM sea-ice models in a second well-studied idealized experiment: the ice arch formation. To this end, we run the SPH model in an idealized domain representing the Nares Strait (see Fig. 6) with an upstream reservoir 5 times the length of the channel (L) to minimize the boundary confinement effect without sacrificing the spatial resolution.

The set of simulations uses a domain with L=60km. The initial conditions for ice thickness, concentration and velocity are h=1m, A=1 and u=0m s−1. The ice is forced with a constant unidirectional wind of −7.5m s−1 in the y^ direction, and ocean current is fixed to uw=0m s−1. The corresponding surface stress is ≈0.04kN m−2, and the total integrated stress at the entry of the channel is slightly smaller than P* (05Lτadx=26.325kN m−1). We use a weaker wind than commonly used in Nares Strait ice arch simulations (≈10m s−1) to limit the ridging phase prior to the formation of the ice arch. In this experiment, we limit ourselves to ice with no tensile strength (kt=0) and a shear strength of 6.875 kN m−2, i.e., an ellipse aspect ratio of 2.

We suspect that the SPH and DEM frameworks have a similar behaviour in certain circumstances even though they have different (implicit) rheologies because of their Lagrangian nature. Indeed, the interpretation of the numerical representation of a particle in SPH as a collection of ice floes is also present in the DEM (Li et al.2014), and the two numerical frameworks compute their quantities with one-to-one interactions. Consequently, we first test whether the SPH approach has the same sensitivity to the relative size of particles with respect to the channel width as in the DEM (Damsgaard et al.2018). Results showed that no stable arch can be formed with the specified forcing for all particle diameter sizes tested (7.5, 5, 3.75 km) (see ice velocity field in Fig. 7). Instead, a “continuous” slow flow of ice is present in the channel. The discontinuity at the entry of the channel is visible in the concentration, thickness and velocity fields (Fig. 7) and can be interpreted as an intermittent (unstable) ice arch formation. Also, we noted that larger particles are not more prone to ice jam than smaller ones. This is contrary to what is known from granular material theory and to results from Damsgaard et al. (2018) that show a transition from stable to no ice arch formation for floe sizes ranging from approximately 1/4 to 1/16 of the strait's width. We explain this difference between SPH and DEM from the continuum description of the ice dynamics equation. In the present model, the constitutive laws prescribe the repulsion of the particles from one another according to the ice strength, which is a function dependent on the ice concentration and mean thickness, not on the particle size. We conclude that to enforce granularity within the SPH framework, the constitutive laws would need to be adapted to account for contact force and particle size which could then reproduce a similar behaviour to that observed in the DEM. However, even though the increase in resolution – or particle size – has no effect on the arch stability, it enables smaller fracture resolutions that are visible at the entrance of the channel (see ϵI and ϵII in Fig. 8). In our SPH model, the stress invariants σI and σII shows oscillation patterns in regions where the ice is in the viscous regime (see the tree-like structure in the normal and shear stress fields in Fig. 8). From our experiments, the “tree-like” peak stresses appear during the transient phase and at steady state. However, the particles never stop moving even in a steady state because viscous deformations are always present. We hypothesized that stress patterns are associated with over-damped elastic waves associated with small movements (but large internal stresses) of the particle in the viscous regime. Those structures are not symmetric, despite symmetrical initial conditions, because of the domino effect between interacting viscous waves. Note that they are absent from the strain-rate fields since viscous deformations are extremely small. They are also absent in sea-ice model based on a continuum approach (Dumont et al.2009; Dansereau et al.2017; Plante et al.2020), but these tree-like structures are qualitatively similar to the stress structure between floes observed in the DEM (e.g., Damsgaard et al.2018, Fig. 5c). Despite the fact that the model solves the same continuum equations as other FDM models, we believe that stress networks can be observed with the SPH method because particles interact in a pairwise fashion according to their relative distance. This can create less-dense ice areas within the medium which can lead to oscillations in the stress field. It is known that SPH can have spurious behaviour in some cases when the stress is solved at the same location as the particle centre (as done here). This can be avoided using stress particles (see Chalk et al.2020, for details). More investigations are required to test whether this behaviour is physical. This is left for future work.

Figure 6Idealized domain of the ice arch experiments. The blue circles represent the ice particles, and the black ones are the boundary particles. The grey arrow shows the wind forcing.


Figure 7Ice concentration, thickness and total velocity (h, A, |ui|) at time t=24h for an initial particle spacing of (a) 7.5, (b) 5 and (c) 3.75 km (8, 12 and 16 particles can fit in the strait respectively), and the initial total integrated surface stress at the entry of the channel is 26.325 kN m−1.


Figure 8Strain rate and stress invariants (ϵ˙I, ϵ˙II, σI, σII) at time t=24h for an initial particle spacing of (a) 7.5, (b) 5 and (c) 3.75 km (8, 12 and 16 particles can fit in the strait respectively), and the initial total integrated surface stress at the entry of the channel is 26.325 kN m−1.


Second, we explored the ability of the model to produce stable ice arches. To this end, we reduce the total integrated surface stress at the entry of the channel to 13.146 kN m−1 (or wind speed of 5.3 m s−1) to a value below the ice compressive strength (P*) to avoid completely ridging the north of the channel and jump immediately to the arch-forming stage. In this case, the results show a clear stable arch (see Fig. 9) with a shape that is qualitatively similar to the one presented by Dansereau et al. (2017), Plante et al. (2020), and West et al. (2022). The formation of a stable arch in an SPH model is possible with the standard shear strength (e=2), in contrast to continuum models that require an increase in shear strength (e.g., Dumont et al.2009; Dansereau et al.2017; Plante et al.2020) – it is important to keep in mind that the domain configurations were different in each of those studies. This suggests that SPH has a different sensitivity of ice arching to the ellipse aspect ratio e and ice thickness h. With a no-slip boundary condition and the same default yield curve (same P* and ellipse aspect ratio e), preliminary results suggest that no arches form – the pack is undeformed – and instead a higher surface wind stress is required to form an arch. Note that in the SPH simulations, only one arch forms close to the outlet. Presumably, the number of arches would increase and location would change if the model was run at higher resolution, with different boundary conditions or in a non-idealized domain geometry. Overall, this shows that SPH is able to capture large-scale features coming from small-scale interactions. The simulation of a stable ice arch (Fig. 9) also shows how SPH can fracture and create discontinuities in the ice field as seen in DEM models. This behaviour is similar to that of the elastic-decohesive sea-ice constitutive model of Schreyer et al. (2006) or the FEM model of Rampal et al. (2016). Finally, in the SPH framework, a lead or polynya can be defined by an absence of particles for leads larger than particle size – akin to DEM – or by particles with reduced concentration for sub-particle size leads – akin to FDM.

Figure 9Thickness field at time t=0, 48 and 168 h for an initial particle spacing of 7.5 km and a total integrated stress at the entry of the channel of 13.146 kN m−1.


4 Discussion and conclusion

In this paper, we have presented a first implementation of the viscous–plastic rheology with an elliptical yield curve and normal flow rule in the framework of SPH with the long-term goal of simulating synoptic-scale sea-ice dynamics. We have described the basics of the SPH approach and how the sea-ice dynamic equations can be formulated in this framework along with the implementation of key components of the numerical method such as the smoothing length, the kernel, the boundaries and the time integration technique. We proposed a different definition of the particle density and showed that the more commonly used density definition involving the ice concentration (Wang et al.1998; Ji et al.2005; Staroszczyk2017) when used together with the average ice thickness leads to erroneous plastic wave speed propagation. A particle density definition independent of the ice concentration corrects this and leads to results that are in line with the VP theory. The SPH model thus developed is in excellent agreement (error of ≈1 %) with an analytical solution of the VP ice dynamics for a simple 1D ridging experiment. The approximations used at the core of the SPH framework result in a dispersive plastic wave speed in the medium – in contrast to its FDM counterpart – which is dependent on the smoothing length (or resolution) and the choice of the kernel. The plastic wave speed is mostly affected for wavelengths 11 times the smoothing length and lower.

From the simple ridging experiment with fixed sea-ice concentration (A=1), we observe nonphysical damped oscillations that propagate in the domain associated with our choice of boundary conditions. The conclusions drawn from our simulations are robust to the choice of boundary conditions. Nevertheless, this behaviour needs to be removed for a proper simulation of sea ice near coastlines. The ridging experiment with an initial ice concentration below 100 % showed that continuity equations for concentration and thickness evolve coherently until a concentration of 85 %. At that point, SPH particles start to ridge locally in the MIZ in addition to the wall where the maximum stress is located. This effect is not observed in the continuum approach and is presumably related to particle collisions in converging motion.

When compared to other numerical frameworks, the SPH model is able to reproduce stable ice arches in an idealized domain of a strait with an ellipse aspect ratio of 2 and a wind forcing of 5.3 m s−1, in contrast to other continuum approaches that require higher material shear strength. However, when using a stronger wind field of 7.5 m s−1, no stable arches are formed when increasing the particles size in the strait (stable arches are only achieved when increasing particle average thickness). We concluded that the number of particles in the strait does not influence the formation of ice arches, in contrast to the DEM, and is analogous to an increase in resolution in a continuum framework: a larger number of particles influence the number of fractures that can form and the resolution of fine-scale structures. The stress fields produced by the SPH model in the channel experiment show a tree-like pattern upstream of the channel where there are low total deformations. This is not observed in FDM experiments, but it is qualitatively similar to the tensile stress network exhibited in DEM simulations (Damsgaard et al.2018) that comes from individual contact force between the ice floes and is hypothesized to be associated with damped viscous sound waves.

Even though we successfully implemented the standard sea-ice viscous–plastic rheology with an elliptical yield curve and a normal flow rule in an SPH framework, the current model does not outperform a classical FDM model. In fact, there are inherent difficulties and instabilities in SPH that do not exist in FDM. It is known that the SPH framework trades consistency – i.e., the ability to correctly represent a differential equation in the limit of an infinite number of points with a null spacing between them – for stability, which gives the SPH a distinct feature of working well for many complicated problems with good efficiency but less accuracy. However, the classical formulation of SPH used and described in the present work does not usually respect zeroth-order consistency because of the unstructured particle position in space (see Belytschko et al.1998, Sect. 3 for a derivation). Nevertheless, consistency can be improved at the expense of computational cost (Chen and Beraun2000; Liu et al.2003) by reformulating the SPH core approximation (Eq. A1). Also, the boundary description has been identified as a weak point of the SPH framework. Prescribing Dirichlet, Neumann or Robin boundary conditions is not as straightforward as in continuum approaches. Moreover, preventing particle penetration through a boundary is still a challenging task (Liu and Liu2010), and the SPH consistency is usually at its worst at the boundary because the support domain is truncated. In the present study, a proper physical representation of the boundary was not adopted, and the boundary treatment was chosen for its numerical simplicity and should be modified in future work. Other major issues with SPH are the zero-energy modes and the tensile instability previously mentioned. The zero-energy modes can be found in FDM and FEM, and they correspond to modes at which the strain energy calculated is erroneously zero (Swegle et al.1995). The tensile instability results in particle clumping or nonphysical fractures in the material. In the present work, we adopted a different kernel from the usual Gaussian spline to avoid those instabilities, but other methods such as the independent stress point (Dyka and Ingel1995; Chalk et al.2020), artificial short-length repulsive force (Monaghan2000), particle repositioning (Sun et al.2018) or adaptive kernel (Lahiri et al.2020) can be used if more stabilization is needed. For example, at a smaller scale, SPH simulation of ice in uniaxial compression was improved by a simplified finite-difference interpolation scheme (Zhang et al.2017). More specifically for sea-ice models, the pressure closure of Kreyscher et al. (2000) is not well-suited for long simulation. Indeed, particles can still move when they are in the viscous state but have low internal ice pressure because of the replacement pressure scheme. Consequently, particles could pass through each other resulting in erroneous locations of the parameters carried out. Finally, using SPH for sea-ice modules in grid-based continuum global climate models (GCMs) complicates the coupling with ocean and atmosphere components since particle quantities need to be converted on a grid and vice versa.

In its current state, the model reproduces very similar behaviour to other FDM continuum models and does not constitute a large improvement. Nevertheless, we believe that SPH enables the possibility to describe sea ice as a continuum at large scale using what is already known from continuum models and to explore some new avenues at small scales, where the continuity approximation is questionable. Indeed, SPH also has interesting properties that could be exploited. For example, SPH can be used with little change for problems involving several fluids, whether liquid, gas or dust fluids (Monaghan2012). This feature could be exploited in the creation of a general approach for all components of a GCM (atmosphere, ocean and sea ice). The method developed is also a proper option for nowcasting sea-ice prediction because only the ice dynamics need to be considered in nowcasting applications and the model has a good ability to carry the ice property in space. SPH can fracture and transition from the continuum to fragments seamlessly since it is not restricted on a grid, which also has the advantage of enabling ice edge shapes independent of it. The ability of SPH to move around particles has the interesting property of concentrating them in converging motion, effectively increasing the spatial resolution of the model in regions under high-stress activity and dispersing particles when the flow is divergent, which decreases the resolution in low-ice-concentration areas. This property should result in higher accuracy than that in typical continuum models. The elastic behaviour assumed for sea ice with a certain rheology can be associated with the weak compressibility inherent in the classical formulation of SPH. Finally, the SPH discretization of the continuum into particles enables the implementation of several new features. For example, angular momentum to individual floes (or pack of floes) can be added to take into account rotation along LKFs. A direct measure of the concentration from the number of particles within a support domain (this takes advantage of the already-computed number of neighbours and helps to ensure the desired number of neighbours in converging flow) can be computed. A subscale parametrization of floe–floe contact force (this short-length repulsive force could also help with the tensile instability) can be implemented. A varying floe size distribution could be incorporated by varying the mass carried by a particle for a given particle density.

For future work and before exploring new features enabled by the SPH numerical framework, a more physical treatment of the boundary conditions should be investigated to properly simulate the grounding of sea ice near the coast enabling the no-slip conditions. Subsequently, the model could be tested against other benchmark problems in an idealized domain to further understand and compare the effect of the SPH method (Flato1993; Hunke2001; Hibler et al.2006; Danilov et al.2015; Mehlmann et al.2021). Also, in order to use the model for pan-Arctic simulations, the Coriolis and sea surface tilt force along with the treatment of the thermodynamics source and sink terms should be implemented in the SPH framework (see preliminary work by Staroszczyk2018). In addition, the parallelization of the code should be improved in order to bring the computational time down to a value comparable to that of an FDM model. Finally, while there still is a significant amount of work to be completed before SPH can be used in large-scale climate simulations, the method shows promise and deserves further investigations and development.

Appendix A: Smoothed particle hydrodynamics basics

The SPH method is at the interface between the finite-element and discrete-element methods. In this framework any function f(r) at a point r is approximated from neighbouring values in the parameter space f(r) using an integral interpolant (see Fig. A1):

(A1) f ( r ) = V f ( r ) W ( | r - r | , l ) d r ,

where W(|r-r|,l) is the interpolating kernel and 𝒱 is the entire space volume. In 2D, the space volume is an area 𝒜, and the kernel has units of per square meter (m−2). This integral interpolant approximation is based on the singular integral mathematical framework of Natanson (1961) and imposes the following restrictions on the kernel:

(A2) A W ( | r - r | , l ) d r = 1


(A3) lim l 0 W ( | r - r | , l ) = δ ( r - r ) ,

where l is the smoothing length of the kernel and δ is the Dirac delta function. Using the particle approximation, Eq. (A1) can be written as a weighted summation over all neighbouring points within the area 𝒜:

(A4) f ( r p ) q = 1 N f ( r q ) W ( | r p - r q | , l p ) Δ A q q = 1 N f ( r q ) W ( | r p - r q | , l p ) m q ρ q ,

where N is the number of points in space referred to as neighbour particles, Δ𝒜q (=m/ρ) is the area associated with the particle q, m represent the mass (kg) and ρ is the 2D density (kg m−2).

Figure A1Graphical representation of the SPH kernel W(|rp-rq|,lp) (solid orange line), the smoothing length lp (red arrow), the particle p, the neighbouring particles q, the support domain 𝒜 (dashed orange line), and the distance between any neighbour particle q and the particle p within the support domain rprq (black arrow). Note that particles are points in space and that their size in this schematic is arbitrary.


From the above approximations, we reformulate differential operators relevant to our study in their discrete SPH forms. We write the divergence of a vector field (V), the divergence of a tensor (T) and the gradient of a vector field (V) (Monaghan2005) as (see Appendix B for complete derivation)


In Eq. (A7), denotes the outer product. pWpq is the gradient of the kernel at the coordinate rprq in the reference frame of particle p and is written as

(A8) p W pq = r p - r q | r p - r q | W ( | r p - r q | , l p ) | r p - r q | .

Note that Wpq is a scalar function and consequently pWpq is a vector, the inner product in Eq. (A5) is a scalar, the inner product in Eq. (A6) is a 2D vector and the outer product in Eq. (A7) is a 2D tensor of rank 2. In addition to Eqs. (A2)–(A3), the smoothing kernel must have the following set of properties to avoid nonphysical behaviour and costly computation (Liu and Liu2003):

(A9)compact support:W(|rp-rq|,lp)=0,  for  |rp-rq|>lp,(A10)positive definite:W(|rp-rq|,lp)0,(A11)monotonically decreasing:W(|rp-rq|,lp)(|rp-rq|)0,(A12)symmetric:W(|rp-rq|,lp)=W(-|rp-rq|,lp),(A13)differentiable:nW(|rp-rq|,lp)(|rp-rq|)n,

where stands for exist. In the above, “differentiable” means that the kernel derivatives exist up to the highest order present in the equations. Finally, to ensure the consistency of the discretization of partial differential equations (PDEs; as defined in Belytschko et al.1998) of the SPH method approximations to the nth order, all kernel moments of order 1 to n need to vanish. In practice, the consistency conditions are satisfied when the number of neighbouring particles is sufficiently large to be evenly distributed in the domain of influence (Fraga Filho2019). Note that, at the boundaries, the domain of influence of the particle is truncated, making it impossible to satisfy the kernel moment equations. This phenomenon is referred to as the particle inconsistency and leads to poorer approximations of physical properties. No clear solutions to this problem have been proposed in the literature yet.

Appendix B: Vector operators in SPH

Vector operators take different forms in the SPH framework because they only operate on the smoothing kernel W and they need to ensure symmetric interactions between particles. The following subsections show the demonstrations to obtain the relevant one to our study.

B1 Divergence of a vector

First, the divergence of vector needs to be changed into a form that can be symmetrized. To do so, we use the identity of the divergence of a scalar function times a vector and chose the scalar function to be the density as follow:

(B1) V = 1 ρ ( ( ρ V ) - V ρ ) .

Now applying the integral interpolant approximation (A1) to the divergence term (∇⋅(ρV)) and to the density (ρ) gives


In the above equations, the prime quantities represent the surrounding values. Note that the kernel is the only function that depends on both primed and non-primed positions as defined in Eq. (A1). Using the divergence theorem, the first term in Eq. (B2) can be cancelled:

(B4) V ( ρ V W ) d r = S ( ρ V W ) d s = 0 ,

since the integration surface 𝒮 encompassing the volume 𝒱 is arbitrary and the kernel W has the compact support property (Eq. A9). Applying the particle approximation (A4) to Eqs. (B2)–(B3), we obtain


where we used the identity p=-q and p and q represent the current particle and neighbour. Finally, substituting the last two Eqs. (B5)–(B6) into Eq. (B1) gives the desired form of the operator:


B2 Divergence of a 2D tensor field

Note that in the following demonstration, the Einstein summation convention is used to simplify the calculation and the tensor representation. We start with the divergence of a 2D tensor divided by the density:

(B9) x i ( T i j ρ ) = 1 ρ T i j x i - T i j ρ 2 ρ x i .

Reorganizing the terms gives

(B10) T i j x i = ρ [ x i ( T i j ρ ) + T i j ρ 2 ρ x i ] .

Now applying the interpolant approximation (A1) to the first term in the bracket leads to


As for the divergence of a vector demonstration (Sect. B1), the first integral above vanishes by using the divergence theorem, and applying the particle approximation gives

(B13) ( x i ( T i j ρ ) ) p = - q ( m q ( T i j ) q ρ q 2 ) W pq ( x i ) q = q ( m q ( T i j ) q ρ q 2 ) W pq ( x i ) p .

Substituting this into Eq. (B10) and using the equality Eq. (B6) we get the following expression:


which is the form presented in Eq. (A6).

B3 Gradient of a vector field

To demonstrate Eq. (A7) we first write

(B17) ( a V ) = a V + V a .

Choosing a=1 and recalling that the zeroth-order moment of the kernel also equals 1, we can substitute it in the last term of the expression (B17) and obtain


Finally using the particle approximation (A4) we get


which is Eq. (A7) and where Einstein summation convention was once again used to simplify the derivation.

Code availability

Our FORTRAN SPH sea-ice model code is public and can be found at (last access: 28 February 2024) and (Marquis2024).

Data availability

Output data from the SPH sea-ice model simulations along with a version of the model used and the analyzing programs are available at (Marquis et al.2022).

Author contributions

OM coded the model, ran all the simulations, analyzed results and led the writing of the manuscript. BT participated in weekly discussions during the course of the work and edited the manuscript. JFL and MI participated in monthly discussions during the course of the work and edited the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


Oreste Marquis is grateful for the support from McGill University, Québec-Océan and Arctrain Canada.

Financial support

This project was partially funded by grants and contributions from the Natural Science and Engineering Research Council Discovery Program (ONR-N00014-11-1-0977) and the National Research Council (A-0038111) awarded to Bruno Tremblay.

Review statement

This paper was edited by Jari Haapala and reviewed by three anonymous referees.


Adcroft, A., Anderson, W., Balaji, V., Blanton, C., Bushuk, M., Dufour, C. O., Dunne, J. P., Griffies, S. M., Hallberg, R., Harrison, M. J., Held, I. M., Jansen, M. F., John, J. G., Krasting, J. P., Langenhorst, A. R., Legg, S., Liang, Z., McHugh, C., Radhakrishnan, A., Reichl, B. G., Rosati, T., Samuels, B. L., Shao, A., Stouffer, R., Winton, M., Wittenberg, A. T., Xiang, B., Zadeh, N., and Zhang, R.: The GFDL Global Ocean and Sea Ice Model OM4.0: Model Description and Simulation Features, J. Adv. Model. Earth Sy., 11, 3167–3211,, 2019. a

Beatty, C. and Holland, D.: Modeling landfast sea ice by adding tensile strength, J. Phys. Oceanogr., 40, 185–198,, 2010. a

Belytschko, T., Krongauz, Y., Dolbow, J., and Gerlach, C.: On the completeness of meshfree particle methods, Int. J. Numer. Meth. Eng., 43, 785–819,<785::aid-nme420>;2-9, 1998. a, b

Bouchat, A. and Tremblay, B.: Using sea-ice deformation fields to constrain the mechanical strength parameters of geophysical sea ice, J. Geophys. Res.-Oceans, 122, 5802–5825,, 2017. a

Bouchat, A., Hutter, N., Chanut, J., Dupont, F., Dukhovskoy, D., Garric, G., Lee, Y. J., Lemieux, J.-F., Lique, C., Losch, M., Maslowski, W., Myers, P. G., Ólason, E., Rampal, P., Rasmussen, T., Talandier, C., Tremblay, B., and Wang, Q.: Sea Ice Rheology Experiment (SIREx): 1. Scaling and Statistical Properties of Sea-Ice Deformation Fields, J. Geophys. Res.-Oceans, 127, e2021JC017667,, 2022. a, b

Cavelan, A., Cabezón, R. M., Korndorfer, J. H. M., and Ciorba, F. M.: Finding Neighbors in a Forest: A b-tree for Smoothed Particle Hydrodynamics Simulations, arXiv [preprint],, 7 October 2019. a

Chalk, C., Pastor, M., Peakall, J., Borman, D., Sleigh, P., Murphy, W., and Fuentes, R.: Stress-Particle Smoothed Particle Hydrodynamics: An application to the failure and post-failure behaviour of slopes, Comput. Meth. Appl. M., 366, 113034,, 2020. a, b

Chen, J. and Beraun, J.: A generalized smoothed particle hydrodynamics method for nonlinear dynamic problems, Comput. Meth. Appl. M., 190, 225–239,, 2000. a

Coon, M., Kwok, R., Levy, G., Pruis, M., Schreyer, H., and Sulsky, D.: Arctic Ice Dynamics Joint Experiment (AIDJEX) assumptions revisited and found inadequate, J. Geophys. Res., 112, C11S90,, 2007. a

Damsgaard, A., Adcroft, A., and Sergienko, O.: Application of Discrete Element Methods to Approximate Sea Ice Dynamics, J. Adv. Model. Earth Sy., 10, 2228–2244,, 2018. a, b, c, d, e, f

Danilov, S., Wang, Q., Timmermann, R., Iakovlev, N., Sidorenko, D., Kimmritz, M., Jung, T., and Schröter, J.: Finite-Element Sea Ice Model (FESIM), version 2, Geosci. Model Dev., 8, 1747–1761,, 2015. a

Dansereau, V., Weiss, J., Saramito, P., and Lattes, P.: A Maxwell elasto-brittle rheology for sea ice modelling, The Cryosphere, 10, 1339–1359,, 2016. a

Dansereau, V., Weiss, J., Saramito, P., Lattes, P., and Coche, E.: Ice bridges and ridges in the Maxwell-EB sea ice rheology, The Cryosphere, 11, 2033–2058,, 2017. a, b, c, d

Dehnen, W. and Aly, H.: Improving convergence in smoothed particle hydrodynamics simulations without pairing instability, Mon. Not. R. Astron. Soc., 425, 1068–1082,, 2012. a, b

Dumont, D., Gratton, Y., and Arbetter, T. E.: Modeling the Dynamics of the North Water Polynya Ice Bridge, J. Phys. Oceanogr., 39, 1448–1461,, 2009. a, b, c

Dyka, C. and Ingel, R.: An approach for tension instability in smoothed particle hydrodynamics (SPH), Comput. Struct., 57, 573–580,, 1995. a

Flato, G. M.: A particle-in-cell sea-ice model, Atmos. Ocean, 31, 339–358,, 1993. a

Fleissner, F., Gaugele, T., and Eberhard, P.: Applications of the discrete element method in mechanical engineering, Multibody Syst. Dyn., 18, 81–94,, 2007. a

Fraga Filho, C. A.: Smoothed Particle Hydrodynamics: Fundamentals and Basic Applications in Continuum Mechanics, Springer, ISBN 978-3-030-00772-0,, 2019. a

Gingold, R. A. and Monaghan, J. J.: Smoothed particle hydrodynamics: theory and application to non-spherical stars, Mon. Not. R. Astron. Soc., 181, 375–389,, 1977. a, b

Girard, L., Bouillon, S., Weiss, J., Amitrano, D., Fichefet, T., and Legat, V.: A new modeling framework for sea-ice mechanics based on elasto-brittle rheology, Ann. Glaciol., 52, 123–132,, 2011. a

Gray, J. and Morland, L.: A Two-Dimensional Model for the Dynamics of Sea Ice, Philos. T. Roy. Soc. B, 347, 219–290,, 1994. a, b

Gray, J. M. N. T.: Loss of Hyperbolicity and Ill-posedness of the Viscous–Plastic Sea Ice Rheology in Uniaxial Divergent Flow, J. Phys. Oceanogr., 29, 2920–2929,<2920:lohaip>;2, 1999. a

Gutfraind, R. and Savage, S. B.: Smoothed Particle Hydrodynamics for the Simulation of Broken-Ice Fields: Mohr–Coulomb-Type Rheology and Frictional Boundary Conditions, J. Comput. Phys., 134, 203–215,, 1997. a, b

Herman, A.: Discrete-Element bonded-particle Sea Ice model DESIgn, version 1.3a – model description and implementation, Geosci. Model Dev., 9, 1219–1241,, 2016. a

Hibler, W., Hutchings, J., and Ip, C.: sea-ice arching and multiple flow States of Arctic pack ice, Ann. Glaciol., 44, 339–344,, 2006. a

Hibler, W. D.: A Dynamic Thermodynamic Sea Ice Model, J. Phys. Oceanogr., 9, 815–846,<0815:adtsim>;2, 1979. a, b, c, d, e

Hopkins, M. A. and Thorndike, A. S.: Floe formation in Arctic sea ice, J. Geophys. Res., 111, C11S23,, 2006. a

Hosseini, K., Omidvar, P., Kheirkhahan, M., and Farzin, S.: Smoothed particle hydrodynamics for the interaction of Newtonian and non-Newtonian fluids using the μ(I) model, Powder Technol., 351, 325–337,, 2019. a

Hunke, E. C.: Viscous–Plastic Sea Ice Dynamics with the EVP Model: Linearization Issues, J. Comput. Phys., 170, 18–38,, 2001. a

Hunke, E. C. and Dukowicz, J. K.: An Elastic–Viscous–Plastic Model for Sea Ice Dynamics, J. Phys. Oceanogr., 27, 1849–1867,<1849:aevpmf>;2, 1997. a

Hutchings, J. K., Jasak, H., and Laxon, S. W.: A strength implicit correction scheme for the viscous-plastic sea ice model, Ocean Model., 7, 111–133,, 2004. a

Hutter, N., Bouchat, A., Dupont, F., Dukhovskoy, D., Koldunov, N., Lee, Y. J., Lemieux, J.-F., Lique, C., Losch, M., Maslowski, W., Myers, P. G., Ólason, E., Rampal, P., Rasmussen, T., Talandier, C., Tremblay, B., and Wang, Q.: Sea Ice Rheology Experiment (SIREx): 2. Evaluating Linear Kinematic Features in High-Resolution Sea Ice Simulations, J. Geophys. Res.-Oceans, 127,, 2022. a, b, c

Ji, S., Shen, H., Wang, Z., Shen, H., and Yue, Q.: A viscoelastic-plastic constitutive model with Mohr-Coulomb yielding criterion for sea ice dynamics, Acta Oceanol. Sin., 24, 54–65, 2005. a, b, c, d, e

Johnson, G. R. and Beissel, S. R.: Normalized smoothing functions for sph impact computations, Int. J. Numer. Meth. Eng., 39, 2725–2741,<2725::AID-NME973>3.0.CO;2-9, 1996. a

Kreyscher, M., Harder, M., Lemke, P., and Flato, G. M.: Results of the Sea Ice Model Intercomparison Project: Evaluation of sea ice rheology schemes for use in climate simulations, J. Geophys. Res.-Oceans, 105, 11299–11320,, 2000. a

Lahiri, S. K., Bhattacharya, K., Shaw, A., and Ramachandra, L. S.: A stable SPH with adaptive B-spline kernel, ArXiv [preprint],, 4 January 2020. a

Lemieux, J.-F. and Tremblay, B.: Numerical convergence of viscous-plastic sea ice models, J. Geophys. Res., 114, C05009,, 2009. a

Li, B., Li, H., Liu, Y., Wang, A., and Ji, S.: A modified discrete element model for sea ice dynamics, Acta Oceanol. Sin., 33, 56–63,, 2014. a

Lilja, V.-P., Polojärvi, A., Tuhkuri, J., and Paavilainen, J.: Finite-discrete element modelling of sea ice sheet fracture, Int. J. Solids Struct., 217-218, 228–258,, 2021. a

Lipscomb, W. H., Hunke, E. C., Maslowski, W., and Jakacki, J.: Ridging, strength, and stability in high-resolution sea ice models, J. Geophys. Res., 112, C03S91,, 2007. a

Liu, G. and Liu, M.: Smoothed Particle Hydrodynamics: A Meshfree Particle Method, World Scientific, ISBN 9789812384560, 2003. a, b

Liu, M. B. and Liu, G. R.: Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments, Arch. Comput. Method. E., 17, 25–76,, 2010. a, b, c, d, e

Liu, M. B., Liu, G. R., and Lam, K. Y.: A one-dimensional meshfree particle formulation for simulating shock waves, Shock Waves, 13, 201–211,, 2003. a

Losch, M., Menemenlis, D., Campin, J.-M., Heimbach, P., and Hill, C.: On the formulation of sea-ice models. Part 1: Effects of different solver implementations and parameterizations, Ocean Model., 33, 129–144,, 2010. a

Lucy, L. B.: A numerical approach to the testing of the fission hypothesis, Astron. J., 82, 1013,, 1977. a

Marquis, O.: McGill-sea-ice/SIMP: Sea Ice Modelling Particles (v1.0.0), Zenodo [code],, 2024. a

Marquis, O., Tremblay, B., Lemieux, J.-F., and Islam, M.: Smoothed Particle Hydrodynamics Implementation of the Standard Viscous-Plastic Sea-Ice Model and Validation in Simple Idealized Experiments, Zenodo [data set],, 2022. a

McPhee, M. G.: The Effect of the Oceanic Boundary Layer on the Mean Drift of Pack Ice: Application of a Simple Model, J. Phys. Oceanogr., 9, 388–400,<0388:teotob>;2, 1979. a

Mehlmann, C., Danilov, S., Losch, M., Lemieux, J. F., Hutter, N., Richter, T., Blain, P., Hunke, E. C., and Korn, P.: Simulating Linear Kinematic Features in Viscous-Plastic Sea Ice Models on Quadrilateral and Triangular Grids With Different Variable Staggering, J. Adv. Model. Earth Sy., 13, e2021MS002523,, 2021. a, b

Monaghan, J.: SPH without a Tensile Instability, J. Comput. Phys., 159, 290–311,, 2000. a

Monaghan, J.: Smoothed Particle Hydrodynamics and Its Diverse Applications, Annu. Rev. Fluid Mech., 44, 323–346,, 2012. a, b

Monaghan, J. and Kajtar, J.: SPH particle boundary forces for arbitrary boundaries, Comput. Phys. Commun., 180, 1811–1820,, 2009. a

Monaghan, J. J.: Smoothed particle hydrodynamics, Rep. Prog. Phys., 68, 1703–1759,, 2005. a

Morland, L. W. and Staroszczyk, R.: A material coordinate treatment of the sea–ice dynamics equations, P. Roy. Soc. Lond. A, 454, 2819–2857,, 1998. a

Morris, J. P., Fox, P. J., and Zhu, Y.: Modeling Low Reynolds Number Incompressible Flows Using SPH, J. Computa. Phys., 136, 214–226,, 1997. a

Natanson, I. P.: Theory of Functions of a Real Variable, vol. 2, Frederick Ungar Publishing Co., New York, ISBN 9780486806433, 1961. a

Overland, J. E., McNutt, S. L., Salo, S., Groves, J., and Li, S.: Arctic sea ice as a granular plastic, J. Geophys. Res.-Oceans, 103, 21845–21867,, 1998. a

Peiró, J. and Sherwin, S.: Finite difference, finite element and finite volume methods for partial differential equations, in: Handbook of materials modeling, 2415–2446, Springer,, 2005. a

Plante, M., Tremblay, B., Losch, M., and Lemieux, J.-F.: Landfast sea ice material properties derived from ice bridge simulations using the Maxwell elasto-brittle rheology, The Cryosphere, 14, 2137–2157,, 2020. a, b, c, d, e

Rabatel, M., Labbé, S., and Weiss, J.: Dynamics of an assembly of rigid ice floes, J. Geophys. Res.-Oceans, 120, 5887–5909,, 2015. a

Rampal, P., Bouillon, S., Ólason, E., and Morlighem, M.: neXtSIM: a new Lagrangian sea ice model, The Cryosphere, 10, 1055–1073,, 2016. a, b

Ranta, J., Polojärvi, A., and Tuhkuri, J.: Limit mechanisms for ice loads on inclined structures: Buckling, Cold Reg. Sci. Technol., 147, 34–44,, 2018. a

Rhoades, C. E.: A fast algorithm for calculating particle interactions in smooth particle hydrodynamic simulations, Comput. Phys. Commun., 70, 478–482,, 1992. a

Ringeisen, D., Losch, M., Tremblay, L. B., and Hutter, N.: Simulating intersection angles between conjugate faults in sea ice with different viscous–plastic rheologies, The Cryosphere, 13, 1167–1186,, 2019. a, b

Salehizadeh, A. M. and Shafiei, A. R.: Modeling of granular column collapses with μ(I) rheology using smoothed particle hydrodynamic method, Granul. Matter, 21, 32,, 2019. a

Schreyer, H. L., Sulsky, D. L., Munday, L. B., Coon, M. D., and Kwok, R.: Elastic-decohesive constitutive model for sea ice, J. Geophys. Res., 111, C11S26,, 2006. a, b

Schulson, E. M.: Compressive shear faults within arctic sea ice: Fracture on scales large and small, J. Geophys. Res., 109, C07016,, 2004. a

Sheikh, B., Qiu, T., and Ahmadipur, A.: Comparison of SPH boundary approaches in simulating frictional soil–structure interaction, Acta Geotech., 16, 2389–2408,, 2020. a

Shen, H. T., Shen, H., and Tsai, S.-M.: Dynamic transport of river ice, J. Hydraul. Res., 28, 659–671,, 1990. a

Staroszczyk, R.: SPH Modelling of Sea-ice Pack Dynamics, Archives of Hydro-Engineering and Environmental Mechanics, 64, 115–137,, 2017. a, b, c, d, e

Staroszczyk, R.: Simulation of Sea-ice Thermodynamics by a Smoothed Particle Hydrodynamics Method, Archives of Hydro-Engineering and Environmental Mechanics, 65, 277–299,, 2018. a

Sulsky, D., Schreyer, H., Peterson, K., Kwok, R., and Coon, M.: Using the material-point method to model sea ice dynamics, J. Geophys. Res., 112, C02S90,, 2007. a

Sun, P., Colagrossi, A., Marrone, S., Antuono, M., and Zhang, A.: Multi-resolution Delta-plus-SPH with tensile instability control: Towards high Reynolds number flows, Comput. Phys. Commun., 224, 63–80,, 2018. a

Sutherland, P. and Dumont, D.: Marginal Ice Zone Thickness and Extent due to Wave Radiation Stress, J. Phys. Oceanogr., 48, 1885–1901,, 2018. a

Swegle, J., Hicks, D., and Attaway, S.: Smoothed Particle Hydrodynamics Stability Analysis, J. Comput. Phys., 116, 123–134,, 1995. a, b

Tsamados, M., Feltham, D. L., and Wilchinsky, A. V.: Impact of a new anisotropic rheology on simulations of Arctic sea ice, J. Geophys. Res.-Oceans, 118, 91–107,, 2013. a

Wang, Z., Shen, H. T., and Wu, H.: A Lagrangian sea ice model with discrete parcel method, in: Proceedings of the 14 th International Symposium on Ice, Potsdam, Germany, 1, 313–320, 1998. a, b, c, d, e, f

Weiss, J., Schulson, E. M., and Stern, H. L.: Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction, Earth Planet. Sc. Lett., 255, 1–8,, 2007. a

Wendland, H.: Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Adv. Comput. Math., 4, 389–396,, 1995. a

West, B., O'Connor, D., Parno, M., Krackow, M., and Polashenski, C.: Bonded Discrete Element Simulations of Sea Ice With Non-Local Failure: Applications to Nares Strait, J. Adv. Model. Earth Syst., 14, e2021MS002614,, 2022. a, b

Williams, J. and Tremblay, L. B.: The dependence of energy dissipation on spatial resolution in a viscous-plastic sea-ice model, Ocean Model., 130, 40–47,, 2018. a

Williams, J., Tremblay, L. B., and Lemieux, J.-F.: The effects of plastic waves on the numerical convergence of the viscous–plastic and elastic-viscous–plastic sea-ice models, J. Comput. Phys., 340, 519–533,, 2017. a, b, c, d, e, f

Xia, X. and Liang, Q.: A GPU-accelerated smoothed particle hydrodynamics (SPH) model for the shallow water equations, Environ. Model. Softw., 75, 28–43,, 2016. a

Yang, E., Bui, H. H., Sterck, H. D., Nguyen, G. D., and Bouazza, A.: A scalable parallel computing SPH framework for predictions of geophysical granular flows, Comput. Geotech., 121, 103474,, 2020. a, b

Zhang, N., Zheng, X., and Ma, Q.: Updated Smoothed Particle Hydrodynamics for Simulating Bending and Compression Failure Progress of Ice, Water, 9, 882,, 2017. a

Short summary
We developed a standard viscous–plastic sea-ice model based on the numerical framework called smoothed particle hydrodynamics. The model conforms to the theory within an error of 1 % in an idealized ridging experiment, and it is able to simulate stable ice arches. However, the method creates a dispersive plastic wave speed. The framework is efficient to simulate fractures and can take full advantage of parallelization, making it a good candidate to investigate sea-ice material properties.