High-resolution regional sea-ice model based on the discrete element method with boundary conditions from a large-scale model for ice drift

Abstract

Understanding sea-ice dynamics at the floe scale is crucial to comprehend polar climate systems. While continuum models are commonly used to simulate large-scale sea-ice dynamics, they have limitations in accurately representing sea-ice behaviour at small scales. DEMs, on the other hand, are well-suited for modelling the behaviour of individual ice floes but face limitations due to computational constraints. To address the limitations of both approaches while combining their strengths, we explored the feasibility of using a DEM within a continuum model, where the latter provides boundary conditions for a rectangular high-resolution DEM domain. This paper presents a feasibility study where a discrete model of a 100 × 100 km² icefield was created using high-resolution optical satellite imagery. Sea-ice dynamics were simulated in the DEM considering environmental forces and integrating large-scale ice-drift velocities as boundary conditions. Model predictions were compared with satellite observations for ice drift and deformation parameters. This numerical approach showed potential for offering accurate, high-resolution predictions of sea ice, particularly in coastal areas and near islands, and may find applications in ice navigation and climate studies. However, further development of the DEM, along with upgrades to the coupled ocean models providing data for the ice component, may be necessary. Additionally, challenges remain to develop a two-way coupling between the DEM and a continuum model, which may be needed to improve the accuracy of large-scale simulations.

1. Introduction

Both climate studies and today's operational sea-ice forecast require information about sea ice at spatial and temporal scales beyond those resolved by common global and regional sea-ice modelling tools. Most of these tools are based on the viscous-plastic (VP) sea-ice model (Hibler, 1979) or its extensions (e.g. Hunke and Dukowicz, 1997). The resolution of the original VP model was on the order of 100 km, below which VP's continuum assumptions may not hold. However, because of deformation processes, various important features, such as leads and pressure ridges, may form in sea ice at much smaller scales. Observations using high-resolution synthetic aperture radar imagery show that nearly all deformation in the Arctic ice cover may be concentrated along so-called linear kinematic features (LKFs), which refer to narrow zones where clusters of leads, ridges and other deformation features of sea ice may be found (Kwok, 2001). To reproduce LKFs in simulations using continuum models that build upon the VP approach, various model modifications have been suggested, which mainly involve using alternative rheology formulations (e.g. Hutchings and others, 2005; Wilchinsky and Feltham, 2006) and/or refining model grids below 10 km or so (e.g. Wang and Wang, 2009; Hutter and Losch, 2020). Notably different from other approaches, new rheological frameworks have also been introduced, for example, the elasto-brittle and Maxwell-elasto-brittle rheologies (Girard and others, 2011; Dansereau and others, 2016). The former has been integrated into neXtSIM and neXtSIM-F (Bouillon and Rampal, 2015; Rampal and others, 2016; Williams and others, 2021), which are sea-ice models running on a Lagrangian (adaptive) grid, as opposed to the traditional VP models based on fixed Eulerian (fixed) grids.

As found from recent comparative studies, all state-of-the-art continuum models have the potential to reproduce some statistical properties of sea-ice deformations observed from satellite imagery at kilometric resolutions (Bouchat and others, 2022). All these models can also produce LKFs, but most of them provide unrealistic LKF distributions (Hutter and others, 2022). Thus, there may still be a need to further increase the resolution of current sea-ice models, which may ultimately require resolving individual floes. Improving high-resolution simulations is important because deformation processes not only create specific narrow features (LKFs) in sea ice but may also influence the sea-ice volume gain in winter through thermodynamic growth in leads and subsequent mechanical redistribution of ice mass under converging conditions (see e.g. Itkin and others, 2018). Furthermore, it is estimated that leads, covering only 1–2% of the central Arctic during winter, may account for more than 70% of the upward heat fluxes in this region, and narrow leads (several metres wide) may be more than twice as efficient in transmitting turbulent heat than those that measure several hundreds of metres (Marcq and Weiss, 2012).

Despite their efficiency on large scales, continuum models cannot accurately resolve elements of sea ice on scales ranging from 10 m to a few kilometres, including floes, leads and ridges, which are relevant for both tactical navigation and climate studies. Today, only DEMs can offer such high resolution. These models consider complex multi-body (floe–floe) interactions, which may also include fracture and ridging processes. The first examples of DEM applications for sea-ice modelling date back to the early 1990s (e.g. Hopkins and Hibler, 1991; Løset, 1994) and represent sea ice as an assembly of 2-D discs in a rectangular control area. It was challenging to realistically model the floe shape as its level of detail directly affected the computational cost of DEM simulations. In general, DEM sea-ice models are computationally demanding, and therefore they typically have been applied to study the interaction of ice floes and ridges with ships and structures at relatively small scales. For these problems, various successful models (including realistic 3-D representations) have been developed, as reviewed by Tuhkuri and Polojärvi (2018). Large-scale discrete models of sea ice have also been considered but often under idealised conditions and with various floe-shape simplifications, for example, using large 2-D floes produced by a randomised Voronoi tessellation (Wilchinsky and others, 2010; Kulchitsky and others, 2016), concave-consolidated floes (Manucharyan and Montemuro, 2022) and circular discrete elements for computational efficiency (Herman, 2016; Turner and others, 2022). However, most of these models have only been used in idealised experiments.

In conventional DEM formulations, the interaction forces are solved by using models mimicking the effect of contact deformation, and the simulations are explicit and advance in very short time steps (Tuhkuri and Polojärvi, 2018). Recently, implicit, non-smooth DEM formulations have also been developed (van den Berg and others, 2018). The major difference between the conventional and non-smooth DEMs lies in the calculation of collision responses, where the non-smooth DEM is formulated by velocities and impulses whereas the conventional method is formulated upon accelerations and forces. Due to this difference, non-smooth DEMs can be stable at relatively large time steps, which broadens the applicability range of the discrete approach, making it potentially suitable for regional sea-ice modelling.

