the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Geometric remapping of particle distributions in the Discrete Element Model for Sea Ice (DEMSI v0.0)
Kara J. Peterson
Dan Bolintineanu
A new sea ice dynamical core, the Discrete Element Model for Sea Ice (DEMSI), is under development for use in coupled Earth system models. DEMSI is based on the discrete element method, which models collections of ice floes as interacting Lagrangian particles. In basinscale sea ice simulations the Lagrangian motion results in significant convergence and ridging, which requires periodic remapping of sea ice variables from a deformed particle configuration back to an undeformed initial distribution. At the resolution required for Earth system models we cannot resolve individual sea ice floes, so we adopt the subgridscale thickness distribution used in continuum sea ice models. This choice leads to a series of hierarchical tracers depending on ice fractional area or concentration that must be remapped consistently. The circular discrete elements employed in DEMSI help improve the computational efficiency at the cost of increased complexity in the effective element area definitions for sea ice cover that are required for the accurate enforcement of conservation. An additional challenge is the accurate remapping of element values along the ice edge, the location of which varies due to the Lagrangian motion of the particles. In this paper we describe a particletoparticle remapping approach based on wellestablished geometric remapping ideas that enforces conservation, bounds preservation, and compatibility between associated tracer quantities, while also robustly managing remapping at the ice edge. One element of the remapping algorithm is a novel optimizationbased flux correction that enforces concentration bounds in the case of nonuniform motion. We demonstrate the accuracy and utility of the algorithm in a series of numerical test cases.
Sea ice, the frozen surface of the ocean at high latitudes, forms an important component of the Earth climate system. Sea ice moderates the exchange of heat, mass, and momentum between the ocean and the atmosphere. The high albedo of sea ice has a significant effect on planetary reflectivity and can help drive the polar amplification of climate change through an albedo feedback mechanism (Ingram et al., 1989), while the rejection of salt during sea ice formation helps drive the thermohaline circulation (Killworth, 1983). The sea ice components of current global climate models use a continuum Eulerian formulation and use either structured grids (e.g., CICE: Hunke et al., 2015, or LIM3: Rousset et al., 2015) or unstructured meshes (e.g., MPASSeaice: Petersen et al., 2019). Continuum Lagrangian sea ice models have also been developed, such as neXtSIM (Rampal et al., 2016), which uses a moving triangulation as its mesh. The discrete element method (DEM) is an alternate Lagrangian method whereby the motion and contacts of finitesized particles are simulated. Several DEM sea ice models have been developed for processscale simulations, such as for sea ice ridging and deformation (Hopkins, 1994), the interaction between sea ice and solid structures (e.g., Tuhkuri and Polojärvi, 2018), the interaction between sea ice and waves (e.g., Xu et al., 2012; Herman, 2017), floe clustering (Herman, 2011), and channel flow (e.g., Gutfraind and Savage, 1998). Basinscale DEM sea ice models have been developed as well (Hopkins, 2004) and have been used to study the formation of the floe size distribution through fracture processes between elements (Hopkins and Thorndike, 2006).
Ice convergence and ridging present a unique set of issues for sea ice DEM models that are not present in traditional DEM applications. During ridging ice area is reduced as ice thickness increases. Capturing the relevant physics of this process using DEM sea ice models has proven to be a challenge. Hopkins (1994) developed a twodimensional crosssectional model of a single pressure ridge, which was used by Hopkins (1996) to derive a contact model for a twodimensional DEM sea ice model representing the plastic deformation produced during ridge formation. This model was then used to derive sea ice yield curves. Hopkins (2004) used the same contact model in sea ice simulations in which the ridging process was represented through a remapping scheme of overlapping converging neighboring elements. In this model, every 24 h the spatial distribution of elements is remapped back to the initial distribution of elements. Elements that overlap thicken after remapping, representing the ridging process. Without this remapping, element overlap would increase indefinitely, producing a highly interpenetrated element distribution for which determination of element contacts would be difficult, as would interpretation of the sea ice state at any particular point. Alternatively, instead of using overlapped elements to represent ridging, elements could be made to shrink instead. Elements with arbitrarily small radii, however, introduce computational challenges with regard to efficient contact searching, time step size, and contact model formulation (Hopkins, 2004; Shire et al., 2020). Both these considerations therefore require longduration DEM simulations of sea ice deformation to periodically perform a remapping of the model elements to an undeformed distribution. A related recent development is the nonsmooth DEM (NDEM) sea ice model developed by van den Berg et al. (2018, 2019), which for floes of any shape reproduces realistic forces for compliant floe–floe interactions.
Computational performance is another important consideration for longduration DEM sea ice simulations. One of the most computationally expensive parts of a DEM model is the detection of contact between neighboring elements. Hopkins (2004) use polygonal elements and so require a computationally expensive algorithm (Preparata and Shamos, 1985) to detect contact between elements and to calculate the contact point. Contact detection between circular elements, on the other hand, is much less computationally expensive since for circular elements only a comparison between the element separation and the element radii is needed to determine if elements are in contact. In this regard, in the work which follows, we explore how to represent sea ice in a Hopkins (2004) type of DEM model using more computationally efficient circular elements. Within this framework, we also investigate how to perform proper remapping of elements, which is necessitated by the ridging process and required for the successful development of a global DEM sea ice model performing longduration simulations.
Considerable research has gone into developing conservative and boundspreserving remapping methods for transferring scalar quantities between two grids. In the geometric approach, a reconstruction of the conserved quantity is integrated over intersections of overlapping cells between the source and target mesh. In the climate modeling community this class of algorithm has been specialized for remapping between structured and unstructured spherical grids (Jones, 1999; Ullrich et al., 2009). Similar remapping algorithms have been developed for use in arbitrary Lagrangian–Eulerian methods (Margolin and Shashkov, 2003) as well as the remapping step of semiLagrangian or incremental remapping schemes for transport (Dukowicz and Baumgardner, 2000; Lipscomb and Hunke, 2001; Lipscomb and Ringler, 2005; Lauritzen et al., 2010). For secondorder and higher remapping schemes, a form of limiting must be done to preserve physical bounds on the remapped field. Bounds preservation can be achieved through limiting the gradient of the reconstruction (Dukowicz and Baumgardner, 2000; van Leer, 1979), applying a flux correction algorithm (Liska et al., 2010), or applying an optimizationbased correction (Bochev et al., 2013, 2014).
In this paper we present a remapping method that builds on standard geometric remapping approaches while addressing the unique challenges associated with remapping for a DEM sea ice model using circular elements. These challenges include defining consistent areas for enforcing conservation, addressing monotonicity errors due to element overlap under nonuniform motion, and enabling accurate reconstructions at the ice edge. This remapping method forms part of a new sea ice dynamical core, the Discrete Element Model for Sea Ice (DEMSI), currently under development for use in coupled Earth system models. The test cases used here to illustrate the method use a greatly simplified DEM contact model, which is sufficient to demonstrate the efficacy of the remapping method. Work is ongoing to develop a contact model for DEMSI representing sea ice dynamics and including a representation of sea ice fracture and ridging. In Sect. 2 we describe the representation of sea ice in our model and introduce an effective element area needed to represent sea ice with 100 % ice concentration. In Sect. 3 we describe our remapping algorithm including a novel optimizationbased flux correction for the case of nonuniform motion and methods for ensuring accurate remap at the ice edge. In Sect. 4 we present numerical examples demonstrating that the method is robust and achieves secondorder accuracy, tracer compatibility, conservation, and bounds preservation.
While using circular elements in DEM models is computationally efficient, they present a unique challenge for the representation of sea ice. Unlike polygonal elements that can be made to completely tessellate a region (e.g., with a Voronoi tessellation), circular elements, in general, cannot completely cover a region, preventing them from directly representing sea ice cover with 100 % concentration. In order to represent 100 % ice cover with circular elements we associate an effective area, e, with each element. This area, which is potentially larger than the geometric area of the circular element, represents the area of sea ice and open water associated with that element. We determine this initial effective element area from a radical Voronoi tessellation (Imai et al., 1985) of the elements. Also called a power or Laguerre–Voronoi diagram, this is a Voronoi tessellation weighted by the element radius and with the element centers as the tessellation generator points. This tessellation results in elements that are contained completely within their polygons for element distributions that do not overlap (see Fig. 1). We take the initial effective area for an element as the area of the radical Voronoi polygon associated with that element. Since the radical Voronoi tessellation covers the whole domain, the effective element area can represent 100 % ice cover. The radical Voronoi polygon associated with an element is carried with it as the element moves.
Over time sea ice deformation and ridging reduce sea ice area (while increasing sea ice thickness and approximately conserving sea ice volume). One way to represent this in a DEM sea ice model is to allow elements to overlap as they ridge and to concomitantly reduce the effective element area as they do. The radical Voronoi polygon associated with an element, P_{i}, can then be handled in either of two ways. Firstly, the size of the polygon associated with the element, ${A}_{{P}_{i}}$, can be kept constant. In this case $e<{A}_{{P}_{i}}$ in general between remappings. Secondly, as the effective element area of an element decreases during ridging the polygon can be decreased in size so that its area remains equal to the effective element area, i.e., $e={A}_{{P}_{i}}$, between remappings. The remapping method described here will work with either methodology. We will examine ridging in DEM sea ice models in a later work but require the remapping scheme described here to periodically remap the element distribution back to the initial Voronoi tessellation to ameliorate the effect of the element overlap associated with ridging during longduration simulations.
Sea ice models used in current global climate models use numerous interdependent tracer fields that form a complex hierarchy with ice concentration, c, or the fractional area of ice in an element, as the root tracer. Sea ice models typically employ an ice thickness category distribution (e.g., Hunke et al., 2015; Rousset et al., 2015), whereby grid cell ice area is divided into a number of categories, each representing sea ice of different thicknesses. The sum of fractional areas in each thickness category is the total element ice concentration. For ice concentration less than 100 %, the remaining area is assumed to be open water. Within each category the ice is further divided into vertical layers, each of which contains tracer fields. Considering both categories and layers, MPASSeaice, for example, utilizes 23 different tracer fields for a typical physics simulation without biogeochemistry. Biogeochemistry uses many more additional tracers. Any remapping method must remap this complex tracer hierarchy in a computationally efficient manner and ensure that the sum of ice concentrations across thickness categories per element after remapping is bounded between 0 and 1.
Another desirable property of the remapping scheme is conservation of the appropriate conserved quantities. The effective area, e, provides a means to define conserved quantities that is consistent for a representation of sea ice with 100 % ice concentration. For example, given ice concentration, c, ice thickness, h, and ice enthalpy, q, in element i quantities that must be conserved during remapping are the total ice area per ice thickness category,
the total ice volume per ice thickness category,
and the total ice energy per ice thickness category and ice layer,
where each thickness category within an element is labeled by the k index and ice layers within a thickness category by the l index. In the remapping implementation we distinguish between primitive tracer variables, such as thickness and enthalpy, and conserved quantities, such as volume and energy.
Several properties are desirable in tracer remapping schemes.

