Ice bridges and ridges in the Maxwell-EB sea ice rheology

This paper presents a first implementation of the Maxwell-EB model on geophysical scales. The model is tested on the basis of its capability to reproduce the complex mechanical and dynamical behaviour of sea ice drifting through a narrow passage. Idealized as well as realistic simulations of the flow of ice through Nares Strait are presented. These demonstrate that the model reproduces the formation of stable ice bridges as well as the stoppage of the flow, a phenomenon occurring within numerous 5 channels of the Arctic. In agreement with observations, the propagation of damage along narrow arch-like kinematic features, the discontinuities in the velocity field across these features dividing the ice cover in floes, the strong spatial localization of the thickest, ridged ice and the opening of polynyas downstream of the Strait are all represented. The model represents different dynamical behaviours linked to an overall weakening of the ice cover and to the shorter lifespan of ice bridges, with implications in terms of increased ice export through narrow outflow pathways of the Arctic. 10


Introduction
The formation of ice bridges is a common phenomenon in the Amundsen Gulf, Bering Strait as well as in many narrow passages of the Canadian Arctic Archipelago (Sodhi, 1977). Commonly referred to as ice arches because of their curved and concave shape, these structures can remain stable for several weeks or months and stop the flow of ice through outflow channels (Kwok, 2005;Kwok et al., 2010;Münchow, 2016). Downstream of ice bridges, expenses of ice-free water and polynyas open, which 15 strongly impacts the local atmosphere-ocean heat exchanges and promotes the generation of new ice (Smith et al., 1990).
Upon breakup of a bridge, the outflow of ice, stored in the basin upstream, drastically increases. Therefore, the capability of representing adequately the complex dynamical behaviour of ice drifting through narrow passages might constitute a key asset when using numerical models to asses the seasonal and interannual variability in the circulation and export of fresh water and sea ice in the Arctic.  Fig. 1a), one of the most extensively studied outflow pathways for seasonal and multi-year Arctic sea ice Kwok, 2005;Kwok et al., 2010;Münchow, 2016). The annual mean ice volume flux through this channel, only a few tens of kilometers wide (30 − 40 km) in some places, is thought to be equivalent to about 7% of the  Hibler, 1979) to increase the area of the yield envelope in the second (or fourth) quadrant in the space of the principal stresses σ 1 and σ 2 (see Fig. 2), corresponding to a cohesive state. This effectively increases the uniaxial unconfined compressive strength σ c (the maximum value of σ 1 that the material can sustain when σ 2 = 0, see Fig. 2, stress state 3). Rasmussen et al. (2010) also performed numerical simulations of the sea ice dynamics in Nares Strait and North Water Polynya, 5 3 The Cryosphere Discuss., doi: 10.5194/tc-2016-276, 2016 Manuscript under review for journal The Cryosphere Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License. using the dynamic-thermodynamic CICE sea ice model, the ice mechanics component of which is also based on a viscousplastic rheology and elliptical yield curve (Hunke et al., 2010a). The authors noted a lack of stability and shorter lifespan of the simulated ice bridges, leading to a slower opening and lower extent of the North Water polynya and to an earlier draining of Nares Strait compared to estimates from satellite imagery. They attributed this deficiency to a too low ice strength in their model, either caused by a too thin ice cover or to the inability of their rheology to reproduce the correct internal strength of sea 5 ice.
Channel flow simulations have not yet been performed using elastic-brittle models. Hence such experiments constitute an interesting test case of their mechanical behaviour. 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 these models, even with a fine spatial resolution (e.g., 0.32 • degrees or 4 km to 10 km between the Lincoln Sea and Baffin Bay in the model of Rasmussen et al. (2010) and 10 0.15 • × 0.04 • degrees in the realistic simulations of Nares Strait of Dumont et al. (2009)), are also able to account for the presence of multiple arch-like leads within and upstream of the channel, as observed from satellite imagery (e.g., see Fig. 1b and Sodhi, 1977). With its capability to represent the extreme localization of damage and deformation , the Maxwell-EB model might be better suited to simulate these fine features.
Ice drift and coverage conditions within a channel moreover represent a severe test of the numerical scheme in terms of 15 handling discontinuities within the simulated fields, as once a stable ice bridge forms, the ice downstream detaches from the bridge and is driven out of the channel without mechanical resistance. At this point, extremely sharp gradients in ice velocity, thickness and concentration are expected to arise (see Fig. 1b and 3a).