Over the past few decades, there has been a movement towards modularity in Earth System Models, which enables combining models developed for different physical processes and scales. This has resulted in integrated systems with various nested components in the next generation of regional sea-ice models, in which a large-scale, continuum model provides boundary conditions for an embedded, higher-resolution domain (e.g. Duarte and others, 2022). However, these embedded ice models are typically continuum models run with horizontal resolutions exceeding 2 km, which is insufficient to resolve floe–floe interactions. In the sea-ice modelling community, it has been speculated that a DEM may be applied near the ice edge or coastlines within a coarse-resolution continuum model, to better capture floe-scale effects at the high resolutions and short timescales needed for navigation (e.g. Hunke and others, 2020). Thereby, the large-scale coverage of the continuum model could be combined with the high local resolution of the DEM. We have further developed this idea through a feasibility study of regional sea-ice modelling involving a nested non-smooth DEM and a continuum framework. Embedding an efficient DEM into a regional sea-ice model may be particularly valuable, for instance, in planning shipping and commercial activities in the marginal ice zone (MIZ). This model may also provide insights into the dynamics of leads and pressure zones in sea ice.

This paper is structured as follows. Section 2 outlines both the DEM and continuum frameworks for modelling sea-ice dynamics and explains the concept of embedding DEM within a larger-scale system. Section 3 describes the data required for our simulations, including atmospheric and ocean forcing data, satellite imagery and methods for reconstructing model ice fields. Here, Section 3.4 details the simulation procedures followed in the presented analysis and describes the boundary conditions used in different simulations. Section 4 focuses on the results of our feasibility study, where a discrete model of a 100 × 100 km² icefield was created using high-resolution optical satellite imagery of sea ice near Svalbard. Sea-ice dynamics were simulated in the DEM considering environmental forces and integrating large-scale ice-drift velocities as boundary conditions. Model predictions are compared with satellite observations for ice drift and deformation parameters. Simulations with different model parameters and boundary conditions are considered to study the variability of predicted ice drift and model errors. Section 5 discusses the results presented in Section 4, focusing on the localisation of ice deformation, the importance of tidal current and ice strength in the simulations. Finally, the conclusive remarks in Section 6 highlight the advantages and shortcomings of the presented approach, suggesting possible directions for future research.

2. Numerical approach

2.1 Sea-ice dynamics in a continuum framework

Traditional continuum models for sea-ice dynamics solve the momentum equation for ice mass per unit area (m) as it flows through a fixed Eulerian mesh driven by winds and ocean currents. Here we used CICE, the Los Alamos Sea Ice Model employed by several Earth System Models (CICE Consortium, 2021). It solves the following 2-D momentum equation obtained by integrating the 3-D equation in the vertical direction through ice thickness (Hunke and Dukowicz, 1997):

(1)$$m\displaystyle{{\partial {\boldsymbol u}} \over {\partial t}} = \nabla \cdot \sigma + {\boldsymbol \tau }_{\rm a} + {\boldsymbol \tau }_{\rm w}-\hat{{\boldsymbol k}} \times mf\,{\boldsymbol u}-mg\nabla H_{\rm o}$$

Here, u is the ice velocity, σ _ij is the internal stress tensor and τ_a and τ_w are the wind and ocean stresses, respectively. The two last terms on the right-hand side are stresses due to Coriolis effects and the sea surface slope. The ocean stress, which is tangential to the ice, is given by

(2)$${\boldsymbol \tau }_{\rm w} = C_{\rm w}\rho _{\rm w}\vert {{\boldsymbol U}_{\rm w}-{\boldsymbol u}} \vert [ {( {{\boldsymbol U}_{\rm w}-{\boldsymbol u}} ) \cos \theta + \hat{{\boldsymbol k}} \times ( {{\boldsymbol U}_{\rm w}-{\boldsymbol u}} ) \sin \theta } ] $$

where ρ _w, C _w and U_w are water density, drag coefficient and the current velocity, respectively; $\hat{{\boldsymbol k}}$ is the vertical unit vector and θ is the turning angle between geostrophic and surface currents. The turning angle is necessary if the top ocean model layers are not able to resolve the Ekman spiral in the boundary layer. If the top layer is sufficiently thin compared to the typical depth of the Ekman spiral, then θ = 0 is a good approximation. Here, as in CICE Consortium (2021), we assume that the top layer is thin enough. Similarly, the wind stress can be expressed via surface wind velocity:

(3)$${\boldsymbol \tau }_{\rm a} = C_{\rm a}\rho _{\rm a}\vert {{\boldsymbol U}_{\rm a}} \vert {\boldsymbol U}_{\rm a}$$

where ρ _a and C _a are air density and drag coefficient, respectively; and U_a is the wind velocity at 10 m.

CICE discretises the momentum equation in time and then solves it, utilising, for example, the elastic-viscous-plastic ice dynamics scheme (Hunke and Dukowicz, 1997). In this scheme, σ _ij is determined from a regularised version of the VP constitutive law, which treats the ice pack as a visco-plastic material that flows plastically under typical stress conditions but behaves as a linear viscous fluid where strain rates are small. At each integration time step (typically, every hour or so), CICE finds velocity u and provides various data including but not limited to ice-drift direction, velocity, concentration and thickness. These parameters can serve as initial and boundary conditions for a nested DEM, as described in Section 2.3.

In our implementation, we used CICE with a mesh of 1° as illustrated in Figure 2. This resolution suffices for integration with the 100 × 100 km² DEM domain considered here.

2.2 Discrete element model

The utilised DEM solves the 3-D rigid-body equations of motion with the same force components as on the right-hand side of Eqn (1). However, since the icefield is now discrete (consisting of convex polyhedra), contact forces between rigid floes need to be modelled instead of the internal stress tensor in Eqn (1), and the fluid forces need to be integrated over the surfaces of the floes. Further, the force term due to the sea surface slope is neglected, as its effect is typically minor compared to the remaining terms (Leppäranta, 2011).

We employed a so-called non-smooth DEM formulated on the level of velocities and impulses, which contrasts with the traditional models formulated upon acceleration and forces. The comparative advantage of the current DEM is its efficiency in resolving a large number of contacts between interacting floes which may involve Coulomb friction and ice crushing (van den Berg and others, 2018). Various analytical solutions for ice fracture due to interaction with icebreaking vessels have been implemented in this DEM. These solutions have been used to model icebreaker performance in ice management operations (e.g. Lubbad and others, 2018). However, we cannot apply these analytical solutions for sea-ice fracture in our large-scale model, since they have been developed specifically for ice–structure interactions, which typically occur within the length of a ship, and may not be valid at larger scales. Thus, in our DEM simulations here, the ice floes cannot break.

