## Abstract

Mechanics plays a key role in the evolution of the sea ice cover through its control on drift, on momentum and thermal energy exchanges between the polar oceans and the atmosphere along cracks and faults, and on ice thickness distribution through opening and ridging processes. At the local scale, a significant variability of the mechanical strength is associated with the microstructural heterogeneity of saline ice, however characterized by a small correlation length, below the ice thickness scale. Conversely, the sea ice mechanical fields (velocity, strain and stress) are characterized by long-ranged (more than 1000 km) and long-lasting (approx. few months) correlations. The associated space and time scaling laws are the signature of the brittle character of sea ice mechanics, with deformation resulting from a multi-scale accumulation of episodic fracturing and faulting events. To translate the short-range-correlated disorder on strength into long-range-correlated mechanical fields, several key ingredients are identified: long-ranged elastic interactions, slow driving conditions, a slow viscous-like relaxation of elastic stresses and a restoring/healing mechanism. These ingredients constrained the development of a new continuum mechanics modelling framework for the sea ice cover, called Maxwell–elasto-brittle. Idealized simulations without advection demonstrate that this rheological framework reproduces the main characteristics of sea ice mechanics, including anisotropy, spatial localization and intermittency, as well as the associated scaling laws.

This article is part of the themed issue ‘Microdynamics of ice’.

## 1. Introduction

Sea ice mechanics, and its consequences in terms of sea ice drift, energy and momentum exchanges between the ocean and atmosphere, or sea ice mass balance, involve a huge range of spatial and time scales, from the thickness scale (approximately in metres) to the scale of the Arctic basin (approx. 10^{6} m), and from seconds to years. Fracturing and crushing of saline ice against offshore structures or ship hulls involve even smaller scales, down to internal microstructural defects, such as grain boundaries or brine pockets (e.g. [1]). A clear understanding of these mechanisms is, of course, necessary to estimate ice loads on these structures [2]. When averaging over much larger scales, say at the scale of the Arctic basin and over decades, a mechanical weakening of the Arctic sea ice cover has been recently evidenced [3], which probably explains the acceleration of both ice drift and deformation over the last decades [4]. This positive trend on sea ice kinematics is not reproduced in global circulation models, thus impacting their estimation of the Arctic sea ice mass balance and possibly explaining why the ice cover declines faster than foreseen [5]. At intermediate spatial scales, i.e. from local (approx. 1 km) to regional (approx. 100 km) scales, and time scales between hours and several days, an accurate modelling of sea ice drift and mechanics is the key milestone towards the development of sea ice operational forecasting platforms [6]. The importance of sea ice mechanics at all scales, though in various contexts (such as climate modelling, regional forecasting and offshore engineering), raises the fundamental question of the linkage of scales. This question is not new [2,7–9], but has been brought again into light during the last decades thanks to the analysis of *in situ* and remote sensing data of the Arctic sea ice cover [10].

This issue is classic in materials science. For the analysis of heterogeneous media, it is generally handled through homogenization procedures (e.g. [11]). In this case, a representative volume element (RVE) can be defined. The mechanical behaviour of the RVE is representative of the behaviour of the medium at any larger scale, hence allowing upscaling. Conversely, the analysis of the interactions between a local entity/heterogeneity (e.g. a grain in a polycrystal and fibres in a composite material) and its environment can be performed by embedding the entity into a representative, averaged medium. Homogenization is, therefore, at the root of multi-scale material modelling [12]. However, two fundamental conditions are required to define an RVE. First, the correlation length associated with the heterogeneities, or with the mechanical fields (such as stress and strain), should be much smaller than the system size, as the RVE scale has to be necessarily larger than this correlation length [13]. Second, the local fluctuations of the mechanical fields or properties should not be fat-tailed, that is, their probability density functions (PDFs) should decay towards large *x*-values faster than *x*^{−3}. This condition ensures that the PDF, *P*(*x*), remains in the Gaussian basin of attraction, hence that the central limit theorem applies [14]. Instead, in the case of a fat-tailed *P*(*x*) ∼ *x*^{−}* ^{κ}* with exponent

*κ*< 3, the moments of the distribution of order equal to or larger than 2, and so the standard deviation, are ill-defined, whereas for

*κ*< 2, the mean itself is ill-defined. The violation of at least one of these two conditions precludes the convergence upon the addition of random variables, meaning that distributions remain fat-tailed whatever the scale of coarse-graining, thus invalidating homogenization at some intermediate scale.

The problem of scale linkage is also crucial in geophysics, which often covers a huge range of space and time scales. However, several geophysical processes and systems, such as atmospheric turbulence [15] or crustal deformation and faulting (e.g. [16]), are characterized by scaling laws that prevent an artificial disconnection between different scales. In this paper, after a brief discussion on small-scale heterogeneities in saline ice, and their role in the variability of mechanical properties (in particular, strength), we will show that sea ice is a paradigm for such scale-invariant systems, characterized by coupled space and time scaling laws for the mechanical fields spanning (at least) hours to months on the one hand, and a few kilometres to the Arctic basin scale on the other. We will then discuss the physical origin of these scaling laws and how they can be used to link scales in sea ice mechanics, in a way that is fundamentally different from homogenization procedures. Finally, we will discuss how these laws constrain developments in the field of sea ice modelling, and present a mechanical framework built in light of these constraints and their underlying physics.

