A Maxwell-Elasto-Brittle rheology for sea ice modelling

GENERAL COMMENTS The literature discussion in the introduction is very much focused on sea ice modeling. Damage mechanics has been a vital topic within the glacier and ice sheet modeling community in recent years though, and – even though length, time and stress scales may differ by several orders of magnitude – the models used in this context are conceptually not fundamentally different from those used for sea ice. Particularly, the pioneering work by Pralong et al. (2006) (and several other articles by the same authors) might be important; but also later contributions, such as Duddu and Waisman (2013), could potentially be relevant.

main difference in the Maxwell-EB framework is that this decrement in damage, called dcrit, is not an arbitrary constant as in these previous progressive damage frameworks, which are based on a sub-iteration loop in which damage is allowed to evolve and stress to be redistributed from over to sub-critical elements until a steady-state is reached (i.e., until all states of stress lie within the yield envelope) before the model is incremented and loading resumes. In the present model, damage evolves in "real" model time. Hence if dcrit was set constant, the stress drop associated with damaging can be such that stresses will still lie outside of the yield envelope after one model time step. As overcritical stresses are not physical, here we make the logical assumption that the decrement of damage dcrit associated with a local damage event should be such that the value of the stress be at most equal to the local critical value.
Concerning anisotropy, we thank the referee for this important point, which was not detailed enough in the initial version of the manuscript. The key point is that, in the present model, an anisotropic damage formulation is not required to generate anisotropy of stress and strain 5elds. Indeed, in an elastic medium submitted to a non-perfectly isotropic loading (i.e., non-perfectly isotropic with respect to the domain geometry or the heterogeneity present in the material), the elastic kernel associated with a damaged "inclusion"  is anisotropic, hence is the redistribution of stresses. Therefore, the combination of (1) small-scale disorder, (2) damage mechanics in an elastic medium, and (3) this anisotropy of the elastic kernel itself is suf5cient to generate anisotropy up to very large space scales through successive elastic interactions between damaged elements. We now discuss this point in more details in the Results section and add a 5gure (6) that shows the anisotropic perturbation in the (Coulomb) stress 5eld (σ1 -qσ 2 ) generated when uniaxial compression is applied to a uniform rectangular plate with an isotropic, circular inclusion and, similarly, when disorder is introduced in the model at the element scale.
Finally, concerning the entropy principle and the healing parameterization, you are right that this remains an open question. As mentioned in section 4.3.2, we used the simplest possible parameterization for healing (using a constant healing rate) in the present uncoupled implementation of the Maxwell-EB model. We did not verify the validity of our approach with respect to entropy. In a dynamic-thermodynamic sea ice model, the rate of healing should logically be a function of the local difference between the temperature of the air near the surface of the ice and the freezing point of seawater below (see section 4.3.2). However we believe entropy is hard, perhaps impossible, to quantify and monitor in the context of an open, dynamically and thermodynamically coupled system such as the Arctic Ocean.