Sea ice often undergoes fracturing processes, wherein ice floes crack or break into fragments or pieces. This can result from a variety of forces, such as tension, compression or shear forces. Although our present model cannot accurately describe the physics of ice fracture, it can be parameterised to predict local deformations of an icefield due to, for example, floe crushing and rafting. In this context, crushing refers to a specific type of fracturing process in which ice is subjected to mechanical compression that breaks it into smaller pieces or causes it to deform. For instance, when two ice floes interact, they may compact or shatter under pressure, resulting in the formation of small ice fragments at the contact zone. While our model does not simulate these small ice fragments, it allows the two floes to penetrate each other within the crushed volume and restore their shapes upon separation. The overlap volume of two interacting ice floes in the numerical simulation represents crushed ice. Figure 1 helps to illustrate this assumption.

Figure 1. Overlap volume of two interacting floes represents crushed ice in the model of van den Berg and others (2018). For clarity, a 2-D sketch is shown here, but the algorithm is implemented fully in 3-D.

The contact model for floe–floe interaction is developed by van den Berg and others (2018) based on the assumption that local ice crushing will occur. In this model, ice crushing is represented by floe overlap (Fig. 1), and the corresponding restoring forces are calculated using exact contact geometries and material properties of interacting floes. Supported by experiments, van den Berg and others (2018) assumed constant energy dissipation per crushed volume or crushing specific energy (CSE). The latter is the amount of energy needed to crush a unit volume of ice, which can be equivalently expressed in J m⁻³ or Pa. The assumption of a constant CSE is equivalent to a constant crushing pressure (or plastic limit stress) during indentation, which is used to calculate the contact force when crushing starts:

(4)$$F_{{\rm cr}} = A_{{\rm proj}}\cdot {\rm CSE}$$

where F _cr is the limit of contact force above which crushing occurs, A _proj is the contact projected area and CSE is the crushing specific energy of sea ice. Below F _cr, the contact force is calculated such that the floes' relative velocity in the normal direction of their contact interface remains zero, or the force is zero in the case of separating floes. Additionally, there are tangential contact forces modelled according to the Coulomb friction law (for details, see van den Berg and others, 2018).

Using this contact model and replacing CSE in Eqn (4) with a tuning parameter for floe ice strength, which may be different from its crushing specific energy, allows us to approximately model various scenarios involving floe rafting, ridging and crushing as plastic deformations of overlapping floes in our DEM. The values of this and other parameters used to model ice dynamics in this paper are given in Table 1.

Table 1. Standard parameters used to model ice dynamics

2.3 Embedding DEM inside a larger sea-ice model

In our current approach, the embedded DEM was implemented to run in parallel with CICE. However, instead of the latter, any alternative tool for predicting large-scale sea-ice motion may be used, as there is currently no coupling between the models. Here, the embedded DEM domain is initially square with four linear boundaries enclosing discrete ice floes. The boundaries are rigid but can move independently from each other according to interpolated ice-drift velocities predicted by the continuum model (or other tools) in its own computational mesh nodes (Fig. 2). As these boundaries drift, they may confine the floes inside, such that no flux of ice through the boundaries is allowed. Thereby, the DEM domain can deform from its initial state while moving over the fixed-continuum model mesh. Note that there is no feedback passed from the DEM into the continuum model in our current approach.

Figure 2. Fixed grid of the continuum model showing in colour an example of ice velocities predicted by CICE (left) and the nested boundaries of the DEM domain northeast of Svalbard with a close-up (right) illustrating how the DEM boundaries can move over the fixed grid. Two islands, Kvit and Victoria, may be located within the DEM domain as it moves and possibly deforms.

The dynamics within the DEM domain are computed based on atmosphere and ocean forcing data (as in Section 3.3), which may have higher temporal and spatial resolutions than the forcing data of the continuum model. The domain size of the adopted DEM in our applications is ~100 km, and the timescales can range from a fraction of a second up to a few days. The resolution of input data must be sufficient to resolve studied processes at these scales. As the duration of our simulations here is ~25 h, no thermodynamic processes are modelled within the DEM domain. Further specific details on the forcing data, boundary conditions, initialisation and execution of the model used in our analysis are given in Section 3.

3. Simulation scenarios and available data

3.1 Satellite imagery of the study site

The study site, which is the region where sea-ice drift is modelled and analysed in this paper, was chosen primarily based on data availability, including satellite optical imagery. We chose the Barents Sea because we also have our own observations and experience from various field experiments in this region (e.g. Tsarau and others, 2017). An interesting region was found using Sentinel-2 satellite imagery, covering the area northeast of Svalbard between Kvit Island and Victoria Island (see Fig. 2) with three successive images taken between 12.06 UTC on 15 April and 13.17 UTC on 16 April 2016 (Fig. 3).

Figure 3. Sentinel-2 optical imagery of the same area northeast of Svalbard at 12.06 and 13.47 UTC on 15 April and at 13.17 UTC on 16 April 2016 (left to right) with Kvit Island and Victoria Island highlighted. The black triangular regions in the southern parts of the images indicate missing data or the presence of clouds.

Source: https://dataspace.copernicus.eu/

Each of the images in Figure 3 depicts the same geographic region of ~100 × 100 km². The GPS coordinates of this square region in decimal degrees are [79.90, 32.18] (bottom left) and [81.01, 38.92] (top right). The images were captured by Sentinel-2 sensors in a visible spectral band with a spatial resolution of 10 m. More details on the Sentinel-2 mission, including sensor characteristics and collected data, are provided by Drusch and others (2012). In this study, it is important that the optical images are nearly cloud-free, enabling the identification of distinct ice floes and island boundaries. The presence of the islands shown in Figure 3 makes the considered case particularly interesting for our modelling because landmasses can cause deformation in pack ice, for example, due to shear, possibly involving breaking, and opening and closing of cracks. The localisation of deformation is an important aspect in our analysis.

3.2 Image analysis

