the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Lagrangian tracking of sea ice in Community Ice CodE (CICE; version 5)
Chenhui Ning
Yan Zhang
Xuantong Wang
Zhihao Fan
Jiping Liu
Sea ice models are essential tools for simulating the thermodynamic and dynamic processes of sea ice and the coupling with the polar atmosphere and ocean. Popular models such as the Community Ice CodE (CICE) are usually based on nonmoving, locally orthogonal Eulerian grids. However, the various in situ observations, such as those from icetethered buoys and drift stations, are subjected to sea ice drift and are, hence, by nature Lagrangian. Furthermore, the statistical analysis of sea ice kinematics requires the Lagrangian perspective. As a result, the offline sea ice tracking with model output is usually carried out for many scientific and validational practices. Certain limitations exist, such as the need for highfrequency model outputs, as well as unaccountable tracking errors. In order to facilitate Lagrangian diagnostics in current sea ice models, we design and implement an online Lagrangian tracking module in CICE under the coupled model system of CESM (Community Earth System Model). In this work, we introduce its design and implementation in detail, as well as the numerical experiments for the validation and the analysis of sea ice deformation. In particular, the sea ice model is forced with historical atmospheric reanalysis data, and the Lagrangian tracking results are compared with the observed buoys' tracks for the years from 1979 to 2001. Moreover, highresolution simulations are carried out with the Lagrangian tracking to study the multiscale sea ice deformation modeled by CICE. Through scaling analysis, we show that CICE simulates multifractal sea ice deformation fields in both the spatial and the temporal domain, as well as the spatial–temporal coupling characteristics. The analysis with model output on the Eulerian grid shows systematic difference with the Lagrangiantrackingbased results, highlighting the importance of the Lagrangian perspective for scaling analysis. Related topics, including the subdaily sea ice kinematics and the potential application of the Lagrangian tracking module, are also discussed.
 Article
(10828 KB)  Fulltext XML

Supplement
(4855 KB)  BibTeX
 EndNote
