Véronique Dansereau
email: veronique.dansereau@lgge.obs.ujf-grenoble.fr

Jérôme Weiss

Pierre Saramito

Philippe Lattes

A Maxwell-Elasto-Brittle rheology for sea ice modelling

A new rheological model is developed that builds on an elasto-brittle (EB) framework used for sea ice and rock mechanics, with the intent of representing both the small elastic deformations associated with fracturing processes and the larger deformations occurring along the faults/leads once the material is highly damaged and fragmented. A viscous-like relaxation term is added to the linear-elastic constitutive relationship together with an effective viscosity that evolves according to the local level of damage of the material, like its elastic modulus. The coupling between the level of damage and both mechanical parameters is such that within an undamaged ice cover the viscosity is infinitely large and deformations are strictly elastic, while along highly damaged zones the elastic modulus vanishes and most of the stress is dissipated through permanent deformations. A healing mechanism is also introduced, counterbalancing the effects of damaging over large time scales. In this new model, named Maxwell-EB after the Maxwell rheology, the irreversible and reversible deformations are solved for simultaneously, hence drift velocities are defined naturally. First idealized simulations without advection show that the model reproduces the main characteristics of sea ice mechanics and deformation: strain localization, anisotropy, intermittency and associated scaling laws.

Introduction

Making reliable predictions of the drift and deformation of sea ice is becoming crucial nowadays for: (1) forecasting the opening of shipping routes across the Arctic, (2) evaluating mechanical constraints on offshore structures and ships and, at larger scales, (3) estimating the future evolution of both the summer and winter sea ice cover in the Arctic and Antarctic to anticipate its short to long-term, regional to global impacts on climate. Current operational modelling platforms, whether assimilating data of not (e.g., TOPAZ4 : [START_REF] Sakov | TOPAZ4: An ocean-sea ice data assimilation system for the North Atlantic and Arctic[END_REF], GIOPS : [START_REF] Smith | Sea ice forecast verification in the Canadian Global Ice Ocean Prediction System[END_REF]), and global 1 climate models including sea ice dynamics (e.g., the Coupled Model Intercomparison Project Phase 5 models involved in the IPCC Fifth Assesment Report [START_REF] Flato | Climate Change 2013: The Physical Science Basis[END_REF]) are based on the same mechanical framework for sea ice developed in the late seventies: the Hibler Viscous-Plastic (VP) model [START_REF] Hibler | A viscous sea ice law as a stochastic average of plasticity[END_REF][START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF]. With this approach, the ice creeps very slowly as a viscous fluid under small stresses and deforms plastically once exceeding a yield criterion. Yet, over the last decade, the viscous hypothesis and other underlying physical assumptions of this VP framework have been revisited and found inconsistent with the observed mechanical behaviour of sea ice [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF][START_REF] Coon | Arctic Ice Dynamics Joint Experiment (AIDJEX) assumptions revisited and found inadequate[END_REF][START_REF] Rampal | Scaling properties of sea ice deformation from buoy dispersion analysis[END_REF]. In the same line of ideas, recent modelling studies have demonstrated that while the VP model can represent with a certain level of accuracy the mean, global (> 100 km) drift of sea ice, it fails at reproducing the observed properties of sea ice deformation and that, especially at the fine scales [START_REF] Lindsay | Sea-ice deformation rates from satellite measurements and in a model[END_REF][START_REF] Kwok | Variability of sea ice simulations assessed with RGPS kinematics[END_REF][START_REF] Girard | Evaluation of high-resolution sea ice models on the basis of statistical and scaling properties of Arctic sea ice drift and deformation[END_REF] relevant for operational modelling, thereby stressing the need to explore alternative rheologies.

Other continuum models have been developed lately with the aim of representing more accurately some important aspects of the mechanical behaviour of sea ice. Considering the discontinuous and anisotropic character of the pack, [START_REF] Schreyer | Elastic-decohesive constitutive model for sea ice[END_REF] have suggested an elastic-decohesive model that explicitly accounts for the deformation arising from discontinuities in displacement across leads, the orientation of which is prescribed. [START_REF] Tsamados | Impact of a new anisotropic rheology on simulations of Arctic sea ice[END_REF] have presented a model based on the rheology of [START_REF] Wilchinsky | Modelling the rheology of sea ice as a collection of diamond-shaped floes[END_REF] that accounts for the subgrid scale anisotropy of the sea ice cover. Their framework incorporates an evolution equation for the orientation of ice floes, for which a diamond shape is assumed. Our present work shares the same objective as these previous initiatives : to build a continuum model for sea ice that is physically consistent with its observed mechanical behaviour. However, we chose to base our approach on a completely isotropic rheology and, by incorporating the relevant brittle mechanics concepts and long-range elastic interactions, aim to develop a model that reproduces the anisotropy and extreme gradients within the sea ice cover naturally, that is, without the need of treating velocity discontinuities explicitly nor prescribing lead orientations or floe shapes.

Early on, sea ice scientists have suspected that the sea ice cover behaves in a brittle instead of a viscous manner, with some strain hardening in compression [START_REF] Nye | Is there any physical basis for assuming linear viscous behavior for sea ice?[END_REF]. Studies of fracture patterns, stresses and strains both in situ and in the laboratory have suggested that the deformation of sea ice is mostly accommodated by a mechanism of multiscale fracturing and frictional sliding [START_REF] Marsan | Scale dependance and localization of the deformation of Arctic sea ice[END_REF][START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF][START_REF] Schulson | Creep and Fracture of Ice[END_REF]Schulson, 2006a). By investigating the dispersion of ice buoys, [START_REF] Rampal | Scaling properties of sea ice deformation from buoy dispersion analysis[END_REF] recently showed that sea ice over the Arctic deforms in a heterogenous and intermittent manner over spatial scales of 300 m to 300 km and time scales of 3 hrs to 3 months. The strong space-time coupling in the scaling laws revealed by their analyses are consistent with (1) a brittle-type material in which permanent deformations are accommodated by displacements along fractures and fault planes over a wide range of scales and (2) long-range elastic interactions, allowing for small, local perturbations to trigger much larger damaging events within the ice pack [START_REF] Marsan | Space/time coupling in brittle deformation at geophysical scales[END_REF].

A close comparison can be made between the deformation of sea ice and that of the Earth crust, in which brittle fracturing and Coulomb stress redistribution also take place and for which scaling properties have been recognized for years [START_REF] Kagan | Spatial distribution of earthquakes: the two-point correlation function[END_REF][START_REF] Kagan | Fractal dimension of brittle fracture[END_REF][START_REF] Kagan | Long-Term Earthquake Clustering[END_REF][START_REF] King | Static stress change and the triggering of earthquakes[END_REF][START_REF] Turcotte | Fractals and Chaos in geology and geophysics[END_REF][START_REF] Stein | The role of stress transfer in earthquake occurrence[END_REF]. Recently, [START_REF] Marsan | Space/time coupling in brittle deformation at geophysical scales[END_REF] established a formal analogy between the mechanical behaviour of sea ice and that of the Earth crust by demonstrating that the space-time coupling in the deformation of sea ice, estimated from continuous displacement fields, is equivalent to a coupled scaling of the discrete ice-fracturing events occurring along the leads, similar to that observed for earthquakes [START_REF] Kagan | Fractal dimension of brittle fracture[END_REF][START_REF] Kagan | Long-Term Earthquake Clustering[END_REF]. The authors suggested that the similarity between sea ice and the Earth crust is attributable to a common cascading mechanism of earth-/ice-fracturing events that extends the influence of local events to longer durations and larger areas than their direct aftershocks.

In the case of rocks, attempts to simulate brittle deformation were first made using random springlike models. Combining local threshold mechanics and long-range elastic interactions, these successfully reproduced the strong localization of rupture in both space and time, the clustering of rupture events along faults and the multifractal properties of strain fields [START_REF] Cowie | Statistical Physics Model for the Spatiotemporal evolution of faults[END_REF][START_REF] Cowie | Multifractal scaling properties of a growing fault population[END_REF]. Building on similar linear-elastic laws and introducing some strain softening at the micro scale, the failure model of [START_REF] Tang | Numerical Simulation of Progressive Rock Failure and Associated Seismicity[END_REF] succeeded in simulating the progressive failure leading to the macroscopic non-linear behaviour of brittle rock, thereby processing discontinuum mechanics by a continuum mechanics method. An analogous approach based on local damage evolution was also taken by [START_REF] Amitrano | From diffuse to localised damage through elastic interaction[END_REF], who combined a linear-elastic constitutive relationship for a continuum solid, a local Mohr-Coulomb criterion for brittle failure, an isotropic progressive damage mechanism for the elastic modulus described by a nondimensional scalar damage parameter, allowing for the redistribution of the stress from overcritical to sub-critical areas of the material, for the triggering of avalanches of damaging events and the for propagation of faults. This rheological framework, named Elasto-Brittle (EB), was recently developed in the context of the Arctic ice pack by Girard et al. (2010a, b) to explicitly introduce brittle mechanics concepts in continuum sea ice models. First implementations of this rheology into short (3-days), no-advection, stand-alone simulations of the Arctic, but using realistic wind forcing from reanalyses, showed that the EB model is able to reproduce the strong localization and the anisotropy of damage within sea ice and agrees very well with the deformation fields estimated from the RADARSAT Geophysical Processor System (RGPS) data (Girard et al., 2010b).