Image processing techniques were employed in this paper to provide our DEM with an initial configuration of ice floes corresponding to the icefield depicted in the first image in Figure 3 (taken at 12.06 UTC on 15 April 2016). We used a procedure that involves the identification of greyscale thresholds for ice floes, image segmentation resulting in regions of ice and water, and subsequently processing the identified floe boundaries and shapes. These steps are commonly employed in image analysis methods used for automatic detection of ice floes (see e.g. Zhang and others, 2013; Zhang and Skjetne, 2015). However, based on our experience, not all these methods may be equally effective when processing large, high-resolution images of a densely packed icefield (as in our case). They may also require manual intervention when floe boundaries are poorly identified, which is typical for highly concentrated icefields.

We achieved a balance between accuracy and efficiency by employing the well-known watershed transform for image segmentation, as inspired by Zhang and others (2013), preceded by the distance transform of the binary (ice-water) image. To prevent over-segmentation, we filtered out any shallow minima that might have resulted from the distance transform before applying the watershed segmentation algorithm. We further applied automatic geometric operations, such as convex hull, trimming and union, to the resulting segments, creating connected convex polygons that approximate the ice floes seen in Figure 3. This procedure simplifies any meanderings of the boundaries (see Fig. 4), preserving the overall shape of the floe, as long as it is nearly convex. Highly irregular, non-convex floes are divided into convex polygons to facilitate DEM modelling. Thus, the result of this image processing can be readily used in our DEM, as demonstrated in Figure 5. Here, a total of 10 737 ice floe with sizes ranging from 50 m to 16 km was modelled.

Figure 4. Result of image segmentation (left) and a close-up view of the ice floes approximated by convex polygons (right).

Figure 5. Icefield in the DEM represented by a total of 10 737 ice floes. The islands are shown in red.

Ice-drift fields can also be extracted from satellite imagery using feature tracking techniques. For instance, Itkin and others (2018) used ridge sails in sea-ice surface topography data obtained by an airborne laser scanner as virtual buoys to identify displacements and drift velocities. Here, we used distinct features, for example, small floes and leads appearing in successive images in Figure 3, as markers and matched them across the images, employing the Computer Vision Toolbox in MATLAB R2022a, to retrieve displacements. These displacements were further converted into drift vectors, which served as the basis for reconstructing the sea-ice drift fields. The identified vector fields of ice motion are presented in Section 4.1.

3.3 Data sources

Optical satellite imagery does not provide data in the ice-thickness space, which is also required to initialise 3-D floes (polyhedra) in the DEM. Here, we used ice-thickness data from CS2–SMOS product (Ricker and others, 2017). This product combines thickness estimates from the European Space Agency CryoSat-2 (CS2) altimeter for thick ice (Ricker and others, 2014) and from the Soil Moisture and Ocean Salinity satellite for thin ice (Tian-Kunze and others, 2014), providing data for Arctic-wide studies on the entire thickness range. In our region of interest, the weekly average ice thickness provided by CS2–SMOS was ~0.3–0.55 m for the period from 9 April until 15 April 2016. We found rafted floe ice in a similar thickness range in the MIZ further south (downstream) from this area 2–3 weeks later (Tsarau and others, 2017), supporting the validity of the above estimates. Unfortunately, CS2–SMOS does not provide data after 15 April (presumably, due to the melt season). Therefore, we used the ice thickness data from the week before but covering a larger region (extending towards the north, from where the ice was drifting). The area-weighted average ice thickness in this region was ~0.45 m, which was used as the floe thickness in the DEM (the same for all floes). If high-resolution ice-thickness data were available on 15 April 2016 and onwards, the model could be initialised with a more realistic thickness distribution.

After initialising an icefield in the nested model, the force terms in Eqn (1) need to be provided. For wind forcing, we utilised the 10 m wind data from the NORA3 reanalysis, a high-resolution numerical mesoscale weather simulation for the North Sea, the Norwegian Sea and parts of the Barents Sea downscaled from the state-of-the-art ERA5 reanalysis (Solbrekke and others, 2021). Hourly wind speeds and directions with 3 km horizontal resolution were obtained for our DEM.

We also used the Hybrid Coordinate Ocean Model (HYCOM), which provides 41-layer ocean data in a 1/12° global grid with an average spacing of 6.5 km (Chassignet and others, 2009). The temporal resolution of the HYCOM non-tidal surface current is 3 h; however, we interpolated it and added a tidal component to update our DEM with hourly estimates of the total surface current. This was done due to the lack of hourly current data from 2016 including both tidal and non-tidal components resolved in the upper ocean layer.

The tidal component was predicted using the Arctic 2 km Tide Model (Arc2kmTM) developed by Howard and Padman (2021). This is a barotropic ocean tide model on a 2 × 2 km² polar stereographic grid, developed using the Regional Ocean Modelling System. Arc2kmTM consists of spatial grids of complex amplitude coefficients for sea surface height and depth-integrated currents for eight principal tidal constituents: four semidiurnal (M2, S2, K2, N2) and four diurnal (K1, O1, P1, Q1). Using the Arc2kmTM model and its bathymetry data in the relevant region (Fig. 6), we obtained depth-averaged tidal currents, which were simply added to the HYCOM current to obtain an approximation of the total current (Fig. 7). The mean values of both current components (together with the wind data) are summarised in Table 2.

Figure 6. Bathymetry map overlaid on top of Sentinel-2 optical imagery of the icefield possibly indicating more open ice in shallow water. The satellite image was taken at 12.06 UTC on 15 April 2016.

Figure 7. Scaled vectors of wind (left) and total current (right) from 12.06 until 13.47 UTC on 15 April (blue) and from 13.47 UTC on 15 April until 13.17 UTC on 16 April 2016 (red) with Sentinel-2 imagery of the area in the background.

Table 2. Mean wind and current data

All directions are measured from north, positive clockwise (e.g. drift towards the east is 90°).

3.4 Simulation procedures