Sea ice floes are inherently Lagrangian points that undergo thermodynamic and dynamic changes throughout their lifetime. Under the dynamic forcings from the atmosphere and the ocean, sea ice drifts, and, as internal stress accumulates, it deforms and undergoes plastic failures. The drift of sea ice floes is associated with constant thermodynamic growth and melt of the sea ice; hence it is fundamental to the energy and ice–water balance in the polar regions (Haas, 2009). Furthermore, highly nonlinear and anisotropic linear kinematic features manifest with the sea ice deformation fields, which prevail from meters to the geophysical scales (Marsan et al., 2004). The accurate observation and modeling of the sea ice drift and deformation are key to our scientific understanding of both the climate system and human activities in the polar region.
Due to the harsh conditions of the polar environment, the longterm direct measurements of the sea ice are usually carried out through in situ deployments of buoys. These autonomous systems, which are usually attached to the sea ice, relay information about their locations, sea ice conditions, and associated atmospheric and oceanic conditions. Since they drift with the sea ice, their locations are also representative of the sea ice floe to which they are attached. For example, the International Arctic Buoy Programme (IABP; https://iabp.apl.uw.edu, last access: 8 September 2024) compiles historical and realtime buoy measurements in the Arctic region. The data are widely used in the study of sea ice dynamics (Rigor et al., 2002) and thermodynamic processes (Perovich et al., 2014), as well as the data assimilation for numerical weather forecasts (Inoue et al., 2009).
The drift and the deformation of sea ice are the result of its dynamic response to the atmospheric and oceanic forcings. Unlike the Newtonian fluids of air and water, sea ice patches undergo multifractal deformation characterized by localized plastic faults. This highly localized and anisotropic deformation typically corresponds to sea ice leads and ridges, manifesting as linear kinematic features (LKFs; Kwok et al., 1998). While sea ice ridging is a major way that thick sea ice is formed, leads serve as important sources of heat and moisture for the polar regions, especially during winter (Rothrock, 1975; Andreas and Cash, 1999). Therefore, sea ice deformation is crucially important for the coupled climate system of the polar region. In order to study how the sea ice deforms, we usually carry out multiscale analysis through the statistics of the deformation rates (i.e., the speed at which the sea ice deforms). Specifically, the deformation rates, denoted $\dot{\mathit{\u03f5}}$'s, can be computed for individual sea ice patches at various temporal and spatial scales. Since the sea ice is constantly drifting, the computing of $\dot{\mathit{\u03f5}}$'s and the related analysis should also be carried out from a Lagrangian perspective. Satellite remotesensingbased datasets such as RGPS (Lindsay and Stern, 2003) utilize synthetic aperture radar (SAR) images collected at various times to produce largescale, highresolution maps of the sea ice drift and deformation. In particular, both the correlationbased and the feature tracking approaches when processing the SAR images ensure that the analysis of the drift/deformation is Lagrangian by nature. These datasets are widely used for the study of multiscale sea ice kinematics (Marsan et al., 2004) and the validation of numerical simulations (Kwok et al., 2008; Rampal et al., 2019).
For the simulation of sea ice and its kinematics, we construct numerical models which involve a layered structure with the following: (1) the mathematical modeling of the physical processes, (2) the numerical treatments for the spatial–temporal discretization and the integration, and (3) the code implementation and the simulation on parallel computers. Popular sea ice models, such as CICE (Community Ice CodE; https://github.com/CICEConsortium/CICE, last access: 8 September 2024) and SI^{3} (Vancoppenolle et al., 2023), are usually based on spatial discretization using locally orthogonal structured grids. Rheology models (Bouchat et al., 2022; Hutter and Losch, 2020) such as viscous–plastic (VP; see also elastic–viscous–plastic (Hunke and Lipscomb, 2008)) are capable of reproducing certain statistics of the observed multiscale sea ice kinematics (Kwok et al., 2008). However, the model's output, including the instantaneous and the average model status at different temporal scales, is typically defined on the model's native Eulerian grid. One notable exception is neXtSIM, which is based on Lagrangian moving mesh and inherently supports the scaling analysis (Rampal et al., 2019). But for CICE and many widely adopted sea ice models, the model output is insufficient, especially for the scaling analysis at large temporal scales, since it inherently requires a Lagrangian perspective. A typical practice to overcome this limitation is to reconstruct Lagrangian tracks with the model's output, such as those derived from daily velocity fields in Bouchat et al. (2022). Certain limitations are still present, however, especially given the evergrowing resolution of current models (Xu et al., 2021; Zhang et al., 2023). Highfrequency model output is needed to reconstruct realistic Lagrangian tracks, entailing substantial data storage and offline computation. Furthermore, the analysis of smallscale sea ice deformation (i.e., minute scale as in Oikkonen et al., 2017) requires even finer spatial and temporal model output and a larger overhead with the offline tracking analysis. Therefore, more flexible Lagrangian diagnostic tools are needed for the scaling analysis of sea ice kinematics and future development of sea ice models.
In this paper we introduce the online Lagrangian tracking of sea ice and its model integration in the model of CICE (version 5). The tracking of sea ice is carried out along with the model's numerical integration, and it supports very high frequency tracking (at the model's time step) and large numbers of Lagrangian points. The model integration is carried out and further validated through the numerical experiments in the coupled framework of Community Earth System Model (CESM, version 2: https://www.cesm.ucar.edu/models/cesm2, last access: 8 September 2024). Specifically, the comparison with observed buoy tracks is carried out with atmospherically forced historical simulations. Furthermore, a highresolution experiment with 7 km resolution in the Arctic region is carried out, and we evaluate the spatial–temporal scaling of wintertime sea ice deformation. In Sect. 2 we introduce the tracking algorithm and the integration in CICE in detail. Section 3 includes the numerical experiments and detailed analysis of simulation results. Finally in Sect. 4, we summarize the article and discuss related topics, including potential applications of the sea ice Lagrangian tracking and the highresolution simulation of sea ice kinematics.
The Lagrangian tracking of sea ice is tightly integrated with the dynamics processes of CICE (version 5). The model grid of CICE is a twodimensional, logically rectangular, structured grid with the size of (nx_global,ny_global) and indexed by (i,j). Typical lateral boundary conditions are supported for global configurations, including the east–west periodic boundary and the tripolar grid boundary (Murray, 1996). For the Lagrangian tracking, each active point has a specific logical location of (x,y) in the twodimensional continuous space, satisfying $\mathrm{0}\le x<$nx_global and $\mathrm{0}\le y<$ny_global. Furthermore, there is a biprojection between each point’s logical location and its corresponding geolocation. The Lagrangian tracking is carried out within the model grid's domain following the tracking algorithm while maintaining the grid's lateral boundary conditions and the land–ocean distribution.
In this section, the Lagrangian tracking algorithm is introduced in detail in Sect. 2.1, which is based on the existing advection framework of transport remapping (Dukowicz and Baumgardner, 2000). Regarding the commonly used functionality of CICE, the Lagrangian tracking also supports parallel computing, tripolar grids and a simple logging system for the record of the tracking results. Detailed implementation and the integration in CICE are covered in Sect. 2.2.
2.1 Lagrangian tracking algorithm
The Lagrangian tracking algorithm is based on the transport remapping advection scheme, which is available in CICE (Dukowicz and Baumgardner, 2000). As a conservative, twodimensional, semiLagrangian scheme, transport remapping operates on the ArakawaB staggered grid. At each dynamics step of CICE, the backward tracking of the corner points is carried out for the advection of tracers onto the Eulerian grid cells (corners of the hollow quadrilateral in Fig. 1a). For the Lagrangian tracking, we utilize the backward tracking vectors and compute the forward tracking for all four corner points (filled grey arrows in Fig. 1a), and they form the new quadrilateral (outlined in blue). Then, based on the previous location of the Lagrangian point (black circle), we can compute its new location (blue circle) as follows.
The tracking of the Lagrangian point is based on the local, normalized coordinate of the grid cell containing the point. Under the local coordinate, the cell's four corners correspond to (0.0,0.0), (0.0,1.0), (1.0,0.0) and (1.0,1.0), respectively. The previous location of the Lagrangian point under this coordinate is (x_{local},y_{local}), satisfying $\mathrm{0}\le {x}_{\mathrm{local}}<\mathrm{1}$ and $\mathrm{0}\le {y}_{\mathrm{local}}<\mathrm{1}$. Note that given the point's logical position in the model grid as (x,y), the following points hold: (1) the cell where the point is present has the index of $(\lfloor x\rfloor +\mathrm{1},\lfloor y\rfloor +\mathrm{1})$, and (2) ${x}_{\mathrm{local}}=x\lfloor x\rfloor $ and ${y}_{\mathrm{local}}=y\lfloor y\rfloor $. After the forward tracking, the four corners have the new locations of $({x}_{i,j},{y}_{i,j})$ for $i\in \mathit{\{}\mathrm{0},\mathrm{1}\mathit{\}}$ and $j\in \mathit{\{}\mathrm{0},\mathrm{1}\mathit{\}}$. Then, the new location of the Lagrangian point, denoted $({x}_{\mathrm{local}}^{*},{y}_{\mathrm{local}}^{*})$, is computed through the bilinear interpolations with x_{i,j}'s and y_{i,j}'s, assuming that its relative location within the new quadrilateral remains (x_{local},y_{local}). In the case that ${x}_{\mathrm{local}}^{*}$ or ${y}_{\mathrm{local}}^{*}$ is larger than 1 or smaller than 0, the Lagrangian point has drifted out of the current grid cell, which we denote as the migration of the point. Since CICE utilizes domain decomposition for parallel computing, the migration potentially is between blocks (subdomains after decomposition) and entails communication between parallel processes. The detailed support is introduced further in Sect. 2.2.
Based on the geolocations of the four cell corners and the Lagrangian point's relative location within the cell, we can compute the geolocation of the Lagrangian point. Specifically, first, the threedimensional locations of the cell corners are computed with their latitudes and longitudes. Through linear interpolations, we locate the threedimensional location of the Lagrangian point. Finally, the latitude and longitude of the Lagrangian point can then be determined.
2.2 Implementation in CICE
2.2.1 Software implementation
For the Lagrangian tracking in CICE, we define the data structure as lagr_point in the Fortran module of ice_transport_driver. It contains necessary fields of information for the Lagrangian point, including its current status, lifetime, current location and other essential information. Furthermore, each of the parallel CICE processes maintains a large, preallocated pool of available Lagrangian points (i.e., instances of lagr_point). When a new Lagrangian point is created in the current process or migrated from another process, an unused slot is claimed from this pool. Similarly, when the point is dead (i.e., due to melting) or migrates out of the block, the slot is reclaimed and recycled in the pool.
The life cycle of a Lagrangian point consists of several stages (Fig. 2) and transitions between the stages, called events. Upon its creation (typeI event), the Lagrangian point is assigned to a specific geolocation and, consequently, a specific block and a specific processor. The Lagrangian point drifts (typeII event), until the sea ice melts (sea ice concentration lower than 5 %, typeV event), or it is automatically deactivated (e.g., exceeding the prescribed maximum lifetime). When the Lagrangian points migrate outside the current block, they will be delivered to the corresponding adjacent block (typeIII event; see also below). This process is carried out through the builtin boundary exchange operations of CICE. A special landing event of Lagrangian points (typeIV) is supported but not possible in the current implementation.
2.2.2 Overall model integration
As shown in Fig. 3, the tracking of all Lagrangian points occurs within each dynamics time step (i.e., step_dynamics) after the computation of the prognostic velocities and the advection process (which computes the backward tracking vectors). To ensure numerically stable integration for the advection, there is strict limitation over the time step and hence the backward/forward tracking vectors. As a consequence, Lagrangian points cannot migrate beyond one grid cell in either direction. When a specific point migrates beyond the current block's boundary (i.e., thick rectangular boundary in Fig. 2), it will be strictly within the outer boundary of the block (i.e., within the thin rectangular boundary in Fig. 2).
The Lagrangian tracking contains four major steps (shown below). For each Lagrangian point, its status is maintained, it is tracked with the forward drifting vectors and its location information is updated. In the case of the Lagrangian point migrating outside the current block, it is recorded for further boundary exchange. Then the Lagrangian points are exchanged between blocks, with newly migrated Lagrangian points recorded for the current block.

For each active Lagrangian point, the following steps should be followed:

Check for deactivation (due to melt, lifespan, etc.), with necessary management of the Lagrangian point slots.

Increase the ge of the Lagrangian point.

Retrieve the local advection information.

Do Lagrangian tracking, and update the logical and the geophysical positions.

Check for potential migration out of the current block.


For a migrating point, the following steps should be followed:

Record its information for boundary exchange.

Carry out management of the Lagrangian point slot.

Carry out boundary exchanges for Lagrangian points.

Activate newly migrated Lagrangian points.

Deactivate Lagrangian points (due to landing, etc.) with necessary management steps.

2.2.3 Miscellaneous
Boundary exchanges for migrating Lagrangian points
Each dynamics time step involves a single Fortran function call to the boundary exchange, carrying out the migration of Lagrangian points between blocks. The maximum exchangeable Lagrangian point count per boundary cell can be changed in the code through the Fortran parameter LAGR_BNDY_SIZE_PARAM. It is worth noting that the single boundary exchange due to the migration of Lagrangian points only incurs a very small computational overhead. Within the same dynamics time step, the elastic–viscous–plastic (EVP) solver usually contains over 100 subcycles of the elastic waves and the corresponding boundary exchanges (Lemieux et al., 2010; Xu et al., 2021). Furthermore, the Lagrangian tracking utilizes the transport remapping scheme's backward tracking. In general, the EVP solver dominates the overall simulation time of the sea ice dynamics, with the Lagrangian tracking consuming less than 3 % of the time for all the numerical experiments in Sect. 3.
Support for tripolar grids
Tripolar grids are commonly used in highresolution global configurations of CICE and CESM (Small et al., 2014). For example, the 0.1° grid of TX0.1 is used for the global ocean mesoscaleresolving simulations of CESM. The commonly used Ufold tripolar grids are supported for our implementation of the Lagrangian tracking in CICE, with the schematics shown in Fig. 4.
When the Lagrangian point drifts beyond the northern boundary of the blocks on the tripolar folding line, it is passed to the corresponding block and grid cell through the boundary exchange process. This corresponds to the case with y larger than ny_global, where ny_global is the grid cell count in the j direction. Then, the x and y positions of the point are modified accordingly to the correct values, following the topology on the tripolar boundary.
Logging system
A simple textbased logging system is implemented to report the results of the Lagrangian tracking. Each parallel processor generates a log file that contains the records of all active Lagrangian points on the processor. The records include the major events, including the creation of the Lagrangian points, the melt events and the migration events, as well as the Lagrangian points' status at regular time intervals. The time step count for reporting Lagrangian points' information is a userprescribed compiletime parameter (see Appendix C for details).
We carry out numerical experiments of the Lagrangian tracking with CICE (version 5) and CESM (version 2). CICE is configured with five ice thickness categories with multiple vertical layers, as well as full thermodynamic and dynamic processes. Key model configurations include the mushylayer vertical thermodynamics, the elastic–viscous–plastic (EVP) rheology model and the deltaEddington radiation scheme. Detailed parameterization schemes and major parameters of CICE are further covered in Appendix A. Two model resolutions of CICE are used: the nominal 1° grid of GX1V6, which is built in in CESM, and the nominal 0.15° grid of TS015 (previously implemented in CESM (Xu et al., 2021)). Notably, the TS015 grid has a horizontal dimension of 2400×1680 globally, and the mean grid resolution is about 7 km in the Arctic region. Prominent, multifractal sea ice deformation can be simulated at this resolution (Xu et al., 2021).
For all the experiments, CICE is coupled to the slab ocean model (SOM) and forced by the COREII dataset through the coupling framework in CESM. The COREII dataset is based on NCEP atmospheric reanalysis and further used in the Ocean Model Intercomparison Projects (Griffies et al., 2016). It contains two separate forcing datasets: the normal year forcing (NYF), which is annually repeating, and the interannual forcing (IAF), which is for the years between 1948 and 2007. While the NYF dataset is usually used for the longterm simulations and the spinup of the ocean–sea ice coupled model, the IAF dataset can be used for the hindcast of the historical states of the ocean and the sea ice (Wang et al., 2016). Both the NYF and the IAF datasets are used for the numerical experiments, covered in the rest of this section.
The simulated Arctic sea ice climatology under NYF is shown in Fig. 5. The seasonality of the sea ice extent at both 1 and 0.15° (not shown) is consistent with the observational climatology (NSIDC). The overestimation of sea ice extent (SIE) during winter months mainly manifests in the peripheral seas of lower latitudes. In particular, in the Atlantic Arctic region, the overestimation of SIE may be due to the absence of ocean heat transport of the SOM. During summer, SIE agrees well with the observation. In terms of the sea ice volume, there is general underestimation compared with PIOMAS (Schweiger et al., 2011). Thick, multiyear sea ice manifests north to Greenland and the Canadian Arctic Archipelago (CAA), with the clockwise circulation in the Canadian Basin as controlled by the Beaufort High. With the generally good agreement of the modeled sea ice states with the observational datasets, we consider it sufficient for further analysis of the Lagrangian tracking. Moreover, further improvements of the model's simulation results are planned for future work, including the coupling with a dynamic ocean model and the tuning of parameterization schemes in CICE.
3.1 Climatology of Lagrangian points under NYF
Under the annually repeating NYF dataset, we carry out the basic test of the Lagrangian tracking with the lowresolution, 1° grid. After the model has reached quasiequilibrium after a spinup of 20 years, we deploy one Lagrangian point at the center of every grid cell with sea ice (sea ice concentration (SIC) >5 %) in the Northern Hemisphere. The initial locations of the points are shown in Fig. 6a. With the summer melt and the transport to lower latitudes, the count of active Lagrangian points decreases with each year. After a full year since the deployment, the Lagrangian points outside the Arctic Basin are all lost due to melting (Fig. 6b). In Fig. 7 we show the 2month tracks of the points for the first whole winter after deployment. The overall distribution of remaining points after a melt season, as well as their tracks in the winter, shows (1) the counterclockwise drift in the Beaufort Gyre at distinctive stages throughout the winter, (2) the convergence of points to the north of Greenland/the CAA, and (3) the transpolar drift from the eastern part of the basin and the outflow in the Fram Strait. After 5 years and 10 years (Fig. 6c and d), the Lagrangian points are generally lost. The surviving points are mainly retained due to the accumulation in the Beaufort Gyre. Among all of the 3599 points that are originally in the basin (80° N and further north), about 2063 (57.3 %) points are lost through melting in the basin within the first 20 years, while 1083 (30.1 %) are lost due to export from the Fram Strait (FS). In later years, the surviving points are mainly lost due to the circulation to the transpolar drift (hence the FS export) or to the eastern part of the basin with a more prominent summer melt. The daily locations of the Lagrangian points are shown by the animation in the Supplement. It is worth noting that the NYF forcing dataset does not support direct comparison with observations. Hence the experiment results, including the simulated circulation of the Lagrangian points and the loss due to melt process, provide us with a basic test of the Lagrangian tracking functionalities.
3.2 Validation of Lagrangian tracking under IAF
We further evaluate the Lagrangian tracking through the IAFbased simulation and compare it with the buoy measurements. Based on the 1° grid (GX1V6) and the model's spinup status under NYF, we carry out the historical simulation for the years between 1979 and 2001. For every day during the winter (December to March) at 00:00 (UTC), we deploy a Lagrangian point in every grid cell where sea ice is present. For each Lagrangian point, we track its location for up to 90 d. The daily instantaneous locations of the points are recorded and further compared with the observations.
The observational dataset from the International Arctic Buoy Programme (IABP; Rigor et al., 2008) during the same years is used for validation (data available at https://iabp.apl.uw.edu/data.html, last access: 10 May 2023). Since we deploy Lagrangian points to the center of the grid cells, we use the following criteria to match the physical buoys' locations of IABP. In total, 621 buoys with 49 004 hourly locations are available for comparison. We further split each buoy's continuous track into subtracks of 14 d for further comparisons. For each buoy's (sub)track, we screen over the newly deployed Lagrangian points on its starting date. The Lagrangian point nearest to the track's starting location is located, and the simulated track is then compared against the observation. Leap years are ignored in the IAF dataset and the simulations; hence the buoys' tracks starting on 29 February are matched with the Lagrangian points on 1 March.
Figure 8 shows the matching between the buoys' and the corresponding Lagrangian points' locations during the winter of 1993–1994 (December to March). All of the subtracks, each covering 14 d of tracking of the buoy, are shown. The two major sea ice drift regimes include the clockwise circulation in the Canadian Basin and the clockwise drift and export of the sea ice in the region east to 60° W and west to 130° E. In both regimes, the modeled Lagrangian points' tracks are highly consistent with the observations.
Furthermore, Fig. 9 shows the distance between the locations of the buoys and the corresponding Lagrangian points for all the subtracks. In total, there are 2849 subtracks for comparison, each corresponding to a buoy's drift within 14 d. For the starting locations of the buoys' subtracks, the average distances to the matching virtual buoys are all within 50 km. The average grid size within the Arctic Basin is about 50 km, with finer (coarser) resolution in the area near (far from) Greenland. For over 75 % of the matching buoy pairs, the initial distance is less than 27 km. After the tracking for 14 d, the distance gradually grows, but for 75 % of all pairs, the distance is within 60 km. Besides, the median (mean) distance is 38 km (43 km). On average, the drift distance of buoys is 87 km at 14 d scale. For fastmoving buoys such as those entering the Fram Strait, the distances between the Lagrangian points and the matching physical buoys are larger, but the relative error always remains low within the 2week tracking period. In general, based on the Lagrangian tracking in CICE, the sea ice drift as simulated by the model matches buoys' observations well. The tracking uncertainties may arise from the limited spatial resolution of the model, the uncertainty of the atmospheric forcing dataset and the sea ice model's dynamics in simulating the observed drifts. Besides, due to the regular deployment of the Lagrangian points, there is no exact match of the buoys' initial locations. Further attribution of the tracking error to the various contributing factors, including the model and the data's uncertainty, as well as the initial location errors, is planned for future study.
3.3 Scaling analysis of sea ice kinematics in the 7 km resolution simulation
We further carry out highresolution simulations of sea ice kinematics and related analysis of the sea ice deformation. In particular, with the 7 km resolution of the TS015 grid in the Arctic region, the model is capable of simulating the multifractal sea ice deformation (Xu et al., 2021). The experiment is based on the NYF dataset and the model's quasiequilibrium status after the spinup process. The Lagrangian tracking is utilized for the diagnostics of the wintertime sea ice deformation and related statistics over multiple spatial and temporal scales. We first introduce the basic framework of scaling analysis and the use of Lagrangian tracking in Sect. 3.3.1. Furthermore, the detailed analysis of the spatial scaling (Sect. 3.3.2) and the temporal scaling and the spatial–temporal coupling (Sect. 3.3.3) are covered.
3.3.1 Lagrangian tracking and scaling analysis
The sea ice deformation rates ($\dot{\mathit{\u03f5}}$'s) are defined with simulated velocities. Since CICE is based on the locally orthogonal grid and ArakawaB staggering, the model's velocity as well as the derived velocity of Lagrangian points are defined as fields of u's and v's. The gridinvariant deformation rates, including the divergence rate (${\dot{\mathit{\u03f5}}}_{\mathrm{div}}$) and the shear rate (${\dot{\mathit{\u03f5}}}_{\mathrm{shear}}$), are then defined as follows. The total deformation rate (${\dot{\mathit{\u03f5}}}_{\mathrm{total}}$), which is subjected to further scaling analysis, can be further defined with ${\dot{\mathit{\u03f5}}}_{\mathrm{div}}$ and ${\dot{\mathit{\u03f5}}}_{\mathrm{shear}}$.
In these equations, u_{x}, u_{y}, v_{x} and v_{y} are the spatial derivatives of u and v in the two orthogonal directions. The direct model output of the instantaneous and/or the temporal mean of the velocity fields, as well as the various deformation rates, is defined on the Eulerian model grid. In Fig. 10 we show the maps of the dailymean ${\dot{\mathit{\u03f5}}}_{\mathrm{total}}$ in the Arctic region for 2 typical days of the experiment with TS015. LKFs manifest with the highly localized, anisotropic sea ice deformation. It is worth noting that the dailymean velocity and the deformation fields (shown in Fig. 10) are inherently Eulerian and hence not appropriate for the scaling analysis. For example, the offline tracking (at hourly scale) should be carried out for further analysis (e.g., Bouchat et al., 2022).
The velocity derivatives at various temporal–spatial scales (i.e., $\mathit{\{}u,v{\mathit{\}}}_{\mathit{\{}x,y\mathit{\}}}$'s) are further computed over the various Lagrangian patches. A set of adjacent Lagrangian points form an enclosed sea ice patch and start with the position on regular Eulerian grid points. Their locations change with the Lagrangian tracking, due to the sea ice drift and deformation. At a certain time (delayed by T), their new locations are recorded and the displacement from their original locations are computed as Δ_{x}'s and Δ_{y}'s. Then the average velocity $\stackrel{\mathrm{\u203e}}{u}$ and $\stackrel{\mathrm{\u203e}}{v}$ for each Lagrangian point can be computed as $\stackrel{\mathrm{\u203e}}{u}=\frac{{\mathrm{\Delta}}_{x}}{{\mathrm{\Delta}}_{t}}$ and $\stackrel{\mathrm{\u203e}}{v}=\frac{{\mathrm{\Delta}}_{y}}{{\mathrm{\Delta}}_{t}}$, respectively. The spatial velocity derivatives of the Lagrangian patch are computed with the patch's area (A) and the line integral of velocity over its outer boundary, following Kwok et al. (2008) and Rampal et al. (2019). Details of the computation are further covered in Appendix B.
Given the total deformation rates of all sea ice patches, we compute the average value of ${\dot{\mathit{\u03f5}}}_{\mathrm{total}}$'s within the similar spatial scale ($L=\sqrt{A}$). In particular, the qth order of ${\dot{\mathit{\u03f5}}}_{\mathrm{total}}$'s, computed as ${\dot{\mathit{\u03f5}}}_{\mathrm{total}}^{q}$'s, along with their average values, is also computed for $q\in \mathit{\{}\mathrm{0.5},\mathrm{1},\mathrm{1.5},\mathrm{2},\mathrm{2.5},\mathrm{3}\mathit{\}}$. By the scaling law (Marsan et al., 2004), under the spatial scale of L and the temporal scale of T, we have
where α(q) and β(q) are the structure functions of the temporal and the spatial scaling. For a specific value of q, we carry out the linear fitting of $\u2329{\dot{\mathit{\u03f5}}}^{q}\left(T,L\right)\u232a$ with respect to T (or L) to estimate the value of α (or β). For the observed multifractal deformation of the sea ice (Marsan et al., 2004), we witness convex structure functions for both α(q) and β(q) with respect to q. A generalized analysis framework with a nonfixed degree of multifractality for the sea ice formation is also available (Weiss, 2008; Bouchat et al., 2022). For comparison, the forms in Eqs. (6) and (7) also assume the underlying multifractal, lognormal multiplicative model. In this study, they are adopted because their quadratic form is sufficient for capturing the convex shape of the structure functions (Marsan et al., 2004; Rampal et al., 2019). Correspondingly, α(q) and β(q) can be fitted as follows, with the fitted parameters of a and c both larger than 0.
In order to evaluate the simulated sea ice deformation, we deploy a Lagrangian point at the center of each grid cell for every model day at the time of 00:00. The maximum allowed lifespan of each Lagrangian point is 30 d. As a consequence, the temporal scaling of up to 30 d can be studied (i.e., the maximum value of T at 30 d). The model regularly reports the location of every Lagrangian point every 6 h (i.e., the minimum value of T at 6 h). Both parameters can be configured at the compile time of the CICE model. The model grid's native resolution in the Arctic is about 7.3 km. Given that the effective resolution is coarser (Xu et al., 2021), we evaluate the deformation at the spatial scale from 30 to 480 km (4 times to 64 times the grid's resolution). For the scaling analysis, we limit the initial locations of all Lagrangian points to be at least 400 km away from land (i.e., within the outlined region in Fig. 10). The regions near the coast and the continental shelf break, where constant deformation features persist, are therefore excluded from the scaling analysis.
After the model's spinup, we output a whole winter's Lagrangian tracking results, together with the daily mean Eulerian fields of the sea ice velocity and the deformation. Furthermore, we also study the difference between the Lagrangiantrackingbased scaling analysis and that based on the temporal–spatial mean fields of the model output at Eulerian locations.
3.3.2 Spatial scaling
We carry out spatial scaling analysis for two typical periods during the winter, focusing on four typical temporal scales: 1, 3, 10 and 30 d. Figures 11 and 12 show the results around 20 December and 6 February, respectively. At q=1 and 1 d scale, the spatial scalings of sea ice deformation are both very shallow: β=0.0623 and 0.121, respectively. This period of 2 d corresponds to different strengths and depths of the Beaufort High and slightly different Arctic Oscillation indices, and consequently the drift and deformation fields are distinct (see Figs. 6 and 11 of Xu et al., 2021). At different orders of momentum, the mean deformation rates all follow the powerlaw scaling. On higher orders (q>1), large deformation events are more dominant in the mean deformation rate, corresponding to steeper scaling (i.e., larger values of β). Furthermore, the value of β grows nonlinearly with q. With the increase in the moment order, the deformation rate decreases faster with respect to the increase in the spatial scale (L). Correspondingly, the convex shape of β(q) indicates that the model simulates multifractal sea ice deformation.
For the period of study surrounding 20 December, at larger temporal scales from 1 to 30 d, the average deformation rates decrease at all spatial scales and moments of order. More importantly, the value and the shape of β(q) also change significantly with the temporal scale. At q=1 and the 1 d temporal scale, the value of β is 0.0623. For larger temporal scales, the mean deformation rate gradually decreases, which also applies to the structure function of β. In particular, the convexity of β(q), indicated by the fitted value of a, decreases drastically from 0.137 (1 d) to 0.064 (30 d). The change in a indicates that there is a strong coupling between the spatial and temporal domain for the sea ice deformation. Similar results hold for the period around 6 February, with the two following major differences: (1) a more convex β(q) around 6 February than 20 December and (2) a faster decrease in a with respect to T (i.e., the temporal scale). The most significant drop in a occurred from the 3 to the 10 d scale for 6 February (from 0.151 to 0.101). For comparison, that for 20 December is between 10 and 30 d (from 0.093 to 0.064). Regarding the differences between the two periods, we conjecture that the processdependent deformations are the major cause, including the strength of the spatial–temporal coupling. Further analysis is needed for the attribution of these differences to various factors, including the sea ice status, its deformation history and the forcings.
For comparison, we also show the results based on Eulerian deformation fields in both Figs. 11 and 12, including the powerlaw fittings and the structure functions of β(q). The method details are introduced in Appendix B. For both periods of study, the scaling analysis using Lagrangian tracking results shows steeper β functions with higher convexity across all temporal scales (i.e., larger values of a). For 6 February, the difference between the two is more pronounced. Note that the scaling analysis should be carried out based on the Lagrangian perspective. The objective of the comparisons is to demonstrate that there are systematic differences in the scaling analysis when using the Eulerian model outputs. The quantitative differences in the deformation rates may arise from (1) the deformation events being misaligned between the Eulerian and the Lagrangian perspectives and (2) the secondary effect of changing shapes and scales when the sea ice deforms within the Lagrangian perspective, which is not captured by the Eulerian framework.
3.3.3 Temporal scaling and spatial–temporal coupling
For the temporal scaling analysis, we focus on the spatial scale of the model's effective resolution (L=22 km). Figure 13 shows the results for the outlined region in Fig. 10 in the month of December (left) and February (right). The analyzed temporal scale is in the range between 1 and 30 d. Similar to the spatial scaling results, the powerlaw scaling is witnessed for the temporal scaling (top panels of Fig. 13). Also, convex structure functions of α(q) are present for both months, with the fitted value of c as 0.063 and 0.059, respectively. These results indicate that the model also simulates multifractal sea ice deformation in the temporal domain.
Evidently, the powerlaw scaling at subdaily scales is much shallower than that between 1 and 30 d. Similar behavior of the sea ice is observed in Oikkonen et al. (2017), which covers the temporal scale from about 10 min. We consider the shallower scaling in the simulations to be qualitatively consistent with the observation. However, although the model's dynamics time step is sufficiently short (see Appendix A), the temporal resolution of the forcing dataset is much coarser, at 6hourly. In order to fully study the simulation of the subdaily sea ice deformation, we need highfrequency forcing datasets or the coupled simulation with the highresolution interactive atmospheric component (Zhang et al., 2023).
Similar to the analysis in Sect. 3.3.2, we also compute the temporalmean Eulerian deformation fields and the equivalent scaling analysis. Since we only output daily model fields, the analysis is limited to the temporal scale between 1 and 30 d (middle row of Fig. 13). Apparent powerlaw scaling is witnessed for both months. However, the structure function of α(q) with the Lagrangian tracking is systematically less convex than that based on Eulerian fields, with lower values of c for both months (bottom row of Fig. 13). This result, together with that in Sect. 3.3.2, shows that the scaling analysis should be carried out in the Lagrangian framework. Although the scaling analysis with the model's Eulerian mean fields yields apparent multifractal sea ice deformation, the results are systematically different from those based on the Lagrangian analysis.
Furthermore, we conduct the initial analysis of the spatial–temporal coupling during the whole winter (December to March). For the region of study (outlined in Fig. 10), we utilize all of the Lagrangian tracking results by forming and tracking Lagrangian patches at various spatial and temporal scales. Specifically, for each combination of the spatial and the temporal scale (i.e., L and T), we form Lagrangian patches that satisfy the following criteria: (1) they start at interleaved Eulerian grid locations separated by $L/\mathrm{2}$ in both directions, and (2) the time difference between their starting time and 0:00 of 1 December is separated by $n\cdot T/\mathrm{2}$, where $n\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots}\mathit{\}}$. With these two criteria, we in effect split the Lagrangian tracking results into sea ice Lagrangian patches with 50 % temporal and spatial overlapping. As a result, not only do we attain full coverage of the study area/period, but also we avoid potential sampling issues associated with the changing weather regimes.
The relationship between the curvature parameter of the spatial (temporal) scaling structure function and the temporal (spatial) scale is shown in Fig. 14a (b). There is evident coupling between the spatial domain and the temporal domain, with decreased curvature of the spatial (temporal) structure function at larger temporal (spatial) scales. In particular, for the temporal scaling, there is a good fit of the curvature parameter c to the power law (Fig. 14b), which is consistent with various estimations based on in situ and satellitebased remote sensing observations (Rampal et al., 2008; Marsan and Weiss, 2010). However, for the spatial scaling, the relationship of a to the temporal scale is much flatter and, in particular, less convex for a powerlaw fit (Fig. 14a). We note that during different periods of the winter, the sea ice drift and deformation patterns are highly heterogeneous (due to changing weather), with the sea ice conditions undergoing significant changes. For example, the spatial scaling exponent shows large temporal variability (Rampal et al., 2019) and is likely linked to the atmospheric forcing patterns (Xu et al., 2021). Therefore, we adopt another fitting with the logquadratic form: $y=p\cdot \mathrm{log}(x{)}^{\mathrm{2}}+q$. This new form yields a much higher fitting to the relationship between the curvature parameter a and the temporal scale T (R^{2} from 0.882 to 0.989). The efficacy of the powerlaw fitting, as well as the full analysis of the spatial–temporal coupling as simulated by CICE, is beyond the scope of this study. In particular, the sensitivity to the sea ice rheology model and other dynamic processes is planned for both CICE on its own and coupled experiments with CESM.
We design and implement the online Lagrangian tracking in the sea ice model of CICE. It is incorporated in the dynamics process of CICE and fully supports the domain decomposition and parallel computing in CICE. Compared with the common practice of offline sea ice tracking (Bouchat et al., 2022; Sumata et al., 2023), the online tracking has several advantages. First, the offline tracking usually requires the storage of highfrequency sea ice drift data, which is not needed for the online tracking. Second, since the online tracking is carried out per dynamics step, the tracking error is minimal. For offline tracking, the tracking uncertainty arises from the relatively coarse temporalmean (e.g., daily) sea ice drift fields, as well as the variable tracking frequency (see examples in Bouchat et al., 2022). Moreover, the study of sea ice dynamics at small temporal scales (such as subdaily motions) with offline tracking requires an even finer resolution of the stored sea ice drift. The study of these processes is readily supported by the online tracking (see Fig. 13 and also Oikkonen et al., 2017). Computationally, the Lagrangian tracking only induces a small overhead to the overall simulation. For example, for the daily deployment of Lagrangian points in the 7 km experiment in Sect. 3, there are on average 30 to 35 million active Lagrangian points, and the extra computational overhead is less than 3 %. Therefore, the current implementation is capable of supporting highresolution simulations with a large quantity of Lagrangian points.
In Sect. 3.3 we carried out the spatial and the temporal scaling analysis with the Lagrangian tracking based on highresolution simulation during winter. As shown, the CICE model simulates multifractal sea ice deformation fields both spatially and temporally; it also simulates the tight coupling between the spatial and the temporal domain. The scaling properties of the simulations are consistent with observed statistics based on highresolution synthetic aperture radar satellite payloads. In particular, we compare the Lagrangianbased scaling statistics with the counterparts based on Eulerian model outputs. Results highlight the importance of using Lagrangianbased diagnosis for the scaling analysis: although the analysis based on Eulerian output also yields convex structure functions, the fitted convexity parameters are significantly different from those based on the Lagrangian perspective.
In order to compare with the observed sea ice drift and deformation (e.g., RGPS) as well as the scaling statistics, we plan to carry out highresolution, atmospherically forced historical simulations with the CICE model or CESM assisted with the Lagrangian tracking. Consequently, fine tuning of the sea ice model parameters is needed for (1) the optimization of the simulation of various sea ice parameters such as sea ice thickness and (2) improved comparability to the observed sea ice deformation. In particular, in Bouchat and Tremblay (2020), the authors proposed a new estimation of deformation rates regarding their uncertainties. Better consistency with model simulations is attained with further incorporation of the model simulation and tracking uncertainties (Bouchat et al., 2022). We plan to incorporate the method and carry out a systematic evaluation of the model's capability to simulate the observed multiscale sea ice deformation. In particular, the sensitivity to the sea ice rheology is a key aspect of the planned future work (Tsamados et al., 2013; Bouchat et al., 2022).
In this first version of the Lagrangian tracking in CICE, virtual buoys are supported, which can be deployed at regular time intervals and in regular locations. Section 3.2 shows that with the pairing between the physical buoys and the respective virtual buoys, the simulated tracks are generally consistent with observations. For future development, we plan to support the deployment of physical buoys in the simulation at prescribed times and locations so that we can make better comparisons with observational datasets such as IABP. Moreover, key Lagrangiantrackingrelated parameters (in Appendix B) are planned to be configured through a subset of the namelist in the future development plan. Finally, the current design of Lagrangian tracking is based on the ArakawaB grid sea ice dynamics and the transport remapping scheme. For the compatibility with ongoing development in CICE6 (https://github.com/CICEConsortium/CICE/releases, last access: 8 September 2024), we also plan to incorporate fluxformbased Lagrangian tracking, which is suitable for ArakawaC grid sea ice dynamics (Lemieux et al., 2024).
The Lagrangian tracking can be used for a variety of science questions and applications related to sea ice. Largescale Lagrangian surveys can be carried out in the simulations, complementary to buoys' measurements, which are limited in terms of spatial and temporal coverage. For example, virtual buoys can be deployed in historical simulations, which enables the systematic study of the thermodynamic and dynamic history of the sea ice (Sumata et al., 2023). Furthermore, CICE is widely adopted for the highresolution sea ice operational systems (Smith et al., 2016; Yang et al., 2020). The online tracking of sea ice in key regions can be carried out for sea ice forecasts, in order to support operations and ship navigation in polar waters.
CICE (version 5.0) is configured with five discretized thickness categories for the ice thickness distribution (ITD) and eight and three layers for sea ice and its snow cover, respectively. The elastic–viscous–plastic (EVP) rheology model is used for all experiments. The ice strength parameterization follows that in Rothrock (1975), which relates the ice strength to the gain of the potential energy during the ridging process. The ridging–rafting parameterization uses an exponential distribution form of the ridged ice in the ITD. For sea ice thermodynamics, the mushylayer physics parameterization is adopted (Turner et al., 2013). The deltaEddington (DE) scheme is used for the radiation processes (Holland et al., 2012). Table A1 lists the key parameters of the model. The parameters are the same across all numerical experiments by default, except the dynamics time step and the EVP subcycling count (NDTE, number of time steps (DT) for the elastic waves). These two parameters are set differently between the experiments with GX1V6 and those with TS015. The shorter time step is adopted for TS015 to ensure numerical stability, and the value of NDTE is further enlarged for the numerical convergence of the EVP rheology (Lemieux et al., 2010; Xu et al., 2021).
We carry out the scaling analysis of the sea ice deformation at various spatial and temporal scales (Rampal et al., 2019; Bouchat and Tremblay, 2020). The Lagrangian tracking results, as well as the model outputs on Eulerian grids, are used to compute the mean sea ice deformation rates at the corresponding scales. A set of Lagrangian points form an enclosed Lagrangian patch, for which we compute the spatial derivatives of velocities (i.e., u_{x}, u_{y}, v_{x} and v_{y}) as follows:
These derivatives are used to compute the deformation rates as in Eqs. (2) to (3). The line integral in Eqs. (B1) is carried out over the Lagrangian patch that originally covers a regular rectangular domain. Then each point on the patch's outer boundary, called a vertex, has its location of $({x}_{i}^{j},{y}_{i}^{j})$ at the time of t_{j} and the new location of $({x}_{i}^{j+\mathrm{1}},{y}_{i}^{j+\mathrm{1}})$ at the time of t_{j+1}. Figure B1 shows the schematics for a Lagrangian patch with four corner points. We can estimate its mean velocity (u_{i},v_{i}) as the ratios of the displacements to the time difference: $\left({x}_{i}^{j+\mathrm{1}}\mathrm{}{x}_{i}^{j}\right)/({t}_{j+\mathrm{1}}{t}_{j})$ and $\left({y}_{i}^{j+\mathrm{1}}\mathrm{}{y}_{i}^{j}\right)/({t}_{j+\mathrm{1}}{t}_{j})$. Furthermore, in order to compute the line integrals at time t, we compute the mean velocity on each edge as follows (for the example in Fig. B1):
where A is the area the set of Lagrangian points covers:
In order to complete the line integral, we set x_{5}=x_{1} and, similarly, the corresponding values of y_{5}, u_{5} and v_{5}. For larger spatial scales, we start with a larger set of Lagrangian points that have their original locations on the Eulerian grid locations and cover a locally rectangular area. Then, the velocities, their spatial derivatives and the line integrals are computed, followed by the computation of the deformation rates. This is equivalent to the aggregation of all the smallest rectangular units of four Lagrangian points that constitute the whole rectangular area. Note that at the beginning of each model day, we deploy the Lagrangian points on all Eulerian grid locations, so that we attain full spatial and temporal coverage, from the daily scale and above. Additionally, we further exclude the cases with the size change by a factor of 2 or larger.
For comparison, we also compute the equivalent deformation rates based on Eulerian fields. We directly obtain the temporally mean velocities on each Eulerian grid location and compute the velocity gradients correspondingly (Fig. B1b). Different spatial scales correspond to the various areas that are associated with different rectangular grid cells, which do not change throughout the simulation.
For both Lagrangian statistics and the Eulerian counterparts, we only compute the deformation field that is contained by the areas covered at the coarsest spatial scale. This practice ensures that there is no preferential sampling of the deformation events.
A wide range of both spatial and temporal scales are adopted for the computation of the deformation rates. The statistics of $\dot{\mathit{\u03f5}}(T,L)$, which is the deformation rate at the spatial scale of L and the temporal scale of T, include its qth order (i.e., ${\dot{\mathit{\u03f5}}}^{q}(T,L)$), and the mean values are computed over the studied area. It is worth noting that the spatial scale of L varies within the studied domain, as well as with time for the Lagrangian tracking. Therefore, both the mean value of ${\dot{\mathit{\u03f5}}}^{q}(T,L)$'s and that of the spatial scales of L's are used for the scaling analysis.
Table C1 shows the major userspecified parameters for running the Lagrangian tracking in CICE. In the current implementation, they are configured at the compile time of the model. In the future, we plan to implement certain parameters to be configurable through the namelist.
The original codebase of CESM (version 2, https://www.cesm.ucar.edu/models/cesm2, UCAR, 2024) and the associated CICE model (version 5) are available through https://github.com/ESCOMP/CESM/tree/releasecesm2.2.2 (last access: 20 January 2024). The CESM code is also archived at https://doi.org/10.5281/zenodo.12616647 (Ning and Xu, 2024). The sea ice concentration data are available at NSIDC through https://noaadata.apps.nsidc.org/NOAA/G02135/ (NSIDC, 2024). The sea ice thickness dataset of PIOMAS was downloaded from https://psc.apl.uw.edu/research/projects/arcticseaicevolumeanomaly/data/ (Polar Science Center, 2024). The buoy tracks from the IABP are downloaded from https://iabp.apl.uw.edu/Data_Products/BUOY_DATA/3HOURLY_DATA/ (International Arctic Buoy Programme, 2024).
The Lagrangian tracking in CICE (version 5) and the sample Lagrangian tracking output are available at https://doi.org/10.5281/zenodo.12616647 (Ning and Xu, 2024). The CICE namelist used in the numerical simulations is also included.
The animation produced with the Lagrangian tracking with the 1° grid under NYF dataset is accessible in the Supplement of this article.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1768472024supplement.
JL and SX conceived the work. SX carried out the design and the implementation of the Lagrangian tracking in CICE. SX, CN and YZ carried out the numerical experiments. CN and SX analyzed the results. All the authors contributed to the writing of the manuscript.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
The authors would like to thank the editors and referees for their invaluable efforts in improving the manuscript. This work is mainly supported by the joint project of INTERAAC cofunded by the National Key R&D Program of China (grant no. 2022YFE0106700) and the Research Council of Norway (grant no. 328957). Jiping Liu is supported by the National Key R&D Program of China (grant no. 2018YFA0605900). Shiming Xu is also partially supported by the National Natural Science Foundation of China (grant no. 42030602), the International Partnership Program of Chinese Academy of Sciences (grant no. 183311KYSB20200015) and the Research Council of Norway (grant no. 328957). The authors would also like to acknowledge the computational and technical support from the National Supercomputing Center at Wuxi for the numerical experiments.
This research has been supported by the National Key Research and Development Program of China (grant no. 2022YFE0106700), the Norges forskningsråd (grant no. 328957), the National Natural Science Foundation of China (grant no. 42030602), the National Key Research and Development Program of China (grant no. 2018YFA0605900), and the International Partnership Program of Chinese Academy of Sciences (grant no. 183311KYSB20200015).
This paper was edited by Christopher Horvat and reviewed by two anonymous referees.
Andreas, E. L. and Cash, B. A.: Convective heat transfer over wintertime leads and polynyas, J. Geophys. Res.Oceans, 104, 25721–25734, https://doi.org/10.1029/1999JC900241, 1999. a
Bouchat, A. and Tremblay, B.: Reassessing the Quality of SeaIce Deformation Estimates Derived From the RADARSAT Geophysical Processor System and Its Impact on the Spatiotemporal Scaling Statistics, J. Geophys. Res.Oceans, 125, e2019JC015944, https://doi.org/10.1029/2019JC015944, 2020. a, b
Bouchat, A., Hutter, N., Chanut, J., Dupont, F., Dukhovskoy, D., Garric, G., Lee, Y. J., Lemieux, J.F., Lique, C., Losch, M., Maslowski, W., Myers, P. G., Olason, E., Rampal, P., Rasmussen, T., Talandier, C., Tremblay, B., and Wang, Q.: Sea Ice Rheology Experiment (SIREx): 1. Scaling and Statistical Properties of SeaIce Deformation Fields, J. Geophys. Res.Oceans, 127, e2021JC017667, https://doi.org/10.1029/2021JC017667, 2022. a, b, c, d, e, f, g, h
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
Griffies, S. M., Danabasoglu, G., Durack, P. J., Adcroft, A. J., Balaji, V., Böning, C. W., Chassignet, E. P., Curchitser, E., Deshayes, J., Drange, H., FoxKemper, B., Gleckler, P. J., Gregory, J. M., Haak, H., Hallberg, R. W., Heimbach, P., Hewitt, H. T., Holland, D. M., Ilyina, T., Jungclaus, J. H., Komuro, Y., Krasting, J. P., Large, W. G., Marsland, S. J., Masina, S., McDougall, T. J., Nurser, A. J. G., Orr, J. C., Pirani, A., Qiao, F., Stouffer, R. J., Taylor, K. E., Treguier, A. M., Tsujino, H., Uotila, P., Valdivieso, M., Wang, Q., Winton, M., and Yeager, S. G.: OMIP contribution to CMIP6: experimental and diagnostic protocol for the physical component of the Ocean Model Intercomparison Project, Geosci. Model Dev., 9, 3231–3296, https://doi.org/10.5194/gmd932312016, 2016. a
Haas, C.: Dynamics Versus Thermodynamics: The Sea Ice Thickness Distribution, chap. 4, pp. 113–151, John Wiley & Sons, Ltd, ISBN 9781444317145, https://doi.org/10.1002/9781444317145.ch4, 2009. a
Holland, M. M., Bailey, D. A., Briegleb, B. P., Light, B., and Hunke, E.: Improved Sea Ice Shortwave Radiation Physics in CCSM4: The Impact of Melt Ponds and Aerosols on Arctic Sea Ice, J. Climate, 25, 1413–1430, https://doi.org/10.1175/JCLID1100078.1, 2012. a
Hunke, E. C. and Lipscomb, W. H.: CICE: the Los Alamos Sea Ice Model Documentation and Software User's Manual Version 4.0 LACC06012, Tech. Rep., Los Alamos National Laboratory, https://www.researchgate.net/publication/237249046_CICE_The_Los_Alamos_sea_ice_model_documentation_and_software_user's_manual_version_40_LACC06012 (last access: 9 September 2024), 2008. a
Hutter, N. and Losch, M.: Featurebased comparison of sea ice deformation in leadpermitting sea ice simulations, The Cryosphere, 14, 93–113, https://doi.org/10.5194/tc14932020, 2020. a
Inoue, J., Enomoto, T., Miyoshi, T., and Yamane, S.: Impact of observations from Arctic drifting buoys on the reanalysis of surface fields, Geophys. Res. Lett., 36, L08501, https://doi.org/10.1029/2009GL037380, 2009. a
International Arctic Buoy Programme: IABP 3hourly bouy data, International Arctic Buoy Programme [data set] https://iabp.apl.uw.edu/Data_Products/BUOY_DATA/3HOURLY_DATA/, last access: 9 September 2024. a
Kwok, R., Schweiger, A., Rothrock, D. A., Pang, S., and Kottmeier, C.: Sea ice motion from satellite passive microwave imagery assessed with ERS SAR and buoy motions, J. Geophys. Res.Oceans, 103, 8191–8214, https://doi.org/10.1029/97JC03334, 1998. a
Kwok, R., Hunke, E. C., Maslowski, W., Menemenlis, D., and Zhang, J.: Variability of sea ice simulations assessed with RGPS kinematics, J. Geophys. Res., 113, C11012, https://doi.org/10.1029/2008JC004783, 2008. a, b, c
Lemieux, J.F., Tremblay, B., Sedlác̆ek, J., Tupper, P., Thomas, S., Huard, D., and Auclair, J.P.: Improving the numerical convergence of viscousplastic sea ice models with the Jacobianfree NewtonKrylov method, J. Comput. Phys., 229, 2840–2852, https://doi.org/10.1016/j.jcp.2009.12.011, 2010. a, b
Lemieux, J.F., Lipscomb, W., Craig, A., Bailey, D. A., Hunke, E., Blain, P., Rasmussen, T. A. S., Bentsen, M., Dupont, F., Hebert, D., and Allard, R.: CICE on a Cgrid: new momentum, stress, and transport schemes for CICEv6.5, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd2023239, in review, 2024. a
Lindsay, R. W. and Stern, H. L.: The RADARSAT Geophysical Processor System: Quality of Sea Ice Trajectory and Deformation Estimates, J. Atmos. Ocean. Tech., 20, 1333–1347, https://doi.org/10.1175/15200426(2003)020<1333:TRGPSQ>2.0.CO;2, 2003. a
Marsan, D. and Weiss, J.: Space/time coupling in brittle deformation at geophysical scales, Earth Planet. Sc. Lett., 296, 353–359, https://doi.org/10.1016/j.epsl.2010.05.019, 2010. a
Marsan, D., Stern, H., Lindsay, R., and Weiss, J.: Scale Dependence and Localization of the Deformation of Arctic Sea Ice, Phys. Rev. Lett., 93, 178501, https://doi.org/10.1103/PhysRevLett.93.178501, 2004. a, b, c, d, e
Murray, R. J.: Explicit Generation of Orthogonal Grids for Ocean Models, J. Comput. Phys., 126, 251–273, https://doi.org/10.1006/jcph.1996.0136, 1996. a
National Snow and Ice Data Center (NSIDC): Sea Ice Index, Version 3, NSIDC [data set], https://noaadata.apps.nsidc.org/NOAA/G02135/, last access: 9 September 2024. a
Ning, C., and Xu, S.: OMARECommunity/LagrangiantrackingofseaiceinCICE (v1.0), Zenodo [code], https://doi.org/10.5281/zenodo.12616647, 2024. a, b
Oikkonen, A., Haapala, J., Lensu, M., Karvonen, J., and Itkin, P.: Smallscale sea ice deformation during NICE2015: From compact pack ice to marginal ice zone, J. Geophys. Res.Oceans, 122, 5105–5120, https://doi.org/10.1002/2016JC012387, 2017. a, b, c
Perovich, D., RichterMenge, J., Polashenski, C., Elder, B., Arbetter, T., and Brennick, O.: Sea ice mass balance observations from the North Pole Environmental Observatory, Geophys. Res. Lett., 41, 2019–2025, https://doi.org/10.1002/2014GL059356, 2014. a
Polar Science Center: PIOMAS (PanArctic Ice Ocean Modeling and Assimilation System) Ice Volume Data, 1979–present, Polar Science Center [data set], https://psc.apl.uw.edu/research/projects/arcticseaicevolumeanomaly/data/, last access: 9 September 2024. a
Rampal, P., Weiss, J., Marsan, D., Lindsay, R., and Stern, H.: Scaling properties of sea ice deformation from buoy dispersion analysis, J. Geophys. Res.Oceans, 113, C03002, https://doi.org/10.1029/2007JC004143, 2008. a
Rampal, P., Dansereau, V., Olason, E., Bouillon, S., Williams, T., Korosov, A., and Samaké, A.: On the multifractal scaling properties of sea ice deformation, The Cryosphere, 13, 2457–2474, https://doi.org/10.5194/tc1324572019, 2019. a, b, c, d, e, f
Rigor, I., ClementeColon, P., and Hudson, E.: The international Arctic buoy programme (IABP): A cornerstone of the Arctic observing network, OCEANS 2008, Quebec City, QC, Canada, 2008, 1–3, https://doi.org/10.1109/OCEANS.2008.5152136, 2008. a
Rigor, I. G., Wallace, J. M., and Colony, R. L.: Response of Sea Ice to the Arctic Oscillation, J. Climate, 15, 2648–2663, https://doi.org/10.1175/15200442(2002)015<2648:ROSITT>2.0.CO;2, 2002. a
Rothrock, D. A.: The energetics of the plastic deformation of pack ice by ridging, J. Geophys. Res., 80, 4514–4519, https://doi.org/10.1029/JC080i033p04514, 1975. a, b
Schweiger, A., Lindsay, R., Zhang, J., Steele, M., Stern, H., and Kwok, R.: Uncertainty in modeled Arctic sea ice volume, J. Geophys. Res.Oceans, 116, C00D06, https://doi.org/10.1029/2011JC007084, 2011. a
Small, R. J., Bacmeister, J., Bailey, D., Baker, A., Bishop, S., Bryan, F., Caron, J., Dennis, J., Gent, P., Hsu, H.m., Jochum, M., Lawrence, D., Munoz, E., diNezio, P., Scheitlin, T., Tomas, R., Tribbia, J., Tseng, Y.h., and Vertenstein, M.: A new synoptic scale resolving global climate simulation using the Community Earth System Model, J. Adv. Model. Earth Sy., 6, 1065–1094, https://doi.org/10.1002/2014MS000363, 2014. a
Smith, G. C., Roy, F., Reszka, M., Surcel Colan, D., He, Z., Deacu, D., Belanger, J.M., Skachko, S., Liu, Y., Dupont, F., Lemieux, J.F., Beaudoin, C., Tranchant, B., Drevillon, M., Garric, G., Testut, C.E., Lellouche, J.M., Pellerin, P., Ritchie, H., Lu, Y., Davidson, F., Buehner, M., Caya, A., and Lajoie, M.: Sea ice forecast verification in the Canadian Global Ice Ocean Prediction System, Q. J. Roy. Meteor. Soc., 142, 659–671, https://doi.org/10.1002/qj.2555, 2016. a
Sumata, H., de Steur, L., Divine, D. V., Granskog, M. A., and Gerland, S.: Regime shift in Arctic Ocean sea ice thickness, Nature, 615, 443–449, https://doi.org/10.1038/s4158602205686x, 2023. a, b
Tsamados, M., Feltham, D. L., and Wilchinsky, A. V.: Impact of a new anisotropic rheology on simulations of Arctic sea ice, J. Geophys. Res.Oceans, 118, 91–107, https://doi.org/10.1029/2012JC007990, 2013. a
Turner, A. K., Hunke, E. C., and Bitz, C. M.: Two modes of seaice gravity drainage: A parameterization for largescale modeling, J. Geophys. Res.Oceans, 118, 2279–2294, https://doi.org/10.1002/jgrc.20171, 2013. a
University Corporation for Atmospheric Research (UCAR): Community Earth System Model 2 (CESM2), UCAR [code], https://www.cesm.ucar.edu/models/cesm2, last access: 11 September 2024. a
Vancoppenolle, M., Rousset, C., Blockley, E., Aksenov, Y., Feltham, D., Fichefet, T., Garric, G., Guemas, V., Iovino, D., Keeley, S., Madec, G., Massonnet, F., Ridley, J., Schroeder, D., and Tietsche, S.: SI3, the NEMO Sea Ice Engine, Zenodo [code], https://doi.org/10.5281/zenodo.7534900, 2023. a
Wang, Q., Ilicak, M., Gerdes, R., Drange, H., Aksenov, Y., Bailey, D. A., Bentsen, M., Biastoch, A., Bozec, A., Boning, C., Cassou, C., Chassignet, E., Coward, A. C., Curry, B., Danabasoglu, G., Danilov, S., Fernandez, E., Fogli, P. G., Fujii, Y., Griffies, S. M., Iovino, D., Jahn, A., Jung, T., Large, W. G., Lee, C., Lique, C., Lu, J., Masina, S., Nurser, A. G., Rabe, B., Roth, C., Salas y Melia, D., Samuels, B. L., Spence, P., Tsujino, H., Valcke, S., Voldoire, A., Wang, X., and Yeager, S. G.: An assessment of the Arctic Ocean in a suite of interannual COREII simulations. Part I: Sea ice and solid freshwater, Ocean Model., 99, 110–132, https://doi.org/10.1016/j.ocemod.2015.12.008, 2016. a
Weiss, J.: Intermittency of principal stress directions within Arctic sea ice, Phys. Rev. E, 77, 056106, https://doi.org/10.1103/PhysRevE.77.056106, 2008. a
Xu, S., Ma, J., Zhou, L., Zhang, Y., Liu, J., and Wang, B.: Comparison of sea ice kinematics at different resolutions modeled with a grid hierarchy in the Community Earth System Model (version 1.2.1), Geosci. Model Dev., 14, 603–628, https://doi.org/10.5194/gmd146032021, 2021. a, b, c, d, e, f, g, h, i
Yang, C.Y., Liu, J., and Xu, S.: Seasonal Arctic Sea Ice Prediction Using a Newly Developed Fully Coupled Regional Model With the Assimilation of Satellite Sea Ice Observations, J. Adv. Model. Earth Sy., 12, e2019MS001938, https://doi.org/10.1029/2019MS001938, 2020. a
Zhang, S., Xu, S., Fu, H., Wu, L., Liu, Z., Gao, Y., Wan, W., Xie, J., Wan, L., Lu, H., Duan, J., Liu, Y., Fei, Y., Guo, X., Wang, Z., Wang, X., Wang, Z., Qu, B., Li, M., Zhao, H., Jiang, Y., Yang, G., Wang, H., An, H., Zhang, X., Zhang, Y., Ma, W., Yu, F., Xu, J., and Shen, X.: Toward Earth System Modeling with Resolved Clouds and Ocean Submesoscales on Heterogeneous ManyCore HPCs, Natl. Sci. Rev., 10, nwad069, https://doi.org/10.1093/nsr/nwad069, 2023. a, b
 Abstract
 Introduction
 Lagrangian tracking in CICE
 Numerical experiments and analysis
 Summary and discussion
 Appendix A: Model configuration and parameters of CICE
 Appendix B: Scaling analysis with Lagrangian tracking and Eulerian fields
 Appendix C: Major parameters of Lagrangian tracking
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Lagrangian tracking in CICE
 Numerical experiments and analysis
 Summary and discussion
 Appendix A: Model configuration and parameters of CICE
 Appendix B: Scaling analysis with Lagrangian tracking and Eulerian fields
 Appendix C: Major parameters of Lagrangian tracking
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement