the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
On the multifractal scaling properties of sea ice deformation
Pierre Rampal
Véronique Dansereau
Einar Olason
Sylvain Bouillon
Timothy Williams
Anton Korosov
Abdoulaye Samaké
In this paper, we evaluate the neXtSIM sea ice model with respect to the observed scaling invariance properties of sea ice deformation in the spatial and temporal domains. Using an Arctic setup with realistic initial conditions, stateoftheart atmospheric reanalysis forcing and geostrophic currents retrieved from satellite data, we show that the model is able to reproduce the observed properties of this scaling in both the spatial and temporal domains over a wide range of scales, as well as their multifractality. The variability of these properties during the winter season is also captured by the model. We also show that the simulated scaling exhibits a space–time coupling, a suggested property of brittle deformation at geophysical scales. The ability to reproduce the multifractality of this scaling is crucial in the context of downscaling model simulation outputs to infer sea ice variables at the subgrid scale and also has implications for modeling the statistical properties of deformationrelated quantities, such as lead fractions and heat and salt fluxes.
The fact that sea ice deformation maps look similar at different scales, with highly localized deformation features intersecting with a wide range of intersection angles (e.g., Hutchings et al., 2005; Wang, 2007; Hutter et al., 2019), suggests scale invariance in the spatial domain (Erlingsson, 1988). We note that scale invariance in space is also observed in sea ice for other deformationrelated quantities, such as floe sizes (Rothrock and Thorndike, 1984; Matsushita, 1985) and keel profiles (Rothrock and Thorndike, 1980). Comprehensive datasets of sea ice drift are now available at different spatial and temporal resolutions, from 50 m every 10 min (Oikkonen et al., 2017), 400 m every 2 d (Thomas et al., 2004, 2007, 2009) to 5–10 km every 3 d (Kwok, 2001; Stern and Moritz, 2002). Analyses of these datasets have confirmed the presence of scale invariance and, in particular, have confirmed that sea ice deformation is highly localized in both space and time.
In the spatial domain, deformation is observed as being concentrated along quasi onedimensional, socalled linear kinematic features (LKFs) organized in “weblike arrays” (Kwok et al., 1998; Thomas et al., 2007) that can be clearly identified over a wide range of space scales (Thomas et al., 2007; Linow and Dierking, 2017; Hutter et al., 2019). Sea ice deformation also appears to be a selfsimilar, highly localized process in the time domain. Isolated, shortduration fracturing events of various intensities occur over a wide range of frequencies. These events also sustain largerscale deformation, maintaining the LKFs as “active” for many days (Coon et al., 2007). The reorganization and formation of new LKFs occur in response to changes in the largescale atmospheric forcing (Kwok, 2001), and permanent deformation with high deformation rates in the ice pack is mainly synchronous with high wind events (Oikkonen et al., 2017).
A quantitative indication of scale invariance in sea ice deformation is given by the shape of the distribution of deformation rate invariants (i.e., shear and divergence) and the total deformation rates, which we refer to here as $\dot{\mathit{\epsilon}}$. These probability density functions (P) have indeed been shown to be “heavytailed”, i.e., dominated by extreme values, following a powerlaw decay of the type
where γ is an exponent larger than 1 that depends on the spatial and timescale considered (Lindsay and Stern, 2003; Marsan et al., 2004; Rampal et al., 2008; Hutchings et al., 2011; Bouillon and Rampal, 2015b). This important characteristic expresses scale invariance, as it is impossible from a powerlaw distribution to determine the scale of a given deformation even by comparing the relative number of deformation events of different sizes.
Localization in the time and space domains is revealed by scaling analysis of the deformation rate invariants. In such analyses, deformation rates are estimated at different spatial and temporal scales, by such methods as “coarsegraining” (see Sect. 3 for more details). The mean deformation rate is estimated using coarsegraining analysis (e.g., Lindsay and Stern, 2003; Marsan et al., 2004; Bouillon and Rampal, 2015b) or dispersion analysis of a pair of buoys (Rampal et al., 2008); the mean has been shown to vary with the spatial scale, L, and temporal scale of observation, T, as
respectively, hence following a power law. The scaling exponents, β and α, are both equal to or greater than zero and quantify the degree of localization of the deformation. In the space domain, β=0 characterizes the homogeneous deformation of an elastic solid or viscous fluid, i.e., a deformation that does not depend on the spatial scale, while β=2, i.e., the topological dimension for the 2Dlike sea ice cover, corresponds to a single “point” concentrating all of the deformation in an otherwise undeformed material (Rampal et al., 2008). Conversely, in the time domain, α=0 corresponds to a homogeneous deformation and a single, temporally isolated deformation event corresponds to the limit of α=1 (Rampal et al., 2008). This scaling has been shown to hold over a very wide range of space scales and timescales (Rampal et al., 2008; Oikkonen et al., 2017; Weiss, 2017), with that α and β larger than 0, even for timescales on the order of the winter season and for space scales on the order of the length of the Arctic basin. This result indicates the absence of a characteristic timescale and/or space scale for mean sea ice deformation over these scales and, as a consequence, that sea ice deformation can not be assumed to be homogeneous over timescales or space scales relevant for simulations of the Arctic system.
The fact that sea ice deformation rates are characterized by a heavytailed statistical distribution, i.e., dominated by extreme events, also indicates that the mean (moment of order 1) is not a sufficient quantity to describe the distribution of deformation rates at a given timescale or space scale. Higher moments of the distribution of deformation rates, such as the variance (order 2) and skewness (order 3), should indeed be explored to better describe the distribution and the associated process of sea ice deformation and be considered in temporal and scaling analyses, as proposed in this study.
While the value of the scaling exponents β and α for the first moment (the mean) describes the rate at which the magnitude of deformation events decreases with the scale of observation, it is the change in the value of β and α with respect to the moment q of the distribution that indicates how the temporal and spatial localization itself changes with the magnitude of deformation events. This change can be described by socalled structure functions of the form
in space and time, respectively. In the case of a linear structure function, i.e., no curvature or equivalent to a=0 or c=0, the amount of localization of large and small deformation events is the same and the scaling is said to be monofractal. In the case where both coefficients a and b or c and d are positive, the structure functions are convex, meaning that the higherorder moments of the distribution therefore increase much faster than the lowerorder moments with a decreasing scale of observation. In other words, large deformation events are more localized in time and space than smaller events, corresponding to the definition of a multifractal scaling (e.g., Kolmogorov, 1962; Lovejoy and Schertzer, 2007). Note that in the literature multifractality is also called intermittency when present in the time dimension and heterogeneity when present in the spatial dimension. The larger the curvature, the stronger the degree of multifractality of the scaling.
Spatial scaling analysis of sea ice deformation retrieved from radar or buoy drift data show a clear multifractal scaling expressed by a powerlaw scaling of the first, second and third moments, ranging from the resolution of the data up to hundreds of kilometers (e.g., Marsan et al., 2004; Rampal et al., 2008; Hutchings et al., 2011; Bouillon and Rampal, 2015a). Recently, Weiss and Dansereau (2017) suggested, based on the combination of all available data, including the ones of Oikkonen et al. (2017), that this multifractality also holds in the time domain over a period of 3 to 160 d. We note that multifractality in space has also been observed for openwater densities (Weiss and Marsan, 2004) and lead fractions and in time for shear stress amplitudes (Weiss and Marsan, 2004) and principal stress directions (Weiss, 2008).
These properties of sea ice deformation imply that observations of these quantities available at large scales can be statistically related, i.e., downscaled, to the same quantities at smaller, unresolved scales. In the case of model simulations, downscaling of outputs could be particularly valuable to infer quantities at the subgrid and/or subtimestep scale. In this context, the capability of these properties to reproduce in monofractal versus multifractal manner becomes very important. Indeed, if one was to estimate the distribution of a variable at the subgrid scale based on a model that would not reproduce the observed multifractality but only a monofractality, the downscaled distribution of this variable would greatly underestimate extreme values.
The multifractal behavior of sea ice has been the subject of a large number of interesting studies and is hypothesized to be of significant importance for sea ice rheology (e.g., Weiss and Dansereau, 2017). Bouillon and Rampal (2015a) and Rampal et al. (2016) showed that previous versions of neXtSIM were capable of reproducing the spatial scaling and multifractal behavior of the ice, with a very weak temporal scaling reported by Rampal et al. (2016). Spreen et al. (2017) and Bouchat and Tremblay (2017) have used scaling analysis to investigate their respective viscous plastic models, without going into the full details of a multifractal analysis or considering the temporal scaling. Hutter et al. (2018), on the other hand, did a full multifractal analysis of the spatial and temporal scaling in a viscous plastic model. Their work shows that with a model running at ∼1 km resolution, they can reproduce reasonably good spatial scaling and multifractality down to the 10 km scale and up to 200 km; it is not shown how well the scaling holds down to the actual model resolution, and their spatial scaling does not hold beyond the 200 km scale. Their results for temporal scaling show some inconsistencies: the authors report a reasonably good temporal scaling when considering the full domain (but do not report on multifractality); however, in the smaller region covered by the Envisat Geophysical Processor System (EGPS) data, the estimated scaling exponents are significantly smaller than for observations. They also only report temporal scaling for up to 30 d. Hutter and Losch (2019) appear to improve on these results, but as this paper is still under review further detailing of their results is premature.
The observed selfsimilarity and multifractality in the deformation and related characteristics of sea ice actually poses great challenges to the development of sea ice models, in particular in the continuum framework. On the one hand, the momentum and evolution equations for sea ice properties are based on mean variables. On the other hand, however, the observed multifractality in sea ice deformation implies that there is not a clear separation of scales between the strain rate due to mesoscale (50–100 m) heterogeneities in the ice (leads, ridges, etc.) and the strain rates at the 10 to 100 km scale. Consequently, no scale appears appropriate for homogenizing sea ice motion and thereby define a mean velocity or deformation rate for model resolutions ranging from 50 to 100 km.
In the absence of a characteristic space scale and timescale for the sea ice deformation, perhaps the best a continuum framework for sea ice modeling can do is to correctly reproduce the statistics of deformation from the smallest scales resolved (the nominal scale) to the largest scale, i.e., from the resolution of the grid in space and the model time step in time to the size of the Arctic basin and the timescale of a season. This is one of the motivation in developing neXtSIM, the numerical model used in the current study. Such localization at the nominal scale is the most faithful representation of the discontinuous nature of sea ice possible in a continuum model. Knowing the importance of essentially discontinuous features, such as leads, for atmosphere–ocean interactions modulated by sea ice (e.g., Smith, 1974; Kozo, 1983; Esau, 2007; Marcq and Weiss, 2012), we can expect the effect of using an ice model that localizes features at the nominal scale to be essential for improving the representation of this interaction in a coupled system.
This paper is the last step in validating neXtSIM against sea ice deformation statistics. While previous works have shown that the model reproduces the observed scaling of sea ice deformation (Bouillon and Rampal, 2015a; Rampal et al., 2016) in space, the temporal scaling and multifractality of both types of scaling have not yet been demonstrated for this model. The comparison performed here is based on satellite observations of sea ice deformation and winterlong simulations over the Arctic Ocean.
Section 2.1 and 2.2 discuss the recent developments of neXtSIM, the simulation setup and the observations. Section 3 describes the methodology used to perform the multifractal scaling analyses on both the model and observational data. Results of these analyses are presented in Sect. 4 and discussed in Sect. 5.
2.1 Model and simulation setup
neXtSIM is a finite elements sea ice model that uses a moving Lagrangian mesh. Its original dynamical component was based on the elastobrittle (EB) mechanical framework of Amitrano et al. (1999), first implemented in the context of sea ice by Girard et al. (2011), to account for brittle fracturing processes and the associated spatial localization of deformation. This framework was later adapted by Bouillon and Rampal (2015b) and Rampal et al. (2016) for longterm simulations of Arctic sea ice cover, including thermodynamical processes and advection, using a Lagrangian treatment of the equations of motion and a dynamical remeshing scheme. Yearlong simulations were presented in Rampal et al. (2016) and evaluated with respect to sea ice area, extent, volume, drift and deformation. In particular, the simulated deformation rates were demonstrated to be in good agreement with observations on the basis of their scaling properties in space.
However, the elastobrittle model does not, by definition, include a physical mechanism for irreversible deformations, as it is based on a strictly linearelastic constitutive law. It therefore cannot represent the transition between the small, elastic deformations associated with the fracturing of the ice cover and the permanent, potentially large, postfracturing deformation that dissipates internal stresses. It is therefore not suited to represent the dynamical behavior of a fractured ice cover over long (> day) timescales and cannot represent fully the properties of sea ice deformation in time.
The recent Maxwell elastobrittle (MEB) rheology addresses this limitation of the EB framework by including a mechanism for the relaxation of internal stresses that depend on the degree of fracturing of the sea ice cover (Dansereau et al., 2016). It is implemented in the current version of neXtSIM, which is used for this study. Another addition to the model is the introduction of a threethicknesscategory scheme that explicitly represents the thin and newly formed ice. The other model components (thermodynamics, slab ocean, etc.) remain unchanged relative to the version presented by Rampal et al. (2016).
All of the relevant equations of the current version of neXtSIM are presented in the appendices (Sect. A1 for the dynamical core and Sect. A2 for the threethicknesscategory scheme and sea ice thermodynamics). The numerics (spatial and temporal discretizations, advection scheme and numerical solvers) are the same as described in Rampal et al. (2016).
The initial mesh is generated in preprocessing over a panArctic region by using the mesh generator presented in Remacle and Lambrechts (2016), with a prescribed mean resolution (i.e., mean length of the edges of the triangular elements) of 10 km. The coasts are defined from the Global Selfconsistent, Hierarchical, Highresolution Geography Database (GSHHS_f_L1.shp, downloaded from: https://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/gshhgshp2.3.51.zip, last access: 1 February 2017). The domain is restricted to the central Arctic by putting open boundaries on the lines cutting through the Bering Strait from (67.7^{∘} N, 166.0^{∘} W) to (65.7^{∘} N, 170.7^{∘} W); cutting through the Canadian Arctic Archipelago from (76.7^{∘} N, 59.0^{∘} W) to (69.5^{∘} N, 121.0^{∘} W); and on the twosegment line cutting through the Greenland, Barents, and Kara Seas by joining the coordinates (77.0^{∘} N, 19.0^{∘} W), (73.0^{∘} N, 11.0^{∘} E), (72.9^{∘} N, 22.0^{∘} E), (76.1^{∘} N, 43.9^{∘} E), (75.7^{∘} N, 75.4^{∘} E), and (73.6^{∘} N, 88.5^{∘} E). We checked that using a larger domain with open conditions much further from the zone of interest does not impact the results presented in this paper.
The atmospheric forcing consists of the 10 m wind velocity, 2 m air temperature, mixing ratio, mean sea level pressure, total precipitation amount and snow fraction, and incoming shortwave and longwave radiation. All of these quantities are provided as 3hourly means and on a 30 km spatial resolution grid from the atmospheric state of the Arctic System Reanalysis^{1} (Bromwich et al., 2016).
The iceocean surface stress is computed from monthly ocean surface geostrophic currents derived, following Armitage et al. (2017), from the Arctic sea surface height data obtained from altimeters by Armitage et al. (2016). The provided surface height fields have a hole of missing data around the North Pole that we filled using a linear interpolation between the northernmost available points and their mean. A smoother is applied to the ocean velocity components in the filled area to avoid spurious oscillations due to the interpolation method. The final ocean current forcing is at a spatial resolution of 50 km. The slab ocean salinity and temperature are nudged towards the daily sea surface temperature and salinity data provided as daily means and at 12.5 km spatial resolution over the Arctic region by the TOPAZ4 reanalysis (available at: http://marine.copernicus.eu/servicesportfolio/accesstoproducts/, last access: May 2018) (Sakov et al., 2012) with a nudging timescale equal to 30 d. TOPAZ4 is a coupled ocean and sea ice data assimilation system for the North Atlantic and the Arctic that is based on the HYCOM ocean model and the ensemble Kalman filter data assimilation method using 100 dynamical members. A 23year reanalysis has been completed for the period 1991–2013 and is the multiyear physical product in the Copernicus Marine Environment Monitoring Service. The ocean depth, H, used for the basal stress parametrization comes from the 1 arcmin ETOPO1 global topography (available at: https://www.ngdc.noaa.gov/mgg/global/, last access: May 2018) (Amante and Eakins, 2009).
Our reference simulation starts on 15 November 2006. The level of damage in the ice cover (see Appendix A1) is initially set to zero where sea ice is present. Initial sea ice concentration and thickness are set from a combination of the TOPAZ4 reanalysis and the OSISAF climate data record (Tonboe et al., 2016) and ICESAT (available at: https://icdc.cen.unihamburg.de/1/daten/cryosphere/seaicethicknesssatobsarc.html, last access: December 2017) (Kwok et al., 2009) datasets, respectively.
2.2 Satellite observations
We use the Lagrangian displacement data produced by the RADARSAT Geophysical Processor System (RGPS) as described in Kwok et al. (1998). This dataset covers the western Arctic for the period 1996–2008 and provides trajectories of sea ice points initially located on a 10 km regular grid (http://rkwok.jpl.nasa.gov/radarsat/lagrangian.html, last access: December 2017). The positions of these points are updated when two successive synthetic aperture radar (SAR) images are available. The time interval between two updates is typically 3 d. For the present analysis we use the data covering the winter period 3 December 2006 to 30 April 2007 from the reprocessed RGPS Lagrangian displacement product, the socalled RGPS Image Pair Product, introduced and used in Bouillon and Rampal (2015b) (Sect. 2.2).
Scaling analyses of sea ice deformation can be performed using two approaches: a socalled coarsegraining method as in, e.g., Marsan et al. (2004) or buoys dispersion method as in Rampal et al. (2008) (using pairs of buoys) or Oikkonen et al. (2017) (using triplets of buoys). We use a similar method as Oikkonen et al. (2017), i.e., computing velocity gradients from the trajectories of triplets of points. The nominal resolution for a scaling analysis is defined as the square root of the surface area of the polygon considered. For example, the minimal spatial resolution that can be achieved with the RGPS dataset when using the threesided polygons obtained from a Delaunay triangulation is about 7.5 km. This also set the nominal spatial resolution of the analysis presented in this study.
Drifters in the model are seeded at the location of the RGPS grid points as of 3 December 2016. The RGPS grid for this initialization is undeformed and the points are regularly spaced by 10 km. The positions of the simulated drifters are updated at each model time step until the end of the simulation or until the ice concentration drops to zero (through melting or opening of a lead). Both the RGPS and simulated trajectories are filtered for the presence of coasts, with a proximity threshold of 100 km. Only the trajectories spanning the same time periods in both the simulation and RGPS dataset are considered in the calculation of the deformation and their statistics. This selection led to discarding only about 1 % of the total trajectory dataset and does not affect the results of the analyses presented in this paper. However, we apply this selection in order to make our comparison between model and observations as consistent and clean as possible.
Triplets of drifting points are defined as the result of Delaunay triangulation of the initial positions of the tracked RGPS points, which ensures that the associated polygons are independent, i.e., nonoverlapping. The exact same triplets of points are considered in the model for the analysis, meaning that we follow the exact same set of triplets of trajectories (or triangles) in the model and in the observations. The polygons after initiation are defined by the positions of their three nodes at any given time. We stress the fact that the simulated trajectories are not reinitialized every 3 d to match the RGPS positions; only the initial positions are identical between the RGPS and the model trajectories.
To perform a spatial scaling analysis of sea ice deformation, one needs to consider triplets of points with different spacing, i.e., different sizes of polygons. In order to obtain sets of polygons of different surface areas, we perform successive Delaunay triangulation through the clouds of points defined by the initial positions of the RGPS points, using increasingly subsampled clouds of these points. Each set of contiguous polygons obtained using this process is associated to a spatial scale, L, defined as the mean of the square root of the polygon surface areas obtained from the triangulation, i.e., from 7.5 to 580 km in this study. We note that due to the finite size of the Arctic basin and the largest spatial scale of 580 km considered here, the number of triplets available for the statistical analyses decreases by a factor of (570/7.5)^{2} as the space scale increases from 7.5 to 580 km.
To perform a temporal scaling analysis of sea ice deformation, one also need to consider the positions of triplets of drifters separated by different times T. The number of available triplets for our analysis in the time domain therefore also decreases as the timescale considered increases due to the finite time covered by our simulation (about 5 months), which is constrained by the fact that we wish to limit our analysis to the winter period, i.e., from early December to midApril.
For each available polygon, the total deformation rate is calculated as follows:
where ${\dot{\mathit{\epsilon}}}_{\mathrm{shear}}$ and ${\dot{\mathit{\epsilon}}}_{\mathrm{div}}$ are the two invariant, shear and divergence, respectively, of the deformation rate. These invariant are estimated using a contour integral calculation as follows: The spatial derivatives of the components u and v of the velocity calculated at a given timescale T are obtained by calculating the contour integrals as in Kwok et al. (2008) and Bouillon and Rampal (2015b) around the boundary of each polygon associated to a given space scale L:
where A is the encompassed area of the polygon approximately equal to L^{2}. For example, u_{x} is approximated by
where n=3 and subscript $n+\mathrm{1}=\mathrm{1}$. The shear rates ${\dot{\mathit{\epsilon}}}_{\mathrm{shear}}$ and divergence rates ${\dot{\mathit{\epsilon}}}_{\mathrm{div}}$ are then computed as follows:
The distribution of total deformation rates is constructed from each given coupled space scale–timescale (L, T), and their first three moments are calculated as $\langle {\dot{\mathit{\epsilon}}}_{\mathrm{tot}}^{q}\rangle $, where $q=\mathrm{1},\mathrm{2},\mathrm{3}$ is the moment order.
Below we discuss some issues that are inherent to the data and our method and their impact in terms of the robustness of the statistics calculated here.

With time, the triangular elements can become too distorted, in which case their length scale, L, is poorly defined. Applying a test for distortion based on the smallest angle of the polygons and discarding the most distorted ones was found to affect the results in terms of the slope of the scaling and the goodness of the fit of the powerlaw fit of the scaling. Hence, here we discard from the analysis the polygons having a minimum angle of 30^{∘} or fewer.

The RGPS trajectories are not sampled at regular time intervals, contrary to the model, due to the irregular interval between the two satellite orbits. The mean sampling is of about 3 d, and 90 % of trajectories are sampled with a frequency between 2.5 and 3 d. Because sea ice deformation depends on the timescale (see results of Sect. 4.2), one should make sure to use similar sampling times for the observations and the model when computing and comparing deformation rate estimates. To deal with this issue, we performed a subsampling of the RGPS trajectory dataset using a nearestneighbor interpolation of the original positions at 3 d intervals but only when the RGPS drifter's position is available within ±6 h around the interpolation target time. The positions simulated by the model, which are outputted every 3 h from midnight to midnight each day, are taken to match the subsampled RGPS time series obtained as described above.

The 3 d RGPS sampling additionally places a lower bound on the timescales one can explore when comparing the simulated and observed deformation rates. In the present analysis, we therefore restrict ourselves to timescales equal to or greater than 3 d.
We find that the relative number of available polygons is what has the largest impact on the deformation statistics. Some facts therefore need to be kept in mind when performing a scaling analysis over a finite period of time. In the time domain in particular this means that sea ice deformation is better sampled, i.e., more triplets are available, for the early than for the late part of the period and for short timescales T than for large timescales. In the present case, the computed statistics are therefore more representative of early than late winter.
Figure 1 shows the maps of the 3 d total, shear and absolute divergence deformation rates simulated by the model and estimated from the RGPS data at the same locations. Note that to obtain a better spatial coverage, these maps show all simulated or observed 3 d deformation rates for a period of 7 d centered on 4 February 2007. The probability density functions of the simulated and observed total, shear, and absolute divergence deformation rates from the snapshots of Fig. 1 are shown in Fig. 2. All distributions exhibit a powerlaw tail with almost identical slopes of about −3, similar to what, e.g., Marsan et al. (2004) found in their study, and with a remarkable agreement between the model and the observations for each invariant of the deformation. From the statistical point of view, this implies that one needs to consider higher moments than the mean of the distributions to fully describe the statistics of the sea ice deformation process (Sornette, 2006). In the scaling analysis presented in the following sections, we thus systematically calculate the three first moments of the distributions of deformation rates.
4.1 Spatial scaling analysis
Figure 3a shows the winter mean of the spatial scaling analysis for the observations and model calculated for a T=3 d temporal scale. We found that both model and observations statistics follow power laws. We used logarithmically spaced bins and applied an ordinary leastsquares method to the binned data in log–log space to get a reasonably accurate estimate of the powerlaw fits (Stern et al., 2018). The deformation rates are very well captured by the model across scales. However, the first and second moments of the distributions are slightly overestimated by the model compared to the observations, whatever the spatial scale considered was. For example, at the nominal scale of 7.5 km, the first and second moments are overestimated by a factor of 2 compared to the observations. This may be due to one or several of the following factors: (1) inaccuracies in the atmospheric forcing, (2) our choice of mechanical parameters values (Weiss and Dansereau, 2017) or (3) the value of the atmospheric drag coefficient (e.g., Bouchat and Tremblay, 2017). It is especially important to note that the simulated deformation rate has not been tuned with respect to the MEB mechanical parameters in the present simulations. We consider such tuning to be out of the scope of this study, which focuses on the ability to reproduce the observed scaling (exponents of the power laws) and, in particular, their multifractal properties (nonlinearity of the structure function). The simulated and observed structure functions (i.e., the dependence of the scaling exponents on the moment order) β(q) are shown in Fig. 3b. The spatial scaling obtained for both the model and the RGPS are clearly multifractal, as their structure functions can be both approximated by a quadratic function as defined by Eq. (4) with coefficients a=0.11 and 0.13, respectively. These values, corresponding to the curvature of the structure functions, are very close to those reported in previous studies (a=0.13–0.14; Marsan et al., 2004; Rampal et al., 2016). It is important to note that the bars on this graph are estimated from the minimal and maximal local scaling exponent values, as in Bouillon and Rampal (2015a) and hence correspond to upperbound estimates rather than to error bars. This good agreement is a relevant indication that the scaling in the simulated deformation fields is consistent with that observed between 7.5 and 580 km.
Using successive and contiguous snapshots throughout the winter, a time series of the value of the spatial scaling exponent β obtained for the mean deformation rates (q=1) is calculated and plotted on Fig. 4. It shows that the spatial scaling exponent varies between −0.1 and −0.34. These exponents are in good agreement with the 1month running means of the scaling exponents calculated by Stern and Lindsay (2009) for the entire period covered by the RGPS dataset (1996–2008). The scaling exponent for the mean is about −0.2 on average over the whole winter period for the simulated and observed total deformation rates, which is also the value found by Marsan et al. (2004) for a snapshot of deformation rates calculated over a 3 d period centered on 6 November 1997. We also note that the model reproduces the observed variability of the scaling exponent throughout the whole period well.
We further characterize the properties of the spatial scaling for both the model and observations by exploring its dependence on the temporal scale, T. We find that the estimated spatial scaling exponent, β, decreases with increasing T, although this behavior is only obvious for the moments of order 2 and 3 (Figs. 5a and 6a). This seems to correspond to the existence of space–time coupling of the scaling properties of sea ice deformation. This property was originally suggested by Rampal et al. (2008) from the result of their scaling analysis of buoy pair dispersion and was further explained in Marsan and Weiss (2010) as being a possible characteristic of brittle deformation at geophysical scales. To our knowledge, this is the first time such coupling is obtained from a sea ice model simulation run at such relatively coarse spatial resolution. The origin of this coupling has been previously proposed to be linked to the complex correlation patterns related to chain triggering of icequakes (Marsan and Weiss, 2010). Further study is, however, needed to explore this hypothesis, which is out of the scope of this paper.
We also note a decrease in the multifractal character of the spatial scaling (i.e., the curvature of β(q)) when increasing the timescales from T=3 to T=96 d (Figs. 5b and 6b). For both the model and the observations, we observe that the multifractality property is present for all scales considered in this study, although a decrease in the degree of multifractality is observed as the scale increases. The curvature values decrease from 0.115 to 0.054 for the model and from 0.13 to 0.063 for the observations following a power law (Fig. 7). While the general behavior of decreasing the degree of multifractality of the spatial scaling as the timescale increases is captured by the model, we note that the degree of multifractality of the deformation is systematically underestimated by the model compared to the observations. This could either be attributed to an inaccurate position, a lack of extreme events in the atmospheric forcing, or an inadequate or insufficiently tuned parameterization of the damage healing in the model. In any case, the reason for this discrepancy should be further explored but is out of scope of the present paper.
4.2 Temporal scaling analysis
The results of the temporal scaling analysis for L=7.5 km is shown in Fig. 8a. We see a robust and very similar powerlaw scaling for the two first moments ($q=\mathrm{1},\mathrm{2}$) for both the model and observations between T=3 d (i.e., the temporal resolution of the observations) and T=96 d. In previous studies based on drifting buoy trajectories, whose positions are sampled at higher frequency, it has been suggested that the temporal scaling for the mean total deformation holds down to 1 h (Hutchings et al., 2011). A recent study based on veryhighresolution ship radar measurements has demonstrated that it holds down even to 10 min (Oikkonen et al., 2017). Here, we obtain a perfect agreement between the slope ($\sim \mathrm{0.3}$) of the temporal scaling for the mean deformation rates estimated by Oikkonen et al. (2017) and those estimated from the RGPS data and the present model simulations.
We note, however, that the third moments of the distributions are slightly underestimated by the model at all timescales. This means that the proportion of extreme deformation events compared to lower ones is too small or that their values are too low in the simulation. This may come from the inaccuracy of the relatively coarse (30 km) atmospheric reanalysis we use to force our model and that is known to poorly resolve the most extreme low pressure systems, a common shortcoming of all the available global or regional atmosphere reanalysis to date. Another explanation could be the fact that we have not tuned the MEB rheology parameters to reproduce the proportion of extreme deformation events versus the less extreme events. In this rheology, the coupling between the damage and the mechanical behavior of sea ice is nonlinear and it is therefore expected that varying parameter values can change the proportion of the simulated extreme events, i.e., the skewness of the distribution of deformation rates.
As in the spatial domain, the temporal scaling is found to be multifractal for the model and observations, and the match is remarkably good. The curvature values (i.e., the coefficient c in Eq. 5) are 0.67 for the model and 0.62 for the observations. This suggests that the multifractal character of the temporal scaling is stronger than the spatial scaling and possibly a robust property of sea ice deformation, at least in the wintertime, independent of the observed change in sea ice cover state and the associated shift of its dynamical regime during the period 1996–2006 (e.g., Rampal et al., 2009a, b). We note that the values of curvature of the structure functions obtained here cannot be directly compared to the one reported in Weiss and Dansereau (2017) since in their paper the authors are plotting the normalized moments $\langle {\dot{\mathit{\epsilon}}}_{\mathrm{tot}}^{q}{\rangle}^{\mathrm{1}/q}$ of the distribution versus the temporal scale instead of the actual moments $\langle {\dot{\mathit{\epsilon}}}_{\mathrm{tot}}^{q}\rangle $ as we do here.
We also investigate the dependence of the temporal scaling on the spatial scale of observation, L, for both the model and RGPS data (Figs. 9a and 10a). We find that the scaling exponent, α, decreases with L. Similar to the spatial scaling analysis performed in Sect. 4.1, we find here the signature of a spacetime coupling in the scaling properties of sea ice deformation. The multifractal character of the temporal scaling holds at all the spatial scales considered here (L=7.5 to L=580 km) and is similar in the model and observations (Figs. 9b and 10b). The curvature values are going from 0.67 down to 0.37 for the model and from 0.63 to 0.35 for the observations following a power law (Fig. 11). The decrease in the degree of multifractality of the temporal scaling as the space scale increases, as seen in the observations, is remarkably well captured by the model.
Our statistical analyses have shown that the neXtSIM model correctly reproduces the distribution of sea ice deformation rates, its scaling properties in both the space and time domains and its multifractal behavior. In particular, it is the first time that multifractality in the time domain is shown to be reproduced in a sea ice model.
The MEB rheology was developed with the aim of improving the representation of the physics of sea ice continuum models by including the ingredients hypothesized by Weiss and Dansereau (2017) to possibly play an important role in the emergence of multifractal heterogeneity and intermittency of sea ice deformation. This hypothesis is based on the analysis of observational data that have highlighted the existence of multifractality of sea ice deformation in space and time (Rampal et al., 2008; Bouillon and Rampal, 2015b; Weiss and Dansereau, 2017) and on a close and arguably sound analogy that can be made with other largescale solids sharing these properties, such as the Earth's crust (Weiss et al., 2009). According to Weiss and Dansereau (2017) the ingredients are a threshold mechanism for brittle fracturing, some disorder that represents the natural heterogeneity of the material at the mesoscale, longrange elastic interactions within the ice cover that promote avalanches of fracturing events through a cascading mechanism, postfracturing relaxation of elastic stresses through viscouslike relaxation and a slow restoring or healing mechanism of the sea ice mechanical properties. We argue that the results obtained here are at least a sound demonstration that a model including these ingredients can indeed reproduce complex aspects of sea ice deformation, such as its multifractality.
We show here that the spatial scaling of sea ice deformation simulated in a realistic setup by neXtSIM holds down to the nominal resolution of the mesh, a result that is in agreement with previous analyses of the MEB model in idealized simulations (Dansereau et al., 2016) and realistic ones (Rampal et al., 2016). It means that neXtSIM does not need to be run at higher spatial resolution in order to reproduce the observed scalings, as, e.g., Hutter et al. (2018) do when running at about 1 km resolution in order to resolve sea ice deformation at scale of about 10 km. Localizing the deformation at the nominal model resolution also means that related quantities, such as ridges, leads and linear kinematic features, are better resolved, although this is not investigated directly here. It is also important to note that using a Lagrangian scheme with adaptive remeshing instead of, for instance, a Eulerian finite volume scheme, helps in preserving such features once formed but plays no role in their formation.
We also show that this spatial localization and the multifractal character of the simulated mean sea ice deformation is resolutionindependent in the model (see Fig. 12). However, despite the fact that the scaling remains multifractal when neXtSIM runs at coarser resolutions (e.g., 15 or 30 km), the level of multifractality is decreasing with decreasing resolution. Indeed, the second and third moments of deformation rates from the 15 and 30 km runs differ from the results obtained from the 7.5 km run (Fig. 12b), which suggests an underestimation of extreme deformation events at the smaller spatial scales with increasing model resolution. Nevertheless, the representation of multifractality at all resolutions implies that neXtSIM could be adequately used to explore a wider range of space scales and timescales than those covered by the currently available observations of the global Arctic. In particular, it could allow us to “zoom in” and explore the statistical properties of sea ice deformation at subsatellite observation scales, which are of increasing interest for both regional to global climate modeling and operational forecasting. A model that could otherwise not represent the multifractal character of sea ice deformation and would only reproduce a monofractal scaling would greatly underestimate extreme deformation events and their impact on sea ice conditions, e.g., the presence (or lack thereof) of leads and ridges, at such scales.
A model that allows for reproducing sea ice deformation and its scaling properties down to its nominal resolution does not preclude the need for appropriate subgridscale parameterizations. On the contrary, we believe that physically sound parameterizations are indeed required and that the knowledge of the distribution of deformation rates at the subgrid scale made possible by neXtSIM could be highly valuable in terms of informing these parameterizations. An appropriate subgridscale parametrization should link the deformation simulated at the scale of the grid cell with the scale at which deformation really occurs within the ice cover, which is the size of individual leads and ridges.
Moreover, we argue that, as sea ice deformation is strongly tied to other model variables, such as drift, lead fraction and thickness distribution, a proper simulation of these variables is a necessary prerequisite to using models for investigating various coupled ocean–ice–atmosphere processes, and their impact on their immediate vicinity and on the polar climate system. For example, the accuracy of neXtSIM in reproducing the observed statistical properties of sea ice deformation as demonstrated in this paper is thought to go handinhand with its capability of representing the observed properties of lead fraction. This is the subject of a parallel study that is in preparation.
In this study we have compared the deformation rates simulated by neXtSIM to those derived from RGPS observations on the basis of their distributions and have shown how these distributions scale in time and space. The conclusions of our analysis are as follows.

The neXtSIM model reproduces the first, second and third moments of the statistical distribution of observed sea ice deformation rates and how it scales in space and time well. In particular, this is the first time the observed scaling invariance in the temporal domain (i.e., intermittency) of sea ice deformation is shown to be reproduced by a model on a realistic panArctic setup over such a large range of scales.

Sea ice deformation rates calculated over a temporal scale of 3 d scale in space from the scale of the model (mesh resolution) and observations up to about 700 km in a multifractal manner.

Sea ice deformation rates calculated over a spatial scale of 7.5 km scale in time over the range of 3 d to 3 months in a multifractal manner.

A space–time coupling in the scaling properties of sea ice deformation is shown to be reproduced by the model. This suggests that neXtSIM could be the proper tool for studying the physical meaning and origin of this coupling, in the context of brittle deformation of geophysical solids.

The simulated mean sea ice deformation rates and their associated scaling invariance characteristics are resolutionindependent, i.e., when running the neXtSIM model at resolutions of 7.5, 15 or 30 km. The most extreme deformation events may be missed, however, if running at coarser resolutions, i.e., the second and thirdorder moments may be underestimated compared to the highresolution run.

As the monofractal versus multifractal character of the scaling of deformation rates is the discriminating factor for the heterogeneity and intermittency of the deformation, we suggest that a multifractal scaling analysis could be considered a meaningful validation step before further analyzing sea ice model outputs that could be influenced by sea ice dynamics.

The good agreement between the model and observations motivates the use of neXtSIM as a tool to further investigate physical processes that are highly sensitive to sea ice deformation.
The datasets used for this study (forcing data, ocean topography, sea ice initial conditions) are all publicly available at the URL addresses mentioned in the main text. Outputs of the simulations analyzed in this study are not publicly accessible but are available upon request. The actual code of the numerical model neXtSIM, despite being hosted on Github (https://github.com/nansencenter/nextsim, last access: February 2019), is however not yet in openaccess, as it is still under development. First release of the code is planned for 2020.
This section presents the dynamical and thermodynamical components of neXtSIM. The waveinice module implemented by Williams et al. (2017) is not included here. Prognostic sea ice variables are listed in Table A1 and all parameter values are listed in Table A2.
A1 Dynamical core
The evolution equation for sea ice velocity comes from vertically integrating the horizontal sea ice momentum equation as follows:
The parameter ρ_{i} is the ice density; H is the mean ice thickness per unit grid cell area; σ is the sea ice internal stress tensor; and τ_{a}, τ_{w}, and τ_{b} are the surface wind, ocean, and basal stresses, respectively, as defined in Rampal et al. (2016). The parameter f is the Coriolis frequency, k is the upward pointing unit vector, g is the acceleration due to gravity and η is the ocean surface elevation. In the region with only thin ice or with thick ice thickness lower than a given threshold (defining our ice edge), the momentum equation is replaced by Laplace's equation so that the velocity linearly decreases from the ice edge to the nearest coast (see Samaké et al., 2017). The additional ice pressure term introduced in Rampal et al. (2016) is not included here.
Following Dansereau et al. (2016), the evolution equation for the internal stress takes the form of the Maxwell constitutive law:
where λ is the relaxation time for the stress, E, is the elastic modulus and $\dot{\mathit{\epsilon}}$, the strain rate tensor, is defined as
Plane stress conditions are assumed and the stiffness tensor K reads
where ν is Poisson's ratio.
As in Dansereau et al. (2016), both the elastic modulus, E, and the relaxation time are functions of the ice concentration, A, and the level of damage, d. The level of damage is a scalar, gridscale variable that represents the density of fractures at the subgrid scale. Its value is 0 for an undamaged and 1 for a “completely” damaged material, which we note is the opposite convention to that of Dansereau et al. (2016). The elastic modulus is a linear function of d:
where E_{0} is the undamaged elastic modulus and f(A) introduces a dependence on the ice concentration via the following exponential function:
where c^{∗} is the ice compactness parameter introduced by Hibler (1979). As in Dansereau et al. (2016), the relaxation time is a power function of d:
where λ_{0} is its undamaged value and α is a constant exponent greater than 1. Here, we use the values α=5 and λ_{0}=10^{7} s (115 d, as in the realistic MaxwellEB simulations of Dansereau et al., 2017) to ensure that the relaxation of stresses is virtually zero over an undamaged ice cover but is significant when the ice is damaged.
The evolution of the damage is controlled by the location of the predicted stress state relative to the failure envelope, which is defined in terms of the principal stress components
with the convention that compressive stresses are positive.
Here, the envelope combines a Mohr–Coulomb failure criterion and maximum tensile and compressive stress criteria. These are given by
where $q={\left[{\left({\mathit{\mu}}^{\mathrm{2}}+\mathrm{1}\right)}^{\mathrm{1}/\mathrm{2}}+\mathit{\mu}\right]}^{\mathrm{2}}$, ${\mathit{\sigma}}_{c}=\frac{\mathrm{2}c}{\left[{\left({\mathit{\mu}}^{\mathrm{2}}+\mathrm{1}\right)}^{\mathrm{1}/\mathrm{2}}\mathit{\mu}\right]}$, μ is the internal friction coefficient, c is the cohesion, σ_{T max} is the maximal tensile strength and σ_{N max} the maximum compressive strength (see Table A2). The cohesion, c, is scaled as a function of the model spatial resolution, as described in Bouillon and Rampal (2015a).
When one of the damage criteria is met, d is modified by multiplying (1−d) with Ψ, or
where
Healing is included here to represent the counteracting effect of refreezing of water within leads on the level of damage of the ice cover. It is implemented via a constant term in the damage evolution equation:
where T_{h} is the characteristic time for healing and T_{d} is the characteristic time for damaging (Dansereau et al., 2016).
A2 Ice thickness redistribution and thermodynamics
neXtSIM includes a multicategory model inspired from Stern and Rothrock (1995), i.e., considering three categories: thick ice, thin ice and open water. In our implementation the thin ice is only newly formed ice, so ice will only be transferred from the thin ice category to thick ice but not in the reverse direction. In addition, we do not apply additional openwater source terms or use the formulation of Gray and Morland (1994) to keep the ice concentration less than 1. (We simply redistribute ice and snow volume if this occurs.) Thin ice is described by its concentration, A_{t}, volume per unit area, H_{t}, and snow volume per unit area, h_{s,t}. Thick ice is described by the concentration, A, volume per unit area, H, and snow volume per unit area, h_{s}. We assume that the thin ice has no mechanical strength and simply follows the motion of the surrounding thick ice.
Note that the total ice concentration and volume per unit area are A+A_{t} and H+H_{t} and the total snow volume per unit area is ${h}_{\mathrm{s}}+{h}_{\mathrm{s},\mathrm{t}}$.
Thin ice thickness is considered to be uniformly distributed with thickness ${h}_{\mathrm{t}}={H}_{\mathrm{t}}/{A}_{\mathrm{t}}$ required to be between h_{min}=5 cm and h_{max}=27.5 cm. The evolution equations for A, H, h_{s}, A_{t}, H_{t} and h_{s,t} have the following form:
where $\frac{D\mathit{\varphi}}{Dt}$ is the material derivative that is defined for any scalar as
Here ∇⋅u is the divergence of the horizontal velocity, Ψ_{ϕ} is a sink or source term due to ridging, and S_{ϕ} is a thermodynamical sink or source term. Volume conservation is imposed by setting ${\mathrm{\Psi}}_{H}={\mathrm{\Psi}}_{{H}_{\mathrm{t}}}$ and ${\mathrm{\Psi}}_{{h}_{\mathrm{s}}}={\mathrm{\Psi}}_{{h}_{\mathrm{s},\mathrm{t}}}$ and an additional constraint is that ${A}_{\mathrm{t}}+A\le \mathrm{1}$.
The evolution of A, H, A_{t} and H_{t} is computed following three main steps (variables updated in each step are denoted with a prime):

Advection. The scheme solves the equation:
$$\begin{array}{}\text{(A18)}& {\displaystyle \frac{D\mathit{\varphi}}{Dt}}=\mathit{\varphi}\mathrm{\nabla}\cdot \mathit{u},\end{array}$$for each conserved scalar quantity (A, H, A_{t}, H_{t}, etc.). For this paper, we use the purely Lagrangian scheme presented in Rampal et al. (2016). After this step the concentration could be larger than 1.

Mechanical redistribution. The scheme imposes the limit ${A}_{\mathrm{t}}+A\le \mathrm{1}$ on the total ice area by following the steps below.

Compute the new openwater concentration as follows:
$$\begin{array}{}\text{(A19)}& {A}_{\mathrm{0}}=max(\mathrm{0},\mathrm{1}A{A}_{\mathrm{t}}),\end{array}$$a source term for the open water could be added here (as in Stern and Rothrock, 1995) to represent subgridscale convergence and divergence.

Compute the new thin ice concentration as follows:
$$\begin{array}{}\text{(A20)}& {A}_{{\mathrm{t}}^{\prime}}=max(\mathrm{0},min(\mathrm{1},\mathrm{1}A{A}_{\mathrm{0}}\left)\right).\end{array}$$ 
Compute the transfer of thin ice if ${A}_{{\mathrm{t}}^{\prime}}<{A}_{\mathrm{t}}$ by setting the following equations:
$$\begin{array}{}\text{(A21)}& {\displaystyle}{H}_{{\mathrm{t}}^{\prime}}& {\displaystyle}={H}_{\mathrm{t}}{\displaystyle \frac{{A}_{{\mathrm{t}}^{\prime}}}{{A}_{\mathrm{t}}}},\text{(A22)}& {\displaystyle}{h}_{\mathrm{s},\mathrm{t}}^{\prime}& {\displaystyle}={h}_{\mathrm{s},\mathrm{t}}{\displaystyle \frac{{A}_{{\mathrm{t}}^{\prime}}}{{A}_{\mathrm{t}}}},\text{(A23)}& {\displaystyle}{H}^{\prime}& {\displaystyle}=H+({H}_{\mathrm{t}}{H}_{{\mathrm{t}}^{\prime}}),\text{(A24)}& {\displaystyle}{h}_{{\mathrm{s}}^{\prime}}& {\displaystyle}={h}_{\mathrm{s}}+({h}_{\mathrm{s},\mathrm{t}}{h}_{\mathrm{s},{\mathrm{t}}^{\prime}}),\text{(A25)}& {\displaystyle}\mathrm{\Delta}A& {\displaystyle}={\displaystyle \frac{{A}_{\mathrm{t}}{A}_{{\mathrm{t}}^{\prime}}}{\mathit{\zeta}}}.\end{array}$$Here, we have transferred ice and snow volume from thin to thick ice in a conservative manner, but we will not try to conserve ice area; ζ is an aspect ratio parameter (tuned to 10), which causes ridging to preferentially increase ice thickness over ice area.

Compute the new thick ice concentration as follows:
$$\begin{array}{}\text{(A26)}& {A}^{\prime}=max(\mathrm{0},min(\mathrm{1},\mathrm{1}{A}_{{\mathrm{t}}^{\prime}}{A}_{\mathrm{0}}+\mathrm{\Delta}A\left)\right).\end{array}$$ 
Apply more ridging if $({A}^{\prime}+{A}_{{\mathrm{t}}^{\prime}})>\mathrm{1}$ by setting ${A}^{\prime}=\mathrm{1}{A}_{{\mathrm{t}}^{\prime}}$.


Growth and melt. The source and sink terms from the thermodynamics are computed by applying the zerolayer Semtner (1976) vertical thermodynamics to the new ice category and that of Winton (2000) for the thick ice, as if the thickness was uniform and equal to H∕A for the thick ice and H_{t}∕A_{t} for the thin ice. Freezing of open water is computed as in Rampal et al. (2016) such that heat loss from the ocean that would cause super cooling is redirected to ice formation. The newly formed ice is transferred to the thin ice category and is assumed to have a thickness equal to h_{min}. The transfer from the thin ice to the thick ice and the lateral melting of thin ice are computed by applying the bounding limit h_{max} – if h_{t}>h_{max}, then we update the variables as follows:
$$\begin{array}{}\text{(A27)}& {\displaystyle}{h}_{{\mathrm{t}}^{\prime}}& {\displaystyle}={h}_{\mathrm{max}},\text{(A28)}& {\displaystyle}{A}_{{\mathrm{t}}^{\prime}}& {\displaystyle}={\displaystyle \frac{{h}_{\mathrm{max}}{h}_{\mathrm{min}}}{{h}_{\mathrm{t}}{h}_{\mathrm{min}}}}{A}_{\mathrm{t}},\text{(A29)}& {\displaystyle}{H}_{{\mathrm{t}}^{\prime}}& {\displaystyle}={A}_{{\mathrm{t}}^{\prime}}*{h}_{{\mathrm{t}}^{\prime}},\text{(A30)}& {\displaystyle}{h}_{\mathrm{s},\mathrm{t}}^{\prime}& {\displaystyle}={\displaystyle \frac{{A}_{{\mathrm{t}}^{\prime}}}{{A}_{\mathrm{t}}}}{h}_{\mathrm{s},\mathrm{t}},\text{(A31)}& {\displaystyle}{H}^{\prime}& {\displaystyle}=H({H}_{{\mathrm{t}}^{\prime}}{H}_{\mathrm{t}}),\text{(A32)}& {\displaystyle}{A}^{\prime}& {\displaystyle}=A({A}_{{\mathrm{t}}^{\prime}}{A}_{\mathrm{t}}),\text{(A33)}& {\displaystyle}{h}_{{\mathrm{s}}^{\prime}}& {\displaystyle}={h}_{\mathrm{s}}({h}_{\mathrm{s},\mathrm{t}}^{\prime}{h}_{\mathrm{s},\mathrm{t}}).\end{array}$$The form of the reduction in thin ice concentration in (Eq. A28) is a little arbitrary, but we wanted to allow the possibility of the thin ice completely changing into thick ice at some point.
This work is the result of a longterm team effort at the Nansen Centre in Norway carried out by, in alphabetical order, SB, VD, EO, PR, AS and TW to develop the new sea ice model neXtSIM. PR led the research and performed the analyses presented in this paper; VD and PR led the writing with contributions from EO, SB, TW and AK. AS implemented the parallel C$++$ version of the model code described in Samaké et al. (2017), with support from SB and EO. VD implemented the MEB rheology in the version of the neXtSIM code used in this study.
The authors declare that they have no conflict of interest.
This paper was prepared thanks to the financial support of the Research Council of Norway (RCN) through the FRASIL project (grant no. 263044). However, the development of the neXtSIM model has been supported since 2013 through several projects from the RCN, in particular the SIMECH (grant no. 231179) and NEXTWIM (grant no. 244001) projects. TOTAL E&P are also thanked for their continuous support over the period 2013–2017 through the KARA project. And finally, we thank Jerome Weiss and David Marsan for fruitful discussions.
This paper was edited by Christian Haas and reviewed by two anonymous referees.
This research has been supported by the Research Council of Norway (grant nos. 263044, 231179 and 244001).
Amante, C. and Eakins, B. W.: ETOPO1 Global Relief Model converted to PanMap layer format, Pangaea, https://doi.org/10.1594/PANGAEA.769615, 2009. a
Amitrano, D., Grasso, J. R., and Hantz, D.: From diffuse to localised damage through elastic interaction, Geophys. Res. Lett., 26, 2109–2112, 1999. a
Armitage, T. W. K., Bacon, S., Ridout, A. L., Thomas, S. F., Aksenov, Y., and Wingham, D. J.: Arctic sea surface height variability and change from satellite radar altimetry and GRACE, 2003–2014, J. Geophys. Res.Oceans, 121, 4303–4322, https://doi.org/10.1002/2015JC011579, 2016. a
Armitage, T. W. K., Bacon, S., Ridout, A. L., Petty, A. A., Wolbach, S., and Tsamados, M.: Arctic Ocean surface geostrophic circulation 2003–2014, The Cryosphere, 11, 1767–1780, https://doi.org/10.5194/tc1117672017, 2017. a
Bouchat, A. and Tremblay, B. L.: Using seaice deformation fields to constrain the mechanical strength parameters of geophysical sea ice, J. Geophys. Res.Oceans, 122, 5802–5825, 2017. a, b
Bouillon, S. and Rampal, P.: On producing sea ice deformation data sets from SARderived sea ice motion, The Cryosphere, 9, 663–673, https://doi.org/10.5194/tc96632015, 2015a. a, b, c, d, e, f, g
Bouillon, S. and Rampal, P.: Presentation of the dynamical core of neXtSIM, a new sea ice model, Ocean Modell., 91, 23–37, 2015b. a, b, c, d, e, f
Bromwich, D. H., Wilson, A. B., Bai, L.S., Moore, G. W. K., and Bauer, P.: A comparison of the regional Arctic System Reanalysis and the global ERAInterim Reanalysis for the Arctic, Q. J. Roy. Meteorol. Soc., 142, 644–658, https://doi.org/10.1002/qj.2527, 2016. a
Coon, M., Kwok, R., Levy, G., Pruis, M., Schreyer, H., and Sulsky, D.: Arctic Ice Dynamics Joint Experiment (AIDJEX) assumptions revisited and found inadequate, J. Geophys. Res., 112, C11S90, https://doi.org/10.1029/2005JC003393, 2007. a
Dansereau, V., Weiss, J., Saramito, P., and Lattes, P.: A Maxwell elastobrittle rheology for sea ice modelling, The Cryosphere, 10, 1339–1359, https://doi.org/10.5194/tc1013392016, 2016. a, b, c, d, e, f, g
Dansereau, V., Weiss, J., Saramito, P., Lattes, P., and Coche, E.: Ice bridges and ridges in the MaxwellEB sea ice rheology, The Cryosphere, 11, 2033–2058, https://doi.org/10.5194/tc1120332017, 2017. a
Erlingsson, B.: Twodimensional deformation patterns in sea ice, J. Glaciol., 34, 118, 1988. a
Esau, I. N.: Amplification of turbulent exchange over wide Arctic leads: Largeeddy simulation study, J. Geophys. Res.Atmos., 112, d08109, https://doi.org/10.1029/2006JD007225, 2007. a
Girard, L., Bouillon, S., Weiss, J., Amitrano, D., Fichefet, T., and Legat, V.: A new modelling framework for sea ice mechanics based on elastobrittle rheology, Ann. Glaciol., 52, 123–132, 2011. a
Gray, J. and Morland, L.: A twodimensional model for the dynamics of sea ice, Philos. Trans. Roy. Soc. Londo A, 347, 219–290, 1994. a
Hibler, W. D.: A dynamic thermodynamic sea ice model, J. Phys. Ocean., 9, 817–846, 1979. a
Hutchings, J. K., Heil, P., and Hibler, W. D.: Modeling linear kinematic features in sea ice, Mon. Weather Rev., 133, 3481–3497, 2005. a
Hutchings, J. K., Roberts, A., Geiger, C. A., and RichterMenge, J.: Spatial and temporal characterization of seaice deformation, Ann. Glaciol., 52, 360–368, 2011. a, b, c
Hutter, N. and Losch, M.: Featurebased comparison of seaice deformation in leadresolving seaice simulations, The Cryosphere Discuss., https://doi.org/10.5194/tc201988, in review, 2019. a
Hutter, N., Losch, M., and Menemenlis, D.: Scaling Properties of Arctic Sea Ice Deformation in a HighResolution ViscousPlastic Sea Ice Model and in Satellite Observations, J. Geophys. Res.Oceans, 117, 1038–1016, 2018. a, b
Hutter, N., Zampieri, L., and Losch, M.: Leads and ridges in Arctic sea ice from RGPS data and a new tracking algorithm, The Cryosphere, 13, 627–645, https://doi.org/10.5194/tc136272019, 2019. a, b
Kolmogorov, A. N.: A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number, J. Fluid Mechan., 13, 82–85, 1962. a
Kozo, T. L.: Initial model results for Arctic mixed layer circulation under a refreezing lead, J. Geophys. Res.Oceans, 88, 2926–2934, https://doi.org/10.1029/JC088iC05p02926, 1983. a
Kwok, R.: Deformation of the Arctic Ocean sea ice cover between November 1996 and April 1997: a qualitative survey, Solid Mechan. Appl., 94, 315–322, 2001. a, b
Kwok, R., Schweiger, A., Rothrock, D., Pang, S. S., and Kottmeier, C.: Sea ice motion from satellite passive microwave imagery assessed with ERS SAR and buoy motions, J. Geophys. Res., 103, 8191–8214, 1998. a, b
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
Kwok, R., Cunningham, G. F., Wensnahan, M., Rigor, I., Zwally, H. J., and Yi, D.: Thinning and volume loss of the Arctic Ocean sea ice cover: 2003–2008, J. Geophys. Res.Oceans, 114, c07005, https://doi.org/10.1029/2009JC005312, 2009. a
Lindsay, R. W. and Stern, H. L.: The RADARSAT geophysical processor system: Quality of sea ice trajectory and deformation estimates, J. Atmos. Ocean. Technol., 20, 1333–1347, 2003. a, b
Linow, S. and Dierking, W.: ObjectBased Detection of Linear Kinematic Features in Sea Ice, Remote Sens., 9, 5, https://doi.org/10.3390/rs9050493, 2017. a
Lovejoy, S. and Schertzer, D.: Scale, Scaling and Multifractals in Geophysics: Twenty Years on, in: Nonlinear Dynamics in Geosciences, 311–337, Springer, New York, NY, 2007. a
Marcq, S. and Weiss, J.: Influence of sea ice leadwidth distribution on turbulent heat transfer between the ocean and the atmosphere, The Cryosphere, 6, 143–156, https://doi.org/10.5194/tc61432012, 2012. a
Marsan, D. and Weiss, J.: Space/time coupling in brittle deformation at geophysical scales, Earth Planet. Sci. Lett., 296, 353–359, 2010. a, b
Marsan, D., Stern, H. L., 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, f, g, h
Matsushita, M.: Fractal Viewpoint of Fracture and Accretion, J. Phys. Soc. JPN, 54, 857–860, https://doi.org/10.1143/JPSJ.54.857, 1985. a
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, d, e, f, g, h, i
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, b, c, d, e, f, g, h, i
Rampal, P., Weiss, J., and Marsan, D.: Positive trend in the mean speed and deformation rate of Arctic sea ice, 1979–2007, J. Geophys. Res.Oceans, 114, c05013, https://doi.org/10.1029/2008JC005066, 2009a. a
Rampal, P., Weiss, J., Marsan, D., and Bourgoin, M.: Arctic sea ice velocity field: General circulation and turbulentlike fluctuations, J. Geophys. Res.Oceans, 114, C10014, https://doi.org/10.1029/2008JC005227, 2009b. a
Rampal, P., Bouillon, S., Ólason, E., and Morlighem, M.: neXtSIM: a new Lagrangian sea ice model, The Cryosphere, 10, 1055–1073, https://doi.org/10.5194/tc1010552016, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m
Remacle, J.F. and Lambrechts, J.: Fast and Robust Mesh Generation on the Sphere? Application to Coastal Domains, Proc. Eng., 163, 20–32, https://doi.org/10.1016/j.proeng.2016.11.011, 25th International Meshing Roundtable, 2016. a
Rothrock, D. A. and Thorndike, A. S.: Geometric properties of the underside of sea ice, J. Geophys. Res., 85, 3955–3963, 1980. a
Rothrock, D. A. and Thorndike, A. S.: Measuring the sea ice floe size distribution, J. Geophys. Res.Oceans, 89, 6477–6486, https://doi.org/10.1029/JC089iC04p06477, 1984. a
Sakov, P., Counillon, F., Bertino, L., Lisæter, K. A., Oke, P., and Korablev, A.: TOPAZ4: An ocean sea ice data assimilation system for the North Atlantic and Arctic, Ocean Sci., 8, 633–662, https://doi.org/10.5194/osd915192012, 2012. a
Samaké, A., Rampal, P., Bouillon, S., and Ólason, E.: Parallel implementation of a Lagrangianbased model on an adaptive mesh in C++: Application to seaice, J. Comput. Phys., 350, 84–96, 2017. a, b
Semtner, A. J.: A model for the thermodynamic growth of sea ice in numerical investigations of climate, J. Phys. Ocean., 6, 379–389, 1976. a
Smith, J.: Oceanographic investigations during the AIDJEX lead experiment, AIDJEX Bull., 27, 125–133, 1974. a
Sornette, D.: Power Law Distributions, 93–121, Springer Berlin Heidelberg, Berlin, Heidelberg, https://doi.org/10.1007/3540331824_4, 2006. a
Spreen, G., Kwok, R., Menemenlis, D., and Nguyen, A. T.: Seaice deformation in a coupled oceanseaice model and in satellite remote sensing data, The Cryosphere, 11, 1553–1573, https://doi.org/10.5194/tc1115532017, 2017. a
Stern, H. L. and Lindsay, R. W.: Spatial scaling of Arctic sea ice deformation, J. Geophys. Res.Oceans, 114, c10017, https://doi.org/10.1029/2009JC005380, 2009. a
Stern, H. L. and Moritz, R. E.: Sea ice kinematics and surface properties from RADARSAT synthetic aperture radar during the SHEBA drift, J. Geophys. Res.Oceans, 107, SHE14–1–SHE14–10, https://doi.org/10.1029/2000JC000472, 2002. a
Stern, H. L. and Rothrock, D. A.: Open water production in Arctic sea ice: Satellite measurements and model parameterizations, J. Geophys. Res., 100, 20601–20612, 1995. a, b
Stern, H. L., Schweiger, A. J., Zhang, J., and Steele, M.: On reconciling disparate studies of the seaice floe size distribution, Elementa Sci. Anthro., 6, 1–16, 2018. a
Thomas, M., Kambhamettu, C., and Geiger, C. A.: Discontinuous NonRigid Motion Analysis of Sea Ice using CBand Synthetic Aperture Radar Satellite Imagery, in: Computer Vision and Pattern Recognition Workshop, 2004, CVPRW '04, Conference on, p. 24, 2004. a
Thomas, M., Geiger, C. A., Kambhamettu, C., Hutchings, J., RichterMenge, J. A., and Engram, M. J.: Nearreal time application of sarderived sea ice Differential motion during aplis ice camp, 2007. a, b, c
Thomas, M., Kambhamettu, C., and Geiger, C. A.: Mapping of large magnitude discontinuous sea ice motion, SIGSPATIAL Special, 1, 45–50, 2009. a
Tonboe, R. T., Eastwood, S., Lavergne, T., Sørensen, A. M., Rathmann, N., Dybkjær, G., Pedersen, L. T., Høyer, J. L., and Kern, S.: The EUMETSAT sea ice concentration climate data record, The Cryosphere, 10, 2275–2290, https://doi.org/10.5194/tc1022752016, 2016. a
Wang, K.: Observing the yield curve of compacted pack ice, J. Geophys. Res.Oceans, 112, c05015, https://doi.org/10.1029/2006JC003610, 2007. 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
Weiss, J.: Exploring the “solid turbulence” of sea ice dynamics down to unprecedented small scales, J. Geophys. Res.Oceans, 122, 6071–6075, 2017. a
Weiss, J. and Dansereau, V.: Linking scales in sea ice mechanics, Philos. Trans. Roy. Soc. A, 375, 20150352, https://doi.org/10.1098/rsta.2015.0352, 2017. a, b, c, d, e, f, g
Weiss, J. and Marsan, D.: Scale properties of sea ice deformation and fracturing, Comp. Rem. Phys., 5, 735–751, 2004. a, b
Weiss, J., Marsan, D., and Rampal, P.: Space and Time Scaling Laws Induced by the Multiscale Fracturing of The Arctic Sea Ice Cover, in: IUTAM Symposium on Scaling in Solid Mechanics, edited by Borodich, F., Vol. 10 of Iutam Bookseries, 101–109, Springer Netherlands, https://doi.org/10.1007/9781402090332_10, 2009. a
Williams, T. D., Rampal, P., and Bouillon, S.: Waveice interactions in the neXtSIM seaice model, The Cryosphere, 11, 2117–2135, https://doi.org/10.5194/tc1121172017, 2017. a
Winton, M.: A reformulated threelayer sea ice model, J. Atmos. Ocean. Technol., 17, 525–531, 2000. a
https://rda.ucar.edu/datasets/ds631.0 (last access: 15 April 2015), ASRv1 30km, formerly ASR final version, Byrd Polar Research Centre and The Ohio State University.