Ice ridges
Sea ice models are most often compared to each other and to observations in terms of the spatial distribution of the simulated 20 ice thickness (e.g., Johnson et al., 2012). An equally important, and perhaps more appropriate, metric to investigate the mechanical behaviour of the sea ice cover is the probability density function (PDF) of the ice thickness, of which some valuable information have been available for some time from drill-hole, sub-marine mounted sonar and airborne electromagnetic sounding measurements. In particular, the tail of PDFs represents the ice that has thickened not only due to thermodynamic but also to mechanical redistribution processes, i.e, the ice incorporated in pressure ridges, which are long, linear rubble piles of ice, 25 meters or tens of meters wide, formed under convergent and shearing motions. A recurrent statistical property of the tail of these PDFs is that it appears to fit a negative exponential function (Wadhams, 1994;Haas, 2009). To this day, it is not clear why it takes this particular form. The important point however is that it indicates the tendency of mechanical redistribution to create extremes (Thorndike et al., 1975), or, in other words, the strong localization of the thickest ice in space.
Over the years, there have been different attempts to represent the formation of pressure ridges in numerical simulation 30 of the sea ice cover (e.g., Parmerter and Coon, 1972;Kovacs and Sodhi, 1980;Hopkins, 1994), but no model is yet capable of describing the entire process. Continuum sea ice models typically have a spatial resolution of a few kilometers to tens of 4 The Cryosphere Discuss., doi: 10.5194/tc-2016-276, 2016 Manuscript under review for journal The Cryosphere Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License. kilometers and hence do not resolve ridges per se. In such models, two main approaches are taken to handle the redistribution of ice thickness associated with ridge building.
The so-called multi-thickness categories scheme based on the pioneer work of Thorndike et al. (1975) and Rothrock (1975) is widely used in current sophisticated sea ice and coupled models. This scheme introduces an areal thickness distribution function which evolves in time due to both thermodynamics and dynamics processes. Ridging is treated "explicitly" by allowing 5 a prescribed thin ice portion of the thickness distribution to be redistributed into thicker ice categories in response to the simulated deformation. The main advantage of this scheme is that it allows accounting for variations in the ice thickness at the sub-grid scale. The relation between the redistribution process and the strength of the ice (often characterized by a pressure, P ) in this modelling framework is however dubiously based on energy conservation principles (the deformational work is equated to the work done in building ridges, which is partitioned between potential energy changes (Thorndike et al., 1975), 10 the frictional dissipation in ridging (Rothrock, 1975) and dissipation in shearing deformation (Pritchard, 1981), all of which being very hard to estimate). As the simulated strain rate tensor does not provide directly the information on the relative amount of opening and ridging (and also sliding) within the ice cover, expressions for the modes of redistribution need to be assumed, the correct form of which remain uncertain to this day (Hunke et al., 2010b). These redistribution functions can be set in an ad-hoc manner (Thorndike et al., 1975), estimated empirically from strain rates observations (Stern et al., 1995) or, in the case 15 of plastic models such as the VP and EVP models, determined based on the prescribed form of the yield criteria and flow rule (Rothrock, 1975;III, 1980;Flato and Hibler, 1995). The multi-categories schemes also introduce several additional parameters (e.g., a frictional dissipation coefficient, a prescribed percentage of the thickness distribution participating in the ridging, ...), which are all badly constrained and to which the simulated thickness distribution and patterns can be highly sensible (Flato and Hibler, 1995;Bitz et al., 2001). Moreover, it necessitates solving an additional evolution equation for the thickness distribution 20 function as well as thermodynamics and transport equations for multiple ice thicknesses, which increases the cost of numerical schemes as the number of ice categories is increased.
The second approach is the simpler two-level model suggested by Hibler (1979) in which the simulated ice cover falls into two thickness categories: the effective thickness, representing the average ice thickness over a model grid cell, and zero thickness, or open water. As opposed to the multi-category scheme, ridge building is treated implicitly (based on a volume 25 conservation principle, see section 3). Known shortcomings of this model are the underestimation of ice thickening in regions of convergent and shearing ice motion compared to the multi-categories scheme, due to the unresolved thin ice that participates in ridge building (e.g., Bitz et al., 2001). Nevertheless, it is still unclear to this day if, when incorporated in viscous-plastic type models, either of the two-level scheme or the multi-categories scheme, even when tuned, is able to reproduce the form of tail of the PDFs calculated from Arctic sea ice thickness measurements (e.g., Flato and Hibler, 1995). 30 3 The Maxwell-EB sea ice model The Maxwell-EB model builds on the continuum Elasto-Brittle (EB) rheology, which has been used to model the fracturing of rocks (e.g., Amitrano et al., 1999) and was implemented for sea ice modelling by Girard et al. (2010) and Rampal 5 The Cryosphere Discuss., doi:10.5194/tc-2016-276, 2016 Manuscript under review for journal The Cryosphere Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License.
(2015); Rampal et al. (2015). As the EB model is based on a linear-elastic constitutive law, it does not solves simultaneously for both the elastic (reversible) deformations associated with the fracturing of the ice pack and the permanent (irreversible) deformations occurring once the ice pack is fractured and ice floes move relative to each other. Therefore, it does not allow estimating ice drift velocities unambiguously. The Maxwell-EB model was developed to deal with these intrinsic shortcomings of the EB framework. This mechanical model has been described in full details by Dansereau et al. (2016). Here we only recall 5 its essential features.
In this augmented rheology, a passage between the small/elastic and large/permanent deformations is made possible by the addition of a viscous-like relaxation term in the linear-elastic constitutive law for a compressible, continuous solid (see Dansereau et al. (2016), Eq. (4)). Associated with the linear elastic term in this constitutive equation is the true elastic modulus of the material, i.e., of sea ice, at the scale of the model grid cell. The viscosity, η, associated with the viscous term, is not the 10 true bulk viscosity of sea ice, but an apparent viscosity that represents the flow resistance of the fractured/fragmented ice cover averaged over the grid cell. The ratio of the two mechanical properties, λ = η E , has the dimension of a time, and sets the rate of dissipation of the stress through permanent viscous-like deformations. Alternatively, it quantifies the capability of the ice pack to retain the memory of elastic deformations. In the Maxwell-EB model, these three mechanical parameters vary with the local level of damage of the material, quantified by a scalar variable d that evolves between 1 for an undamaged and 0 for a 15 "completely damaged" ice cover.
Damage occurs when the state of stress becomes overcritical with respect to a Mohr-Coulomb failure criterion or a tensile cut-off. The combined criteria are represented in Fig. 2 (thick solid lines) in the principal stresses space, with the convention that compressive stresses are positive. With such criterion, the uniaxial (unconfined) compressive strength, σ c , and maximum tensile strength, σ t , are directly related to a cohesion parameter C as In the model, the damage criterion varies spatially, as some disorder in C is introduced at the local scale to represent the natural heterogeneity of the ice cover associated with the presence of various defects at the sub-grid scale Weiss and Dansereau, 2016). 25 A healing mechanism counterbalances the effects of damaging over much larger time scales and represents the refreezing of leads. This mechanism is distinct and independent from pure thermodynamic growth as it allows the level of damage to re-increase and recover at most the undamaged value of d = 1. Both processes, damaging and healing, are combined in a single equation for the evolution of damage (see Dansereau et al. (2016), eq. 17).
The coupling of E, η and λ with d is such that the mechanical strength, as well as the capability of the material to retain  (20)). It is this coupling of the mechanical properties and level of damage of the ice cover that allows the model to dissipate internal stresses in large, permanent deformations along faults, or leads, once the ice pack is highly 6 The Cryosphere Discuss., doi: 10.5194/tc-2016-276, 2016 Manuscript under review for journal The Cryosphere 2 + σ 2 12 . The thin solid lines represent the damage criterion in the case of C = 0. The ellipses represents the yield criterion in the standard VP model of Hibler (1979) for different aspect ratios (< 2, outer; 2, center; > 2, inner ellipse). The numbers 1 to 4 indicate the states of (1) uniaxial tension, (2) biaxial tension and compression, (3) uniaxial compression and (4) biaxial compression.
. damaged, while reproducing the small deformations associated with the fracturing process and retaining the memory of elastic deformations over relatively low damage areas.
Analyses of the deformation and damage fields simulated using idealized geometries, simple forcing conditions and mechanical parameters values consistent with sea ice on geophysical scales  have demonstrated that the Maxwell-EB rheological framework successfully reproduces the anisotropy of sea ice deformation as well as the strong strain 5 localization in both space and time and associated spatial scaling laws (Stern et al., 1995;Kwok, 2001;Marsan et al., 2004;Rampal et al., 2008Rampal et al., , 2009Weiss, 2008). The observed spatial and temporal coupling between these scalings is also represented . Sensitivity analyses on the damage parameter α (see Weiss and Dansereau (2016)), which sets the rate of viscous dissipation of the internal stress as a function of the increasing level of damage of the ice cover, have shown that the model, with few independent variables, can represent a large range of mechanical behaviours: from a regular, 10 7 The Cryosphere Discuss., doi: 10.5194/tc-2016-276, 2016 Manuscript under review for journal The Cryosphere Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License. predictable stick-slip with a single damaging frequency related to the prescribed rate of healing, to a marginally stable, unpredictable deformation with temporal correlations in the damaging activity at all time scales below the material's healing time . Over a range of values of this parameter, the model reproduces both the persistence of creeping leads in the ice cover and the activation of new leads with different shapes and orientations Weiss and Dansereau, 2016). 5 In the channel flow simulations presented here, this new rheological model is implemented in a continuum modelling framework typical of regional and global sea ice models. In such framework, the ice cover is assimilated to a 2-dimensional plate.
Hence plane-stresses are assumed. The motion of the ice is described by the following Navier-Stokes type equation: with u, the ice velocity, σ, the stress tensor, ρ, the ice density, τ a and τ w , the air and water drags, −ρhf k × u, the Coriolis 10 pseudo-force with f , the Coriolis parameter and k, the upward unit vector normal to the ice surface and −ρhg∇H, the force due to gradients in the sea surface dynamic height, H, with g, the gravitational acceleration, and which can be expressed alternatively in terms of the geostrophic ocean current velocity, u w .
In this 2-dimensional momentum equation, the variables h and A represent respectively the mean ice thickness and ice concentration over a model grid cell. The mean thickness is assumed to be the weighted sum of two ice categories: thick ice, 15 with thickness h thick = h A , and open water, with thickness h thin = 0 m (Hibler, 1979). Mass conservation is ensured by the following evolution equation for h: where the term S h represents the ice thickening/thinning through thermodynamic processes and mechanical redistribution, i.e., through pressure ridge formation. 20 In the present implementation of the Maxwell-EB model, the mechanical redistribution is accounted for in a very simple manner. This is done on purpose, to test the input of the new rheology in the representation of the thickness distribution.
Following the parametrizations of Hibler (1979), Thomson et al. (1988) and Lietaer et al. (2008), an evolution equation for the ice concentration is first solved