The developed nested model was executed in different testing modes, which differed mainly in the way the DEM boundary conditions were provided. As presented in Section 2, a trial set-up of the nested model was developed to run in parallel with CICE; however, any alternative tool for predicting large-scale sea-ice motion may be used. We considered this trial set-up, running both models, CICE and DEM, locally to make sure that the nested DEM concept is feasible. A standard version of CICE 6.4.0 on a 1° global grid was utilised here following the preparation and execution steps described by CICE Consortium (2021). The DEM set-up was as described in Section 2.2, with its own wind and ocean data as presented in Section 3.3. Both models have their own numerical integration schemes with different time stepping. A typical time step in CICE is ~1 h, while the DEM solves sea-ice dynamics at 0.2 s intervals. The larger of these two time steps was used as a communication interval to update the boundary conditions in the DEM with CICE data. Communication with the DEM was implemented via the TCP/IP network interface, allowing data exchange between different machines on the network. In this simulation, we were only sending data from CICE to the DEM; however, it may also be possible updating some data fields in CICE based on DEM predictions (e.g. information about ice concentration, leads and floe-size distribution) for coupling. Fully coupled simulations may be considered in future. Since CICE with its large time step may run much faster than the DEM, the predicted continuum ice velocities were stored in a buffer. From there, the DEM TCP/IP module was reading the hourly data, interpolating it to the DEM boundaries, and then updating the DEM boundary conditions. This simulation is depicted in Figure 2, where the DEM domain moves over the 1° grid with velocities predicted by CICE. Here, the centres of the rigid DEM boundaries move with ice velocities in the nearest nodes of CICE.

To compare the DEM performance against observed data and investigate the effect of boundary conditions, we considered another simulation procedure with prescribed boundary velocities. To avoid unnecessary efforts validating CICE simulations against satellite observations and rather focus on the performance of the DEM, we obtained possible sets of boundary conditions directly from the available satellite images. This was done based on the same feature tracking technique described in Section 4.1. The results are presented in Table 3 providing time-dependent (N1) and constant (N2) boundary velocities, which are the same for the northern, southern, eastern and western boundaries in both cases here. Additionally, we implemented a set-up with no confinement along the boundaries (N3), which means no internal force was prescribed on the floes at the boundaries. We used different boundary conditions in our numerical experiments to investigate how the accuracy and different update frequencies (constant vs variable) of boundary velocities impact the accuracy of DEM simulations. Among the sets of boundary conditions presented in Table 1, N1 is considered the most accurate, while N3 is an extreme simplification.

Table 3. Boundary conditions used in the DEM simulations

In N1 and N2, ice-drift velocities obtained from Sentinel-2 images were prescribed on the rigid DEM boundaries (the same for the northern, southern, eastern and western boundaries here). In N3, no rigid boundaries were modelled representing a scenario with no confinement.

These sets of boundary conditions (Table 3) were used to update our DEM in the simulations presented in Section 4. Otherwise, the procedure for executing the model was the same as described above for concurrent CICE-DEM simulations. Standard (or rather typical) parameters for modelling ice dynamics were used, as listed in Table 1.

In all our simulations here, the ice floes were generated using the satellite image taken at 12.06 UTC on 15 April (in Fig. 3) as described in Section 3.2. For practical reasons, all floes were initialised with the velocities of the current in the DEM. In some scenarios, more accurate initial conditions could potentially be obtained, for example, by interpolating the DEM boundary velocities. In simulations with CICE, continuum velocity predictions could be used. However, we did not explore these options.

4. Results

This section provides a concise presentation of our modelling results, focusing on ice drift (Section 4.1) and deformation rates (Section 4.2), and compares them with the corresponding parameters observed from satellite images. Additionally, Section 4.3 presents supplementary results from a sensitivity study, where we varied certain model parameters and explored different boundary conditions. These results offer insights into how different model parameters may affect ice drift and deformation but are not intended for a thorough quantitative analysis. Furthermore, Section 4.4 offers rough estimates of uncertainties in both the input data for the model and the observations. All data and simulation results presented here will be further discussed in Section 5.

4.1 Ice drift

Ice-drift fields were extracted from the satellite imagery in Figure 3 using feature tracking techniques as described in Section 3.2. The results are depicted in Figures 8, 9, showing the mean ice drift within two temporal intervals defined by the time stamps in Figure 3: from 12.06 until 13.47 UTC on 15 April, and 13.47 UTC on 15 April until 13.17 UTC on 16 April 2016. For brevity, we will refer to all average parameters identified within these two intervals as the ‘2 h’ and ‘24 h’ means, respectively. All graphical results presented in this paper were plotted in the same reference frame as the original satellite images shown in Figure 3, to facilitate comparison between the model and observations.

Figure 8. Ice-drift vectors between 12.06 and 13.47 UTC on 15 April 2016 obtained from satellite imagery (left) and the DEM (right).

Figure 9. Ice-drift vectors between 13.47 UTC on 15 April and 13.17 UTC on 16 April 2016 obtained from satellite imagery (left) and the DEM (right).

As the positions of all ice floes in the DEM can be tracked at any time, the model ice-drift vector fields are readily available. The results obtained by executing our model with its standard input parameters (Table 1) and using the time-dependent boundary conditions N1 from Table 3 are plotted in Figures 8, 9 for both temporal intervals considered. Additionally, we calculated the spatial means for both the model and the satellite data, as summarised in Table 4. As seen, these average drift speed estimates are reasonably captured by the model, deviating from the satellite data by no more than 20%. However, there are greater differences when considering all data points separately and computing their RMSE, as presented in Section 4.4. It can also be observed (e.g. in Fig. 9) that there is a substantial spatial variation in the predicted drift vectors, with seemingly more accurate predictions in the central part of the model domain.

Table 4. Mean ice drift and deformation rates with percentage difference between the predicted and observed parameters in parentheses

All directions are measured from north and are positive clockwise (e.g. drift towards east is 90°).

When comparing model results with observational data, it is important to also consider the observational error. In Section 4.4, we estimated a possible error of 10% or so in the ice-drift velocities extracted from the satellite images and used as observational data here. Sections 4.3 and 4.4 provide also results that help to understand the effect of boundary conditions on the predicted drift speeds.

As Table 4 shows, the predicted drift directions are accurately captured only when considering 24 h means. The short term, 2 h, drift directions in the model deviate significantly from the satellite data, by up to 10–20°. Possible reasons for this discrepancy are investigated and discussed in Section 5.

4.2 Deformation rates

Deformation rates were derived from the ice motion data obtained in Section 4.1 and described by a velocity field (u, v) in two dimensions as follows (e.g. Kwok, 2001; 2005; Bouchat and others, 2022):

(5)$$div = \displaystyle{{\partial u} \over {\partial x}} + \displaystyle{{\partial v} \over {\partial y}}$$