## 2. Sea ice microstructure and the small-scale variability of mechanical properties

Our goal here is not to give a comprehensive review of sea ice microstructural features, or of their impact on sea ice mechanical properties, as this can be found elsewhere [17–19]. It is sufficient to say that sea ice is a multi-phase material composed of crystalline ice, gas, brines and solid salts [19], with a complex microstructure characterized by tortuous grain boundaries, brine pockets and channels, as well as columnar or more isotropic grains [17]. At the scale of the sea ice cover thickness, stratification can be important, especially in multi-year ice, inducing heterogeneity. These structural heterogeneities have an impact on mechanical properties, such as elastic moduli or strength. For instance, at the scale of the polycrystal, brine pockets can play the role of stress concentrators, and so of microcrack initiators, leading to a decreasing average flexural strength with increasing porosity [19]. At the scale of the ice thickness, the columnar structure and stratification lead to significant strength anisotropy [18,20].

However, instead of focusing on the role of the microstructure on *mean* mechanical properties, here we briefly discuss its impact on the spatial as well as temporal (as sea ice is known to age [19]) *variability*. If the literature on mean mechanical properties from laboratory experiments is vast, datasets allowing a sound estimation of associated variability are scarce. Kamio *et al*. [21] reported a statistical analysis of tensile and bending strengths of first-year fast sea ice cored in the Okhotsk Sea. They fitted the data with a Weibull's distribution, suggesting a weakest-link mechanism, plausible for tensile fracturing [22], and linked this variability of the mechanical strength to that of the microstructure. However, to ascertain the relevance of the weakest-link hypothesis, a comparison of datasets obtained for different sample sizes is necessary [23]. As a matter of fact, the strength distributions reported in [21] are close to Gaussian distributions, with standard deviations of the order of the means. It indicates that very low local strengths (less than 100 kPa) are not unlikely. Note that a significant variability of the compressive strength has also been observed for fresh-water artificial granular (isotropic) ice [24]. In this case, the strength variability is normally distributed as well [23], though with a smaller standard deviation/mean ratio compared with the tensile or flexural strength of sea ice [21], as its only possible microstructural origin is the polycrystalline structure.

Therefore, microstructural complexity of sea ice gives rise to a significant variability in *local* mechanical strength. It is difficult however, from the existing literature, to estimate a correlation length *ξ*_{m} associated with either these microstructural characteristics or strength properties. Existing data [21] from samples taken from the same ice plate but exhibiting a large variability of strength suggest a small *ξ*_{m}, that is below the ice thickness *h*. At much larger (regional) scales, one may expect some very long wavelength variability of sea ice properties, such as brine content and porosity, as the result of large-scale environmental conditions. However, at scales larger than a few ice thicknesses, it seems reasonable to assume the relevance of homogenization procedures to define the average characteristics of sea ice, such as salinity, porosity or elastic moduli. If so, does this mean a smooth mechanical response of the sea ice cover? As shown below, the answer is ‘no’, as the small-scale variability (or disorder) of strength, coupled with long-ranged elastic interactions and brittle mechanics, gives rise to scale invariance.

## 3. Scaling of sea ice mechanics and deformation

Historically, the first evidence of scaling in sea ice mechanics came from the analysis of sea ice floe size distributions [25,26]. Qualitatively, this spatial scale invariance is expressed by the impossibility of estimating the resolution of an aerial or satellite image of the sea ice cover without a scale of reference (e.g. an island or a ship). Quantitatively, a power-law tail for the PDF of the floe sizes *s*, *P*(*s*) ∼ *s*^{−}* ^{η}*, with the exponent

*η*in the range 2.0–3.0, is generally observed (e.g. [10,27] for a review). When combining datasets, such scaling has been shown to hold from the scale of the ice thickness up to at least 50 km, meaning that no characteristic scale can be defined over this range. Simple hierarchical models can be proposed to explain these observations conceptually, all of which are based on scale-invariant breaking rules [28,29]. Combined with more detailed analyses of fracturing patterns [30], this argues for an underlying scale-independent fracturing and fragmentation process for pack ice, from the thickness of the ice cover to regional scales, and probably beyond. This scale invariance can break down near the ice edge, where the action of the ocean swell, associated with a well-defined wavelength, sets an upper bound to scaling as the result of bending failure of ice floes [31,32].

Away from this influence of waves, the observed scale invariance implies that small scales cannot be arbitrarily mechanically disconnected from large ones, and that the definition of an RVE is elusive. However, these observations are episodic, and say little about sea ice drift, deformation and time patterns. Hence, they can hardly be used to constrain continuum mechanics sea ice models. Two datasets, Lagrangian trajectories of ice-tethered buoys, and satellite-derived velocities [33,34], available since the early 1980s, allowed a much more comprehensive view of sea ice motion and deformation over most of the Arctic (Antarctic datasets are more scarce [35]), from spatial scales ranging from approximately a few kilometres to approximately 1000 km and time scales from hours (buoys) to years. After the pioneering work of Thorndike and co-workers inspired from the analysis of turbulent Lagrangian trajectories [36,37], the last decade has seen a renewing interest in the scaling analysis of these datasets. Other datasets, such as internal stress time series [38], exhibit an intermittent behaviour with scale invariance in the time domain, up to a time scale of a few months (100 days) [39]. Although fully consistent with the intermittency of sea ice deformation, these data obtained from local sensors do not allow an analysis in the spatial domain. Hence, here we focus on sea ice deformation, for which space and time scaling relationships can be analysed concomitantly. Section 4 will discuss the implication of these scaling properties in terms of the underlying physical and mechanical processes, whereas §5 presents how they can be used to constrain the development of a new sea ice mechanical framework.

