Geometric remapping of particle distributions in the Discrete Element Model for Sea Ice (DEMSI v0.0)

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 basin-scale 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 5 the sub-gridscale 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. 10 In this paper we describe a particle-to-particle remapping approach based on well-established 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 optimization-based flux correction that enforces concentration bounds in the case of non-uniform motion. We demonstrate the accuracy and utility of the algorithm in a series of numerical test cases. 15

the sub-gridscale 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. 10 In this paper we describe a particle-to-particle remapping approach based on well-established 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 optimization-based flux correction that enforces concentration bounds in the case of non-uniform motion. We demonstrate the accuracy and utility of the algorithm in a series of numerical test cases. 15

Introduction
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 20 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.
Over time sea-ice deformation and ridging reduces 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 95 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 Pi , can be kept constant.
In this case e < A Pi 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 Pi between remappings. The remapping method described here will work with either methodology. We will examine ridging in DEM sea 100 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 long duration 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 fractional area of ice in an element, as the root tracer. Sea ice models typically employ an ice 105 thickness category distribution (e.g. Hunke et al., 2015;Rousset et al., 2015), where 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 contain tracer fields. Considering both categories and layers, MPAS-Seaice, for example, utilizes 23 different tracer fields for a typical physics simulation without biogeochemistry.

110
Biogeochemistry uses many more additional tracers. Any remapping method must remap this complex tracer hierarchy in a computationally efficient manner, as well as ensuring 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 con-115 centration. 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 120 and the total ice energy per ice thickness category and ice layer where each thickness category within an element is labelled 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.

Geometric remapping implementation
Several properties are desirable in tracer remapping schemes: -Conservation: A remapping scheme should preserve conservation of the appropriate quantities. Conservation of mass and energy are very important for long-term global climate simulations.
-Accuracy: Since first order schemes are very numerically diffusive, the method should be at least second order accurate 130 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 bounds-preserving.

135
-Computational efficiency: The method should be computationally efficient for many tracers.
The geometric 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 MPAS-Seaice (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 peri-140 odically 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 145 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, are transferred to a destination polygon ensuring conservation of Z.

150
For a first order method, where 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 Pi denotes the area of the P i polygon, and recall that e i is the effective area of element i based on tessellation of the 155 initial element configuration. For a second order 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, form 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.

160
The major steps of the remapping method are described below: 1. Determine overlap polygons, P ij , between the source (P i ) and destination (P j ) polygons and remap effective element area.
2. 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 165 uniform motion.
3. Integrate the conserved variables over the intersection polygons using the linear reconstructed tracer fields and aggregate in the destination polygons.
4. Enforce bounds preservation for remapped effective area and tracers using an optimization-based flux correction.
These steps are described in more detail in the following sections.

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 link-cell method (Hockney et al., 1974;Plimpton, 1995).

175
Unlike for the discrete element application described here, when computing intersections between two well-defined 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 non-uniform there may be gaps and overlaps between the source polygons used to compute the intersections. Additionally, 180 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 Pj and in some cases may 185 exceed the destination element effective area. The optimization-based flux algorithm described in Section 3.4 is used to correct the remap effective area to ensure e j ≤ A Pj 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 Pi and e j will only differ from A Pj along the ice edge.

190
As the first step in remapping the tracers, a mean-preserving linear reconstruction of the tracer fields, t p (r) :::: 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 195 and q p p (r) = q + α q ∇q · (r −r) where 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 second-order spatial accuracy of the remapping. To satisfy conservation 200 of A, V, and Q, the reconstructed tracer fields must equal the known cell-averaged tracer values (c, h, and q) when integrated over the source polygons so that 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 215 :: as ::::: being :::::::::: neighboring ::::: source :::::::: elements. : For the n m neighbouring elements (including the element itself) the tracer gradients are given by and , 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 and where t min i and t max i 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 neighbours of those source elements, as in the case for the incremental remapping algorithm. In addition, in a similar way gradients 250 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].

275
Instead of directly calculating the integrals over 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 large tracer number.
For example equation 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).