Conservation. A remapping scheme should preserve conservation of the appropriate quantities. Conservation of mass and energy is very important for longterm global climate simulations.

Accuracy. Since firstorder schemes are very numerically diffusive, the method should be at least secondorder accurate in space.

Monotonicity. The method should preserve monotonicity of the tracer fields so no new extrema are generated in these fields by remapping.

Compatibility. The method should ensure that no new extrema are created when the primitive tracer field is diagnosed from the remapped conserved field; that is, both the primitive and conserved variables are consistently boundspreserving.

Computational efficiency. The method should be computationally efficient for many tracers.
The geometrically based remapping scheme proposed here satisfies all these requirements. The method proposed is based partly on the incremental remapping transport algorithm used for sea ice transport in CICE and MPASSeaice (Dukowicz and Baumgardner, 2000; Lipscomb and Hunke, 2001). The scheme for remapping a source element distribution (indexed by i) to a destination element distribution (indexed with j) is described below. This remapping method is used within a sea ice simulation to periodically remap a deformed sea ice element (source) distribution to an undeformed (destination) distribution, after which the simulation is restarted. Source element i has its effective area polygon, P_{i}, advecting with it (including both translation and rotation), while the fixed destination element j has its effective area polygon, P_{j}, associated with it. Both polygons are convex. The remapping method is based on the polygon, P_{ij}, associated with the geometric overlap between the P_{i} and P_{j} polygons, which is also convex, and consists of transferring quantities associated with P_{ij} to P_{j}. Hence, a conserved quantity Z_{i} on the source polygons is remapped to P_{j} as
where Z_{ij} is the part of Z_{i} in the intersection polygon P_{ij} and the sum is over the source element polygons, P_{i}, that overlap the destination polygon, P_{j}. Since the destination polygons tessellate the whole domain, every part of all the Z_{i} values, represented in the Z_{ij} values, is transferred to a destination polygon ensuring conservation of Z.
For a firstorder method in which a tracer, t_{i}, is assumed constant within the P_{i} polygons,
where e_{ij} is the effective area associated with the overlap polygon P_{ij} and is given by the fractional area of the overlap as
where ${A}_{{P}_{i}}$ denotes the area of the P_{i} polygon; recall that e_{i} is the effective area of element i based on tessellation of the initial element configuration. For a secondorder method, with the tracer t_{i}(x) represented as a linear function of position,
The set of source elements does not necessarily fill the entire domain, as some portions of the domain can be made up of open water. The set of destination polygons, however, forms a complete tessellation of the domain without gaps or overlap, since for conservation every part of the P_{i} polygons must overlap with a P_{j} polygon exactly once.
The major steps of the remapping method are described below.