### (a) Sea ice deformation

Owing to the huge aspect ratio, lateral extension/thickness of the sea ice cover, non-horizontal mechanics can be safely neglected down to at least the kilometre scale. To fully calculate a two-dimensional strain-rate tensor and the associated invariants, such as shear rate, , or divergence rate, , over a given domain, at least three sea ice trajectories have to be recorded simultaneously. This strongly limits the use of buoy trajectories to analyse the scaling of sea ice deformation, except in the case of dense arrays [35]. However, a strain-rate proxy , excluding solid rotation but unable to differentiate shear from divergence, can be defined from the dispersion rate of pairs of buoys, making a statistically significant analysis of scaling possible [40]. Satellite synthetic aperture radar (SAR) imagery allowed a much denser sampling of the sea ice velocity field, and mapping of the strain-rate tensor over a large part of the Arctic at a spatial resolution of 10 km, over several years [34]. This revealed that most of the sea ice deformation is accommodated along linear-like structures [41], concentrating both shear and divergence/convergence [42]. These structures can be considered as active Coulombic shear faults [30,43], giving a first indication of the brittle character of sea ice mechanics. The main limitation of this satellite-derived dataset is its time resolution, of the order of 3 days, which prevents the exploration of sea ice kinematics at small time scales.

If Thorndike and co-workers were the first, in the 1980s, to give a mechanical interpretation of the dispersion of sea ice trajectories [37], a more comprehensive analysis of the scaling properties of sea ice deformation was performed only over the last 12 years [44] from both buoy data and SAR imagery. The main results are summarized below.

At the scales of the satellite-data resolution (10 km, 3 days), the PDFs of shear rates exhibit power-law tails, , with *κ* = 2.4 ± 0.1 < 3 [42,44]. As noted above, this precludes the convergence upon the addition of random variables, and thus represents a first argument against homogenization in sea ice mechanics. In the case of divergence rates, an ambiguity remains due to the presence of artificial noise in RADARSAT Geophysical Processor System (RGPS) data [45]. The power-law exponent *κ* increases with increasing space and/or time scale, constituting a first indication of strong spatial and temporal correlations of these mechanical fields [10]. However, a direct scaling analysis of these PDFs is complex. To explore these aspects, multi-fractal analyses of the corresponding moments were, therefore, privileged [44]. These revealed both space and time scaling laws characterizing, respectively, the localization and intermittency of sea ice deformation.

In the spatial domain, the moments of the distributions of shear (), absolute divergence () and total () strain rates, measured at the RGPS time resolution (3 days), scale as
3.1
over scales *L* ranging from the RGPS resolution (10 km) to the scale of the Arctic basin as a whole (approx. 1000 km) [44–46]. The structure function *β*(*q*) is quadratic, *β*(*q*) = *aq*^{2} + *bq*, and convex (*a* > 0, *b* > 0). This spatial multi-fractality of the ice cover quantitatively characterizes its scale-dependent localization. Indeed, the convexity of *β*(*q*) implies that the spatial heterogeneity of sea ice deformation increases towards small scales, in full agreement with the scale dependence of the power-law exponent *κ*. As *β*(1) > 0, the mean strain-rate invariants are scale-dependent themselves, i.e. sea ice deformation is non-conservative. It has been shown that *β*(1) depends on both the season (i.e. it is larger in summer (approx. 0.25) than in winter (approx. 0.15–0.2)) and the age of the sea ice (larger for first-year than for multi-year ice) [47]. This is probably reminiscent of a different behaviour between a highly cohesive sea ice pack and a more fragmented, granular-like cover [10]. However, the main consequence of the observed scaling (3.1) is the presence of a correlation length *ξ*_{M} of sea ice deformation of the order of the Arctic basin scale. This raises the fundamental question of the underlying mechanisms that allow a very small correlation length *ξ*_{m} for the material's strength (see §2) to be translated into a more than 1000 km correlation length *ξ*_{M} for the mechanical fields, and represents a second impediment to homogenization procedures.

Figure 1 shows that the spatial multi-fractal scaling (3.1) is mirrored in the time domain, over time scales ranging in winter from the RGPS time resolution (3 days) to a season (approx. 150 days), such that
3.2
with strain rates estimated at a spatial scale of 10 km. As for spatial scaling, the structure function *α*(*q*) is quadratic and convex. In a way similar to that discussed above, relation (3.2) thus characterizes the strong intermittency of sea ice deformation, i.e. the time clustering of the transient deformation episodes. The associated correlation time *T*_{M} is of the order of 5 months in winter, meaning that sea ice mechanics loses its memory of previous deformation history very slowly, in agreement with the observation that most of the deformation is accommodated along long-lasting active faults [48]. This signs the link between sea ice deformation and the state of fracturing (or ‘damage’, see §5). Available SAR-derived data make an analysis over longer time scales difficult. In particular, checking whether this long-term memory survives the summer season in multi-year ice, in particular, is tricky, if not impossible.