Optimization-based flux correction to remapping 290
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 Section 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 flux based optimization correction to the effective area and tracer fields that 295 is applied after the remapping discussed above. This optimization scheme determines a minimal flux, f j1j2 , 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 j1 and e j2 , are less than the destination polygon areas, A Pj 1 and A Pj 2 for all destination polygons. For a given destination polygon P j1 we require that the corrected effective area, e j1 , obey The corrected effective area is given by where B j1j2 is 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 j1j2 is into or out of the polygon j 1 ). Thus the determination of the set of fluxes, f j1 , can be cast in the following 305 inequality constrained quadratic program: where f is the vector of edge fluxes, f i1j2 , P is the identity matrix, B is the matrix B j1j2 , 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(Balay et al., , 2019. The fluxes, f j1j2 , 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 j1j2 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 320 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 325 larger than their geometric area, the flux correction algorithm ameliorates this effect.

Effect of open water
One difference between particle and Eulerian methods that causes potential issues for remapping is the treatment of open water.
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 pseudo-element 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 . In two dimensions we estimate E i based on a hexagonal close packing of elements, so that E i 6Ā Pj , whereĀ Pj is the mean area of destination polygons. In one 350 dimension each element is only surrounded by two elements so for one dimensional tests we use E i 2Ā Pj .

One dimensional uniform motion
The first test we perform is a simple one dimensional 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 385 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 km and 200 km in the x direction. Initial ice concentration, c 0 (x), is given by where x 1 and x 2 are 100 km and 200 km, respectively.
Tracer gradients are set to zero so that the remapping algorithm becomes a low order method. The final situation ("Higher order" in Figure 2) is the same as the "Low order" situation except tracer gradients are not set to zero, making the remapping 400 a more accurate higher order 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 "Higher order" method produces significantly less numerical diffusion of the concentration profile than the "Low order" 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 section 3.5. We rerun the "Higher order" simulation with and without the effects of open water taken into account 405 (see Figure 3). As can be seen in Figure 3 the effect of the absence of open water 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 Figure 3b).
Similar results were found with the two dimensional simulations described next. Based on these results we neglect the effect of the absence of open water elements for the remainder of this work.
The next one dimensional test case is based on that found in Lipscomb and Hunke (2001), and is used to test how well the 410 remapping algorithm preserves the compatibility of the particle tracers. The test case is the same as the previous one dimensional 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 Figure 4a,b) and provide a harder test of compatibility. Initial ice concentration, c 0 (x), and thickness, h 0 (x), are given by where x 1 , x 2 , x 3 , x 4 , and x 5 are given by 100 km, 112.5 km, 150 km, 187.5 km, and 200 km respectively. As before we preform three separate simulations: One with no remapping (Figure 4a,b), one with low order remapping (Figure 4c,d), and one with higher order remapping (Figure 4e,f). As before both the low and higher order remapping schemes are capable of remap-420 ping the initial particle distribution, with the higher order scheme producing much less numerical diffusion. Figures 4d,f also clearly show that the thickness tracer preserves monotonicity as well as the concentration and volume, clearly demonstrating compatibility of the method.

Two dimensional uniform motion
The next test case we examine involves a two dimensional domain of size 300 km in the x direction and 100 km size in the 425 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) like 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 Figure 1 430 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  where r = (x − x 0 ) 2 + (y − x 0 ) 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 435 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 440 in the x direction with speed 500 m s −1 and examined after 200 seconds ::::::: 100 km :: of ::::::::: translation. As before we perform three simulations, one with no remapping, one with low order remapping and one with higher order remapping, with remapping performed every second ::::: 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 Figure 5. The higher order remapping scheme is able to preserve the approximate shape of the initial concentration distributions while high numerical diffusion with the low order remapping scheme leaves the final 445 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 Figure 6. As can be seen almost no trace of the slot in the slotted cylinder is present for the low order scheme, while the slot is preserved well for the higher order method.
Next we study the error scaling properties of the remapping method. To do this we create various initial element distribu-450 tions 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 0 i is the ice concentration for no remapping. In Figure 7, the values of the L 2 error norm are plotted against the mean 455 radius of the element distribution for both the cosine bell and slotted cylinder initial distributions and the low and higher order remapping methods. As can be seen from the figure the higher order method has significantly lower errors than the low order 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 second order accurate for the higher order method.
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 'L' shaped island 485 consisting of fixed coastal elements. Initial elements are coastal or sea-ice elements according to if (x > x 0 and x <= x 1 and y > y 0 and y <= y 1 ) and not (x > x 0 and x <= x 2 and y > y 0 and y <= y 2 ) then    where (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 non-coastal elements having a concentration of one 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 open water wake behind it.

505
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 one 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, again, is stationary but the wind field, U a , takes the form of a vortex with U a (r) = min ω|r|, where r is the position vector relative to the center of the vortex at (250 km, 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 days and then the element distribution is remapped. As for the island test case, Figure 12 shows that the flux correction scheme is capable of removing excess concentration caused by more realistic non-uniform motion of the elements for this test case.

515
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 is 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 geometricbased remapping algorithm designed to address the unique challenges of accurately remapping sea ice tracer fields for circular 520 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, bounds preserving, 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 525 the method inherits the properties of the incremental remapping algorithms. In the case of non-uniform motion, where gaps and overlaps occur between the deformed effective element areas, monotonicity can be lost. To address this, we developed a novel flux-based correction method, which maintains conservation while correcting for bounds violations due to overlapping elements. Numerical examples are provided that demonstrate the second-order accuracy of the method, bounds preservation for both uniform and non-uniform motion, and compatibility between primitive and conserved tracer quantities.

530
Code availability. The code and test cases used in this paper can be found at https://doi.org/10.5281/zenodo.5800083.