(6)$$shr = \sqrt {{\left({\displaystyle{{\partial u} \over {\partial x}}-\displaystyle{{\partial v} \over {\partial y}}} \right)}^2 + {\left({\displaystyle{{\partial u} \over {\partial y}} + \displaystyle{{\partial v} \over {\partial x}}} \right)}^2} $$

(7)$$D = \sqrt {shr^2 + div^2} $$

where ∂u/∂x, ∂u/∂y, ∂v/∂x and ∂v/∂y are the spatial gradients in ice motion computed on a regular, square grid (with a side of ~5 km here); and div, shr and D are the divergence, shear and total deformation rates, respectively.

The 2 h deformation rates, derived from both the satellite data and the model, are presented in Figure 10 next to each other; and the 24 h means are presented in Figure 11. Additionally, the area-averaged values of D, shr and drift speeds are tabulated in Table 4 to facilitate the comparison. All model results presented here were obtained by executing the DEM with its standard parameters (Table 1) and the time-dependent boundary conditions (N1 in Table 3), while results with alternative model configurations will be presented in Section 4.3.

Figure 10. Total deformation, shear and divergence between 12.06 and 13.47 UTC on 15 April 2016 obtained from satellite imagery (top row) and the DEM (bottom row). Because of missing data in the southern region of the satellite image captured at 13.47 UTC on 15 April, observational drift vectors are not depicted in this area.

Figure 11. Total deformation, shear and divergence between 13.47 UTC on 15 April and 13.17 UTC on 16 April obtained from satellite imagery (top row) and the DEM (bottom row). Because of missing data in the southern region of the satellite image captured at 13.47 UTC on 15 April, observational drift vectors are not depicted in this area.

The deformation rates obtained from both the model and the satellite imagery are of the same order of magnitude, and both datasets demonstrate that sea-ice deformation depends on temporal scales. For instance, Table 4 shows that the mean hourly deformation rates (‘2 h means’) are roughly twice as large as their daily (24 h) counterparts. This observation appears consistent with the results of research on the temporal scaling of sea-ice deformation, for example, Hutchings and others (2011), where a log–log linear scaling relationship was found between deformation and timescale, T: log (D) ∝ log (1/T).

Further, both the model and satellite imagery demonstrate that ice deformation is localised near the landmasses, with the highest deformation rates observed within a distance of 20–40 km from the islands (see Figs 10, 11). Despite all these apparent similarities between the model and the satellite data, the DEM with its standard model parameters tends to substantially overestimate (by approximately a factor of two) the average values of total deformation, which seems to be shear dominated in this case. We discuss these results and possible strategies to improve the modelling of ice deformation in our DEM in the following sections.

4.3 Alternative model parameters

As presented in Section 4.1, running our DEM with its standard parameters (as listed in Table 1) and time-dependent boundary conditions (Table 3) yielded reasonably accurate results in terms of ice-drift speeds. However, when compared with the satellite data, the predicted deformation rates were not highly accurate (Section 4.2). To gain a better understanding of the sensitivity of ice deformation to various model parameters and input data, we considered alternative DEM configurations.

In two alternative simulations, we reduced CSE from the standard 2 MPa, typically used to model sea-ice crushing, to values of 200 and 20 kPa. These adjustments were made to account for reduced strengths resulting from possible ice fracture and ridging processes. In another simulation, we focused on the tidal component, which is considered one of the most uncertain forces in the ice dynamics model. In this simulation, we excluded the tidal component from the total current model (see Table 2). In two additional simulations, we considered different boundary conditions described in Table 3: boundaries with a constant velocity (N2) and a set-up with no confinement along the boundaries (N3). These simulations allowed us to explore how boundary conditions may affect the DEM results.

The results obtained with these alternative model configurations include 2 h mean ice-drift parameters and deformation rates, as summarised in Table 5. Additionally, selected graphical data are presented in Figures 12 through 15. We briefly describe the main findings and our interpretation of these results below.

Table 5. Mean 2 h ice drift and deformation rates predicted with different DEM parameters and model configurations

Figure 12. Total deformation, shear and divergence between 12.06 and 13.47 UTC on 15 April 2016, as predicted by the DEM without the tidal-current component.

Excluding the tidal current has the strongest effect on the prediction of short-term drift speeds and deformation rates in our analysis. This can be seen when comparing the results in Table 5. Additionally, Figure 12 shows that excluding the tidal component leads to considerably less deformation concentrated around the islands (where the water is shallow) than in Figure 10 predicted with the tide. We delve into the analysis of these results and their implications in Section 5.2.

The effect of CSE appears relatively less significant than, for example, the tidal effects in our simulations, except in the case with the lowest ice strength considered here (CSE = 20 kPa). Too low ice strength in the model renders the predicted ice drift almost insensitive to the presence of landmasses. This is evidenced by the drift vectors in Figure 13 predicted with CSE = 20 kPa, which show nearly straight trajectories close to the islands. However, both the observational data and simulations with much stronger ice (e.g. CSE = 2 MPa) suggest that the presence of the islands strongly affects the ice drift, as the ice trajectories gradually bend around the islands at a distance from them. This is shown by the drift vectors in Figure 10. The differences between the drift vectors in Figures 10, 13 arise mainly near Kvit Island, where the observed drift vectors in Figure 10 tend to align with the boundaries of the island, whereas the predicted vectors in Figure 13 show nearly uniform velocities. Considering this behaviour, the drift vectors predicted at CSE = 2 MPa appear more realistic than those at CSE = 20 kPa. Additional discussions on the role of this ice strength parameter in our model are provided in Section 5.3.

Figure 13. Total deformation, shear and divergence between 12.06 and 13.47 UTC on 15 April 2016, as predicted by the DEM with CSE = 20 kPa.

Appropriate boundary conditions are important for coupled DEM–continuum model simulations. Here, we experimented with different prescribed boundary conditions in the DEM, without the necessity of implementing different coupling solutions. This allowed us to test our set-up without enduring the computational burden of running extensive coupled simulations. Nevertheless, the insights gained from these experiments may be equally useful for both assessing the performance of our DEM and establishing requirements for boundary conditions in a coupled system, for example, with respect to the accuracy of boundary velocities and their update frequency.