Buoy-derived analyses [35,40,49] are fully consistent with these results. Using a strain-rate proxy defined from buoy dispersion, Rampal *et al*. [40], though computing only the first-order moments , went one step further in analysing the coupling between space and time scaling laws. They showed that the spatial scaling is time-dependent, whereas the temporal scaling is spatial-scale-dependent, with exponent *β*(1) (respectively, *α*(1)) decreasing as a function of *t* (respectively, *L*):
3.3

Relation (3.3*b*) is consistent with the power spectral analysis of sea ice deformation of Hutchings *et al*. [49], who reported a power-law behaviour with *f* = 1/*t* and *θ* an increasing function of *L*, as theoretical considerations imply *θ*(*L*) = 1 − 2*α*(1,*L*) [10]. As *α*(1) can be considered as a measure of intermittency, varying between 0 (non-intermittent ‘flow’) and 1 (extreme intermittency), relation (3.3*b*) expresses a more continuous sea ice mechanics with increasing spatial scale. Similarly, the measure of spatial localization, *β*(1), is bounded between 0 (a perfectly homogeneous field) and 2 (a point-like process), and equation (3.3*a*) reflects a progressive smoothing as one cumulates deformation events decreasingly correlated with increasing time. The unique expression compatible with (3.3*a*,*b*) is the following coupled equation [50]:
3.4
which summarizes the complex space/time coupling of sea ice deformation, with *c* = 0.10 a coupling constant. Interestingly, a similar scaling symmetry is present in the brittle deformation of the Earth's crust, arguing for similar deformation mechanisms [50].

As already stressed above, the scaling properties of sea ice deformation rule out the use of spatial homogenization procedures to link scales in sea ice mechanics. Long-term correlations (at least 5 months in winter) raise similar problems in the time domain. However, these scaling laws provide tools for upscaling or downscaling in the statistical sense. Indeed, knowing the strain-rate PDF at some given scales (e.g. at the SAR resolution), the multi-fractal properties of the associated field (equation (3.1) or (3.2)) can be used to extrapolate the PDF to either larger (upscaling) or smaller (downscaling) scales [44]. This can be done in either the space or the time domain, separately, as a coupled space/time multi-fractal analysis has not yet been performed. Assuming that the scaling laws extrapolate down to scales smaller than the resolution of satellite imagery (3 days, 10 km), the PDF can be downscaled and compared to the strain rate defining the transition from a ductile to a brittle behaviour for bulk (undamaged) saline ice [18,51]. In this way, it has been shown that, if considering spatial scales of the order of the sea ice thickness (approximately metres), all deformation occurs at rates larger than , i.e. in a brittle manner [10,44].

## 4. The physical origin of scaling

The analyses summarized above revealed (i) that sea ice mechanics is characterized by a correlation length of the order of the system as a whole and a correlation time at least of the order of a season, (ii) a multi-fractal scaling signing an extreme intermittency and spatial heterogeneity of deformation, and (iii) in conjunction with other sources [30,52], a brittle mechanics of the ice cover. Said differently, and raising again the analogy with Earth's crustal deformation, sea ice responds to the action of winds and currents through episodic fracturing/faulting events spanning a wide range of scales and whose complex mechanical interactions result in the observed scaling laws [50].

Many other physical systems, such as ferromagnets [53], plastically deformed crystals [54], individual cracks propagating in heterogeneous materials [55], or, as noted above, the brittle crust, respond to external forcing through such jerky dynamics. As all these systems share common scale-free statistics and scaling properties, universal underlying ingredients have been identified to give rise to this ‘crackling’ dynamics [56]. These cornerstones are summarized below in the context of sea ice mechanics, and will be used in §5 to constrain the development of physically sound sea ice mechanical models.

(1) A

*threshold mechanics*should be at work at the elementary scale, meaning that damage, fracturing, slip and/or deformation is triggered when and where the driving ‘force’ (e.g. stress) overcomes a ‘strength’. In sea ice, faulting is triggered once the Coulomb stress reaches a cohesion threshold*τ*_{0}, where*τ*is the shear stress resolved on the fault's plane,*σ*_{n}the normal stress and*μ*an internal friction coefficient [30,43,52].(2)

*Disorder*should be associated with this threshold. If the friction coefficient weakly depends on temperature [57], the value*μ*= 0.7–0.8 is generally considered as being scale-independent and rather universal [43,58]. We thus expect disorder to be mainly associated with the cohesion*τ*_{0}, in relation to microstructural effects and expressing the variability of mechanical strength discussed in §2. As noted above, most of this variability is short-range-correlated and thus cannot explain by itself the emergence of scaling.(3)

*Long-ranged interactions*are the key ingredient to translate this short-range-correlated disorder into long-ranged space and time correlations of the mechanical fields.*Elastic interactions*can play this role in the spatial domain [55]: in an elastic medium, once the strength is locally overcome and the material fractures, damages and/or deforms inelastically, elastic stresses redistribute in 1/*r*, where^{δ}*r*is the distance and*δ*> 0 depends on the topology of the problem and the geometry of the failed ‘inclusion’ [59]. This long-ranged stress redistribution can trigger failure in the neighbourhood of the initial event, and possibly avalanches of failure throughout the entire system from a cascading mechanism.(4) For such a cascading effect to take place, a

*slow driving condition*must be fulfilled [60], meaning that the external driving rate must be much smaller than the rate of elastic stress redistribution in the medium. This is naturally met in the case of the sea ice cover, as the synoptic time scale associated with wind forcing (of the order of week) is orders of magnitude larger than the time needed for an elastic wave to travel across the Arctic basin (of the order of 10^{3}s, taking a shear wave speed of*c*_{s }≈ 1000 m s^{−1}[61]).

The combination of these four ingredients gives rise to jerky dynamics, with the system reaching a critical point upon increasing the external loading. In the vicinity of this critical point, the stored elastic energy is released through power-law-distributed avalanches. The archetype of this process is the so-called depinning transition [62]: Consider an elastic manifold pushed by a spatially homogeneous force *F* against a field of obstacles with variable but short-range-correlated pinning strengths. For small *F*, after some local depinning events from the weakest obstacles, the manifold is rapidly stopped. As *F* increases towards a critical value *F*_{c}, a jerky dynamics develops, characterized by power-law-distributed depinning avalanches and a correlation length *ξ* scaling with the system size. For *F* ≥ *F*_{c}, the manifold moves continuously. Far above *F*_{c}, it flows as if there were no disorder, i.e. in a non-jerky way and with short-range interactions only [63].

A formal analogy has been proposed recently between this depinning critical transition and the damage and failure of brittle media [23]. In this case, the ‘manifold’ is the spatially varying field of damage and macroscopic failure occurs for a critical value of the applied stress . A significant difference remains, however, as the classical depinning transition does not induce localization. In addition, if this scenario can explain the emergence of jerkiness and scaling on approaching the critical stress [64,65], it is unable to reproduce several important features of the deformation of geophysical objects, such as the sea ice cover or the Earth's crust. In particular, long-range time correlations, expressed, for example, by aftershock series in crustal seismicity or by the scaling relation (3.2) for sea ice, are not explained. More generally, instead of evolving towards a critical *point*, these systems remain in a marginally stable state characterized by strong jerkiness but stable power-law statistics and scaling relations. For these features to emerge, additional ingredients are needed.

(5)

*Viscous-like relaxation*of elastic stresses should be considered [66]. This introduces a viscosity*η*, hence a relaxation time λ =*η*/*E*(where*E*is an elastic modulus) that sets the characteristic time scale of direct aftershock sequences. Combined with cascading effects, this can lead to long-ranged time correlations of the kind of equation (3.2) [50].(6) A

