M Philippe Lattes

Anne-Marie Côté

Elizabeth Côté

Collège Rivier

Pascaline Daniel Léveillé

Pauline Quentin

Yoann

Nelly

Stéphane

Mathias

David

Jonathan Lucas

Fred Bruno

Xue Meng

Keywords: sea ice deformation, Maxwell-elasto-brittle rheology, brittle mechanics, elastic interactions, viscous stress relaxation, heterogeneity, intermittency, anisotropy, finite elements and discontinuous Galerkin methods déformation de la banquise, rhéologie Maxwell-élasto-fragile, mécanique fragile, interactions élastiques, relaxation visqueuse des contraintes, hétérogénéité, intermittence, anisotropie, méthodes de calcul

In recent years, analyses of available ice buoy and satellite data have revealed the strong heterogeneity and intermittency of the deformation of sea ice and have demonstrated that the viscous-plastic rheology widely used in current climate models and operational modelling platforms does not simulate adequately the drift, deformation and mechanical stresses within the ice pack.

A new alternative rheological framework named Maxwell-Elasto-Brittle (Maxwell-EB) is therefore developed here in the view of reproducing more accurately the drift and deformation of the ice cover in continuum sea ice models at regional to global scales. The model builds on an elasto-brittle framework used for ice and rocks. A viscous-like relaxation term is added to a linear-elastic constitutive law together with an e↵ective viscosity that evolves with the local level of damage of the material, like its elastic modulus. This framework allows for part of the internal stress to dissipate in large, permanent deformations along the faults/leads once the material is highly damaged while retaining the memory of small, elastic deformations over undamaged areas. A healing mechanism is also introduced, counterbalancing the e↵ects of damaging over large time scales. The numerical scheme for the Maxwell-EB model is based on finite elements and variational methods. The equations of motion are cast in the Eulerian frame and discontinuous Galerkin methods are implemented to handle advective processes. Idealized simulations without advection are first presented. These demonstrate that the Maxwell-EB rheological framework reproduces the main characteristics of sea ice mechanics and deformation : the strain localization, the anisotropy and intermittency of deformation and the associated scaling laws. The successful representation of these properties translates into very large gradients within all simulated fields. Idealized numerical experiments are conducted to evaluate the amount of numerical di↵usion associated with the advection of these extreme gradients in the model and to investigate other limitations of the numerical scheme. First large-deformation simulations are carried in the context of a Couette flow experiment, which allow a comparison with the result of a similar laboratory experiment performed on fresh-water ice. The model reproduces part of the mechanical behaviour observed in the laboratory. Comparison of the numerical and experimental results allow identifying some numerical and physical limitations of the model in the context of large-deformation and laboratory-scale simulations. Finally, the Maxwell-EB framework is implemented in the context of modelling the drift and deformation of sea ice on geophysical scales. Idealized simulations of the flow of sea ice through a narrow channel are presented. The model simulates the propagation of damage along arch-like features and successfully reproduces the formation of stable ice bridges.

A Maxwell-Elasto-Brittle model for the drift and deformation of sea ice

Résumé

Depuis quelques années, des analyses statistiques de données satellitales et de bouées dérivantes ont révélé le caractère hautement hétérogène et intermittent de la déformation de la banquise Arctique, démontrant de ce fait que le schéma rhéologique visco-plastique traditionellement utilisé en modélisation climatique et opérationnelle ne simule pas adéquatement le comportement dynamique des glaces ainsi que les e↵orts mécaniques en leur sein.

Un cadre rhéologique alternatif, baptisé Maxwell-Élasto-Fragile (Maxwell-EB) est donc ici développé dans le but de reproduire correctement la dérive et la déformation des glaces dans les modèles continus de la banquise à l'échelle régionale et globale. Le modèle se base en partie sur un cadre de modélisation élasto-fragile utilisé pour les roches et la glace. Un terme de relaxation visqueuse est ajouté à la relation constitutive d'élasticité linéaire avec une viscosité e↵ective, ou "apparente", laquelle évolue en fonction du niveau d'endommagement local du matériel simulé, comme son module d'élasticité. Ce cadre rhéologique permet la dissipation partielle des contraintes internes par le biais de déformations permanentes, possiblement grandes, le long de failles (ou "leads") lorsque le matériel est fortement endommagé ainsi que la conservation de la mémoire des contraintes associées aux déformations élastiques dans les zones où le matériel reste relativement peu endommagé. Un mécanisme de cicatrisation permettant de représenter le regel le long des failles formées est également inclus.

The schéma numérique du modèle Maxwell-EB est basé sur des méthodes de calcul variationnel et par éléments finis. Une représentation Eulérienne des équations du mouvement est utilisée et des méthodes dites Galerkin discontinues sont implémentées pour le traitement des processus d'advection.

Une première série de simulations idéalisées et sans advection est présentée, lesquelles démontrent que la rhéologie Maxwell-Élasto-Fragile reproduit les caractéristiques principales du comportement mécanique de la banquise, c'est-à-dire la localisation spatiale, l'anisotropie et l'intermittence de sa déformation ainsi que les lois d'échelle qui en découlent. La représentation adéquate de ces propriétés se traduit par la présence de très forts gradients au sein des champs de contrainte, de déformation et du niveau d'endommagement simulés par le modèle. Des tests visant à évaluer la di↵usion numérique découlant de l'advection de ces gradients extrêmes ainsi qu'à identifier certaines contraintes numériques du modèle sont ensuite présentés. De premières simulations en grandes déformations, incluant les processus d'advection, sont réalisées, lesquelles permettent une comparaison aux résultats d'une expérience de Couette annulaire sur de la glace fabriquée en laboratoire. Le modèle reproduit en partie le comportement mécanique observé. Par ailleurs, les di↵érences entre les résultats des simulations et ceux obtenus en laboratoire permettent d'identifier certaines limitations, numériques et physiques, du modèle en grandes déformations. Finalement, le modèle rhéologique est utilisé pour modéliser la dérive et la déformation des glaces à l'échelle de la banquise Arctique. Des simulations idéalisées de l'écoulement de glace dans un chenal étroit sont présentées. Le modèle simule une propagation localisée de l'endommagement, définissant des failles en forme d'arche, et la formation de ponts de glace stables. The ice that covers the polar oceans moves and deforms under the action of the winds and ocean currents.

Table of Contents

Making reliable predictions of the drift and deformation of this thin floating sheet of ice is becoming crucial nowadays for: (1) forecasting the opening of shipping routes across the Arctic, (2) evaluating mechanical constraints on o↵shore structures and ships and, at larger scales, (3) estimating the future evolution of both its summer and winter extent in the Arctic and Antarctic to anticipate its short to long-term, regional to global impacts on climate.

Early aerial pictures of the Arctic ice cover have revealed a highly heterogeneous, "densely fractured material" [START_REF] Coon | Modeling the pack ice as an elastic-platic material[END_REF], divided in plates, called "floes", of sizes ranging from a few meters to several kilometres (see figure 1.1). The relative motion of these floes translates into the opening of cracks, joining along larger features called "leads", the shearing deformation along opened cracks, the closing of leads and the formation of pressure ridges [START_REF] Coon | Modeling the pack ice as an elastic-platic material[END_REF][START_REF] Feltham | Sea ice rheology[END_REF]Weiss, 2013, and others). In more recent years, satellite remote sensing data such as the RADARSAT Geophysical Processor System sea ice motion products have allowed observing the strong localization of the deformation of sea ice along linear features, termed "linear kinematic features" [START_REF] Kwok | Deformation of the Arctic Ocean sea ice cover: November 1996 through[END_REF], which corresponds to active faults within the ice cover and the associated discontinuities in the drift velocities (see figure 1.2).

From the modelling point of view, the approach that appears the most natural to represent the dynamics of the ice cover at the scale of an aggregate of a few floes (see figure 1.1) is one in which individual ice plates are resolved and floe-to-floe interactions are treated explicitly. Discrete element-type methods have indeed been recently developed to model the mechanics and kinematics of sea ice at these small scales (e.g. [START_REF] Hopkins | A discrete element lagrangian sea ice model[END_REF][START_REF] Herman | Molecular-dynamics simulation of clustering processes in sea-ice floes[END_REF][START_REF] Wilchinsky | Modelling the reorientation of sea-ice faults as the wind changes direction[END_REF][START_REF] Rabatel | Dynamics of an assembly of rigid ice floes[END_REF]. In the context of long-term and global scale simulations of the sea ice cover, such models however remain computationally too expensive. As in general they do not handle fracturing processes within the floes, they are also better suited for the representation of a low concentration ice cover (< 80%) than of a dense ice pack. Current regional and global climate models as well as most operational modelling platforms are instead based on a continuum mechanics description of sea ice. Besides computational e ciency, this approach presents the advantage of being more suitable for a coupling with atmosphere and ocean models. In this case the motion of the thin ice cover is described by a with ⇢ the density, u the mean velocity, h the mean thickness of the ice and F ext , the external forces on the ice cover, typically the air and ocean drags, the Coriolis force and the sea surface tilt. The last term in this momentum balance equation stands for the mean internal force that arises from the sum of all mechanical interactions between ice floes. Without a treatment of the kinematics of the discrete floes or individual leads, the relationship between the corresponding internal stress tensor, , and the macroscopic deformation of the ice cover must be prescribed through a rheological law. Perhaps the biggest challenge of the continuum modelling approach lies in the formulation of this constitutive law [START_REF] Feltham | Sea ice rheology[END_REF], which must allow representing adequately through mean quantities a material with a inherently discontinuous character.

Space (m)

Current operational modelling platforms, whether assimilating data or not (e.g., TOPAZ4 : Sakov et al.

(2012), GIOPS : [START_REF] Smith | Sea ice forecast verification in the Canadian Global Ice Ocean Prediction System[END_REF]), and global climate models including sea ice dynamics (e.g., the Coupled Model Intercomparison Project Phase 5 models involved in the IPCC Fifth Assessment Report [START_REF] Flato | Climate change 2013: The physical science basis. contribution of working group i to the fifth assessment report of the intergovernmental panel on climate change[END_REF]) are based on the same rheological 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 modelling frameworks 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 displacements 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 continuum 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.

The Elasto-Brittle approach

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]Schulson, 2006a;[START_REF] Schulson | Creep and Fracture of Ice[END_REF][START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF][START_REF] Weiss | Coulombic faulting from the grain scale to the geophysical scale: lessons from ice[END_REF]. By investigating the dispersion of ice buoys over the Arctic, [START_REF] Rampal | Scaling properties of sea ice deformation from buoy dispersion analysis[END_REF] recently showed that sea ice 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]. The form of the scaling laws itself implies that the spatial heterogeneity and intermittency characterizing sea ice deformation decrease very slowly when increasing the scale of observation and therefore indicates that even at large scales, averaging the mechanical properties and deformation rates of the Arctic sea ice does not give rise to a smooth, possibly viscous-like behaviour as suggested by [START_REF] Hibler | A viscous sea ice law as a stochastic average of plasticity[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 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 spring-like 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 di↵use to localised damage through elastic interaction[END_REF], who combined

• a linear-elastic constitutive law for a continuum solid,

• a local Mohr-Coulomb criterion for brittle failure,

• an isotropic progressive damage mechanism for the elastic modulus described by a non-dimensional scalar damage parameter, allowing for the redistribution of the stress from over-critical to sub-critical areas of the material, for the triggering of avalanches of damaging events and the propagation of faults.

Their model was shown to reproduce macroscopic brittleness and scaling laws in the distribution of damage events, consistent with earthquake observations [START_REF] Amitrano | Brittle-ductile transition and associated seismicity: Experimental and numerical studies and relationship with the b value[END_REF].

In light of the analyses of sea ice deformation of [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] Rampal | Scaling properties of sea ice deformation from buoy dispersion analysis[END_REF] and others and of the similarity between the mechanical behaviour of sea ice and that of the Earth crust emphasized by [START_REF] Marsan | Space/time coupling in brittle deformation at geophysical scales[END_REF], this rheological framework, named Elasto-Brittle (EB) was recently developed in the context of the Arctic ice pack by Girard et al. (2010b) 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 pack, but using realistic wind forcing from reanalyses, showed that a Lagrangian EB model is able to reproduce the strong spatial localization and the anisotropy of damage within sea ice and to simulate deformation fields in good agreement with that reconstructed from the RADARSAT Geophysical Processor System (RGPS) ice motion data (Girard et al., 2010b).

1.0.3 The small and large deformations of sea ice

In the context of longer-term simulations of ice conditions and coupling to an ocean component, a suitable sea ice model needs to represent not only the small deformations associated with the fracturing of the ice pack, but also the permanent deformations occurring once the pack is fragmented and 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 point is an important and intrinsic limitation of the EB framework. To illustrate this limitation, let us represent schematically a linear-elastic and damageable material by a spring, with one free end and an initially undamaged elastic modulus E 0 , as shown in figure 1.3a. Let us further consider that the maximum stress this material (i.e., the spring) can sustain without being structurally damaged is c . If a stress, , smaller than c , is applied at time t, the resulting deformation of the spring, " is purely elastic and directly proportional to the stress, as shown on figure 1.3b (dotted line). Conversely, if > c , the spring becomes damaged, which translates into a diminution of its elastic modulus E and a non-linear deformation (figure 1.3b, solid black curve). and a progressive damage mechanism. At time t, an overcritical stress, > c, is applied to the system. It is removed at time t + t. In the all-elastic deformation assumption case (1), the damaged spring goes back to its initial position when the loading is removed. In the all-permanent deformation assumption case (2), it keeps its final position. In both cases, the damaged spring has a degraded elastic modulus, E < E 0 . (b) Stress-strain diagram for the linear elastic (dotted black line) and EB (solid and dashed lines) model. The black dot indicates the onset of damaging in the EB model, after which the model behaviour diverges from the linear-elastic case. The dashed red unloading path corresponds to the all-elastic deformation limit (1) and the dashed purple unloading path, to the all-permanent deformation limit (2).

1 t t + t E < E 0 2 E 0 E < E 0 " E 0 E < E 0
Because its structure has been damaged, if the spring is released at time t + t, it is not expected to go back to its exact initial position. Instead, part of its deformation is expected to remains permanent. Yet, if the deformation of the spring is modelled using a linear-elastic constitutive law, i.e., Hooke's law, as is the Elasto-Brittle framework [START_REF] Amitrano | From di↵use to localised damage through elastic interaction[END_REF]Girard et al., 2010b), the model will solve for its total deformation " : it will not distinguish between the elastic (recoverable) and the potentially permanent part of its deformation. Hence to estimate the deformation rate (i.e., the velocity) of a damaged elastic material in such a model, assumptions about the amount of reversible versus irreversible deformation must be made. The partitioning is bounded by two limit cases, which can be illustrated using the same linear spring system (see figure 1.3a).

1. If all of the deformation of the damaged material (spring) is assumed elastic, whatever its damage level, the material goes back to its initial position when unloaded (figure 1.3a, case 1) as represented by the red, dashed unloading path on figure 1.3b. The elastic modulus of the material is degraded, but its velocity in this case is zero. This assumption was made by Girard et al. (2010b), who neither updated the position of mesh element nodes nor estimated advection processes in their Lagrangian, short-term simulations of the Arctic sea ice cover.

2. If all of the deformation of the material is considered permanent, the material always keeps its final position when unloaded (see figure 1.3a, case 2 and figure 1.3b, purple dashed unloading path). In this case, velocities and deformation rates can be trivially estimated as the ratio of the total deformation and of the time associated with the loading process ( t).

In an Elasto-Brittle model for sea ice, the second assumption might be justified by the fact that elastic deformations within an undamaged ice pack are small compared to the permanent deformations associated with the opening, closing, and shearing along leads. Considering 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], 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 a time scale for these stress measurements of 1 day, upper bound values for daily shear strain rates in a one meter thick, purely elastic ice pack would be on the order of 10 5 day 1 . This is less than the lower bound deformation rate estimates from RGPS data (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 assumption is taken in the recently developed neXtSIM sea ice model, which is based on the EB rheology and does estimate advective processes over the Arctic [START_REF] Bouillon | Presentation of the dynamical core of neXtSIM a new sea ice model[END_REF][START_REF] Rampal | neXtSIM: a new Lagrangian sea ice model[END_REF]. However in the all-permanent deformations limit, internal stresses are instantaneously dissipated in the damaged material.

Hence the memory of the stresses associated with elastic deformations is erased whenever the applied loading is removed or reset. Without carrying this history of previous stresses, the model cannot reproduce the intermittency intrinsic to the mechanical behaviour of the material. In the context of sea ice modelling, this implies that the model only reproduces the part of the intermittency that is inherited from fluctuations of the wind forcing and cannot represent the observed long-range temporal correlations and extreme fluctuations in the deformation that emanate from elastic interactions within the ice cover [START_REF] Weiss | Intermittency of principal stress directions within Arctic sea ice[END_REF][START_REF] Rampal | Scaling properties of sea ice deformation from buoy dispersion analysis[END_REF]. In order to represent the observed properties of the deformation of sea ice and estimate adequate drift velocities, a suitable rheological model must therefore have the capability to distinguish between reversible and irreversible deformations.

The goal of this work is to develop such a model allowing a transition from the small/elastic to the 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 law takes the form of the Maxwell viscoelastic model. The all-important di↵erence 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]. Instead, it is an "apparent" viscosity which, like the elastic modulus in the Elasto-Brittle framework, is coupled to the level of damage of the ice cover and allowed to evolve in both space and time according to a progressive damage mechanism. The coupling is designed so that strains remain elastic 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 by the existence of similar approaches to model lithospheric faulting. Active faults in the Earth crust have been known to deform in two distinct manners: either abruptly, causing earthquakes, or in a 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 ismetpa?a segment of the north anatolian fault (turkey) from persistent scatterer insar[END_REF]. Similar to sea ice, co-seismic fracturing activates aseismic creep, leading to deformations that can be much larger than that associated with the fracturing itself [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 ismetpa?a segment of the north anatolian fault (turkey) from persistent scatterer insar[END_REF] and to slip rates that decrease progressively over years to decades due to various healing processes [START_REF] Gratier | Postseismic pressure solution creep: Evidence and time-dependent change from dynamic indenting experiments[END_REF]. Recent studies on aseismic faults (e.g., the Izmit fault) have suggested that creep relaxes a significant amount of elastic strain along the fault, retarding stress accumulation along some portions (and concentrating stresses on other locked portions), thereby demonstrating that this dissipative process needs to be included in earthquakes models [START_REF] Cakir | Onset of aseismic creep on major strike-slip faults[END_REF][START_REF] Gratier | Postseismic pressure solution creep: Evidence and time-dependent change from dynamic indenting experiments[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 figure 2.5), rocks along active faults are highly fragmented. Sheared granular media flow in a viscous manner when inertial e↵ects 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 closed-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 di↵erent stages of geological faulting, from subcritical crack growth to increasing crack concentration and material degradation, macroscopic brittle failure, post failure deformation and healing. Their progressive damage mechanism involved a scalar variable similar to that used in the EB and our Maxwell-EB framework, but 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-viscoplastic 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 (with the damage still evolving according to energy conservation) to account for the gradual accumulation of irreversible strain observed in typical rock mechanics experiments. Doing so they obtained better agreement with the observed deformation near the peak loading stress preceding the macroscopic failure of rock samples. The use of a viscous relaxation term in their rheological framework therefore had a fundamentally di↵erent purpose than in the present approach in that it was intended for the representation of 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 in the context of sea ice modelling.

A finite element sea ice model

In terms of formulating the discrete approximation of the equations of motion, finite element (FE) are known for the ease with which they can handle structured or unstructured meshes, complex geometries, boundary conditions, interface conditions and material properties. Although the relevance of these methods in the context of sea ice modelling were recognized as early as the 1970's Arctic Ice Dynamics Joint Experiment (AIDJEX) [START_REF] Becker | Finite element method in AIDJEX[END_REF], finite di↵erences (FD) have been common practice among the sea ice community due to the practicality of regular grids and because of their widespread use in ocean and climate models.

In late years, significant e↵orts have been put into incorporating FE methods in sea ice and ocean general circulation models in order to take advantage of multi-resolution mesh grids and better resolve complex flows with local phenomenons. In 2004, Wang and Ikeda presented a FE formulation for a dynamical sea ice model based on the Hibler VP rheology and proved their numerical scheme to be both stable and e cient. [START_REF] Danilov | A finite-element ocean model: principles and evaluation[END_REF][START_REF] Wang | The finite element sea ice-ocean model (FESOM) v.1.4: formulation of an ocean general circulation model[END_REF][START_REF] Danilov | Finite-element sea ice model (FESIM), version 2[END_REF]. The current version of this model can use both the standard VP rheology and its "numerically-elastic" version, the Elastic-Viscous-Plastic rheology of [START_REF] Hunke | An elastic-viscous-plastic model for sea ice dynamics[END_REF]. The FE approach has also been taken in the development of new rheological frameworks for sea ice, such as the elatic-decohesive [START_REF] Schreyer | Elastic-decohesive constitutive model for sea ice[END_REF] and the Elasto-Brittle model (Girard et al., 2010b;[START_REF] Bouillon | Presentation of the dynamical core of neXtSIM a new sea ice model[END_REF][START_REF] Rampal | neXtSIM: a new Lagrangian sea ice model[END_REF].

In building the Maxwell-EB model, we choose to take advantage of the recent expansion of finite element methods in the fields of sea ice and climate modelling. The definition of the advection scheme, in particular, is key to the numerical development of the model. The Maxwell-EB rheology aims, and succeeds, in reproducing the extreme gradients in the fields of sea ice deformation. A suitable advection scheme must therefore handle the transport of these gradients. The use of finite elements as opposed to finite di↵erences leaves the option of a Eulerian, Lagrangian or mixed approach. With the idea of taking the path that is the most compatible with the numerical habits of the climate modelling community and foreseeing an ineluctable coupling with an atmosphere and ocean component, we decide on an Eulerian approach. We implement discontinuous Galerkin methods for handling advection as well as other non-linear terms in the Maxwell-EB constitutive law. To our knowledge, such methods have not yet been employed in a sea ice model based on an elastic rheology. [START_REF] Lietaer | The e↵ects of resolving the Canadian Arctic Archipelago in a finite element sea ice model[END_REF] for instance, have used in their finite element sea ice model a finite volume upwind-weighted scheme with constant-by-element approximations for the ice thickness and concentration fields, which coincides with the discontinuous Galerkin method for FE approximation of order 0, together with linear nonconforming approximations for the ice velocity. Their mechanical framework however, was based on the viscous-plastic rheology, which is known to be defective in reproducing the extreme gradients within the ice cover [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)), and their numerical scheme did not include transport terms for the internal stress tensor. This makes the implementation and testing of discontinuous Galerkin methods within our Maxwell-EB framework an interesting experiment in itself.

Composition of the dissertation

This thesis relates the entire process of building the Maxwell-EB model for sea ice, from the laying of the theoretical grounds to the validation the physical approach, the building of the numerical scheme and the implementation of the mechanical framework in the context of sea ice modelling. The report is structured as follow:

• Chapter 2 presents the Maxwell-EB framework, and in particular the concept of stress relaxation. The progressive damage and healing mechanisms are described, as well as the coupling between these mechanisms and the mechanical properties of the material, which constitute a fundamental feature of the rheological model. Follows a discussion of the important characteristic times and numbers involved in the Maxwell-EB framework and of their relative values in the context of sea ice modelling. First idealized simulations are analyzed that use a simplified, small-deformation version of the model in which advection processes are not yet included. These simulations provide a validation of the mechanical framework in terms of its capacity in reproducing a highly heterogeneous, anisotropic and intermittent deformation.

Finally, a sensitivity analysis to the value of one key, yet poorly constrained, model parameter is presented that allows exploring the range of mechanical behaviours reproduced by the model.

• Chapter 3 describes the transition between the small and the large-deformation Maxwell-EB model.

In terms of the numerics, this transition is twofolds. It consists in (1) introducing advective processes and (2) treating the full, objective internal stress tensor derivative arising in the Maxwell-EB constitutive equation. In the first part of this chapter, the discontinuous Galerkin approximation of the advection term for transported quantities and of the constitutive equation are introduced. Simple tests are performed to explore the limitations of the numerical scheme. In the second part, the numerical scheme for the full, large-deformation Maxwell-EB dynamical model is presented. Large-deformation simulations are carried that use a typical Couette flow geometry. These are compared to the results of a laboratory Couette experiment performed on a thin plate of fresh-water ice and are discussed in terms of the simulated mechanical behaviour.

• In chapter 4, the Maxwell-EB framework is implemented in the context of modelling sea ice on regional to global scales. In particular, evolution equations for the ice thickness and concentration are introduced in the model and their simple coupling with the Maxwell-EB mechanical parameters is described. Simulations of the flow of ice through a narrow, idealized channel with dimensions consistent with straits over the Arctic, for instance in the Canadian Arctic Archipelago, are presented. These are analyzed in terms of the simulated velocity fields, internal stress states and thickness distributions and with a particular emphasis on the capacity of the mechanical framework to represent the formation of arch-like leads and stable ice bridges.

• advective processes are neglected is described in section 2.3. The numerical scheme for these small-deformation model experiments is presented in section 2.3.1. These simulations are analyzed in section 2.4 on the basis of the macroscopic behaviour and convergence properties of the model, and in particular, of the heterogeneity, anisotropy and intermittency of the simulated deformation. Section 2.5 presents an analysis of the sensitivity of the mechanical model on the value of its least constrained parameter. Conclusions for this chapter are summarized in section 2.6.

The building of a viscous-elastic-brittle continuum model

Viscoelastic materials, whether fluid or solid, exhibit both viscous and elastic properties at the microscopic scale. Elasticity is inhered from the reversible stretching of bonds along crystallographic planes or from the flexibility of their molecules. In continuum models, this property is characterized at the macroscopic scale by an elastic modulus, E. The viscous properties of viscoelastic materials are attributable to rearrangements of their molecular structure, to the motion of defects (vacancies, dislocations), to the frictional sliding between grains, or to the friction between polymer and solvent molecules. These irreversible plastic deformations at the microscopic scale result in energy dissipation when viscoelastic materials are subjected to an external stress, as opposed to purely elastic materials. At the macroscopic scale, these viscous e↵ects translate into a dependence of the deformation of the material on time, characterized in continuum models by a constant or stress-dependent dynamic viscosity ⌘.

In 1867, James Clerk Maxwell presented a linear model suitable for the macroscopic behaviour of a viscoelastic material, typically an incompressible fluid, undergoing small deformations [START_REF] Maxwell | On the dynamical theory of gases[END_REF]. This model is represented schematically in one dimension by a spring and dashpot connected in series (see figure 2.1b). The elastic component of this system obeys the linear elastic constitutive relation

= 2E" (2.1)
where is the elastic (shear) stress, E is the elastic (shear) modulus and " = D(U) = 1 2 (rU + rU T ) is the strain tensor, defined in terms of the displacement U. The viscous component of the system is described by a linear Newtonian fluid constitutive law = 2⌘ "(u)

(2.2) with ⌘, the material's (stress-independent) viscosity and "(u) = D(u) = 1 2 (ru + ru T ), the rate of deformation tensor, defined in terms of the velocity, u.

When, at a time t, a stress is suddenly applied to this Maxwell system, the resulting deformation is split between the instantaneous, reversible, deformation of the spring, " E , and the permanent deformation of the dashpot, " , increasing linearly with time, such that the total deformation of the system, ", is given by the sum " = " + " E .

(2.3) (see figure 2.1b). If after a time t + t the stress is removed, the spring goes back to its initial position, but the dashpot retains its deformation. Conversely, when a deformation is applied and the system is held to a given strain rate, the resulting internal stress dissipates exponentially with time.

The Maxwell constitutive law is deduced by summing the viscous and elastic components of the deformation given by equations 2.1 and 2.2 :

1 E D Dt + 1 ⌘ = 2 "(u)
This relationship can be written equivalently in the from

D Dt + E ⌘ = 2E "(u)
where the ratio ⌘ E is interpreted as a relaxation time, setting the rate of dissipation of the applied stress through the permanent deformation of the dashpot and hence characterizing the capacity of the viscoelastic material to retain the memory of reversible deformations. As both the viscous and elastic components of the deformation are solved for simultaneously in the Maxwell model, deformation is defined naturally in terms of the strain rate.

(a) a uniform material with elastic (shear) modulus E and viscosity ⌘ (top panels) and associated stress-strain diagrams (bottom panels). At time t, a stress is applied on each material. It is removed at time t + t. In the linear elastic case, the material goes back to its initial position when the stress is removed and no energy is dissipated in the loading-unloading process. In the Maxwell viscoelastic case, the total deformation is split between the deformation of the elastic (spring) and of the viscous (dashpot) components. When the stress is removed, the elastic part of the deformation is recovered, but viscous deformations remain permanent. In the Maxwell model, the energy dissipated through permanent deformations is represented by the area between the loading (blue) and unloading path (red) on the stress-strain diagram.

" E " (b) " " " = " E + " ⌘ " E t t + t " = " E E = 2E" E = 2⌘ "

Towards a Maxwell-Elasto-Brittle rheology

In developing the Maxwell-EB rheology, we extend the underlying concepts of the Maxwell viscoelastic model to an isotropic, compressible solid, in particular, the idea of stress relaxation and the partionning of the deformation into an elastic, recoverable, and a permanent component. As in the Elasto-Brittle framework, the elastic deformation of the ice cover obeys Hooke's law for an isotropic material. In tensor form, this linear elastic constitutive law reads

= EK : " E (2.4)
with the Cauchy stress tensor, E, the elastic (Young) modulus, K a sti↵ness tensor defined in terms of ⌫, the Poisson's ratio, such that for all symmetric tensor ✏ = ✏ ij 8 i, j;

1  i, j  3, (K : ✏) ij = ⌫ (1+⌫)(1 2⌫) tr(✏) ij + 2 1 2(1+⌫) ✏ ij .
As in the Maxwell model, we consider the total deformation to be the linear sum of an instantaneous, elastic and of a time-dependent, permanent contribution, such that " = " + " E . Substituting for " E in terms of the total deformation in equation 2.4, the constitutive law writes, = E(K : " K : " ) or, in terms of the rate of deformation tensor

D Dt = E (K : "(u) K : "(u) ) .
Following the idea of stress relaxation behind the Maxwell model, we consider the time-dependent part of the deformation to be a linear function of the applied stress of the form

= ⌘K : "(u) (2.5)
where ⌘ sets the rate of increase of the permanent deformation with time and has the dimensions of a viscosity.

Substituting for "(u) , the constitutive equation therefore becomes

D Dt = E ✓ K : "(u) 1 ⌘ ◆ , or D Dt + E ⌘ = EK : "(u),
or again, in "Maxwell-like" form:

1 E D Dt + 1 ⌘ = K : "(u). (2.6)
As in the Maxwell model, the ratio ⌘ E has the dimension of a time, and sets the rate of dissipation of the internal stress through the permanent deformation of the material. In the following, we refer to this ratio as the relaxation time for the stress, . The equation 2.5 for the viscous deformation of the compressible material implies that this rate is assumed equal for both the volumetric and deviatoric components of its deformation.

The mechanical parameter ⌘ therefore di↵ers intrinsically from the dynamic viscosity of a fluid. This point, which is further discussed in section 2.1.2, is easier to comprehend when expending K : "(u) and introducing the (adimensional, i.e., normalized by

E) Lame parameters ⇤ = ⌫ (1+⌫)(1 2⌫) and G = 1 2(1+⌫) , D Dt + E ⌘ = E [⇤tr ( "(u)) I + 2G "(u)] .
(2.7)

If the material was incompressible, div(u) = tr ( "(u)) = 0, the constitutive law would then coincide with that of the original Maxwell model.

Bridging small and large deformations

Besides compressibility, an important distinction between the standard Maxwell and the Maxwell-EB model pertains to the scale the represented deformations. On the one hand, the linear Maxwell model, based on a constant elastic modulus and viscosity, is usually suitable for the representation of the small deformations of viscoelastic materials [START_REF] Maxwell | On the dynamical theory of gases[END_REF]. The Maxwell-EB framework on the other hand was conceived with the particular intention of allowing the representation of both small and large deformations of a continuum solid. The transition from small to large deformations requires introducing some non-linearity in the rheological model. In the Maxwell-EB model, non-linearity is twofold.

First, non-linear terms are introduced in the Maxwell-EB constitutive law through the material derivative of the internal stress tensor. The material derivative D Dt of any scalar quantity ' writes

D' Dt = @' @t (u • r)'
and can be approximated as @' @t in the limit of small deformations, hence neglecting the non-linear advection term. In the case of large deformations, advective processes cannot be neglected and this form of the material derivative is not objective in the case of vectorial or tensorial quantities. In the Maxwell-EB constitutive law, the material derivative of the internal stress tensor therefore takes the form of the objective Gordon-Schowalter derivative [START_REF] Saramito | Complex fluids. Modelling and algorithms[END_REF])

D Dt = @ @t + (u • r) + a (ru, )
where the additional term a accounts for the e↵ects of rotation and deformation on the evolution of the stress tensor and expresses as

a (ru, ) = W (u) W (u) a ( D(u) + D(u) ) with D(u) = ru+ru T 2
the symmetric and W (u) = ru ru T 2 the anti-symmetric part of the velocity gradient.

The cases of a = 0, 1 and 1, represent the Jaumann, upper convected and lower convected objective derivative respectively.

With E and ⌘ constant in the standard Maxwell constitutive equation (2.6), the model is linear with respect to (but induces a coupling between and u as u is also an unknown in the conservation of momentum equation). When a load, i.e., a stress, is applied to the material, permanent deformations increase linearly with time. Therefore if the loading is applied for a long enough time, irreversible deformations can accumulate and become potentially large. Yet, if the viscosity, ⌘, was taken on the order of the dynamic viscosity for undamaged sea ice [START_REF] Duval | Rate-controlling processes in the creep of polycrystalline ice[END_REF], the simulated strains would remain small on time scales relevant for the drift and deformation of the ice cover. In order for the total deformation to become large compared to elastic deformations on the time scales considered in regional and global sea ice models, the mechanical parameters E and ⌘ in the Maxwell-EB framework are allowed to evolve both spatially and temporally according to its local characteristics such as its level of damage, which in the original EB model [START_REF] Amitrano | From di↵use to localised damage through elastic interaction[END_REF], represents the weakening of the material when stressed beyond its mechanical resistance (see section 2.2.3).

Another property impacting the mechanical strength of the material is its compactness, referred to as the packing fraction in the case of granular materials and concentration in the specific case of sea ice (ratio of the model grid cell covered by ice as opposed to open water). This aspect is further discussed in chapter 4. This dependence of the mechanical parameter ⌘ and E on the physical characteristics of the material e↵ectively introduces additional non-linearity in the model, as the constitutive equation becomes coupled to evolution equations for these characteristics.

The stress-strain diagram of figure 2.2 contrasts the behaviour the linear-elastic (dotted black curve), original EB (dashed red, blue and purple curves), original Maxwell (solid red and blue curves) and Maxwell-EB models (solid black curves) in terms of their capacity in partionning the elastic and permanent deformations and the relative size of both types of deformation.

"

" E E < E 0 E 0 " Figure 2
.2: Loading-unloading paths for a material with initial elastic modulus E 0 in the linear-elastic (dotted black line), Elasto-Brittle (dashed colored), linear Maxwell (solid colored) and non-linear Maxwell-EB (solid black curves) model. Unlike the purely elastic and EB models, the viscoelastic models allow partitioning the total deformation into a permanent and an elastic contribution. The non-linear Maxwell-EB model allows for permanent deformations to become much greater than in the standard, linear Maxwell model. It is important to note that the diagram is still 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.

The Maxwell-EB model

The previous section has introduced the Maxwell-EB constitutive law and, in particular, the dissipative mechanism for the internal stress, which allows part of the deformations of a Maxwell-EB material to become large and permanent. The present section presents the remaining ingredients of the Maxwell-EB framework, which are:

1. the criterion for brittle failure, 2. the disorder introduced in the failure criterion, representing the material's heterogeneity, 3. the damage mechanism, 4. the healing mechanism, 5. the coupling between the Maxwell-EB mechanical parameters and the damage and healing mechanism.

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][START_REF] Weiss | Coulombic faulting from the grain scale to the geophysical scale: lessons from ice[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 (2.8) (or 2 = q 1 + c
, by symmetry of the criterion along the the 1 = 2 axis -see figure 2.3). In the following, symmetry is always implied and only the lower branch of the damage criterion is discussed). The slope of the envelope in the principal stresses plane, q, is expressed in terms of the internal friction coe cient µ as

q = h (µ 2 + 1) 1/2 + µ i 2 .
(2.9)

In the model, µ 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]. The intercept c of the MC criterion with the 1 axes (see figure 2.3), interpreted as the uniaxial (unconfined) compressive strength, is given by