In the context of longer-term simulations of ice conditions and coupling to an ocean component, a suitable sea ice model however needs to represent not only the small deformations associated with the fracturing of the pack, but also the permanent deformations occurring once it is fragmented and undamaged ice floes move relative to each other along open leads, as these much larger deformations set its overall drift patterns and advective processes. This last point is an important and intrinsic limitation of the EB framework, since the linear-elastic constitutive law does not allow solving for the elastic (reversible) and permanent deformations of the simulated material separately.

Hence to estimate the material's velocity, assumptions about the amount of reversible versus irreversible deformation need to be made in the EB model. The partitioning is bounded by two limit cases. (1) If a loading stress is applied to the damaged material (see Fig. 1b, dashed blue loading path) and all of the resulting deformation is assumed elastic, the material goes back to its initial position if unloaded and its velocity is zero (red dashed unloading path). This assumption was made in the no-advection simulations of Girard et al. (2010b). ( 2) Alternatively, if all of the resulting deformation is considered permanent, the material keeps its final position if unloaded (Fig. 1b, purple dashed unloading path) and the velocity is trivially estimated as the ratio of the total deformation and of the time associated with the loading. In the case of sea ice, the second assumption might be justified by the fact that elastic deformations within an undamaged pack are small compared to the permanent deformations associated with the opening, closing, and shearing along leads. Considering the maximum in-situ values of shear stress of 10 5 Pa reported by [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF] and an undamaged elastic modulus between 1.0 and 10.0 • 10 9 Pa [START_REF] Timco | A review of the engineering properties of sea ice[END_REF], upper bound values for shear strains in a one meter thick elastic ice pack would be on the order of 10 -5 . On daily time scales, these are at the lower bound of RGPS deformation rate estimates (between 10 -4 and 10 0 day -1 , for [START_REF] Marsan | Scale dependance and localization of the deformation of Arctic sea ice[END_REF]; [START_REF] Girard | Evaluation of high-resolution sea ice models on the basis of statistical and scaling properties of Arctic sea ice drift and deformation[END_REF]), suggesting a dominant contribution of irreversible deformations. This second assumption is taken in the recently developed neXtSIM sea ice model, which is based on the EB rheology and does represent advective processes over the Arctic [START_REF] Bouillon | Presentation of the dynamical core of neXtSIM a new sea ice model[END_REF]. However in this all-permanent deformations limit, internal stresses are immediately dissipated, hence the memory of the stresses associated with elastic deformations is erased whenever the applied loading is removed or reset. Without carrying the history of previous stresses, the model cannot exhibit the intermittency intrinsic to the mechanical behaviour of sea ice, i.e., not directly inhered from the wind forcing [START_REF] Rampal | Arctic sea ice velocity field: General circulation and turbulent-like fluctuations[END_REF]. In order to estimate adequate drift velocities, a suitable rheological model must therefore have the capacity to distinguish between reversible and irreversible deformations.

The goal of this work is to develop such a model allowing a passage between the small/elastic and large/permanent deformations and with the capability of damage mechanics models to reproduce the observed space and time scaling properties of sea ice deformation. Our approach consists in introducing a viscous relaxation term into the linear-elastic constitutive law of the original EB framework. The new constitutive relationship takes the form of the Maxwell viscoelastic model. The all-important difference with respect to the Maxwell framework however is that the viscosity associated with the stress dissipation term is not meant to represent the viscoplastic creep of bulk ice [START_REF] Duval | Rate-controlling processes in the creep of polycrystalline ice[END_REF], but instead is an "apparent" viscosity that depends on the local level of damage and concentration of the ice cover. As the elastic modulus, this mechanical parameter is coupled to the progressive damage mechanism through a scalar variable d representing the time and space-evolving level of damage of the ice pack. The coupling is designed so that stresses induce elastic strains over undamaged portions of the ice and are dissipated through permanent deformations where the pack is highly fractured.

The use of a viscoelastic rheology and apparent viscosity in the case of sea ice can be supported again by the similarity between the mechanical behaviour of the ice pack and that of the Earth crust and the existence of similar approaches to model lithospheric faulting. Active faults in the Earth crust have been known to deform in two distinct ways: either abruptly, causing earthquakes, or in an transient, aseismic manner [START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF][START_REF] Gratier | Postseismic pressure solution creep: Evidence and time-dependent change from dynamic indenting experiments[END_REF][START_REF] Cakir | Onset of aseismic creep on major strike-slip faults[END_REF][START_REF] Cetin | Extent and distribution of aseismic slip on the Ismetpasa segment of the North Anatolian Fault (Turkey) from Persistent Scatterer InSAR[END_REF].

Similar to sea ice, co-seisimic fracturing activates aseismic creep, leading to deformations that can be much larger than that associated with the fracturing itself and to the relaxation of a significant amount of elastic strain [START_REF] Cakir | Onset of aseismic creep on major strike-slip faults[END_REF][START_REF] Cetin | Extent and distribution of aseismic slip on the Ismetpasa segment of the North Anatolian Fault (Turkey) from Persistent Scatterer InSAR[END_REF]. A further justification of the introduction of such pseudo-viscosity comes from the rheology of granular media. As sea ice along leads (see Fig. 3), rocks along active faults are highly fragmented. Sheared granular media flow in a viscous manner when inertial effects can be neglected [START_REF] Jop | A constitutive law for dense granular flows[END_REF] with an apparent viscosity diverging as the packing fraction approaches the close-packed limit [START_REF] Aranson | Patterns and collective behavior in granular media: Theoretical concepts[END_REF]. This last point will justify the dependence of our apparent viscosity on sea ice concentration.

Viscous-elastic rheological models using apparent viscosities have already been used to model the deformation of rock-like materials. [START_REF] Lyakhovsky | Distributed damage, faulting and friction[END_REF] built a viscoelastic damage rheology model with the intent of representing the different stages of geological faulting, from subcritical crack growth to increasing crack concentration and material degradation, macroscopic brittle failure, post failure deformation and healing. However, the evolution of damage in their model was derived from energy conservation principles rather than from a brittle failure criterion and was coupled to the elastic modulus only. [START_REF] Frederiksen | Numerical modelling of strain localisation during extension of the continental lithosphere[END_REF] successfully simulated strain localization during lithospheric extension using an elasto-visco-plastic model together with an ad hoc viscosity.

As their work was concerned with the ductile rather than the brittle deformation regime, strain softening in their model did not involve a progressive damage mechanism but instead was achieved by coupling the viscosity to the accumulated strain and the elastic modulus of the material was kept constant. [START_REF] Hamiel | A viscoelastic damage model with applications to stable and unstable fracturing[END_REF] modified the coupled linear elasticity and progressive damage rheological framework of [START_REF] Lyakhovsky | Distributed damage, faulting and friction[END_REF] with a non-linear damage-elastic moduli relation and by adding a damage-dependent Maxwell-like viscous term to account for the gradual accumulation of irreversible strain observed in typical rock mechanics experiments. The addition of this term had however a fundamentally different purpose than in the present approach in that it was intended for the representation of the small pre-macroscopic brittle failure deformations, not to bridge between small and large deformations.

To our knowledge, it is therefore the first time a viscoelastic Maxwell constitutive law is coupled to a progressive damage (and healing) mechanism through both the elastic modulus and an apparent viscosity with the intent of reproducing the small deformation associated with brittle fracturing and the large, permanent post-fracture deformation of geomaterials. It is certainly the first time such a rheological model has been adapted for sea ice modelling.

The paper is structured as follow: the Maxwell-EB rheological framework is described in section 3. A dynamical Maxwell-EB sea ice model is presented in section 4 along with its adimensional version and a discussion of the important non-dimensional numbers involved in the model. The numerical scheme employed in the case of small-deformation experiments is presented and idealized model simulations are described in section 5. In section 6, these simulations are analyzed and discussed on the basis of the macroscopic behaviour and convergence properties of the model and of the heterogeneity, anisotropy and intermittency of the simulated deformation. Conclusions are summarized in section 7.

3 The Maxwell-EB model

Constitutive relationship

The Maxwell rheology describes the behaviour of a continuum material exhibiting both elastic and viscous properties and combines a Newtonian viscous fluid-like damper and a linear elastic term, typically represented by a spring and dashpot connected in series (see Fig. 1a). Considering the material, typically an incompressible fluid, as being isotropic at the elementary scale for both elastic and viscous properties and assuming plane stress conditions, the Maxwell constitutive relationship reads

1 G Dτ Dt + 1 η τ = 2ε (1)
with τ , the deviatoric part of the Cauchy stress tensor, G and η the (shear) elastic modulus and viscosity of the material associated to the spring and dashpot components respectively, ε the strain rate tensor, defined in terms of the velocity u. The objective derivative of the stress tensor is given by

Dτ Dt = ∂τ ∂t + (u • ∇)τ + β a (∇u, τ ) (2) with β a (∇u, τ ) = τ W (u)-W (u)τ -a (τ D(u) + D(u)τ ), D(u) = ∇u+∇u T 2 and W (u) = ∇u T -∇u 2
the symmetric and anti-symmetric parts of the velocity gradient and a = 0, 1 or -1 if using the Jaumann, upper convected or lower convected objective derivative.