*restoring/healing mechanism*, characterized by a healing time*t*_{h}, must be present to allow the highly damaged zones to slowly recover their mechanical properties, hence to avoid irreversibility of the process. In the case of sea ice, refreezing of cracks and faults is an obvious mechanism.(7) These relaxation and healing processes introduce feedback mechanisms that allow the system to remain all the time in the vicinity of the critical state [55] and therefore to exhibit a marginally stable state mimicking self-organized criticality [60]. For this self-adjustment to operate, a

*slow relaxation and healing condition*must be fulfilled, i.e. these processes should be much slower than elastic stress redistributions. Field studies of the refreezing within sea ice ‘leads’ (faults) showed that the time for 1 m of ice to grow within a lead opening of 10 cm for an air temperature of −15°C is of the order of a few days. The slow healing condition is, therefore, observed.

The ingredients for the emergence of scaling in mechanical systems being relevant in the case of sea ice, they provide a compelling physical framework to interpret the observations summarized in §3, and can, therefore, constrain the development of physically sound sea ice rheological/mechanical models.

## 5. Modelling the scale-invariant brittle mechanics of sea ice within a continuum mechanics framework

When considering small (less than 10 km) scales and/or regions of sea ice concentration smaller than approximately 85%, such as the marginal ice zone, the discontinuous nature of an ice cover made of an assembly of individual floes can hardly be ignored [67]. However, when pack ice and/or larger scales are concerned, as, for example, in climate models, sea ice mechanical models are all built from a continuum mechanics framework [68]. Some works parametrized a continuous representation of the mechanical response of floe assemblies from homogenization procedures [69,70], but this approach requires the maximum floe size to be significantly smaller than the grid scale of the continuum model.