Determine overlap polygons, P_{ij}, between the source (P_{i}) and destination (P_{j}) polygons and remap effective element area.

Compute linear reconstructions of average tracer fields on source elements based on tracer values in neighboring elements. The gradients in the reconstruction are limited to ensure monotonicity of the remapped fields in the case of uniform motion.

Integrate the conserved variables over the intersection polygons using the linear reconstructed tracer fields and aggregate in the destination polygons.

Enforce bounds preservation for remapped effective area and tracers using an optimizationbased flux correction.
These steps are described in more detail in the following sections.
3.1 Polygon intersections and remapped area
In the first step of the algorithm, the intersection polygons, P_{ij}, are calculated using the algorithm of Preparata and Shamos (1985). In order to avoid the large computational cost of calculating the intersection polygon for every source and destination polygon pair, only destination polygons that are close enough to source polygons to have potential overlap are included in the search. This is implemented with a linkcell method (Hockney et al., 1974; Plimpton, 1995).
Unlike for the discrete element application described here, when computing intersections between two welldefined grids that cover the same domain, the sum of the intersections will be equal to the domain area. Each destination grid cell will also be entirely covered by source grid cells so that the sum of the source intersections of a destination grid cell will equal the destination grid cell area. For our discrete element application, if the motion leading to the deformed element distribution is nonuniform there may be gaps and overlaps between the source polygons used to compute the intersections. Additionally, near the ice edge there will be destination polygons that are only partially covered by source polygons. To account for this, we compute a remapped effective area, e_{j}, given by
Recall that e_{i} may not be equivalent to the effective polygon area associated with element i due to area changes from ridging. The remapped effective area e_{j} will, in general, not equal the destination element effective area ${A}_{{P}_{j}}$ and in some cases may exceed the destination element effective area. The optimizationbased flux algorithm described in Sect. 3.4 is used to correct the remapped effective area to ensure ${e}_{j}\le {A}_{{P}_{j}}$ and correct tracer values with the computed area fluxes to enforce monotonicity while maintaining conservation. In the case of uniform motion there are no area changes due to ridging, and the remapped effective area simplifies such that e_{i} is equal to ${A}_{{P}_{i}}$ and e_{j} will only differ from ${A}_{{P}_{j}}$ along the ice edge.
3.2 Linear tracer reconstruction
As the first step in remapping the tracers, a meanpreserving linear reconstruction of the tracer fields, t^{p}(r), is made in each polygon of the source element distribution. For the concentration, c, thickness, h, and enthalpy, q, the reconstructions are
and
where $\mathit{r}=(x,y)$ is the position vector within the element polygon, ∇c, ∇h, and ∇q are estimates of the tracer gradients for the c, h, and q tracers in the source polygons, and α_{c}, α_{h}, and α_{q} are limiting coefficients for the c, h, and q tracers that enforce monotonicity. A linear tracer reconstruction ensures secondorder spatial accuracy of the remapping. To satisfy conservation of 𝒜, 𝒱, and 𝒬, the reconstructed tracer fields must equal the known cellaveraged tracer values (c, h, and q) when integrated over the source polygons so that
and
This requires that (Lipscomb and Hunke, 2001)
and
Tracer gradients for a source element are calculated as a multivariate linear regression of the tracer values in that element and the neighboring source elements. If two source polygons jointly overlap with a particular destination polygon they are defined as being neighboring source elements. For the n_{m} neighboring elements (including the element itself) the tracer gradients are given by
and
where
Neighboring element m has position (x_{m},y_{m}) and tracer value t_{m}.
Tracer gradients are limited with a form of van Leer limiting (van Leer, 1979) to preserve monotonicity of the tracer fields. Tracer gradients are limited so that the extrema values of the linear reconstructed tracer field, t_{p}(r), within a source polygon are within the range of tracer values for the surrounding neighbor source elements. These neighbor elements are defined in the same way as for the gradient calculation above. Since the extremal values of a linear reconstructed tracer field for a polygon are at the corners of that polygon, we use the minimum and maximum of the corner polygon tracer values of the reconstructed field, ${t}_{i}^{\text{cmin}}$ and ${t}_{i}^{\text{cmax}}$, to perform the limiting. The gradients for a source element i are limited by multiplying them by a limiting coefficient, α_{i}, which lies in the range [0,1] and whose value is given by
where
and
Here, ${t}_{i}^{min}$ and ${t}_{i}^{max}$ are the minimum and maximum, respectively, of the tracer values of the surrounding neighbor source elements to source element i. This method restricts reconstructed tracer values of a given destination element, t_{j}, to be within the range of source tracer values for source elements that overlap with the destination element as well as source element neighbors of those source elements, as in the case for the incremental remapping algorithm. In addition, in a similar way gradients of the concentration fields for all the ice thickness categories are limited so that the sum of the reconstructed concentration over thickness categories lies within [0,1].
3.3 Integration over intersection polygons
With tracer gradients calculated and limited, conserved quantities can now be calculated from integrals of the reconstructed tracer fields across the intersection polygons. For destination polygon j,
and
where the sum is over source element distribution i, and the integral is over the ij intersection polygon area. e_{ij} is nonzero only when source polygon i overlaps destination polygon j. Finally, new remapped tracer values for the destination elements are determined with
and
Only destination elements with a remapped effective area greater than zero then need to be kept for the continuing DEM simulation. To ensure consistency for partially icefilled destination elements, before the final calculation of tracer values in Eqs. (31) to (33), the effective area is set to the polygon area ${A}_{{p}_{j}}$, decreasing the concentration of sea ice calculated in Eq. (31). This is compatible with Eq. (8) since the values of effective area and concentration are not unique for representing the physically relevant quantity of sea ice area within an element, given by 𝒜_{i}=c_{i}e_{i}. Provided the area of sea ice within an element remains fixed, $\mathrm{0}\le {c}_{i}\le \mathrm{1}$, and $\mathrm{0}\le {e}_{i}\le {A}_{p}$, the concentration and effective area can vary together without affecting the physical representation.
Instead of directly calculating the integrals over the source and intersection polygon area required by the remapping method, the integrands of these integrals are expanded so that only integrals of geometric quantities of the source polygons are required. These quantities are calculated once for all tracers, improving computational efficiency for simulations with a large tracer number. For example, Eq. (16) once expanded becomes
and
and only source polygon integrals of x, y, x^{2}, y^{2}, and xy are required. The integrals are calculated by breaking polygons into triangles by connecting the polygon centroid to the polygon vertices and using the Gaussian triangular quadrature rules of Dunavant (1985).
3.4 Optimizationbased flux correction to remapping
For uniform motion of the elements there is no relative motion (either rectilinear or rotational) between elements, so their polygons do not overlap during their motion. For general motion, however, relative motion between elements results in element polygons overlapping. For the formulation of effective area from Sect. 2, this would potentially result in the remapped effective area exceeding the available geometric area of the destination polygons, resulting in remapped concentrations greater than 1. To ameliorate this issue we implement a fluxbased optimization correction to the effective area and tracer fields that is applied after the remapping discussed above. This optimization scheme determines a minimal flux, ${f}_{{j}_{\mathrm{1}}{j}_{\mathrm{2}}}$, between polygons j_{1} and j_{2} of the destination element distribution that conservatively removes excess effective area from destination polygons, so after the flux is applied to the remapped effective area and tracer fields the effective areas, ${e}_{{j}_{\mathrm{1}}}$ and ${e}_{{j}_{\mathrm{2}}}$, are less than the destination polygon areas, ${A}_{{P}_{{j}_{\mathrm{1}}}}$ and ${A}_{{P}_{{j}_{\mathrm{2}}}}$, for all destination polygons. For a given destination polygon ${P}_{{j}_{\mathrm{1}}}$ we require that the corrected effective area, ${e}_{{j}_{\mathrm{1}}}^{\prime}$, obey
The corrected effective area is given by
where ${B}_{{j}_{\mathrm{1}}{j}_{\mathrm{2}}}$ is the j_{1},j_{2} component of a sparse matrix with entries of either 0 (if polygons j_{1} and j_{2} do not share an edge), −1, or 1 (signifying whether the flux ${f}_{{j}_{\mathrm{1}}{j}_{\mathrm{2}}}$ is into or out of the polygon j_{1}). Thus, the determination of the set of fluxes, ${f}_{{j}_{\mathrm{1}}{j}_{\mathrm{2}}}$, can be cast in the following inequality constrained quadratic program:
where f is the vector of edge fluxes, ${f}_{{i}_{\mathrm{1}}{j}_{\mathrm{2}}}$, P is the identity matrix, B is the matrix with components ${B}_{{j}_{\mathrm{1}}{j}_{\mathrm{2}}}$, and u and l are vectors defined by
and
The quadratic program is solved in parallel by the PermonQP library (Hapla et al., 2016; Kružík et al., 2020; Hapla et al., 2021), which uses the PETSc library (Balay et al., 1997, 2019). The fluxes, ${f}_{{j}_{\mathrm{1}}{j}_{\mathrm{2}}}$, are then used to conservatively transfer effective area and tracer quantities between destination polygons, removing excess effective area. Since the fluxes are fluxes of effective area, a conserved quantity, Z, is corrected according to
and
for all values of f, where ${f}_{{j}_{\mathrm{1}}{j}_{\mathrm{2}}}$ flows from element j_{1} to j_{2}.
Another issue ameliorated by the flux correction is the overlapping of source elements with coastal elements. Here we represent coastlines through a series of immovable elements coincident with land. DEM models, however, model compressive interaction through an overlap of elements, so sea ice elements pushed by winds onto fixed coastal elements will overlap slightly with those elements. During the remapping described previously some sea ice will be remapped onto the coastal elements. This can be fixed by moving this remapped ice to the nearest destination element that is not a coastal element before the flux correction is performed. While moving this ice has the potential to increase the effective area of these elements to be larger than their geometric area, the flux correction algorithm ameliorates this effect.
3.5 Effect of open water
One difference between particle and Eulerian methods that causes potential issues for remapping is the treatment of open water. In Eulerian methods, empty Eulerian cells, representing open water, can be included in the calculation of tracer gradients and monotonicity limits. For DEM sea ice models, elements only exist where sea ice exists, with elements that lose all their ice being removed from the simulation. One advantage of this is that fewer elements are needed than Eulerian cells (which must cover the whole domain), increasing computational performance for the DEM model. Consequently, there are no empty elements representing open water to be included in tracer gradient calculations or in the neighboring element tracer extrema used to limit those gradients. A possible effect of this is that tracer values for ice concentration and thickness, t_{i}, for elements at the ice pack edge are likely to be the minimum of the neighbor tracer values so that ${t}_{i}={t}_{ij}^{min}$. From Eqs. (25) to (27) this would mean the gradients would be limited to zero in these ice edge elements, and the remapping here would revert to a highly diffusive firstorder method. To investigate the effect of this issue we implement a method to account for open water in the tracer gradient calculation and limiting. Firstly, for the neighboring source elements, N_{i}, of a given source element (not including the element itself) we calculate the total effective area, E_{i}, and the centroid of that effective area, ${\mathit{r}}_{{E}_{i}}$. These are given by
and
If the E_{i} is below some limit we assume that there is some open water surrounding the source element and include an additional openwater neighbor pseudoelement into the gradient and gradientlimiting calculations. This pseudoelement has zero concentration and a position on the opposite side of the element to the centroid of the surrounding effective area ($\mathrm{2}{\mathit{r}}_{i}{\mathit{r}}_{{E}_{i}}$). The addition of a neighbor concentration value of zero would prevent the gradient at the ice edge from being limited to zero. We include the zero concentration pseudoelement if the sea ice surrounding the element of interest is below a fraction (taken as 0.75) of an estimate of the expected value of E_{i} for a fully covered pack, ${E}_{i}^{\prime}$. In two dimensions we estimate ${E}_{i}^{\prime}$ based on a hexagonal close packing of elements so that ${E}_{i}^{\prime}\simeq \mathrm{6}{\overline{A}}_{{P}_{j}}$, where ${\overline{A}}_{{P}_{j}}$ is the mean area of destination polygons. In one dimension each element is only surrounded by two elements, so for onedimensional tests we use ${E}_{i}^{\prime}\simeq \mathrm{2}{\overline{A}}_{{P}_{j}}$.
We now test the remapping algorithm described in the previous section with a series of idealized test cases. Motion of the elements between remapping is determined by the LAMMPS (Plimpton, 1995) molecular dynamics code, which has the capability to run DEM simulations, wherein particle interaction forces are computed based on contacts and the equations of motions are integrated in time. LAMMPS can either explicitly enforce the motion of elements (such as in uniform rectilinear motion) or solve the equations of motion of the elements with a modified velocity Verlet integrator (Swope et al., 1982). For test cases involving solving the equations of motion, we solve a simplified sea ice momentum equation, which in a DEM context reduces to a set of coupled firstorder equations for all elements:
where m_{i} is the mass of each element, u_{i} is the element velocity, F_{ij} is the contact force on element i from neighboring element j, and τ_{a} and τ_{w} are the atmospheric and ocean stresses, respectively, evaluated at the location of the element x_{i}. To determine element rotation we also solve an angular momentum equation given by
where I_{i} is the moment of inertia of element i, ω_{i} is the element angular velocity, (r_{ij}×F_{ij}) is the torque on element i produced by the contact force with element j, and r_{ij} is the contact position between elements i and j relative to the center of mass of element i. The surface stresses have a quadratic form given by
where e is the effective element area, c_{a,w} is a drag coefficient, ρ_{a,w} is the fluid density, U_{a,w} is the fluid velocity, subscript a refers to air, and subscript w refers to water. This momentum equation is discretized in time as
where superscripts n and n−1 signify the current and previous time steps, respectively, and the ocean velocity is treated partly implicitly for stability. For the test cases considered here a simple Hookean elastic repulsion contact model is employed with velocity damping in both the normal and tangential contact directions. While the Hookean contact model has several limitations, such as sensitivity to the time step and to the time stepping method used, it is sufficient to demonstrate the efficacy of the remapping method. The air drag coefficient, c_{a}, and air density, ρ_{a}, are given by 0.0012 and 1.3 kg m^{−3}, respectively (Hibler, 1979), while the ocean drag coefficient, c_{w}, and water density, ρ_{w}, are given by 0.00536 and 1026 kg m^{−3} (Kim et al., 2006), respectively.
4.1 Onedimensional uniform motion
The first test we perform is a simple onedimensional test case with uniform motion and a “top hat” initial ice concentration distribution. The destination element distribution consists of uniform elements with radius, r, of 500 m arranged in the x direction across a domain 1000 km long. The domain accommodates a single element in the y direction and is 1 km wide. The initial source distribution is a subset of the destination distribution consisting of the 100 elements between 100 and 200 km in the x direction. Initial ice concentration, c_{0}(x), is given by
where x_{1} and x_{2} are 100 and 200 km, respectively.
These elements have an initial effective area of a square of side 2r, a sea ice concentration of 1, and a sea ice thickness of 1 m. The elements undergo uniform motion in the positive x direction for a total of 100 km of translation, whereupon the final source distribution is examined. Figure 2 shows this final distribution for three situations. The first situation (“no remap” in Fig. 2) has no remapping performed during the uniform motion so that the final distribution of ice concentration and thickness is the same as the initial one, other than being translated 100 km in the positive x direction. For the second situation (“low order” in Fig. 2) the source distribution elements move half their diameter (500 m) before remapping and undergo a total of 200 remappings, resulting in the 100 km total translation. Tracer gradients are set to zero so that the remapping algorithm becomes a loworder method. The final situation (“higher order” in Fig. 2) is the same as the loworder situation except tracer gradients are not set to zero, making the remapping a more accurate higherorder method. The no remap situation is useful for demonstrating the final distribution that a perfect remapping scheme, with no numerical diffusion, would produce. Figure 2 shows that the higherorder method produces significantly less numerical diffusion of the concentration profile than the loworder method.
We now use this test case to explore the effect of taking account of open water in the tracer gradient and limiting, as described in Sect. 3.5. We rerun the higherorder simulation with and without the effects of open water taken into account (see Fig. 3). As can be seen in Fig. 3 the effect of the absence of openwater elements appears negligible with only a very minor amount of extra numerical diffusion present at the very edge of the ice due to this absence (see inset in Fig. 3b). Similar results were found with the twodimensional simulations described next. Based on these results we neglect the effect of the absence of openwater elements for the remainder of this work.
The next onedimensional test case is based on that found in Lipscomb and Hunke (2001) and is used to test how well the remapping algorithm preserves the compatibility of the particle tracers. The test case is the same as the previous onedimensional test case, except for the initial distribution of ice concentration, thickness, and volume. These initial distributions are more complex than the previous test case (see Fig. 4a and b) and provide a harder test of compatibility. Initial ice concentration, c_{0}(x), and thickness, h_{0}(x), are given by
and
where x_{1}, x_{2}, x_{3}, x_{4}, and x_{5} are given by 100, 112.5, 150, 187.5, and 200 km, respectively. As before we perform three separate simulations: one with no remapping (Fig. 4a and b), one with loworder remapping (Fig. 4c and d), and one with higherorder remapping (Fig. 4e and f). As before both the low and higherorder remapping schemes are capable of remapping the initial particle distribution, with the higherorder scheme producing much less numerical diffusion. Figure 4d and f also clearly show that the thickness tracer preserves monotonicity as well as the concentration and volume, clearly demonstrating compatibility of the method.
4.2 Twodimensional uniform motion
The next test case we examine involves a twodimensional domain of size 300 km in the x direction and 100 km size in the y direction. The destination and initial source element distributions consist of a random arrangement of circular elements with a distribution of element radius, with mean radius ∼ 513 m. To create the tightly packed initial element distribution we use the method of Liang et al. (2015). This method uses a Lloyd (1982) type of iteration of the radical Voronoi tessellation of an initial random set of generator points each with an assigned radius drawn from a distribution. During each iteration the generator locations are moved to the centers of the largest inscribed circles within the radical Voronoi polygons (see Fig. 1 for an example element distribution generated with this method). We explore two different initial source element distributions discussed in Lipscomb and Ringler (2005). The first is a cosine bell distribution with an initial ice concentration, c_{0}(x,y), given by
where $r=\sqrt{(x{x}_{\mathrm{0}}{)}^{\mathrm{2}}+(y{x}_{\mathrm{0}}{)}^{\mathrm{2}}}$, (x_{0},y_{0}) is the center of the cosine bell at (25 km, 50 km), and r_{0} is the radius of the cosine bell with size 15 km. The second is a slotted cylinder with an initial ice concentration, c_{0}(x,y), given by
where x_{0}, y_{0}, and r_{0} are the same as for the cosine bell distribution. Elements on the edge of the slotted cylinder are assigned an initial concentration according to their areal fraction within the slotted cylinder. The element distributions are moved uniformly in the x direction and examined after 100 km of translation. As before we perform three simulations (one with no remapping, one with loworder remapping, and one with higherorder remapping), with remapping performed every 500 m of translation for a total of 200 remappings. Spatial maps of the element concentration at the end of the simulation are shown in Fig. 5. The higherorder remapping scheme is able to preserve the approximate shape of the initial concentration distributions, while high numerical diffusion with the loworder remapping scheme leaves the final element distribution more spread out with any structure of the slotted cylinder lost. The difference in numerical diffusion of the two remapping methods is evident in cross sections of final ice concentration through the middle of the distributions, as shown in Fig. 6. As can be seen almost no trace of the slot in the slotted cylinder is present for the loworder scheme, while the slot is preserved well for the higherorder method.
Next we study the error scaling properties of the remapping method. To do this we create various initial element distributions with different mean particle radii but the same relative distribution of radii. We rerun the simulations for these element distributions and record the L_{2} error norm for each simulation. We calculate the L_{2} error norm as
where the sum is over the final element distribution, e_{i} and c_{i} are the effective area and concentration of the final distribution, and ${c}_{i}^{\mathrm{0}}$ is the ice concentration for no remapping. In Fig. 7, the values of the L_{2} error norm are plotted against the mean radius of the element distribution for both the cosine bell and slotted cylinder initial distributions as well as the low and higherorder remapping methods. As can be seen from the figure the higherorder method has significantly lower errors than the loworder method for both initial distributions. The gradient of the slope for the smooth cosine bell distribution is also much higher than the slotted cylinder with its sharp discontinuities. It is also apparent that for the smooth cosine bell distribution the remapping method is approximately secondorder accurate for the higherorder method.
4.3 Flux correction tests
Having demonstrated in Sects. 4.1 and 4.2 that the higherorder remapping method for uniform motion is capable of remapping the particle distribution with low numerical diffusion while preserving compatibility, we next turn to examining the effect of nonuniform motion and testing the ability of the flux correction algorithm described in Sect. 3.4 to correct the issues created by this motion. The first test we examine is a variation of the twodimensional test used in the previous section. Instead of making the elements undergo uniform motion, we fix their positions and randomize the orientations of their associated polygons before each remapping. This introduces the overlaps and gaps between source polygons that might be expected from a general motion of the elements. As before, a total of 200 remappings are performed before the concentration distributions are examined. Figure 8b and e show the effect of these randomized orientations on the remapped ice concentration for the higherorder remapping scheme. As can be seen, while the general shape of the initial distributions is preserved, the randomization results in ice concentrations significantly larger than 1 for some elements (where elements overlapped) and ice concentration less than that expected for others (from gaps in the polygon distribution), with these effects taking on a random nature. Clearly, without correction this effect would have a seriously detrimental effect for any sea ice simulation. Figure 8c and f show the same final ice concentration distribution but generated with the flux correction scheme described in Sect. 3.4 performed after remapping. Now, ice concentrations are bounded correctly by [0,1], and the final distribution is smooth with results similar to those obtained in Sect. 4.2. Exactly the same result is not expected as the distribution of particles for remapping is different and the remapping method tested in Sect. 4.2 does not require the flux correction scheme. However, the similarity in the results indicates that the flux correction scheme works effectively and does not introduce significantly more numerical diffusion. This is also evident in the cross section of ice concentration shown in Fig. 9. Repeating the L_{2} error scaling analysis in the previous section (Fig. 10) shows larger L_{2} error norms for the uncorrected simulations, while the fluxcorrected simulations show significantly lower errors and also show secondorder convergence with mean particle size for the smooth cosine bell test case.
To demonstrate the effect of the flux correction in more realistic situations we present two more test cases. In the first, adapted from Lipscomb et al. (2007), we have a square domain of size 1000 km by 1000 km with a random distribution of elements covering the domain generated by the algorithm of Liang et al. (2015). In the middle of the domain is an “Lshaped” island consisting of fixed coastal elements. Initial elements are coastal or sea ice elements according to the following. Here, (x,y) is the position of the element center, x_{0} and y_{0} equal 400 km, x_{1} and y_{1} equal 600 km, and x_{2} and y_{2} equal 550 km. Elements have a mean element radius of ∼ 10 km, with noncoastal elements having a concentration of 1 everywhere. A constant wind forcing is applied to the elements with the x and y components of wind velocity both equal to 10 m s^{−1}, while the ocean is stationary. During the simulation the elements flow around the island, leaving an openwater wake behind it. Since the contact model used here consists of elastic repulsion without a representation of ridging, ice does not build up on the windward side of the island. The simulation is run for 2 d after which the element distribution is remapped to the destination element distribution with the higherorder scheme and either with or without the flux correction algorithm. Figure 11a and b show the effect on ice concentration of remapping without and with the flux correction, respectively. Figure 11c and d show the same but with a value of 1 subtracted from the ice concentration to highlight the excess ice concentration generated by the remapping when the flux correction is not used. These figures clearly show elements with concentration greater than 1 after remapping without flux correction and also show how the flux correction removes this excess concentration and results in a concentration field closer to the concentration field before remapping with a concentration of 1 most places within the pack. Elements on the edge of the ice pack and the edge of the island wake have concentrations less than 1, since during the remapping these elements were not fully covered by source elements.
The final test case we examine is that from Flato (1993). Here the domain consists of a 500 km by 500 km closed square with the upper half (y > 250 km) initially filled with elements with a concentration of 1 and a mean radius of ∼ 2.5 km. The elements have the same air and ocean stress formulation and the same elastic repulsion contact model as the previous island test case. The ocean is again stationary, but the wind field, U_{a}, takes the form of a vortex with
where r is the position vector relative to the center of the vortex at (250, 200 km), ω is 0.5 × 10^{−3} s^{−1}, λ is 8 × 10^{5} m^{2} s^{−1}, and k is the unit vertical vector. The simulation is run for 2 d and then the element distribution is remapped. As for the island test case, Fig. 12 shows that the flux correction scheme is capable of removing excess concentration caused by more realistic nonuniform motion of the elements for this test case.
Modeling sea ice dynamics with the discrete element method has the potential to enable the capture of the anisotropic deformation and fracture seen in observational data that are difficult to reproduce in continuum models. The DEM models, however, present a challenge due to convergence and overlap of the Lagrangian elements in simulations. To manage this element clustering, a deformed element distribution is periodically remapped to an undeformed distribution. This paper presents a geometrically based remapping algorithm designed to address the unique challenges of accurately remapping sea ice tracer fields for circular discrete elements. In particular, the method includes a representation of effective element area defined with a radical Voronoi tessellation for the undeformed element configuration that enables an accurate definition of area conservation in the case of 100 % sea ice concentration. Our remapping approach builds on conservative, boundspreserving, and compatible algorithms developed for incremental remapping and applied to sea ice (Lipscomb and Hunke, 2001; Lipscomb and Ringler, 2005). In the case of uniform motion, the effective element area enables a complete coverage of the ice domain without overlap, and the method inherits the properties of the incremental remapping algorithms. In the case of nonuniform motion in which gaps and overlaps occur between the deformed effective element areas, monotonicity can be lost. To address this, we developed a novel fluxbased correction method, which maintains conservation while correcting for bounds violations due to overlapping elements. Numerical examples are provided that demonstrate the secondorder accuracy of the method, bounds preservation for both uniform and nonuniform motion, and compatibility between primitive and conserved tracer quantities.
The code and test cases used in this paper can be found at https://doi.org/10.5281/zenodo.5800083 (Turner et al., 2021).
Input test case data used in this paper can be found at https://doi.org/10.5281/zenodo.6226411 (Turner et al., 2022).
AKT developed the methodology of the remapping scheme and the remapping software. AKT also implemented the test cases and performed model validation and analysis. KJP contributed to the remapping methodology and software development, as well as model analysis. DB integrated the LAMMPS model into the remainder of the code. AKT and KJP wrote this paper with DB contributing to the editing of the paper. All authors contributed to funding acquisition.
The contact author has declared that neither they nor their coauthors have any competing interests.
This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the U.S. Government.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors wish to thank the PermonQP team for help in using the PermonQP library and Andrew Roberts for contributing to the writing of the SciDAC DEMSI proposals.
This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research and Office of Biological and Environmental Research, Scientific Discovery through Advanced Computing (SciDAC) program. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration (contract no. DENA0003525).
This paper was edited by Alexander Robel and reviewed by two anonymous referees.
Balay, S., Gropp, W. D., McInnes, L. C., and Smith, B. F.: Efficient Management of Parallelism in Object Oriented Numerical Software Libraries, in: Modern Software Tools in Scientific Computing, edited by: Arge, E., Bruaset, A. M., and Langtangen, H. P., Birkhäuser Press, Boston, MA, pp. 163–202, 1997. a
Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W. D., Karpeyev, D., Kaushik, D., Knepley, M. G., May, D. A., McInnes, L. C., Mills, R. T., Munson, T., Rupp, K., Sanan, P., Smith, B. F., Zampini, S., Zhang, H., and Zhang, H.: PETSc Users Manual, Tech. Rep. ANL95/11 – Revision 3.11, Argonne National Laboratory, Lemont, IL, 2019. a
Bochev, P., Ridzal, D., and Shashkov, M.: Fast optimizationbased conservative remap of scalar fields through aggregate mass transfer, J. Comput. Phys., 246, 37–57, https://doi.org/10.1016/j.jcp.2013.03.040, 2013. a
Bochev, P., Ridzal, D., and Peterson, K.: Optimizationbased remap and transport: A divide and conquer strategy for featurepreserving discretizations, J. Comp. Physics, 257, Part B, 1113–1139, https://doi.org/10.1016/j.jcp.2013.03.057, 2014. a
Dukowicz, J. K. and Baumgardner, J. R.: Incremental Remapping as a Transport/Advection Algorithm, J. Comput. Phys., 160, 318–335, https://doi.org/10.1006/jcph.2000.6465, 2000. a, b, c
Dunavant, D. A.: High degree efficient symmetrical Gaussian quadrature rules for the triangle, Int. J. Numer. Meth. Eng., 21, 1129–1148, https://doi.org/10.1002/nme.1620210612, 1985. a
Flato, G. M.: A particleincell seaice model, Atmos.Ocean, 31, 339–358, https://doi.org/10.1080/07055900.1993.9649475, 1993. a
Gutfraind, R. and Savage, S. B.: Flow of fractured ice through wedgeshaped channels: smoothed particle hydrodynamics and discreteelement simulations, Mech. Mater., 29, 1–17, https://doi.org/10.1016/S01676636(97)000720, 1998. a
Hapla, V., Horák, D., Pospíšil, L., Čermák, M., Vašatová, A., and Sojka, R.: Solving Contact Mechanics Problems with PERMON, in: High Performance Computing in Science and Engineering, vol. 9611 of Lecture Notes in Computer Science, edited by: Kozubek, T., Blaheta, R., Šístek, J., Rozložník, M., and Čermák, M., Springer International Publishing Switzerland, Cham, Switzerland, https://doi.org/10.1007/9783319403618_7, pp. 101–115, 2016. a
Hapla, V., Horák, D., Kružík, J., Pecha, M., Pospíšil, L., Sojka, R., Vašatová, A., Čermák, M., Dostál, Z., and Markopoulos, A.: PERMON Web page, http://permon.vsb.cz/ (last access: 3 February 2021), 2021. a
Herman, A.: Moleculardynamics simulation of clustering processes in seaice floes, Phys. Rev. E, 84, 056104, https://doi.org/10.1103/PhysRevE.84.056104, 2011. a
Herman, A.: Waveinduced stress and breaking of sea ice in a coupled hydrodynamic discreteelement wave–ice model, The Cryosphere, 11, 2711–2725, https://doi.org/10.5194/tc1127112017, 2017. a
Hibler, W. D.: A Dynamic Thermodynamic Sea Ice Model, J. Phys. Oceanogr., 9, 815–846, https://doi.org/10.1175/15200485(1979)009<0815:ADTSIM>2.0.CO;2, 1979. a
Hockney, R., Goel, S., and Eastwood, J.: Quiet highresolution computer models of a plasma, J. Comput. Phys., 14, 148–158, https://doi.org/10.1016/00219991(74)900102, 1974. a
Hopkins, M.: A discrete element Lagrangian sea ice model, Eng. Computation., 21, 409–421, https://doi.org/10.1108/02644400410519857, 2004. a, b, c, d, e
Hopkins, M. A.: On the ridging of intact lead ice, J. Geophys. Res.Oceans, 99, 16351–16360, https://doi.org/10.1029/94JC00996, 1994. a, b
Hopkins, M. A.: On the mesoscale interaction of lead ice and floes, J. Geophys. Res.Oceans, 101, 18315–18326, https://doi.org/10.1029/96JC01689, 1996. a
Hopkins, M. A. and Thorndike, A. S.: Floe formation in Arctic sea ice, J. Geophys. Res.Oceans, 111, C11S23, https://doi.org/10.1029/2005JC003352, 2006. a
Hunke, E. C., Lipscomb, W. H., Turner, A. K., Jeffery, N., and Elliott, S.: CICE: the Los Alamos sea ice model documentation and software user's manual version 5.1, Tech. rep., Los Alamos National Laboratory, Los Alamos, NM, 2015. a, b
Imai, H., Iri, M., and Murota, K.: Voronoi Diagram in the Laguerre Geometry and Its Applications, SIAM J. Comput., 14, 93–105, https://doi.org/10.1137/0214006, 1985. a
Ingram, W. J., Wilson, C. A., and Mitchell, J. F. B.: Modeling climate change: An assessment of sea ice and surface albedo feedbacks, J. Geophys. Res.Atmos., 94, 8609–8622, https://doi.org/10.1029/JD094iD06p08609, 1989. a
Jones, P. W.: First and SecondOrder Conservative Remapping Schemes for Grids in Spherical Coordinates, Mon. Weather Rev., 127, 2204–2210, https://doi.org/10.1175/15200493(1999)127<2204:FASOCR>2.0.CO;2, 1999. a
Killworth, P. D.: Deep convection in the World Ocean, Rev. Geophys., 21, 1–26, https://doi.org/10.1029/RG021i001p00001, 1983. a
Kim, J. G., Hunke, E. C., and Lipscomb, W. H.: Sensitivity analysis and parameter tuning scheme for global seaice modeling, Ocean Model., 14, 61–80, https://doi.org/10.1016/j.ocemod.2006.03.003, 2006. a
Kružík, J., Horák, D., Čermák, M., Pospíšil, L., and Pecha, M.: Active set expansion strategies in MPRGP algorithm, Adv. Eng. Softw., 149, 102895, https://doi.org/10.1016/j.advengsoft.2020.102895, 2020. a
Lauritzen, P. H., Nair, R. D., and Ullrich, P. A.: A conservative semiLagrangian multitracer transport scheme (CSLAM) on the cubedsphere grid, J. Comput. Phys., 229, 1401–1424, https://doi.org/10.1016/j.jcp.2009.10.036, 2010. a
Liang, G., Lu, L., Chen, Z., and Yang, C.: Poisson disk sampling through disk packing, Comp. Visual Media, 1, 17–26, https://doi.org/10.1007/s4109501500037, 2015. a, b, c
Lipscomb, W. H. and Hunke, E. C.: Modeling Sea Ice Transport Using Incremental Remapping, Mon. Weather Rev., 132, 1341–1354, https://doi.org/10.1175/15200493(2004)132<1341:MSITUI>2.0.CO;2, 2001. a, b, c, d, e
Lipscomb, W. H. and Ringler, T. D.: An Incremental Remapping Transport Scheme on a Spherical Geodesic Grid, Mon. Weather Rev., 133, 2335–2350, https://doi.org/10.1175/MWR2983.1, 2005. a, b, c
Lipscomb, W. H., Hunke, E. C., Maslowski, W., and Jakacki, J.: Ridging, strength, and stability in highresolution sea ice models, J. Geophys. Res.Oceans, 112, C03S91, https://doi.org/10.1029/2005JC003355, 2007. a
Liska, R., Shashkov, M., Váchal, P., and Wendroff, B.: Optimizationbased synchronized fluxcorrected conservative interpolation (remapping) of mass and momentum for arbitrary Lagrangian–Eulerian methods, J. Comput. Phys., 229, 1467–1497, https://doi.org/10.1016/j.jcp.2009.10.039, 2010. a
Lloyd, S.: Least squares quantization in PCM, IEEE T. Inform. Theory, 28, 129–137, https://doi.org/10.1109/TIT.1982.1056489, 1982. a
Margolin, L. G. and Shashkov, M.: Secondorder signpreserving conservative interpolation (remapping) on general grids, J. Comput. Phys., 184, 266–298, https://doi.org/10.1016/S00219991(02)000335, 2003. a
Petersen, M. R., AsayDavis, X. S., Berres, A. S., Chen, Q., Feige, N., Hoffman, M. J., Jacobsen, D. W., Jones, P. W., Maltrud, M. E., Price, S. F., Ringler, T. D., Streletz, G. J., Turner, A. K., Van Roekel, L. P., Veneziani, M., Wolfe, J. D., Wolfram, P. J., and Woodring, J. L.: An Evaluation of the Ocean and Sea Ice Climate of E3SM Using MPAS and Interannual COREII Forcing, J. Adv. Model. Earth Sy., 11, 1438–1458, https://doi.org/10.1029/2018MS001373, 2019. a
Plimpton, S.: Fast Parallel Algorithms for ShortRange Molecular Dynamics, J. Comput. Phys., 117, 1–19, https://doi.org/10.1006/jcph.1995.1039, 1995. a, b
Preparata, F. P. and Shamos, M. I.: Computational Geometry: An introduction, SpringerVerlag, Berlin, https://doi.org/10.1007/9781461210986, 1985. a, b
Rampal, P., Bouillon, S., Ólason, E., and Morlighem, M.: neXtSIM: a new Lagrangian sea ice model, The Cryosphere, 10, 1055–1073, https://doi.org/10.5194/tc1010552016, 2016. a
Rousset, C., Vancoppenolle, M., Madec, G., Fichefet, T., Flavoni, S., Barthélemy, A., Benshila, R., Chanut, J., Levy, C., Masson, S., and Vivier, F.: The LouvainLaNeuve sea ice model LIM3.6: global and regional capabilities, Geosci. Model Dev., 8, 2991–3005, https://doi.org/10.5194/gmd829912015, 2015. a, b
Shire, T., Hanley, K. J., and Stratford, K.: DEM simulations of polydisperse media: efficient contact detection applied to investigate the quasistatic limit, Computational Particle Mechanics, 8, 653–663, https://doi.org/10.1007/s40571020003612, 2020. a
Swope, W. C., Andersen, H. C., Berens, P. H., and Wilson, K. R.: A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters, J. Chem. Phys., 76, 637–649, https://doi.org/10.1063/1.442716, 1982. a
Tuhkuri, J. and Polojärvi, A.: A review of discrete element simulation of ice–structure interaction, Philos. T. R. Soc. A, 376, 20170335, https://doi.org/10.1098/rsta.2017.0335, 2018. a
Turner, A. K., Peterson, K. J., Bolintineanu, D., Kuberry, P. A. , Clemmer, J. T., and Nikolov, S.: DEMSI release v0.0.1 (0.0.1), Zenodo [code], https://doi.org/10.5281/zenodo.5800083, 2021. a
Turner, A. K., Peterson, K. J., Bolintineanu, D., Kuberry, P. A., Clemmer, J. T., and Nikolov, S.: DEMSI v0.0.1 input data (0.0.1), Zenodo [data set], https://doi.org/10.5281/zenodo.6226411, 2022. a
Ullrich, P. A., Lauritzen, P. H., and Jablonowski, C.: Geometrically Exact Conservative Remapping (GECoRe): Regular latitude–longitude and cubed–sphere grids, Mon. Weather Rev., 137, 1721–1741, https://doi.org/10.1175/2008MWR2817.1, 2009. a
van den Berg, M., Lubbad, R., and Loset, S.: An implicit timestepping scheme and an improved contact model for ice–structure interaction simulations, Cold Reg. Sci. Technol., 155, 193–213, https://doi.org/10.1016/j.coldregions.2018.07.001, 2018. a
van den Berg, M., Lubbad, R., and Løset, S.: The effect of ice floe shape on the load experienced by verticalsided structures interacting with a broken ice field, Mar. Struct., 65, 229–248, https://doi.org/10.1016/j.marstruc.2019.01.011, 2019. a
van Leer, B.: Towards the ultimate conservative difference scheme. V. A secondorder sequel to Godunov's method, J. Comput. Phys., 32, 101–136, https://doi.org/10.1016/00219991(79)901451, 1979. a, b
Xu, Z., Tartakovsky, A. M., and Pan, W.: Discreteelement model for the interaction between ocean waves and sea ice, Phys. Rev. E, 85, 016703, https://doi.org/10.1103/PhysRevE.85.016703, 2012. a