MINOR COMMENTS
p.2 l.6 (and many other places): "constitutive relationship:" Even though admittedly a constitutive relationship sounds very romantic, I guess the more common terminology in this context is "constitutive relation." We replace this formulation by "constitutive law". p.3 1.25 "physically consistent with observed behavior:" Is it physically consistent, or consistent with observed behavior, or both? In either case, that should be made clear. It is naturally both, as something that is physically consistent should necessarily be consistent with observation. We replace this sentence by: ''consistent with its observed mechanical behaviour''.
p. 5 l.9 "processing... method:" I am not sure whether "processing" is the right word for this. "Processing" is synonym of "treating" here, or "representing". We rephrase: thereby treating discontinuum mechanics with a continuum mechanics method.
p.5 l.25-p.6 l.2: Maybe make two sentences out of this in order to make it more readable. OK.
p.6 l.5: "material's velocity" rather "the velocity of the material" (as the material is not a person). OK.
p.6 l.13f: "the ration of ... and of the time": drop second "of" OK.
p.6 l.25: add commma after "However" OK. p.7 l.1: "i.e. .. forcing" what do you mean by this? Here we make the distinction between the part of the intermittency that is attributable to the forcing applied on the material and the part that is inherited from its mechanical behaviour. This is an important difference in the context of sea ice. The turbulent wind forcing itself exhibits some intermittency, which is "transmitted" to the ice cover. This part of the intermittency is therefore expected to be reproduced in most sea ice models. However, it was shown that the wind forcing is less intermittent than the deformation of sea ice, hence that the "extra" intermittency of sea ice deformation is attributable to its mechanical behaviour (Weiss, 2008;. It is this part of the intermittency of the deformation that cannot be adequately reproduced in sea ice model if elastic interactions are not accounted for and if the memory of elastic stresses is not adequately retained. p.7 l.5 what does "such a model" refer to? It refers to the previous sentence: to a rheological model that has "the capacity to distinguish between reversible and irreversible deformations" and that "allows a passage between the small/elastic and large/permanent deformations". We try improving the reading by adding "that" after "such a model". p.7 l.20: I do not quite understand why the use of a viscoelastic constitutive relation has to be justi5ed from rock mechanics? There is a large literature about viscoelasticity of ice. Once again, here we do not mean to model the classical bulk viscoelasticity of ice (see response to major comment above). Viscous creep is negligible in the deformation of sea ice, which is essentially brittle (e.g., Weiss et al., 2007). The apparent viscosity introduced here is intended to represent the slow relaxation of elastic stresses through permanent, large deformations of the damaged material: this is why we make a parallel with models of lithospheric faulting rather than with viscoelastic models for glacier ice.
p.8 l.26: Even if I feel bad about this self promotion, but a very similar model for glacier ice has been proposed by . See above.
p.9 l.23: I guess the strain rate tensor is the symmetric gradient of the velocity? Please state this explicitly. Furthermore, I think something should be said about how strain rates and strains are related (the notation suggests that the strain rate tensor is the rate of the strain tensor, which is generally not the case). Yes, the strain rate tensor is the rate of the strain tensor and hence is the symmetric gradient of the velocity. We clarify this point by de5ning the strain rate tensor in equation (1) as the symmetric gradient of the velocity. Large deformations in the Maxwell-EB model are accounted for by introducing the objective derivative of the stress tensor in the constitutive law (eqn. 1). To clarify this point further, we include the de5nition of the objective material derivative for the stress tensor after introducing the constitutive relation and develop the βa term, which accounts for the effects of rotation. It is important to recall however that all simulations presented in the paper are in the small-deformation regime (no advection and rotational effects neglected), hence these distinctions do not apply here.
Besides, we recognize that the notations used to describe the standard Maxwell model at the beginning of section 4.1were rather confusing, especially when referring to 5gure 1(a) and (b) to introduce these concepts. For instance, G was used for the elastic modulus and τ was used for the stress tensor in the text and in 5gure 1a, instead of E and σ in 5gure 1b. Hence we slightly reformulated the description of the standard Maxwell model and changed the notation in 5gure 1(a). We believe this improves the reading and presents the transition from the standard to the Maxwell-EB model more clearly.
p.10 l.20: Maybe write this as equation. This de5nition of the (adimensional) elastic stiffness matrix is quite standard. Hence we do not think this ought to be listed as a separate equation. But we reformulate this de5nition in index notation, which makes it somewhat easier to grasp.
p.11 l.4: The principal stresses are eigenvalues, not components. OK.
p.13 2nd paragraph: The physical interpretation of the damage variable crucially depends on how it affects the stress-strain relation. Therefore, I think this should be discussed right when de5ning the damage measure.
We are not sure we understand this comment. On the one hand, in the elastic regime, the effective stress is given by σ d =E 0 ε with 0<d ⩽1 , which is consistent with the typical interpretation of damage in progressive damage models for elastic materials. On the other hand, in the viscous regime, the effective stress is given by which does contrast with the classical de5nition of the effective stress due to the introduction of the damage parameter (exponent α). This viscous term does not represent a true viscous =ow of the bulk material, but instead an apparent viscosity aimed at slowly relaxing the elastic stresses within a damaged material. Hence the notion of effective stress is not entirely relevant here. As the parameter α is larger than 1, our formulation means that damage plays a stronger role on the apparent viscosity than the effective stress concept would do. We agree that this point might have been confusing, especially as we state about the damage variable "This variable is interpreted as a measure of sub-grid cell defects or crack density (Kemeny and Cook, 1986) and is allowed to evolve (...)" (p. 15, lines 14-16). Hence we reformulate this passage as "This variable is interpreted as a measure of sub-grid cell defects".
p.13 l.22f: Please de5ne h. Thanks for catching this: h stands for the ice thickness.
p.14 l.22: "Time steps" and "elements" are concepts of numerical methods. One should be careful with using them for motivating the governing equations of a model. In physics, space and time are not discrete.
Here the space and time discretizations do not "motivate" the governing equations of the model, although they need to be taken into consideration when writing the continuous form of the damage equation (that is, the part of the evolution equation for d pertaining to the damaging process). As already discussed in this section, this arises because of our treatment of the damage mechanism, which is similar to that of linear-elastic progressive damage models (e.g . The governing equation for this discrete process is written as a recursive relation (formerly on line 20, p. 15 and now numbered).
In the linear-elastic damage mechanics model on which the Maxwell-EB model is based, time does not enter the governing equations (e.g., Girard, 2010), hence the formulation of an evolution equation for this mechanism is not an issue and the damage equation is simply written in this recursive form. Here evolution equations are written and as the damage process is tied to the space and time discretizations, the order of the time scheme and spatial and temporal resolutions need to be taken into account when writing the recursive equation in continuous form. We agree that this might constitute a limitation of the current Maxwell-EB model, especially as it requires the use of an explicit time-stepping scheme for the damage evolution and the use of a small model time step.
p.15 l.1: Equation numbers missing? Please give a clear de5nition of epsilon (particularly in the context mentioned above, strains and strain rates). OK, we numbered this equation and de5ned epsilon, the strain tensor.
p.15 l.6: Again, principal stresses are not stress components (they do not have the transformation properties of tensor or vector components). OK.
p.17 l.10: It may be overly rigorous, but I think it is not a proper use of the Landau notation to write O(some constant number) We rewrite these formulations in words instead.
As the 5rst author is Canadian, Canadian english conventions were used (which combines that of American and British spelling).
p.18 l.23: "... are entirely de5ned by ... :" Either they are well-de5ned or not, there is nothing in between, so no need for this emphasis. Maybe better rephrase to "... only depend on" OK.
p.18 l.24: "the constitutive equation..." please specify which one (the constitutive framework consists of more than one equation). OK.
p.19 Eqs. 12-13: I think this regularization technique would be worth a more concise discussion: Is this purely a method to ensure convergence of the numerics? (In this case I'd suggest not to mix this with the physical governing equations). Or is this a conceptual problem? Finally, this is a known problem in continuum damage mechanics that the more damage is increasing the weaker its physical interpretation becomes. As mentioned on lines 13 to 15, p. 19, this formulation is introduced to ensure mathematical consistency, i.e., so that the constitutive equation be de5ned in the limit of d = 0. It is not meant to handle the physical interpretation of damage in the limit of a "completely damaged" material. We mentioned this point with the intend of being as rigorous as possible, as it is how the damage equation is written in our numerical scheme. However, as also pointed out, the use of a ''regularization technique'' for d = 0, as well of this speci5c technique (as opposed to introducing a minimum value on d instead of on η) has no impact on our results, as the level of damage d never reaches zero values (we set d 0 > 0 in all simulations). Hence we agree that introducing this level of precision within the description of the governing equations might not be necessary in the present paper. We suggest to remove this entire paragraph. To take care of any ambiguities, we specify the condition 0 < d 0 = d(t = 0)≤ 1 in the equations for E and η as a function of d (end of former page 18). We modify 5gure 3 accordingly and reformulate lines 13 and 14, p. 13, as "The level of damage is equal to 1 for an undamaged material and approaches the value of 0 in the case of a "completely damaged" material." p.19 l.17f: "we take this approach.... but it had really no impact on our results.... :" This sounds very sloppy (and the switch from present to past tense is somewhat random). This sentence was cut (see above).
p.20 l.1: "a 2-dimensional plate ... and a constant healing rate ..." What are the consequences of the two-dimensional plate assumption? Probably a simpler velocity 5eld? Apart from that, it sounds odd to me to squeeze those two assumptions (concerning completely different parts of the model) into one sentence. The 2-dimensional plate assumption here is equivalent to the plane-stress assumption, in the sense that we assume no stresses in the z direction (σ13, σ31, σ23, σ32, σ33 = 0). In this case, the (adimensional) elastic stiffness tensor K is de5ned such that for all symmetric tensor However, one difference is that, as in all regional and global sea ice models, the z-components of the deformation are not taken into account here (the sea ice momentum balance equation is 2-dimensional). Hence yes, the velocity 5eld is calculated in the horizontal plane (ux, uy) only. This assumption was listed independently of the constant healing rate assumption. We reformulate this sentence as: "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 plane stresses are assumed. A constant healing rate is used".
p.20 l.6: Not sure whether "assimilate" is the right word for this. It is now replaced by "represents".
p.20 l.9f: What does "internal stress" refer to, Cauchy stress? And shouldn't it rather be distributed "over" the depth? What is the advantage of keeping the entire stress tensor? Yes, the internal stress refers to the in-plane stress, σ. In virtually all sea ice models based on the continuum assumption, it is the common name for the stress arising from the sum of the mechanical interactions between ice =oes ( Weiss, 2013, Feltham, 2008, and many others). It is thereby distinguished from external stresses, which is the formulation used for the skin (air and ocean) drags per unit area on the ice cover.
Here, the entire stress tensor instead of the deviatoric part is used, as the simulated material is a compressible, elastic solid.
The ice concentration is the fraction of a model grid cell covered by ice, or in other words, the surface of ice (as opposed to open water) per unit area. In the present model, it has the same de5nition as in virtually all continuum sea ice models developed since 1979 (i.e., since the Hibler viscous-plastic model). Hence we do not believe a lengthly de5nition of this variable is necessary.
p.20 l.17: Rather make two sentences out of this. OK.
p.21 l.3: What is the reason not to write A = 1? OK.
p.21 l.6: "c*" is the * a typo? No, we did intent to use this symbol. This constant is usually denoted C . We use c* to distinguish it from the cohesion in the Maxwell-EB and EB frameworks.
p.22 l.11: "...parameters must evolve within..." so, they evolve while the model is running? This would change the dynamical equations. Otherwise, rather rephrase this. We rephrase: "In order for the Maxwell-EB model to represent the intended physics, the value of these parameters must lie within a certain range".
p.23 l.13: "This separation.... calculations:" Somewhat weird semantics and ambiguous syntax in this sentence. The content absolutely makes sense though. Maybe rephrase this (and split into two sentences). Agreed. The sentence has been split and reformulated as "This separation of scales ensures that elements cannot recover by healing more strength than they have lost by damaging within one time step. In the case of sea ice for instance, excess healing would effectively entail a net growth, or thickening, of the pack, a process that should instead be accounted for by thermodynamic balance considerations. However, considering the estimates of the speed of elastic waves and of the healing rate of leads aforementioned, the sea ice cover naturally meets the condition Th << Td".
p.25 l.5: "... transmitting the damage information within the material:" Not sure whether "within" is the appropriate preposition. We do think "within" works here.
p.25 l.6: I don't quite understand this logic... If the waves are not resolved, why should they be 5ltered out of the solution? Furthermore, replace "the model's solution" by "the solution of the model" (A model can be a person, just in this case for sure it is not ;-)) Agreed: this paragraph is somewhat confusing, and not strictly necessary. Hence we cut these two sentences. The main point here is expressed in the two remaining sentences of this paragraph, which is that inertial effects and the effect of the propagation of viscoelastic waves on the stress and deformation 5elds can be safely neglected in most sea ice implementations of the model.
p.26/27, description of the test geometry: I don't quite understand the boundary conditions. Is the velocity on the lower (short) edge set to (0, 0) or may there be a non-zero component in x-direction? Concerning the lateral boundary, in the text it says "no con5nement is applied on the lateral sides" (that is, the boundary may move freely?), whereas in the sketch (Fig. 4) it says u(x,y) = 0, thus the velocity is 5xed. Please clarify these ambiguities! Furthermore, if the lateral boundary is 5xed, what happens to the in=owing ice mass if the ice thickness is kept constant?
The 5gure that was included in this submitted version of the paper was not the good one, and yes, was somewhat confusing. The current, corrected version indicates that: (1) The velocity on the lower and upper edges of the plate is set to 0 in the y direction only : the x-component of the velocity can be non-zero.
(2) The velocity is strictly zero only on the lower-left corner of the plate (indicated by the triangle).
(3) No con5nement is applied on the lateral sides, hence these boundaries may move freely . However, as these smalldeformation simulations are run for a short time (until the cumulative applied deformation is of at most 10% -there was an error in the text [l.4, p. 27] so we corrected it-of the size of an average, single model element), the position of grid nodes is not updated in time (hence the FE spatial discretization is de5ned based on the initial mesh grid, i.e., is not updated in time). The lateral boundaries therefore do not "move" in the simulations presented here. We now mention this point in section 6. (4) As deformations are small, mass-conservation is not prescribed (see page 27, 5rst paragraph).
p.27 l.26f: "so that to be representative:" the conjunction "so that" should not be followed by an in5nitive. OK. p.28 l.1: drop the comma. OK.
p.28 l.10-26: Maybe move this to the literature discussion? Agreed, we move this to the Introduction section.
p.29 l.8: "... is its strong anisotropy" would this 5nding not call for the use of a tensor damage variable? Please explain why a scalar damage model is suf5cient. Again, a scalar damage model is suf5cient to generate large-scale anisotropy in an elastic medium. The kernel of elastic interactions and the resulting perturbation of the stress 5eld in such a material becomes anisotropic as soon as some spatial heterogeneity in its mechanical strength, or if the applied forcing is not purely isotropic with respect to this heterogeneity or to the geometry of the experiment (see response to major comment above).
p.30 l.19: "spatial distribution of the damage criteria" isn't the damage criterion the same everywhere? Only its parameters vary. If the parameters (e.g., C) of the damage criterion vary, then the damage threshold itself varies. Here the local damage threshold varies with the cohesion, C. We replace "criteria" by "threshold" here, to make this point clearer.
p.32 l.6: "we use the output of strain rate 5elds from simulations..." isn't it rather the output of simulations, not that of strain rate 5elds? We removed "the". p.32 l.24: "....and clusters in space ..." This seems somehow syntactically lost in the sentence? Unfortunately, we do not understand this comment. We replaced this formulation by "localizes again". p.34 l.9: "of the discrete failure events:" drop the article. OK.
p.34 l. 25: "... the of damage rate are anti-correlated" something missing here? Yes: increments of the damage rate.
p.36 l.7: "permanent deformations within the material" There are no deformations outside of the material, so no need to state this. OK.
p.36 l.9 "show the Maxwell-EB model simulates...." maybe make the beginning of the subordinate clause clear by using "that". Or even better split into two sentences. OK.
p.36 l.26: "internal stress within the material" why not just "the Cauchy stress"? Or even "the stress?" I guess there are neither external stresses, nor stresses outside of the material. See previous comment: "internal stress" is the term used extensively by the sea ice modelling community.
p.50, Fig. 6b,c: Excessive use of colored plots. Except for the yellow one, they all look more or less the same to me. I am probably not the only one who will have trouble distinguishing those plots: statistically, you can expect that about one out of ten male readers has a similar color vision de5ciency. So it probably makes sense to use dash patterns instead, or to add plot labels. Agreed: we switched to black and white and added markers wherever possible. p.51, Fig. 7a (and various other 5gures): Are the units dimensionless? If so, this should be made clear, e.g. with the tilde notation used in the text. Thank you for this comment: this was indeed not made clear since we forgot to mention that we dropped the "~" notation for all adimensional variables in the Results section. We now make this point clearer by adding a mention to this effect and by rearranging the order of a few sentence in section 6 (formerly section 4).

Answers to reviewer 2 GENERAL COMMENTS
The thoroughness and detail of the manuscript are commendable, although in some places the presentation was somewhat confusing as a result. I would recommend splitting the Introduction into separate Introduction and Background/Theory sections. A shorter Introduction section that clearly outlines the motivation for the new model, the most relevant context, and the general approach would bene5t the reader. As it is written, the Introduction currently is very long and detailed but in some places confusing in terms of both the writing and the relevance of this level of detail. In other places, relevant details seem to be left out, and it seems to be left to the reader to be familiar with all of the references in order to understand certain points. It seems to be an excellent review of the state of sea ice modelling, and demonstrates that the authors have a good grasp of the 5eld, but as a reader I found myself a little "lost in the weeds" at times. We followed your suggestion and divided the introduction into an Introduction and a Background section. In particular, we also integrated the suggestion of reviewer 1 of moving the discussion of heterogeneity, intermittency and anisotropy (at the beginning of the Results section) to the Background section. We believe this helps presenting our motivations in developing this new rheological framework more clearly, as we indeed aim to create a model that is able to reproduce these (3) allimportant characteristics of sea ice deformation.
Even in a fully viscous model of ice deformation, the stress balance can be non-local (for instance, in a glaciology context, viscous ice stream or ice shelf deformation is described by a stress balance that is inherently nonlocal, (e.g. MacAyeal, 1989)). You seem to be implying in several places that "long-range" interactions must come from elastic deformation (e.g. p 4, l 14). Am I reading this wrong, or are you indeed stating that long-range interactions can only be accounted for by elastic interactions? This is indeed not what we implied. You are totally right on the point that stress balance can be non-local in fully viscous models, in other words that viscous deformations can redistribute stresses in a non-local manner. In our case, elastic interactions by nature redistribute the stress over long distances and need to be accounted for in order for the heterogeneity of deformation to be adequately represented in models of the ice pack.
The results of the model are mesh-dependent, as damage localizes to the scale of an individual element. This is often viewed as a negative result, because the results of the model thus depend on how the user sets it up. However, many different approaches for nonlocal regularization of damage models have been proposed and adopted in a variety of settings (e.g. Bazant and Jirasek, 2002;Borstad and McClung, 2011). In these approaches, the stresses/strains/constitutive relation are computed by integrating over an intrinsic length scale related to the scale of heterogeneity of fracturing of the material. As long as the element size is smaller than this intrinsic length scale, the results of the model are independent of the resolution of the mesh. I think the authors should mention this type of approach in the manuscript, and discuss whether it might be feasible to produce mesh-objective results. Thank you very much for this comment: this is a very important point that we did not discuss, mainly to keep the paper as short as possible. Introducing an intrinsic length scale, i.e., a correlation length, ξ, for damage that is larger than the model grid cell would indeed allow the model solution to converge. However, in the context of sea ice modelling, introducing a correlation length ξ > Δx would not be physical, as the scale of natural heterogeneities (thermal cracks, brine pockets, etc) within the ice cover that serves as stress concentrators is much smaller than the typical model spatial resolution (on the order of a few to several kilometres). Furthermore, invariance of sea ice fracturing, as revealed from =oe sizes distributions, holds down to the meter scale . Hence here, disorder is introduced at the smallest available scale: that of the mesh element (through the 5eld of cohesion, C). We believe this is the most rigorous choice in terms of representing the physics behind the deformation of the ice cover, but as a result the model solution is mesh-dependent and does not converge locally. We agree that this point deserves a more extensive discussion in the paper. Hence we somewhat reformulated lines 15-22, page 11 and 1-7, page 12 to clarify the role of the spatial noise introduced in the model through the cohesion variable, C, and added a few line discussing non-local damage and convergence in section 5.1 (former lines 19-27, page 30, and 1-11, page 31).
The discussion of the damage formulation is a bit confusing. You mention in the text (p 14, l 24-26) that stresses outside the failure envelope are non-physical. However, unless I am missing something, you seem to be calculating your damage variable according to the distance beyond the failure envelope. It seems, then, that damage is a sort of constitutive postprocessing to "correct" the stress level such that it lies directly on the failure envelope. Some clari5cation is needed in the text on this point, since your schematic representation of the failure envelope in Figure 2 seems to contradict what you state in the text. Damage is based on the distance of the stress state outside of the envelope, and yet a stress state outside the envelope is non-physical... See response to reviewer 1's comment above. The damage of an element is calculated such that, just after a local damage event, the stress state of the damaged element lies on the envelope. I'm confused as to why a separate term for the ice concentration (A) is needed, as this seems at least partially redundant with damage. Why is it necessary to have both a concentration term and a damage term that modify both the elastic modulus and the viscosity? Isn't there some redundancy here, as a damaged fault/lead will necessarily have a reduced ice concentration? The ice concentration term seems to be simply added to the model at the very end of the model discussion, without much explanation. We indeed did not discuss sea ice concentration in length in the paper, since the de5nition of the ice concentration is the same in the Maxwell-EB model as in typical continuum sea ice models (it is the ice-covered surface by unit area, hence it varies between 0 and 1). Its equation of evolution is therefore also similar to that employed in these models. The coupling between the ice concentration and both mechanical parameters (see equations 17, 18) here was inspired from the coupling of A to the ice strength in compression parameter (P) suggested in the VP model of  as well as from the coupling between the elastic modulus E and A suggested by Girard et al. (2011) and Bouillon and Rampal (2015). In the case of the elastic modulus, E, this formulation represents the fact that when the ice concentration drops below about 90%, internal stresses become negligible and the ice is essentially in free drift. In the case of η, this dependence on ice concentration is consistent with the rapid decay of the apparent viscosity of granular media when decreasing their packing fraction from the close-packed limit . While concentration might indeed seem partially redundant with damage, we emphasize that the two variables represent different things, as the ice pack might be densely fractured, hence not withstand large stresses, but still retain a high concentration (for instance under convergent motion). Nevertheless, we agree that this parameterization, which we chose here so that to be consistent with the approach taken in previous sea ice models, could eventually be re5ned.

SPECIFIC COMMENTS
25: "or" OK. p 4, l 6-9: how have these VP hypotheses been found inconsisent? Can you summarize these for the reader? The VP hypotheses have been found inconsistent in many respects, which are discussed in length in the studies cited here (Weiss et al., 2007, see lines 8-9, page 3). Although discussing these inconsistencies could be pertinent here, it would also make the introduction signi5cantly longer. Hence we chose to add only a short list of the main points at the end of this sentence (line 9, page 3) with the references to the relevant studies.
check English spelling throughout the document, a number of words are misspelled (looses instead of loses, euqation, it's instead of its, dependance, assymptotes,...) OK.
p 10, l 20: there is some inconsistency in the text formatting of the different versions of "I" for the identity tensors We reformulated the de5nition of K using index notation. p 13, l 22-23: I was confused about what h is here Thank you for catching this: h refers to the ice thickness and its de5nition in now de5ned. p 19, l 18: might there be other contexts or model setups (e.g. realistic domains) for which the minimum value of d might come into play? The results of a damage model can be quite sensitive to this choice. Problems, both conceptual and numerical, arise when the value of d is zero, which according to the evolution equation for d (line 20, p. 15 or equation 9) occurs if, and only if, the initial value of d (d 0 ) is set to zero. In the uniaxial experiment presented here, this problem does not arise since the experiment is started from a homogeneous, undamaged state (d 0 = 1). A regularization technique, such as the one presented in this paragraph, can be used to avoid numerical problems when d 0 is allowed to take the value of 0. Alternatively, a cutoff minimum value for d could be used. As this section brings essentially useless complexity in the context of the simulations presented in this paper, we decided to remove it in the revised manuscript (see response to reviewer 1's comments). p 20, l 11-13: Some motivation or explanation is needed for why you choose to write the momentum equation in terms of internal stress rather than the vertically integrated stress tensor, especially if you are departing from what is more commonly done in the sea ice modelling community. We chose to incorporate the vertical integration of the internal in-plane stress in the momentum equation instead of de5ning σ as the vertically-integrated stress tensor to avoid confusion when evaluating the distance to the damage criterion, i.e., when comparing the local state of stress σ to the critical tensile stress σt and critical stress with respect to the Mohr-Coulomb criterion σc, both de5ned in Pascals (Nm -2 ) as the cohesion C of the material. De5ning σ as the vertically-integrated stress tensor would indeed necessitate rede5ning σt and σc in terms of an effective cohesion, C x h. This approach is taken in the recently developed NeXtSIM model, based on the Elasto-Brittle rheology of Girard et al. (2011), and as this point was explained in the paper by Bouillon et al. (2015), we do not feel it needs to be discussed in length here, but we do agree that a clear reference to this paper is needed. Hence we reformulate the lines 10 to 13 as follows to explain this point more clearly: "We assume the internal stress to be homogeneously distributed over the depth h and write the momentum equation in terms of the internal stress rather than the vertically integrated stress tensor more commonly used in the sea ice modelling community. This approach was also taken in the Elasto-Brittle model of Bouillon et al. (NeXtSIM, 2015), as it allows a direct comparison between the local state of stress and the critical stress (here σt or σc) when estimating the distance to the damage criterion." p 23, l 6-9: why not perform a sensitivity analysis as you describe then? We did perform sensitivity analyses on the relative values of td and Δt in the Maxwell-EB model. The results of such analyses however, depends on the applied forcing, domain geometry, etc, and hence is speci5c to each numerical experiment. Presenting the details of these analyses is therefore beyond the scope of the present paper. The remark of lines 6-9 was instead meant as a warning one should be careful and perhaps carry sensitivity analyses when using a model time step that is larger than the prescribed time of propagation of damage ( td). To ensure numerical stability and allow for highest resolution of the elastic interactions in the Maxwell-EB simulations presented in this paper, we set td = Δt. As this also ensures the one to one correspondence between the progressive damage mechanism as described by the recursive relation of line 20, page 15 and the continuous form of equation 20, this appears to be the most logical choice. Hence to avoid giving the impression thatΔt could be set arbitrarily larger than td in the model (at the cost of a sensitivity analysis) without any loss of physical rigor, we removed this sentence and slightly reformulated this paragraph. p 26, l 5: you previously described the inertial term as being negligible, so why is it here? some clari5cation is needed. We replace "inertial term" by the "time derivative" of the stress tensor.
The 5rst part of the Results section is not really results, but background. It is now included in the background section. p 29, l 23-24: this would be helpful to state also in the 5gure caption OK.
p 30, l 26-27: well, the localization scale is the element scale, so the choice of resolution dictates the localization of damage We totally agree with this remark, and this is what we mean by "there is no physical scale associated with the localization of damage. Through elastic interactions, damage and deformation tend to localize at the 5nest scale (the mesh element)". p 31, l 6-11: but you didn't introduce disorder initially, so you cannot claim this here. Perhaps we do not fully understand this comment, but we indeed introduced disorder (in the damage threshold) through the 5eld of cohesion, set at the beginning of each simulation. We try rephrasing this line as "the initial disorder introduced in the model".
p 31, l 27: "...has already been investigated in depth..." is another example of the reliance on the reader to be familiar with all of the literature you are citing. It would be more helpful to summarize the 5ndings. What did these investigations 5nd? Thank you for this comment. Here we meant that damage models based on a linear-elastic constitutive law and without healing have been demonstrated to reproduce a highly heterogeneous deformation (the 5nding that is most relevant for our purpose). However, as these frameworks neither include a healing mechanism nor a slow relaxation of elastic stresses, the post macro-failure behaviour of these models is physically inconsistent, and only the path to the 5rst rupture was analyzed. Here we seek to establish if the Maxwell-EB model, based on a viscoelastic constitutive law, has the same capability of reproducing a highly heterogeneous deformation in a partially damaged material (over subsequent healingfracturing cycles). We hence reformulate these sentences to make this point clearer. p 34, l 25: some word has been omitted here Yes, "increments", thank your for catching this. Figure 5: it doesn't look like the elements are getting smaller from the top row to the bottom row of panel (b), but isn't the resolution supposed to be getting 5ner moving down in the 5gure? Also, the damage rate axis in panel (a) is missing a numerical scale other than the zero. Yes, resolution is increasing from top to bottom on this 5gure. However, you might get the impression that the resolution is the same between the experiments, especially at the beginning of the simulation (panels 1 and 2,) because the initial 5eld of cohesion is the same in all simulations. That is, we use the 5eld of C prescribed in the N = 10 simulation (lowest resolution, top row). This 5eld is interpolated onto the 5ner resolution grids (see lines 16 to 19, p. 29) in the other three simulations. Perhaps this point was not made clear enough, hence we slightly reformulate the description of these experiments, at the beginning of section 5.1. As for the right y-axis on panel (a), thank your for catching this. Figure 6 is presented in the discussion of heterogeneity, the dependence on the spatial scale of observation. It's still not clear to me how this is represented in the 5gure, which only shows one realization of the experiment at one resolution. Indeed, this 5gure shows one realization of the experiment. However, repeating the same scaling analysis for other realizations of the experiment gave very similar results. In this sense, averaging over several realizations did not give more insight. This is what we meant by "Repeating the procedure for subsequent healing and damaging cycles and for multiple realizations of the experiment initialized with different cohesion 5elds 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., β = 0.15 ± 0.02 ; Girard et al., 2010a)." The 5gure below shows for instance the result of averaging ⟨ε tot ⟩ l computed as a function of the spatial scale, at 5ve equidistant stages (as indicated on 5gure 7a, formerly 5gure 6a) between the minimum in macroscopic stress that follows the propagation of a fault (red) and the maximum that precedes the next macro-rupture (purple), over the second cycle of stress build-up and macro-rupture for 5 different realizations of the uniaxial compression experiment with a resolution of N = 100. The total deformation rate still shows a clear power law decrease with increasing spatial scale, with β largest for the post-macro-rupture (red) and pre-macro-rupture (purple) stages. The averaging could alternatively be done over multiple stress build-up/macro-rupture cycles of a single experiment. However, we believe that giving an appropriate description of the method for partitioning of the results in different stress build-up/macro-rupture cycles and averaging ⟨ε tot ⟩ l over these multiple cycles so that to present an "average" 5gure such as the one below instead of 5gure 7 would make the reading of this part of the Results section more dense than insightful.
Moreover, we do not believe that varying the resolution would be relevant here. Increasing the resolution could indeed be interesting, in the sense that it would allow performing the scaling analysis over a larger range of space scales. Obviously it would also be computationally more expensive. Here, the analysis and the power law obtained already spans almost 2 orders of magnitude in l.

Abstract
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
A close comparison can be made between the deformation of sea ice and that of the Earth crust, 90 in which brittle fracturing and Coulomb stress redistribution also take place and for which scaling properties have been recognized for years (Kagan and Knopoff, 1980;Kagan, 1991;Kagan and Jackson, 1991;King et al., 1994;Turcotte, 1992;Stein, 1999). Recently, Marsan and Weiss (2010) 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 95 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 (Kagan, 1991;Kagan and Jackson, 1991). 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.

100
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 (Cowie et al., , 1995. Building on similar linear-elastic laws and introducing some strain softening at the micro scale, the failure 105 model of Tang (1997) 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 :::::: treating :::::::::::: discontinuum ::::::::: mechanics :::: with ::::::::: continuum ::::::::: mechanics :::::::: methods. An analogous approach based on local damage evolution was also taken by , who combined 110 a linear-elastic constitutive relationship for a continuum ::: law ::: for : a :::::::::: continuous 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 115 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. (2011) 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 120 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., 2011).
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 125 the fracturing of the pack, but also the permanent deformations occurring once it is fragmented and undamaged ::::::: fractured ::: and :::::::::: fragmented. :::::: When ice floes move relative to each other along open leads, as these much larger deformations set the advective processes and overall drift pattern of the ice cover.
This last point is an important and intrinsic limitation of the EB framework, since the linear-elastic 130 constitutive law does not allow solving for the elastic (reversible) and permanent deformations of the simulated material separately. :::::::: separately. : Hence to estimate the material's velocity :::::: velocity :: of ::: the ::::::::: simulated ::::::: material, 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 dam-135 aged 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. (2011).
(2) Alternatively, if all of the resulting deformation is considered permanent, the material keeps its does represent advective processes over the Arctic (Bouillon and Rampal, 2015). 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., ::: that :: is, ::: the :::: part :: of ::: the ::::::::::: intermittency :::: that :: is not 155 directly inhered from the wind forcing (Rampal et al., 2009). In order to estimate adequate drift velocities, a suitable rheological model must therefore have the capacity :::::::: capability : to distinguish between reversible and irreversible deformations.
The goal of this work is to develop such a model allowing ::: that :::::: allows a passage between the small/elastic and large/permanent deformations and with the capability of damage mechanics models of bulk ice (Duval et al., 1983), 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 170 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 :: by : 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 175 transient, aseismic manner (Scholz, 2002;Gratier et al., 2014;Cakir et al., 2012;Cetin et al., 2014).
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 (Cakir et al., 2012;Cetin et al., 2014). A further justification of the introduction of such pseudo-viscosity comes from the rheology of granular media. As sea ice along leads (see Fig.   180 3), rocks along active faults are highly fragmented. Sheared granular media flow in a viscous manner when inertial effects can be neglected (Jop et al., 2006) with an apparent viscosity diverging as the packing fraction approaches the close-packed limit . This last point will justify the dependence of our apparent viscosity on sea ice concentration.
For metals and rocks, the MC theory was shown to be defective in the case of tension (Paul,305 1961), 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 310 N ⇠ 50 kPa (Weiss et al., 2007). 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 (Paul, 1961), as shown on Fig. 2. The tensile strength cutoff therefore takes the form: 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 (Weiss et al., 2007) and in the lab (Schulson, 2006a 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 (Weiss et al., 2007). Besides, imposing a truncation was shown to have little impact on the 325 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  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 pa-330 rameters 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 devolving between . :::: The :::: level ::: of :::::: damage :: is ::::: equal ::: to 1 (undamaged ) and :: for ::: an ::::::::: undamaged ::::::: material :::: and ::::::::: approaches ::: the ::::: value ::: of 0 ( :: in ::: the :::: case :: of :: a "completely dam-335 aged" material). This variable is interpreted as a measure of sub-grid cell defects or crack density 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 strength- the ::: ice :::::::: thickness) and effective apparent viscosity (⌘ ⇥ h) of the ice, healing is distinguished from pure thermodynamic growth or dynamically-driven thickness redistribution (e.g., Rothrock, 1975) in that it applies only where and when the material has been damaged. It therefore allows d, E and 345 ⌘ 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.

350
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., Tang, 1997;Hamiel et al., 2004, and others). Because of the long-range interactions within the elastic medium, local drops in E imply a stress redistribution 355 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 Girard et al., 2010). As the elastic perturbation generated by such events is anisotropic , this propagation mechanism naturally leads to the emergence of both spatial het-erogeneity and anisotropy in the stress and strain fields, i.e., to the formation of linear-like faults (see 360 section 7).
In 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 365 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 5.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 370 of stress lying just on the failure envelope. With these assumptions, the following equality holds for each damaged element: where the (12) 375 ::::: where : ✏ :: is ::: the ::::: strain :::::: tensor ::: and ::: the : superscript 0 denotes the post-damage state of deformation and stress. In terms of the principal stress components :::::: stresses, the change in level of damage of a given element is given by which implies that as the level of damage varies, all stress components vary in the same proportions.

380
Hence the state of stress 0 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 0 , 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: 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 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:

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 400 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. (2010), which follows parameterizations of the vertical growth of sea ice (Maykut, 1986). 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, 405 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 : Combining both the damaging and healing mechanisms (Eq. 14, 16 and 17), the complete evolution equation for d is @d @t Although the two processes apply simultaneously on the level of damage in the model, they are 410 inherently distinct. On the one hand, damaging is a discrete threshold mechanism applying only where and when :::: when ::: and :::::: where the state of stress becomes overcritical. As mentioned in sections 4.2 and 4.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 (Marsan et al., 2011). For spatial resolutions between 415 that of current global climate and high resolution regional sea ice models ( x = 1 :::::: x ⇡ 1 : to 100 km), the characteristic time for damaging :::::: damage :::::::: evolution, t d therefore varies between O(1) and O(10 2 ) s : a :::: few ::::::: seconds :: to :: a ::: few ::::::::: hundreds :: of ::::::: seconds. 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 420 within an opening of 10 cm under atmospheric temperatures of T a = 15 C is of O(100) hours or O(10 5 ) seconds : a ::::::: hundred ::::: hours : (between 10 5 and 10 6 seconds, Petrich et al., 2007). 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.

425
The coupling between the Maxwell-EB constitutive relationship ::: law : 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 430 and a very large apparent viscosity ⌘ 0 . In this case, the viscous term in (4) is negligible and a linear-elastic constitutive relationship ::: law : 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 435 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 440 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 (Agnon and Lyakhovsky, 1995;Schapery, 1999), we use the simplest parameterization 445 and set such that , with :::::::::::::: 0 < d(t = 0)  1 :::: such ::: that : 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 ⌘ 455 and E are entirely defined by their initial value , a constant, and by ::::: depend :::: only ::: on :::: their :::::::::: undamaged :::: value :::: and :: on : the level of damage variable d. 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 Substituting for ⌘ in the expression for the relaxation time, the elastic modulus then becomes Using such a cutoff on ⌘, the elastic term in the Maxwell-EB constitutive 465 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 470 imposing a minimum value of d, but it had really no impact on our results since in the simulations presented here d > d(E min ) at all times.
3. The constitutive relationship ::: law : (4) with :::::::::::::::: , :::::::::::::::: , where f 1 and f 2 represent the functional dependance ::::::::: dependence : on the level of damage of the ice d, given by (19) and (21) respectively. The exponential function of the ice concentration allows the internal stress term to be maximal when A = 100% ::::: A = 1 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 . Here the non-500 dimensional parameter c⇤ :: 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 (18) with the damage criterion defined by Eq. (6) and Eq. (9) and q, c and t given by (7), (8), (10) in terms of the cohesion variable C and of 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 simulated velocity field. Table 1 lists all model variables and parameters.

Characteristic numbers and times
where The temporal resolution that is optimal in terms of capturing all elastic interactions within the simulated material ::: and ::::::: ensuring ::::::::: numerical ::::::: stability is therefore t = t d . In the model experiments :: all ::::::::: simulations : 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 545 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 . :: In ::: the :::: case :: of ::: sea ::: ice ::: for ::::::: instance, : excess healing would effectively entail a net growthof the material, :: or :::::::::: thickening, :: of ::: the :::: pack, a process that is not intended by this parameterization and should instead be accounted for by thermodynamic , balance calculations.

We
The Weissenberg number, We, defined as the dimensionless product of the viscous relaxation time 555 for the stress and of time T = L U characterizing the deformation process: 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 560 level of damage as We = We 0 d ↵ 1 with We 0 , its maximum value.
As viscous dissipation should be insignificant over :: in 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 ::: law (4) 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 565 second, due to a numerical scheme failure known in the field of viscoelastic flow computations as the high Weissenberg number problem (Keunings, 1986;Kupferman, 2004, 2005;Saramito, 2014). 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 (5). With We 0 (or equivalently, 0 ) too low, simulations can run for a time t ⇠ 0 and unphysical viscous dissipa-570 tion can occur over :: in undamaged parts of the simulated material. To get round this problem, :::: This :::: issue ::: can ::: be :::: dealt :::: with ::: by ::::::::: multiplying : the viscous term in the Maxwell constitutive relationship can be multiplied ::: law 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 575 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 , : is : 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 585 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 590 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 595 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 600 time becomes very small whatever the damage level (see section 4.3.3). This means that elastic deformations are almost immediately dissipated after damaging, that is, the model becomes purely elasto-plastic. ::: The ::::::::: sensitivity :: of ::: the :::::: model :: to ::: the ::::: value ::: of ::: this ::::::::: parameter :::: was :::: kept ::: for : a :::::::: separate ::::: paper, ::::: hence :: ↵ :: is ::: not ::::: varied :::: here. : For the experiments presented here : in ::::::: section : 7, we find that ↵ = 4 allows representing both the brittle behaviour and the relaxation of the internal stress within : in : a 605 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).

Numerical scheme and experiments
The objective time derivative for :: of : the Cauchy stress in the Maxwell-EB constitutive relationship 610 ::: law (4) is composed of an inertial term : a :::: time :::::::: derivative, 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 :::: time :::::::: derivative : while neglecting the advection and a terms allows retaining a 615 Lagrangian scheme, similar to the original EB model (Girard et al., 2011). 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 620 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.

625
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 set small enough to ensure a low driving rate (i.e., slow compared to time scale of damage propagation, ).

630
In the present implementation, the model is not yet coupled to a thermodynamic component, hence 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 ::: the ::::::::: cumulative ::::::::::: deformation :: is ::::::  10% :: of ::: the :::: size ::: of : a :::::: single ::::: model ::::::: element), dynamicsinduced variations (through convergence-divergence) of the ice thickness and concentration are not 635 accounted for and hence the mechanical parameters E, ⌘ and C are not yet coupled to h or A.

Results
In this section we analyze the mechanical behaviour of the Maxwell-EB model. In particular, we evaluate its capacity :::::::: capability to reproduce the main characteristics of sea ice deformation , which are its spatial -::: its heterogeneity, intermittency and anisotropy , following the methodology developed 675 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, 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 680 total deformation rate <" tot > l invariant and the corresponding averaging scale l of the form 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.

685
This coarse-graining calculation was later extended to ice buoy data 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, 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 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 -:::::: which :::: have :::: been ::::::: recently used as benchmarks to validate (or invalidate) sea ice models (e.g., Girard et al., , 2010Bouillon 695 and Rampal, 2015).
7.1 Spatial resolution, convergence and dependence on the initial conditions ::::::::: anisotropy In a first time, we analyze the overall, macroscopic behaviour of the Maxwell-EB model , its convergence properties and the dependance :: as :::: well :: as ::: the ::::::::::: convergence ::: and :::::::::: dependence : of the so-705 lution 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 simulationsas well as : . ::: So : is : the field of cohesion, which is defined at the lowest resolution (N = 10) and interpolated onto the higher resolution mesh grids.

810
To quantify the heterogeneity ::::: spatial :::::::::: localization of the simulated deformation, we follow Marsan et al. (2004) 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. (2010). 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 815 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   7 (a and b) shows the total deformation rate <" tot > l as a function of the space scale l at 5 equallyspaced 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 820 next macro-rupture, as indicated in Fig. 7(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.
(1) 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 825 mechanical strength by healing. Deformation 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 ::::::: localizes ::::: again. :::: The exponent ::::::: therefore re-increases towards its post macro-rupture value (blue and purple curves).

830
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., 2010, = 0.15 ± 0.02).
However, an important difference between the present results and that of Girard et al. (2010) is the 835 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 :: its : 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 840 RGPS-derived deformation rates of Marsan et al. (2004) and Stern and Lindsay (2009), 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 8 First, the evolution of the macroscopic stress is clearly characterized by cycles of slow stress 850 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. 8a).
Second, as identified on the same time series, some periods (e.g., the interval delimited by the 855 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. 8b) 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 860 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 865 realizations of the compression experiment initialized with different fields of C over domains with N = 40. Fig. 9(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 un-870 correlated 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 instan-875 taneous 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 :::::: damage :::: rate ::::::::: increments : 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. Finally, we analyze the dependance ::::::::: dependence : of the simulated deformation on the time scale 885 of observation using a temporal coarse-graining method (e.g., . 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 ::: rate : <" 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 890 chosen spatial averaging scale is that of the highest deformation rate, which as shown in section 1 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 9(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 :::: grey lines) centred on different t 0 for a simula-895 tion 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 900 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 ::::::::: asymptotes to a value corresponding to the prescribed forcing. Simulations over larger systems using nonhomogeneous 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 905 spans.

Conclusions
In this paper we have presented a new mechanical framework suited for modelling the brittle behaviour of the sea ice cover (Weiss et al., 2007) while keeping a continuum description. A relaxation term for the internal stress is added to the original Elasto-Brittle constitutive relationship ::: law : and 910 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 ::::: based :: on ::: the local level of damage :: of ::: the ::::::: material.
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 pat-915 terns, in close agreement with observations of the Arctic sea ice cover. Anisotropy in the simulated damage and deformation fields arises naturally from :: the :::::::::: small-scale ::::::: disorder ::: and 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 ori-entations, in agreement with the observed deformation of sea ice . Analyses of 920 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 de-925 formation 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 930 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 :: of ::: the ::::::: material. 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.

935
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, 940 the partitioning between the simulated brittle and creep-like behaviour as well as the degree of localization of the deformation (Frederiksen and Braun, 2001) 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.

945
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 ::: out numerical experiments over periods of several days to months and over realistic domains of regional to global scales. At these spatial and temporal scales, deforma-950 tions 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 955 advection and is both efficient and practical in view of dynamic-thermodynamic and fully coupled ocean-sea ice-atmosphere simulations is underway.