One of the earliest sea ice rheological models was the AIDJEX elastic–plastic rheology [71]. Sea ice is considered as a spatially homogeneous (in terms of rheological properties) isotropic elastic medium as long as the stress state remains inside a yield envelope, but flows plastically following a normal flow rule once it reaches the envelope. As this requires keeping track of the ice strain state indefinitely, thus raising challenges for Eulerian numerical schemes, this model was tested for short-term simulations only. To solve this, Hibler proposed a viscous–plastic (VP) rheology in which sea ice is considered as a nonlinear viscous compressible fluid, coupled to a plastic yield curve defined by a concentration-dependent yield strength *P* [72]. This VP rheology became the standard sea ice mechanical modelling tool for regional to large-scale simulations. Several key ingredients listed in §4 are absent in this VP rheology. If a threshold mechanics is introduced through a plastic yield criterion, disorder is not considered. Without long-ranged elastic interactions, this precludes the emergence of spatial scaling and extreme statistics [42]. The dependence of the strength on ice concentration induces, in an Eulerian scheme, a weakening mechanism through a possible advection of ice to neighbouring cells. However, this mechanism remains local, so is the redistribution of stresses. Similarly, a plastic rheology implies by construction a vanishing relaxation time once the stress state reaches the yield envelope. This implies that the condition (7) of §4 is not fulfilled, thus preventing time correlations and intermittency to emerge as well.

These fundamental limitations, as well as other drawbacks of this framework, such as a normal flow rule in contradiction with a Coulomb's faulting mechanism of sea ice [43,52], called for a change of paradigm in sea ice mechanical modelling. Inspired by the modelling of brittle failure of heterogeneous media [73], Girard *et al*. proposed an elasto-brittle (EB) model that assimilates sea ice to a heterogeneous, damageable, elastic medium [74]. In this framework, ice is purely elastic everywhere, though with a variable modulus *E*(*x*,*t*) that deteriorates locally (i.e. at the grid scale) each time the stress state reaches a Coulomb's failure envelope. This deterioration is represented through a damage variable 0 < *d* ≤ 1 such that *E*(*x*,*t*) = *E*_{0}*d*(*x*,*t*), with *E*_{0} Young's modulus of intact ice, i.e. *d* = 1 for intact ice. Disorder is introduced at the local scale on the cohesion *τ*_{0} of the Coulomb's criterion. In this framework, a local damage event (*d* decreasing locally) induces an elastic redistribution of stresses that can trigger additional damage of the material on much larger scales through a cascading process. This model thus combines the first three ingredients for scaling listed in §4. The first implementation of this EB rheology into short (3 days), no-advection, stand-alone pan-Arctic simulations forced by realistic wind forcing from reanalyses was performed with an infinitely slow driving condition (external driving is artificially frozen as long as a damage cascading process proceeds). In agreement with observations, deformation was observed to strongly localize along linear-like features, to be characterized by power-law PDFs, and to follow a scaling relation . Although ice is considered as being elastically isotropic at the grid scale, stress and strain anisotropy spontaneously emerge at larger scales. This results from the anisotropy of the elastic kernel describing the Coulomb stress redistribution after a local damage event. Therefore, anisotropy does not need to be introduced at the grid element scale in an ad hoc manner. Recently, this EB model was implemented in a Lagrangian finite-element scheme and coupled to a thermodynamic model which also provides a parametrization of damage healing [75]. Stand-alone simulations forced by atmospheric and oceanic reanalyses were performed over durations ranging from 10 days [6] to 1 year [75] and exhibited realistic velocity fields and a multi-fractality of strain rates (equation (3.1)). This pleads for the pertinence of this model.

However, one key ingredient listed in §4, namely a viscous-like relaxation, is lacking in this seminal EB framework. This raises conceptual as well as technical difficulties. (i) As ice in this model remains *strictly elastic*, it is not consistent with permanent, large deformations, since, by definition, an isolated piece of material retrieves its initial shape when unloaded, whatever its level of damage. (ii) If *permanent* deformations were nonetheless allowed in this framework by assuming the shape of the elastic material to remain ‘frozen’ when un- or re-loaded, the model would become essentially elasto-plastic. In this case, long-ranged time correlations intrinsic to the mechanical behaviour (i.e. not inherited from the forcing) would not take place. To solve these fundamental issues, we recently revised the seminal EB model by introducing a viscous-like relaxation term in the otherwise linear-elastic constitutive relationship, and coupled this viscous component to the progressive damage mechanism. The model, named Maxwell–elasto-brittle (Maxwell-EB) after the classical Maxwell rheology, is described in full detail elsewhere [76]. Here, we recall its main characteristics in the context of the above discussion. The constitutive relationship reads
5.1
where D*σ*/D*t* is the objective time derivative of the total Cauchy stress tensor, *E* Young's modulus and **K** an adimensional stiffness matrix depending on Poisson's ratio *ν*. In this framework, *η* is an effective viscosity expressing how the fractured/damaged material loses its capacity to retain the memory of recoverable (i.e. elastic) deformations. The associated relaxation time λ = *η*/*E* should be theoretically infinite for an undamaged, purely elastic material, and decay with increasing level of damage, allowing elastic stresses to relax and irreversible deformations to take place in damaged zones. In this way, the accommodation of most of sea ice deformation along Coulombic faults can be accurately described in a continuum mechanics framework. This implies the need to couple both *E* and *η* to the damage variable *d* as follows:
5.2
with *α* > 1 such that decreases with increasing damage (that is, decreasing *d*). A healing mechanism is also introduced, which counterbalances damaging over large time scales (approx. *t*_{h}). Instead of using an infinitely slow driving, an explicit description of damage propagation is introduced: elastic perturbations propagate at the shear wave speed *c*_{s}, thus defining a discretization time scale for the numerical model, *t*_{d} = Δ*x*/*c*_{s} where Δ*x* is the model grid cell scale. Taking *c*_{s} ≈ 1000 m s^{−1} and, for example, Δ*x* = 10 km, then *t*_{d} ≈ 10 s; hence *t*_{d} is orders of magnitude smaller than synoptic time scales or typical healing times. Considering finally a very large relaxation time for undamaged ice, the Maxwell-EB framework henceforth incorporates all the ingredients listed in §4 and required for the emergence of the scale-invariant jerky dynamics of sea ice.