When a stress τ is applied to the Maxwell system, the resulting deformation ε total is split between two components: the instantaneous, reversible, deformation of the spring, ε E , and the permanent deformation of the dashpot, ε η , increasing linearly in time (see Fig. 1a). For a given total deformation applied to the system, the rate of dissipation of the associated stress through the permanent deformation of the dashpot is determined by the ratio, η G , of the viscosity of the dashpot and of the (shear) elastic modulus of the spring, G. This ratio can be interpreted as a characteristic memory time for elastic deformations : as it decreases, the material looses its capacity to retain the memory of recoverable deformations.

Here we apply the idea of stress dissipation to a compressible, elastic continuum solid and formulate the following constitutive equation by adding a Maxwell-like viscous damper term to the linear elasticity (i.e., Hooke's law) constitutive relationship :

1 E ∂σ ∂t + (u • ∇)σ + β a (∇u, σ) + 1 η σ = K : ε(u) ( 3
)
where σ is the total Cauchy stress tensor, E is the elastic (or Young) modulus and K is the (adimensional) elastic stiffness matrix, which in terms of the Poisson ratio

ν writes K = ν (1+ν)(1-2ν) I ⊗ I + 2 2(1+ν)
I with I the second rank identity tensor and I the symmetric part of the fourth rank identity tensor. In this rheological framework, the mechanical parameter η is not the true dynamic viscosity of the material but rather is an "apparent" viscosity. The related relaxation time, λ = η E , characterizing the rate at which internal stresses dissipate into permanent deformations, is assumed equal for both the volumetric and deviatoric components of the deformation of the compressible material.

Damage criterion

In agreement with in-situ stress measurements [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF], and as in the original EB model, the damage criterion in the Maxwell-EB rheology is based on the Mohr-Coulomb (MC) theory of fracture. In terms of the principal stress components σ 1 and σ 2 , and using the rock mechanics convention that compressive stresses are positive, the MC criterion reads

σ 1 = qσ 2 + σ c (4)
(or σ 2 = qσ 1 +σ c , by symmetry of the criterion along the the σ 1 = σ 2 axis -see Fig. 2). The slope of the envelope in the principal stresses plane, q, is expressed in terms of the internal friction coefficient

µ as q = (µ 2 + 1) 1/2 + µ 2 . ( 5
)
The intercept σ c of the MC criterion with the σ 1 axes (see Fig. 2), interpreted as the uniaxial (unconfined) compressive strength, is given by

σ c = 2C (µ 2 + 1) 1/2 -µ . (6)
" with the cohesion C setting the local resistance of the material to pure shear. Disorder is introduced in the damage criterion through the spatial distribution of C. This noise represents the material's natural heterogeneity that causes progressive failure behaviour (e.g., [START_REF] Amitrano | From diffuse to localised damage through elastic interaction[END_REF]Herrmann and Roux, 1990;[START_REF] Tang | Numerical Simulation of Progressive Rock Failure and Associated Seismicity[END_REF] under homogeneous forcing conditions and is associated with structural defects at the sub-grid scale, thermal cracks in sea ice for instance, serving as stress concentrators [START_REF] Schulson | Creep and Fracture of Ice[END_REF]. No correlation length is associated to these heterogeneities, hence their spatial scale corresponds to the spatial resolution of the model, ∆x [START_REF] Hutchings | Modeling Linear Kinematic Features in Sea Ice[END_REF][START_REF] Bouillon | Presentation of the dynamical core of neXtSIM a new sea ice model[END_REF], and the value of C over each model element is drawn randomly from a uniform distribution of values spanning estimates from in-situ stress measurements in Arctic sea ice [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF]. The internal friction coefficient µ is set to 0.7, a value seemingly scale-independent and consistent with laboratory experiments on Coulombic shear faults in fresh ice (Schulson et al., 2006b;[START_REF] Fortt | The resistance to sliding along Coulombic shear faaults in ice[END_REF][START_REF] Weiss | Coulombic faulting from the grain scale to the geophysical scale: lessons from ice[END_REF] and also common for geomaterials [START_REF] Byerlee | Friction of rocks[END_REF][START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF].

" E " ⌘ " total E < E 0 E 0 " ⌘ ⌘ G ⌧ (a) (b)
For metals and rocks, the MC theory was shown to be defective in the case of tension [START_REF] Paul | A modification of the Coulomb-Mohr theory of fracture[END_REF], as the mechanism of tensile failure is intrinsically different to that of compressive failure and, in general, does not involve friction. In the case of σ 1 , σ 2 < 0, fracture occurs whenever σ 1 or σ 2 reaches a critical value. However, in-situ stress measurements in Arctic sea ice have revealed that pure tensile failure does not significantly modify the Coulombic-like failure envelope of pack ice and that Coulomb branches well describe this envelope even under large tensile stresses, up to at least σ N ∼ 50 kPa [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF]. Here, we therefore extend the Mohr-Coulomb criterion to tensile stresses and for practical reasons, set the critical value to the ultimate tensile stress σ t , defined as the intersection of the Mohr-Coulomb criterion with the σ 2 axis [START_REF] Paul | A modification of the Coulomb-Mohr theory of fracture[END_REF], as shown on Fig. 2. The tensile strength cutoff therefore takes the form:

σ 1 < 0; σ 2 = σ t , (7)
where

σ t = - σ c q = -2C (µ 2 + 1) 1/2 + µ . (8)
This gives a ratio of the ultimate tensile stress and uniaxial compressive stress of σt σc ≈ 0.27, which might slightly overestimate the tensile strength for sea ice as measured on the field [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF] and in the lab (Schulson, 2006a) (σ t ≈ 0.2σ c ). However, as such large values of tensile strength are rarely obtained in the Maxwell-EB model simulations, this choice does not significantly affect our results.

No truncation to the MC criterion is used to close the envelope towards biaxial compression (i.e., beyond σ c ) as instances of large biaxial compressive stresses are seldom encountered in Arctic sea ice [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF]. Besides, imposing a truncation was shown to have little impact on the simulation results. The damage criterion combining the MC envelope and the tensile strength cutoff is represented in Fig. 2 in the principal stresses plane and has the same shape as deduced by [START_REF] Coon | Arctic Ice Dynamics Joint Experiment (AIDJEX) assumptions revisited and found inadequate[END_REF] from measurements in undamaged pack ice.

Progressive damage mechanism and healing

The Maxwell-EB rheology differs from the standard Maxwell rheology in that the mechanical parameters E, η and λ are not constant but all coupled to the spatially and temporally evolving level of damage of the material, which controls its local degradation and re-increase in strength. Consistent with previous damage rheological frameworks, the level of damage is represented by a nondimensional, scalar parameter d evolving between 1 (undamaged) and 0 ("completely damaged" material). This variable is interpreted as a measure of sub-grid cell defects or crack density (Kemeny and Cook, 1986) and is allowed to evolve through two competing mechanisms : damaging and healing. On the one hand, damaging represents fracturing and the opening of faults, or "leads" in the case of sea ice, occurring when and where the internal stress exceeds the mechanical resistance of the material and which leads to its weakening. Healing on the other hand represents the reconsolidating and strengthening of the damaged material through sintering or, in the case of sea ice, refreezing within open leads. Although this mechanism also contributes to the increase in elastic stiffness (E × h) and effective apparent viscosity (η × h) of the ice, healing is distinguished from pure thermodynamic growth or dynamically-driven thickness redistribution (e.g., [START_REF] Rothrock | The energetics of the plastic deformation of pack ice by ridging[END_REF]) The calculation of the distance to the damage criterion dcrit, defined by the intersection (σ 1 , σ 2 ) of the line relating the state of stress (σ1, σ2) of a given element to the origin of the principal stress plane, is represented in red in the case of exceeding the Mohr-Coulomb criterion and in purple, the tensile strength criterion.

4 3 2 1 1 1 2 2 3 3 ( 1 , 2 ) ( 0 1 , 0 2 ) c t c t 2 ( 1 , 2 ) ( 0 1 , 0 2 ) 1 = 2 4 1 q 2 = 0 2 q 1 = 0
in that it applies only where and when the material has been damaged. It therefore allows d, E and η to re-increase at most to their undamaged value; d 0 = 1, E 0 and η 0 respectively. Because the two processes operate simultaneously within the simulated material, an evolution equation for d needs to include both mechanisms. In the following damaging and healing are first treated separately and then combined in a single equation for d.

Damaging

Contrary to typical sea ice modelling frameworks, no plastic (i.e., normal) flow rule is prescribed when the damage criterion is reached in the Maxwell-EB model. Instead, when the stress locally exceeds the critical stress, the elastic modulus is allowed to drop, leading to local strain softening (e.g., [START_REF] Amitrano | From diffuse to localised damage through elastic interaction[END_REF][START_REF] Cowie | Statistical Physics Model for the Spatiotemporal evolution of faults[END_REF][START_REF] Tang | Numerical Simulation of Progressive Rock Failure and Associated Seismicity[END_REF]Hamiel et al., 2004, and others). Because of 290 the long-range interactions within the elastic medium, local drops in E imply a stress redistribution that can in turn induce damaging of neighbouring elements. By this process, "avalanches" of damaging events can occur and damage can propagate within the material over long distances [START_REF] Amitrano | From diffuse to localised damage through elastic interaction[END_REF]Girard et al., 2010a). As the elastic perturbation generated by such events is anisotropic [START_REF] Eshelby | The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems[END_REF], this propagation mechanism naturally leads to the emergence of both spatial heterogeneity and anisotropy in the stress and strain fields, i.e., to the formation of linear-like faults (see section 6).