As far as the accuracy of the model boundary velocities is concerned, its effect on the DEM predictions is important, as the motion of the boundaries reflects how other floes surrounding the region affect the dynamics of the floes within the DEM. For example, the high deformation rates shown at the southern boundary of the model domain in Figure 10 (and also at the northern boundary in Fig. 11) are presumably due to the prescribed velocity of the boundaries, which may somewhat deviate from the actual ice velocities in those locations. In contrast, in Figures 14, 15, obtained from the simulation with no confinement along the boundaries (i.e. free boundary conditions with no effect from external floes), no artificial deformation can be seen. Although employing free boundary conditions here represents a significant simplification, comparing this scenario with other simulations helps identify the potential drawbacks associated with introducing artificial boundaries in the DEM. We further analyse the role of boundary condition in our modelling in Section 4.4, where an analysis of the error propagation from the boundaries into the model domain is presented.

Figure 14. Total deformation, shear and divergence between 12.06 and 13.47 UTC on 15 April 2016, as predicted by the DEM without confinement.

Figure 15. 24 h total deformation, shear and divergence between 13.47 UTC on 15 April and 13.17 UTC on 16 April, as predicted by the DEM without confinement.

4.4 Errors and uncertainties in the modelled and observational data

It is difficult to carry out a proper uncertainty analysis of all modelled results presented above because the errors in the input data are unknown. To run our DEM, we used multiple data sources but, unfortunately, we do not have control of their accuracy. Therefore, we can only provide some ideas about the size of possible errors in the wind and ocean forcing data, as described below. Additionally, we estimated the observation error of the ice-drift data extracted from satellite images and evaluated possible ice-drift errors in the DEM boundary conditions. This information about the uncertainties in both the model's input and the observational data may be useful for interpreting the numerical results presented above.

Solbrekke and others (2021) compared several NORA3 wind datasets with observational data collected at sites for possible offshore wind power installations. They found that the model underestimated the wind speed at all sites, and the largest difference between the observed and simulated mean wind speeds was 8.9%.

The accuracy of the HYCOM current dataset used here is unknown. However, utilising these data to approximate the total current at short intervals (e.g. 2 h) without accounting for the tidal current could introduce an uncertainty range of ±40% or so, as can be seen in Table 2. We tried to reduce this uncertainty by also considering the tide, but there were no observational data available to compare with our simulated total current and therefore no estimates of its possible error could be provided.

We were able to estimate the accuracy of the observed ice-drift fields extracted from satellite imagery. This was done based on the known resolution of the used Sentinel-2 images (10 m) and our experience that at least 5 pixels were needed to represent a distinct feature that could be matched with certainty across two images. Thus, assuming a maximum displacement error of 50 m, we could measure the mean drift speed with an accuracy of ±0.01 m s⁻¹ on a 2 h interval. However, this accuracy was probably compromised when interpolating the calculated drift data based on tracking irregularly spaced features into a regular 5 km mesh, which was necessary for comparing drift fields and their analysis in Section 4.1. The errors introduced because of this interpolation could be ~0.02 m s⁻¹, based on the variation of the observed drift data on 5 km intervals. Thus, the errors in the observed drift data presented in Section 4.1 could be of a size of 10% or so.

The variation of the observed ice drift along the 100 km-long boundaries of the DEM was considerable (with RMS estimates of 0.04–0.07 m s⁻¹), implying potential errors of ~20% in the drift speeds used as the DEM boundary conditions (both N1 and N2 in Table 3). It is important to note that providing very accurate ice velocities at the DEM boundaries was not assumed in our approach. However, large errors certainly could negatively affect the simulation results, especially closer to the boundaries, as demonstrated in Table 6. In this table, we computed RMS differences between the observed and simulated drift speeds within the largest possible region located inside the DEM domain at a given minimum distance from its boundaries. For example, more internal regions show smaller prediction errors, which indicates that errors in the boundary conditions may propagate into the model domain but decrease with distance from the boundaries.

Table 6. Mean absolute and relative (in parentheses) differences between the observed and simulated ice-drift speeds at different distances from the model boundaries

5. Discussion

5.1 Localisation of deformation zones

Sea-ice deformation processes, such as compression, tension and shear, may lead to the formation of so-called linear kinematic features (LKFs), including leads and pressure ridges. Remote-sensing techniques suggest that LKFs concentrate along islands and coastlines (e.g. Hutter and Losch, 2020). This is an indication of significant ice deformation in these areas. Our analysis of satellite data also reveals high deformation rates near landmasses (see Figs 10, 11). As presented in Section 4.1, the DEM appears to accurately capture the localisation of these deformation zones. This demonstrates the potential value of using DEMs as nested components in a global or regional continuum model, especially when improved local resolution near coastlines and islands is required, for example, for navigation planning.

5.2 Temporal variation

Temporal variation in ice deformation is evident in the results presented in Section 4.1, where it was found that the mean hourly and daily deformation rates may differ by a factor of 2 or so. In the presented scenario, this temporal variation may be correlated with the forcing data due to ocean currents. Among all model components listed in Table 2, the tidal-current component has the largest relative difference (94%) between the 2 and 24 h intervals, accounting for a major part of the total current variation. This may explain the observed temporal variation in both ice-drift speeds and the deformation rates, implying that the tidal-current component may be of crucial importance for accurately modelling sub-daily ice motion. The effect of tidal currents may be particularly strong in shallow waters, as Figure 6 indicates, showing more open ice in areas of shallow water.

Watkins and others (2023) analysed sea-ice dynamics based on ice-buoy data collected during the MOSAiC Expedition in the Fram Strait. The authors found that, in shallow seas, strong tidal currents affect ice drift, resulting in repeated opening and closing of the ice. According to their analysis, boundary currents increase ice-drift speeds near the shelf edge, causing shear in the ice pack. Based on these observations, the authors speculate that sea-ice models that disregard small-scale ocean currents will underestimate ice deformation.

In our case, the DEM overestimates ice deformation, but it may be due to inaccurate modelling of the surface currents, especially the tidal component (see Section 3.3). Further research is necessary to investigate this issue. Nevertheless, it is clear that to improve sea-ice predictions on short temporal scales (e.g. hourly), it is insufficient to modify the sea-ice model alone. Upgrades to the ocean components of coupled models and their coupling mechanisms are also essential.

5.3 Ice strength in DEM