The damage parameter *α* controls the rate at which the relaxation time λ decreases with increasing damage, hence the mechanical behaviour itself. It is the only truly ad hoc parameter of the Maxwell-EB rheology and perhaps the least physically constrained. However, the following upper and lower bounds can be identified:

— For

*α*= 0, the apparent viscosity*η*is independent of*d*and hence constant in both space and time. Coupled with a very large viscosity for the initially undamaged material, , this implies a vanishing viscous dissipation term in the constitutive equation (5.1) whatever*d*, i.e. all deformations are elastic. Besides its physical inconsistency (see above), this case violates the principle of elastic stress relaxation, our ingredient (5) listed in §4.— For

*α*= 1,*η*decreases with increasing*d*within the same proportions as the material's elastic modulus. Hence λ is constant in both space and time, implying that the rate of dissipation of the internal stress in permanent deformations does not depend on the level of damage.— For

*α*> 1,*E*,*η*and λ all decrease with increasing damage, consistent with a degradation in both the mechanical strength and the capacity of the material to retain the memory of elastic deformations.— In the limit of , the relaxation time drops to zero at the onset of damage and the local memory of elastic deformations is lost. In this case, the Maxwell-EB rheology becomes elasto-plastic, and the condition (7) is violated.

We present the results of a set of small-deformation simulations (i.e. in which advection is neglected) representing the uniaxial, strain-driven and constant-rate compression of a two-dimensional rectangular ice plate, for different values of *α*. In terms of the prescribed forcing and boundary conditions, this constitutes perhaps the simplest test case that allows (i) investigating the range of mechanical behaviours reproduced by the model and (ii) demonstrating that, over a certain range of *α*-values, both heterogeneity and intermittency emerge naturally, i.e. from the physics included in the Maxwell-EB framework, and hence do not need to be accounted for in an ad hoc manner. All simulations are started from an initially undamaged plate with a uniform elastic modulus, chosen to be representative of sea ice on geophysical scales, and a uniform viscosity *η*_{0} set such that the undamaged relaxation time is as large as possible to accurately represent the limit of strictly elastic deformations (while ensuring the convergence of the numerical scheme). To mimic the small-scale, short-range-correlated variability of the sea ice strength (see §2), the (quenched) field of cohesion *τ*_{0}, identical for all simulations, is drawn randomly over each grid element from a uniform distribution spanning *in situ* stress measurements in Arctic sea ice [52]. This simulation set-up is described briefly in the electronic supplementary material, and in full detail elsewhere [76] along with a spatio-temporal characterization of the model's behaviour in terms of *deformation* fields. Hence, here we follow an approach taken for fracture-type models that record the number of broken fibres, ruptured bonds, depinning events, etc., and analyse the mechanical behaviour in terms of the spatial distribution and time series of *failure,* i.e. *damage* events. Animations of the evolution of the damage rate (defined as the number of damaged model elements per time step multiplied by their distance to the damage criterion [76]) and of the total strain-rate for the different simulations are available as the electronic supplementary material.

Figure 2 shows the time series of (i) the macroscopic stress and (ii) the damage rate for simulations with *α* = 1, 3, 4, 10 and 100, and an identical initial relaxation time. After an initial elastic loading phase, the mechanical behaviour differs widely. For *α* = 1, viscous dissipation is hindered as λ is large and does not deteriorate with damage. No abrupt drop of the macroscopic stress is observed after the onset of damage. Instead, the input loading is balanced by the stress relaxations associated with damage, and the states of stress remain always near-critical. Both damage and deformation occur within a dense network of linear-like features (figure 3), but this structure is quenched. The correlation integral for the damage events, *C*(*r*), defined as the probability for a pair of events to be separated by less than a distance *r*, scales as , where the correlation dimension *D* is close to the topological dimension of the system (*D* ≈ 2); whatever the duration of the time window over which the damage activity is considered (the departure from this scaling towards large *r* results from a finite-size effect). This indicates an essentially spatially homogeneous distribution of damage, i.e. a weak localization, whatever the time scale considered.