25
where S A accounts for both the thermodynamics and mechanical sources and sinks of ice concentration. If as a result of convergent ice motion, the ice coverage A over a given model element exceeds unity, the excess concentration, is used to increment the thickness of thick ice over that element, and the ice concentration is reset to 1. The mechanical contribution (sink) to S A in this case reads: 30 8 The Cryosphere Discuss., doi: 10.5194/tc-2016-276, 2016 Manuscript under review for journal The Cryosphere Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License. and the associated mechanical contribution (source) to S h in eq. 4 is Finally, the mechanical parameters E and η are coupled to the ice concentration A, as follows : 5 with f 1 and f 2 , the functional dependence of E and η on the level of damage of the ice cover d given respectively by Eq.
(18) and (19) in Dansereau et al. (2016) and c * a non-dimensional, constant parameter. The functional dependence on the ice concentration follows the one suggested by Hibler (1979) and widely employed in VP models for the pressure term (P ) setting the ice strength in compression. The exponential function of A simply allows the value of both E and η to be maximal when the ice concentration is of 100% (A = 1) and to decrease rapidly when leads open and A drops. It is employed to characterize the 10 dependence of the elastic modulus (or effective elastic stiffness) on A in the elasto-brittle models of Girard et al. (2010) and of Bouillon and Rampal (2015). In the case of η, this parametrization is compatible with the rapid decay of the apparent viscosity of granular media when decreasing their packing fraction from the close-packed limit (Aranson and Tsimring, 2006). In the present implementation of the Maxwell-EB model, this simple parametrization as well as the value of the non-dimensional parameter c * is the same for both mechanical parameters, but this could eventually be refined in future developments of the 15 rheology.