c = 2C ⇥ (µ 2 + 1) 1/2 µ ⇤ .
(2.10)

where C represents the cohesion of the material and sets its resistance to pure shear.

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 di↵erent 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]. In the Maxwell-EB model, the Mohr-Coulomb criterion is therefore extended to tensile stresses. For practical reasons, the critical value is set 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] (see figure 2.3). The tensile strength cuto↵ therefore takes the form: (2.11) where

1 < 0; 2 = t ,
t = c q = 2C h (µ 2 + 1) 1/2 + µ i .
(2.12)

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 of 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 a↵ect our results.

No truncation to the MC criterion is used to close the envelope towards biaxial compression (i.e., beyond c ) as instances of such 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 cuto↵ is represented in figure 2.3 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.

Disorder

In the Maxwell-EB model, a noise is introduced in the damage criterion via the spatial distribution of the cohesion of the material, C (see equations 2.10 and 2.12). This noise represents the natural heterogeneity of a real material that is associated with various structural defects at the sub-grid scale like vacancies, deviations in the crystalline order or thermal cracks in the case of sea ice [START_REF] Schulson | Compressive shear faults within arctic sea ice: Fractures on scales large and small[END_REF][START_REF] Schulson | Creep and Fracture of Ice[END_REF], and which serve as stress concentrators. This heterogeneity causes the progressive failure of the material even under homogeneous forcing conditions (e.g., [START_REF] Herrmann | Statistical Models for the fracture of disordered meria[END_REF][START_REF] Amitrano | From di↵use to localised damage through elastic interaction[END_REF][START_REF] Tang | Numerical simulation of progressive rock failure and associated seismicity[END_REF]. In the case of the ice pack, the correlation length associated to the natural heterogeneities is likely much smaller than the typical spatial resolution of models ( 1 km). Hence heterogeneity in the Maxwell-EB model is introduced at the smallest resolved scale, that is, the mesh element size x. The value of C over each model element is therefore 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]. Although not physically strictly tied to the cohesion of the simulated material, heterogeneity is a necessary ingredient in the model, at least when using idealized domain geometries and homogeneous forcing and initial conditions [START_REF] Bouillon | Presentation of the dynamical core of neXtSIM a new sea ice model[END_REF][START_REF] Rampal | neXtSIM: a new Lagrangian sea ice model[END_REF]. In other damage mechanics model it is alternatively introduced through other mechanical resistance parameters, for instance in the spatial distribution of the elastic modulus (e.g., [START_REF] Cowie | Statistical physics model for the spatiotemporal evolution of faults[END_REF]. If the evolution of the sub-grid scale heterogeneity is considered much slower than the propagation of fractures within the material, the disorder is time-independent, or "quenched", and the field of C is set once at the beginning of a model simulation. Conversely, if disorder is considered to evolve with the fracturing and damaging of the material (e.g., [START_REF] Amitrano | From di↵use to localised damage through elastic interaction[END_REF], it is said "annealed" and the local value of C is then time-dependant.

4 3 2 1 1 2 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 an Eulerian framework, a transport equation of the form

@C @t + (u • r)C = 0 (2.13)
is necessary to handle the advection of the field of C with the simulated velocity field.

Progressive damage mechanism and healing

As introduced in section 2.1.2, the Maxwell-EB rheology di↵ers fundamentally from the standard Maxwell rheology in that the mechanical parameters E, ⌘ (and, as discussed below, ) are not constant but coupled to the spatially and temporally evolving level of damage of the simulated material, which controls its local degradation and re-increase in strength.

Consistent with previous damage rheological models, this property is represented in the Maxwell-EB framework by a dimensionless scalar parameter, d (e.g., [START_REF] Amitrano | From di↵use 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), and is interpreted as a measure of sub-grid cell defects or crack density [START_REF] Kemeny | E↵ective moduli, non-linear deformation and strength of a cracked elastic solid[END_REF]. The value of d evolves between 1 (undamaged) and 0 ("completely damaged" material).

The level of damage 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 these mechanisms also contribute to the increase in elastic sti↵ness (E ⇥ h, with h, the ice thickness) and e↵ective 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] 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, while respecting the following constraints:

1. d decreases with damaging, which occurs only when the stress state becomes over-critical, 2. d increases with healing, 3. the increase in d due to healing cannot o↵set its decrease due to damaging, i.e., d cannot increase beyond its undamaged value of 1.

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 (and EB) model. Instead, when the stress locally exceeds the critical stress, the level of damage d and the elastic modulus, are allowed to drop, leading to local strain softening [START_REF] Amitrano | From di↵use 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]. Because of the long-range interactions in space using a (here, structured) mesh of triangular elements and the damage criterion is set through the field of cohesion C, in which noise is introduced at the scale of the element. Forcing, homogeneous or not, is applied on the simulated material as an external force or prescribed displacement and the internal stress within the simulated material builds up such that at time t it exceeds the damage criterion over a given element (circled in red). At the next time step, t + t, the element is by this fact weakened and its post-damaging level of damage is d 0 < d. Because it is weaker, the part of the internal stress it can no longer withstand is redistributed to its close neighbours through elastic interactions. By this redistribution, some of the neighbouring elements (circled in red) become in turn over-critical. At time t + 2 t, the weakening of these elements triggers further damaging events.

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 in the following by 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. This assumption holds when considering not a single model element, but a matrix of elements. In this case, when damaging occurs locally and neighbouring elements remain undamaged, the first e↵ect is a stress redistribution between these neighbouring elements. The second assumption is that, for a su ciently small model time step t, i.e. very small compared to the viscous relaxation time (see section 2.2.4), a negligible part of the stress is dissipated into viscous deformations. 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 damaging acts so that after being damaged, an element has its state of stress lying just on the failure envelope. With these assumptions made, the following equality holds for each damaged element:

" 0 = " ! K 1 0 E ⇥ d crit = K 1 E , (2.14)
where the superscript 0 denotes the post-damage state of deformation and stress (lying on the failure envelope).

Using Hooke's law (given by 2.4 and here under plane stress conditions), the equality (2.14) can can be written for each component of the stress tensor and principal stress, such that the change in level of damage of a given element can be written as

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

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 figure 2.3).

Two cases must be distinguished when calculating 0 , hence d, depending on which of the Mohr-Coulomb or tensile criteria has been exceeded.

1. In the first case, 1 > 0 and 1 q 2 > 0 and the post-damage state of stress is the intersection of the two lines of equations

0 1 = q 0 2 + c , 0 1 = 1 2 0 2 .
Hence the increment in damage is

d crit = 8 < : min(1, c 1 q 2 ) if 1 q 2 > 0, 1 otherwise.
(2.16) .

2. In the second case, 1 < 0 and 2 < 0 and the post-damage state of stress is given by the intersection of

0 1 = 1 2 0 2 . 0 2 = t ,
The increment of damage is

d crit = 8 < : min(1, t 2 ) if 1 2 , 1 otherwise.
(2.17)

Combining (2.16) and (2.17), d crit is evaluated simultaneously over all mesh elements of the model domain as:

d crit = min  1, t 2 , c 1 q 2 .
(2.18)

According to equations (2.10) and (2.12) for c and t , the local change in d is hence defined entirely in terms of the local state of stress and the local value of C.

Following previous progressive damage models (e.g., [START_REF] Amitrano | From di↵use to localised damage through elastic interaction[END_REF]Girard et al., 2010a), the level of damage, d, 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 (2.19)
with d n+1 crit and d n+1 , the distance to the damage criterion and level of damage at the (n + 1) th model time step and d n the level of damage at the previous time step. 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, c, of propagation of elastic waves 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 time-discretized equation for d is obtained

d n+1 d n t = d n+1 crit 1 d n 1 t d .
Taking the limit of t ! 0 the intrinsically discrete propagation mechanism can therefore be processed through the continuous equation

Dd Dt = d crit 1 t d d.
(2.20)

It is important to note that the direct correspondence between the discrete recursive sequence (2.19) and the continuous equation (2.20) relies on the specific use of an explicit scheme of order 1. This point is central in the development of the Maxwell-EB numerical scheme. It will be discussed more extensively in chapters 3 and 5.

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 sti↵ness (d = 1), which in a dynamic-thermodynamic sea ice model would depend on the local di↵erence between the temperature of the air near the surface of the ice and the freezing point of seawater below. In the uncoupled, dynamic model described here, t h is set constant in both space and time.

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. In this case the evolution equation for d reads

Dd Dt = 1 t h 1 d , 0  d  1.
However as there is no physical evidence that the healing rate should depend on the level of damage, in the following description and implementation of the Maxwell-EB model we use an even simpler parameterization that assumes a constant healing rate set to 1

t h : Dd Dt = 1 t h , 0  d  1. (2.21)
Combining both the damaging and healing mechanisms (Eq. 2.18, 2.20 and 2.21), the complete evolution

equation for d is Dd Dt = ✓ min  1, t 2 , c 1 q 2 1 ◆ 1 t d d + 1 t h , 0  d  1.
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, e↵ective only where and when the state of stress becomes overcritical. The characteristic time for this process, t d , is tied to the speed of propagation of (shear) elastic waves as well as to the scale of heterogeneities within the modelled material. 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], which is significantly smaller than the theoretical shear wave speed for bulk saline ice [START_REF] Gammon | Elastic constants of artificial and natural ice samples by brillouin spectroscopy[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 di↵erence between t h and t d therefore implies that the two processes are intrinsically decoupled in the case of the ice pack.

Coupling between d, E and ⌘

The coupling between the Maxwell-EB constitutive 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 and a very large apparent viscosity ⌘ 0 . In this case, the viscous term in (2.6) is negligible and a linear-elastic constitutive law is recovered (figure 2.5, 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, d ! 0, elastic interactions are hindered and deformations are strictly irreversible (figure 2.5, left panel). In this case, ! t d and a soft elasticplastic behaviour is recovered in which the memory of the elastic stresses is totally lost.

• As damaged areas heal, E, ⌘ and all re-increase, up to their initial undamaged values.

Di↵erent 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 di↵use 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 (2.22)
Undamaged! pack ice small, elastic deformations ! 0

Fractured,! damaged zones large, permanent deformations Here ↵ is a constant greater than 1, introduced to fulfil the constraint that the relaxation time for the stress decreases with increasing damage and increases with healing, as the material respectively looses and recovers the memory of reversible deformations. In the following, we refer to this constant as the damage parameter.

⌘ ! ⌘ min 1 E ! 0 1 ⌘ ! 0 E = E 0
Using this formulation, both ⌘ and E are entirely defined by their initial value, a constant, and by the level of damage. Hence the Maxwell-EB constitutive law (2.6) can be written as

1 E 0 d D Dt + 1 ⌘ 0 d ↵ = K : ".
However in this form, 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 cuto↵, ⌘ min ⌧ ⌘ 0 , on the value of the apparent viscosity can be introduced and the expression for ⌘(d) modified as

⌘ = (⌘ 0 ⌘ min )d ↵ + ⌘ min = 8 < : ⌘ 0 for d = 1, ⌘ min for d = 0.
(2.23) Substituting for ⌘ in the expression for the relaxation time (2.22), the elastic modulus then becomes

E = ⌘ 0 ⌘ min ⌘ 0 E 0 d + ⌘ min ⌘ 0 1 d ↵ 1 E 0 .
(2.24)

In this case, E ⇡ E 0 for d = 1, E decreases with d until a minimum at d(E min ) = h ⌘min ⌘ 0 ⌘min (↵ 1) i 1 ↵ and E ! 1 for d ! 0.
The functions E(d) and ⌘(d) above are represented in figure 2.6 for di↵erent values of ↵. Using such a cuto↵ 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 is however has no physical grounds in the context of a progressive damage model for a continuum solid and is not used here to represent the flow of a Newtonian fluid. It is rather introduced to insure the mathematical consistency and existence of a solution while retaining a continuous function for the level of damage. In the following implementation of the Maxwell-EB rheology, we take this approach instead of imposing a minimum value for d. This choice has no significant e↵ect on our

Characteristic numbers and times

In this section, the equations constituting the dynamical Maxwell-EB model are formulated in an adimensional form. This allows describing the framework in terms of important characteristic numbers and times and to explore its sensitivity to a reduced set of parameters. In order for the model to represent the intended physics, these parameters must evolve within a certain range of values, which is discussed at the end of the section.

Here and in the following, the Maxwell-EB rheology is implemented in the context of modelling the drift and deformation of a continuum material with horizontal extent very large compared to its vertical extent. The balance of forces is described by a momentum equation of the form 1.1. The model is made adimensional with respect to the horizontal extent L, the thickness H, the flow velocity U and internal stress ⌃ within the material. The time characterizing the deformation process is then T = L U . The superscript ' ˜' is used for all dimension-less variables and operators, which are listed in table 2.1. In terms of these adimensional variables and operators, the Maxwell-EB system of equations reads:

(the momentum equation)

⇢ hH U 2 L Dũ D t = ⌃ H L Fext + ⌃H L r • ( h˜ ) ⇢ U 2 ⌃ h Dũ D t = Fext + r • ( h˜ ) (the constitutive equation) ⌘ 0 E 0 d ↵ 1 U L ⌃ D˜ D t + ⌃˜ = U L ⌘ 0 ✓ (1 ⌘ min ⌘ 0 )d ↵ + ⌘ min ⌘ 0 ◆ K(⌫) : "(ũ) ⌘ 0 E 0 U L d ↵ 1 D˜ D t + ˜ = ⌘ 0 ⌃ U L ✓ (1 ⌘ min ⌘ 0 )d ↵ + ⌘ min ⌘ 0 ◆ K(⌫) : "(ũ)
In the following and throughout the text, we use the notation

d 0 ↵ = h⇣ 1 ⌘min ⌘ 0 ⌘ d ↵ + ⌘min ⌘ 0 i .
(the damage equation)

1 T Dd D t = ✓ min  1, t ⌃ 1 ˜ 2 , c ⌃ 1 ˜ 1 q ˜ 2 1 ◆ 1 t d d + 1 t h , 0 < d  1 Dd D t = ✓ min  1, t ⌃ 1 ˜ 2 , c ⌃ 1 ˜ 1 q ˜ 2 1 ◆ 1 t d /T d + 1 t h /T , 0 < d  1
In the limit of small deformations, all velocity gradients are small such that advection, rotation and deformation terms can all be neglected. If neglecting the e↵ect of the (small) elastic deformations on the material's density as well, mass conservation does not need to be imposed and the thickness of the material h remains constant. In the case of large deformations, additional equations are introduced that handle mass conservation (see chapter 4) and, in the case of quenched disorder, the transport of the field of cohesion C (see equation 2.13).

Variables, dimensions and operators Non-dimensional equivalent

Horizontal dimension x x = x L Time t t = t T Velocity u ũ = u U Internal stress ˜ = ⌃ External forcing F ext Fext = Fext ⌃ L H

Level of damage d d

Thickness

h h = h H Del Operator r r = Lr
Table 2.1: Dimensional model variables and operators and their adimensional counterpart.

In viscoelastic models, two di↵erent scalings can be considered for the characteristic stress : either ⌃ scales as the elastic modulus of the material or as the product of the material's viscosity and of the characteristic time for the deformation process, ⌘ ⇥ T . As the Maxwell-EB framework is meant to represent the deformation of a continuum solid undergoing progressive damage rather than the flow of a fluid material, scaling ⌃ as E 0 , the (undamaged, constant) elastic modulus, seems the most appropriate choice here. In this case, the non-dimensional set of equations reads:

Ca 0 h  @ ũ @ t + (ũ • r)ũ = Fext + r • ⇣ h˜ ⌘ (2.25) We 0 d ↵ 1  @ ˜ @ t + (ũ • r)˜ + a ( rũ, ˜ ) = We 0 d 0 ↵ K(⌫) : "(ũ) (2.26) @d @ t + (ũ • r)d = ✓ min  1, ⌃ t 1 ˜ 2 , ⌃ c 1 ˜ 1 q ˜ 2 1 ◆ 1 T d d + 1 T h , 0 < d  1. (2.27)
In this form, the model involves 8 characteristic numbers and time scales: Ca 0 , the (undamaged) Cauchy number, We 0 , the (undamaged) Weissenberg number, ⌫, Poisson's ratio, ⌃ t = t E 0 , the dimensionless critical tensile stress, ⌃ c = c E 0 , 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 parameter. 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.

The characteristic time of damaging, T d

As mentioned in section 2.2.3, the (adimensional) characteristic time for the propagation of damage, T d = t d T , is determined by the speed of propagation of elastic waves, c, within the simulated material and therefore is strongly tied to the mean spatial resolution of the model, as t d should be on the order of x c . In turn, T d 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 greater or equal to 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 and (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 characteristic time of healing, T h

In order for healing not to o↵set damaging in the rate of change of d, the (adimensional) time for healing, T h = t h T , must be set 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 e↵ectively 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 2.2.3), pack ice naturally meets this condition.

The Weissenberg number, We

The Weissenberg number, We, defined as the dimensionless product of the viscous relaxation time for the stress and of the deformation rate or, equivalently, the ratio of the viscous relaxation time and of the time characterizing the simulated deformation process

We = ⌘ E U L = T , (2.28)
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 locally according to the level of damage as We = We 0 d ↵ 1 with We 0 , its maximum, undamaged value.

On the one hand, as viscous dissipation should be insignificant over undamaged and strictly elastic portions 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 Maxwell constitutive law (2.6) e↵ectively 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 values of We approaching 1, numerical instabilities arise in Maxwell-type models due the presence of deformation source terms ( a ) in the transport equation for the stress tensor. 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 law can be multiplied by a Heaviside-type function d ⇤ that e↵ectively sets 1 ⌘ to 0 when and where d d c , with d c , a chosen threshold value, and leaves the constitutive equation unchanged otherwise. In this case the modified constitutive equation reads:

We 0 d ↵ 1 D˜ D t + d ⇤ ˜ = We 0 d 0 ↵ K(⌫) : "(ũ) with 8 < : d ⇤ = 0 if d d c d ⇤ = 1 if 0  d < d c (2.29)
When using a constant heal rate parameterization, d c can be simply set to 1. For small values of We 0 (⌧ 1) and/or large model time steps, a more regular version of d ⇤ could be a more suitable alternative in order to avoid the generation of numerical instabilities. In small-deformation simulations, i.e., in which the advection, rotation and deformation of the internal stress tensor are not accounted for and hence the value of We is not limited, or in simulations 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.

On the other hand, 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. In the case of a constant heal rate parameterization, this implies We

0 d ↵ 1 ⌧ T h (or 0 d ↵ 1 ⌧ t h ) where d is small.
Conversely, for the Weissenberg number to remain within a range of values ensuring the convergence of the numerical scheme, i.e., We < 1 for most schemes, deformation rates within the material should be smaller than 1 0 where the material is undamaged, i.e., strictly elastic. This point is discussed further in chapter 3, section 3.3.

Poisson's ratio, ⌫

Poisson's ratio is a measure of the compressibility of a homogeneous material. "Perfectly compressible" materials, i.e., showing no lateral expansion when compressed longitudinally, have ⌫ ! 0. Incompressible materials, Newtonian viscous fluids for instance, have ⌫ = 0.5. In this case, the first (adimensional) Lame parameter is ⇤ ! 1 and hence, r • u ! 0 (see equation 2.7). Sea ice is somewhat compressible, with a dynamic Poisson ratio of 0 < ⌫ < 0.5 depending on its temperature [START_REF] Timco | A review of the engineering properties of sea ice[END_REF]. As mentioned in section 2.1.1, this is an important di↵erence between the Maxwell-EB and the standard Maxwell viscoelastic model. In the limit of ⌫ = 0.5 the Maxwell-EB constitutive equation (2.7) takes the form of the standard Maxwell constitutive law.

In the current implementation of the Maxwell-EB rheology, ⌫ = 0.3 and is constant in both space and time, whatever the level of heterogeneity of the simulated material. Yet, rock mechanics experiments have shown that ⌫ e↵ectively increases with damaging [START_REF] Heap | The evolution of elastic moduli with increasing crack damage during cyclic stressing of a basalt from mt. etna volcano[END_REF], consistent with a dilatation of the fractured material [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF][START_REF] Martin | The progressive fracture of Lac du Bonnet granite[END_REF]. In other progressive damage models of elastic materials, Poisson's ratio is therefore also allowed to evolve as a function of d (e.g. [START_REF] Hamiel | A viscoelastic damage model with applications to stable and unstable fracturing[END_REF]. However, tests performed with the EB model of [START_REF] Amitrano | From di↵use to localised damage through elastic interaction[END_REF] suggest that accounting for this additional level of complexity has no significant impact on the simulated mechanical behaviour. The e↵ect of coupling ⌫ and d could be further investigated in future refinements of the Maxwell-EB rheology.

The Cauchy number, Ca

The dimensionless number that arises in place of the Reynolds number when adimensionalizing stresses in the momentum equation with respect to the elastic modulus as opposed to ⌘ ⇥ T , is the Cauchy number

Ca = ⇢U 2 E (2.30)
which interprets as the ratio of inertial to elastic forces.

If inertial forces are comparable to elastic forces and Ca ⇠ 1, the e↵ect 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. In order for the wave contribution not to have a significant e↵ect on simulated deformation and stress fields, Ca must 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 e↵ects can be safely neglected. For simulated ice velocities U < 1 ms 1 and ↵ > 2, inertial e↵ects in the Maxwell-EB model remain small when damage becomes important. In the case of faster flows the value of ↵ should be set > 2 in order for inertial e↵ects to be negligible. However, in the context of the ice pack, drift velocities larger than 1 ms 1 on space scales relevant for continuum sea ice models are seldom, perhaps never, encountered.

The damage parameter, ↵

The damage parameter ↵ is introduced in the Maxwell-EB framework in a completely ad-hoc manner, to control the rate at which the apparent viscosity drops and the material looses its elastic properties with damaging.

As mentioned in previous sections, one physical constraint on the value of ↵ is that the relaxation time should decrease with damaging : hence ↵ should be set greater than 1. The requirements that (1) the viscous relaxation time drops well below the time for healing over highly damaged areas and (2) inertial e↵ects remain negligible for high deformation rates somehow place a constraint on the minimum value of ↵. Conversely, there is no theoretical upper bound for the value of ↵. However, one can think that for ↵ large, the relaxation time becomes very small at the onset of damage, whatever the damage level (see section 2.2.3). In this case elastic deformations are readily dissipated after damaging, which can a↵ects the mechanical behaviour of the Maxwell-EB model over time. Sensitivity analysis are required in order to determine a range of adequate values for this poorly constrained parameter. This point is discussed in more details in section 2.5.

Small-deformation simulations

Each of the three terms composing the objective time derivative of the Cauchy stress tensor in the constitutive law -the inertial, the advection, and the rotation and deformation term a , see equation (2.1.2) -implies a di↵erent level of mechanical and numerical complexity in the Maxwell-EB model. In developing the numerical scheme, our approach was to introduce each term separately in order to evaluate their relative importance and impact on the simulated mechanical behaviour.

On the one hand, neglecting the advection and a terms allows retaining a Lagrangian scheme, similar to the original EB model (Girard et al., 2010a,b). Without any remeshing of the domain, the model is then suitable for small-deformation simulations only. In this minimal form, its mechanical behaviour can be analyzed in terms of the statistical and scaling properties of the simulated damage and deformation fields. On the other hand, when permanent deformations accumulate over long simulation times, the advection term is no longer negligible and the rotation and deformation of the stress tensor terms become potentially important.

In the remaining of this chapter, the results of small-deformation simulations performed with a highly idealized configuration for the domain geometry, the applied loading and boundary conditions are presented. These will demonstrate that the principal features of the Maxwell-EB model (the spatial heterogeneity, anisotropy and intermittency) naturally emerge from the underlying physics and do not need to be implemented in an ad hoc manner. The treatment of advection and the implementation of rotation and deformation terms is presented in chapter 3.

(b) The small-deformation simulations presented here represent a uniaxial compression experiment and are conceived to be consistent and easily comparable to similar numerical experiments conducted by Girard et al. (2010a) with the original EB model. In terms of the prescribed forcing and boundary conditions, it is perhaps the simplest test case that could be performed.

• n = 0 ⌦ L • n = 0 (a) R L T B u(0, 0) = 0 u y (x, 0) = xy (x, 0) = 0, x > 0 xy (x, L) = 0, x 0 u y (x, L) = U, x 0
The domain is rectangular with dimensions of L 2 ⇥L (see figure 2.7a). 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. Hence Fext = 0. The velocity U is set 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] and to be consistent with a small deformation regime. The Cauchy number (equation 2.30) being very small over both undamaged and highly damaged portions of the material throughout the experiment, the inertial term is neglected in the momentum equation. As advection is also 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 volume are not accounted for and conservation of mass is therefore not imposed. The ice thickness h is constant and set to 1 m.

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 geophysical 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 to accurately represent the limit of strictly elastic deformations within an undamaged material while the numerical scheme converges in a reasonable runtime. The field of cohesion is set once at the beginning of each simulation, as in quenched disorder, by drawing randomly the value of C over each grid element from a uniform distribution spanning 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].

In all simulations presented in section 2.4, the damage parameter ↵ is set to 4 and the characteristic healing time to 10 5 s. These values allows representing both the brittle behaviour and the relaxation of the internal stress within a material with mechanical parameters set as above. All numerical and mechanical parameter values for these uniaxial compression simulations are listed in table 2.2.

In these simulations, the model is made adimensional with respect to the length of the rectangular plate, L, the prescribed velocity U on the top boundary and the undamaged elastic modulus E 0 . The system of equations is given by (2.25) to (2.27), where D ũ D t = 0, h = 1 and Fext = 0 and is solved for the 6 unknowns : ũ (2 components), ˜ (3 components) and d. Finite elements and variational methods are used to solve the time-discretized problem on a Lagrangian grid within the C++ environment RHEOLEF [START_REF] Saramito | E cient C++ finite element computing with Rheolef[END_REF] http://cel.archives-ouvertes.fr/cel-00573970). As cumulative deformations are small, the deformation of the mesh is not calculated and the position of grid nodes, not updated in time. Meshes with triangular elements are built using the Gmsh grid generator [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF] (see figure 2.7b). These are chosen unstructured so that to avoid introducing preferential orientations for the propagation of the damage and be consistent with our fully isotropic rheological framework (see figure 2.7b). The average spatial resolution, x, is set by choosing the number N of elements along the short side of the domain, such that x = L 2N . The model time step is set equal to the characteristic time for damage propagation, which is estimated based on this mean element size, x. The numerical scheme is described in section 2. Table 2.2: Model variables, parameters and domain dimensions for the uniaxial compression experiment of section 2.4.

Numerical scheme : the small-deformation Maxwell-EB model

This section presents the time and space discretizations as well as the numerical algorithm employed to solve the Maxwell-EB system of equations in the limit of small deformations and in the particular case of the uniaxial compression experiment described in 2.3.

The (adimensional) system of equations (2.25) to (2.27), with D ũ D t = 0, h = 1 and Fext = 0 is solved over the closed domain ⌦ 2 R 2 (see figure 2.7a) with boundary partitioned as Simulations are started form rest and from a homogeneous, undamaged state. The strain-driven uniaxial compression problem writes :