For *α* > 1, the macroscopic stress exhibits abrupt relaxations corresponding to peaks in the damage activity and to the formation of faults, or ‘leads’, where deformation rates are orders of magnitude higher than in undamaged areas of the plate (figure 3). For *α* = 3, a macroscopic fault remains activated throughout the entire simulation and the internal stress is almost always near-critical. In this case, when considering short time windows, *C*(*r*) significantly departs from the homogeneous scaling (figure 3*b*, black curve), consistent with the strong concentration of damage within these linear structures. Similar to creeping faults in the Earth's crust showing a continuous micro-earthquake activity [77], the structure and orientation of the fault are maintained by low-amplitude and high-frequency fluctuations of the damage activity within the fault. An episodic, spatially spread damage activity also occurs (see the electronic supplementary material, movie_3). This explains that damage heterogeneity decreases as the considered time scale increases (figure 3*b*, grey curves), a first indication of space/time coupling in the simulated dynamics. At all frequencies *f* beyond that associated with the prescribed healing time *t*_{h}, the power spectral density (PSD) of the damage-rate time series shows a power-law decay , reminiscent of temporal correlations and consistent with a marginally stable state (figure 4). For *f* < 1/*t*_{h}, the PSD is flat, indicating that the memory of previous damage events is lost within the Maxwell-EB material when healing is completed.

For larger values of α (e.g. *α* = 4), the mechanical behaviour is partitioned between (i) periods of creep-like deformation during which the *same* main fault remains activated and (ii) cycles of slow stress build-ups due to the healing of existing faults and rapid relaxations associated with the activation of *new* faults, with different shapes and locations (see the electronic supplementary material, movie_4). This translates into a broad peak on the PSD of the damage-rate time series and in a rapid spatial homogenization of damage when cumulating several macro-rupture events. Consequently, the correlation dimension shifts from *D* ≈ 1.2 when considering a small time window (and so a single macro-rupture event) to *D* ≈ 2 at much larger time scales (figure 3*b*). This loss of spatial correlation, as one increases the time scale, is fully consistent with the space/time coupling of sea ice deformation described in §3.

Both types of deformation coexist for a range of *α*-values. Beyond this range (e.g. *α* = 10), episodes of creep-like deformation disappear and damage clusters in time, i.e. there is almost no precursory damage activity to the propagation of macro-ruptures (see the electronic supplementary material, movie_10). The *same* fault is formed at each stress build-up and relaxation cycle, equivalent to re-starting the experiment from its initial undamaged state, thereby suggesting the total suppression of the memory of previous damage events within one cycle. Consistent with this observation, the PSD progressively flattens as *α* increases, indicating the loss of the long-term memory of damage. Thus *C*(*r*) does not exhibit a clear power law in this case, as a randomly distributed (uncorrelated) damage activity precedes each re-activation of the fault.

For *α* = 100, deformation over damaged areas of the plate becomes readily permanent, whatever their level of damage: as expected, the model becomes essentially elasto-plastic. The PSD is characterized by a single damage frequency, corresponding to the healing frequency. At the onset of damage, the numerical scheme does not always fully converge within the convergence rules we prescribed, impacting the spatial distribution of damage, which is, therefore, not analysed here.

## 6. Concluding remarks

The scaling analyses performed over the last years on sea ice mechanical observables revealed a brittle mechanics that accommodates deformation imposed by the forcing fields (such as winds and currents) through episodic fracturing and faulting events over a wide range of space and time scales. This mechanics, which shares many similarities with other geophysical systems (e.g. the Earth's crust), is characterized by a complex coupling between space and time scaling laws that prohibits the use of homogenization procedures to link scales.

Based on these observations as well as theoretical considerations, we built a novel sea ice rheological model called Maxwell–elasto-brittle [76]. This framework considers fast elastic stress redistributions following brittle damage as well as a slow, post-damage viscous-like relaxation. The competition between these mechanisms is controlled by the level of damage of the ice. By varying the parameter *α*, the proposed rheological framework can represent widely different mechanical behaviours, from pathological ones (*α* ≤ 1) to elasto-plasticity (). In this last limit, the memory of elastic deformations is lost and the macroscopic behaviour becomes predictable, with a periodic (re)activation of a single fault pattern. In this case, any intermittency in the system would be directly inherited from a complex applied forcing (e.g. winds), not from the mechanical framework itself. Conversely, with intermediate *α*-values, the model integrates all the cornerstones for the emergence of a completely unpredictable (in the deterministic sense), marginally stable state characterized by extreme heterogeneity, intermittency, anisotropy and space/time coupling. In other words, it simulates ‘by itself’ all the characteristics of sea ice mechanical complexity. If the simulations presented here under very simple loading configurations demonstrate the ability of the model to reproduce the mechanical and dynamical properties of the sea ice cover, the next step will be the realization of, for example, pan-Arctic simulations introducing a coupling with a thermodynamic sea ice model as well as realistic forcing fields (winds and ocean currents). This will allow one to better constrain the model parameters and particularly *α*. This is left for future work.

## Authors' contributions

J.W. supervised the project and wrote the paper. V.D. developed the Maxwell-EB model and performed the simulations. All authors contributed to the discussions and revised the paper.

## Competing interests

We have no competing interests.

## Funding

This work has been partly supported by Total E&P Recherche Developpement.

## Acknowledgements

We thank Pierre Saramito for his expertise on the modelling part of the work, and the use of the modelling platform RHEOLEF.

## Footnotes

One contribution of 11 to a theme issue ‘Microdynamics of ice’.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3579836.

- Accepted July 9, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.