Channel flow simulations
The drift of sea ice within Nares Strait is thought of being primarily driven by the prevalence of sturdy northerly winds associated with the strong pressure gradient between the Lincoln sea to the north and Baffin Bay to the south, and which are orographically channelled by the steep coastal topography of Ellesmere Island and Greenland (Ingram et al., 2002;Gudmand-20 sen, 2004; Samelson and Barbour, 2008;Münchow, 2016). The most recurrent location for an ice bridge is in the southern Kane Basin (see Fig. 1c) (Kwok et al., 2010)), where a stable arch is observed to form almost every year between November and March ) as a result of the convergence of a mixture of first and multi-year ice into Kane Basin (Kwok et al., 2010;Gudmandsen, 2004). Disintegration of the ice cover downstream of this arch leads to the opening of the North Water polynya (see Fig. 1c) in Smith Sound Ingram et al., 2002). 25 In both the idealized and realistic simulations presented here, we aim to reproduce the formation of such ice bridge. In the idealized case, the domain consists in a 120 km wide rectangular basin that converges into a 40 km wide, 40 km long channel (see Fig. 3b). The geometry, similar to that used in the idealized simulations of Dumont et al. (2009), is conceived to be roughly consistent with the shape of the constriction between Kane Basin and Smith Sound (dashed box, Fig. 3a). This simple configuration, symmetric with respect to the y−axis, facilitates the analysis of the dynamical behaviour of the model 30 9 The Cryosphere Discuss., doi: 10.5194/tc-2016-276, 2016 Manuscript under review for journal The Cryosphere Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License. and, in particular, of the simulated states of stress. The dynamics described in the next section is however not specific to this geometry : simulations were also performed over different domains (narrower, longer channels, smaller basins) and produced similar results. The realistic domain covers the entire Strait down to northern Smith Sound (see Fig. 1c, grey shading). This configuration allows investigating the formation of secondary arches in various locations as well as other phenomena related to the presence of topographic features such as islands and fjords. In this case the mesh was built using the Gmsh mesh generator 5 (Geuzaine and Remacle, 2009) and the GSHHG high resolution shoreline data (http://www.soest.hawaii.edu/pwessel/gshhg/) between 73 and 85 • N and 280 and 320 • W. This data is available in geodetic longitude/latitude on the WGS-84 ellipsoid. Here it is converted in Cartesian coordinates using a stereographic projection. The projection is centered on the North Pole, with the Greenwich meridian aligned on the positive x-axis and the Strait is roughly oriented along the y-axis (see Fig. 3a).
The prescribed wind forcing is made as simple as possible to facilitate the analysis. Consistent with observations of oro-10 graphic channelling, a uniform along-channel, i.e., northerly, wind stress, τ a , is applied. The stress is increased steadily between 0 and 1 Nm −2 over a period of 24 hours, and then held constant, to simulate the passage of a storm (Kwok, 2006;Samelson et al., 2006;Samelson and Barbour, 2008). Considering a simplified quadratic wind drag based on the absolute, instead of the relative, wind speed of the form 15 with u a , the wind velocity ρ a = 1.3 kg m −3 , the surface air density and C da , the air drag coefficient, commonly set to 1.2·10 −3 in sea ice models following Hibler (1979), this corresponds to a maximum wind speed of ∼ 22 ms −1 (∼ 82 km h −1 ). The ocean is at rest (u w = 0), hence the oceanic drag is given by the following quadratic formula where ρ w = 1027 kg m −3 is the density of sea water and C dw is the drag coefficient, set to 5.5 · 10 −3 (McPhee, 1980). the Coriolis term is also set to zero. As the ocean is at rest, the force associated with gradients in the sea surface dynamic height is also zero. All other quantities (σ, d, C, h, A) are advected with the ice flow. Thermodynamics processes are not accounted for, hence S h and S A in Eqn. 4 and 5 represent only the mechanical redistribution of ice thickness and concentration. 30 Mechanical parameter values are based on measurements within sea ice. An undamaged elastic modulus (Young's modulus) of E 0 = 5.0 · 10 8 Pa is considered. Poisson's ratio is set to ν = 0.3 (Timco and Weeks, 2010). The undamaged relaxation time,

10
The Cryosphere Discuss., doi: 10.5194/tc-2016-276, 2016 Manuscript under review for journal The Cryosphere Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License. λ 0 , is set to 10 7 s (∼ 115 days), a value that allows the numerical scheme to converge, while also ensuring that non-physical viscous dissipation over low damage areas of the ice cover is insignificant . The characteristic time for the healing process, which corresponds to the time required for the local level of damage d to re-increase from 0 to 1, is set to 5·10 5 s (∼ 5.7 days) based on estimates of ice growth within open leads : Petrich et al. (2007) reported a time for the growth of 1 m of ice within an opening of 10 cm under air temperatures of −15 • C between 10 5 and 10 6 seconds. As healing is meant to 5 represent the local recovery of the ice mechanical strength within refreezing leads, not the thermodynamic growth of ice within polynyas, healing in the simulation is turned off as soon as the ice concentration locally drops below 75%, which occurs when and where the ice detaches from a stable ice bridge. Table 1 summarizes all model parameter values employed the simulations presented here.
As a determining factor for the formation of stress-free surface is the cohesive nature of the material, simulations with a 10 different range of values of cohesion were compared. The field of C was set as in quenched disorder as follows. First, its spatial distribution for all simulations using the idealized or realistic domain was obtained by randomly by drawing a value in the non-dimensional interval [1, 2] over each model element. This noise was then multiplied by a minimum value of cohesion, C min , such that C ∈ [C min , 2 × C min ]. This minimum cohesion was varied between 2, 5, 10, 20, 30 kPa. Hence, the same spatial distribution of C was used in all simulations, but the magnitude of C was varied between simulations. The largest 15 values of C employed are consistent with in-situ stress measurements in the Beaufort reported by Weiss et al. (2007) (see Fig.   8a). According to the presumed scale effect on shear strength, set by the size of the defects (thermal cracks, brine pockets, ...)

11
The Cryosphere Discuss., doi: 10.5194/tc-2016-276, 2016 Manuscript under review for journal The Cryosphere Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License. present in the ice cover (Schulson, 2004;Weiss et al., 2007), lower values of C are consistent with larger defect sizes and a lower shear strength.

Parameters Values
Internal friction coefficient µ 0.7 Ice density ρ 900 kg m −3 Undamaged elastic modulus Characteristic time for healing t h 5 · 10 5 s c * 20 Air drag coefficient Water density ρw 1027 kg m −3 Table 1. Model parameters for the idealized and realistic channel flow simulations.
In all simulations, the ice is initially undamaged (d = 1) and has a uniform thickness of 1 m. Although a higher value of h would perhaps be more representative of the mean thickness of the mixture of first year and multi-year ice reported in the Strait (0.9 to 1.6 meters, Haas et al., 2006), we find that the evolution of the simulated fields does not depend on this prescribed 5 initial value of the ice thickness, as long as the coverage is initially uniform. The domain is initially completely covered with undamaged ice (A = 1, d = 1), so that the location of ice bridges is not prescribed. Simulations are started from rest. A no-slip condition (u = 0) is applied at the lateral boundaries Γ left , Γ right , representing the coasts (see Fig. 3a  The model is entirely developed within the C++ environment RHEOLEF (Saramito, 2013a). The numerics is based on finite elements and variational methods. The equations of motion are cast in the Eulerian frame and discontinuous Galerkin methods are used to handle advective processes (Saramito, 2013b), as well as the non-linear terms arising in the objective derivative of 15 the stress tensor (see ). Polynomial approximations of degree 1 are used for the velocity field and all The numerical scheme is divided into linear sub-problems. A semi-implicit scheme is used to linearize the momentum and constitutive equations. The internal stress term in the momentum equation is discretized using an implicit scheme while the water drag term is linearized as τ w = ρ w C dw |u n |u n+1 and the ice thickness and concentration are taken at the n th time step.  15 in order to simulate the propagation of damage over the ice cover with the highest possible resolution (see Dansereau et al. (2016)). Similarly, the ice thickness and concentration could also be solved for using an implicit scheme and be updated as part of the fixed point iteration. However, variations in the mechanical parameters associated with changes in h and A are small and slow compared to the changes due to damaging. Hence we find that solving for these equations explicitly does not impact the solution. 20

Dynamical behaviour
Here we investigate the formation of arches and stable ice bridges in the Maxwell-EB model and analyze the evolution of the simulated level of damage, ice velocity and internal stresses. Idealized and realistic simulations with a different range of values of C (see section 4) were compared. 25 The evolution of the applied wind forcing (dashed line) and of the damage rate (solid grey line), defined as the number of damaged elements per model time step times their respective distance to the Mohr-Coulomb or tensile damage criterion (see Dansereau et al. (2016)), is represented in the case of one idealized simulation using C min = 20 kPa in Fig. 4a. The spatial distribution of the level of damage is shown at three different stages of the forcing (numbered 1, 2, 3) for this idealized simulation in Fig. 4b and at the corresponding stages in a realistic simulation using C min = 20 kPa in Fig. 6a. We first discuss 30 the idealized simulations results and then comment on the realistic case.

Idealized simulations
In all simulations using different values of C min , deformation rates first concentrate along narrow, concave damaged features that form in the interior and downstream of the channel (Fig. 4b, panel 1). The profile of the y-component of the ice velocity, u y (red curve), plotted along the central meridional axis (see Fig. 5a, stage 1) shows that after the onset of damaging (indicated by the first peak in damage rate on Fig. 4a), the ice over this portion of the domain is set in motion while the undamaged ice 5 upstream remains motionless. The sharp no-flow transition coincides with the constriction point that defines the entrance to the channel, across which an ice arch has clearly formed. It is important to note that here, "no-flow" or, as later mentioned, "flow stoppage" is not defined as a zero drift speed, but rather as a drift velocity on the order of that associated with strictly elastic deformations within an undamaged ice cover.
As the wind forcing is further increased, damage propagates upstream of the channel, which corresponds to the second peak Up to that point, the overall spatial and temporal evolution of damage is similar between the simulations using different ranges of cohesion. However, as the local value of C sets the local damage criterion in the Maxwell-EB model, i.e., the local value of σ c and σ t , the minimum value of cohesion over the domain controls the timing of the onset of damaging in the simulations, with damaging occurring sooner as C min is lower. Furthermore, the width of the distribution of C impacts the rate 20 of propagation of the damage, with the propagation being more progressive for a larger distribution.  Fig. 5a, panel 3). 5 We further investigate the mechanism behind the formation of the ice arches by analyzing the simulated states of stress for the idealized simulation with C min = 20 kPa. Similar results were obtained for C min = 30 kPa. The symmetry of the idealized domain facilitates the analysis, but the same general conclusions also apply to the realistic simulations. Figure 5b represents the instantaneous profiles along the central meridional axis of the principal stresses σ 1 and σ 2 at the times indicated by the numbers 1, 2 and 3 on Fig. 4a. 10 Internal stresses within the initially undamaged ice cover are compressive (σ 1 , σ 2 > 0) over the basin, change sign at the middle of the channel (at y = 0) and are tensile (σ 1 , σ 2 < 0) downstream of the channel (not shown). At the onset of damage, σ 2 becomes negative (i.e., tensile, see Fig. 2) over the converging part of the basin (y < 50 km) and in the interior of the channel. Within and downstream of the channel, multiple minima in σ 2 are observed. These minima are collocated with either a maximum (positive) value of σ 1 (shearing state of stress) or a minimum (negative) value of σ 1 < 0 (tensile state of stress) and 15 correspond to the location of arch-like features (see Fig. 4b, panel1). Figure 5c, panel 1 shows the elements that have exceeded either the local Mohr-Coulomb (in blue) or tensile (in red) criterion and confirms that these features were formed by a shear, tensile or a combination of both failure mechanisms.
Later in the simulation, the ice cover fails preferentially in shear upstream of the channel, while tensile failure progressively becomes predominant within and near the exit of the channel (see Fig. 5b and c, panels 2 and 3). Once the stable ice bridge 20 is formed, its location is clearly indicated by collocated extrema of the value of the principal stresses (σ 1 > 0, σ 2 < 0, on Fig. 5b, panel 3). Both the profiles of the principal stress components and the spatial distribution of the over-critical elements downstream of the bridge give evidence that this structure sustains biaxial tensile stresses (see Fig. 4b and 4c, panel 3). This is an important point, as standard viscous-plastic sea ice models do not account for pure uniaxial or biaxial tensile strength and hence would not be able to reproduce the formation of a stable ice arch with self-obstruction to flow under the conditions simulated 25 here. Consistent with the formation of a stable stress-free surface and with the detachment of the ice from the stable bridge, both principal stress components decrease to zero near the exit and outside of the channel (see Fig. 4b, panel 3). Velocity profiles (Fig. 5a, panel 3) are effectively uniform downstream of the no-flow transition with u x ≈ 0 and u y ≈ 0.44 ms −1 corresponding to the free drift velocity 30 for τ a = 1 Nm −2 .

Realistic simulations
The temporal evolution of damage, velocity and internal stress in the realistic simulations is similar to the idealized case.
The ice first fails near the exit of Kane Basin and rapidly weakens in Smith Sound (see Fig. 6a, t = 6 hours). Damage then progressively propagates within and upstream of Kane Basin along closed arches up to the entrance of Robeson Chanel and eventually along open arches and more linear features upstream of the Strait (Fig. 6a, t = 20 and 72 hours). The opening of 5 multiple arch-like leads under the effect of the wind forcing appears clearly on the corresponding fields of ice concentration (Fig. 6b). Comparison of these fields and corresponding snapshots of the ice drift velocity over the domain indicates that these leads divide the ice into relatively undamaged plates, or floes. Drift speeds are piecewise constant over the floes and discontinuous across the floes (see Fig. 7, t = 6 hours), a feature that is also reproduced in the idealized simulations, as shown by the "stair-case" like profiles of u y (see Fig. 5a). This behaviour is consistent with the motion of the Arctic ice cover as 10 revealed by Synthetic Aperture Radar imagery analysis and RGPS motion products (Kwok, 2001;Moritz and Stern, 2001).
Also evident from the fields of ice drift speed are regions where the ice remains landfast even under strong wind forcing, either enclosed in bays and narrow fjords or between islands and the nearby shore (see Fig. 7, t = 72 hours), consistent with the observation of landfast ice along the coast of Kane Basin and Smith Sound in winter and spring (Mundy and Barber, 2001).
The sharp gradients in ice velocity between these motionless regions and the fast flowing ice are well simulated. 15 For C min ≥ 20 kPa, stable ice arches form in the realistic simulations. Consistent with observations, the main bridge is located at the constriction point between Kane Basin and Smith Sound (see Fig. 6b and 7, t = 72 hours). Downstream of this bridge, the ice concentration rapidly drops, as the ice detaches and is driven out the channel at the free drift speed, leading to the opening of the North Water polynya (see Fig. 6b, t = 72 hours). The model also simulates the detachment of the ice from the leeward side of islands (see Fig. 7, t = 72 hours) and, as the wind forcing is further increased, from a secondary bridge at 20 the entrance of Kane Basin (see Fig. 6b, t = 72 hours), which is also observed on satellite imagery (Kwok et al., 2010).
Internal stresses in the realistic simulations also evolve similarly to the idealized case. After the onset of damage, the simulated instantaneous states of stress, represented in the principal stress plane, appear in good agreement with in situ stress measurements from the Beaufort Sea reported by Weiss et al. (2007) (see Fig. 8). Consistent with these observations, biaxial tensile stresses occur in a large number both prior and after the formation of stable ice bridges in the Strait. This again supports 25 the relevance of accounting for some resistance of the ice in pure tension in sea ice models. Also consistent with in-situ measurements, large biaxial compressive stresses do not appear so frequent compared to pure tensile and biaxial tensilecompressive states (Weiss et al., 2007;Weiss and Schulson, 2009), although convergent ice motion occurs over a significant portion of the domain. It is 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/contraction (Richter-Menge et al., 2002), which 30 are not represented in the present Maxwell-EB framework. 16 The Cryosphere Discuss., doi: 10.5194/tc-2016-276, 2016 Manuscript under review for journal The Cryosphere Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License.

Consequences in terms of sea ice export
While studies have highlighted the strong interannual variability of the flux of ice through Nares Strait (Kwok et al., 2010), observations have also shown a tendency for ice bridges to form later and break earlier in later years than in the 1980s , suggesting than an increased ice outflow through the Strait is expected in a future climate. The year 2007 for instance was characterized by the absence of flow stoppage (Münchow, 2016), which resulted into a record areal and volume 5 outflow, equivalent to twice their average value over the 1997 to 2009 period (Kwok et al., 2010). The following question arises: is this tendency towards shorter lifespan of ice bridges and associate increased ice export through straits the consequence of a thinning and/or a mechanical weakening of the ice cover in recent years?
In a recent study, Gimbert et al. (2012b)  of the fields of |u| in the summer and winter cases moreover suggests that refreezing of leads may play a significant role not only in supporting stable ice bridges, but also in maintaining the ice landfast along portions of the coast.
It is important to note that these numerical experiments are by no means an attempt to determine the physically appropriate value of the cohesion of the ice cover in either of these scenarios, as the presence or not of a stable ice bridge also depends on the magnitude and form of the applied wind forcing as well as on the prescribed initial ice thickness. Instead, this exercise 5 serves to show that by varying few mechanical strength parameters in a range of values which, based on observations, seems physically appropriate, the Maxwell-EB model is able to reproduce widely different dynamical behaviours. Conversely, the simulations suggest that mechanical weakening, independent of thinning of the ice cover, can play a role in increased ice exports through narrow outflow pathways of the Arctic. 10 We now investigate the ice thickness distribution simulated by the Maxwell-EB model and, in particular, the distribution of ridged ice, i.e. the ice that has thickened through mechanical redistribution. In the simulations analyzed for this purpose, a stable ice bridge does not form: the flow of ice upstream the channel therefore slows but does not stop, such that the ice is more readily redistributed under the applied wind forcing.  1 m), which at all times is well described by a negative exponential (the coefficient of determination for the goodness of the fit varies between 96% and 98%). This function is of the form

Ice thickness distribution
with h * increasing in time until the convergence and thickening of ice at the entrance of the channel starts slowing significantly 25 the ice flow (see Fig. 10). The flattening of the PDF is a signature of the thickening of the ice along these highly localized features.
The part of the distribution with h < 1 is characterized by a second mode near h = 0 associated with the presence of open water (e.g., Wadhams, 1981Wadhams, , 1994Haas et al., 2006;Haas, 2009), which becomes more important as ice is driven out of the domain. At t = 3 days (orange curve), the ice downstream of the channel has nearly left the domain and a large patch of ice is the Galerkin discontinuous method coincides with a finite volume scheme with upwinding, which is known to be diffusive. As expected, numerical diffusion is most important at the edge of the detached ice, where mechanical stresses vanish, gradients between the ice and opening water are the strongest and drift velocities the highest. It is important to note however that the present continuum model is not meant to represent the dynamical behaviour of sea ice in regions of low ice concentration dominated by free drift conditions such as the marginal ice zone, but rather that of the ice pack. A Lagrangian model would 5 perhaps be more suitable to simulate the edge of the detached ice, although some diffusion would be unavoidable in free drift mode as some remeshing would be required. Where ice concentration is high A > 90%, mechanical stresses are significant and constantly redistributed under damaging. In the Maxwell-EB model, the associated deformation is highly localized in both space and time, which acts to mitigate numerical diffusion and re-increase gradients in all fields.
The realistic simulations show a similar evolution of the ice thickness in the channel. In the funnel-like entrance to Nares 10 Strait, the converging part of Kane Basin and upstream of coasts and islands, the ice builds-up along narrow and oriented features (see Fig. 11a). The effect of numerical diffusion in smoothing these features increases with ice drift velocities downstream of the domain. Nevertheless, at all times the simulated probability density function is strongly asymmetric, consistent with thickness distributions estimated for sea ice with little history of melting (e.g., Haas, 2009). As in the idealized case, the strong localization of the ridged ice translates into an exponential tail for P (h) of the form of (13) with h * increasing in time 15 (see Fig. 11b).
According to Eqn. (4) and to the simple redistribution scheme employed here, the evolution of h is a function of the flux ∇ · (hu) = u · ∇h + h · ∇u and of the mechanical redistribution term given by Eqn. (7), which itself is a function of the flux of ice concentration ∇ · (Au). The spatial distribution of h in the present Maxwell-EB model therefore depends essentially on the simulated velocity field. Furthermore, both simulations, idealized and realistic, show a similar exponential decrease for the 20 tail of P (h). This suggests that the localization of the thickest ice does not arise from the complexity of the domain geometry.
The shape of the tail of P (h) is also conserved when using a lower spatial resolution, for instance, with ∆x = 4 and 8 km in idealized simulations (orange, dotted and dotted-dashed curves on Fig. 11b), suggesting that this property of P (h) is not resolution-dependant either. The strong localization of ridged ice in the model therefore appears to be only the consequence of its capability to reproduce the extreme localization of ice deformation and associated sharp gradients in the ice velocity field. 25 6 Summary and conclusion In this paper we have presented the results of a first implementation of the Maxwell-EB rheology for modelling sea ice on geophysical scales. Idealized and realistic simulations of the flow of ice through Nares Strait have shown that the Maxwell-EB sea ice model is able to reproduce the formation of multiple arch-like leads within the ice cover downstream, in the interior and upstream of the channel, 30 in agreement with observations of ice conditions in narrow straits of the Arctic (Sodhi, 1977;Kwok et al., 2010). As these features appear in high numbers in both highly idealized and realistic simulations, they are not attributable to the
complexity of the domain geometry or boundary conditions but to the mechanical behaviour of the Maxwell-EB model itself. In coupled thermodynamic and dynamic models, a high density of leads is expected to impact the simulated heat fluxes between the atmosphere, the ice and the ocean (Smith et al., 1990).
the formation of a stable, concave ice bridge, the associated stoppage of the ice flow upstream and the opening of a polynya downstream of the bridge. In the realistic simulations of Nares Strait, the location of the main bridge is 5 consistent with the stable arch observed to form seasonally downstream of Kane Basin.
regions of relatively undamaged ice with uniform, plate-like motion, and extreme gradients in ice velocity across the leads that delimit these regions, consistent with RGPS observations (see Kwok (2001); Moritz and Stern (2001)). In particular, the model reproduces the very sharp gradients in ice velocity, concentration and thickness associated with the edge of the ice bridge while remaining numerically stable.

10
the presence of landfast ice along portions of the coast of Nares Strait, in agreement with observations (Mundy and Barber, 2001). In the model, this landfast ice resists the strong applied wind forcing. The associated sharp gradients in ice velocity, thickness and concentration between this stagnant and the nearby fast flowing ice are well represented.
states of stress that are overall in good agreement with in-situ measurements in Arctic sea ice. In particular, the simulations have shown that pure (i.e., uniaxial or biaxial) tensile states of stress are rather frequent, hence pointing out the 15 importance of accounting for some resistance of the ice in tension in sea ice models.
the strong localization of ridged ice and an associated thickness distribution with an exponential tail, in agreement with probability density functions calculated from sea ice thickness measurements.
This last point is an important outcome of this first implementation of the Maxwell-EB rheology in a realistic context, as the treatment of the mechanical redistribution of the ice thickness, in particular of the ridging process, has been the subject 20 of numerous studies since the 1970's and continues to challenge the sea ice modelling community to this day. Although the excessively simple redistribution scheme used here does not include multiple ice categories nor allows prescribing the thickness of the ice involved in the ridging, the Maxwell-EB model seems to allow ice to thicken in regions of strong convergence and shear, a process that is known to be underestimated in VP models using a two-level scheme and which has been more adequately simulated at the cost of increasing the number of ice thickness categories in VP (or EVP) models using a multi-25 categories redistribution scheme (e.g., Bitz et al., 2001). In the Maxwell-EB model, this capability of accounting for a sufficient thickening of the ice as well as the spatial localization of extreme thickness values arises from the appropriate description of extreme strain localization. On a mechanical point of view, this may therefore question the relevance of using multi-categories redistribution schemes.
In terms of the simulated thermodynamic fluxes and ice energy balance however, a number of studies have highlighted thin ice using a two-level versus a multi-categories model (e.g., Walsh et al., 1985). As opposed to mechanical redistribution processes, thermodynamic processes are expected to "seek the mean" (Thorndike et al., 1975) and smooth the distribution by allowing ice to grow over open water, thinner ice to thicken faster than thick ice and thicker ice to melt comparatively faster (Haas, 2009 Besides numerical efficiency, other advantages of using a simple redistribution scheme such as the one employed here is that no thickness redistribution function needs to be assumed and the redistribution is not directly tied to the prescribed failure strength of the ice. In the Maxwell-EB model, the assumed strength is instead based on in-situ stress measurements, which point to a Mohr-Coulomb failure criterion and directly provide information on the order of magnitude of the shear 10 strength. In particular, prescribing a cut-off for biaxial compressive strength (equivalent to the pressure, P , in VP models) appears unnecessary. In the model, the strength of the ice is set through a single parameter, the cohesion C. The channel flow simulations performed here showed that varying this parameter in the limit of physically acceptable values and in agreement with the suggested evolution of the failure strength of the ice pack in later years (Gimbert et al., 2012a), while keeping the initial thickness unchanged, is sufficient to reproduce different dynamical behaviours of the ice cover (the presence of stable

J Weiss and E M Schulson
Arctic sea ice cover over a 3-day 97), obtained at a scale of llite images s/radarsat.html). (a) Shear km 2 cells with a shear strain-rate ed, revealing linear faults where ugh these cells represent only ta, they concentrate 54% of the r strain-rate below this threshold, rate (in 1 day −1 ): the 10 × 10 km 2 han 0.02 day −1 are indicated, in (a). Although these cells ered by the data, they concentrate cells are plotted in uniform grey. g scale of measurement (see t al (2008)) and section 3.3). ns for Coulombic faulting are ale-independent friction and in section 3.1 offer strong e-independent fracture physics (a) (b) 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  at different times (coloured lines). A bin width of 0.1 meters is used. Note the semi-logarithmic axis. The tail of the distribution (h ≥ 1.1 meters) is fitted using a linear regression to estimate the slope, − 1 h * . The percentage of the variance (R) explained by this fit is > 95%. The insert shows the values of h * calculated at different times using a linear regression of ln P (h) for h ≥ 1.1 m.