In

ε = ε ←→ K -1 σ E × d crit = K -1 σ E ,
where the superscript denotes the post-damage state of deformation and stress. In terms of the principal stress components, the change in level of damage of a given element is given by

d crit = σ 1 σ 1 = σ 2 σ 2 , (9)
which implies that as the level of damage varies, all stress components vary in the same proportions.

Hence the state of stress σ after each damaging event is given by the intersection of the failure envelope and of the line connecting the pre-damage state of stress (σ 1 , σ 2 ) with the origin, in the principal stress plane (see Fig. 2). Two cases must be distinguished when calculating σ , depending on which of the Mohr-Coulomb or tensile criterion has been exceeded. Combining the two, d crit is evaluated simultaneously over all mesh elements of the model domain as:

d crit = min 1, σ t σ 2 , σ c σ 1 -qσ 2 . ( 10
)
Following progressive damage models, the level of damage of a given element in the Maxwell-EB model at any given time is determined by both its instantaneous distance to the damage criterion d crit , i.e., its current state of stress, and its previous damage level. This implies that the variable d carries the entire history of damage of model elements and, if discretizing time as t n = n∆t, n ≥ 0, translates into the discrete recursive equation

d n+1 = d n+1 crit d n , 0 < d 0 ≤ 1.
A continuous evolution equation for d can be obtained by considering that the time characterizing the redistribution of stress between model elements is intrinsically tied to the speed of propagation of elastic waves, c, in the material, which carry the damage information. Using a Backward explicit scheme of order 1, and setting the model time step to ∆t = t d with t d = ∆x c , the exact time of propagation of an elastic wave with speed c over a distance ∆x, the following equation arises:

Dd Dt = d crit -1 t d d.
(11)

Healing

By healing, the simulated material is allowed to regain some strength. The characteristic time for this process is designated in the following by t h . It corresponds to the time required for a completely damaged element (d = 0) to recover its initial stiffness (d = 1), which in a dynamic-thermodynamic sea ice model would depend on the local difference between the temperature of the air near the surface of the ice and the freezing point of seawater below. Healing schemes of varying level of complexity could be used in the Maxwell-EB model. One possibility is the one employed in the EB sea ice model of Girard et al. (2010a), which follows parameterizations of the vertical growth of sea ice [START_REF] Maykut | The surface heat and mass balance[END_REF]). An underlying assumption is that the rate of healing is inversely proportional to the level of damaging of the ice. However as there is no physical evidence for this assumption, in the following, uncoupled, implementation of the Maxwell-EB model we use an even simpler parameterization that implies a constant healing rate, 1 t h :

Dd Dt = 1 t h , 0 ≤ d ≤ 1. (12)
Combining both the damaging and healing mechanisms (Eq. 10, 11 and 12), the complete evolution equation for d is

∂d ∂t + (u • ∇)d = min 1, σ t σ 2 , σ c σ 1 -qσ 2 -1 1 t d d + 1 t h , 0 < d ≤ 1, (13)
Although the two processes apply simultaneously on the level of damage in the model, they are inherently distinct. On the one hand, damaging is a discrete threshold mechanism applying only where and when the state of stress becomes overcritical. As mentioned in sections 3.2 and 3.3.1, the characteristic time for this process, t d , is tied to the speed of propagation of (shear) elastic waves and to the model's spatial resolution. In the case of an heterogeneous ice pack, an average value for c is on the order of 500 m/s [START_REF] Marsan | Low frequency bursts of horizontally-polarized waves in the Arctic sea-ice cover[END_REF]. For spatial resolutions between that of current global climate and high resolution regional sea ice models (∆x = 1 to 100 km), the characteristic time for damaging, t d therefore varies between O(1) and O(10 2 ) s. Healing on the other hand is a continuous process acting on all model elements, independently of the local distance to the damage criteria. Studies on the refreezing within leads in sea ice showed that the time for 1 meter of ice to grow within an opening of 10 cm under atmospheric temperatures of

T a = -15 • C is of O(100)
hours or O(10 5 ) seconds [START_REF] Petrich | Formation and structure of refrozen cracks in land-fast first-year sea ice[END_REF]. The orders of magnitude of difference between t h and t d therefore imply that the two processes are intrinsically decoupled in the case of the ice pack.

Coupling d with E and η

The coupling between the Maxwell-EB constitutive relationship and the progressive damage mechanism constitutes one of the main features of this new modelling framework. It is defined such that:

deformations within an undamaged medium are small and reversible, i.e., strictly elastic.

Hence undamaged portions of the simulated material have a maximum elastic modulus E 0 and a very large apparent viscosity η 0 . In this case, the viscous term in ( 3) is negligible and a linear-elastic constitutive relationship is recovered (Fig. 3, right panel),

deformations can accumulate over highly damaged areas of the material to become arbitrarily large. These deformations are permanent and dissipate most of the the stress applied to the material within a short relaxation time. Hence the elastic modulus, viscosity and relaxation time drop locally over damaged areas. In the limit of a completely damaged material, elastic interactions are hindered and deformations are strictly irreversible (Fig. 3, left panel). In this case, λ -→ t d and a soft elastic-plastic behaviour is recovered in which the memory of the elastic stresses is totally lost (narrow-dashed blue line on Fig. 1).

as damaged areas are allowed to heal, E, η and λ all re-increase, up to their initial undamaged values.

Different functions could be used to express the dependence of E, η and λ on d that meet these criteria. In the absence of physical evidences for a higher level of complexity, and consistent with the relationship between the elastic modulus and crack density used in damage models of rocks [START_REF] Agnon | Damage distribution and localization during dyke intrusion[END_REF][START_REF] Amitrano | From diffuse to localised damage through elastic interaction[END_REF][START_REF] Schapery | Nonlinear viscoelastic and viscoplastic constitutive equations with growing damage[END_REF], we use the simplest parameterization and set

E(t) = E 0 d(t) η(t) = η 0 d(t) α , such that λ(t) = η 0 E 0 d(t) α-1
Undamaged! pack ice small, elastic deformations ! 0

Fractured,! damaged zones large, permanent deformations with α a constant greater than 1 introduced to fulfil the constraint that the relaxation time for the stress also decreases with increasing damage and re-increases with healing, as the material respectively looses and recovers the memory of reversible deformations. Using this formulation, both η and E are entirely defined by their initial value, a constant, and by the level of damage variable d.

⌘ ! ⌘ min 1 E ! 0 1 ⌘ ! 0 E = E 0
However, the constitutive equation becomes undefined in the limit of d → 0. This problem can be handled by imposing a fixed minimum value d min > 0 for the level of damage. Alternatively, a cutoff η min η 0 on the value of the apparent viscosity can be introduced and the expression for η(d)

modified as η = (η 0 -η min )d α + η min =      η 0 for d = 1, η min for d = 0. ( 14
)
Substituting for η in the expression for the relaxation time, the elastic modulus then becomes

E = η 0 -η min η 0 E 0 d + η min η 0 1 d α-1 E 0 , (15)
such that E ≈ E 0 for d = 1, E decreases with d until a minimum at d(E = E min ) = ηmin η 0 -ηmin (α -1)

1 α
and E → ∞ for d → 0. Using such a cutoff on η, the elastic term in the Maxwell-EB constitutive equation therefore vanishes in the limit of a "totally" damaged material and the rate of viscous dissipation is then set by the minimum viscosity η min . It is important to note that this limit has no physical significance in the context of a progressive damage model for a continuum solid and is rather introduced to insure mathematical consistency while retaining a continuous function for the level of damage. In the following implementation of the model, we take this approach instead of imposing a minimum value of d, but it had really no impact on our results since in the simulations presented here d > d(E = E min ) at all times.

A Maxwell-EB sea ice model

In this section, the Maxwell-EB rheology is implemented in the context of sea ice modelling. As in regional and global sea ice models, the ice cover is considered as a 2-dimensional plate due to its very large aspect ratio and a constant healing rate is assumed. In this case, the complete dynamical model is given by the following system of equations:

1. The momentum equation:

ρh ∂u ∂t + (u • ∇)u = F ext + ∇ • (σh), (16)
with u the velocity, h the thickness and ρ the density of sea ice. F ext assimilates all external forces on the sea ice cover, which in regional and global sea ice models are typically the air and ocean drags and the forces associated with the Coriolis acceleration and gradients in sea surface height. We assume the internal stress to be homogeneously distributed within the depth h and following [START_REF] Bouillon | Presentation of the dynamical core of neXtSIM a new sea ice model[END_REF] and [START_REF] Sulsky | Using the material-point method to model sea ice dynamics[END_REF], we write the momentum equation in terms of the internal stress rather than the vertically integrated stress tensor commonly used in the sea ice modelling community.

2. Conservation equations for the ice concentration A and ice thickness h:

∂h ∂t + ∇ • (hu) = S h , (17)
∂A ∂t + ∇ • (Au) = S A , (18)
where S h and S A represents thermodynamic source and diffusion terms and elastic compressibility effects are assumed negligible relative to dynamic variations of the ice volume in the conservation of the mass of the sea ice cover.

3. The constitutive relationship (3) with

E = f 1 (E 0 , η 0 , η min , d) exp [-c * (1 -A)], (19)
η = f 2 (η 0 , η min , d) exp [-c * (1 -A)], (20)
where f 1 and f 2 represent the functional dependance on the level of damage of the ice d,

given by ( 15) and ( 14) respectively. The exponential function of the ice concentration allows the internal stress term to be maximal when A = 100% and to decrease rapidly when leads open and A drops. It is of the same form as that used for the pressure term (P , or ice strength in compression) in the VP rheology of [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF]. Here the non-dimensional parameter c * characterizing this dependence on the ice concentration has the same (constant) value for both mechanical parameters, but could be set different in a refined parameterization.

4. The equation for the evolution of damage (13) with the damage criterion defined by Eq. ( 4)

and Eq. ( 7) and q, σ c and σ t given by ( 5), ( 6), (8) in terms of the cohesion variable C and of the constant internal friction coefficient µ. In the case of "quenched disorder" (i.e., when the field of C is set at the beginning of a model simulation), an additional equation arises that handles the advection of the field of cohesion with the 425 simulated velocity field. Table 1 lists all model variables and parameters.

Characteristic numbers and times

Neglecting all thermodynamic effects and variations in ice thickness and concentration (considering h = 1 and A = 100%) as well as external forcings and adimensionalizing with respect to the ice velocity U , the horizontal extent of the model domain L, the thickness of the ice H and the undamaged elastic modulus E 0 , the dynamical system of equations read

Ca 0 Dũ D t = ∇ • σ (21) We 0 dα-1 Dσ D t + σ = We 0 dα K(ν) : ε(ũ) (22) D d D t = min 1, Σ t 1 σ2 , Σ c 1 σ1 -q σ2 -1 1 T d d + 1 T h , 0 < d ≤ 1. ( 23
)
where dα = 1 -ηmin η 0 dα + ηmin η 0

and the superscript '˜' is used for all non-dimensional variables and operators.

In this form, the model involves 8 characteristic numbers and time scales, some of which are constant and some, evolving with the local level of damage d: Ca 0 , the (undamaged) Cauchy number, We = We 0 dα-1 , the Weissenberg number, with We 0 its undamaged value, ν Poisson's ratio, Σ t the dimensionless critical tensile stress, Σ c the dimensionless critical stress with respect to the Mohr-Coulomb criterion, T d the characteristic time for damaging, T h the characteristic time for healing and α the damage constant. In order for the Maxwell-EB model to represent the intended physics, the value of these parameters must evolve within a certain range of values. In the following we elaborate on the absolute and relative values of those numbers which are the most critical in the context of sea ice modelling.

T d

As mentioned in the previous section, the (adimensional) characteristic time for the propagation of damage, T d = t d T with T = L U , is determined by the speed of propagation of elastic waves within the simulated material and is strongly tied to the mean spatial resolution of the model, as t d should be of O( ∆x c ). In turn, this time places a strong constraint on the Maxwell-EB model time step. Setting ∆t < ∆x c is indeed unphysical, as the time associated to one model iteration would then be too short for the stress to be redistributed from one overcritical element to its direct neighbour. For the model to resolve the propagation of damage, the time step must therefore be ≥ t d .

No strict upper bound to ∆t is imposed by the damage mechanism. One the one hand, choosing ∆t > t d could be interesting in terms of reducing computational costs. Physically, it implies that damage is allowed to propagate beyond the first neighbour barrier and over larger distances within one model time step. On the other hand, increasing ∆t with respect to t d also implies (1) a decrease in the resolution of damaging, as the model might miss important intermediate damage events that trigger additional interactions between neighbouring elements (2) larger local drops in the level of damage, inducing large stress perturbations and, potentially, numerical instabilities in the model. Sensitivity analyses on the propagation of the damage should therefore be performed when choosing ∆t > t d . The temporal resolution that is optimal in terms of capturing all elastic interactions within the simulated material is therefore ∆t = t d . In the model experiments presented in the following, this is the choice we make.

T h

In order for healing not to offset damaging in the rate of change of d, the (adimensional) time for healing, T h = t h T , must be much larger than the (adimensional) time for damage propagation. This separation of scales ensures that elements cannot recover by healing more strength than they have lost by damaging within one time step, as excess healing would effectively entail a net growth of the material, a process that is not intended by this parameterization and should instead be accounted for by thermodynamic, balance calculations. Considering the estimates of the speed of elastic waves and of the healing rate of leads aforementioned (section 3.3), pack ice naturally meets this condition.

We

The Weissenberg number, We, defined as the dimensionless product of the viscous relaxation time for the stress and of time T = L U characterizing the deformation process:

We = η E U L = λ T , (24)
sets the viscous versus elastic character of the flow of a viscoelastic material. In the original Maxwell model, We = 0 represents the limit of zero elastic stresses, while a very large We characterizes a strictly elastic solid. In the Maxwell-EB model, the Weissenberg number evolves according to the level of damage as We = We 0 d α-1 with We 0 , its maximum value.

As viscous dissipation should be insignificant over undamaged and strictly elastic areas of the material, We 0 should be chosen very large, representing the limit of 1 η 0 -→ 0. In this case the viscous term in the constitutive relationship (3) effectively vanishes and a linear elastic rheology is recovered. In practice, the value of We 0 is however limited, first, by the machine precision and second, due to a numerical scheme failure known in the field of viscoelastic flow computations as the high Weissenberg number problem [START_REF] Keunings | On the high Weissenberg number problem[END_REF]Fattal andKupferman, 2004, 2005;[START_REF] Saramito | On a modified non-singular log-conformation formulation for Johnson-Segalman viscoelastic fluids[END_REF]. For large values of We, numerical instabilities arise in Maxwell-type models due the presence of deformation source terms (β a ) in the transport equation for the stress tensor (2). With We 0 (or equivalently, λ 0 ) too low, simulations can run for a time t ∼ λ 0 and unphysical viscous dissipation can occur over undamaged parts of the simulated material. To get round this problem, the viscous term in the Maxwell constitutive relationship can be multiplied by a Heaviside function d * that effectively sets 1 η to the limit value of 0 when and where d ≥ d c , with d c a chosen threshold value (e.g., d c = 1 when using a constant heal rate parameterization) and leaves the constitutive equation unchanged (d * = 1) otherwise. In small-deformation experiments, i.e., run for a time t λ 0 , viscous dissipation over undamaged parts of the material is not significant and the inclusion of such a function, unnecessary.

Conversely, where damage becomes important, the viscous relaxation time λ should decrease significantly below the characteristic time for healing to allow for internal stresses to "have time" to dissipate and deformations to become large.

Ca

The dimensionless number that arises when adimensionalizing stresses in the momentum equation with respect to the elastic modulus is the Cauchy number, defined as the ratio of inertial to elastic forces (Ca = ρU 2 E ). If inertial forces are comparable to elastic forces and Ca ∼ 1, the effect of the propagation of viscoelastic waves in the material cannot be neglected. Yet, setting ∆t ≥ t d , that is ∆t at least equal to the period of shear elastic waves, implies that the model does not resolve these waves, but only their consequence of transmitting the damage information within the material.

Hence the wave signal cannot be properly filtered out of the model's solution. In order for the wave contribution not to have a significant effect on the simulated deformation and stress fields, Ca must therefore be 1. Dimensional analysis indicates that over an undamaged ice pack with velocity ranging between 0.001 and 1 m/s, Ca 0 is in the range [10 -12 -10 -6 ]. Hence inertial effects can be safely neglected. For simulated ice velocities U < 1 m/s, and α > 2, inertial effects in the Maxwell-EB model remain negligible when damage becomes important.

α

The damage parameter α controls the rate at which the apparent viscosity decreases and the material looses its elastic properties with damaging. As mentioned in previous sections, it should be set greater than 1 in order for the viscous relaxation time to decrease with damaging. The requirements that (1) the viscous relaxation time drops well below the time for healing over highly damaged areas and (2) inertial effects remain negligible for high deformation rates (i.e., large velocities) can also place a constraint on the minimum value of α. Conversely, for large values of α, the relaxation time λ becomes very small whatever the damage level (see section 3.3.3). This means that elastic deformations are almost immediately dissipated after damaging, that is, the model becomes purely elasto-plastic. For the experiments presented here, we find that α = 4 allows representing both the brittle behaviour and the relaxation of the internal stress within a material with mechanical parameters in the range of the values suitable for sea ice. For α larger than about 7, memory effects become insignificant and the experiment instead exhibits a stick-slip behaviour with a well-defined characteristic frequency (not shown).

• n = 0

⌦ L • n = 0 @⌦ u(x, y) = 0 u(y) = U u(y) = 0

Numerical scheme and experiments

The objective time derivative for the Cauchy stress σ in the Maxwell-EB constitutive relationship

(3) is composed of an inertial term, an advection term and of a sum of rotation and deformation (β a )