As discussed above, differences in the 2 h deformation rates between the model and satellite data are expected, as the modelled tidal current may not accurately represent conditions in the upper ocean layer. However, for the 24 h case, there must be other reasons for the overestimated deformation. It can be argued that this DEM, which does not explicitly model ice fracturing, cannot accurately predict long-term, large ice deformations, although small deformations may still be captured on shorter time intervals. For example, the high deformation rates seen in Figure 11 between Victoria Island and the eastern model boundary indicate that large rigid floes might become stuck in that location as the distance between the island and the boundary here is ~20 km. In the satellite images in Figure 3, large floes are seen northeast of Victoria Island as of 12.06 and 13.47 UTC on 15 April but not at 13.17 UTC on 16 April 2016, suggesting that these large floes were likely broken into pieces when pushed against the landmass. However, the same floes in the model could become stuck between the island and the eastern model boundary because this fracture behaviour cannot be modelled.

The scenario shown in Figure 11 was simulated using the crushing strength of ice (CSE = 2 MPa), rendering the entire icefield relatively rigid, which explains the observed behaviour with floes becoming stuck. In Section 4.3, we also made attempts to model more compliant ice fields with strength values of 200 and 20 kPa. Reducing the CSE parameter had some effect on the predicted deformation rates (which were gradually reduced for lower CSE), showing a more plastic behaviour of the icefield. Thus, CSE may potentially be tuned in the employed DEM to model a reduced strength due to, for example, ridging and rafting processes, but it does not offer any solution for handling tensile fracture processes, which are important, for example, for predicting leads and cracks. This issue emphasises the necessity of implementing techniques for explicit fracture modelling within the DEM framework, which could substantially improve deformation predictions on long timescales. Such techniques have been developed for smaller scales, for example, for ice–structure interaction problems, based on analytical solutions (Lubbad and Løset, 2011; Lu and others, 2015). However, there is still a need to develop appropriate DEM solutions for ice fracture on regional scales.

5.4 Model limitations and future possibilities

In Section 4.4 it was found, the model drift field can be sensitive to errors in the boundary velocities, but using a sufficiently large simulation domain can substantially mitigate the propagation of these errors into more internal fields. If sufficiently accurate boundary conditions cannot be provided using, for example, CICE or other prediction tools for large-scale sea-ice drift, the DEM domain can be expanded, which may help when simulating sea ice in the MIZ. However, for pack ice, this strategy may be ineffective or impractical because a too large simulation domain may be required. Field studies indicate that the effect of coastal boundaries on the drift of pack ice is felt within ~400 km of the coast (Thorndike and Colony, 1982), and the compression of the pack against a coastal boundary can be measured by stress sensors located over 500 km offshore (Richter-Menge and others, 2002). In our simulations, we used a model domain with a size of ~100 km. Using much larger domains with the same floe sizes may be unfeasible.

The main limitation of our approach is probably due to the DEM's computational efficiency. In the presented case study, we simulated a 25 h process involving ~11 000 interacting floes with various sizes and shapes in 3-D. Without taking advantage of parallel computation and high-performance central processing units, we conducted our DEM simulations on an average laptop, which resulted in computational speeds 5–10 times slower than processing in real time, that is, simulating a real time of 1 h may take up to 10 h. Although we did not focus on computational efficiency in this feasibility study, our simulations provided insight into possible model optimisation for faster computations. We found that the DEM behaviour in our simulation scenarios was essentially 2-D, with the ice deformations primarily occurring in a horizontal plane. If out-of-plane ice deformation (e.g. rafting and ridging) can be approximated in 2-D, there may be no need to use a 3-D DEM. Instead, a 2-D non-smooth discrete model would be more suitable, potentially offering a significant improvement in computational efficiency. We leave this aspect for future development.

As our study shows, there are still challenges that need to be addressed within the DEM before it can be fully integrated with a continuum model for large-scale analyses. Therefore, we explored only the potential of one-way coupling and focused mainly on the performance of our DEM with additional forcing from boundary conditions within a large-scale context. The ability of this set-up to replicate some high-resolution patterns in sea-ice drift and deformation rates observed in satellite imagery is promising. However, to fully exploit the capabilities of DEM, for example, for improving continuum-model simulations, integrating feedback from the DEM into the continuum model may be necessary. This feedback may include updated information on the floe-size distribution, lead areas and shapes, rafted ice and other parameters that may be useful for parameterising continuum rheology models. Hopefully, future research will uncover what other data may be useful for improving these models.

6. Conclusions

The need for high-resolution sea-ice predictions, which cannot be provided by today's Earth System Models, has motivated us to investigate the feasibility of employing a nested DEM in the framework of a continuum approach commonly used for global and regional sea-ice modelling. We developed an approach where a non-smooth DEM can utilise boundary conditions obtained from a continuum sea-ice model with no feedback from the DEM, that is, so-called one-way coupling. In the presented case study, the employed DEM was capable of resolving contact interactions between thousands of ice floes that were identified in satellite images covering a region of 100 × 100 km² at 10 m resolution. This set-up was used in a regional reanalysis of sea-ice motion spanning ~25 h, which allowed evaluating the ice drift and deformation parameters on hourly and daily intervals. These parameters predicted by the DEM reasonably aligned with satellite data, although the model overestimated ice deformation.

Resolving floe–floe interactions without explicitly modelling ice fracture processes, the nested model was able to provide kilometric details of ice deformation zones, which were concentrated near landmasses, in agreement with satellite images. High-resolution modelling of ice deformation offers new capabilities needed in Earth System Models, such as predicting the formation and dynamics of leads. However, this advancement may require further development of the DEM to explicitly model ice fracture processes.

It was also found that in shallow seas, tidal currents have a strong effect on sea-ice dynamics at short temporal scales. However, ocean models may not accurately predict this current component in the upper ocean layer, where it interacts with sea ice. Moreover, very few ocean models provide current data on sub-hourly intervals, which may be needed for high-resolution DEM simulations. Therefore, the development of new sea-ice models towards high-resolution forecasting should be accompanied by upgrades to the coupled ocean models that provide forcing data for the ice component.

Overall, the results of this feasibility study suggest that combining the DEM and continuum model approaches is a valid, flexible option for future developments of high-resolution sea-ice models, which may also benefit from utilising satellite optical imagery for initialising DEM simulations. Two-way coupling mechanisms between a DEM and continuum model were not explored in this study but may offer additional improvements in large-scale sea-ice modelling.