@⌦ = T [ L [ B [ R . Let
(P ) : Find ũ, ˜ , d, defined in ⌦⇥]0, 1[, such that r • ˜ = 0 in ⌦⇥]0, 1[, We 0 d ↵ 1 @˜ @ t + ˜ = We 0 d 0 ↵ K : D(ũ) in ⌦⇥]0, 1[, @d @ t = ✓ min  1, ⌃ t 1 ˜ 2 , ⌃ c 1 ˜ 1 q ˜ 2 1 ◆ 1 T d d + 1 T h , 0 < d  1 in ⌦⇥]0, 1[, with initial conditions ũ( t = 0) = 0 in ⌦, ˜ ( t = 0) = 0 in ⌦, d( t = 0) = 1 in ⌦
and boundary conditions

ũy ( t) = 0 on B ⇥]0, 1[, ũy ( t) = 1 on T ⇥]0, 1[, ũ( t) = 0 on L \ B ⇥]0, 1[, ˜ ( t) • n = 0 on L ⇥]0, 1[ and R ⇥]0, 1[, ˜ xy ( t) = 0 on T ⇥]0, 1[ and B ⇥]0, 1[.

Time discretization

Advection and a terms being neglected, non-linearity in the small-deformation Maxwell-EB model arises only through the coupling of the constitutive relationship and of the damage evolution equations. Let us discretize time t such that t n = n t, with t > 0 and n = 0, 1, 2, .... Problem (P ) is solved using a backward Euler (implicit) scheme of order 1 for the momentum and constitutive equations, and a semi-implicit scheme for the damage evolution equation. The discretized system of equations reads:

r • ˜ n+1 = 0, We 0 d n+1 ↵ 1 ˜ n+1 ˜ n ˜ t + ˜ n+1 = We 0 ⇣ d 0 n+1 ⌘ ↵ K : D(ũ n+1 ), d n+1 d n ˜ t = ✓ min  1, ⌃ t ˜ 2 n+1 , ⌃ c ˜ 1 n+1 q ˜ 2 n+1 1 ◆ 1 T d d n + 1 T h , 0 < d n+1  1.
For simplicity, in the following description of the numerical scheme we drop the superscript 0 ˜0 for adimensional variables.

The discretized problem is solved using a fixed-point (FP) algorithm which allows linearizing the system of equation. In this scheme, the field of d is initialized with its value at the previous time step and the momentum and constitutive equations are first solved simultaneously for n+1 and u n+1 (P1). The level of damage is then updated using the field of n+1 (P2). The algorithm iterates between these two steps until the residual of the linearized constitutive equation drops below a chosen tolerance (tol) ensuring the convergence of the solution.

Using the superscript k = 0, 1, 2, ... for sub-iterations, FP the algorithm reads:

For k = 0, let n+1,0 , u n+1,0 , d n+1,0 = ( n , u n , d n ),
For k 0,

• (P1) Find n+1,k+1 and u n+1,k+1 such that r • n+1,k+1 = 0, (2.32)

We 0 d n+1,k ↵ 1 n+1,k+1 n t + n+1,k+1 = We 0 ⇣ d 0 n+1,k ⌘ ↵ K : D(u n+1,k+1
).

(2.33)

and with

u n+1,k+1 y = 0 on B , u n+1,k+1 y = 1 on T , u n+1,k+1 = 0 on L \ B , n+1,k+1 • n = 0 on L and R , n+1,k+1 xy = 0 on T and B .
(2.34)

• (P2) Find d n+1,k+1 , such that 0 < d n+1,k+1  1 and d n+1,k+1 d n t = min " 1, ⌃ t n+1,k+1 2 , ⌃ c n+1,k+1 1 q n+1,k+1 2 # 1 ! 1 T d d n + 1 T h . (2.35) • Stopping criterion : compute res = We 0 d n+1,k+1 ↵ 1 n+1,k+1 n t + n+1,k+1 We 0 ⇣ d 0 n+1,k+1 ⌘ ↵ K : D(u n+1,k+1 ) , If res < tol STOP , and set n+1 , u n+1 , d n+1 = n+1,k+1 , u n+1,k+1 , d n+1,k+1 .
Subproblem P1 is solved by first substituting for n+1,k+1 as a function of u n+1,k+1 n+1,k+1

= K  We 0 d n+1,k ↵ 1 1 t n + We 0 ⇣ d 0 n+1,k ⌘ ↵ ⇤div(u n+1,k+1 ) + 2GD(u n+1,k+1 ) (2.36) where K = h We 0 d n+1,k ↵ 1 1 t + 1 i 1
and ⇤ and G are the dimensionless Lamé coe cients, into the momentum equation (2.32). The internal stress n+1,k+1 is then computed explicitly by substituting the velocity u n+1,k+1 back into (2.36). Subproblem P1 make use of variational methods. The level of damage d n+1,k+1 is solved for exactly by estimating the distance of the known field of n+1,k+1 to the local damage criterion.

Variational formulation and finite element approximations

To obtain the weak form of the momentum equation, let us introduce a test function v 2 V with V, the functional space

V = {v = (v 1 , v 2 ) 2 H 1 (⌦) 2 ; v 2 = 0 on T [ B , v = 0 on L \ R }.
and H 1 the Hilbert space of functions whose gradient is square integrable. Multiplying equation (2.32) by this virtual displacement and integrating over ⌦ :

Z ⌦ r • n+1,k+1 • v dx = 0, 8v 2 V.
Let us now make use of the vectorial Green formula [START_REF] Saramito | Complex fluids. Modelling and algorithms[END_REF], and write the momentum equation in the form

Z ⌦ r • n+1,k+1 • v dx = Z ⌦ n+1,k+1 : D(v) dx + Z @⌦ ( n+1,k+1 • n) • v ds, 8v 2 V. (2.37)
where the boundary integral vanishes over L and R because of the Neumann condition and on T and B because of the Dirichlet conditions (equations 2.34).

Substituting for n+1,k+1 (equation 2.36) into (2.37) leads to a linear expression of the momentum equation in terms of u n+1,k+1 :

Z ⌦ KWe 0  d n+1,k ↵ 1 1 t n : D(v) + ⇣ d 0 n+1,k ⌘ ↵ ⇤div(u n+1,k+1 ) • div(v) + 2GD(u n+1,k+1 ) : D(v) dx = 0
The variational formulation of this subproblem therefore reads:

(P 1.1) : Find u n+1,k+1 2 V such that a(u n+1,k+1 , v) = l (n+1,k+1) (v), 8v 2 V,
where the bilinear form a() and linear form l() are defined for all u, v 2

H 1 (⌦) 2 as a(u, v) = Z ⌦ KWe 0 d 0n+1,k ↵ [⇤div(u) • div(v) + 2GD(u) : D(v)] dx l (n+1,k+1) (v) = Z ⌦ KWe 0 d n+1,k ↵ 1 1 t n : D(v) dx
and where d n+1,k , d 0 n+1,k and n are known quantities, calculated at the n th iteration and k th subiteration.

The computation of n+1,k+1 (P1.2), is done by substituting u n+1,k+1 into equation (2.36). The level of damage d n+1,k+1 is computed explicitly as well from equation (2.35) and by comparing the field of n+1,k+1

to the local damage criteria (i.e., to ⌃ c and ⌃ t ).

The 2-dimensional domain ⌦ is discretized in space as ⌦ h and its boundary as h using a family of triangulations T h made of triangular finite elements, with h > 0 is the discretization parameter representing the greatest diameter of a triangle in T h . The subscript h is used throughout to denote the finite element approximations a continuum quantity. We introduce 1. X h , the finite element (FE) space for the velocity field u h and by V h the space for the associated test function v h :

X h = {v h 2 H 1 (⌦) 2 ; v h/K 2 (P k ) 2 , 8K 2 T h } V h = X h \ V
with k 1, the degree of the polynomial functions on each mesh triangle.

2. T h , the FE space for the field of internal stress h

T h = {⌧ h 2 L 2 (⌦) 2⇥2 ; ⌧ h = ⌧ T h and ⌧ h/K 2 (P k 1 ) 2⇥2 , 8K 2 T h , } (2.38)
As the stress tensor is a function of the velocity gradient, the approximations in the space T h is of degree k 1 and discontinuous at inter-element boundaries.

3. S h the space for the level of damage, mechanical parameters and other scalar quantities.

S h = {' h 2 L 2 (⌦); ' h/K 2 P k 1 , 8K 2 T h } (2.39)
which, in the case of d h , is bounded and defined on the interval [0, 1].

In the Maxwell-EB framework, the degree k of the polynomial functions for each of these finite element spaces is chosen based on the degree of the polynomial approximations for the field of cohesion, C h , in which disorder is introduced to represent the material's heterogeneity. As mentioned in previous sections, in the case of sea ice the correlation length of heterogeneities, ⇠, is much smaller than the mean model resolution. Hence disorder is introduced at the smallest resolved scale: the scale of the mesh element, x. This sets the degree of the polynomial approximations for the field of cohesion to 0, i.e., constant by element. The FE approximations for scalar mechanical parameters, for the level of damage d h , the distance to the damage criterion and the local damage criteria (⌃ ch and ⌃ th ), all defined in S h , are therefore necessarily P 0 . In turn, the degree of polynomial approximations in T h is also 0, which makes V h for the velocity field P 1 , piecewise linear. In finite element models, higher levels of precision in the computed solutions are usually obtained by increasing the degree of the polynomial functions of the FE spaces. With the degree of the FE approximation fixed in the current implementation of the Maxwell-EB rheology, the structure and convergence of the model solution is analyzed in section 2.4 only in terms of the spatial resolution. This point is further discussed in details in chapter 3.

Using these temporal and spatial discretizations, the strain driven uniaxial compression problem becomes:

(P ) h : Initialization (n = 0) u n h = 0 in ⌦ h , n h = 0 in ⌦ h , d n h = 1 in ⌦ h ,
For n 0

• For k = 0, ⇣ n+1,0 h , u n+1,0 h , d n+1,0 h ⌘ = ( n h , u n h , d n h ).
• For k 0

(P 1.1) h : With n h and d n+1,k h known, find u n+1,k+1 h 2 V h such that a h (u n+1,k+1 h , v h ) = l (n+1,k+1) h (v h ) 8v h 2 V h
with the bilinear and linear forms a h and l h defined for all u h , v h 2 V h by

a h (u h , v h ) = Z ⌦ h KWe 0 ⇣ d 0 n+1,k h ⌘ ↵ [⇤div(u h ) • div(v h ) + 2GD(u h ) : D(v h )] dx l (n+1,k+1) h (v h ) = Z ⌦ h KWe 0 ⇣ d n+1,k h ⌘ ↵ 1 1 t n h : D(v h ) dx (P 1.2) h : With u n+1,k+1 h , n h and d n+1,k h known, compute n+1,k+1 h 2 T h explicitly from (2.36). (P 2) h : With n+1,k+1 h known, compute d n+1,k+1 h 2 S h explicitly from (2.

35).

Compute

res h = We 0 ⇣ d n+1,k+1 h ⌘ ↵ 1 n+1,k+1 h n h t + n+1,k+1 h We 0 ⇣ d 0 n+1,k+1 h ⌘ ↵ K : "(u n+1,k+1 h ) If res h < tol STOP, Set n+1 h , u n+1 h , d n+1 h = ⇣ n+1,k+1 h , u n+1,k+1 h , d n+1,k+1 h ⌘ .
To ensure that the cumulative deformation of the simulated material is small enough so that the small deformation approximation is valid, the uniaxial compression simulations are run up to a maximum of n 20% ⇥ 1 ˜ t 1 N model time steps (i.e., until the total, cumulative deformation prescribed on the top edge of the domain is of maximum 20% of the size of an average mesh element).

At each time step, the model exits the fixed-point iteration (in k) when the residual, res h , evaluated as the L 2 norm | • | of the sum of the residuals of the constitutive equation on ⌦ h , drops below a chosen tolerance. The log-linear plot of figure 2.8b shows the evolution of res h as a function of the number of fixed-point iterations k at di↵erent time steps (identified by di↵erent colour curves) of a uniaxial compression simulation with N = 10. Figure 2.8a shows the evolution of the macroscopic stress (dashed-dotted curve) and damage rate (solid grey curve) for this simulation. During periods of stress-buildup and low damaging activity, the residual drops to the machine precision within a few iterations of the fixed point algorithm. The convergence is slower and a larger number of sub-iterations are required to reach convergence when damage rates are high, such as during macro-rupture events. As these events are highly localized in time, the fixed-point numerical scheme remains computationally e cient for relatively small systems. In the uniaxial compression simulations, the tolerance on res h (tol) is fixed to 10 12 .

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 spatial 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 ice motion product. They estimated deformation rate invariants (shear, divergence and total deformation rates) at di↵erent spatial scales by averaging the components of the 3-days strain rate tensor over square boxes of length 10 km  l  1000 km. Doing so, they obtained a power-law relationship between the total deformation rate < "tot > l and the corresponding averaging scale l of the form

< "tot > l ⇠ l (2.40)
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]. With a higher temporal resolution than the RGPS motion products, ice buoys positions 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 (2.41) 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][START_REF] Rampal | neXtSIM: a new Lagrangian 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[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 di↵erent 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. In this section and the remaining of this chapter, results are expressed in terms of adimensional variables and parameters. However for the sake of simplicity, we drop the superscript 0 ˜0 for all quantities. 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 di↵use to localised damage through elastic interaction[END_REF]Girard et al., 2010a). The Maxwell-EB model simulates 1. a 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 figure 2.9 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 figure 2.9 b2).

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. In all four cases, the macroscopic behaviour is not linear in the early stage of the stress build-up phase, consistent with the significant contribution of viscous-like dissipation within the highly damaged material. The behaviour then becomes nearly linear as the material heals. The associated slope of the macroscopic stress-strain curve indicates an "e↵ective" elastic modulus E < E 0 , i.e., lower for the partially damaged than for the initially undamaged 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 locations of the first damage events are the same at all resolutions, as shown by the maps of the instantaneous level of damage d near the onset of damaging (figure 2.9 b1). However, soon after these first failure events, model solutions do not converge (figure 2.9 b2-4) and fractures form with a shape and orientation di↵ering 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 di↵erent redistribution of the stress between neighbouring elements at di↵erent 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] Sornette | Statistical physics of fault patterns self-organized by repeated earthquakes[END_REF][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 2.4.2) represented on figure 2.9

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 di↵use 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 law.

Heterogeneity

As shown in the previous section, when simulations are started from an undamaged state, the simulated mechanical behaviour of the material is intrinsically di↵erent between the first and subsequent loading and damaging cycles. The path to the first rupture in "irreversible damage" (i.e., models without healing) elastobrittle models has already been investigated in depth (e.g. [START_REF] Tang | Numerical simulation of progressive rock failure and associated seismicity[END_REF][START_REF] Amitrano | From di↵use to localised damage through elastic interaction[END_REF]Girard et al., 2010a). Hence in the following 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 Marsan et al. ( 2004) and estimate deformation rates over two orders of magnitude in space scales using a coarse-graining procedure. The calculation is similar to that performed by Girard et al. (2010a). This analysis uses 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 invariants and total deformation rates are first computed at the largest scale, l = L, by averaging the components of the strain rate tensor over the whole domain. Then, deformation rates are estimated at di↵erent space scales l = L 2n with 1  n  N/2 by dividing the domain into 2 ⇥ n 2 square boxes of equal size l, the smallest scale therefore corresponding approximately to two times the width of a mesh grid element. At each scale l, the components of the strain rate tensor are first averaged spatially over each individual box by finding all elements that have their centre lying inside that box, with their contribution weighted by their actual area. The shear, divergence and total deformation rates are computed over each box as

"div = @u @x + @v @y , (2.42) "shear =  @u @x @v @y 2 +  @u @y + @v @x 2 ! 1 2 , (2.43) "tot = "2 div + "2 shear 1 2 , (2.44)
with the overline indicating the space-averaged components. The three invariants are then averaged over all boxes of same size. The mean total deformation rates is denoted by < "tot > l .

The dependence of the deformation rates on the spatial scale of observation is investigated at di↵erent stages of the healing-damaging cycle. Figure 2.10 (a and b) 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 figure 2.10(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 equation (2.40) 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. 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 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 di↵erent 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 di↵erence between the present results and that of Girard et al. (2010a), who performed a similar analysis for the path to the first macro-rupture simulated with their irreversible-damage EB model, 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 spatial 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 cuto↵ 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. 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 figure 2.11a).

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 (figure 2.11b) 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 slowly dissipates all of the input loading. Here, creep is defined by analogy with creeping faults in the Earth crust showing a continuous micro-earthquake activity [START_REF] Amelung | Earthquake scaling laws for creeping and non-creeping faults[END_REF], not in terms of the dislocation-mediated creep of ice. This period of creep is interrupted by a macro-failure event as evident by the abrupt and complete cessation of damaging, and the formation of a new fault outside the system (last panel of figure 2.11b). The partitioning between the two types of mechanical behaviour depends on the rate of healing and of viscous dissipation set in the model through through the parameters t h and ↵.

Sensitivity analysis on the value of ↵ are presented in section 2.5. 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 the damage rate.