terms, each of which implies a different level of numerical complexity. In developing the model, our approach is to introduce each of these terms separately in order to evaluate their respective contribution to the simulated mechanical behaviour. On the one hand, introducing the inertial term while neglecting the advection and β a terms allows retaining a Lagrangian scheme, similar to the original EB model (Girard et al., 2010b). Without any remeshing of the domain, the model is then suitable for short-term, small-deformation simulations only. On the other hand, when permanent deformations accumulate over a long time, the advection term is no longer negligible and β a terms become potentially important.

In the following, we present small-deformation numerical experiments that allow analyzing the mechanical behaviour of the Maxwell-EB model in terms of the statistical and scaling properties of the simulated damage and deformation fields. Performed with a highly idealized configuration for the domain geometry, the applied loading and boundary conditions, these will demonstrate that the main characteristics of sea ice deformation (spatial heterogeneity, anisotropy, intermittency) naturally emerge from the underlying physics and do not need to be implemented in an ad-hoc manner.

The simulations represent the uniaxial compression of a (2-dimensional) rectangular ice plate with dimensions L 2 × L (see Fig. 4a). Compression is applied by prescribing a constant velocity U on the upper short edge of the plate with the opposite edge maintained fixed in the direction of the forcing. No confinement is applied on the lateral sides. The velocity U is small enough to ensure a low driving rate (i.e., slow compared to time scale of damage propagation, [START_REF] Cowie | Statistical Physics Model for the Spatiotemporal evolution of faults[END_REF]).

In the present implementation, the model is not yet coupled to a thermodynamic component, hence

S A = S h = 0.
As advection is neglected and simulations are run for a short enough time such that the macroscopic and local deformations within the ice cover remain small (∼ 1% of the area of model elements), dynamics-induced variations (through convergence-divergence) of the ice thickness and concentration are not accounted for and hence the mechanical parameters E, η and C are not yet coupled to h or A. Conservation of mass is therefore not imposed in these small-deformation simulations, equivalent to assuming a uniform, constant thickness (1 m) and ice concentration (100%). In this case, the system of equations reduces to Eq. ( 21) to ( 23 In all simulations, the time step is set equal to the characteristic time for damage propagation. A semi-implicit scheme is used that linearizes the system, in which the momentum and constitutive equations are first solved simultaneously using a backward Euler scheme of O(1) and the value of d at the previous model time step. The level of damage is updated in a second time using the estimated ũ and σ and an explicit scheme of O(1). A fixed-point algorithm iterates between these two steps until the residual of the linearized constitutive equation drops below a chosen tolerance, ensuring the convergence of the solution. Finite elements and variational methods are used to solve the timediscretized problem on a Lagrangian grid within the C++ environment RHEOLEF (Saramito, 2013: http://cel.archives-ouvertes.fr/cel-00573970). An unstructured mesh with triangular elements is used and the average spatial resolution is set by choosing the number N of elements along the short side of the domain.

All simulations are started from an initially undamaged ice cover with uniform elastic modulus and viscosity. Undamaged mechanical parameter values are chosen so that to be representative of sea ice on regional to global scales (c = 500 ms -1 and ν = 0.3). The undamaged elastic modulus is given by the relation E 0 = 2c 2 (1 + ν)ρ and the undamaged viscosity η 0 is set such that the initial relaxation time λ 0 is as large as possible while the maximum Weissenberg number, We 0 is small (< 1). All model variables and parameters are listed in table 1. Parameter values are not varied in any of the experiments presented here as a sensitivity study is kept for a separate paper.

Results

In this section we analyze the mechanical behaviour of the Maxwell-EB model. In particular, we evaluate its capacity to reproduce the main characteristics of sea ice deformation, which are its spa-tial heterogeneity, intermittency and anisotropy, following the methodology developed in previous observational studies of the deformation and drift of the Arctic ice pack.

One signature of the strong heterogeneity of sea ice deformation is the emergence of a spatial scaling in the deformation fields over a wide range of scales. Using a coarse-graining procedure, [START_REF] Marsan | Scale dependance and localization of the deformation of Arctic sea ice[END_REF] performed a scaling analysis of the deformation of sea ice over the Arctic using the 3-days, 10 km×10 km gridded RGPS deformation product. Doing so, they obtained a power-law relationship between the total deformation rate < εtot > l invariant and the corresponding averaging scale l of the form

< εtot > l ∼ l -β (25)
with a constant exponent β > 0, indicating correlations in the deformation fields over 2 orders of magnitude in l and an increase in the mean strain rate with decreasing scale of observation, in agreement with a strong spatial localization of the deformation.

This coarse-graining calculation was later extended to ice buoy data (e.g., [START_REF] Rampal | Scaling properties of sea ice deformation from buoy dispersion analysis[END_REF][START_REF] Hutchings | Spatial and temporal characterization of sea-ice deformation[END_REF] which, with a higher temporal resolution than the RGPS data, allowed performing scaling analyses of Arctic sea ice deformation in the temporal dimension as well. Using the dispersion rate of buoys as a proxy for the strain rate, Rampal et al. ( 2008) obtained a power-law relationship between the total deformation rate < εtot > t computed at a chosen space scale and the time scale of observation t

< εtot > t ∼ t -γ (26)
with a constant exponent γ > 0 over 2 orders of magnitudes in t (3 hours to 3 months), indicating an increase of strain rates with decreasing temporal scale consistent with an intermittent deformation process. Recently, these temporal and spatial scaling properties have been used as benchmarks to validate (or invalidate) sea ice models (e.g., [START_REF] Girard | Evaluation of high-resolution sea ice models on the basis of statistical and scaling properties of Arctic sea ice drift and deformation[END_REF]Girard et al., , 2010a;;[START_REF] Bouillon | Presentation of the dynamical core of neXtSIM a new sea ice model[END_REF].

An additional and all-important characteristic of the deformation of sea ice that is not captured by these scaling analyses is its strong anisotropy. This property has been made evident since the availability of satellite imagery-derived ice motion products (e.g. [START_REF] Stern | Open water production in arctic sea ice: Satellite measurements and model parameterizations[END_REF], which showed that high strain rates concentrate along oriented, linear-like faults, or leads, often termed "linear kinematic features" [START_REF] Kwok | Deformation of the Arctic Ocean sea ice cover: November 1996 through April 1997[END_REF].

Spatial resolution, convergence and dependence on the initial conditions

In a first time, we analyze the overall, macroscopic behaviour of the Maxwell-EB model, its convergence properties and the dependance of the solution on the prescribed initial conditions. To do so, a set of four uniaxial compression simulations is run using different spatial resolutions, with N = 10, 20, 40 and 80. The values of the initial, undamaged mechanical parameters are identical between the simulations as well as the field of cohesion, which is defined at the lowest resolution (N = 10) and interpolated onto the higher resolution mesh grids. and damaging sequence suggests that the mechanical behaviour is similar to that obtained with other elasto-brittle models (e.g., [START_REF] Tang | Numerical Simulation of Progressive Rock Failure and Associated Seismicity[END_REF][START_REF] Amitrano | From diffuse to localised damage through elastic interaction[END_REF]Girard et al., 2010a). The Maxwell-EB In the Maxwell-EB model, this last stage is characterized by a drop in the Weissenberg number (i.e., in λ) localized along the main fault (not shown), where strain rates are orders of magnitude higher than over undamaged parts of the material. Then, as damaged areas heal, stress builds up again within the material. At all spatial resolutions, the model simulates cycles of slow stress build ups (healing phase) and rapid stress relaxations (damaging phase).

Because the simulations use the same spatial distribution of the damage criteria (i.e., of C) the location of the first damage events is the same at all resolutions, as shown by the maps of instantaneous level of damage d near the onset of damaging (Fig. 5 b1). However, soon after these first failure events, model solutions do not converge (Fig. 5 b2-4) and fractures form with a shape and orientation differing between simulations. This divergence between the post-damage solutions illustrates an all-important and intrinsic characteristic of the Maxwell-EB framework arising from the fact that there is no physical scale associated with the localization of damage in the model. Through elastic interactions, damage and deformation tend to localize at the finest scale (the mesh element), resulting in a different redistribution of the stress between neighbouring elements at different spatial resolutions and hence a non-identical propagation of the damage. Put another way, the divergence of the solutions indicates that while the disorder in C sets the location of the first damage events, the heterogeneities introduced in the stress field by these events prevail in setting the location and timing of subsequent events. This result is consistent with previous elasto-brittle model simulations which have shown that the number of active faults as well as the degree of localization of the deformation over long time scales do not depend systematically on the disorder initially introduced in the model [START_REF] Cowie | Statistical Physics Model for the Spatiotemporal evolution of faults[END_REF] and that once formed, faults produce their own stress field which dominates further fracture growth [START_REF] Tang | Numerical Simulation of Progressive Rock Failure and Associated Seismicity[END_REF].

Another important property of the deformation made evident by this set of experiments is its strong anisotropy. The fields of d and of the total deformation ( εtot , see section 25) represented on Fig. 5 indeed show that at all spatial resolutions, the simulated damage and deformation are both highly localized and oriented along linear features. This is an important result, as no anisotropy is introduced at the local scale on either the elastic or viscous properties, or in the damage parameterization.

This property arises naturally due to elastic interactions within the material and without the need to prescribe fault orientations. It was reproduced by the original EB model [START_REF] Amitrano | From diffuse to localised damage through elastic interaction[END_REF]Girard et al., 2010a, b) and is not lost when including a viscous dissipation term for the stress in the Maxwell-EB constitutive relationship.

Heterogeneity

As shown in the previous section, when simulations are started from an undamaged state, the simulated mechanical behaviour of the material is intrinsically different between the first and subsequent loading and damaging cycles. The path to the first rupture in "irreversible damage" (i.e., models without healing) elasto-brittle models has already been investigated in depth (e.g. Girard et al., 2010a;[START_REF] Tang | Numerical Simulation of Progressive Rock Failure and Associated Seismicity[END_REF][START_REF] Amitrano | From diffuse to localised damage through elastic interaction[END_REF]. Hence we focus our analysis of the spatial dependence of the Maxwell-EB model strain rate fields on the post macro-rupture behaviour.

To quantify the heterogeneity of the simulated deformation, we follow [START_REF] Marsan | Scale dependance and localization of the deformation of Arctic sea ice[END_REF] and estimate deformation rates over two orders of magnitude in space scales using a coarse-graining procedure. The calculation is described in details by Girard et al. (2010a). For this analysis we use the outputs of strain rate fields from simulations with N = 100, averaged over a time interval corresponding to the time of propagation of an elastic shear wave with speed c through the width of the domain ( L 2 1 T ×c = N time steps). The dependence of the deformation rates on the spatial scale of observation is investigated at different stages of the healing-damaging cycle. Figure 6 (a andb) shows the total deformation rate < εtot > l as a function of the space scale l at 5 equally-spaced steps along the path towards a given macroscopic failure event, that is, between the minimum in macroscopic stress that follows the propagation of a fault and the maximum that precedes the next macro-rupture, as indicated in Fig. 6(a).

Deformation rates are normalized by < εtot > at the smallest averaging scale (L/N ). At the first stage, just following the rupture (red curve), the total deformation rate shows a clear power law decrease with increasing spatial scale of the form of Eq. ( 25) over nearly two orders of magnitude of l, consistent with a strong localization of the deformation. At the subsequent stages (yellow and green curves), damaged elements progressively recover their mechanical strength by healing. Deforma-tion rates decrease along the main fault and re-increases over undamaged areas, hence deformation homogenizes over the domain and the rate of decrease of < εtot > l with l is reduced. Then, as healing allows stress to build up within the material, damaging resumes and clusters in space and the exponent β re-increases towards its post macro-rupture value (blue and purple curves).

Repeating the procedure for subsequent healing and damaging cycles and for multiple realizations of the experiment initialized with different cohesion fields showed a similar evolution of the rate of decrease of < εtot > l with l between macro-ruptures events, with values of β in the vicinity of the rupture consistent with previous EB model analyses (e.g., Girard et al., 2010a, β = 0.15 ± 0.02).

However, an important difference between the present results and that of Girard et al. (2010a) is the absence of a clear cross-over scale for which < εtot > l becomes independent of l and which implies a finite correlation length of damage events. This suggests that the Maxwell-EB system progressively looses the memory of it's initial homogeneous, undamaged state and that an elasto-brittle material experiencing both healing and damaging enters a "marginally stable" state with scale invariance spanning the size of the system. This result is consistent with the scale-dependence analysis of RGPS-derived deformation rates of [START_REF] Marsan | Scale dependance and localization of the deformation of Arctic sea ice[END_REF] and [START_REF] Stern | Spatial scaling of Arctic sea ice deformation[END_REF], in which no cutoff scale was observed for l varying between 10 and 1000 km, suggesting that Arctic sea ice is most often in a near-critical state.

Intermittency

In this section we characterize the temporal behaviour of the Maxwell-EB model. Figure 7 First, the evolution of the macroscopic stress is clearly characterized by cycles of slow stress build-ups and very fast relaxations. The strong asymmetry of the signal in time is confirmed by a high (negative) skewness (-6) of the distribution of the macroscopic stress increments ∆σm ∆t (not shown). Associated with these cycles is a succession of progressive increases in damage events and very sharp drops, after which damaging stops momentarily (red arrow on Fig. 7a) .

Second, as identified on the same time series, some periods (e.g., the interval delimited by the dashed red box) are characterized by a continuous damage activity and by both low amplitude and low frequency fluctuations of the stress. This contrasted behaviour translates into a significantly more symmetric (skewness of -1.9) distribution of ∆σm ∆t . Inspection of the spatial distribution of damage (Fig. 7b) and strain rate fields (not shown) over this time interval indicates that the same system of interacting faults remains activated, with not much damaging activity over the rest of the domain and therefore suggests that creep-like deformation along this system dissipates all of the input loading.

Following the approach taken for fracture-type models which record the number of broken fibres, ruptured bounds, depinning events, etc., we investigate the time-dependence of the simulated damage activity by analyzing time series of the discrete failure events. We estimate the power spectral density (PSD) of damage rate time series. The resulting squared Fourier coefficients are averaged over 5 realizations of the compression experiment initialized with different fields of C over domains with N = 40. Fig. 8(a) represents the spectral density estimated by averaging the power over a 5 values window centred on each frequency f . We checked that using a smaller averaging window does not affect the shape of the PSD discussed below.

At low frequencies, the PSD is almost flat, suggesting that the number of damage events is uncorrelated in time. As these frequencies are lower than 1 T h , this is consistent with the fact that the Maxwell-EB material entirely looses the memory of previous damage events when allowed to heal completely. At higher frequencies, the PSD shows a decrease with increasing f reminiscent of a temporal correlation of damaging events in the material. This expresses as a power law decay with PSD(f ) = 1/f γ . At intermediate frequencies, we estimate a slope γ = 2, suggesting that the instantaneous damage rate is correlated in time but increments of the damage rate are uncorrelated. At the highest frequencies, γ > 2, indicating that the damage rate is correlated in time and the of damage rate are anti-correlated. The break in the slope occurs around f = 10 6 , a frequency that we relate to the minimum propagation time of a macro-rupture, i.e., the time of propagation of damage (i.e., of an elastic shear wave with speed c) across the width L 2 of the domain (N time steps). The transition between the flat and power law decaying parts of the PSD is marked by a clear peak spanning the range of frequencies corresponding to the cycles of healing and damaging, the red dashed line indicating the frequency of such a cycle, as identified by the double arrow on Fig. 7(a).

Finally, we analyze the dependance of the simulated deformation on the time scale of observation using a temporal coarse-graining method (e.g., [START_REF] Rampal | Scaling properties of sea ice deformation from buoy dispersion analysis[END_REF]. Components of the strain rate at a given spatial scale are averaged over a time window of duration t to compute the mean total deformation < εtot > t . The window is centred on an arbitrary time t 0 and has a size t = 2n×(N ∆t) with n = 1, 2, 3, ... and with the smallest averaging time scale corresponding to the time of propagation of an elastic shear wave with speed c across the width L 2 of the domain. The chosen spatial averaging scale is that of the highest deformation rate, which as shown in section 25 is of L N . The domain is therefore divided in square boxes of equal size l = L N and the calculated deformation invariants are averaged over all available boxes. Figure 8(b) shows the total deformation rate < εtot > t as a function of the time of observation t (thick black line) averaged over 20 realizations of the coarse graining calculation (thin, coloured lines) centred on different t 0 for a simulation with N = 40.

Consistent with the localizing of the deformation and an intermittent process, < εtot > t decreases with increasing t over almost two orders of magnitudes of t. The observed scaling is however altered in two ways, which relate to the specific geometry, loading and boundary conditions used in the present simulations. First, as one main fault always dominates the deformation in the system, curves of < εtot > t are strongly modulated by a succession of peaks associated with the cycles of stress build-up and macro-rupture, the amplitude of which decreases with the scale of observation t. Second, at large t, the scaling assymptotes to a value corresponding to the prescribed forcing.

Simulations over larger systems using non-homogeneous surface forcing should allow for multiple macroscopic scale faults to be active simultaneously and hence to observe a clearer scaling of the simulated deformation over larger time spans.

Conclusions

In this paper we have presented a new mechanical framework suited for modelling the brittle behaviour of the sea ice cover [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF] while keeping a continuum description. A relaxation term for the internal stress is added to the original Elasto-Brittle constitutive relationship and both the linear and viscous components are coupled to a progressive damage mechanism to allow partitioning between the reversible and permanent deformations within the material based on its local level of damage.

Highly idealized simulations using forcing conditions homogeneous in both space and time show the Maxwell-EB model simulates a complex temporal and spatial evolution of the deformation patterns, in close agreement with observations of the Arctic sea ice cover. Anisotropy in the simulated damage and deformation fields arises naturally from elastic interactions, although the material's properties are fully isotropic at the element scale. The model also reproduces both the persistence of creeping leads and the activation of new leads with different shapes and orientations, in agreement with the observed deformation of sea ice [START_REF] Coon | Arctic Ice Dynamics Joint Experiment (AIDJEX) assumptions revisited and found inadequate[END_REF]. Analyses of the simulated damage and deformation fields reveal 1. a highly heterogeneous deformation, translating into a power law decrease of the deformation rate with increasing spatial scale. The associated exponent varies periodically: it is highest in the vicinity of macro-rupture events and decreases between events as the material partially heals. The disappearance after a few "spinup" rupture events of a cross-over scale at which deformation rates become independent of the scale of observation suggests that the Maxwell-EB model, including both damaging and healing processes, successfully reproduces a "marginally stable" state, as observed for Arctic sea ice.

2. an intermittent deformation, manifested by the highly asymmetric temporal evolution of the internal stress within the material, which shows a succession of slow build-ups and very rapid relaxation phases. This intermittency is supported by the existence of a temporal correlation in the rate of damage at all timescales below the material's characteristic healing time. A temporal scaling of the deformation rate is also obtained but due to the specific setup of the simulations analyzed here, it is modulated by the cycles of stress build-up and relaxation and its span is limited by the prescribed forcing.

Considering the highly idealized setup of the simulations analyzed here, these temporal and spatial scaling properties in the deformation fields cannot possibly be inherited from the prescribed forcing.

Instead, their emergence is a signature of the mechanical behaviour of the Maxwell-EB model itself.

The next logical step in the development of a Maxwell-EB sea ice rheology consists in analyzing the sensitivity of the simulated deformation and damage fields to the model parameters. In particular, the partitioning between the simulated brittle and creep-like behaviour as well as the degree of localization of the deformation [START_REF] Frederiksen | Numerical modelling of strain localisation during extension of the continental lithosphere[END_REF] might depend on the rate of decrease of the viscous relaxation time with increasing level of damage (parameter α) and on the characteristic time for healing and associated healing parameterization, all of which are poorly constrained in the case of the ice pack.

Further validation of the Maxwell-EB framework and the determination of the range of model parameters values suitable for sea ice call for a thorough comparison of the scaling properties of the simulated deformation rates with that estimated from the available ice buoy and RGPS data. Such analysis necessitates carrying numerical experiments over periods of several days to months and over realistic domains of regional to global scales. At these spatial and temporal scales, deformations within the sea ice cover become large. Hence advective processes cannot be neglected. As the Maxwell-EB rheology effectively reproduces very strong spatial gradients within the velocity, strain and stress fields, its use in large-deformation experiments requires the implementation of a robust advection scheme in order to limit diffusion and retain the strong localization of damage and deformation rates. The development of a numerical scheme for the the Maxwell-EB model that includes advection and is both efficient and practical in view of dynamic-thermodynamic and fully coupled ocean-sea ice-atmosphere simulations is underway.

Figure 1 .

 1 Figure 1. (a) Schematic representations of the Maxwell model for a continuum material with elastic (shear) modulus G and viscosity η. At time t, a stress is applied on the system. It is removed at time t + ∆t : the spring goes back to its initial position but the dashpot retains its deformation εη. (b) Loading-unloading paths for a material with initial elastic modulus E 0 in the linear-elastic (dotted), EB (dashed) and Maxwell-EB (solid lines) model. The black dot indicates the onset of damaging in the EB and Maxwell-EB models. Unlike the EB model, the Maxwell-EB model allows partitioning the total deformation into a permanent and an elastic contribution, indicated by the red arrows along the deformation axis. The diagram is not to scale in the context of modelling the lithosphere or sea ice : in these geomaterials, permanent deformations can become much greater than elastic deformations as damage events accumulate over time.

Figure 2 .

 2 Figure 2. Damage criterion of the Maxwell-EB model in the principal stresses plane (solid line) combining the Mohr-Coulomb and tensile stress criteria. The thick dashed line represents a biaxial compression truncation that closes the envelope but is not applied in the present model. Compression is taken positive and the dotted line indicates the σ1 = σ2 axis. Numbers indicate the states of (1) uniaxial tension, (2) biaxial tension and compression, (3) uniaxial compression and (4) biaxial compression and their location relative to the envelope.

  285

Figure 3 .

 3 Figure 3. Dependence of the apparent viscosity (η) the elastic modulus (E) and the relaxation time (λ) on the level of damage in the Maxwell-EB sea ice model. The image is a SPOT satellite aerial picture of a 59 km ×59 km portion of the Arctic sea ice cover centred around 80.18 N, 108.55 W.

Figure 4 .

 4 Figure 4. (a) Domain and boundary conditions for the uniaxial compression experiment.

  ) with D ũ D t = 0 and the 6 unknowns ũ (2 components), σ (3 components) and d. The model is made adimensional with respect to the length of the rectangular plate, L, the prescribed velocity U on the top boundary, the undamaged elastic modulus E 0 .

Figure 5

 5 Figure 5 shows the (adimensional) macroscopic stress, σ m (normal stress integrated on the upper boundary of the domain), as a function of the adimensional macroscopic strain, ε m , set by the prescribed displacement of the upper boundary, for these four simulations. The dotted line represents the damage rate (the number of damaged elements per model time step times their distance to the damage criterion, 1d crit ) for the simulation with N = 40. Inspection of the initial loading

  strictly linear-elastic behaviour at the initial stage of the experiment, as the material is initially undamaged, 2. a deviation from the linear-elastic behaviour after the onset of damage (marked by the red dot 1), indicative of macroscopic strain softening, with damage distributed homogeneously throughout the material (see Fig. 5 b1), 3. the formation of clusters of damaged elements, non-interacting at first, then joining along linear features. This stage is marked by a rapid increase in the number of damaged elements, 4. a sharp stress drop associated with the macroscopic failure of the sample and propagation of a main fault spanning the entire domain (see Fig. 5 b2).

  (a) represents the simulated macroscopic stress as a function of time (black dashed-dotted line) along with the corresponding damage rate (grey solid line) record for one realization of the uniaxial compression experiment with N = 40. Inspection of both temporal series reveals two types of mechanical behaviour of the Maxwell-EB material.

Figure 5 .

 5 Figure 5. (a) Macroscopic stress versus macroscopic strain (solid lines) for four uniaxial compression simulations with different spatial resolutions and damage rate (dashed grey line) for the simulation with N = 40. All simulations are initialized with the same values of mechanical parameters and cohesion field C defined at the lowest spatial resolution (N = 10). (b) Fields of the instantaneous damage (left panels) and of the order of magnitude of the total deformation rate (log 10 ( εtot), right panels) at the four different times indicates on Fig. (a) and for the four simulations (resolution increasing from top to bottom).

Figure 6 .Figure 7 .Figure 8 .

 678 Figure 6. (a) Macroscopic stress as a function of the macroscopic strain for one realization of the uniaxial compression experiment with N = 100. (b) Total deformation rate as a function of the spatial scale l (l = L 2n with 1 ≤ n ≤ N/2), normalized at the smallest scale L/N , at the five stages indicated on panel (a). (c) Zoom into panel (b) for the second, third and fourth stages. (d) Corresponding fields of the order of magnitude of the total deformation rate (log 10 ( εtot)) normalized by the maximum value of εtot.

  the Maxwell-EB model, the change in level of damage corresponding to a local damage event is determined as a function of the distance of the damaged model element to the yield criterion. Three important assumptions are made when calculating this distance, denoted d crit . The first is that the

deformation of each model element is conserved during a damaging event, i.e., at initiation, damage modifies only the local state of stress, not strains. The second is that, for a sufficiently small model time step ∆t, i.e. very small compared to the viscous relaxation time λ (see section 4.1), a negligible part of the stress is dissipated into viscous deformation. A third constraint is based on the fact that stresses outside the failure envelope are not physical because brittle failure would occur before the material could support them. Hence we consider that after being damaged, an element has its state of stress lying just on the failure envelope. With these assumptions, the following equality holds for each damaged element:

Table 1 .

 1 Model variables, parameters and domain dimensions for the uniaxial compression experiment.

	Parameters		Values
	Poisson's ratio	ν	0.3
	Internal friction coefficient	µ	0.7
	Ice density	ρ	900 kg m -3
	Shear wave propagation speed	c	500 ms -1
	Undamaged elastic modulus	E 0	2c 2 (1 + ν)ρ Pa
	Undamaged apparent viscosity	η 0	10 15 Pa s
	Minimum apparent viscosity	ηmin	10 4 Pa s
	Cohesion	C	(25 -50) • 10 3 Pa
	Damage parameter	α	4.0
	Undamaged relaxation time	λ 0	η 0 E 0 s
	Characteristic time for damage	t d	∆t s
	Characteristic time for healing	t h	10 5 s
	Dimensions of compression experiment		Values
	Length of the ice plate	L	200 • 10 3 m
	Prescribed velocity of forced edge	U	10 -3 ms -1
	Number of elements along short edge	N	10, 20, 40, 80, 100
	Mean model resolution	∆x	L 2N m
	Model time step	∆t	∆x c s
	Ice thickness	h	1 m
	Ice concentration	A	100%
	Variables		Non-dimensional equivalent
	Horizontal dimension	x	x = x L
	Time	t	t = tU L
	Ice velocity	u	ũ = u U
	Internal stress	σ	σ = σ E 0
	Level of damage	d	d = d
	Ice thickness	h	h = h H

Acknowledgements. The financial support of TOTAL EP RECHERCHE DEVELOPPEMENT is gratefully acknowledged. V. Dansereau has been supported by ANRT. A. Audibert-Hayet, E. Coche and K. Riska are thanked for valuable suggestions and support on this work. V.D. acknowledges support from the National Sciences and Engineering Research Council of Canada and from the Fonds Québécois de la Recherche sur la Nature et les Technologies.