Time series are truncated by removing the first macro-rupture event and the squared Fourier coe cients are averaged over 5 realizations of the compression experiment initialized with di↵erent fields of C over domains with N = 40. Figure 2.12(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 a↵ect the shape of the PSD discussed below.

At low frequencies, the PSD is almost flat, suggesting that the damage rate 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 that increments of the damage rate are uncorrelated. At the highest frequencies, > 2, indicating that the damage rate is correlated in time and increments of the damage rate are anti-correlated. The break in the slope occurs around f = 10 6 , a frequency Power Spectral Density of the damage rate! (cumulative number of damaged elements !

x ! distance to damage criteria) 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 figure 2.11(a).

T 1 h 2 T ⇥ c L > 2 (a) PSD, (1 d crit ) f t 10 -6 10 -5 10 -4 10 -3 < tot > t 10 0 10 1 T h (b) 1 2 L T ⇥ c
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 2.4.2 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 2.12(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 di↵erent 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. The fact that the average of the 20 realizations of the coarse-graining calculation is rather smooth indicates that these cycles are not associated to a single frequency but instead span a range of frequencies, consistent with the PDS of the damage rate signal. 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. For t . 10 5 , no scaling is observed. This time corresponds to the average time of macroscopic failure, as measured from fields of the deformation rate. It is about one order of magnitude larger than the minimum time of propagation of damage across the width of the domain (for which a break in the scaling is observed on the PSD of the damage rate), consistent with the fact that deformation evolves as a function of the local state of damage d, not of the damage rate (⇠ 1 d crit ). The average of the 20 realizations of the coarse graining calculation also indicates a change of scaling exponent (equation 2.41) associated with the period of stress build-up/relaxation cycles (red dotted line). Beyond this period, is much reduced but the scaling is not entirely lost, in agreement with the point discussed at the beginning of this section that creeping faults can remain activated over an extended period of time.

Sensitivity analyses

In this section, we investigate the dependence of the Maxwell-EB model behaviour on the value of the damage parameter ↵, which is the only truly ad-hoc parameter of the Maxwell-EB rheological framework and perhaps the least constrained, both physically and numerically. Based on the coupling between ⌘, E, and the level of damage d described in section 2.2.3, we can distinguish and comment on the following values of this parameter.

• For ↵ = 0, the apparent viscosity of the material, ⌘ is independent of d and hence constant in both space and time. If the initial undamaged value of ⌘ 0 is chosen very large, such that 1 ⌘ 0 ! 0, the viscous dissipation term in the constitutive equation (2.6) vanishes and all deformations are elastic.

• For ↵ = 1, the apparent viscosity decreases with d in the same proportions as the material's elastic modulus. Hence the relaxation time is constant in both space and time, implying that the rate of dissipation of the internal stress in permanent deformations does not depend on the level of damage.

• For ↵ > 1, the elastic modulus, apparent viscosity and relaxation time all decrease with increasing level of damage, consistent with a degradation in both the mechanical strength and the capacity of the material to retain the memory of elastic deformations.

• In the limit of ↵ ! 1, the relaxation time drops to zero at the onset of damaging and the local memory of elastic deformations is lost.

Here we analyze a set of uniaxial compression experiments in which the value of the damage parameter is varied over the interval [0,100]. The value of the characteristic healing time and of all other mechanical parameters are identical to that used in section 2.3 (see table 2.2). The domain is as shown in figure 2.7 and the mean mesh resolution is x = 1 N with N = 40. Simulations are started from an undamaged state and initialized with the same field of cohesion.

We first briefly investigate the cases of ↵ = 0 and ↵ = 1, which both violate the principle of stress relaxation in the context of the Maxwell-EB framework. With ↵ = 0, the elastic modulus decreases with damaging but the apparent viscosity remains constant. Hence increases with damaging. As 0 is set very large over the initially undamaged material to represent the limit of pure elasticity (see section 2.2.4), there is virtually no viscous dissipation and stress relaxation occurs only through damaging. No abrupt drop of the macroscopic stress is observed after the onset of damaging in this case. Instead, the damage and deformation slowly cluster along linear features (see 2.15) and the macroscopic stress stabilize after some time (see figure 2.13a, solid black curves), suggesting a balance is reached between the rate of loading and the rate of damaging. This is can be explained by the fact that in the Maxwell-EB model, the part of the stress relaxation due to elastic deformation is "prescribed" by calculating the change in level of damage of any given element according to its exact distance to the damage criterion, d crit (see section 2.2.3). Without viscous dissipation, all damaged elements have their state of stress lying just on the damage envelope and all undamaged elements, just inside the envelope, near the maximum uniaxial compressive strength (⌃ c ) state at which the first damaging events preferentially occur under uniaxial compression forcing conditions (see figure 2.16a). As the maximum value of , or equivalently the value of We, is unbounded in this case, numerical instabilities arises as the level of damage drops over the domain and the simulations stop after a short runtime.

For ↵ = 1, viscous stress relaxation is e↵ective. However, as the rate at which it operates is independent of the level of damage, the model behaviour depends strongly on the chosen value of 0 , which sets the elastic versus viscous character of the deformation. Here we vary 0 , or equivalently, the adimensional Weissenberg number (We 0 = 0 U L ) for the simulations, over four orders of magnitudes. Figures 2.13a and 2.13b show respectively the evolution of the macroscopic stress and of the damage rate for We 0 = 0.05 (dotted lines), which is the value otherwise used in all simulations presented in this chapter and which approaches the limit of pure elasticity for the undamaged material, We 0 = 5 • 10 4 (dashed-dotted lines) and We 0 = 5 • 10 6 (dashed lines). For the highest value of We 0 , the behaviour is similar to the ↵ = 0 case : the input loading is balanced by the stress relaxation associated to small, elastic deformations and states of stress are always near-critical.

Damage does concentrate along linear features (see figure 2.15). However the structure of these features does not evolve in time and no main fault emerges that dominates the deformation. When lowering the value of We 0 and allowing for more significant viscous dissipation (We 0 = 5 • 10 4 ), deformation localizes along one main fault (see 2.15). However this localization happens very slowly, not in an abrupt, brittle-like manner, and the damage rate still shows no significant fluctuations over time. For even smaller values of We 0 , the behaviour is essentially viscous-like: no damaging occur and all of the applied loading is readily dissipated into permanent deformations.

For values of ↵ greater than 1 (and We 0 = 0.05), time series of the stress and damage rate reveal a range of di↵erent mechanical behaviours (see figures 2.14a and 2.14b). The case of ↵ = 2 (darkest blue curves) is similar to ↵ = 0 and ↵ = 1 and We 0 = 0.05 : no brutal stress relaxation is observed after the initial elastic loading and damage rates stabilize after a long period of softening.

For ↵ = 3 (blue curves), the macroscopic stress shows a small but rapid drop associated to the propagation of a first system of faults across the domain. The stress subsequently builds up and stabilizes. States of stress are scattered near the damage envelope (see figure 2.16b). The damage rate time series shows low amplitudes and high-frequency fluctuations of the damaging (see figure 2.14c), which concentrates along the initially activated faults. Some small macroscopic drops in the stress are observed that are associated with higher damaging activity and to the reactivation of a portion of the same system of faults. Figures 2.14a and 2.14b represent only 10% of the time series for this simulation. During the entire simulation, the same feature remained activated, slowly dissipating the applied loading through creep-like deformation.

The case of ↵ = 4 (clear blue curves) was described in section 2.4.3. The mechanical behaviour is partionned between extended periods of creep-like deformation and cyclic, sharp stress relaxations associated with the activation of new features. The behaviour is similar for ↵ = 5 (green curves) and ↵ = 6 (not shown), suggesting that there is a range of values of ↵ for which both types of deformation can coexist.

When increasing ↵ further, episodes of creep-like deformation disappear. After each macro-rupture event, states of stress drop well inside the damage envelope (see figure 2.16c). Damaging clusters in time (see 2.14b, yellow curve), and somewhat localizes in space (see figure 2.15). For ↵ = 10, the same fault is formed again at each stress-build up and relaxation cycle, equivalent to re-starting the experiment from its initial undamaged state, thereby suggesting the total suppression of the memory of previous damaging events within one stressbuild up cycle. At ↵ = 100 (orange curves), there is almost no precursory damaging activity to the propagation of macro-ruptures (see figure 2.14b and 2.14c). The stress drops to near-zero values at each failure event. The fluctuations in the stress near peak values are due to the fact that the numerical scheme does not converge in the set maximum value of iterations of the fixed-point algorithm.

Figure 2.17 shows power spectral densities calculated for the damage rate time series in the case of ↵ = 2, 3, 4, 5, 10 and 100 (time series were truncated by removing the first macro-rupture event). An averaging window of 15 frequencies centred on any given frequency is used. For ↵ = 3, the PSD shows no distinct peak. A power law decay with a slope of ⇡ 2 is observed for frequencies larger than that associated with the prescribed time of healing, T h . As ↵ is increased beyond 3, a characteristic frequency appears that is associated with the cycles of stress build-ups and relaxation. The PSD flattens starting at the lowest frequencies, consistent with the loss of the long-term memory of damaging events. At ↵ = 100, the PSD is flat, except around a peak centred on the frequency of healing (which appears broader due to the failure of the numerical scheme near the peak macroscopic stress), indicating a complete loss of temporal correlations of the damage rate, at all time scales.

To summarize, by varying only the value of the damage parameter the Maxwell-EB model can represent widely di↵erent mechanical behaviours, even under homogeneous and highly idealized forcing conditions. The cases of ↵ = 0 and ↵ = 1 are pathological, as expected from the misrepresentation of the dependence of the relaxation time for the stress on the level of damage, central to the rheological framework. With ↵ = 3, the scale invariance observed all time scales below the prescribed healing time is consistent with a marginally-stable state and the maintenance of a stable structure and orientation of the fault by short-term damage fluctuations.

In this case, a scaling analysis of the deformation in space would likely reveal stable statistical properties, i.e.

a spatial scaling exponent , (see section 2.4.2) that is constant in time. For ↵ 10, the model is essentially elasto-plastic. The memory of elastic deformations is lost just after a few model time steps. Deformation over damaged elements becomes readily permanent, whatever their level of damage. The macroscopic behaviour is characterized by a single fault pattern and a single damaging frequency, corresponding to the healing frequency, and is hence completely predictable. In this elasto-plastic limit, the model could reproduce the activation of faults with di↵erent shapes and orientations and over a range of frequencies, provided a non-homogeneous and su ciently complex forcing is used. However, any intermittency in the system would then be inherited from the applied forcing, not from the mechanical framework itself. In the context of sea ice modelling, this implies that the model would not reproduce the extreme fluctuations and long-range temporal correlations that characterize the intermittency of deformation of the ice pack [START_REF] Rampal | Scaling properties of sea ice deformation from buoy dispersion analysis[END_REF][START_REF] Weiss | Intermittency of principal stress directions within Arctic sea ice[END_REF]. One question remains unanswered at this point. Consistent with the work of [START_REF] Frederiksen | Numerical modelling of strain localisation during extension of the continental lithosphere[END_REF] who used a strain-dependent apparent viscosity to model strain localization during lithospheric extension with an elasto-visco-plastic mechanical framework, we note that increasing the rate of change of the viscosity ⌘ with the level of damage, i.e., increasing ↵, leads to a somewhat higher degree of localization of the deformation.

↵ = 0 ↵ = 1 0 = 10 5 s 0 = 10 7 s ↵ = 3 ↵ = 10
However, figure 2.15 (↵ = 10) suggests that in the limit of ↵ large, damage concentrates along a single feature and is otherwise homogeneously distributed over the domain. Hence we can wonder if spatial scale-invariance is still maintained in the elasto-plastic limit, or if characteristic scales eventually arise for the deformation.

This point remains to be clarified through appropriate spatial scaling analyses.

Finally, it is important to note that the range of values of ↵ for which a creep-like, elastic-brittle or elastic-plastic dynamics is simulated with the Maxwell-EB model is expected to depend on the value of other parameters such as the healing time, t h , and the initial, undamaged relaxation time for the stress 0 . It also depends on the magnitude, temporal and spatial structure of the applied forcing and, possibly, on the level of heterogeneity initially introduced in the model (i.e., the width of the distribution of C). Hence, in order to determine the values of ↵ that are appropriate to represent either kind of mechanical behaviour under specific conditions, or for a specific material, sensitivity analysis should be performed with the model. 3.1 From small to large deformations in the Maxwell-EB framework

From the continuum mechanics point of view, sea ice is both (1) a solid whose mechanical properties depend on the history of its deformation and (2) a material experiencing high strains. In this sense, modelling the deformation and drift of the ice pack lies between a solid mechanics (small deformations) and a fluid dynamics (large deformations) problem. The numerical experiments carried in chapter 2 suggest that in the smalldeformation regime, the Maxwell-EB rheology succeeds in representing the mechanical behaviour of sea ice and, in particular, in reproducing the linear and highly localized features that concentrate its deformation. The advection scheme that handles the large deformations along these features once formed and their transport with the mean flow must allow retaining the strong gradients and the high degree of heterogeneity arising within all simulated fields.

On the one hand, the most common approach to model the deformation of continuum solids is to cast the equations of motion in Lagrangian form. In this case computational points, i.e, mesh grid nodes, are tied to the material and therefore keep a historical record of its mechanical properties. Strain and strain rates are computed from the deformation of the mesh. Hence if deformations within the material become large, restructuring of the grid is necessary to preserve a homogeneous sampling and maintain a numerically workable grid. Remeshing inevitably results in interpolation-related di↵usion of the recorded information, which can be somewhat minimized if one is willing to pay the computational price and call the remeshing scheme frequently.

Remeshing also introduces some inaccuracy in the model solution, as the restructured fields might not respect the local force balance or yield criterion.

On the other hand, typical fluid flow problems make use of an Eulerian advection scheme and fixed model grid, which allows dealing with arbitrary large strain rates, best defined in terms of the velocity field. This approach is less suitable however if the mechanical properties of the material are allowed to evolve in both space and time, as transport implies interpolating the properties onto the grid point positions where velocity is estimated, which lead to numerical di↵usion of gradients and to a deterioration of the material's history over time.

Most climate models including sea ice dynamics use an Eulerian-type advection scheme. As these are based on a viscous-plastic rheology [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF][START_REF] Hunke | An elastic-viscous-plastic model for sea ice dynamics[END_REF][START_REF] Tsamados | Impact of a new anisotropic rheology on simulations of arctic sea ice[END_REF], their mechanical framework is indeed compatible with a fluid-type numerical approach. Besides, Eulerian sea ice models present the advantage of being easily coupled with an ocean and atmospheric components.

One exception in the sea ice community is the recently developed NeXtSIM model [START_REF] Bouillon | Presentation of the dynamical core of neXtSIM a new sea ice model[END_REF][START_REF] Rampal | neXtSIM: a new Lagrangian sea ice model[END_REF], which is based on the solid mechanics-type Elasto-Brittle rheological framework of Girard et al. (2010b) and uses a Lagrangian scheme together with a remeshing method. With their remeshing criterion (based on a threshold for the smallest angle of a triangular mesh element), the model grid needs to be updated on average every hour (in simulation time). So far this model has been used for relatively short (10 days to 2 months) simulations over the Arctic. Even if adaptive remeshing techniques are becoming more and more computationally e cient, these might still not be practical for long-term, global climate simulations.

Between the Lagrangian and Eulerian approaches, alternative methods have been developed to e ciently handle large deformations within geomaterials in continuum models. [START_REF] Braun | Dynamical lagrangian remeshing (dlr): A new algorithm for solving large strain deformation problems and its application to fault-propagation folding[END_REF] have developed a Dynamical Lagrangian Remeshing (DLR) scheme, based on a Lagrangian finite element method and the Delaunay triangulation, in which the mesh was defined using the same material particles throughout the entire simulation but the connectivity between the particles was updated at the end of each time step to accommodate the material's deformation. The DLR method was shown to successfully simulate fault propagation and the large deformations in an elasto-plastic model of the Earth crust. However, one downside of this method is that while it eliminates the need to interpolate scalar values associated with the mesh nodes, some interpolation is still implied for values defined over mesh elements, like stress tensors. Particle-in-cell (PIC, [START_REF] Evans | The particle-in-cell method for hydrodynamic calculations[END_REF] methods combine the Eulerian and Lagrangian approaches by tying material properties to non-interacting particles that can move through an Eulerian mesh grid with fixed nodal points. Their main advantage is that they allow casting the equations of motion in the frame of reference of the moving particles, which limits the numerical di↵usion of gradients. These methods, originally designed for hydrodynamics problems, have been used in sea ice models based on the Hibler viscous-plastic rheology, in particular to advect the ice thickness and compactness (i.e, concentration) and take advantage of a Lagrangian representation of the ice edge (e.g., [START_REF] Flato | A particle?in?cell sea?ice model[END_REF].

PIC approaches were later extended to solids with history-dependent variables [START_REF] Sulsky | A particle method for history-dependent materials[END_REF][START_REF] Sulsky | Application of a particle-in-cell method to solid mechanics[END_REF] by allowing for the evaluation of the constitutive equation on the material points rather than on the grid nodes. These so-called Material Point Methods (MPM) were implemented in the context of sea ice modelling by [START_REF] Sulsky | Using the material-point method to model sea ice dynamics[END_REF], which demonstrated their robustness for the transport of material properties and the representation of large deformation within a continuum material. [START_REF] Moresi | A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials[END_REF] developed a variant of the MPM approach for viscoelastic material, called the Lagrangian Integration Point -finite element method, with the aim of solving convection problems in geodynamics. A general limitation of these methods however arises for certain configurations of the flow, for instance under local convergence or divergence, which can leave some grid elements empty of material points and others significantly overepresented. Then, the density of points need to be made more uniform over the computational domain by somehow splitting, merging, creating or destroying particles [START_REF] Moresi | A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials[END_REF]. In the case of sea ice, this can be dealt with through thermodynamics considerations: material points can be created through refreezing within diverging leads and destroyed where the ice melts [START_REF] Sulsky | Using the material-point method to model sea ice dynamics[END_REF]. However, if thermodynamic processes are not (or not yet) accounted for, other, possibly non-trivial, artefacts need to be implemented that deal with this issue.

Here, we base our choice of an advection scheme for the Maxwell-EB model on its flexibility and its capacity in handling the advection of strong gradients while keeping in mind an eventual coupling of the rheological framework to an ocean and atmospheric component. We cast the equations of motion in the Eulerian frame and write their variational formulation using discontinuous Galerkin methods. Doing so, we therefore retain a finite element framework and trade continuous trial and test functions for discrete polynomial approximations and the associated numerical fluxes at element interfaces.

Discontinuous Galerkin methods have experienced a significant expansion in recent years owing to their many appealing properties [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]. First, as all finite element-based approaches, they allow dealing with complex geometries and boundary conditions. Second, these methods, which can be viewed as a higher order upwind finite volume scheme in which the solution is approximated over each mesh element as a polynomial function instead of being piecewise constant, are known to be very stable. This makes them wellsuited for problems where advection is dominant. DG methods were also shown to be e cient for high Reynolds flows [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF][START_REF] Kanschat | Discontinuous Galerkin Methods for Viscous Incompressible Flow[END_REF][START_REF] Saramito | E cient C++ finite element computing with Rheolef[END_REF] and for discontinuous approximation of stress tensors in viscoelastic [START_REF] Saramito | On a modified non-singular log-conformation formulation for johnson-segalman viscoelastic fluids[END_REF] and elasto-visco-plastic [START_REF] Cheddadi | A new operator splitting algorithm for elastoviscoplastic flow problems[END_REF] flow problems.

In the context of ocean modelling, continuous finite elements and finite volume schemes (with upwinding) were found to generate unphysical oscillations and have problems handling steep gradients. Their discontinuous and non-conforming counterparts were found to perform better [START_REF] Hanert | Advection schemes for unstructured grid ocean modelling[END_REF]. Finally, discontinuous

Galerkin methods also represent a practical framework for the development of high-order accurate numerical schemes and are known to be readily parallelizable, which in the context of modelling sea ice on large scales, can eventually prove to be a highly valuable property.

As we introduce the full objective internal stress tensor to deal with the large deformations of a Maxwell-EB material, we also become concerned with numerical constraints a↵ecting the representation of its elastic character. In practice, purely elastic deformations cannot be exactly represented by standard viscoelastic models due to the failure of their numerical schemes when the value of the adimensional number characterizing the elastic versus the viscous nature of the flow, the Weissenberg number We, becomes large, a limitation widely known and referred to in the applied mathematics literature as the high Weissenberg number problem [START_REF] Keunings | On the high Weissenberg number problem[END_REF]. The origin of these numerical instabilities has been traced to the rotation and deformation terms in the transport equation for the stress tensor. [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] attributed them to the inaccurate representation of the exponential growth of the stress tensor by polynomial FE approximations.

In the rheological framework presented here, the deformation of an undamaged medium, for instance an unfractured ice pack, is expected to be strictly elastic. The numerical limitation on the value of the Weissenber number however can lead to unphysical viscous dissipation within the material. Hence in developing the Maxwell-EB framework, the high Weissenberg number problem represents "more than a mathematical issue for purist" [START_REF] Keunings | On the high Weissenberg number problem[END_REF], and its impact on the range of applications of the rheological model must be assessed.

For the first tests of the large-deformations Maxwell-EB model, we choose to tackle a Couette flow problem.

The annular Couette setting allows (1) to test the model in a fundamentally di↵erent configuration than that used in the small-deformations simulations presented in chapter 2, i.e., in simple shear as opposed to uniaxial compression, (2) to maintain an idealized, easily interpretable simulation setup while having the possibility of applying arbitrarily large strains and (3) to compare the model's solution to known analytical solutions of standard viscoelastic models as well as to the results of a laboratory experiment performed on ice.

The chapter is structured as follow. Section 3.2 presents the discontinuous Galerkin approximation of the advection term for transported scalar quantities in the Maxwell-EB framework and assesses the capacity of the DG method in handling the advection of the strong gradients and discontinuous fields simulated by the model.

Section 3.3 presents the DG approximation of the objective derivative for the stress tensor in the Maxwell-EB constitutive law in the context of a homogeneous, undamaged material and evaluates the maximum allowed value of the Weissenberg number that ensures the convergence of the numerical scheme. Section 3.4 presents the numerical scheme for the full, large-deformation dynamical Maxwell-EB model, compares the results of a numerical annular Couette flow experiment to that of a similar laboratory experiment and addresses some numerical and physical limitations of the rheological framework.

Advection and di↵usion

In this section, we evaluate the ability of a discontinuous Galerkin-based numerical scheme of handling the x . We consider a passive scalar field transported by this velocity field. With the time t defined over the interval [0, 1], the problem reads:

(P ) : Find in ⌦ ⇥ [0, 1], such that @ @t + (u • r) = 0. (3.2)
We compare the accuracy of the numerical scheme between three fields of with di↵erent levels of regularity.

1. In the first case, is a cosine function 2 H 2 (⌦) \ C 1 ( ⌦) 2 given by

1 (✓) = 8 < : 1 2 (cos(4✓) + 1) on ✓ 2 [ ⇡ 4 , ⇡ 4 ] 0 on ✓ / 2 [ ⇡ 4 , ⇡ 4 ] (3.3)
t > 0 and n = 0, 1, 2, ... . The domain ⌦ is discretized in space using a two-dimensional mesh, T h , made of triangles. The transport equation is written in weak form by introducing a scalar test function ' h for the field and the discontinuous finite element space

S h = {' h 2 L 2 (⌦); ' h|K 2 P k , 8K 2 T h } (3.5)
with discontinuous elements P k of degree k > 0. We obtain

Z ⌦ h 1 t n+1 Z ⌦ h 1 t n + Z ⌦ h (u • r h ) n+1 h ' h dx = 0.
Broken gradients are involved when computing the gradient of the discretized field h . The discontinuous Galerkin approximation of the advection term reads

X K Z K (u • r h ) n+1 h ' h dx + X K X S2@K Z S ✓ # 2 |u • n| [[ n+1 h ]][[' h ]] (u • n)[[ n+1 h ]]{{' h }} ◆ ds, 8' h 2 S h
where # = 1 for an upwinding scheme, [[•]] and { {•} } denote respectively the jump and the average of discontinuous sets across inter-element boundaries and n is the normal on an oriented element side (Di Pietro and Ern (2012), chapter 2, Saramito (2013), chapter 1). Since in this case there is no flux of in or out of the domain and hence, no contribution from the domain boundaries, the second sum representing the jumps in across all element sides S can be expressed alternatively as a single sum over the sides belonging to the set

L (i)
h of internal sides of the mesh T h . With these definitions and using a backward scheme of order 1, the discretized problem (P ) n,h reads:

For n 0 and with n h known, find

n+1 h 2 S h such that b h ( n+1 h , ' h ) = o n+1 h (' h ), 8' h 2 S h
with the bilinear form b h () and linear form o h () given by

b h ( n+1 h , ' h ) = 1 t Z ⌦ h n+1 h ' h dx + X K Z K (u • r) n+1 h ' h dx + X S2L (i) h Z S ✓ 1 2 |u • n| [[ n+1 h ]][[' h ]] (u • n)[[ n+1 h ]]{{' h }} ◆ ds o n+1 h (' h ) = 1 t Z ⌦ h n h ' h dx
For a backward scheme of order 2, and with n h and n 1

h known, b h ( n+1 h , ' h ) = 3 2 t Z ⌦ h n+1 h ' h dx + X K Z K (u • r) n+1 h ' h dx + X S2L (i) h Z S ✓ 1 2 |u • n| [[ n+1 h ]][[' h ]] (u • n)[[ n+1 h ]]{{' h }} ◆ ds, o n+1 h (' h ) = 4 2 t Z ⌦ h n h ' h 1 2 t Z ⌦ h n 1 h ' h dx
and in the case of the order 3, with n h , n 1 (Süli and Mayers, 2003, p. 349).

h and n 2 h known, b h ( n+1 h , ' h ) = 11 6 t Z ⌦ h n+1 h ' h dx + X K Z K (u • r) n+1 h ' h dx + X S2L (i) h Z S ✓ 1 2 |u • n| [[ n+1 h ]][[' h ]] (u • n)[[ n+1 h ]]{{' h }} ◆ ds, o n+1 h (' h ) = 18 6 t Z ⌦ h n h ' h 9 6 t Z ⌦ h n 1 h ' h dx + 2 6 t Z ⌦ h n 2 h ' h dx.
Model simulations are run with (adimensional) time steps t 2 { 1 10 , 1 100 , 1 1000 , 1 10000 , 1 100000 }. Unstructured meshes with mean (adimensional) element size x = 1 N are used, with N 2 {10, 20, 40, 80}. Figure 3.1b shows an example of a mesh with x = 1 10 . The simulations are compared on the basis of the error e , defined as the Euclidian (L 2 ) norm of the di↵erence between the model solution, h , and the projection of the exact solution onto the space of the FE approximation after one complete revolution (at t = 1)

e = k h (t = 1) ⇧ S h ( )(t = 1)k L 2 (⌦) , (3.6)
with ⇧ S h representing the Lagrange interpolation operator. Under some reasonable assumption on the mesh family and the regularity of the solution, we expect the error to converge optimally with t and x as

e ( x, t) = C 1 t p + C 2 x k+1 , (3.7)
where p is the order of the temporal discretization scheme and k the order of the FE polynomial approximation, and C 1 and C 2 are some constants [START_REF] Richter | An optimal-order error estimate for the discontinuous galerkin method[END_REF].

P 0 approximations

The di↵usivity of the numerical scheme is first evaluated for finite element approximations of order k = 0. The function h is then constant by element, as are the fields of , d and C in the implementation of the Maxwell-EB framework described in the previous chapter (see section 2.3.1). In this case the DG method coincides with a finite volume approximation with upwinding. The first sum in the bilinear form b h (), equation (3.31), disappears as the gradient of h is zero inside elements. The only remaining contribution to the gradient of h is that of the jumps in h across the set of internal sides of the mesh elements, L

h . Figure 3.2 gives a taste of the amount of di↵usion to expect in the model for P 0 FE approximations. The initial fields of 1 , 2 and 3 are compared to their values after one quarter and a complete revolution, in the case of a relatively high temporal and spatial resolution ( t = 1 10000 , x = 1 40 ). For the most regular function,

1 , the di↵erence in amplitude of the initial and final approximations of 1 is of 25%-35% (note that the solution appears very noisy due to the fact that we plot the value of as a function of ✓ for all elements in ⌦ h , i.e., not along a constant radius). In the case of the step function 2 , initially discontinuous at ✓ = ± ⇡ 4 , di↵usion is such that the peak of the function decreases while its width increases. After one complete revolution, the function looks like a smooth sinusoid. In the case of the field of noise, as we are not interested in the specific spatial structure of 3 , we instead evaluate the numerical di↵usion based on the evolution of its distribution. Figure 3.2c shows that after 1/4 of a revolution, 3 is almost completely homogenized. In the context of the Maxwell-EB model, this implies an essentially complete suppression of the spatial noise in the damage criteria.

Figure 3.3 shows the error, e , as a function of t, for simulations run at di↵erent spatial resolutions x (coloured lines) and using a time discretization scheme of order 1 (solid lines), 2 (dashed lines) and 3 (dotted lines). In the case of 1 (see figure 3.3a) and 2 (figure 3.3b), the error becomes independent of t for x t between 10 2 and 10 3 and that, for all orders of the time discretization scheme. In the case of sea ice on geophysical scales, this ratio is of the same order of magnitude as the speed of propagation of elastic waves (c ⇡ 500 ms 1 , see chapter 2) which sets the relative values of the spatial and temporal resolutions in the Maxwell-EB framework (i.e., when using P 0 approximations for the damage criterion). This suggests that, in a Maxwell-EB sea ice model using P 0 approximations for the fields of C, d and (1) the value of the time step would be such that the error depends essentially on the spatial resolution and (2) the use of a higher order temporal discretization scheme would therefore not be e cient in reducing the numerical di↵usion of gradients.

In the case of 1 (figure 3.4a), the rate of convergence of the error, minimized with respected to t, e ( t ! 0), with decreasing spatial resolution is O( x 1 2 ), consistent with the DG method error estimates of [START_REF] Johnson | An analysis of the discontinuous galerkin method for a scalar hyperbolic equation[END_REF] for polynomial approximations of order 0. This is smaller however than the O( x k+1 ) convergence rate demonstrated by [START_REF] Richter | An optimal-order error estimate for the discontinuous galerkin method[END_REF] assuming some regularity on the solution and often observed for common mesh families. This could be due to the fact that here, 1 2 H 2 as the second derivative is locally discontinuous at ✓ = ± ⇡ 4 , and hence is not a su ciently smooth function. For the step function 2 , the error is larger than for 1 but does converge as O( x 1 4 ) (see figure 3.4b). The error calculated for the field of noise, 3 , does not converge neither in space nor in time and for any order of the time discretization scheme (see figures 3.3c and 3.4c). For P 0 approximations of d and , we can therefore expect the numerical 63 h using FE approximations of order k = 0. (Lower panels) Corresponding P0 approximation of the fields of 1 and phi2 at t = 0 (black curve), t = 1/4 (blue curve) and t = 1 (red curve) and distribution of phi3 at t = 0, t = 1/4 and t = 1. di↵usion of gradients in the Maxwell-EB model to be important and to decrease slowly with increasing spatial resolution, that is, at a rate of at most O( x 1 4 ).

Higher order FE approximations

We repeat the error calculations for higher order FE approximations and compare the di↵usivity of the numerical scheme to the P 0 case. Similar to figure 3.2, figures 3.5 and 3.6 (upper panels) compare the initial fields of 1 , 2 and the initial distribution of 3 to their values after one quarter and one complete revolution in the case of P 1 and P 2 discontinuous approximations respectively. The order of the time discretization scheme is increased to p = 2 in the P 1 and to p = 3 in the P 2 case and t = 1 10000 and x = 1 40 . For the most regular function 1 2 H 2 , increasing the order of the FE approximations by 1 degree is su cient to practically eliminate numerical di↵usion (see figure 3.5a). An optimal convergence rate of the error of O( x 2 ) is recovered. For the step function 2 , the numerical scheme preserves much sharper gradients around the initial location of the discontinuities and presents a higher rate of convergence of the error, that is, of O( x 1 3 ) for P 1 approximations and of O( x 1 2 ) for P 2 approximations. Increasing the order of the polynomial approximations does decrease the rate at which the field of noise, 3 , homogenizes, but does not increase the rate of convergence of the error (see figure 3.5c). For P 2 approximations, the error somewhat increases at the highest resolution. One important point is that in the case of an implicit time discretization scheme, the DG method results in non-monotonic dispersion and generates values of outside the interval over which it is initially defined. This e↵ect is more important around local discontinuities and as the advected gradients are strong (compare for instance figures 3.5a and 3.6a with figures 3.5b and 3.6b around ✓ = ± ⇡ 4 ). Yet, in the Maxwell-EB framework, some of the advected fields must evolve within a certain range of values, for instance the level of damage d, which is defined on the interval [0, 1]. As will be discussed in chapter 4 in the context of sea ice modelling, it is also the case for the ice thickness (strictly positive) and the ice concentration, defined on the interval [0%, 100%]. While these variables could be "manually" reset within their prescribed interval in the model, the generation of spurious oscillations and unphysical values could also be avoided by replacing the implicit time discretization scheme used here by a Strong Stability-Preserving Runge-Kutta scheme [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF] for which the dispersion of the DG method is known to be strictly monotonic (e.g., [START_REF] Hanert | Advection schemes for unstructured grid ocean modelling[END_REF]. When using a time discretization scheme of order 1, the implicit method presented above can be replaced by a Euler explicit scheme, in which case the advection terms in the bilinear form b h () simply move to the linear form o h ().

For su ciently small time steps, i.e., on the order of that set for sea ice (i.e, t ⇡ x c with c ⇡ 10 2 10 3 ), we find that the magnitude of the error associated with the numerical di↵usion is comparable between the Euler implicit and Euler explicit schemes (results not shown).

These simple tests suggest that the DG method could be e cient in handling the advection of the strong gradients generated by the Maxwell-EB model, provided that the order of the FE element approximations for the advected fields can be of order k > 0. However, using polynomial approximations of order k > 0 pose a consequential problem when evaluating the distance of the internal stress to the local damage criterion, which in the current implementation of the framework varies at the scale of the mesh element.

To allow increasing the order of the FE approximations, disorder needs to be introduced in the model by a more regular function, varying on scales larger than x. This amounts to assuming a correlation length

t = 0 t = 1/4 t = 1 -2 0 2 0 0.2 0.4 0.6 0.8 1 t = 0 t = 1/4 t = 1 (a) (b) (c) 1/3 2 Figure 3
.5: (Upper panels) P1 discontinuous FE approximations of the fields of (a) 1 and (b) 2 at t = 0 (black curve), t = 1/4 (blue curve) and t = 1 (red curve) and (c) of the distribution of 3 at t = 0, t = 1/4 and t = 1. (Lower panels) corresponding L 2 norm of the error for t = 1 10000 (i.e., the minimum error with respect to time) as a function of the spatial resolution x for the implicit scheme of order 2. Note the di↵erent extent of the vertical axis between panel (a) and panels (b) and (c).

⇠ > x for heterogeneities within the Maxwell-EB material. Doing so, the characteristic time for propagation of damage would then be set as ⇠ c . Due to elastic interactions, deformation would tend to localize at a scale < ⇠, but > x. In the case of sea ice modelling on regional to global scales, this assumption is however not physical, as the true correlation length of associated with strong fluctuations within the ice cover is likely much smaller than that of the model grid cell [START_REF] Schulson | Compressive shear faults within arctic sea ice: Fractures on scales large and small[END_REF].

When choosing the order of the FE approximations for the Maxwell-EB model, there is therefore a tradeo↵ between the mitigation of di↵usion and the resolution of the progressive damage mechanism. In the following implementation of the rheology, we persist in introducing disorder at the scale of the model element and hence use P 0 approximations for the fields of C, d and (i.e., all advected fields in small Cauchy number flow problems). Given the current formulation of the damage mechanism (equation 2.27), we believe this is the most rigorous choice in terms of respecting the underlying physics. This also allows for the highest resolution of the damaging process. In the simulations presented here, levels of numerical di↵usion comparable to a finite volume scheme with upwinding are therefore to be expected.

At this point, the problem of the homogenization of the field of cohesion in the case of quenched disorder (i.e., when noise is introduced in the field of C once, at the beginning of the simulations) remains unsolved.

In order for the mechanical framework to reproduce an heterogeneous and intermittent mechanical behaviour when deformations within the material are large and advection processes are important, or when simulation times become long, spatial disorder needs to be regularly "re-injected" in the model. One way to handle this is to redraw the local value of cohesion from its initial distribution every time a given element becomes overcritical. This annealed-type disorder can be thought as to represent the heterogeneity in the strength of an elasto-brittle material during the damage increase [START_REF] Amitrano | From di↵use to localised damage through elastic interaction[END_REF] and does not a↵ect the macroscopic (Upper panels) P2 discontinuous FE approximations of the fields of (a) 1 and (b) 2 at t = 0 (black curve), t = 1/4 (blue curve) and t = 1 (red curve) and (c) of the distribution of 3 at t = 0, t = 1/4 and t = 1. (Lower panels) corresponding L 2 norm of the error for t = 1 10000 (i.e., the minimum error with respect to time) as a function of the spatial resolution x for the implicit scheme of order 3. Note the di↵erent extent of the vertical axis between the three panels.

behaviour nor the scaling laws characterizing the deformation of the material [START_REF] Herrmann | Statistical Models for the fracture of disordered meria[END_REF].

Rotation and deformation

In this section we present the numerical scheme for solving the Maxwell-EB constitutive equation in which the convection, rotation and deformation terms are all included and we evaluate the capacity of this scheme to handle large Weissenberg numbers. To do so, we solve a simplified version of the Maxwell-EB system of equations in which the mechanical parameters E, ⌘ and are kept constant, that is, in which the level of damage does not evolve. With d = 1 everywhere and at all times, the (adimensional) Maxwell-EB system of equations reduces to:

Ca  @ ũ @ t + (ũ • r)ũ = Fext + r • ˜ (3.8) We  @ ˜ @ t + (ũ • r)˜ + a ( rũ, ˜ ) + ˜ = We ⇣ ⇤ div(ũ) • I + 2G D(ũ) ⌘ (3.9)
with Ca = U 2 E and We = ⌘ E U L , the Cauchy and Weissenberg numbers introduced in section 2.2.4 and a ( rũ, ˜ ) = ˜ W (ũ) W (ũ)˜ a (˜ D(ũ) + D(ũ)˜ ), with a 2 { 1, 0, 1}. In our implementation of the Maxwell-EB framework, we use the upper-convected objective derivative for and hence set a = 1.

The geometry and boundary conditions of the problem solved here are that of a typical annular Couette flow experiment (see figure 3.7). In such an experiment, a material is being sheared in between two concentric circular boundaries, the angular velocity of which is prescribed. In the present case, the Maxwell-EB material is 2-dimensional, homogeneous and undamaged, the external boundary is fixed (u = 0), the internal boundary has a velocity u = U e ✓ and F ext = 0. The axial symmetry in this particular geometry is such that the solutions of u and , when cast in polar coordinates, are expected to vary only in the radial direction. With Ca small, as estimated in the case of an ice pack drifting on large scales (see section 2.2.4), and the inertial term neglected in the momentum equation, the only adimensional parameters characterizing the Couette flow problem are the Weissenberg number, We = U L = T , with L = R ext R int the distance between the internal and external boundaries, and the adimensional model time step, ˜ t = t T . In the following, the convergence towards a stationary solution is investigated for a semi-implicit and implicit numerical scheme and for di↵erent values of We and ˜ t.

Numerical scheme

Equations (3.8) and (3.9) with Ca = 0 and Fext = 0 are solved over an unstructured Eulerian grid with triangular elements of constant mean size (see figure 3.1b). We consider the domain ⌦ ⇢ R 2 represented in figure 3.7a and delimited by the boundary @⌦ = ext [ int , with Rint Rext = 0.4, and the time t 2 [0, +1[. Forcing is applied on the system by prescribing the angular velocity on the boundaries. The Dirichlet conditions are ũ(t) = 0 on ext and ũ(t) = 1e ✓ on and int . There are no flux in or out of the domain. The problem reads:

(P ) : Find ũ and ˜ defined in ⌦⇥]0, +1[, such that

0 = r • ˜ We  @ ˜ @ t + (ũ • r)˜ + a ( rũ, ˜ ) + ˜ = We ⇣ ⇤ div(ũ) • I + 2G D(ũ) ⌘ with initial conditions ũ( t = 0) = 0 in ⌦, ˜ ( t = 0) = 0 in ⌦,
and boundary conditions

ũ( t) = 0 on ext ⇥]0, +1[ ũ( t) = 1e ✓ on int ⇥]0, +1[
In the following, we drop the superscript 0 ˜0 for the adimensional variables and operators.

Non-linearity in this problem arises through the advection, rotation and deformation of the stress tensor

. Discretizing time such that t n = n t, with t > 0 and n = 0, 1, 2, ..., the constitutive equation can be linearized by using either (1) a semi-implicit scheme in which the advection, rotation and deformation terms are estimated using the field of velocity and internal stress at the previous (n th ) model time step or (2) an implicit scheme together with a fixed-point iteration for the non-linear terms. We present the details of the semi-implicit version of the numerical scheme. The passage from the semi-implicit to the implicit formulation is then straightforward. With a semi-implicit scheme, the time-discretized system of equations reads

r • n+1 = 0 (3.10) We  n+1 n t + (u n • r) n + a (ru n , n ) + n+1 = We K : D(u n+1 ) (3.11)
With the constitutive equation linearized, problem P is solved in two steps as in the small-deformation, uniaxial compression case (see section 2.3.1). First, the momentum equation is solved for u n+1 by substituting for the following expression of n+1 into equation (3.10) :

n+1 = KWe  1 t n [(r • u n ) n + a (ru n , n )] + K : D(u n+1 ) (3.12)
where K = ⇥

We t + 1 ⇤ 1 . Second, the internal stress, n+1 , is computed explicitly by substituting back for u n+1 into (3.12). In the following, we present the continuous and discrete variational formulations of the problem.

Introducing the continuous functional spaces

V 0 = {v 2 H 1 (⌦) 2 ; v = 0 on } (3.13) and V = {v 2 H 1 (⌦) 2 ; v = 1e ✓ on int and v = 0 on ext } (3.14)
the weak form of the momentum equation writes:

Z ⌦ r • n+1 v dx = Z ⌦ n+1 : D(v) dx + Z @⌦ n+1 • n v ds = 0, 8v 2 V 0 ,
with the boundary integral vanishing because the test function v vanishes on the boundary. Substituting for n+1 given by (3.12) leads to a linear problem expressed in terms of u n+1 only:

(P.1) : Find u n+1 2 V such that a(u n+1 , v) = l(v), 8v 2 V 0
where the bilinear form a() and linear form l() are defined for all u, v 2 H 1 (⌦) 2 as

a(u, v) = Z ⌦ KWe (K : D(u)) : D(v) dx l(v) = Z ⌦ KWe  (u n • r) n + a (ru n , n ) 1 t n : D(v) dx
and where u n and n are known from the n th model time step. In terms of the dimensionless Lamé coe cients ⇤ and G, the bilinear form a() writes

a(u, v) = Z ⌦ KWe [⇤div(u) • div(v) + 2GD(u) : D(v)] dx. (3.15)
The problem is discretized in space by introducing the two-dimensional mesh T h and the functional space of P 1 elements

X h = {v h 2 H 1 (⌦) \ C 0 ( ⌦) 2 ; v h|K 2 P 2 1 , 8K 2 T h } (3.16)
for the velocity field. Let us define another test function ⌧ h 2 T h with T h , the functional space

T h = {⌧ h 2 (L 2 (⌦)) 2⇥2 ; ⌧ h = ⌧ T h and ⌧ h|K 2 (P 0 ) 2⇥2 , 8K 2 T h } (3.17)
for the stress tensor. Then the discontinuous Galerkin approximation of the convective term for the stress tensor leads to the following approximation of the linear form l():

l h (v) = Z ⌦ h  ⌧ n h + a (ru n h , n h ) 1 t n h : D(v h ) dx (3.18)
where n h 2 T h is the piecewise approximation of n , u n h 2 X h is the approximation of u n and ⌧ n h 2 T h is an approximation of (u n • r) n defined for all ⌧ h 2 T h by

Z ⌦ ⌧ n h : ⌧ h dx = X K2T h Z K (u n h • r) n h : ⌧ h dx + X K2T h X S2@K Z S [[ n h ]] : ✓ |u n h • n| 2 [[⌧ h ]] (u n h • n){{⌧ h }} ◆ ds (3.19)
with n, the normal on an oriented element side. In the case of a P 0 FE approximation for the stress tensor, (3.19) simplifies as the sum of the gradient of h is zero inside elements and the contribution from jumps in h across element sides is non-zero only over the set of the internal sides of the mesh T h , L (i) h . The discrete version of expression (3.12) for the stress tensor reads:

n+1 h = KWe  1 t n h [⌧ n h + a (ru n h , n h )] + K : D(u n+1 h ) . (3.20)
With these definitions, the finite dimensional Couette problem for an homogeneous, or undamaged, material writes:

u n h = 0 in ⌦ h n h = 0 in ⌦ h
For n 0

• (P.1) h : With u n h and n h known, find u n+1 h 2 V h such that a(u n+1 h , v h ) = l n+1 h (v h ), 8v h 2 V 0,h where V 0,h = {v h 2 X h , v h = 0 on } V h = {v h 2 X h , v h = 1e ✓ on int and v h = 0 on ext }
and with the bilinear form a() and linear form l h () defined by (3.15) and (3.18) respectively.

• (P.2) h : With u n+1 h known, compute n+1 h 2 T h explicitly from (3.20).
This scheme can be enhanced into an implicit version by using a fixed-point algorithm to update the values of ⌧ n h and of a (ru h , h ). In this case, the momentum and constitutive equations are solved iteratively until the residual of the constitutive equation drops below a chosen tolerance, tol. Using the superscript k = 0, 1, 2, ... for the sub-iterations, the discretized fully implicit scheme writes:

(P ) h : Initialization (n = 0)

u n h = 0 in ⌦ h n h = 0 in ⌦ h
For n 0

• For k = 0, set ⇣ n+1,0 h , u n+1,0 h ⌘ = ( n h , u n h ),
• For k 0

-(P.1) h : With u n h and n h known, find u n+1 h 2 V h such that a(u h , v h ) = l (n+1,k+1) h (v h ), 8v h 2 V 0,h ,
with with the bilinear form a() given by (3.15) and the linear form l h () given by

l n+1,k+1 h (v h ) = Z ⌦ h KWe  ⌧ n+1,k h + a (ru n+1,k h , n+1,k h ) 1 t n h : D(v h ) dx
and with

Z ⌦ ⌧ n+1,k h : ⌧ h dx = X K2T h Z K ⇣ u n+1,k h • r ⌘ n+1,k h : ⌧ h dx + X S2L (i) h Z S [[ n+1,k h ]] : |u n+1,k h • n| 2 [[⌧ h ]] (u n+1,k h • n){{⌧ h }} ! -(P.2) h : With u n+1,k+1 h known, compute n+1,k+1 h 2 T h explicitly from n+1,k+1 h = KWe  1 t n h ⇣ ⌧ n+1,k h + a (ru n+1,k h , n+1,k h ) ⌘ + K : D(u n+1,k+1 h ) -Compute res h = We " n+1,k+1 h n h t + ⇣ ⌧ n+1,k+1 h + a (r h u n+1,k+1 h , n+1,k+1 h ⌘ # + n+1,k+1 h We K : D(u n+1,k+1 h ) , If res h < tol STOP , • Set n+1 h , u n+1 h = ⇣ n+1,k+1 h , u n+1,k+1 h ⌘ .
For the simulations presented in this section that use this implicit scheme, the tolerance was set to the machine precision.

Stationary solution

We first establish the convergence of both the semi-implicit and implicit numerical schemes towards a stationary solution for We 2 {0.001, 0.01, 0.1, 1} and t 2 {1, 1 10 , 1 100 , 1 1000 , 1 10000 , 1 100000 }. For We = 1, neither the semiimplicit nor the implicit scheme converges, and that, for any value of t. For We = 0.1, a stationary solution is obtained for t  1 100 using both schemes, with a practically equivalent rate of convergence (see figure 3.8). For lower values of We, both schemes converge for all t. Comparing simulations using mesh grids with x = 1 10 , 1 20 , 1 40 , 1 80 (not shown) indicates that the convergence towards the stationary solution does not depend on the spatial resolution [START_REF] Keunings | On the high Weissenberg number problem[END_REF]. int, as a function of the simulation time in the explicit and implicit scheme (black curves) and di↵erence between the solutions for both schemes (red curve) for We 0 = 0.1 and t = 1 100 . The di↵erence is of maximum 0.3% of the calculated torque between the explicit and implicit schemes.

We therefore conclude that (1) the semi-implicit and implicit schemes presented above are equivalent in terms of handling the convection, rotation and deformation terms for high We values and (2) the highest value of the Weissenberg number that can be used to ensure the stability and convergence of our numerical scheme in this geometry is of about 0.1.

In light of these tests, we therefore use an explicit formulation of the discretized advection, rotation and deformation terms in the following implementation of the rheology. Compared to the implicit formulation, the explicit approach reduces computational costs considerably since the piecewise approximations ⌧ h and a (ru h , h ) defined above are then calculated once per model time step, instead of an indeterminate number of times within a fixed-point iteration. In order to ensure the convergence of our numerical scheme in the flow geometries explored here, we set the undamaged (i.e., maximum) relaxation time 0 in all model simulations started from an undamaged state to at most 0 = 0.1 ⇥ U L , with U L the rate of deformation characterizing the simulated process on the large scale. Considering the small time-stepping imposed by our discrete-like formulation of the progressive damage mechanism, the requirement that the (adimensional) model time step be t < 1 100 for the convergence of the numerical scheme with We 0 = 0.1 is naturally met in the case of regional to global sea ice models. The dependance of We on the deformation rate U L implies that the maximum allowed value of the relaxation time for the internal stress 0 decreases as the simulated strain-rates are large, hence that the capacity of the model to represent the limit of strictly elastic deformations within an undamaged material decreases in the case of large-deformation simulations. When and where the Maxwell-EB material becomes damaged, deformation rates are expected to increase as the material weakens. However, as the local value of the relaxation time, = 0 d ↵ 1 drops like d ↵ 1 (with ↵ > 1), the damaging process acts so that to keep the local value of We small.

In the case of the ice cover, we can deduce the maximal allowed value of the undamaged relaxation time 0 for the current implementation of the Maxwell-EB model based on observational deformation rate estimates. [START_REF] Marsan | Scale dependance and localization of the deformation of Arctic sea ice[END_REF] and [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], among others, have estimated deformation rates between 10 4 and 10 0 day 1 from the 10 km RGPS motion products. Assuming for instance that the lower bound of this interval corresponds to small elastic deformations over low-damaged portions of the ice pack, this implies a maximum value of 0 in simulations with spatial resolution of ⇠ 10 km and intended for the representation the deformation of the ice pack on daily time scales, of about 10 3 days (⇠ 3 years). As this is much larger than the length of the polar winter (about 180 days), during which large portions of the ice pack are expected to remain relatively undamaged and hence to retain the memory of elastic deformations, unphysical viscous dissipation in the model is not expected to be significant. Conversely, where deformation rates are large, i.e. on the order of 10 0 day 1 , the relaxation time should drop with d to values below ⇠ 0.1 day to ensure numerical stability. Appropriate scaling laws should be used when interpolating these estimates onto di↵erent time and space scales [START_REF] Rampal | Scaling properties of sea ice deformation from buoy dispersion analysis[END_REF].

We can now further comment on the form the solution obtained in these Couette simulations. Analytical solutions have been derived for the laminar Couette flow problem in the case of the (upper and lower-convected) Oldroyd and Bingham models (e.g., [START_REF] Saramito | Complex fluids. Modelling and algorithms[END_REF][START_REF] Cheddadi | Numerical modelling of foam couette flows[END_REF][START_REF] Cheddadi | Modelisation numerique d'ecoulements de mousse[END_REF] for viscoelastic fluids. The Maxwell model is a special, simpler case of these models, without a plastic threshold for the deformation or a solvent viscosity term. With the divergence of the velocity being zero in incompressible fluids, the radial velocity in this problem is, by construction, also zero. For the upper and lower-convected derivatives (a ± 1), the equations are linear and the system can be solved explicitly (e.g., [START_REF] Cheddadi | Modelisation numerique d'ecoulements de mousse[END_REF].

With the domain geometry and boundary conditions as shown on figure 3.7 and a = 1, the analytical solution is given by:

u ✓ = UR int R 2 ext R 2 ext R 2 int  1 r r R 2 ext , u r = 0, rr = 0, r✓ = 2WeG UR int R 2 ext R 2 ext R 2 int 1 r 2 , ✓✓ = 2G  2We UR int R 2 ext R 2 ext R 2 int 1 r 2 2 .
As here we consider a compressible material and do not impose mass conservation, the solution to (3.8) and

(3.9) with Ca = 0 and Fext = 0 cannot be derived explicitly in the case of a Maxwell-EB material. Figure 3.9 shows the P 0 approximation of the (adimensional) stationary solution for the velocity and stress tensor components in polar coordinates, as a function of the radial distance to the centre of the domain, r, for a mesh resolution of x = 1 40 . It is important to note that, as in section 3.2, the solutions appear noisy due to the fact that we plot the value of the fields over all elements (or nodes) of the unstructured mesh, instead of along a radial cut though the domain. We otherwise checked that the solutions converge with increasing spatial resolution. Comparing the model and analytical solution for an incompressible Maxwell fluid (black curves), we indeed find the largest di↵erence for the normal component rr and for We = 0.1. With a Poisson's ratio of ⌫ = 0.3 and We = 0.001, 0.01, 0.1, we find that the model solution for the shearing stress r✓ and the ✓✓ component is close to the analytical solution for an incompressible Maxwell fluid (black curves). This suggests that r✓ and ✓✓ in the Maxwell-EB model vary approximately as We and We 2 respectively.

If neglecting the rotation and deformation terms, i.e., the a term, in the Maxwell-EB constitutive equation, both rr and ✓✓ would then be 0. According to the solution with We ⇡ 0.1, with the Weissenberg number taken as large as possible to represent the limit of a purely elastic solid, we can expect these normal stresses to be important in an undamaged material where velocity gradients are large. However, because the largest velocity gradients in the Maxwell-EB material arise were the material is highly damaged and the e↵ective Weissenberg number, We, varies as d ↵ 1 with ↵ > 1, the magnitude of these normal stresses is expected to decrease as d 2(↵ 1) and hence become negligible.

The full, large-deformation Maxwell-EB model

In this section, we investigate the behaviour of the full Maxwell-EB model in the context of large-deformation simulations. With the progressive damage mechanism included, there is obviously no analytical solution available to us for the validation the Maxwell-EB framework. Hence here we compare the mechanical behaviour of the model to the results of a laboratory Couette flow experiment that was performed on a thin sheet of fresh-water ice with the objective of investigating the competing e↵ects of fracturing and healing within a fault formed naturally by shear failure.

Laboratory Couette experimental setup

The experimental setup is represented in figure 3.10. A circular tank of radius R ext = 0.5 m was filled with fresh water. A circular disk of radius R int = 0.2 m was partially immersed at the middle of the tank. The apparatus was placed in a cold room at subzero temperatures so that to form a layer of ice about 4 mm (±1 mm) thick. The temperature of the water was regulated by placing heating pads at the bottom of the tank.

The thin ice plate was deformed by prescribing the constant angular velocity of the inner rotating disk. The associated torque, (t), was measured using a torquemeter with precision of 1 Nm 1 . The circular Couette geometry implies a tangential shear stress in an homogeneous, elastic material that varies in 1 r 2 , with r the radial distance from the centre of the apparatus. Hence the shear stress in the initially undamaged ice plate was expected to be maximum at R int . In order to avoid for the ice to fail right along the rotating disk, it was sanded to enable the ice to better "stick" and grow a narrow meniscus. This e↵ectively allowed for the main shearing fault to form a few centimetres away from the disk.

Numerical Couette experiment setup

In the numerical simulations designed to represent this annular Couette flow experiment, the thin ice plate formed at the surface of the water tank is assimilated to a 2-dimensional continuum Maxwell-EB material. The domain is shown in figure 3.7a and has dimensions R ext and R int as in the laboratory experiment. Loading is applied on the material by prescribing the angular velocity on the inner boundary, while u = 0 on the outer boundary. No other external forcing is applied on the ice (F ext = 0) and the drag from the underlying water, probably negligible considering the applied rotation speeds, is not represented. Thermodynamic e↵ects and dynamic variations of the ice volume are not accounted for either and mass conservation is not enforced, as in the small-deformation experiments. Simulations are run until the cumulated displacement along the inner boundary is of at least a few degrees. Hence advection, rotation and deformation terms are all potentially important and the complete constitutive equation (2.6) is solved. We checked that the introduction of this step-like function does not generate numerical instabilities and does not significantly impact the simulations results presented here.

Based on the analytical solution of the annular Couette flow problem of an homogeneous Maxwell material and as in the laboratory experiments, the shear stress in the initially undamaged Maxwell-EB material is expected to be maximum at R int , hence deformation and damaging are expected to concentrate at the inner domain boundary. In order to resolve the propagation and width of shearing faults, the critical strength of the simulated ice plate is increased near the rotating inner boundary by adding an increment to the field of cohesion, C, defined as 4C min (1 [r R int ]) 40 with C min the lower bound of the uniform distribution of C.

The form of this completely ad hoc function ensures that this increment is restricted to a few grids cells near int . As the ice meniscus in the laboratory experiments, it allows for the simulated macroscopic fault to form a few grid cells away from the inner domain boundary.

In these simulations, the Maxwell-EB model is made adimensional with respect to the undamaged elastic modulus of the (fresh) ice plate E 0 , the lower bound of the distribution of C, C min , the distance L = R ext R int , between the inner and outer domain boundaries, the prescribed velocity on the inner boundary, U = !R int , and the corresponding shearing deformation rate U L = 1 T . The velocity on the inner boundary is small enough so that Ca 0 ⌧ 1 and the inertial term is neglected in the momentum equation.

r • ˜ = 0 We 0 d ↵ 1  @˜ @ t + (ũ • r)˜ + a ( rũ, ˜ ) + d1 de˜ = We 0 d 0 ↵ K : "(ũ) @d @ t + (ũ • r)d = min " 1, ⌃ t C ˜ 2 , ⌃ c C ˜ 1 q ˜ 2 # 1 ! 1 T d d + 1 T h , 0 < d  1 @ C @ t + (ũ • r) C = 0 (3.21)
Table 3.1 lists the variables, dimensions and parameters relevant to the laboratory and numerical experiments.

As deformation and damaging are expected to concentrate along a circular fault near the inner boundary, leaving a large portion of the domain motionless, we do not set the mean mesh element size constant here but instead let it vary as x(r) = Rext N r 2 , with Rext N the spatial resolution at the outer boundary (see figure 3.7b) and N = 10. This choice was made to allow increasing the spatial resolution where required while maintaining computational costs low. We verified the convergence of the algorithm towards a stationary solution over this mesh grid in the no-damaging case and for the values of (adimensional) model time step and We 0 used here. Section 3.4.3 describes the time and space discretizations and gives the details of the discontinuous Galerkin methods-based numerical scheme employed to solve the system of equations (3.21).

Numerical scheme : the large-deformation Maxwell-EB model

The computational domain ⌦ ⇢ R 2 with boundary @⌦ = ext [ int is represented in figure 3.7a. Dirichlet conditions apply on both ext and int . There is strictly no flux in or out of the domain. Simulations are started from rest, from a homogeneous undamaged state (d = 1) and run for a time t 2 [0, +1[. With these definitions, the annular Couette flow problem for a damageable material reads:

(P ) : Find ũ, ˜ , d and C, defined in ⌦⇥]0, +1[, such that

r • ˜ = 0 in ⌦⇥]0, +1[, We 0 d ↵ 1  @˜ @ t + (ũ • r)˜ + a ( rũ, ˜ ) + d1 de˜ = We 0 d 0 ↵ K : D(ũ) in ⌦⇥]0, +1[, @d @ t + (ũ • r)d = min " 1, ⌃ t C ˜ 2 , ⌃ c C ˜ 1 q ˜ 2 # 1 ! 1 T d d + 1 T h , 0 < d  1 in ⌦⇥]0, 1[, @ C @ t + (ũ • r) C = 0, in ⌦⇥]0, +1[, where
⌃ c = C min E 0 2 ⇥ (µ 2 + 1) 1/2 µ ⇤ , in ⌦⇥]0, +1[, ⌃ t = 2 C min E 0 h (µ 2 + 1) 1/2 + µ i , in ⌦⇥]0, +1[,
and with initial conditions

ũ( t = 0) = 0 in ⌦, ˜ ( t = 0) = 0 in ⌦, d( t = 0) = 1 in ⌦, C( t = 0) = C0 in ⌦,
and boundary conditions

ũ( t) = 0 on ext ⇥]0, +1[, ũ( t) = 1e ✓ on int ⇥]0, +1[.
Here C0 represents the initial field of cohesion, in which spatial noise is introduced at the scale of the element.

To facilitate the reading, the superscript 0 ˜0 for adimensional variables is dropped in the remaining of this section.

Time discretization

Non-linearity in the full, large-deformation Maxwell-EB model arises through advection, through the rotation and deformation of the stress tensor in the constitutive equation, through the coupling of the cohesion trans-port and damage equations and, as in the small-deformation model, through the coupling of the constitutive and damage equations. The system of equations is discretized in time as follow.

• A Euler implicit scheme of order 1 is used for the momentum equation.

• A semi-implicit scheme is used for the constitutive equation, in which the advection, rotation and deformation terms are estimated using the field of velocity and internal stress computed at the previous (n th ) model time step.

• A semi-implicit scheme is used for the damage equation. The advection term is explicit and uses the velocity at the n th time step. The distance to the damage criterion (d crit ) is evaluated using the field of internal stress at the current time step (n + 1) and the local damage criterion (i.e., C) at the previous time step.

• An explicit scheme with the value of the velocity field at the n + 1 model iteration is used for the cohesion transport equation.

Let us discretize t such that t n = n t, with t > 0 and n = 0, 1, 2, .... The system of equations then writes:

r • n+1 = 0 We 0 d n+1 ↵ 1  n+1 n t + (u n • r) n + a (ru n , n ) + d1 d n+1 e n+1 = We 0 d 0n+1 ↵ K : D(u n+1 ) d n+1 d n t + (u n • r)d n = ✓ min  1, ⌃ t C n n+1 2 , ⌃ c C n n+1 1 q n+1 2 1 ◆ 1 T d d n + 1 T h , 0 < d n+1  1 C n+1 C n t + (u n+1 • r)C n+1 = 0
The problem is solved using a fixed-point algorithm, which divides it in three linear subproblems. Initializing d with its previous time step value, the constitutive and momentum equations are first solved for n+1 and u n+1

(subproblem P1) as in the small-deformation model. Then the damage evolution equation is solved using the fields of n+1 and u n (subproblem P2). The value of d is updated and these computations are iterated until the residual of the constitutive equation drops below a chosen tolerance. The cohesion transport equation is then solved using u n+1 and the value of the local damage criterion (i.e., ⌃ c C and ⌃ t C) is updated (subproblem P3). Using the superscript k = 0, 1, 2, ... for the sub-iterations, the fixed point algorithm reads:

For k = 0, let n+1,0 , u n+1,0 , d n+1,0 = ( n , u n , d n ),
For k 0

• (P1) Find n+1,k+1 and u n+1,k+1 such that r • n+1,k+1 = 0 (3.22) We 0 d n+1,k ↵ 1  n+1,k+1 n t + (u n • r) n + a (ru n , n ) + d1 d n+1,k e n+1,k+1 = We 0 d 0n+1,k ↵ K : D(u n+1,k+1 ) (3.23)
and with u n+1,k+1 = 0 on ext , (3.24)

u n+1,k+1 = 1e ✓ on int .
(3.25)

• (P2) Find d n+1,k+1 , such that 0 < d n+1,k+1  1 and

d n+1,k+1 d n t + (u n • r)d n+1 = ⇣ min h 1, ⌃tC n 2 n+1,k+1 , ⌃cC n n+1,k+1 1 q n+1,k+1 2 i 1 ⌘ 1 T d d n + 1 T h . (3.26) • Stopping criterion : compute res = We 0 (d n+1,k+1 ) ↵ 1  n+1,k+1 n t + (u n • r) n + a (ru n , n ) + d1 d n+1,k+1 e n+1,k+1
We 0 d 0n+1,k+1 ↵ K : D(u n+1,k+1 )

If res < tol STOP , and set n+1 , u n+1 , d n+1 = n+1,k+1 , u n+1,k+1 , d n+1,k+1 .

As in the small-deformations and no-damaging case (see sections 2.3.1 and 3.3), subproblem P1 is solved in two steps. The momentum equation (3.22) is first solved for u n+1,k+1 by substituting for the following expression of n+1,k+1 :

n+1,k+1 = K " We 0 d n+1,k ↵ 1 t n We 0 d n+1,k ↵ 1 [(u n • r) n + a (ru n , n )] + We 0 d 0n+1,k ↵ K : D(u n+1,k+1 ) i (3.27) with K = h We 0 d n+1,k ↵ 1 1 t + d1 d n+1,k e i 1
. The internal stress n+1,k+1 is then calculated by substituting u n+1,k+1 back into (3.27). All subproblems make use of a variational formulation and of discontinuous Galerkin finite element methods.

Variational formulation and discontinuous Galerkin FE approximation

Subproblem P1

Introducing the continuous functional spaces V 0 and V defined by (3.13) and (3.14) respectively, and substituting for (3.27) in the weak form of equation (3.22), the momentum balance subproblem P1.1 reads:

(P1.1) : Find u n+1,k+1 2 V such that a(u n+1,k+1 , v) = l(v), 8v 2 V 0 with a() and l() defined for all u, v 2 H 1 (⌦) 2 as a(u, v) = Z ⌦ KWe 0 ⇣ d 0 n+1,k ⌘ ↵ (K : D(u)) : D(v) dx (3.28) l(v) = Z ⌦ KWe 0 d n+1,k ↵ 1  (u n • r) n + a (ru n , n ) 1 t n : D(v) dx
Here, d n+1,k , d 0 n+1,k , u n and n are known from the k th subiteration and from the n th time step.

Problem P1 is discretized in space as in section 3.3, by introducing a two-dimensional mesh T h made of triangular elements and the functional spaces X h and T h defined by (3.16) and (3.17). Using discontinuous Galerkin methods, the convective term for the stress tensor is approximated by

⌧ n h 2 T h such that Z ⌦ ⌧ n h : ⌧ h dx = X K2T h Z K (u n h • r) n h : ⌧ h dx + X S2L (i) h Z S [[ n h ]] : ✓ |u n h • n| 2 [[⌧ h ]] (u n h • n){{⌧ h }} ◆ ds, 8⌧ h 2 T h
and the linear form l() reads:

l h (v) = Z ⌦ h KWe 0 ⇣ d n+1,k h ⌘ ↵ 1  ⌧ n h + a (ru n h , n h ) 1 t n h : D(v h ) dx. (3.29)
With these definitions, the finite dimensional version of P1.1 writes

(P 1.1) h : Find u n+1,k+1 h 2 V h such that a(u n+1,k+1 h , v h ) = l n+1,k+1 h (v h ), 8v h 2 V 0,h where V 0,h = {v h 2 X h , v h = 0 on } V h = {v h 2 X h , v h = 0 on ext [ int }
and the computation of n+1,k+1 h is explicit using the following discrete version of (3.27) :

n+1,k+1 h = K 2 6 4 We 0 ⇣ d n+1,k h ⌘ ↵ 1 t n h We 0 ⇣ d n+1,k h ⌘ ↵ 1 [⌧ n h + a (ru n h , n h )] + We 0 d 0 h n+1,k ↵ K : D(u n+1,k+1 h ) (3.30) Subproblem P2
The damage evolution equation is written in weak form by introducing a scalar test function ' h and the discontinuous finite element space S h defined by equation 3.5 (see section 3.2). Using the broken gradient, the damage evolution equation reads

Z ⌦ h 1 t (d n+1,k+1 h d n h )' h dx + X K Z K (u n h • r)d n h ' h dx + X K X S2@K Z S ✓ 1 2 |u n h • n| [[d n h ]][[' h ]] (u n h • n)[[d n h ]]{{' h }} ◆ ds = Z ⌦ h min " 1, ⌃ t C n h 2 n+1,k+1 h , ⌃ c C n h 1 n+1,k+1 h q 2 n+1,k+1 h # 1 ! 1 T d d n h ' h dx + Z ⌦ h 1 T h ' h dx, 0 < d n+1,k+1 h  1
As d h is defined constant by element (P 0 ) and there is no flux of d across the domain boundary, the sum of the gradient of d inside elements vanishes and the jumps of d across element boundaries are non-zero only over the subset S 2 L (i) h of internal mesh element sides. The discrete variational formulation of the damage subproblem reads:

(P 2) h : Find d n+1,k+1 h 2 S h , 0 < d n+1,k+1 h  1, such that b h (d n+1,k+1 h , ' h ) = m n+1,k+1 h (' h ), 8' h 2 S h , with b h () and m h () given by b h (d n+1,k+1 h , ' h ) = Z ⌦ h 1 t d n+1,k+1 h ' h dx (3.31) m n+1,k+1 h (' h ) = Z ⌦ h " 1 t + min " 1, ⌃ t C n h 2 n+1,k+1 h , ⌃ c C n h 1 n+1,k+1 h q 2 n+1,k+1 h # 1 ! 1 T d # d n h ' h dx + Z ⌦ h 1 T h ' h dx X S2L (i) h Z S ✓ 1 2 |u n h • n| [[d n h ]][[' h ]] (u n h • n)[[d n h ]]{{' h }} ◆ ds (3.32) Subproblem P3
As the damage evolution equation, the cohesion transport equation is written in weak form using a scalar test function ' h :

Z ⌦ h 1 t (C n+1 h C n h )' h dx + X K Z K (u n+1 h • r)C n h ' h dx + X K X S2@K Z S ✓ 1 2 u n+1 h • n [[C n h ]][[' h ]] (u n+1 h • n)[[C n h ]]{{' h }} ◆ ds = 0
where the sum of the gradient of C h over the elements vanishes and the jumps in C h at inter-element boundaries are non-zero only on the subset S 2 L (i) h of internal mesh element sides. The discrete formulation of this problem writes:

(P 3) h : With u n+1 h known, find C n+1 h 2 S h such that c h (C n+1 h , ' h ) = o n+1 h (' h ), 8' h 2 S h , with c h () and o h () given by c h (C n+1 h , ' h ) = Z ⌦ h 1 t C n+1 h ' h dx (3.33) o n+1 h (' h ) = Z ⌦ h 1 t C n h ' h dx X S2L (i) h Z S ✓ 1 2 u n+1 h • n [[C n h ]][[' h ]] (u n+1 h • n)[[C n h ]]{{' h }} ◆ ds (3.34)
Using this FE spatial discretization, the complete annular Couette flow problem, (P ) h leads to the following scheme:

u n h = 0 in ⌦ h , n h = 0 in ⌦ h , d n h = 1 in ⌦ h , C n h = C 0 in ⌦ h ,
For n 0

• For k = 0, ⇣ n+1,0 h , u n+1,0 h , d n+1,0 h ⌘ = ( n h , u n h , d n h ),
• For k 0

-(P 1.1) h : With n+1,k h , u n+1,k h , d n+1,k h known, find u n+1,k+1 h 2 V h such that a(u n+1,k+1 h , v h ) = l (n+1,k+1) h (v h ), 8v h 2 V 0,h ,
with the bilinear and linear forms a() and l h () defined by equations (3.28) and (3.29) respectively and with

V 0,h = {v h 2 X h , v h = 0 on } V h = {v h 2 X h , v h = 1e ✓ on int and v h = 0 on ext }.
- -

(P 1.2) h : With u n+1,k+1
(P 2) h : With n+1,k+1 h known, find d n+1,k+1 h 2 S h , 0 < d n+1,k+1  1, such that b h (d n+1,k+1 h , ' h ) = m n+1,k+1 h (' h ), 8' h 2 S h ,
with b h () and m h () given by (3.31) and (3.32).

-Compute

res h = We 0 (d n+1,k+1 h ) ↵ 1 " n+1,k+1 h n h t + (⌧ n h + a (r h u n h , n h ) # + d1 d n+1,k+1 h e n+1,k+1 h We 0 ⇣ d n+1,k+1 0 h ⌘ ↵ K : D(u n+1,k+1 h ) , If res h < tol STOP . -Set n+1 h , u n+1 h , d n+1 h = ⇣ n+1,k+1 h , u n+1,k+1 h , d n+1,k+1 h ⌘ . • (P 3) h : With u n+1 h known, find C n+1 h 2 S h such that c h (C n+1 h , ' h ) = o n+1 h (' h ), 8' h 2 S h ,
with c h () and o h () given by (3.33) and (3.34).

In the simulations presented in this section, the tolerance, tol, was set to 10 3 times the machine precision, thereby ensuring that res h drops by several orders of magnitude with respect to its initial (k = 0) value. While using an implicit formulation for the advection, rotation and deformation terms and updating their value within the fixed-point loop would not increase the level of complexity of the numerical scheme, the explicit formulation used here allows building the transport form for the stress tensor (⌧ n h + a (ru n h , n h )) outside this loop, i.e., once per model time step. This reduces the numerical cost of the algorithm compared to the implicit case, especially over large systems. As the model time step is always very small, it does not impact the stability of the scheme (see section 3.3) and has no significant e↵ect on the simulated mechanical behaviour.

A comparison of Couette experiments

In the laboratory, Couette experiments were started from an unbroken ice place with an (approximately) uniform thickness. Within minutes, a circular fault a few millimetres wide was formed about 1 to 2 centimetres away from the inner rotating disk, concentrating most of the shearing deformation (see figure 3.11). The subsequent mechanical behaviour emanated from a competition between fracturing and refreezing of the fault.

The relative importance of the two processes was set primarily by (1) the prescribed angular frequency of the inner disk (!) and ( 2) the temperature of the cold room, which determined the rate of refreezing within the fault. These parameters were varied in a set of 25 realizations of the experiment: rotation rates ! between 1 and 80 revolutions per day (rpd) and the cold room temperature between 5 and 15 C. At high rotation rates, the mean (temporally averaged) shear stress along the fault was low. Time series of the shear stress averaged along the fault showed high frequency and low asymmetry fluctuations, similar to typical shear stress records obtained for dry granular materials (J. Weiss et al., in prep). This suggests that healing within the fault was not significant. For relatively low temperatures and slow rotation rates, shear stress time series were characterized by a succession of slow build-ups and rapid drops, which were associated with the re-activation of large portions of the fault. These strongly asymmetric cycles are visible on the records of the measured torque. Figure 3.12a shows such a time series for a realization of the experiment with ! = 1 rpd and the temperature of the cold room set to 10 C. In between large rupture events, some periods of high mean stress values and low frequency, low asymmetry fluctuations were also observed, as indicated by the dashed red box on figure 3.12, consistent with a slow dissipation of the applied loading through creep-like deformation along the fault [START_REF] Amelung | Earthquake scaling laws for creeping and non-creeping faults[END_REF]. It is important to note that the positive trend observed on this record is related to an unavoidable thickening of the ice during the experiment.

Looking at this time series, we cannot help but notice the similarity with the macroscopic stress-strain curves discussed in chapter 2 in the context of the uniaxial compression simulations, which exhibit a succession of slow stress build-ups and rapid relaxation associated with cycles of healing and damaging as well as periods of creep-like deformation with lower amplitude and asymmetry of the stress fluctuations. Yet, the mechanical behaviour observed in the laboratory experiments is highly complex and the possibilities of analyses, multiple.

Here we do not aim to perform a thorough comparison of the Maxwell-EB model and Couette experimental data. Instead, we simply seek to establish if the large-deformation Maxwell-EB model is capable of reproducing the cycles of slow stress build-up and rapid relaxation emanating from the competition between healing and fracturing observed in the laboratory Couette experiment under the simple shear loading conditions described above and if not, why.

Our choice of model parameter values for the Couette simulations is based on one realization of the experiment performed at a room temperature of 10 C and with a rotation rate of ! = 1 rpd during which the refreezing of the gouge within the fault was indeed significant. The torque record for this particular experiment was described in the previous paragraph and is represented on figure 3.12a. A characteristic healing time for this experiment is estimated as follow. Assuming a linear variation of the temperature within a thin ice layer [START_REF] Maykut | The surface heat and mass balance[END_REF], a rate of heat exchange that is proportional to the temperature di↵erence between the air and the surface of the ice plate and estimating an average coe cient for both the sensible and latent heat exchanges from records of the temperatures and ice thickness during the experiment, a simple growth model was obtained for the ice plate (J. Weiss et al., in prep). According to this model, the thickening of the very thin plate was almost linear in time. The time required to grow a 4 mm plate over the initially ice-free basin at 10 C was of about t h = 10 4 s. This estimate of the characteristic healing time is in agreement with the power spectrum density analysis performed in section 2.4.3 in the context of the uniaxial compression Maxwell-EB simulations.

In these analyses, we indeed identified a range of frequencies for the healing-damaging cycles (equivalent to the stress build-up and relaxation cycles) that is 1 to 3 times larger than the frequency associated with the prescribed characteristic healing time (see figure 2.12a). From figure 3.12c, we can estimate visually the frequency of stress build-up and relaxation cycles: we count 4 such cycles in ⇠ 6 hours, equivalent to a frequency of 2.0 • 10 4 s 1 . This is e↵ectively twice as high as the inverse estimated healing time (10 4 s 1 ). This "backof-the-envelope" calculation increases our confidence in our assumption that the stress build-up and relaxation cycles observed in the Couette experiment are indeed controlled by the rate of healing of the shearing fault and that healing can be adequately represented in the Maxwell-EB model using the simple parameterization introduced in section 2.2.3.

With the angular frequency of the inner boundary set to ! = 1 rpd, the shearing deformation rate within the initially undamaged ice plate is of U L = 2⇡Rint 24•3600 1

Rext Rint ⇡ 5.0 • 10 5 s 1 . In order to approach the limit t (hours) of a strictly elastic undamaged ice plate, while ensuring the convergence of the numerical scheme, we set the initial Weissenberg number to We 0 = 0.1, based on the results of the sensitivity analyses of section 3.3. This gives 0 = 0.1 ⇥ L U ⇡ 2 • 10 3 s. The damage parameter ↵ is set to 4 as in the uniaxial compression experiments. The order of magnitude of the values of cohesion is chosen based on the duration of the initial elastic loading in the laboratory experiment (⇠ 20 minutes, see figure 3.12). The local value of C is re-drawn from the same uniform distribution of values every time a given element becomes overcritical, as in annealed disorder. The model time step is set as follow.

T (Nm) (a) (b) (c)
Considering for instance an unstructured mesh with constant element size x (i.e., independent of r) and N = 10, such that x = Rext N = 0.05 m and a speed c of propagation of shear elastic waves of 2 • 10 3 ms 1 , consistent with estimations of the Young modulus of laboratory fresh ice (e.g., Gammon et al., 1983, 10 GPa), the characteristic time of propagation of damage in the ice plate, t d , would then be of at most x c ⇡ 2.5 • 10 5 s. If choosing t = t d , on the order of 10 8 model time steps would be required to simulate one hour in real time. With a rotation rate of ! = 1 rpd, one quarter of a revolution of the inner disk would represent 10 9 time steps. Considering instead a mesh grid with resolution x varying as r 2 , as used here, one quarter of a revolution would require on the order of 10 10 model time steps and would represent several months of runtime.

Here, we therefore choose a model time step based on a reasonable runtime for our simulations. We set t = 1 s. Doing so, the characteristic time for the propagation of damage, t d , can be set in di↵erent ways. We explore 2 possibilities : 1. t d = t with t = 1 s, which amounts to assuming a slow propagation of the damage within the ice plate,

2. t d = c
x , with x the size of the smallest element of our mesh grid and c representative of the speed of propagation of shear elastic waves in laboratory fresh ice. In this case, t d ⇡ 1.0 • 10 6 s and hence All model simulations are started from an initially undamaged ice plate with uniform elastic modulus and viscosity. Consistent with the behaviour observed in the laboratory experiment, the time series of the (adimensional) torque on the inner boundary, , shows an initially linear-elastic loading phase, followed by an abrupt drop of the stress (see figure 3.13). In the case of t d ⌧ t, the numerical scheme does not converge at this point and simulations stop. This can be trivially explained by the large separation of scales between t and t d . According to the time-discretized equation (3.26) for the level of damage, if t t d (or equivalently, ˜ t T d ), d drops to near zero values over damaged elements within very few time steps, hence rapidly generating numerical instabilities. In this case, the model does not resolve the progressive propagation of the damage. Here we therefore only discuss "slow damage" simulations, in which t d is set equal to t and the model does resolve the progressive damage mechanism. In these simulations, the abrupt stress relaxation is associated with a sharp increase in the damage rate and with a propagation of damage along a roughly circular feature a few grid cells away from the inner boundary which concentrates virtually all of the simulated deformation (see figure 3.13a and 3.13b). Consistent with the Couette experiment, the stress within the partially damaged material subsequently build-ups to values lower than the pre-macro-rupture maximum. The torque time series shows some macroscopic fluctuations which correlate with small increases in the damage rate and are associated with the reactivation of small portions of the fault (see figure 3.13c).

t d ⌧ t. (a) (b) (c)
As the spatial extent of these damaging events is relatively small, the torque in the Maxwell-EB Couette simulations never drops to low values after the first rupture event, as opposed to the torque measured in the laboratory (see figure 3.12). Another obvious di↵erence between the simulation and experiment results pertains to the symmetry of the torque signal. In the model, no abrupt stress relaxation is simulated after the first macro-rupture event. The time series of shows almost symmetric increases and decreases. The distribution of torque fluctuations, t after the first rupture (not shown), is characterized by a (negative) skewness of -1.2, which is indeed low compared to the (negative) skewness of -6 computed for the distribution of the macroscopic stress fluctuations, m t , in the uniaxial compression simulations presented in 2 (see figure 2.11, section 2.4.3). The frequency of the macroscopic stress build-up and relaxation cycles also appears much higher in the simulations than in the laboratory experiment. The average power spectral density of the damage rate time series, calculated for 5 simulations initialized with di↵erent fields of cohesion (see figure 3.14a) shows a factor of 5 separation between the frequency of these cycles (indicated by the dashed red line) and that associated with the prescribed value of t h , suggesting a decoupling between the small-extent damage events and the healing mechanism. Beyond the frequency of these cycles and up to f ⇡ 10 3 , the PSD shows a power-law decrease consistent with a long-term memory of the damage rate within the material. One can think of two possible adjustments to the model parameters that could result in a more pronounced asymmetry of the simulated macroscopic behaviour and more widespread damaging events. The first consists in increasing the value of the poorly constrained damage parameter, ↵, which in the case of small-deformation simulations was shown to result in more abrupt drops of the internal stress at the onset of damaging and to a higher temporal localization of damaging events (see section 2.5). In the present case however, using a larger ↵ does not help obtaining a better agreement with the observed macroscopic behaviour. The second adjustment consists in decreasing the prescribed time of healing, so that to allow for a widespread refreezing and a rapid increase of the internal stress within the simulated fault. In the context of the laboratory experiment, this amounts to lowering the cold room temperature, which was indeed shown to induce larger fracturing events. Here, we try simulating a more rapid freeze-up of the shearing fault by decreasing our estimated value of t h by a factor of 10 (t h = 10 3 s). Figure 3.16 indicates that this does amplify the simulated fluctuations of both the torque and damaging activity and somewhat increases the temporal asymmetry of the torque signal (skewness of the t distribution of -2). Fields of the level of damage and deformation rates indicate that in this case the highly damaged zone almost completely heal between damaging events, which "forces" the system to enter phases of extensive damaging activity and to form a new fault. The power spectral density of the damage signal (see figure 3.14b) shows a clear and narrow peak at a frequency that is about twice as high as the prescribed healing time, consistent with the fact that widespread healing is the main driving mechanism of these healing and damaging cycles.

Discussion

Here we elaborate on the possible reasons why the Maxwell-EB model with the present set of numerical and physical parameters does not fully reproduce the mechanical behaviour observed in the laboratory.

• Considering the adimensional form of the constitutive and damage equations (2.26 and 2.27), the five parameters controlling the simulated mechanical behaviour in these experiments are (1) the e↵ective Weissenber number, We = We 0 d ↵ 1 , (2) the (adimensional) model time step ˜ t, (3) the characteristic time of damage, T d , (4) the damage parameter, ↵, and (5) the characteristic healing time, T h .

As already discussed in chapter 2 (section 2.2.4), for We T d the relaxation time for the internal stress is much larger than the time of propagation of the damage : permanent deformations are small and elastic interactions are e↵ective in redistributing the internal stress within the material. For We ⇠ T d , deformations are essentially permanent, elastic interactions and the propagation of damage are hindered.

The rate at which the transition between elastic and permanent deformation occurs depends on the initial separation between the values of T d and We 0 .

On the one hand, in the present simulations the time step is chosen conveniently large and T d is set equal to ˜ t so that to resolve the propagation of damage. On the other hand, the value of We 0 is limited by the numerics. Therefore, the initial separation between We 0 and ˜ t is such that the e↵ective value of We readily drops to ⇠ T d as the material becomes damaged. Elastic interactions are by this fact inhibited and the extent of the damaging, restricted. Consistent with this idea, increasing the damage parameter ↵, which controls the rate of decrease of We with the level of damage, does not lead to larger damaging events but instead to an even "softer", viscous-like behaviour.

Put another way, the propagation of damage in these Couette simulations is too slow. Considering the value chosen for t d (1 s), the maximal resolution of the mesh grid, x ⇡ Rext 10 r 2 at r = R int (near the inner boundary) and the e↵ective speed of propagation of damage near the inner boundary c = x t d , we can estimate the minimum time for propagation of damage along the circular fault as 2⇡Rint c ⇡ 600 s. This time is smaller than the undamaged value of relaxation time for the internal stress ( 0 ⇡ 2 • 10 3 s), consistent with the fact that the model does represent the initial propagation of the macro-rupture with large damage rates and a first abrupt relaxation of the stress. However, considering an average level of damage of the material of d = 0.5 in in the vicinity of a partially healed fault near the point of rupture and ↵ = 4, the e↵ective relaxation time, 0 d ↵ 1 (⇡ 300 s) is then less then the time for the propagation of damage around the inner boundary. Viscous dissipation therefore dominates the mechanical behaviour before damage has time to propagate along the fault, hence explaining the small extent of damaging events after the first macro-rupture.

To better resolve the brittle character of ice deformation, a larger separation of scales between the time for the propagation of damage and the relaxation time for the internal stress in the model is necessary.

As the maximum value of We 0 is set by the construction (convergence) of our numerical scheme, in the current implementation of the Maxwell-EB model this implies using a smaller time step and hence long, unpractical simulation times. The representation of laboratory-scale experiments with the Maxwell-EB model therefore calls for a reformulation of the "discrete-like" progressive damage mechanism in a truly continuous form. This would (1) eliminate the requirement of an explicit scheme for the damage evolution equation and the necessity of a small model time step, (2) allow exploiting fully-implicit and higher-order numerical methods and (3) perhaps exploring numerical schemes that were designed to handle larger values of the Weissenberg number in viscoelastic flow problems, for instance those based on a log-conformation formulation (e.g., [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF][START_REF] Saramito | On a modified non-singular log-conformation formulation for johnson-segalman viscoelastic fluids[END_REF].

• Another important point related to the choice of model parameters and characteristic numbers in these Couette simulations concerns the relative value of the (adimensional) healing time T h and of the Weissenberg number, We. With t h = 10 4 s, T h is always larger than We, which means that the rate of relaxation of the internal stress in permanent deformations is always larger than the rate of healing of the material. Using this value of the characteristic healing time, healing does not seem e cient in reincreasing the mechanical strength over damaged areas and slowing deformation rates. With t h = 10 3 s, the rate of healing is about two times lower than the rate of viscous dissipation within the undamaged material ( 0 ⇡ 2 • 10 3 s). The simulations presented here demonstrate that healing is then e↵ective in lowering levels of damage and allow the stress to significantly build-up within the material so that to trigger subsequent large damaging events. Although the balance between the rate of healing and damaging also depends on the applied driving rate, a necessary condition to represent the competition between the fracturing and refreezing within the ice plate might be t h < 0 or, equivalently, T h < We 0 .

• In the case of large deformations simulations and relatively high driving rates, we can think that the di↵usion of the strong gradients and local weakening of the material through advection can also smooth elastic interactions and the propagation of damage and hence lead to a "softer" macroscopic behaviour.

However, as the e↵ect of numerical di↵usion is hard, perhaps impossible, to quantify here, we can only suppose it is of second order in the present problem.

• Finally, some physical processes might also be missing in our rheological framework that prevent the Maxwell-EB model from representing adequately the mechanical behaviour of a shearing fault. Based on the Couette simulations presented here, it appears that healing, at a rate consistent with the laboratory conditions, alone is not su cient for the stress to build up to critical values over large enough portions of the fault so that to trigger large-scale damaging events. In the laboratory, friction along the fault acts to slow the movement and hence might allow for more e cient refreezing. Because of the irregular geometry of the fault, frictional sliding can also concentrate the stress locally by locking some portions of the fault, thereby initializing fracturing events. In the Maxwell-EB framework, the friction associated with the shearing deformation along activated faults is not represented. The inclusion of a threshold on the value of the internal stress for the deformation of damaged elements might be a necessary ingredient in simple shear conditions. This point is further discussed in chapter 5.

Concluding remarks

In this chapter a numerical scheme for the large-deformation Maxwell-EB model based on discontinuous Galerkin methods was presented. Three series of numerical experiments with di↵erent level of complexity were performed in the context of large deformations simulations.

In the first series of tests, a solid body rotation problem was solved, hence only advective processes were represented. These have allowed evaluating the capability of the discontinuous Galerkin method in handling the advection of fields presenting extreme gradients and a high level of discontinuity, such as the fields of the level of damage, d, and of the stress, , simulated by the model. Using polynomial approximations of order 0 for the advected fields, the DG scheme coincides with a finite volume method with upwinding. In this case the numerical scheme is di↵usive and di↵usion decreases slowly with increasing spatial resolution. In the current implementation of the Maxwell-EB model, all advected fields are approximated constant by element. Di↵usion is therefore expected to be important, and most important when and where the level of damage is high and internal stresses are low. Conversely, where internal stresses are high, it might be mitigated by the tendency of the model to "re-localize" deformation and damaging at the scale of the model element. This point is further discussed in the next chapter.

In the second series of simulations, the advection, rotation and deformation terms for the stress tensor were all accounted for and the Maxwell-EB system of equations was solved in the case of a homogeneous, undamaged material under simple shear conditions. These have allowed evaluating the highest value of the Weissenberg number that the numerical scheme can handle while ensuring the convergence towards a stationary solution and to discuss its impact in terms of unphysical viscous dissipation of the stress in the Maxwell-EB model.

Finally, the full, large-deformation Maxwell-EB framework was implemented in the context of a Couette flow experiment which allowed comparing the simulated mechanical behaviour to that observed in the laboratory for a naturally formed shearing fault in fresh-water ice. The model was found to simulate adequately the failure of the ice plate and the formation of the shearing fault. The subsequent healing and fracturing within the fault was however underrepresented in the model: in particular, the spatial extent of damage events was limited and the amplitude of the associated stress relaxation was small. The simulated behaviour was explained in terms of numerical and possible physical limitations of the Maxwell-EB model. An important, perhaps the most important, limitation in the context of these experiments stems from our treatment of the progressive damage mechanism, which requires using an explicit time stepping scheme and setting the model time step to a value comparable to that of the characteristic time of propagation of damage t d in the simulated material. In the case of ice at the laboratory scale, this requirement translates into unpractical simulation times. A slower propagation of damage in the material therefore needs to be assumed, which results in a too small separation between the rate of damage propagation and the rate of viscous relaxation of the internal stress within the material where the level of damage is low and a too soft, viscous-like behaviour.

In the case of sea ice on geophysical scales, setting the model time step equal to the characteristic time of propagation of damage places the time step between a couple of seconds and a couple of minutes, for spatial resolutions of 1 km to 100 km respectively. Although small, such a time step is comparable to that of current ocean models and translates into reasonable simulation times. In this case, the separation between t d and the viscous relaxation time within a relatively undamaged ice cover is large enough for the Maxwell-EB model to simulate adequately the brittle character of sea ice deformation. This point is further discussed, and demonstrated, in the next chapter.

Coulombic granular models and obtained deformation patterns over Bering Strait and the Amundsen Gulf that were in good agreement with satellite (Landsat) imagery. Stable ice bridges have also been successfully

reproduced by plastic-type sea ice models, providing the prescribed plastic yield criteria allowed some cohesion or tensile strength [START_REF] Ip | Numerical investigation of di↵erent rheologies on sea-ice dynamics[END_REF][START_REF] Hibler | Sea-ice arching and multiple flow states of arctic pack ice[END_REF]. Recently, [START_REF] Dumont | Modeling the dynamics of the north water polynya ice bridge[END_REF] were able to simulate the formation of ice arches in both an idealized and a realistic representation of Nares Strait using the Viscous-Plastic model of [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF]. As their model was based on the widely used elliptical yield curve, it did not account for uniaxial or biaxial tensile strength of the ice. Stable ice bridges and flow stoppage were obtained by decreasing the ellipticity of the yield curve below the standard value (2, [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF], hence by increasing the area of the yield envelope in the second and fourth quadrants in the principal stress space, corresponding to a cohesive state.

Channel flow simulations have not yet been performed using elastic-brittle models. Hence such experiments constitute an interesting test case of their mechanical behaviour. Here we therefore establish the capability of the Maxwell-EB model to simulate the formation of stable ice bridges and self-obstruction to flow under strong along-channel wind conditions. Moreover, while other rheological models have been shown to simulate both the occurrence of ice bridges and flow stoppage, it is not at all clear if they are also able to represent the presence of multiple, fine arch-like leads within and upstream of the channel [START_REF] Sodhi | Ice arching and the drift of pack ice through restricted channels[END_REF], as observed from satellite imagery (e.g., see figure 4.1b). With its ability to represent the propagation of damage along highly localized faults, the Maxwell-EB model might be more suited to simulate these features. Although using a realistic configuration of Nares Strait would not represent much technical di culty in the present model, here

we retain an idealized, symmetric domain geometry and simple, homogeneous forcing conditions to facilitate the analysis of the dynamical behaviour and, in particular, of the states of stress within the channel.

Once a stable ice bridge forms, the ice downstream is expected to detach from the bridge and be driven out of the channel without any mechanical resistance. Extreme gradients are therefore expected to arise in the ice velocity as well as in the ice thickness field (see figures 4.1b and 4.2a). These represent severe conditions for testing the limits of the numerical scheme in terms of its di↵usivity and capability in handling discontinuities within the simulated fields.

The complete picture

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.

In this case, the complete dynamical model is given by the following system of equations:

1. The momentum equation:

⇢h  @u @t + (u • r)u = A(⌧ a + ⌧ w ) ⇢hf k ⇥ u ⇢hgrH + r • (h ), (4.1)
with ⇢, the ice density, h, the mean ice thickness, A, the ice concentration, defined as the fraction of a model grid cell covered by ice, ⌧ a and ⌧ w , the air and water drags, ⇢hf k ⇥ u, the Coriolis pseudo-force with f , the Coriolis parameter and k, the upward unit vector normal to the ice surface, and ⇢hgrH, the force due to gradients in the sea surface dynamic height, H, with g the gravitational acceleration.

This last term is often expressed alternatively in terms of the geostrophic ocean current velocity, u w , as ⇢hf k ⇥ u w . We assume the internal stress to be homogeneously distributed throughout the ice depth h and, as in the models of [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], 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 and mean ice thickness:

@h @t + r • (hu) = S h , h 0, (4.2) @A @t + r • (Au) = S A , 0  A  1. (4.3)
Mass conservation is ensured by this first equation for the evolution of the mean ice thickness, hence neglecting the compressibility e↵ect on the ice density, ⇢. The mean thickness is assumed to be the weighted sum of two ice categories: thick ice, with thickness h thick = h A , and open water, with thickness h thin = 0 m [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF]. The terms S h and S A represent the mechanical redistribution and thermodynamic sources and sinks of ice. In the following implementation of the model, thermodynamics processes are not represented.

Mechanical redistribution is accounted for following [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF] and [START_REF] Thomson | Short-term ice motion modeling with application to the beaufort sea[END_REF]. If over a given element, the ice concentration A exceeds unity, the excess concentration, max[0, (1 A)], is used to increment the thickness of thick ice over that element, and the ice concentration is reset to 1. With this approach, the mechanical source (h + ) and sink (A ) terms of h and A read:

h + = max[0, (1 A)]h (4.4) A = max[0, (A 1)] (4.5)
Again, our formulation of the mechanical redistribution of ice neglects any compressibility e↵ects. Based on the order of magnitude di↵erence between elastic and permanent deformations within sea ice (see chapter 1), the error associated with this approach is expected to be small. Comparison of channel flow simulations in which equations (4.2) and (4.3) were solved using the total versus the permanent part of the deformation (i.e., of u) only indeed showed no significant di↵erence in terms of the simulated spatial distribution of the ice thickness and concentration.

3. The constitutive law:

1 E  @ @t + (u • r) + a (ru, ) + 1 ⌘ = K : "(u).
As in chapter 3, the upper convected form of the objective tensor derivative is used here (a = 1). An additional relation of E and ⌘ in terms of the ice concentration A is introduced, such that :

E = f 1 (E 0 , ⌘ 0 , ⌘ min , d) exp [ c ⇤ (1 A)], (4.6) ⌘ = f 2 (⌘ 0 , ⌘ min , d) exp [ c ⇤ (1 A)], (4.7)
with f 1 and f 2 , the functional dependance of E and ⌘ on the level of damage of the ice cover d given by (2.24) and (2.23) respectively and c ⇤ , a non-dimensional parameter characterizing the dependence of the mechanical parameters on the ice concentration. This exponential function of the ice concentration allows the value of both mechanical properties to be maximal when the ice concentration is of 100% (A = 1) and to decrease rapidly when leads open and A drops. This function 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] and the e↵ective elastic sti↵ness in the elasto-brittle models of Girard et al. (2010b) and of [START_REF] Bouillon | Presentation of the dynamical core of neXtSIM a new sea ice model[END_REF]. In the case of ⌘, this dependence on ice concentration is compatible with the rapid decay of the apparent viscosity of granular media when decreasing their packing fraction from the close-packed limit [START_REF] Aranson | Patterns and collective behavior in granular media: Theoretical concepts[END_REF]. In the present implementation of the Maxwell-EB model, the value of the parameter c ⇤ is constant and set to 20 as in the model of [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF] for both E and ⌘, but this parameterization could eventually be refined.

4. The equation for the evolution of damage:

@d @t + (u • r)d = ✓ min  1, t 2 , c 1 q 2 1 ◆ 1 t d d + 1 t h , 0 < d  1, (4.8)
in which a constant healing rate is assumed. The damage criterion is defined by equations (2.8) and

(2.11), with q, c and t given by equations (2.9), (2.10), (2.12) in terms of the cohesion, C, and of the constant internal friction coe cient, µ.

5. The transport equation (2.13) handles the advection of the passive field of C.

Channel flow simulations

The domain for the channel flow simulations presented here represents a 120 km wide rectangular basin, converging into a 40 km wide, 40 km long channel (see figure 4.2b). It is symmetric with respect to the y axis.

The geometry and dimensions are similar to that used in the idealized simulations of [START_REF] Dumont | Modeling the dynamics of the north water polynya ice bridge[END_REF] and conceived to be roughly consistent with the shape of the constriction between Kane Basin and Smith Sound, in Nares Strait. The dynamics described in the next section is not specific to this geometry. Indeed, simulations were also performed over di↵erent domains (narrower, longer channels, smaller basins) and produced similar results.

A reduced form of the momentum equation (4.1) is solved here, in which the advection and acceleration terms are neglected. To retain symmetry in forcing conditions, the Coriolis term is also set to zero. The prescribed wind forcing is made as simple as possible to facilitate the analysis. A uniform northerly wind channels, hoppers, and silos. This last problem has been extensively studied using the granular theory for bulk solids [START_REF] Richmond | Limiting spans for arching of bulk materials in vertical channels[END_REF]Walters 1973;Savage and Sayed 1981). In contrast to fluids, bulk solids can transmit shear stresses while at rest. Cohesive materials (e.g., damp sand), which are able to support higher static shear stresses, are capable of forming a self-obstruction to flow [START_REF] Walker | An approximate theory for pressures and arching in hoppers[END_REF]. Accordingly, ice arching has been observed and modeled assuming sea ice to be a plastic (discontinuous flow under continuous forcing) and cohesive (with the ability to maintain its integrity while submitted to tensile forces) material. [START_REF] Sodhi | Ice arching and the drift of pack ice through restricted channels[END_REF] compared ice arches forming in Bering Strait and in the Amundsen Gulf with a cohesive Mohr-Coulomb granular rheology and found a good correspondence between the modeled ice arch form and the observed profile. He followed the analytical analysis of Morrison and Richmond (1976) that relates the arching process with the cohesive strength of the granular material. [START_REF] Ip | Numerical investigation of di↵erent rheologies on sea-ice dynamics[END_REF] studied the problem of ice arching in converging channels with different rheologies and showed that plastic yield curves lying on or crossing the principal stress axes allow ice arch formation, given the ade-quate loading and thickness. His main conclusion is that only cohesive materials are able to form arches. The elliptical yield curve for sea ice, first introduced by [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF], has become the most widely used dynamic sea ice model in climate and process studies. It represents a cohesive material because a part of it covers the second and fourth quadrants of the principal stress space (Fig. 2). It does not allow pure tension, which occurs when both principal stress components are positive. Shear strength relative to compressive strength is scaled by the major to minor axis ratio parameter e. Cohesion increases along with shear strength as e decreases. A value of e 5 2 was originally chosen by [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF] to approximately match the ratio of energy dissipation in shear to energy dissipation in sea ice ridging calculated by [START_REF] Rothrock | The energetics of the plastic deformation of pack ice by ridging[END_REF]. [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF] showed that decreasing e stiffens the ice flow throughout the Arctic Basin but does not significantly alter spatial patterns of modeled ice thickness. On the other hand, such a change has a drastic impact on both the flow and spatial distribution of ice thickness in smaller enclosed areas. Kubat et al. (2006) used the ellipse and showed that modifications of the shear strength (cohesion) may stress, ⌧ a , is applied and increased steadily between 0 and 1 Nm 2 over a period of 24 hours, then held constant. Considering a quadratic wind drag of the form

⌧ a = ⇢ a C da |u a |u a , (4.9)
with ⇢ a = 1.3 kg m 3 , the surface air density and C da , the air drag coe cient, commonly set to 1.2 • 10 3 in sea ice models following [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF], this corresponds to a maximum wind speed of ⇠ 22 ms 1 (⇠ 82 km h 1 ). These conditions are thought as to represent a strong, along-strait, winter wind episode in Nares Strait [START_REF] Samelson | Evidence for atmospheric control of sea-ice motion through nares strait[END_REF]. The ocean is at rest (u w = 0), hence the force associated to gradients in the sea surface dynamic height is also zero. The ocean drag is given by the following quadratic formula

⌧ w = ⇢ w C dw |u|u, (4.10)
where ⇢ w = 1027 kg m 3 is the density of sea water and C dw is the drag coe cient, set to 5.5 1980). Following [START_REF] Gray | A two-dimensional model for the dynamics of sea ice[END_REF], both drag terms are weighted by the local ice concentration, A, to account for the fraction of open water within a grid cell.

• 10 3 (McPhee,
The ice is initially undamaged (d = 1) and has a uniform thickness of 1 m, representative of thin first year ice in Nares Strait [START_REF] Haas | Comparison of the sea-ice thickness distribution in the lincoln sea and adjacent arctic ocean in 2004 and 2005[END_REF]. Contrary to [START_REF] Dumont | Modeling the dynamics of the north water polynya ice bridge[END_REF], the domain is initially completely covered with ice (A = 100%), so that to avoid prescribing the location of ice bridges. Mechanical parameter values are based on measurements within sea ice. A speed of propagation of shear elastic waves of c = 500 ms 1 is considered [START_REF] Marsan | Low frequency bursts of horizontally-polarized waves in the arctic sea-ice cover[END_REF], such that the undamaged elastic modulus,

E 0 = 2c 2 (1 + ⌫)⇢,
is of approximately 5.0 • 10 8 Pa. Poisson's ratio is set to ⌫ = 0.3 [START_REF] Timco | A review of the engineering properties of sea ice[END_REF]. The undamaged relaxation time is set to 10 7 s (⇠ 115 days), a value that allows the numerical scheme to converge, while also ensuring that unphysical viscous dissipation over low damage areas of the ice cover is insignificant. The time for healing, t h , is set to 5 • 10 5 s based on the ice growth rate estimates of [START_REF] Petrich | Formation and structure of refrozen cracks in land-fast first-year sea ice[END_REF]

0 = A (⌧ a ⇢ w C dw |u|u) + r • (h ) in ⌦⇥]0, +1[, 0 d ↵ 1  @ @t + (u • r) + a (ru, ) + d1 de = ⌘ 0 d 0 ↵ exp[ c ⇤ (1 A)]K : D(u) in ⌦⇥]0, +1[, @d @t + (u • r)d = ✓ min  1, c 2 , t 1 q 2 1 ◆ 1 t d d + 1 t h , 0 < d  1, in ⌦⇥]0, +1[, @C @t + (u • r)C = 0 in ⌦⇥]0, +1[, @h @t + (u • r)h = h(r • u) in ⌦⇥]0, +1[, @A @t + (u • r)A = A(r • u) in ⌦⇥]0, +1[,
with c and t defined in terms of C as (2.10) and (2.12) respectively,

d 0 ↵ = (1 ⌘min ⌘0 )d ↵ + ⌘min ⌘0
and with initial conditions

u(t = 0) = 0 ms 1 in ⌦, (t = 0) = 0 Pa in ⌦, d(t = 0) = 1 in ⌦, h(t = 0) = 1 m in ⌦, A(t = 0) = 1 in ⌦, C(t = 0) = C 0 (Pa) in ⌦,

Time discretization

Non-linearity in this version of the model is seven-folds: it arises through (1) the advection terms, (2) the rotation and deformation terms, a , in the constitutive equation, (3) the coupling between the damage and constitutive equations, (4) the coupling between the momentum and the mass (thickness) and ice concentration conservation equations, (5) the coupling between the constitutive and ice concentration equations, (6) the coupling between the cohesion transport and damage equations ( 7) the quadratic formulation of the water drag term in the momentum equation. The system of equations is discretized in time as follow.

• A semi-implicit scheme is used for the momentum equation. The internal stress term is discretized using an implicit scheme while the water drag term is linearized as ⌧ w = ⇢ w C dw |u n |u n+1 . The ice thickness and concentration are both taken at the n th time step.

• As in the Couette flow simulations, a semi-implicit scheme is used for the constitutive equation in which the advection, rotation and deformation terms are estimated using the fields of velocity and internal stress at the n th model time step. The ice concentration is also taken at the previous model time step.

• The damage evolution equation is discretized using a semi-implicit scheme as in the Couette simulations.

• The cohesion, ice thickness and concentration equations use an explicit scheme and the value of the ice velocity at the (n + 1) th time step.

Discretizing the time, t, such that t n = n t, with t > 0 and n = 0, 1, 2, ..., the system of equations reads:

A n ⌧ a ⇢ w C dw |u n |u n+1 + r • (h n n+1 ) = 0 0 (d n+1 ) ↵ 1  n+1 n t + (u n • r) n + a (ru n , n ) + d1 d n+1 e n+1 = ⌘ 0 (d 0n+1 ) ↵ exp[ c ⇤ (1 A n )]K : D(u n+1 ) d n+1 d n t + (u n • r)d n = ✓ min  1, n t n+1 2 , n c n+1 1 q n+1 2 1 ◆ 1 t d d n + 1 t h , 0 < d n+1  1 C n+1 C n t + (u n+1 • r)C n = 0 h n+1 h n t + (u n+1 • r)h n = h n (r • u n+1 ) A n+1 A n t + (u n+1 • r)A n = A n (r • u n+1 )
The channel flow problem (P) is divided into four linear subproblems. A fixed point algorithm is used in which the momentum and constitutive equations are first solved simultaneously, as in the small deformation and Couette flow numerical schemes (subproblem P1). The damage equation is then solved exactly for d n+1 using the fields of n+1 and u n (subproblem P2). The value of d is updated in the constitutive equation and the two computations are iterated until the residual of the constitutive equation drops below a prescribed tolerance.

The conservation equation for the ice thickness and concentration are solved using u n+1 and both fields are adjusted for mechanical redistribution (subproblem 4). The cohesion transport equation is then solved using u n+1 and the local damage criterion (i.e., c and t ) is updated (subproblem P5).

Using the superscript k = 0, 1, 2, ... for the sub-iterations, the fixed point algorithm for solving the momentum, constitutive and damage equations reads :

For k = 0, let n+1,0 , u n+1,0 , d n+1,0 = ( n , u n , d n ),
For k 0

• (P1) Find n+1,k+1 and u n+1,k+1 such that

A n ⌧ a ⇢ w C dw |u n |u n+1,k+1 + r • (h n n+1 ) = 0 (4.11) 0 (d n+1 ) ↵ 1  n+1 n t + (u n • r) n + a (ru n , n ) + d1 d n+1 e n+1 = ⌘ 0 (d 0 n+1 ) ↵ exp[ c ⇤ (1 A n )]K : D(u n+1 ) (4.12)
and with ,k+1 , such that 0 < d n+1,k+1  1 and

u n+1,k+1 = 0 on left [ right , n+1,k+1 • n = 0 on in [ out , d n+1,k+1 = d on , h n+1,k+1 = h on , A n+1,k+1 = A on , n+1,k+1 = on , C n+1,k+1 = C on . • (P2) Find d n+1
d n+1,k+1 d n t + (u n • r)d n = ⇣ min h 1, n t 2 n+1,k+1 , n c n+1,k+1 1 q n+1,k+1 2 i 1 ⌘ 1 t d d n + 1 t h .
(4.13)

• Stopping criterion : compute

res = 0 (d n+1,k+1 ) ↵ 1  n+1,k+1 n t + (u n • r) n + a (ru n , n ) + d1 d n+1,k+1 e n+1,k+1 ⌘ 0 (d 0 n+1,k+1 ) ↵ exp[ c ⇤ (1 A n )]K : D(u n+1,k+1 )
If res < tol STOP , and set n+1 , u n+1 , d n+1 = n+1,k+1 , u n+1,k+1 , d n+1,k+1 .

While the value of the ice thickness and concentration could also be updated as part of the fixed point iteration, modulations of these fields and their impact on the momentum balance are slow-evolving compared to that associated with damaging. Hence variations of h and A do not a↵ect the stability of the numerical scheme nor the rate of convergence of the solution and can be estimated outside of the FP iterations together with the transport of the field of cohesion. This approach has also no significant impact on the model results.

As in previous numerical schemes, subproblem P1 is solved in two steps. The momentum equation (4.12)

is solved for u n+1,k+1 by substituting for the following expression of n+1,k+1 :

n+1,k+1 = K " 0 d n+1,k ↵ 1 t n 0 d n+1,k ↵ 1 [(u n • r) n + a (ru n , n )] +⌘ 0 d 0n+1,k ↵ exp[ c ⇤ (1 A n )]K : D(u n+1,k+1 ) i (4.14) with K = h 0 d n+1,k ↵ 1 1 t + d1 d n+1,k e i 1
. The internal stress n+1,k+1 is calculated exactly by substituting u n+1,k+1 into (4.14).

Variational formulation and discontinuous Galerkin FE approximation

Subproblem P1

Introducing the continuous functional space V = {v 2 H 1 (⌦) 2 ; v = 0 on left and right }, the weak form of the time-discretized momentum equation writes:

Z ⌦ ⌧ a ⇢ w C dw |u n |u n+1,k+1 • v dx Z ⌦ h n n+1 : D(v) dx + Z @⌦ h n n+1 • n v ds = 0,
with the boundary integral vanishing because v = 0 on the lateral boundaries ( left,right ) and because of the Neumann condition • n = 0 on the top and bottom boundaries ( in , out ). Substituting for n+1,k+1 given by (4.14), the variational formulation of the linear momentum subproblem reads:

(P1.1) : Find u n+1,k+1 2 V such that a(u n+1,k+1 , v) = l(v), 8v 2 V,
with the bilinear form a() and linear form l() defined for all u, v 2 H 1 (⌦) 2 as

a(u, v) = Z ⌦ K⌘ 0 d 0n+1,k ↵ exp [ c ⇤ (1 A n )]h n [K : D(u)] dx + ⇢ w C dw Z ⌦ |u n |u • v dx, (4.15) l(v) = Z ⌦ ⌧ a • v dx + Z ⌦ K 0 d n+1,k ↵ 1 h n  (u n • r) n + a (ru n , n ) 1 t n : D(v) dx
and with d n+1,k , d 0 n+1,k , h n , A n , u n and n known from the previous fixed point or model iteration.

The problem is discretized in space by introducing the two-dimensional mesh T h , the functional space X h of P 1 elements for the velocity field given by (3.16) (section 3.3), the functional space T h of P 0 discontinuous elements for the internal stress tensor given by (3.17) and the functional space S h of P 0 discontinuous elements for the damage level field and other mechanical properties given by (3.5). We define the discontinuous Galerkin approximation of the convective term for the stress tensor, ⌧ n h 2 T h , as

Z ⌦ h ⌧ n h : ⌧ h dx = X S2L (i) h Z S [[ n h ]] : ✓ |u n h • n| 2 [[⌧ h ]] (u n h • n){{⌧ h }} ◆ ds + Z @⌦ h max(0, u n h • n) n h : ⌧ h , 8⌧ h 2 T h .
where the last term handles the flux of through the upstream part of the boundary [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF][START_REF] Saramito | E cient C++ finite element computing with Rheolef[END_REF]. With these definitions, the finite dimensional version of P1.1 reads:

(P 1.1) h : Find u n+1,k+1 h 2 V h such that a(u n+1,k+1 h , v h ) = l n+1,k+1 h (v h ), 8v h 2 V h where V h = {v h 2 X h , v h = 0 on left [ right }
and the discretized form of l() is

l n+1,k+1 h (v) = Z ⌦ h ⌧ a • v h dx + Z ⌦ h K 0 ⇣ d n+1,k h ⌘ ↵ 1 h n h  ⌧ n h + a (ru n h , n h ) 1 t n h : D(v h ) dx Z @⌦ h h n h max(0, u n h • n) : D(v h ) ds (4.16)
with n h 2 T h , the piecewise approximation of n , u n h 2 X h , the approximation of u n and d n+1,k h 2 S h , the piecewise approximations of d n+1,k .

The computation of n+1,k+1 h is once again explicit. The discrete version of (4.14) writes:

n+1,k+1 h = K 2 6 4 0 ⇣ d n+1,k h ⌘ ↵ 1 t n h 0 ⇣ d n+1,k h ⌘ ↵ 1 [⌧ n h + a (ru n h , n h )] + ⌘ 0 ⇣ d 0 h n+1,k ⌘ ↵ exp[ c ⇤ (1 A n h )]K : D(u n+1,k+1 h
) .

(4.17)

Subproblem P2

The weak form of the damage evolution equation is obtained by introducing a scalar test function ' h 2 S h .

The discontinuous Galerkin approximation of the advection term leads to the following discretized form of the damage evolution subproblem:

(P 2) h : Find d n+1,k+1 h 2 S h , 0 < d n+1,k+1 h  1, such that b h (d n+1,k+1 h , ' h ) = m n+1,k+1 h (' h ), 8' h 2 S h , with b h () and m h () given by b h (d n+1,k+1 h , ' h ) = Z ⌦ h 1 t d n+1,k+1 h ' h dx (4.18) m n+1,k+1 h (' h ) = Z ⌦ h " 1 t + min " 1, t n h 2 n+1,k+1 h , c n h 1 n+1,k+1 h q 2 n+1,k+1 h # 1 ! 1 t d # d n h ' h dx + Z ⌦ h 1 t h ' h dx X S2L (i) h Z S ✓ 1 2 |u n h • n| [[d n h ]][[' h ]] (u n h • n)[[d n h ]]{{' h }} ◆ ds Z @⌦ h max(0, u n h • n) (d n h d )' h ds (4.19)
where the last terms in both (4.18) and (4.19) handle the flux of ice through the upstream boundaries.

Subproblem P3

The transport equation for the cohesion is written in weak form using the test function ' h 2 S h . The discretized subproblem writes :

(P 3) h : With u n+1 h known, find C n+1 h 2 S h such that c h (C n+1 h , ' h ) = o n+1 h (' h ), 8' h 2 S h , with c h () and o h () given by c h (C n+1 h , ' h ) = Z ⌦ h 1 t C n+1 h ' h dx (4.20) o n+1 h (' h ) = Z ⌦ h 1 t C n h ' h dx X S2L (i) h Z S ✓ 1 2 u n+1 h • n [[C n h ]][[' h ]] (u n+1 h • n)[[C n h ]]{{' h }} ◆ ds Z @⌦ h max(0, u n+1 h • n) (C n h C )' h ds (4.21) Subproblem P4
The conservation equations for the ice thickness and concentration are also written in weak form using the scalar test function ' h 2 S h . Both problems read similarly as

(P 4.1) h : With u n+1 h known, find A n+1 h 2 S h such that c h (A n+1 h , ' h ) = p n+1 h (' h ), 8' h 2 S h , (P 4.2) h : With u n+1 h known, find h n+1 h 2 S h such that c h (h n+1 h , ' h ) = p n+1 h (' h ), 8' h 2 S h ,
with the bilinear form c h () and linear form p h () given by

c h (A n+1 h , ' h ) = Z ⌦ h 1 t A n+1 h ' h dx (4.22) p n+1 h (' h ) = Z ⌦ h 1 t A n h ' h dx X S2L (i) h Z S ✓ 1 2 |u n+1 h • n|[[A n h ]][[' h ]] (u n+1 h • n)[[A n h ]]{{' h }} ◆ ds Z @⌦ h max(0, u n+1 h • n) (A n h A )' h ds Z ⌦ h A n h (r • u n+1 h ) dx (4.23)
and similarly for h h . After solving these equations, mechanical redistribution is applied to both quantities such that :

h 0 h = max[0, (1 A n+1 h )]h n+1 h + h n+1 h A 0 h = A n+1 h max[0, (A n+1 h 1)]
and the fields of h n+1 h and A n+1 h are set respectively to h 0 h and A 0 h .

Combining all these definitions leads to the following scheme for the complete discretized channel flow problem, (P ) h :

(P ) h : Initialization (n = 0)

u n h = 0 ms 1 in ⌦ h n h = 0 Pa in ⌦ h , d n h = 1 in ⌦ h , C n h = C 0 (Pa) in ⌦ h , A n h = 1 in ⌦ h , h n h = 1 m in ⌦ h , For n 0 • For k = 0, ⇣ n+1,0 h , u n+1,0 h , d n+1,0 ⌘ = ( n h , u n h , d n h ),
• For k 0

-(P 1.1) h : With n+1,k h , u n+1,k h , d n+1,k h , h n h and A n h being known, find u n+1,k+1 h 2 V h such that a(u n+1,k+1 h , v h ) = l (n+1,k+1) h (v h ), 8v h 2 V h ,
with the bilinear and linear forms a() and l h () defined by (4.15) and (4.16) respectively and with

V h = {v h 2 X h , v h = 0 on left [ right } -(P 1.2) h : With u n+1,k+1 h and d n+1,k h known, compute n+1,k+1
h 2 T h explicitly from (4.17).

108

-(P 2) h : With n+1,k+1 h known, find d n+1,k+1 h 2 S h , 0 < d n+1,k+1 h  1, such that b h (d n+1,k+1 h , ' h ) = m n+1,k+1 h (' h ), 8' h 2 S h ,
with b h () and m h () given by (4.18) and (4.19).

-Compute

res h = 0 (d n+1,k+1 h ) ↵ 1 " n+1,k+1 h n h t + (⌧ n h + a (ru n h , n h ) # + d1 d n+1,k+1 h e n+1,k+1 h ⌘ 0 (d 0 h n+1,k+1 ) ↵ exp[ c ⇤ (1 A n h )]K : D(u n+1,k+1 h ) If res h < tol STOP . -Set n+1 h , u n+1 h , d n+1 h = ⇣ n+1,k+1 h , u n+1,k+1 h , d n+1,k+1 h ⌘ . • (P 3) h : With u n+1 h known, find C n+1 h 2 S h such that c h (C n+1 h , ' h ) = o n+1 h (' h ), 8' h 2 S h , with c h (C n+1 h
, ' h ) and o h () given by (4.20) and (4.21).

•

(P 4.1) h : Find A n+1 h 2 S h such that c h (A n+1 h , ' h ) = p n+1 h (' h ), 8' h 2 S h , (P 4.2) h : Find h n+1 h 2 S h such that c h (h n+1 h , ' h ) = p n+1 h (' h ), 8' h 2 S h , with c h (A n+1 h , ' h ), c h (h n+1 h , ' h )
given by (4.22) and p h () given by (4.23) respectively. Apply mechanical redistribution, such that 0  A n+1 h  1 and 0  h n+1 h .

Results

Dynamical behavior

We first investigate the formation of arches and stable ice bridges by analyzing the evolution of damage and ice velocity over the channel. As a determining factor for the formation of free arches is the cohesive nature of the material, here simulations with a di↵erent range of values of C are compared, with the lowest value, Apart for the timing and level of localization of the damage, the simulated dynamics is similar in between the simulations (this is also true for the case of C min = 30000 Pa, not shown).

Near the onset of damaging, indicated by dot number 0 on the wind forcing curves, panels a, the (downstream) ice velocity at the centre of the channel is on the order of 10 5 ms 1 (not shown), consistent with small and strictly elastic deformations. In all cases, the first peak in the damage rate corresponds to the propagation of damage in the interior and downstream of the channel. Damage is observed to concentrate along curved and concave features (see figures 4.3b to 4.6b, panel 1), which appear more localized as the value of C min is increased. After this first damaging event, the flow of ice upstream of the channel stops in all simulations (see figures 4.3c to 4.6c, panel 1). The profile of u y indicates a jump in the ice velocity that coincides with the location of a concave arch, meeting the channel boundaries at the constriction point. Downstream of this arch, u y varies between 10 2 to 10 1 ms 1 and is hence several orders of magnitude larger than within the initially undamaged ice cover, consistent with an important contribution of viscous stress relaxation into large deformations. In all cases, the profile of u y is "stair-case"-like, indicating some regions of piecewise rigid motion separated by strong discontinuities. Comparison of fields of the magnitude of u and of the level of damage (see figure 4.7, panel 1 and figure 4.6b, panel 1) indeed indicate that narrow, highly damaged areas divide the ice into relatively undamaged plates, themselves moving like solid bodies, consistent with the motion of the Arctic sea ice cover as revealed by Synthetic Aperture Radar imagery analysis and RGPS motion products [START_REF] Kwok | Deformation of the Arctic Ocean sea ice cover: November 1996 through[END_REF][START_REF] Moritz | Relationship between geostrophic winds, ice strain rates and the piecewise rigid motions of pack ice[END_REF], see figure 1.2, chapter 1).

Damaging resumes as the wind forcing is increased beyond this point. First, it slowly propagates upstream of the channel along closed, concave arches, along with the no-flow transition (see figures 4.3b, c and 4.4b, c, panel 2). We observe that closed arches form somewhat further northward of the constriction point as the range of values of C is increased. Damage then propagates rapidly along open arches in the interior and along nearly strait features towards the exterior of the domain, leaving low damage areas in the centre and near the domain boundaries in the converging part of the basin (see figures 4.3b to 4.6b, panel 3). These features define a converging channel within the basin, with some regions of fast ice attached to the domain boundaries (see figure 4.7, panel 3). Similar results were also obtained by [START_REF] Dumont | Modeling the dynamics of the north water polynya ice bridge[END_REF]. At the centre of this inner channel, the southward flow of ice is almost uniform (see figures 4.3c to 4.6c, panel 3).

From that point on, the dynamics di↵ers between the simulations. In the case of C min = 2000 Pa, the flow of ice keeps increasing as the wind stress is further intensified. Inflowing ice gets damaged as it enters the basin. In the case of C min = 5000 Pa, the flow of ice both upstream and within the channel subsequently slows down and almost, but not quite, stops (see 4.4 c, panel 4). At that point the damage rate time series shows a third, smaller peak indicating the partial collapse of an ice bridge within the channel, after which the ice velocity re-increases (see 4.4 c, panel 5). In the simulations with C min 10000 Pa, the flow of ice upstream and within the channel decreases and e↵ectively stops (see figures 4.5c and 4.6d, panels 4). The ice over the basin and channel remains motionless while u y subsequently increases with the wind forcing downstream of the channel, suggesting the formation of a stable ice bridge near the opening of the channel.

We further investigate the evolution of the internal stress in the simulations with C min 10000 Pa. the principal (i.e., the maximum and minimum normal) stresses, defined as

1 = ( 11 + 22 ) 2 + s  11 22 2 2 + 2 12 , 2 = ( 11 + 22 ) 2 s  11 22 2 2 + 2 12 ,
with the convention that compressive stresses are positive. The principal stresses within the initially undamaged ice cover (panels 0) are compressive ( 1 , 2 > 0) over the basin, change sign at the middle of the channel (y = 0) and are tensile ( 1 , 2 < 0) downstream of the channel (not shown). After the first damaging event, 2 becomes negative (tensile) over the converging part of the basin, with 1 > 0 (compressive), hence indicating a stress state dominated by shear (panels 1). A maximum in 1 collocated with a minimum in 2 coincides with the position of the ice arch at the mouth of the channel and with the no-flow transition, suggesting this arch was e↵ectively formed by Coulomb failure.

Downstream of this arch, other local maxima in

1 are collocated with local minima in 2 , consistent with the presence of multiple arch-like features. The subsequent stress profiles (panels 2) suggest a similar failure mechanism for the arches forming upstream of the channel. After the second damaging event, both stress components show large fluctuations over the basin, with stress states on average compressive for y > 0 and tensile for y < 0. Panels 4 corresponds to the stopping of the flow over the the basin and formation of the stable bridge near the opening point of the channel. In the case of C min = 10000 Pa, a clear maximum in 1 , a state of biaxial tension. This is an important point, as most viscous-plastic sea ice models, based on an elliptical yield curve suggested by [START_REF] Hibler | A viscous sea ice law as a stochastic average of plasticity[END_REF], do not account for pure (uniaxial or biaxial) tensile strength. This suggests that these models would not be able to reproduce a stable ice bridge and self-obstruction to flow under the conditions simulated here.

As the wind stress is further increased, both principal stresses goes to zero downstream of the ice bridge, consistent with the formation of a stable stress-free surface. or

u y = s ⇢ a C da ⇢ w C dw u a ⇡ 2% u a , (4.25)
with u a ⇡ 22 ms 1 in this case.

We further comment on the nature, i.e., tensile, compressive, biaxial tensile-compressive (i.e. sheardominated), of the internal stress within the simulated ice cover and compare the simulated stress states to the in-situ stress measurements over the Beaufort Sea, Arctic Ocean, reported by [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF] and [START_REF] Weiss | Coulombic faulting from the grain scale to the geophysical scale: lessons from ice[END_REF]. Instantaneous stress states are plotted for each model element after the second damaging event (figure 4.8b, panel number 3), that is, before the formation of a stable ice bridge, and after the formation of the bridge and stoppage of the ice flow (panel number 5). In agreement with stress measurements over the open Arctic Ocean, biaxial tensile stresses occur in a large number before the formation of the ice bridge and large biaxial compressive stresses are not so frequent compared to pure tensile and biaxial tensile-compressive states [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF]. After the formation of the bridge, biaxial compressive stresses become more frequent, consistent with the convergence of the flow over the basin and within the channel. Tensile stresses on the other hand are less frequent, but still present in a significant number. This again supports the relevance of accounting for some resistance of the ice in pure tension in sea ice models. It is also important to note that isotropic stresses ( 1 = 2 ) are frequent in the observations and absent in the model. Theses are associated with thermal processes, i.e., thermal expansion and fracturing [START_REF] Richter-Menge | Relating arctic pack ice stress and deformation under winter conditions[END_REF], which are not represented in the present Maxwell-EB framework.

Ice thickness and concentration

In this section we briefly describe the evolution of the fields of ice concentration and thickness simulated by the The presence of multiple arch-like features within and downstream of the channel appears clearly on the fields of ice concentration (see figure 4.9, panel number 3). The value of A rapidly drops to  0.9 at the opening point of the channel and along the domain boundaries downstream of the channel. As our functional dependance of the mechanical parameters on the ice concentration, based on the dependance of the ice strength in compression (P ) in the VP model of [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF] (see equations 4.6 and 4.7), allows for the value of E and ⌘ to drop by almost 90% for A = 0.9 and c ⇤ = 20, this is consistent with the detachment of the ice at the opening point of the channel, i.e., where the velocity field becomes divergent. Observational studies of sea ice in narrow straits however have emphasized a preferential location of stable ice bridges near constriction points (e.g. [START_REF] Kwok | Large sea ice outflow into the nares strait in 2007[END_REF]. As in most channel geometries explored in the present simulations (not shown), an ice arch does form near the constriction point but ice preferentially detaches along an arch located at the opening point of the channel, this might indicate that a lower value of c ⇤ would be more appropriate in the Maxwell-EB model, in which case a sensitivity study of the ice bridge position on the value of c ⇤ would be needed. Also evident from these fields is the di↵usivity of the numerical scheme. While the ice edge is clearly defined at the location of the bridge, the edge of the detached, free drifting ice is much less clearly defined. After a short time, gradients in ice concentration associated with the presence of arch-like features downstream of the bridge are also significantly smoothed due to di↵usion. Strong di↵usion also appears obvious from the fields of the level of damage (see figures 4.3b to 4.6b). It is most important within the free-drifting ice. Upstream of the ice bridge, damaging tends to refine and re-increase gradients in all fields. It is also strongly asymmetric, consistent with thickness distributions estimated for first and second-year ice at the end of the growth season, i.e., for an ice cover with little history of melting (e.g., [START_REF] Haas | Dynamics Versus Thermodynamics: The Sea Ice Thickness Distribution[END_REF]. The tail of the distribution (h > 1 m), appears to fit a negative exponential function. The associated slope decreases with time as the ice in the converging part of the domain thickens. This statistical property has been observed for the tail of probability density functions (PDF) calculated for in-situ sea ice thickness measurements [START_REF] Wadhams | Sea Ice Thickness Changes and their Relation to Climate[END_REF][START_REF] Haas | Dynamics Versus Thermodynamics: The Sea Ice Thickness Distribution[END_REF], indicating the strong spatial localization of the thick, ridged ice.

In the Maxwell-EB model, this property seems to arise naturally from the capacity of reproducing the strong spatial localization of sea ice deformation.

The part of the frequency histogram for h < 0.8 m corresponds to the ice at the opening and downstream of the channel. As ice drifts out of the domain, a second mode appears and progressively moves towards h = 0, consistent with the increasing presence of open water. The distribution of the ice concentration, figure 4.11a, shows a similar behaviour. This mode however, is not nearly as important as that normally obtained for PDF of sea ice thickness measurements in the presence of open water (e.g., [START_REF] Wadhams | Sea-ice topography of the arctic ocean in the region 70 degrees w to 25 degrees e[END_REF][START_REF] Wadhams | Sea Ice Thickness Changes and their Relation to Climate[END_REF][START_REF] Haas | Comparison of the sea-ice thickness distribution in the lincoln sea and adjacent arctic ocean in 2004 and 2005[END_REF][START_REF] Haas | Dynamics Versus Thermodynamics: The Sea Ice Thickness Distribution[END_REF]. The fact that this mode is small here and the frequency histogram (approximating a PDF here) almost flat between h = 0 and h ⇡ 0.8 m is mostly attributable to the numerical di↵usion of gradients within and especially at the edge of the free drifting ice.

Concluding remarks

In this chapter, the Maxwell-EB rheological framework was used for modelling the dynamics of sea ice on geophysical scales. The simulations described here represent the fracturing and drift of sea ice within a narrow strait under the action of orographically channeled winds. The simulations have shown that the model is able to reproduce:

• The formation of concave arch-like faults within, downstream and, when increasing the value of C, further upstream of the channel. As opposed to previous viscous-plastic model simulations, both idealized and realistic, of the flow through Nares strait and in agreement with observations, the Maxwell-EB model reproduces a high density of arch-like faults within the ice cover, with the associated narrow zones of low ice concentration.

• The formation of a stable ice bridge and the stoppage of the flow of ice upstream of the channel, provided that the ice has su cient cohesion, and the opening of the basin downstream of the bridge. It is important to note that the simulations presented here are by no means an attempt to determine a physically appropriate range of values of the cohesion and of the critical tensile strength for the sea ice cover in the Maxwell-EB model, as the presence of a stable ice bridge is expected to depend on the applied wind forcing, on the domain geometry, the channel width, etc. Adequate values of C at the scale of the model element should rather be determined from in-situ stress measurements within sea ice and from the appropriate scaling laws [START_REF] Weiss | Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction[END_REF][START_REF] Schulson | Compressive shear faults within arctic sea ice: Fractures on scales large and small[END_REF].

• Discontinuities in the ice velocity fields between regions of uniform, plate-like motion, consistent with RGPS observations (see figure 1.2, left panels). In particular, the model reproduces the very sharp gradients in ice velocity, concentration and thickness associated with the edge of the bridge.

• States of stress that are overall in good agreement with in-situ measurements in Arctic sea ice.

• The formation of ice ridges, with a strong localization of the thickest ice and an associated thickness frequency distribution that fits a negative exponential, in agreement with probability density functions calculated from thickness measurement within sea ice. Due to the high di↵usion of the numerical scheme, the mode associated with the presence of open water is underrepresented. This last point is an important outcome of this first implementation of the Maxwell-EB rheological framework in a realistic context, as the adequate representation of the mechanical redistribution of the ice thickness, in particular the formation of ice ridges, has been the subject of numerous studies since the 1970's and continues to challenge the sea ice modelling community to this day. Multi-thickness categories schemes based on the pioneer work of [START_REF] Thorndike | The thickness distribution of sea ice[END_REF] and [START_REF] Rothrock | The energetics of the plastic deformation of pack ice by ridging[END_REF] are widely used in models to allow accounting for the presence of di↵erent ice thicknesses at the subgrid scale and determining the thickness of the sea ice involved in the ridging process. Such schemes handle an evolution equation for the thickness distribution over each model grid cell and use redistribution functions based on the simulated strain rates. As the strain rate tensor itself does not directly provide the information on the relative amount of opening and ridging (and also sliding) within the ice cover, expressions for these modes of redistribution need to be prescribed in such models, the form of which can be set in an ad-hoc manner [START_REF] Thorndike | The thickness distribution of sea ice[END_REF], estimated empirically from strain rates observations [START_REF] Stern | Open water production in arctic sea ice: Satellite measurements and model parameterizations[END_REF] or determined based on of the prescribed form of the yield criteria and flow rule in the case of plastic models [START_REF] Rothrock | The energetics of the plastic deformation of pack ice by ridging[END_REF].

Although the excessively simple redistribution function used here does not include multiple ice categories nor allow prescribing the thickness of the ice involved in the ridging, the Maxwell-EB model seems to be able to reproduce the observed distribution of extreme thickness values corresponding to ridged ice because of its adequate representation of the mechanical processes and, in particular, of the strong spatial localization of strain rates within the sea ice cover. Of course, larger scale and longer term simulations would be required in order to investigate more extensively the behaviour of the Maxwell-EB model in terms of the mechanical redistribution of sea ice. A thorough comparison of simulated and measured thickness distributions also calls for the representation of thermodynamic processes in the model, which contrary to mechanical redistribution, are expected to "seek the mean" [START_REF] Thorndike | The thickness distribution of sea ice[END_REF] and smooth the distribution by allowing ice to grow over open water in the winter, thinner ice to thicken faster than thick ice and thick ice to melt in the summer [START_REF] Haas | Dynamics Versus Thermodynamics: The Sea Ice Thickness Distribution[END_REF].

Chapter 5

Conclusions 5.1 A brief summary

This thesis presents a new continuum mechanics framework suited for modelling the brittle behaviour and deformation of the sea ice cover that is built from the following ingredients:

1. a constitutive law combining linear elasticity and a dissipative mechanism for the stress, 2. a threshold mechanism for brittle failure, in the form of a Mohr-Coulomb and tensile damage criterion, 3. disorder in the local damage criterion, representing the material's natural heterogeneity and expressing the spatial variability of its mechanical strength, 4. an isotropic progressive damage mechanism allowing for the elastic redistribution of the stress from overto sub-critical areas of the material and involving a non-dimensional scalar damage parameter that reflects the integrated history of the damaging, 5. a recovery, or healing, mechanism for the material's mechanical properties.

The leading novelty of this mechanical framework, named Maxwell-Elasto-Brittle, with respect to the original Elasto-Brittle [START_REF] Amitrano | From di↵use to localised damage through elastic interaction[END_REF]Girard et al., 2010a) and to the recent EB-based sea ice model NeXtSIM [START_REF] Bouillon | Presentation of the dynamical core of neXtSIM a new sea ice model[END_REF][START_REF] Rampal | neXtSIM: a new Lagrangian sea ice model[END_REF] lies in the formulation of the constitutive law, that is, in the inclusion of a viscous-like relaxation term for the elastic stresses. Like the elastic modulus in the EB model, the associated apparent viscosity was coupled to the local level of damage, with the idea of letting internal stresses dissipate into large, permanent deformations where the ice cover is highly fractured while allowing for relatively undamaged areas of the ice to retain the memory of the stresses associated with elastic deformations.

The following paragraphs summarize the main conclusions drawn from the analysis of the Maxwell-EB simulations presented in this manuscript.

A first series of conclusions relates to the mechanical behaviour reproduced by the Maxwell-EB model.

The set of small-deformations simulations presented in chapter 2 has allowed demonstrating that the model represents:

• the strong localization of deformation in both space and time, indicative of its heterogeneity and intermittency, and the associated scaling laws,

• the anisotropy of the deformation due to the elastic kernel describing the Coulomb stress redistribution at the onset of damaging,

• a competition between healing and damaging within the material, which translates into two main types of behaviour : (1) a succession of rapid damaging events with the associated stress relaxation and slower recovery periods, characterized by the activation of new faults with di↵erent shapes and orientations and

(2) creep-like deformation and the persistence of activated faults,

• the emergence of a marginally stable state in which the material remains near criticality, with long range temporal correlations in the damage rate that translate into an invariance of scales at all time scales below the characteristic time of healing of the material.

The last point in particular is the consequence of allowing for the relaxation of the internal stresses, at a rate that depends on the local level of damage. A sensitivity analysis on the rate of decrease of the relaxation time with increasing level of damage, set by the damage parameter ↵, has revealed that a too brutal relaxation of the internal stress at the onset of damaging e↵ectively results in the loss of temporal correlations in the damage rate. In this case the model becomes essentially elasto-plastic. The deformation is still highly localized in space, although the emergence of spatial scale invariance remains to be established, but exhibits no intermittency.

Chapter 3 described the transition from a small-deformation to a large-deformation Maxwell-EB model and, in particular, the treatment of the advection processes and of the rotation and deformation of the internal stress tensor terms in the constitutive equation. A numerical scheme for the large-deformation dynamical model was presented, and used in the context of simulating a laboratory Couette experiment performed on a thin sheet of ice. The main outcomes of this chapter are:

• As the model reproduces the very strong gradients observed within sea ice, our current discontinuous Galerkin approximation of advection terms using finite element approximations of order 0 is expected to lead to significant numerical di↵usion within the advected fields.

• The Couette flow simulations allowed identifying some numerical limitations of the model as well as some important requirements on the relative value of the model parameters. An important numerical limitation pertains to the treatment of the progressive damage mechanism and the implicit requirement that the time step be comparable to that of the propagation of damage, t d , within the simulated material.

In the context of ice (fresh ice or sea ice) at the laboratory scale, this requirement results in unpractical simulation times. Hence a slower propagation of damage was assumed. This assumption translated into an insu cient separation of scales between the rate of damage propagation and the rate of viscous relaxation of the internal stress within the material, itself limited by the numerics, i.e., by the "large Weissenberg number problem". Consequently the simulated deformation was too viscous-like. One important requirement for the adequate representation of the brittle character of the deformation in the Maxwell-EB model hence consists in that the relaxation time for the stress over low damage areas of the ice be orders of magnitude higher than the characteristic time for the propagation of damage.

• Another important requirement for the proper representation of the competition between fracturing and refreezing within the ice cover when large deformations occur is that the relaxation time for the internal stress over low damage areas be larger than the characteristic time for healing. This condition ensures that healing is e cient in allowing damaged areas to recover and the stress the material to significantly build up, so that to trigger subsequent extensive damaging events. 123 Finally, chapter 3 presented a first implementation of the Maxwell-EB sea ice model in a realistic context.

The main conclusions drawn from the channel flow simulations are:

• The Maxwell-EB rheological framework simulates the propagation and localization of damage along narrow, concave arch-like features and the formation of a stable ice bridge within the channel with the associated stoppage of the ice flow.

• The model represents the formation of ice "plates", or floes, and their piecewise rigid motion, with clear discontinuities in the velocity field between the floes, consistent with observations.

• The simulated states of stress are in good agreement with measurements within the Arctic ice pack.

Uniaxial and biaxial (i.e., "pure") tensile stresses were shown to be frequent, hence emphasizing the importance of accounting for some resistance in pure tension in sea ice models.

• The model reproduces the formation of ice ridges, with the associated thickness distribution fitting a negative exponential, also in agreement with observations. Numerical di↵usion however leads to an underrepresentation of the mode corresponding to open water.

Another point worth mentioning is that in these realistic simulations, setting the time step to the characteristic time for damage, t d , estimated from the speed of propagation of shear elastic waves within the ice cover and from the model spatial resolution (2 km), lead to reasonable simulation times. The relaxation time for the internal stress, , was set to 10 7 seconds over undamaged ice, allowing the convergence of the numerical scheme.

The large separation between t (4 seconds) and 0 lead to a proper representation of the brittle character of the deformation. This value of the undamaged relaxation time corresponds to about 115 days (almost four months), and hence is comparable to the duration of the polar winter (about 5 months). Therefore, unphysical viscous dissipation is expected to be negligible in these realistic simulations of the Arctic sea ice cover.

5.2 Some perspectives 5.2.1 Towards a fully continuous formulation of the Maxwell-EB rheology Some 95 pages ago, in chapter 2, a continuous expression was presented for the damage term in the damage evolution equation (see equation 2.20). At that point, it was mentioned that this formulation relies on the specific use of a backward explicit time discretization scheme of order 1. The progressive damage mechanism in the Maxwell-EB framework is indeed an intrinsically discrete process, tied to the local value of the damage criterion and to the spatial discretization of the model. The current representation of this mechanism constitutes a limiting factor in the numerical development of the model in many respects.

• First, the precise representation of this mechanism implies setting the model time step to the characteristic time for damage, t d . In some applications of the model, in particular when simulating the deformation of ice at the laboratory scale, this leads to unpractical simulation times. This issue can be handled in two ways. First, the model time step can be increased with respect to t d . However for time steps significantly larger than t d , elastic interactions within the simulated material are hindered and numerical instabilities arise. Second, for the sake of computational e ciency the time of propagation of damage in the material can be assumed longer. However, this assumption might result in an inadequate representation of the brittle character of the deformation. Luckily, in some sense, for regional and global sea ice model spatial resolutions, t d , and hence t, falls between a couple of seconds and a couple of minutes, leading to reasonable simulation times. However, the dependance of the model time step on the spatial resolution and speed of propagation of damage within the simulated material still represents a limitation of the current Maxwell-EB framework as it restricts its range of applications.

• Second, the dependance of t d (hence of t) on the spatial resolution of the model precludes the use of multi-resolution mesh grids, otherwise one of the main advantage of the finite element approach. On the one hand, setting t d based on the largest element of a multi-resolution grid would indeed result in a loss in the resolution of the damage propagation over higher resolution areas of the domain. On the other hand, setting t d based on the smallest model element would not be physically rigorous as it implies a speed of propagation of damage that is locally (i.e., over the largest elements) higher the speed of elastic waves in the material.

It is also important to note that the introduction of disorder in the damage criterion at the smallest resolved scale, combined to the tendency of the model to localize deformation and damaging at the scale of the mesh element, whatever the element size, also precludes the use of adaptive remeshing techniques with the Maxwell-EB model. Increasing spatial resolution where gradients are steep would e↵ectively lead to finer features, but would not help "capturing" the solution. Instead, the solution would be ever-evolving towards an increased localization. In a model based on a finite correlation length for the damaging and allowing the convergence of the solution at a finite scale, adaptive remeshing would otherwise be easily implemented in our discontinuous Galerkin numerical scheme.

• Finally, our current treatment of the damage evolution, based on the direct calculation of the distance of the local state of stress to the local damage criterion, is consistent with finite element approximations of order 0 for the fields of cohesion, internal stress and level of damage.

These numerical limitations prompts for a truly continuous mathematical representation of the damage mechanism in the Maxwell-EB model. This would allow the use of implicit, higher-order time-stepping schemes, of a larger model time step, of multi-resolution mesh grids and of higher order finite element approximations for the internal stress and damage fields, which as demonstrated in section 3.2, would significantly help reducing the numerical di↵usion stemming from advection.

It is also important to note that the use of low order finite element approximations is not restricted to the present mechanical framework. To our knowledge, in no Eulerian finite element scheme for sea ice modelling the degree of polynomial approximations have been increased beyond 1 for the velocity (e.g., [START_REF] Lietaer | The e↵ects of resolving the Canadian Arctic Archipelago in a finite element sea ice model[END_REF][START_REF] Terwisscha Van Scheltinga | A finite element sea ice model of the canadian arctic archipelago[END_REF][START_REF] Danilov | Finite-element sea ice model (FESIM), version 2[END_REF] and the ice thickness and concentration (e.g., [START_REF] Terwisscha Van Scheltinga | A finite element sea ice model of the canadian arctic archipelago[END_REF][START_REF] Danilov | Finite-element sea ice model (FESIM), version 2[END_REF] and beyond 0 for the internal stress and ice strength (i.e., the pressure P in viscous-plastic type models, e.g., [START_REF] Lietaer | The e↵ects of resolving the Canadian Arctic Archipelago in a finite element sea ice model[END_REF][START_REF] Danilov | Finite-element sea ice model (FESIM), version 2[END_REF]. [START_REF] Lietaer | The e↵ects of resolving the Canadian Arctic Archipelago in a finite element sea ice model[END_REF], who have used the two-levels thickness model of [START_REF] Hibler | A dynamic thermodynamic sea ice model[END_REF]) and a mechanical redistribution scheme consistent with a constant-by-element thickness and concentration similar to the one employed here, have also used order 0 FE approximations for the ice thickness and concentration. Increasing the order of the finite element approximations would probably result in important improvements of model performances. Discontinuous Galerkin methods would represent a compact, numerically e cient and flexible way to do so.

On the physical point of view however, a critical question remains unanswered: are we dealing with an inherently discrete problem? In other words, can we represent the mechanical behaviour simulated with the Maxwell-EB framework with a truly continuous model? The same question was asked in the context of modelling the Earth crust [START_REF] Main | Statistical physics, seismogenesis, and seismic hazard[END_REF]. This point might constitute the main challenge in the future development of EB-type sea ice models. On the one hand, the classical viscous-plastic framework might be numerically more convenient, especially in the context of a coupling with an atmosphere and ocean component, but it raises issues about the representation of the brittle character of the ice deformation. The present Maxwell-EB framework on the other hand, reproduces the brittle character of the deformation, with the observed strong spatial gradients within the ice velocity, strain and stress fields, but this performance comes at the cost of numerical limitations.

Future developments and additions to the rheology

Besides numerical improvements, we can also think of relevant additions to the Maxwell-EB rheological framework itself. Here some of these possible enhancements are briefly discussed.

• Coupling Poisson's ratio to the level of damage:

In the present implementation of the Maxwell-EB rheology, Poisson's ratio, ⌫, is taken constant. As mentioned in section 2.2.4, rock mechanics experiments have shown that ⌫ e↵ectively increases within a fractured material consistent with its dilatation [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF][START_REF] Martin | The progressive fracture of Lac du Bonnet granite[END_REF][START_REF] Heap | The evolution of elastic moduli with increasing crack damage during cyclic stressing of a basalt from mt. etna volcano[END_REF]. This additional level of complexity could be easily accounted for in the Maxwell-EB model by allowing Poisson's ratio to evolve with the level of damage, equivalent to letting the sti↵ness tensor K be a function of d. In this case, an appropriate form for the coupling between ⌫ and d would need to be determined and tested through sensitivity analyses. In particular, the dependance of ⌫ on d is expected to impact the orientation of the simulated sets of faults.

• Coupling the cohesion and level of damage:

Another possible refinement to the rheology would consist in letting the cohesion of the material vary as a function of the level of damage, consistent with a local decrease of the failure threshold over weakened elements. In general, this is expected to result in more abrupt fracturing events as elements are susceptible to get damaged over successive time steps. By the same fact, it can result in a more rapid loss of the memory of elastic stresses within the fault. The form of this additional coupling, as well as its impacts on the spatial and temporal localization of the deformation and associated scaling laws needs to be determined.

• A threshold for permanent deformations:

In the Maxwell-EB framework, a damage threshold is introduced in the form of a Mohr-Coulomb and tensile criterion to represent terminal failure within a brittle material. The model does not include a deformation threshold. Besides mechanical strength recovery through healing, no mechanism represents impediments to permanent deformations nor "true" friction along faults once formed. Such mechanisms could be roughly accounted for by prescribing a minimum value of the internal stress for permanent deformations. This is represented schematically on figure 5.1a by the sliding friction element, connected in parallel to the dashpot and in series with the spring. One possible formulation for this threshold is represented on figure 5.1b by the branch of equation 1 q 2 = 0.

(5.1) 126

Combining both this threshold and the damage criterion, the simulated deformation would be partionned as follow:

-The deformation of undamaged elements would be strictly elastic, -The deformation of damaged elements with 1 > q 2 would be split between small, elastic and permanent, potentially large, deformations, -The deformation of damaged elements with 1 < q 2 would again be strictly elastic. Finally, the next logical step in the development of the Maxwell-EB sea ice model consists in a complete validation of the mechanical framework by comparison of the simulated deformation rates and associated scaling properties to that estimated from the available ice buoy and RGPS Arctic ice motion data. These analyses will necessitate carrying numerical experiments with both realistic wind forcing and ocean currents over periods of days to months and over domains of regional to global scales. Over such long time scales, thermodynamics e↵ects will also need to be accounted for. Parameterizations of heat fluxes between the atmosphere, the ice cover and the ocean, as well as a refined formulation of the healing rate in terms of the local di↵erence between the surface air temperature and underlying water temperature will need to be included. Conversely, the representation of thermodynamic processes in the model will open multiple perspectives of analyses. For instance, it will allow assessing the impact of the Maxwell-EB rheology on the simulated thickness distribution of sea ice, as well as on the simulated heat fluxes and ice melt and growth rates of over the Arctic.

Figure 1 . 1 :

 11 Figure 1.1: Areal pictures of sea ice (from bottom to top) at the scale of 100 meters, ⇠ 5 km and 60 km. At the smallest time and space scales, the discontinuous nature of the ice cover cannot be ignored when modelling its deformation and drift. At large scales ( 10 4 m), the ice cover is constituted by a large number of individual ice floes of di↵erent shapes and sizes and is described by mean quantities in continuum models.

Figure 1 . 2 :

 12 Figure 1.2: Weekly fields of divergence (upper left), vorticity (middle left), shear (lower left) and sea ice motion (x-component, y-component and magnitude) for the period of February 8 to February 14 2008 over the Arctic basin, showing the concentration of deformation within linear features and clear discontinuities in the velocity field. The fields are estimated from the 3-days RADARSAT Geophysical Processor System sea ice motion products (http://rkwok.jpl.nasa.gov/graphics/index.html)

"

  Figure 1.3: (a) Schematic representation of the Elasto-Brittle model, which combines a linear elastic constitutive law

  The final chapter restates the novelties of the Maxwell-EB framework with respect to the current rheological approaches for modelling the drift and deformation of sea ice. Possible solutions to the current limitations of the Maxwell-EB model, both physical and numerical, are presented and future additions to the framework are suggested. ., Weiss, J., Saramito, P. and Lattes, P. (2015), A Maxwell-Elasto-Brittle rheology for sea ice modelling, The Cryosphere Discussion, and with results from Weiss, J. and Dansereau, V. (2015), Linking scales in sea ice mechanics, submitted to Philosophical Transactions of the Royal Society A. The main objective of this chapter is to introduce the theoretical bases for a viscous-elasto-brittle model for sea ice and to describe in details the Maxwell-EB rheological framework. Section 2.1 relates and contrasts the original Maxwell and the Maxwell-EB constitutive laws. Section 2.2 presents the remaining ingredients of the Maxwell-EB model : the damage criterion, the progressive damage and healing mechanisms and the coupling of both mechanisms with the material's mechanical properties. Section 2.2.4 discusses the important non-dimensional numbers involved in the rheological framework as well as their absolute and/or relative values in the context of sea ice modelling. The setup of a first set of highly idealized model simulations in which

Figure 2 . 1 :

 21 Figure 2.1: Schematic 1-dimensional representations of the linear elastic model (a) and of the Maxwell model (b) for

Figure 2 . 3 :

 23 Figure 2.3: 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 would close 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. The calculation of the distance to the damage criterion dcrit, defined by the intersection ( 0 1 , 0 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.

  Figure 2.4: Schematic representation of the propagation of damage in the Maxwell-EB model. The domain is discretized

Figure 2 . 5 :

 25 Figure 2.5: 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.

d

  Figure 2.6: Elastic modulus (dashed lines) and apparent viscosity (plain lines) as a function of the level of damage, d, as described in equations (2.24) and (2.23), for di↵erent values of the constant ↵. The right panel zooms into the low ranges of values of d.

Figure 2

 2 Figure 2.7: (a) Domain and boundary conditions for the uniaxial compression experiment. (b) Example of finite element mesh with N = 10.

  us consider the time t T = t 2 [0, 1[. Let us further denote the velocity ũ = (ũ x , ũy ), the strain rate tensor "(ũ) = D(ũ) = 1 2 ( rũ + rũ T ) and the stress mixed: compression is applied by prescribing the Dirichlet conditions ũy = 1 on T and ũy = 0 on B . No confinement is applied on the lateral sides, hence the Neumann condition ˜ • n = 0 with n the unit normal vector on L and R . The top and bottom edges of the plate are allowed to deform in the x direction, implying ˜ xy = 0 on T and B . At the lower left corner of the plate ( B \ B ), ũ = 0.

Figure 2

 2 Figure 2.8: (a) Macroscopic stress (dashed-dotted line) and damage rate (solid grey line) as a function of time for a uniaxial compression experiment with N = 10. (b) Residual of the constitutive equation (res h ) as a function of the number of fixed-point sub-iterations k for a few model time steps, equally spaced along the interval of simulation delimited by the colorbar in (a).

Figure 2 .

 2 Figure 2.9 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, 1 d crit ) for the simulation with N = 40. Inspection of the initial loading and damaging sequence suggests that the mechanical behaviour

Figure 2

 2 Figure 2.9: (a) Macroscopic stress versus macroscopic strain (solid lines) for four uniaxial compression simulations with di↵erent 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 di↵erent times indicates on figure (a) and for the four simulations (resolution increasing from top to bottom).

Figure 2

 2 Figure 2.10: (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.

  Figure 2.11(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 2

 2 Figure 2.11: (a) Macroscopic stress (black dashed-dotted line) and damage rate (solid grey line) as a function of time for one realization of the uniaxial compression experiment with N = 40. The dashed red box indicates an interval of uninterrupted damaging activity, during which deformation is accommodated by a persisting system of interacting faults. (b) Instantaneous fields of level of damage at the five times indicated by blue dots on the macroscopic stress curve, showing the formation of the system of faults (first panel), which remains active for some time (three following panels), until the propagation of a new, non-interacting fault (last panel).

Figure 2

 2 Figure 2.12: (a) Average power spectral density of the damage rate time series for 5 realizations of the uniaxial compression experiment initialized with di↵erent fields of C and with N = 40. Blue dashed lines indicate, from left to right, the frequency associated with the characteristic time for healing, the inverse time of propagation of damage across the width of the domain and 1 2 ⇥ the frequency associated with the characteristic time for damage. The red dashed line indicates the frequency of the healing and damaging cycle marked with an arrow on figure 2.11(a). (b) Total deformation rate < "tot >t as a function of the observation time t, for 20 realizations of the coarse graining calculation centred on di↵erent arbitrary times t0 along a uniaxial compression experiment with N = 40 (coloured lines) and average of the 20 realizations (thick black line).

Figure 2

 2 Figure 2.13: Time series of (a) the macroscopic stress and (b) the damage rate in uniaxial compression experiments with ↵ = 0 (solid curves) and ↵ = 1 with We 0 = 5 • 10 2 (dotted line), We 0 = 5 • 10 4 (dashed-dotted line) and We 0 = 5 • 10 6 (dashed line).

Figure 2 . 14 :

 214 Figure 2.14: Time series of (a) the macroscopic stress and (b) damage rate in uniaxial compression experiments with ↵ > 1. The initial elastic loading is identical between all simulations, as all other mechanical parameter values are the same, and time series hence start at the onset of damaging. (c) Zoom in on the damage rate time series.

Figure 2 . 15 :

 215 Figure 2.15: Instantaneous fields of the order of magnitude of the total deformation rate (log 10 ( "tot)) normalized by the maximum value of "tot (upper panels) and of the level of damage (lower panels) at the instant indicated by the dot on figures 2.13a and 2.14a and for ↵ = 0, ↵ = 1 (We 0 = 5 • 10 4 and We 0 = 5 • 10 2 ), ↵ = 3 and ↵ = 10.

Figure 2 Figure 2 . 17 :{3

 2217 Figure 2.16: Principal states of stress at the instant indicated by the dot on figures 2.13a and 2.14a for the uniaxial compression experiments with (a) ↵ = 1 and We 0 = 5 • 10 2 , ↵ = 3 and ↵ = 10. The case of ↵ = 0 is identical to the case of ↵ = 1 and We 0 = 5 • 10 2 . The solid black lines represent the damage envelope corresponding to the minimum and maximum values of the cohesion C over the domain.

Figure 3

 3 Figure 3.1: (a) Domain and boundary conditions for the solid body rotation problem. (b) Example of low resolution mesh with (adimensional) mean element size of x = 1 10 .

Figure 3

 3 Figure 3.2: (Upper panels) Fields of (a) 1, (b) 2 and (c) 3 at t = 0 interpolated onto T

Figure 3 Figure 3

 33 Figure3.3: L 2 norm of the error between the simulated and exact (interpolated) solutions of the rigid body rotation problem as a function of t calculated over meshes with spatial resolution x = 1/10, 1/20, 1/40, 1/80, using a backward (Euler) finite di↵erence scheme of order 1 (solid lines), 2 (dashed lines) and 3 (dotted lines) for (a) 1, (b) 2 and (c) 3

  Figure 3.6: (Upper panels) P2 discontinuous FE approximations of the fields of (a) 1 and (b) 2 at t = 0 (black

Figure 3

 3 Figure 3.7: (a) Domain and boundary conditions for the annular Couette flow problem, in both the no-damage (section 3.3) and full Maxwell-EB (section 3.4) simulations. (b) Finite element mesh used in the full Maxwell-EB Couette flow simulations (section 3.4). The element size is a function of the radius r and is set as x(r) = R ext N (r) 2 , with N = 10.

Figure 3 . 8 :

 38 Figure 3.8: Torque on the inner domain boundary,

Figure 3 . 9 :

 39 Figure 3.9: Stationary solution for (upper panels) the velocity components ur (red circles) and u ✓ (purple circles) and the stress components (lower panels) rr (dark blue circles), r✓ (light blue circles) and ✓✓ (green circles) for (a) We 0 = 0.1 (b) and We 0 = 0.01 and (c) We 0 = 0.001.

Figure 3

 3 Figure 3.10: Schematic representation of the Couette experimental setup, showing (1) the ice plate, (2) the heating mat at the bottom of the tank, (3) the rotating disk, (4) the torquemeter and (5) the engine. Adapted from V. Pellisier, Ph.D. thesis.

Figure 3 .

 3 Figure 3.11: Photography of the ice layer during one realization of the laboratory Couette flow experiment showing the main fault (2) of mean width ⇡ 5 mm, between a thin region of "fast" ice (1) along the inner rotating disk and a large motionless region (3). Credit: V. Pellisier.

Figure 3 .

 3 Figure 3.12: (a) Time series of the measured torque during one realization of the laboratory Couette flow experiment with ! = 1 rpd. Arrows along the time axis indicate 4 cycles of stress build-up and rupture. The dashed red box indicate a period of low frequency and low asymmetry stress fluctuations. (b) Photograph during the initial elastic loading showing the intact ice plate. (c) Photograph indicating the position and shape of the shearing fault (green shading). Credit: V. Pellisier.

Figure 3 .

 3 Figure 3.13: (a) Time series of the (adimensional) torque on the inner boundary (black line) and of the damage rate (blue line) for a Couette simulation with ! = 1 rpd , and t h = 10 4 s. The time is expressed in terms of the fraction (percentage) of one complete revolution (25% of a revolution represents 6 hours in real time). Fields of the order of magnitude of the (adimensional) shearing deformation rate (log 10 ( "S), left panels) and of level of damage (d, right panels) just after (b) the first macro-rupture event and (c) a subsequent stress relaxation phase.

Figure 3

 3 Figure 3.14: Power spectral density of the damage rate time series in 5 realizations of the Couette simulations with ! = 1 rpd, ↵ = 4 and the characteristic healing time set to (a) t h = 10 4 s and (b) t h = 10 3 s. The dashed blue line indicates the (adimensional) frequency associated to the (adimensional) healing time in both cases. An averaging window of 5 frequencies is used to compute the PSD in both cases.

  Figure 3.15 represents for instance the torque and damage rate time series for one realization of the Couette simulation using ↵ = 6 which shows lower amplitude fluctuations of the damaging activity compared to the ↵ = 4 case and more symmetric variations of (skewness of -0.8).

Figure 3 .

 3 Figure 3.15: (a) Time series of the (adimensional) torque on the inner boundary (black line) and of the damage rate (blue line) for a Couette simulation with ! = 1 rpd, t h = 10 4 s and ↵ = 6.

Figure 3

 3 Figure 3.16: (a) Time series of the (adimensional) torque on the inner boundary (black line) and of the damage rate (blue line) for a Couette simulation with ! = 1 rpd, t h = 10 3 s and ↵ = 4. Fields of the instantaneous level of damage (d, upper pannels) and order of magnitude of the shearing deformation rate (log 10 ( "S), lower pannels), at the four (preand post-rupture) instants indicated on the time series.

Figure 4

 4 FIG. 1. Moderate Resolution Imaging Spectroradiometer (MODIS) image of the North Water polynya fully opened by 25 May 2001. The ice bridge prevents ice from drifting in the polynya: (a) Kane Basin; (b) Smith Sound; (c) Greenland; and (d) Ellesmere Island, Canada.

  ) boundaries with the Neumann condition • n = 0. Simulations are started from rest, from a homogeneous, undamaged ice cover with d = 1, A = 1, h = 1 m, = 0 Pa, from undamaged mechanical parameters E 0 , ⌘ 0 , 0 and from an initial field of cohesion, C 0 . The numerical scheme here includes the use of a Heaviside-type function d ⇤ = d1 de that kills the pseudo-viscous term in the constitutive equation for d = 1, as in the Couette flow simulations 3.4. As the value of the undamaged relaxation time, 0 , is very large compared to the duration of the present simulations (115 days versus 2 days in simulation time), including this threshold has no significant impact on the simulation results. The value of the transported quantities , d, C, A and h are all prescribed on the upstream part of the top and bottom boundaries, , defined as = {x 2 ⌦; u(x) • n(x) < 0}. These are chosen so that to represent inflowing undamaged ice, with d = 1, A = 1, h = 1 m, = 0 and C randomly drawn from the same uniform distribution prescribed as initial condition. Simulations are run for a time t 2 [0, +1[ and the channel flow problem reads: (P ) : Find u, , d, C, A, and h defined in ⌦⇥]0, +1[, such that

C

  min 2 2000, 5000, 10000, 20000, 30000 Pa and the maximum value set to 2 ⇥ C min . Besides impacting the formation and stability of ice arches, the local value of C sets the local damage criterion in the Maxwell-EB model, i.e., the local value of both the uniaxial compressive strength and tensile strength, c and t . Hence the minimum value of C obviously controls the timing of the onset of damaging in these simulations and the width of the distribution of C is expected to impact the rate at which damage propagates over the domain as well as the extent of the damaging, with the propagation being more progressive for a larger distribution of C, leading to a higher localization of the damage.

Figures 4 .

 4 Figures 4.3 through 4.6 represent the evolution of the wind forcing (blue curve) and of the damage rate (grey curve), together with the spatial distribution of the level of damage and the profiles of the x and y components of the ice velocity (u x , u y ) along the central meridional axis of the channel (i.e., at x = 0) at di↵erent times during the simulations, in the case of C min = 2000, 50000, 10000 and 20000 Pa respectively.

Figures 4 .Figure 4

 44 Figures 4.5d and 4.6d represent the instantaneous profiles along the central meridional axis of the domain of

Figure 4

 4 Figure 4.4: (a) Time series of the wind forcing (blue curve) and of the damage rate (grey curve) in the channel simulation using Cmin = 5000 Pa. (b) Instantaneous spatial distribution of the level of damage at the times indicated by dots 1 to 5 on the wind forcing time series. (c) Corresponding profiles of the x and y components of the ice velocity (ux, uy) along the central meridional axis of the channel.

Figure 4 . 5 :Figure 4 . 6 :

 4546 Figure 4.5: Time series of the wind forcing (blue curve) and of the damage rate (grey curve) in the channel simulation using Cmin = 10000 Pa. (b) Instantaneous spatial distribution of the level of damage at the times indicated by dots 1 to 5 on the wind forcing time series. Corresponding profiles of (c) the x and y components of the ice velocity and (d) of the principal stresses, 1 and 2, along the central meridional axis of the channel.

  Figure 4.7: (a) Instantaneous fields of the magnitude of the ice velocity, |u|, for the simulation with Cmin = 20000 Pa, at the times 1, 3, 4, indicated on the corresponding time series of the wind forcing (figure 4.6a).

  Fields of |u| show a sharp discontinuity in the drift velocity at the location of the bridge (see figure 4.7, panel 4), with near zero velocities to the north and uniform drift speed south of the bridge, indicating the detachment of the ice downstream of the bridge. Velocity profiles indeed suggest a free drift state (i.e., r • (h ) = 0) south of the ice edge.

  Figure 4.6d, panel 5, for instance indicates a value of u y of ⇠ 0.44 ms 1 for a northerly wind forcing of ⌧ a = 1 Nm 1 , consistent with the free drift ice velocity, given by u y = r ⌧ a ⇢ w C dw (4.24)

  Figure 4.8 shows the results of the simulation with C min = 20000 Pa. Similar results are obtained with lower and higher values of C min .

Figure 13 .Figure 4

 134 Figure 13. Stress states recorded during the SHEBA experiment in the Beaufort sea at one stressmeter, from mid-October, 1997 to end of June, 1998 (1 measure per hour). (a) plotted in a principal stress space. The dotted line represents σ1 = qσ2 + σc with q = 5.2 and σc = 250 kPa; (b) plotted in a τ versus σN graph. The dotted line represents |τ | = |τ0| + µσN, with |τ0| =40 kPa and µ = 0.7. If one assumes that the fault orientation (unknown) corresponding to the stress states along the dotted lines maximizes the Coulomb stress, the internal friction coefficient is given by relation (4a), µi = 0.9. (From Weiss et al (2007).) and Elder 1998, Richter-Menge et al 2002). The stress sensors determine the stress acting on a point in the horizontal plane of the ice cover by measuring changes in the radial deformation of a cylindrical annulus, allowing one to estimate the associated maximum (σ 1 ) and intermediate (σ 2 ) principal stresses (Cox and Johnson 1983). Figure 13 shows in principal stress space the stress states recorded at one sensor from October 1997 to July 1998 (one measure/hour) in the Beaufort sea during the SHEBA experiment (Richter-Menge et al 2002, Weiss et al 2007). Data recorded on different sensors or during different field campaigns give similar results. In 14

  and 4.10 show respectively the instantaneous fields of A and h at the times 3 to 7 indicated on the time series of the damage rate and wind forcing, figure 4.5a, for the simulation with C min = 10000 Pa. Simulations using other values of C min show a similar evolution of both A and h.

Figure 4 . 9 :Figure 4

 494 Figure 4.9: Fields of the instantaneous ice concentration, A, at the times indicated by numbers 3 to 7 on figure 4.5a.

Figure 4

 4 Figure 4.11: Log-normal frequency distribution of (a) the ice concentration and (b) the ice thickness corresponding to the fields of A and h shown in figures 4.9 and 4.10.

Figure 5

 5 Figure 5.1: (a) Schematic representation of the Maxwell-EB model including a threshold for permanent deformations.(b) Damage criterion of the Maxwell-EB model in the principal stresses plane, combining the Mohr-Coulomb and maximum tensile stress criteria (solid line), and the criteria for the deformation in the viscous regime, with q (i.e., µ) set to the same value as in the Mohr-Coulomb criteria (dashed black lines). With the inclusion of a threshold on the permanent deformations, the deformation of elements with a state of stress lying inside the grey shaded area would be strictly elastic.

  approach and investigated the sensitivity of the simulated thickness of Arctic sea ice to the resolution of the narrow, intricate straits of the Canadian Arctic Archipelago. A first global and coupled sea ice-ocean FE model, FESOM (The Finite Element Sea Ice-Ocean Model), was recently developed for large-scale climate simulations

[START_REF] Lietaer | The e↵ects of resolving the Canadian Arctic Archipelago in a finite element sea ice model[END_REF]

developed a large-scale FE sea ice model, operational for climate studies and destined to be coupled to the FE Second-generation Louvain-la-Neuve Ice-ocean Model (SLIM, http://www.climate.be/SLIM). They tested their unstructured grid

  The complete system of equations Table3.1: Dimensions, variables and parameters for the numerical and laboratory Couette experiments.

	Mechanical and model parameters		Values
	Poisson's ratio	⌫	0.3
	Internal friction coe cient	µ	0.7
	Undamaged elastic modulus	E 0	10 • 10 9 Pa
	Distribution of cohesion values	C	4 8 • 10 8 Pa
	Undamaged relaxation time Minimum (adimensional) apparent viscosity ⌘min 0 ⌘0	⇠ 2.0 • 10 3 s 10 12
	Damage parameter	↵	4
	Characteristic time for damage	t d	-
	Characteristic time for healing	t h	10 4 s
	Mean model resolution	x	Rext 10 r 2
	Model time step	t	1 s
	Dimensions of the domain/experiment		Values
	Radius of the tank	R ext	0.5 m
	Radius of the rotating disk	R int	0.2 m
	Angular frequency of the rotating disk	!	2⇡ rad day 1 (1 rpd)
	Angular velocity of the rotating disk	U e ✓	!R int ms 1
	therefore reads:

  (see section 2.2.3). The field of cohesion is set once at the beginning of each simulation as in quenched disorder, to allow a direct comparison between simulations run with a di↵erent range of values of C. All simulations use the same spatial distribution of the noise in C, set by randomly drawing a value over each model element over the non-dimensional interval [1, 2]. The local value is multiplied by a minimum value of the cohesion, C min , such that C 2 [C min , 2 ⇥ C min ], and the value of C min is varied between simulations. Table4.1 summarizes all model parameter values employed in these simulations.An unstructured mesh with triangular elements is used. The average spatial resolution, x, is of 2 km.Hence the characteristic time for damage, t d , and the model time step, t, are set to x c = 4 seconds. The numerical scheme is described in the following section.The computational domain, ⌦ ⇢ R 2 , for the channel flow simulations is shown in figure 4.2b. The domain boundary @⌦ is partionned as @⌦ =

	4.4 Numerical scheme
	boundaries (	left ,	in [ right ). The channel is open at its top ( left [ out [	right . A no-slip condition is applied at the lateral in ) and bottom ( out
		Parameters		Values
		Poisson's ratio	⌫	0.3
		Internal friction coe cient	µ	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 7 ⇥ E 0 Pa s
		Undamaged relaxation time	0	10 7 s
		Minimum apparent viscosity	⌘ min	10 4 Pa s
		Minimum cohesion	C min 2000, 5000, 10000, 20000, 30000 Pa
		Damage parameter	↵	4
		Characteristic time for damage	t d	t s
		Characteristic time for healing	t h	5 • 10 5 s
		Mean model resolution	x	2000 m
		Model time step	t	4 s
				c ⇤	20
		Air drag coe cient	C da	1.5 • 10 3
		Air density	⇢ a	1.3 kg m 3
		Water drag coe cient	C dw	5.5 • 10 3
		Water density	⇢ w	1027 kg m 3
			Table 4.1: Model parameters for the channel flow simulations.

Acknowledgements -Remerciements

I am writing this page as the last of this thesis. Yet I feel it is the most important one.

I would like to thank Jérôme, first for o↵ering me the opportunity to come to Grenoble and work on this Ph.D. project. This opportunity came as a surprise; unplanned. The work, as the life here, has been challenging in many ways. But most importantly, it has been even more stimulating than I hoped it to be. Hence saying that I am thankful for the role he played in this unexpected turn of events is simply an understatement.

Second, I would like to thank him for his support. And by support I mean more than correcting me for the inadequate use of the word "constraint" a thousand times.

The second person I would like to thank is Pierre, for his precious advices, his contagious enthusiasm and his never-failing patience in answering my multiple, slightly bewildered or completely confused, math questions.

I am also grateful for the support I received and stimulating contacts I had with many people at Total. I especially thank Annie Audibert-Hayet, Edmond Coche, Philippe Lattes, Philippe Ricoux and Kaj Riska for their valuable suggestions on this work and their many encouragements. I thank my reviewers, Anne Mangeney and Vincent Legat, as well as my examinators, Annie Audibert-Hayet and François Renard for their time and their numerous, useful advices.

Then let me just make a list of the people I would like to express my gratitude to, although simply listing their names does not really pay justice to my appreciation of their presence and help:

Concluding remarks

In this chapter, we have presented all the components of the Maxwell-EB model. Analyses of small-deformation simulations using the simplest, most idealized domain, forcing and boundary conditions have revealed :

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 temporal correlations 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 be inherited from the prescribed forcing. Instead, their emergence is a signature of the rheological framework itself. Moreover, the simulations have demonstrated that anisotropy in the Maxwell-EB model arises naturally from elastic interactions, although the material's properties are fully isotropic at the element scale. This is likely due to the anisotropy of the elastic kernel representing the redistribution of stresses within the material after a damaging event.

Furthermore, sensitivity analysis on the damage parameter ↵ have demonstrated that the Maxwell-EB dynamical framework, with few numbers of independent variables, can represent a large range of mechanical behaviours, from a regular, predictable stick-slip-like behaviour with a single damaging frequency corresponding to the prescribed healing rate, to a marginally stable, unpredictable, creep-like deformation with temporal correlations in the damaging activity at all time scales below the material's healing time. Our simulations show that over a limited range of values of ↵, the model can reproduce both the persistence of creeping leads and the activation of new leads with di↵erent shapes and orientations and that, under completely homogeneous forcing conditions. Ultimately, comparisons to temporal scaling analysis based on measurements of the stress and/or deformation rates would be required in order to determine the values of ↵ most relevant for sea ice.

2. In the second case, is a locally discontinuous step function in H 0 \ C 1 ( ⌦) defined as

(3.4)

3. In the third case, 3 2 L 2 represents a field of uncorrelated noise. It is defined locally such that the value of 3 over each model element is drawn randomly from a uniform distribution spanning the range [0, 1]. It is thus discontinuous everywhere and non-di↵erentiable (see figure 3.2c).

Here

3 is of the same form as the field of cohesion described in chapter 2, which represents the material's heterogeneity in the current implementation of the Maxwell-EB framework. The regularity of the fields of level of damage (d) and of the internal stress ( ) simulated by the model, which are characterized by the presence of highly localized linear features arising from the propagation of damage (see for instance, figure 2.9), can be considered to lie somewhere between that of 3 and 2 : as these linear features introduce spatial correlations, at least in one dimension, d and are more regular than pure uncorrelated noise, but, due to the tendency of the model to localize damage and deformation at the scale of the element, they show more discontinuities than the 2-dimensional step function defined by (3.4). Di↵usion of the strong gradients within these fields is therefore expected to be somewhat less important than for 3 , but more than for 1 and 2 .

Problem (P) is solved using a Euler implicit scheme and discretizing the time t such that t n = n t, with

Chapter 4

A Maxwell-EB sea ice model

After the highly idealized, small-deformations and large-deformations simulations presented in chapters 2 and 3, we finally use the Maxwell-EB rheological framework in the context of modelling the sea ice cover on geophysical scales. The problem presented here represents the flow of sea ice through a narrow channel and constitutes a relevant test case in multiple respects.

A realistic case study

The formation of ice bridges in narrow passages is a common phenomenon over the Arctic, occurring for instance in the Amundsen Gulf and Bering Strait as well as in many channels of the Canadian Arctic Archipelago (CAA).

Because ice bridges stop the flow of ice and can remain stable for several weeks, maybe months in the northern CAA, understanding the conditions under which they form is important in terms of assessing the variability of the Arctic ice circulation and export through these pathways [START_REF] Kwok | Large sea ice outflow into the nares strait in 2007[END_REF]. Expenses of ice-free

water and polynyas open downstream of these bridges and strongly impact the local atmosphere-ocean heat exchanges, while also leading to the recurrent generation of new ice. As the geometry of these bridges is always concave and arch-like, they are also commonly referred to as ice arches. important and extensively studied outflow pathway for thick, multi-year Arctic sea ice [START_REF] Ingram | An overview of physical processes in the North Water[END_REF][START_REF] Kwok | Large sea ice outflow into the nares strait in 2007[END_REF]. In the case of Nares Strait, ice bridges have been known to form preferentially throughout the winter, at the entrance of the Strait and south of Kane Basin [START_REF] Gudmandsen | Lincoln sea and Nares Strait[END_REF], leading to the opening of Smith Sound and to the formation of the North Water polynya (NOW). The ice breakup is thought of being primarily mechanically driven, as opposed to thermally (i.e., ocean upwelling) driven [START_REF] Ingram | An overview of physical processes in the North Water[END_REF][START_REF] Gudmandsen | Lincoln sea and Nares Strait[END_REF]. The opening and the drift of sea ice within the channel has been associated to the prevalence of northerly winds, associated with strong pressure gradients between the Lincoln sea, to the north, and Ba n Bay, to the south and orographically channeled by the steep coastal topography of Ellesmere Island and Greenland [START_REF] Gudmandsen | Lincoln sea and Nares Strait[END_REF][START_REF] Samelson | Evidence for atmospheric control of sea-ice motion through nares strait[END_REF].

Granular materials have been known for a long time to form concave stress-free surfaces and exhibit selfobstruction to flow under certain conditions (e.g., [START_REF] Richmond | Limiting spans for arching of bulk materials in vertical channels[END_REF][START_REF] Walker | An approximate theory for pressures and arching in hoppers[END_REF]. By assimilating sea ice as a 2-dimensional continuum material that obeys Coulomb's failure criterion, [START_REF] Sodhi | Ice arching and the drift of pack ice through restricted channels[END_REF] showed that the formation and shape of stable ice bridges in the Arctic ice cover can be explained by two-dimensional