neXtSIM: a new Lagrangian sea ice model

. The Arctic sea ice cover has changed drastically over the last decades. Associated with these changes is a shift in dynamical regime seen by an increase of extreme fracturing events and an acceleration of sea ice drift. The highly non-linear dynamical response of sea ice to external forcing makes modelling these changes and the future evolution of Arctic sea ice a challenge for current models. It is, however, increasingly important that this challenge be better met, both because of the important role of sea ice in the climate system and because of the steady increase of industrial opera-tions in the Arctic. In this paper we present a new dynam-ical/thermodynamical sea ice model called neXtSIM that is designed to address this challenge. neXtSIM is a continuous and fully Lagrangian model, whose momentum equation is discretised with the ﬁnite-element method. In this model, sea ice physics are driven by the combination of two core components: a model for sea ice dynamics built on a mechanical framework using an elasto-brittle rheology, and a model for sea ice thermodynamics providing damage healing for the mechanical framework. The evaluation of the model performance for the Arctic is presented for the period September 2007 to October 2008 and shows that observed multiscale statistical properties of sea ice drift and deformation are well captured as well as the seasonal cycles of ice volume, area, and extent. These results show that neXtSIM is an appropriate tool for simulating sea ice over a wide range of spatial and temporal scales.


Introduction
Sea ice dynamics are very complex and share many characteristics with earth crust dynamics, such as dynamical triggering and clustering of deformation events or earth/ice quakes. Both sea ice and the earth's crust are geophysical solids that can be viewed from the mechanical point of view as two-dimensional plates due to their very small geometrical aspect ratio (O(10 −6 )). These plates then experience planar internal stresses under the action of winds and ocean currents in the case of sea ice and magmatic currents of the mantle in case of the earth's crust. Similarly to earth crust dynamics, sea ice dynamics are controlled by processes interacting and evolving over a wide range of spatial and temporal scales. Mechanical processes like fracturing and faulting are important as they both drive large-scale sea ice drift and deformation patterns (e.g. . These processes are the expression of the mechanical damage of the ice pack which, as a result, may look more like an assembly of plates (> O(1 km)) and floes (< O(100 m)) than an intact solid plate. In addition to the damaging processes the formation of new ice is also important. Indeed, new ice formed in fractures and leads can fuse together broken ice and thus contribute to an effective mechanical strength recovery, or "healing". The observed complex dynamical behaviour of the sea ice cover therefore emerges from the interplay of these dynamical and thermodynamical processes (see for example Korsnes et al., 2004). As an example of this complexity, recent studies showed that the statistical properties of sea ice deformation are characterised by a coupled space-time scaling invariance (Marsan et al., 2004;Rampal et al., 2008), similar to what is observed for earthquakes (e.g. Kagan, 1991;Kagan and Jackson, 1991;Marsan and Weiss, 2010), and which is a fingerprint of the presence of long-range elastic interactions within P. Rampal et al.: neXtSIM: a new Lagrangian sea ice model the ice cover. In this paper we present a new general sea ice model, called neXtSIM, which has been recently developed and designed to correctly capture this complexity. Coon et al. (1974) introduced the first realistic dynamical sea ice model, based on observations from the Arctic Ice Dynamics Joint Experiment (AIDJEX). In their model, sea ice was described primarily as a plastic material that deforms irreversibly once a critical internal stress state is reached. When the stress is subcritical, however, the ice is modelled as an elastic material that deforms under stress, but returns back to its original shape when the stress is removed. Hibler (1979) replaced the elastic response of the AIDJEX model by a viscous one, producing the viscous-plastic model (VP). This made his model easier to solve numerically and easier to couple to ocean general circulation models. Hunke and Dukowicz (1997) suggested adding an elastic term to the VP model of Hibler (1979), producing the elastic-viscousplastic model (EVP). This modification was based purely on numerical considerations, making the model easier to parallelise, but offers no additional physical insights. Virtually all modern sea ice models use either the VP or EVP formulation, combined with a thermodynamics model (e.g. Semtner, 1976;Bitz and Lipscomb, 1999;Vancoppenolle et al., 2009;Turner and Hunke, 2015) and variously detailed subgrid-scale parameterisations (for commonly used large-scale sea ice models as for instance CICE , LIM3 (Vancoppenolle et al., 2009), MITgcm (Adcroft et al., 2016) or MPI-ESM (Notz et al., 2013)).
One of the reasons to redesign or replace the VP and EVP rheologies is that classical models give a poor representation of ice drift and deformation statistics and scaling, compared with satellite observations (Girard et al., 2009). Girard et al. (2011) introduced the elasto-brittle rheology and showed that this has the potential to accurately reproduce the aforementioned statistics and scaling. Bouillon and Rampal (2015b) then introduced the dynamical core of the new sea ice model presented in this paper, using the elasto-brittle rheology. In their paper they described a preliminary stand-alone version of the model used to simulate the sea-ice-damaging process and the associated ice cover deformation over short timescales (up to 10 days), while neglecting the thermodynamical processes and feedbacks (e.g. on the sea ice mechanical strength). This sea ice model was capable of reproducing one of the complex aforementioned characteristics of sea ice dynamics, i.e. the multifractal spatial scaling of sea ice deformation, revealed by satellite observations analysis and reported for the first time in Marsan et al. (2004).
The main goal of the neXtSIM development is to reproduce the mechanical behaviour and state of the Arctic sea ice cover on seasonal timescales (over 1 year in this paper). In particular, we wish to simulate realistic sea ice drift and deformation statistics and annual cycle, as well as sea ice volume and extent seasonal cycles. Addressing such a temporal scale required some developments from the first simplified version of neXtSIM presented in Bouillon and Rampal (2015b), such as an adapted rheological framework to take care of post-damage sea ice motion and permanent deformation, and a thermodynamical model capable of producing the necessary feedback on the sea ice mechanical behaviour over a seasonal timescale.
This paper presents the first comprehensive version of the neXtSIM model, a fully Lagrangian dynamical/thermodynamical sea ice model. A generic presentation of the model is made in Sect. 2. The remeshing/remapping scheme is described in Sect. 2.4. The model set-up, tuning of parameters, and evaluation of the model performance are described in Sect. 3. Section 4 presents the summary and conclusion.

Model description
The sea ice variables used in neXtSIM are the following: h and h s are the effective sea ice and snow thickness respectively; A is the sea ice concentration; d is the sea ice damage; u is the horizontal sea ice velocity vector; and σ is the ice internal stress tensor. These variables are listed in Table 1. The model has two ice thickness categories: ice and open water. As in Bouillon and Rampal (2015b), scalar and tensorial variables are defined at the centre of the elements of the mesh, whereas vectors are defined at the vertices.

Evolution equations
The evolution equations for h, h s , and A (here denoted φ) have the following generic form: where Dφ Dt is the material derivative of φ, ∇ · u is the divergence of the horizontal velocity, and S φ is a thermodynamical sink/source term. The actual form of S h , S A , and S h s are defined in Sect. 2.3. An additional constraint for the concentration is that A ≤ 1. The evolution of d is given by where d is the damage source term (defined below) and S d is a thermodynamical sink term (defined in Sect. 2.3). The damage is an ice volume tracer and is equal to 0 for newly formed ice. When new ice is formed, the new damage is decreased and calculated as a volume-weighted average over the old and new ice, meaning that the sea ice partly recovers its mechanical strength. The evolution of sea ice velocity derives from the vertically integrated sea ice momentum equation: where m is the inertial mass, P is a pressure term, τ a is the surface wind (air) stress, τ w is the ocean (water) stress and τ b is the basal stress in case of grounded ice. All these terms are defined in Sect. 2.2. The other symbols in Eq.
(3) are f , the Coriolis parameter; k, the upward pointing unit vector; g, the gravity acceleration; η, the ocean surface elevation; and σ , the internal stress tensor. The evolution of the internal stress is computed as in Bouillon and Rampal (2015b) in two steps that correspond to where superscripts n and n + 1 correspond to the previous and current time steps, respectively. The first step accounts for the elastic deformation without considering the damaging process and gives a first estimate of the internal stress, σ , by where E(A, d) is the effective elastic stiffness, ν is Poisson's ratio, and˙ is the strain rate tensor defined by˙ = 1 2 ∇u + (∇u) T . The effective elastic stiffness is parameterised as where Y is the sea ice elastic modulus (Young's modulus) and f (A) is a decreasing function of A parameterising the effect of the compactness (defined in Sect. 2.2). The second step accounts for the damaging process. With this estimate of the internal stress, the failure criteria are checked. For the elements where the estimated internal stress σ falls outside the failure envelope (defined in Sect. 2.2), the damage factor is set to the value for which the stress state, is scaled back onto the failure envelope. For the elements for which the estimated internal stress σ is inside the failure envelope, is simply set to 1. The damage source term d corresponding to the decrease of σ has been derived in Bouillon and Rampal (2015b) as

Dynamical component
In this subsection, we detail each term of the sea ice momentum equation (Eq. 3) and the underlying parameterisations. The inertial mass m depends on the assumption made about the motion of the water present in leads: either the water in the leads moves as the ocean below or as the ice around it (Connolley et al., 2004). We choose the second hypothesis for our model, as we think it is more relevant for high-resolution models and when the rheology allows for sharp transitions of sea ice concentration within the ice cover. Using this approach, the inertial mass m corresponds to the mass of ice and snow plus the mass of the water in the leads as m = ρ i h + ρ s h s + ρ w h w , where ρ w is the reference density of seawater and h w is the volume of water from the base of the ice to the sea surface per unit area. By isostasy, , and m is then given by The failure envelope is defined as in Weiss et al. (2007)  respectively, following the convention that compressive stresses are positive. The envelope represents a combination of a Mohr-Coulomb criterion, a tensile stress criterion, and a compressive stress criterion, defined by where q and σ c are defined as in Weiss and Schulson (2009) σ c = 2c where µ is the internal friction coefficient and c is the cohesion. σ N min < 0 and σ N max > 0 are the maximal tensile stress and the maximal compressive stress, respectively. The friction coefficient µ for sea ice is chosen equal to 0.7, which is a common value for geo-materials (Amitrano et al., 1999) and seems to be scale-independent (Weiss and Schulson, 2009). The values of the cohesion c, σ N min , and σ N max seem to be inversely proportional to the square root of the spatial scale (Weiss et al., 2007). Here we use this scaling relationship to define their values at 10 km as equal to c = 8 kPa, σ N min = −9.52 kPa, and σ N max = 150 kPa from values measured in the field (Weiss et al., 2007) and in lab experiments (Schulson, 2009). P is a vertically integrated sea ice pressure term that is set to avoid excessive convergence when the ice concentration in a cell is at 100 % and at the same time highly damaged. Without this term, one obtains unrealistically large local thickness of the ice cover, for example north of Greenland and the Canadian Archipelago. This term implies no memory effect, meaning that it cannot be included in the evolution equation of the internal stress. This term is parameterised as where P * is the pressure parameter, f (A) is the same function as in Eq. (6), and determines if the pressure term is active or not. In our case, is defined as a function of the divergence rate at the previous time step and is computed as where min is a parameter set to a small value to regularise the transition when the divergence rate is close to 0. The quadratic dependence on the mean thickness and the value of P * used in this study (P * = 12 kPa) comes from Hibler (1986) and corresponds to the redistribution scheme of Thorndike et al. (1975) when it is applied to only one ice thickness category.
The effect of the concentration on the mechanical response of sea ice is parameterised here by a decreasing exponential function of the concentration: where α is the compactness parameter (see Bouillon and Rampal, 2015b, for more details). The air stress τ a is computed following the quadratic expression: where u a is the wind velocity, ρ a the air density, θ a the atmospheric turning angle, and c a the atmospheric drag coefficient.
The oceanic stress τ w is also computed following a quadratic expression, namely where u w is the ocean velocity, ρ w the reference density of seawater, θ w the water turning angle, and c w the water drag coefficient.
The basal stress τ b is a term added to simulate grounded fast ice. It is parameterised as in Lemieux et al. (2015) by the expression: where k 2 is the maximum basal stress parameter, u 0 is the basal stress velocity parameter, and C b is the basal stress concentration parameter. The critical thickness from which the parameterisation starts acting is defined as h c = A H +η k 1 , where k 1 is the critical thickness parameter and H is the ocean depth.

Thermodynamical component
In neXtSIM damaged sea ice recovers its mechanical strength (i.e., decrease of the damage) through time via two processes: formation of new ice in open water and leads and thermodynamical healing. Sea ice melting is supposed to have no direct impact on the damage. New ice formation is naturally treated by updating the value of the local damage as a volume-weighted average over the old and new ice. When sea ice volume in a cell increases due to ice formation, the damage then automatically decreases as new ice is supposed to have a damage equal to zero. For the thermodynamical healing process, more assumptions need to be established. We assume here that the thermodynamical healing process is driven by the local temperature gradient between the bottom of the ice and the snow-ice interfaces and decreases the effective compliance, defined as C = 1 1−d , at a constant rate. This is based on the fact that the cooler the environment, the faster the ice will freeze, so presumably low temperatures result in fast healing and warm temperatures in slower healing, with no healing occurring for temperatures over the freezing point. The damage relaxation term S d of Eq. (2) is then computed by with C (n+1) given by and where T C is the compliance relaxation time. This relaxation time is assumed here to be inversely proportional to the temperature difference T between the bottom and snowice interface, which is given by where T b is the temperature at the ice base, T s the temperature at the ice or snow surface, k i is the heat conductivity of ice, and k s is the heat conductivity of snow. The compliance relaxation time is then defined as where T d is the parameter controlling the healing rate. We use the constant 1000 and 40 K, which are typical values of effective compliance and temperature difference given by the model during winter, so that T d can be interpreted as the time needed to heal the ice in winter conditions. The sensitivity to this parameter is discussed in Sect. 3. We also limit T C to be positive so that melting conditions alone (T s > T b ) cannot damage the ice. The other components of the thermodynamical model are similar to those in classical sea ice models. There are three thermodynamical source and sink terms corresponding to S φ in Eq. (1), one for each of ice volume, concentration, and snow volume. The source/sink term for the ice volume stems from the conservation of mass and can be written as where h is the change in level ice volume and h ow is ice formation in open water. The change in A is calculated by assuming a given thickness for the ice forming over open water ( h ow ). We use a constant, h 0 , for this thickness, giving a source/sink term for A as Ice formation in the open water portion of the grid cell, h ow , is calculated such that heat loss from the ocean that would cause supercooling is redirected to ice formation. This is an adaptation of the form suggested by Hibler (1979) in that he uses prescribed growth rates, but we calculate those depending on the atmosphere and ocean states, as described below. The source/sink term for snow thickness, S h s , also stems from the conservation of mass and is The first term accounts for snowfall and snowmelt, the second term removes snow due to the lateral melt of ice, and the third term converts snow into ice whenever the ice freeboard, D, falls below the water surface due to snow loading. Energy needed to melt snow due to the lateral melt of ice is removed from the ocean as an additional heat flux. Thickness changes in level ice and snow are calculated using the zero-layer model of Semtner (1976), using the same parameter values, unless otherwise stated. This is arguably the simplest usable thermodynamic sea ice model, and it has some well-known deficiencies (most notably Semtner, 1984). It does, however, suffice for short runs with a stand-alone ice model, like the ones discussed in Sect. 3.4. In this model the incoming radiative fluxes are interpolated from the forcing data, applying constant albedos to the incoming short-wave radiation of α i = 0.64 for the ice and α s = 0.85 for the snow (Maykut and Untersteiner, 1971). The turbulent heat fluxes are calculated using bulk formula for the sensible heat flux: and the latent heat flux: Here, c p is the atmospheric heat capacity and L v is the latent heat of sublimation. The temperature difference, T a − T s is taken between the ice surface and the atmosphere at 2 m, the same as the difference in specific humidity, q a − q s . We calculate the specific humidity using the formulation of Buck (1981). The drag coefficients C t and C q are set to 1.3 × 10 −3 based on Ebert and Curry (1993). Fluxes between the ice and ocean, Q oi , are calculated assuming the ocean underneath the ice must always be at the freezing point. In order to produce realistic heat fluxes through the ice, the thermodynamical ice model must be coupled to an ocean model. Here we use a simple slab ocean model that consists of a single ocean layer with one temperature and salinity www.the-cryosphere.net/10/1055/2016/ The Cryosphere, 10, 1055-1073, 2016 point per grid cell. The flux of heat at the ocean surface is calculated as a weighted average of the ocean-ice and oceanatmosphere fluxes: where Q oa is the ocean-atmosphere flux. The oceanatmosphere heat flux is calculated using the bulk formulas (Eqs. 30 and 31), with C t = 0.83×10 −3 and C q = 1.5×10 −3 (Gill, 1982) for the turbulent heat fluxes, while the radiative heat fluxes are read in from the atmospheric forcing, applying a constant albedo of 0.07 to the short-wave flux. The change in ocean salinity is calculated assuming the total salt content of the ice-ocean system is conserved, resulting in a change in salinity of for a change in ice and snow area mean thickness of h and h s respectively, and where S o and S i are the ocean salinity and ice salinity, respectively, ρ i and ρ s are the ice and snow density, respectively, F fw is freshwater flux at the surface, H ml the mixed layer depth, and t the model time step. We assume a constant ice salinity of S i = 5 psu.
When using a slab ocean, simulated ocean temperature and salinity have to be artificially maintained at realistic values because of the missing representation of both vertical and horizontal heat and salt exchanges within the ocean. Here, we use Newtonian nudging to relax the simulated ocean temperature and salinity towards the values of the uppermost ocean layer from a full ocean model (in this case the TOPAZ4 system, see Sect. 3.1). The local mixed layer depth from the full ocean model is also used as the mixed layer depth (H ml ) of the slab ocean model.

Remeshing and remapping
Most sea ice models use an Eulerian approach for the advection. However, we believe that a purely Lagrangian approach as in Wang and Ikeda (2004) may be more suitable to preserve highly localised features (i.e. one cell wide ridged or open water areas) generated by the model. Continuous, purely Lagrangian schemes require unstructured meshes and a procedure for mesh adaptation. Local mesh modifications can be done in parallel and introduce very low numerical dissipation (Compère et al., 2009). They also show local conservation (Compère et al., 2008).
In the purely Lagrangian approach, the vertices of the element (i.e. the nodes of the grid) move with the sea ice velocity u. The material derivative is then simply equal to the temporal derivative ∂φ ∂t X relative to the moving mesh so that the quantities are naturally transported with the ice. The sea ice thickness and concentration, for example, are simply updated by and where S (n) and S (n+1) are the surface of the element at time steps n and n + 1. The variables defined at the nodes do not need to be updated.
In this approach the model mesh deforms as the ice cover itself deforms. When the mesh becomes too distorted the results of the finite element method are no longer reliable and the mesh must be adjusted, a process referred to as remeshing. In the current implementation the mesh is considered too deformed when the smallest angle of any triangle of the mesh is smaller than 10 • . Using this criterion and the set up we use here, the mesh needs to be adapted on average every model hour.
In order to save computational time the forcing fields are only interpolated onto the model grid after remeshing, or when new forcing fields are required. This means that even though the model grid drifts and deforms in-between remeshings the forcing seen by the nodes and elements of the model does not change. We checked that this method gives virtually identical results as when we interpolate the forcing fields every time step. Indeed as the remeshing criterion is global the error in the position of the forcing field is in practice never larger than a single model element. Given the high resolution of the model grid in our tests, the forcing fields are too smooth and too coarsely resolved for this error to have any substantial effect.
The approach used for the slab-ocean is similar to that used for the forcing in that the fields are only interpolated after remeshing. The slab-ocean model resides on the same mesh as the ice model, but the relative displacement of ice and ocean is ignored between remeshing steps. When the model mesh becomes too deformed and therefore needs to be remeshed, the temperature and salinity are interpolated from the old onto the new mesh using a linear interpolation and ignoring the displacement of the old mesh. This ensures that the temperature and salinity fields do not drift with the ice as the ice-model mesh moves.
The new mesh is created by a version of the BAMG mesh generator by Hecht (1998) taken from the Ice Sheet System Model (Larour et al., 2012). This mesh generator can be instructed to preserve as many of the nodes from the old mesh as possible. The mesh is thus only modified in a limited number of locations, hereafter called cavities, at each remeshing. Doing this allows the model to track large expanses of drifting ice that is deforming very little without any artificial diffusion, since it is only necessary to interpolate values from the new grid to the old one inside the cavities. Outside the cavities the tracer values are not affected by the remeshing. For the variables defined at the nodes (i.e. the sea ice velocities, etc.), a non-conservative linear interpolation is performed for the new nodes.
For the quantities that are defined at the centre of the elements, a conservative remapping scheme is applied to each cavity independently. The cavities are defined as the smallest partitions of the mesh for which the external edges are the same before and after the remeshing. We implemented an algorithm that uses the information provided by BAMG (i.e. node-element connectivity and previous numbering of the preserved nodes) to efficiently detect the cavities and the intersections between the triangles of the old and new meshes. For each intersection, the corresponding integrated quantities are transferred from the old mesh to the new one. The process is fully conservative and generates only limited numerical diffusion.

Simulation set-up
For the model evaluation we force our model with the ocean state of the TOPAZ4 reanalysis (see Sakov et al., 2012), and with the atmospheric state of the Arctic System Reanalysis, Interim version (ASR-Interim hereafter) (http://rda.ucar. edu/datasets/ds631.4/, Byrd Polar Research Centre/The Ohio State University (2012). Accessed 1 January 2014.) TOPAZ4 is a coupled ocean-sea ice system combined with a stateof-art ensemble Kalman filter data assimilation scheme of both ocean and sea ice variables, running at an average resolution of 12.5 km in the Arctic. The TOPAZ4 bathymetry is based on the 1 arc minute GEBCO bathymetry (Jakobsson et al., 2008) and the coastline is derived from the 5 m isobath. The main benefit of using the TOPAZ reanalysis in this context is its accuracy in simulating the location of the ice edge, and therefore to provide realistic forcing for ocean temperature and salinity. The ASR-Interim is a high-resolution atmospheric reanalysis (30 km) known to reproduce the nearsurface wind fields particularly well (Bromwich et al., 2015).
In order to simplify the forcing of the slab ocean with TOPAZ4, the domain of our model is defined from TOPAZ4 coastlines and open boundaries. The resulting mesh covers the Arctic and North Atlantic oceans, extending from an open boundary at 43 • N in the North Atlantic to an open boundary in the Bering Strait (see Fig. 1). The resolution of the finite element mesh is about 11 km, where the resolution is defined as the square root of the mean element area. Note that the ocean depth H used for the basal stress parameterization comes from the 1 arc minute ETOPO1 global topography (Amante and Eakins, 2009).
The oceanic forcing variables are sea surface height, velocity at 30 m depth, and sea surface temperature and salinity, which are provided as daily means by the TOPAZ4 system. We interpolate these quantities temporally and spatially onto the model mesh at run time using linear and bilinear interpolation methods, respectively. The slab-ocean temperature and salinity are nudged towards TOPAZ4 and we found 30 days to be an appropriate nudging timescale for this set-up. This value allows our model to reproduce the location of the ice edge well without unduly affecting heat fluxes in the central Arctic. The heat flux resulting from the nudging is usually slightly below 0.5 W m −2 in the central Arctic in midwinter, which is a reasonable value (see e.g. Sirevaag et al., 2011).
The atmospheric forcing consists of the 10 m wind velocity, the 2 m temperature and mixing ratio, mean sea level pressure, total precipitation and the fraction of that which is snow, and the incoming short-wave and long-wave radiation. All of these quantities are provided as 3-hourly means by the ASR-Interim. Similarly to the oceanic variables, we interpolate them temporally and spatially onto the model grid at run time using a linear and bilinear interpolation method, respectively.
In order to prevent an initial shock to the system when the model is started, the strength of applied wind and ocean currents is increased linearly from zero to full strength over the period of 1 day.
Our reference simulation starts on 15 September 2007 and ends on 9 October 2008. We choose this year because this is the only winter when the GlobICE (http://www. globice.info) and RGPS (Kwok et al., 1990) data sets overlap. The values of the model parameters that are used for this simulation are listed in Table 2. Damage is initially set to zero where sea ice is present. Initial sea ice concentration and thickness are interpolated from the sea ice of the TOPAZ4 reanalysis. The modelled ice thickness in TOPAZ4 is known to be unrealistically low on average compared to other Arctic ice-ocean coupled systems (Sakov www.the-cryosphere.net/10/1055/2016/ The Cryosphere, 10, 1055-1073, 2016  Warren et al. (1999) climatology, but we limit the snow thickness to 20 % of the ice thickness so that we do not get unrealistically thick snow on the relatively thin ice that was present at the start of the simulation.

Drag coefficient optimisation and evaluation of simulated ice drift
In neXtSIM, as in most classical sea ice models, the air and water drags depend on four parameters: c a the air drag coefficient, θ a the air turning angle, c w the water drag coefficient, and θ w the water turning angle. The value of these four parameters has to be calibrated depending on which atmospheric and oceanic forcing are being used. The purpose of this section is to present how we proceed with this calibration for the present study.

Basics of the method
By performing a scale analysis, it can be shown that the sea ice momentum equation (Eq. 3) is actually dominated by three terms: the ice internal stress, the surface wind drag, and the surface ocean drag. This equation can therefore be written as To prevent the rheology from affecting the optimisation process of the drag parameters, we only consider situations for which sea ice is in "free drift", i.e. situations where the internal stress term can be neglected in Eq. (36). By using Eqs. (20) and (21) and assuming that |u| |u a |, the solution of Eq. (36) then becomes where Na = √ ρ a c a /ρ w c w is the Nansen number. The first estimate of this number (Na ≈ 2 %) was made by Fridtjof Nansen during the Fram expedition (1893-1896) by comparing the drift of his boat, while trapped in sea ice, to local wind and ocean velocities. The air and water density being considered as constant, the Nansen number only depends on the ratio between the two drag coefficients, c a and c w . In the free drift mode, calibrating the Nansen number is then equivalent to calibrating one of the drag parameters while keeping the other one constant.
The calibration method uses a full winter of sea ice drift data (here from the GlobICE data set, http://www.globice. info), a reference run of 1 year starting in September 2007, and a series of short simulations restarted every 9 days from the reference run with the model being in free drift mode. To perform the simulation in free drift mode, we set the Young's modulus Y and the pressure parameter P * to 0. For each drift vector from the observation data set, we compute the corresponding simulated drift vector from the 6-hourly Lagrangian sea ice displacement fields produced by the two sets of experiments, the reference run and the 9-day free drift run. The simulated drift from the reference run is selected for the optimisation analysis only if it differs by less than 10 % from the drift simulated by the free drift run. As in McPhee (1979), we also restrict the analysis to the range of ice speeds going from 7 to 19 km day −1 . As a result, about 15 000 vectors are selected from the 20 analysed periods of 9 days (from 31 October 2007 to 28 April 2008). As shown on Fig. 2, the identified free drift events are mainly located in the Transpolar Drift, in the Beaufort Gyre, and near the ice edge (i.e. in Greenland, Barents, and Kara seas). Note that the number of identified free drift events also depends on the observation coverage, which is indeed high in the areas just mentioned.

Results of the method
By optimising the air drag parameters for these selected free drift vectors, we find an optimal value of c a = 0.0076, corresponding to a Nansen number equal to Na = 4.2 %. Here, the optimal value is found to be higher than the classical values, which is consistent with the negative bias documented for ASR-Interim surface winds (Bromwich et al., 2015). Doing the same exercise for ERA-Interim winds (Dee et al., 2011), which are frequently used to force large-scale sea ice models, the optimal air drag coefficient is found to be c a = 0.0023 (i.e. Na = 2.3 %), which is much closer to classical values. The scatter plots for each component of the selected vectors are shown in Fig. 3. For each component, the correlations between the simulated and observed free drift vectors are equal to 0.94 and 0.92, respectively, and the root mean square errors (RMSEs) are equal to 2.5 and 2.3 km day −1 . The cumulative probability of the error in velocity is shown in Fig. 3, along with the median and mean error, which are equal to 2.3 and 2.8 km day −1 , respectively. Note that the RMSE and median and mean errors, when using ERA-Interim with its optimal air drag coefficient, were found to be about 50 % larger than the ones found when using ASR-Interim.

Evaluation of simulated sea ice drift and sensitivity to the healing timescale
The quality of the simulated sea ice drift is evaluated by comparing simulated velocity vectors to the ones provided by the RGPS and GlobICE data sets between 31 October 2007 and 28 April 2008. The high spatial and temporal resolution of these data sets (about 3 days and 10 km) make it possible to compare a very large number of simulated and observed drift vectors, as shown on the scatter plots in Fig. 4, and in-crease the robustness and level of confidence of our model evaluation. The correlation for each component is slightly lower than in the free drift analysis (0.85 and 0.87, respectively, compared to 0.92 and 0.94) but is still very good.
The RMSE values are similar to these of the free drift analysis (2.5 and 2.2 km day −1 , respectively) which is remarkable knowing that here no selection nor restriction has been applied to the data. The median and the mean velocity errors are remarkably low, 1.9 and 2.4 km day −1 respectively. These results are in good agreement with observations, which may be attributed to the quality of the atmospheric forcings and to a proper representation of sea ice drift in the pack. The sensitivity of the correlation, RMSE, and velocity errors to the healing timescale parameter is presented on Fig. 5. The reference simulation is obtained with a healing timescale parameter equal to 28 days. Using larger healing timescales does not improve the correlation with either the observations, the RMSE, or the mean and median errors. Using shorter healing timescales decreases the skills of the model at reproducing the observed ice velocity.
To identify potential systematic errors, we also look at the mean sea ice drift by averaging modelled and observed drift over the whole season on a mesh grid of 21 by 21 km, covering the whole observation domain (see Fig. 6). For each cell, we also show the total number of observations to indicate the areas where data coverage is poor. As shown on Fig. 7, the largest differences between the observed and simulated mean ice drift are located in the Beaufort Gyre and Fram Strait and in some areas of the Kara, East Siberian, and Chukchi seas. These systematic errors may be partly explained by errors in the oceanic surface currents of TOPAZ, especially for the Beaufort Gyre. In the rest of the domain, the error on the mean winter drift is remarkably low, i.e. < 1 km day −1 .
The mean drift speed, taken over the central Arctic, correlates closely with the mean wind speed taken over the same area. This is to be expected, since the wind is the main driver of ice drift. We do, however, expect to see a significant difference between the ice response to wind in summer and in winter, due mainly to changes in ice concentration and thickness. In order to assess this effect, Fig. 8 shows the ratio  of drift speed to wind speed for the reference run, a model run forced with ERA-Interim, and the ratio of the climatology of the International Arctic Buoy Programm (IABP) buoy drift speed to ERA-Interim wind speed climatology. The two model runs considered here show a clear seasonal cycle in the drift speed to wind speed ratio. This is highest in summer, decreasing steadily from August to January, when it plateaus until April and then starts increasing again. Ólason and Notz (2014) found that over the 33 years they considered, the observed drift speed depends on concentration when concentration is low, thickness when concentration is high, and is related to an increased number of active fractures in April-May. It is not clear how well our model captures this relationship since the results for only 1 year can be heavily influenced by the timing and intensity of storms passing through the region. However, it is clear that the general shape of the observed time series for the drift speed to wind speed ratio is reproduced by the model, indicating that it captures the transition between freely drifting ice and pack ice correctly. In terms of magnitude, the drift speed to wind speed ratio for the run forced with ERA-Interim is slightly higher than the climatology. This is to be expected since both are based on ERA-Interim wind; the slight magnitude shift between the two is likely to be caused by the positive trend in Arctic sea ice drift speed that was originally revealed from the analysis of the IABP data set and reported by Rampal et al. (2009). Additionally, one can note a significant difference between the ratio time series of the reference run (in cyan) and that of the run forced with ERA-Interim (in blue). This can be explained by the fact that the winds in ASR-Interim (used as forcing in the reference run) are weaker than in ERA-Interim, but this effect is counteracted in the model by tuning the drag coefficient, as discussed earlier.

Evaluation of simulated sea ice deformation and sensitivity to the healing timescale
One of the main differences of neXtSIM compared to other sea ice models is the rheology, which defines the link between internal stress and deformation. For the internal stress, only a few observations are available and cannot be directly used for a complete evaluation. However, since we have calibrated the two other dominant terms of the momentum equation (i.e. the oceanic and atmospheric drag terms), then we can use an evaluation of the overall drift and deformation of the ice as an evaluation of the internal tress. A good way to evaluate the new rheology is then to compare the simulated deformation fields to the large amount of data available from satellite products. The data The Cryosphere, 10, 1055-1073, 2016 www.the-cryosphere.net/10/1055/2016/ used here are produced from the RGPS sea ice drift data set with the method proposed by Bouillon and Rampal (2015a). An interesting specificity of sea ice deformation is its strong localisation in space (see Fig. 9) and in time.
This makes a comparison of the geographical location of the observed and simulated deformation features impractical, since small errors in the applied forcing are bound to result in significant changes in the simulated location, compared to the extent of the features. Instead we compare simulated and observed deformation in a statistical sense using, among others, the multi-scale metrics introduced by Marsan et al. (2004) and Rampal et al. (2008). Ratio (%) Figure 8. Ratio of drift speed over wind speed for the reference simulation forced with ASR-Interim (cyan) and a simulation forced with ERA-Interim (blue). As a reference the same ratio is shown for the IABP buoys drift speed climatology over the ERA-Interim wind speed climatology (green).
The comparison with observation is focussed on the period January-April 2008, which has been identified as the period for which deformation is typically lower than during the rest of the year (see the annual cycle presented later in the paper). Figure 10 shows time series of the observed and simulated mean shear rate from January to April 2008. The deformation rates are computed at a spatial scale of 20 km and for 13 periods of 9 days. As in the previous sea ice drift evaluation, the model data are built to match the observations spatially and temporally. The correlation between the observed and simulated mean shear values is satisfactory, but we note that the model systematically underestimates the mean shear rate during this period.
Spatial scaling properties of sea ice deformation (or the degree of heterogeneity of sea ice deformation) can be studied from the analysis of Lagrangian trajectories, as e.g. in Marsan et al. (2004) who applied it to the trajectories of the RGPS data set. Here we perform this analysis for both the model and the satellite observation, following the method used in Bouillon and Rampal (2015a), which is very simi- lar to the one of Marsan et al. (2004), and which also gives an estimate of the error on the spatial scaling exponent. The spatial scaling analysis has been first applied to all the 13 snapshots corresponding to the 9-day periods between 1 January and 28 April 2008. Four snapshots had to be discarded because the power-law model fit was not significant in the least squared sense due to the excessive noise in the observed deformation fields. The values of the first-order moment of the shear rate (here denoted˙ s ) distribution obtained when gathering the selected snapshots is shown in Fig. 11 for different spatial scales. The power law ˙ s ∼ L −β that fits the data (grey and black lines) corresponds to a scaling exponent β of −0.16 for the observations and −0.11 for the model. The departure from the power-law fit at L < 20 km and L > 500 km comes from finite size effects (model resolutions and size of the Arctic basin, respectively). The scaling exponents β(q) for the other moment orders of the shear distribution fit well with a quadratic function of q, whose curvature is equal to 0.13 for the observation and 0.07 for the model. This shows that like observed deformations (Marsan et al., 2004;Rampal et al., 2008), the deformation simulated by the model is characterised by a multi-fractal spatial scaling.
Temporal scaling properties of sea ice deformation (or the degree of intermittency of sea ice deformation) can be studied from the dispersion of passive tracers (Rampal et al., 2008). Note that this approach is inspired by a classical methodology developed originally to study fluid turbulence (Richardson and Stommel, 1948;Batchelor, 1952). Here we perform the temporal scaling analysis for both the model and the observations following the same method as in Rampal et al. (2008), using pairs of vertices of the model mesh and pairs of tracking points of the RGPS trajectory data set, respectively. Indeed, in a Lagrangian modelling framework, each vertex of the mesh can be considered as a passive tracer of sea ice, and directly compared with tracking points of the RGPS data set. For each pair of vertices/RGPS points initially separated by a distance L of ∼ 30 km on average, a proxy of sea ice deformation D is measured by looking at the relative variation L/L of the distance between the two vertices/RGPS points for different time intervals t. The deformation rateḊ is estimated asḊ = L L t . For the model the intervals are 0. 25, 0.5, 1, 2, 4, 8, 16, 32, and 64 days, whereas the analysis for the RGPS trajectories only starts from a time interval of 2 days due to the limited temporal resolution of these data. The number of analysed measurements is similar to Rampal et al. (2008) and decreases from 630 000 for the 0.25 day interval to 2600 for the 64 days interval. Figure 12 shows the mean value of |Ḋ| for the different timescales for the model and the observations. A power-law model |Ḋ| ∼ T −α from T = 2 to 64 days fits both the observed and simulated data very well, with the same scaling exponent α = 0.3. Previous studies based on buoy data indicate that the scaling should also hold for smaller scales. This is not the case for the model data and cannot be verified from the RGPS data used in this study. The right panel of Fig. 12 indicates that the model only gives the right scaling exponent for healing timescales equal and larger than 7 days. We note that the order of magnitude of the healing timescale obtained here is consistent with the optimal value obtained in the previous section when analysing the sensitivity of the  Figure 11. Scaling analysis of sea ice deformation performed for the period January-May. The left panel shows the mean shear rate (i.e., the first-order moment of the distribution) computed for spatial scales ranging from ∼10 to 1000 km. A power law ˙ q ∼ L −β(q) is fitted to the data sets for each moment order q ranging from 0.5 to 3 (grey and black lines on the left panel for q = 1). The departure from the powerlaw fit at L < 20 km and L > 500 km comes from finite size effect (model/data resolutions and size of the Arctic basin, respectively). The scaling exponents β(q) for the other moment orders of the shear distribution are shown in the right panel. These values fit remarkably well with a quadratic function, which reveals the multi-fractal character of the scaling. The error bars correspond to the minimum and maximum exponents computed over two consecutive scales. . Temporal scaling analysis of sea ice deformation performed for the period January-May. Left panel shows the mean deformation rate (i.e., the 1 st order moment of the distribution) computed for temporal scales ranging from 6 h to 60 days. The deformation rate is here defined as in Rampal et al. (2008). A power-law fit ˙ q ∼ T −α(q) is calculated for the data sets for each moment order q ranging from 0.5 to 3 (grey and black lines on the left panel for q = 1) and the scaling exponents are found to be very similar between the model and the observations. This shows the model captures the observed intermittency of sea ice deformation. The right panel shows the sensitivity of the scaling exponent for the mean (q = 1) to the healing timescale. The green square corresponds to the exponent obtained for the reference run (shown on left panel). model with respect to sea ice drift. The low sensitivity to the healing time parameter for values larger than 14 days may indicate that the thermodynamical healing term is not needed and that the healing due to new ice formation is sufficient. However, as this may not be true for all model configurations, we prefer keeping the thermodynamical healing term in the description of the model.
The simulated mean value and spatial scaling exponent of the 3-day deformation evolves during the year of simulation (see Fig. 13), with much lower mean deformation between January and April, and a more negative scaling exponent in summer than in winter. We note that this behaviour, as well as the high variability, compares well with the results found by Stern and Lindsay (2009)

Evaluation of simulated sea ice extent and volume seasonal cycles and sensitivity to thermodynamical parameters
We now consider the modelled seasonal cycle in total ice volume and area. This section is intended to demonstrate that the model produces a reasonable seasonal cycle and to explore briefly its sensitivity to key parameters. An in-depth evaluation and tuning of these aspects of the model's behaviour would require several multi-decadal runs, which we consider outside the scope of this paper. For this purpose results from three runs, in addition to the reference run are shown: a run with fixed albedos of α i = 0.7 and α s = 0.9 (high albedo case), a run with temperature-dependent albedos , and a run forced with the ERA-Interim reanalysis results. Figure 14a shows the modelled total ice volume compared to monthly mean outputs from PIOMAS and observations from ICESat (Kwok et al., 2004, data downloaded from http://rkwok.jpl.nasa.gov/icesat/download.html on 6 March 2015). The estimates are calculated for the region . Both uncertainty estimates are probably upper bounds according to their authors. It is difficult to assess model performance in terms of total ice volume due to the lack of reliable observations. The uncertainty on the ICESat observations is substantial and the quality of the October-November estimate in particular is suspect. Because of this lack of data we chose to also plot the total volume from the PIOMAS model, but this should also only be considered a reference and not an accurate representation of the state of the ice cover. With these caveats in mind we see that the performance of the reference run is acceptable when it comes to ice volume. The melt rate can also be substantially affected by tuning the albedos, as expected. Figure 14b and c show the modelled sea ice extent and area respectively, compared to satellite observations. The area is simply the total area covered by sea ice, while the extent is the total area of grid cells covered with more than 15 % of sea ice. The observations shown are the mean values and extremes for daily observations using the ASI algorithm AMSR-E (Kaleschke et al., 2001;Spreen et al., 2008, data obtained from the Integrated Climate Date Center, University of Hamburg, Germany, http://icdc.zmaw.de), OSI-SAF (EUMETSAT Ocean and Sea Ice Satelitte Application Facility, 2015), NASA Team (Cavalieri et al., 1996), and bootstrap (Comiso, 2000) products. Using these four products gives a good idea of the uncertainty involved in the satellite observations .
In terms of extent, the results of the neXtSIM model are within the limits for the uncertainty estimates for the observations until the start of May. At this point the modelled melt is considerably more rapid than the observed one, leading to a difference of about 1.5 million km 2 at the beginning of June. From then, however, the observed melt becomes much more rapid and the end result is that the modelled extent at the extent minimum is within the limits of the observations uncertainties. Changing the forcing or albedos does affect the melt substantially, but it does not alter the fact that the model fails to capture the two-phased melt observed, a slow phase from early April to early June and a rapid phase from early June to early September.
In terms of total ice area, the model slightly overestimates the ice area during the freeze-up, but is in good agreement with observations for the rest of the model run. This is, however, not the case when using the ERA-Interim forcing or the temperature-dependent albedos since in those cases the melt is too rapid, resulting in total ice area that is about 1.5 million km 2 smaller than in the reference run.
There is, therefore, a discrepancy between the modelled extent and area when compared to the observations, in that the modelled extent is too low during melt but the modelled area is correct. This seems to indicate that as the ice concentration is reduced during melt, the ice compacts too easily, resulting in the correct area but too low extent. This issue is currently under investigation.
The spatial distribution of concentration is shown in Fig. 15. The concentration map shows that the sea ice distribution at minimum extent is well captured. The largest differences between the modelled and observed ice extent occur in the regions where the modelled ice concentration is low and the ice is easily influenced by the wind. Kauker et al. (2009) have shown the shape of the ice extent minimum to be heavily influenced by the wind, but we have not investigated this in our model.
In addition to concentration, Fig. 16 shows the spatial distribution of ice thickness at the beginning of the simulation, in midwinter and at the sea ice minimum. These distributions clearly show the high degree of heterogeneity that appears in the model, despite very smooth initial conditions. The midwinter map shows substantial amounts of fracturing and ridge formation in the Beaufort Gyre and the Transpolar www.the-cryosphere.net/10/1055/2016/ The Cryosphere, 10, 1055-1073, 2016 Drift Stream in particular. This heterogeneity persists until the end of the melt season, even if the melt does smooth it out somewhat. Overall, the model performs well in terms of total volume, area, and extent. This behaviour is largely controlled by the atmospheric and oceanic forcing. However, a poorly tuned or conceived ice model is still free to diverge considerably from the observed state, and it is reassuring to see that this is not the case here. The only genuine discrepancy between the model results and observations is that the model does not capture the two phases of melting observed in the extent. The model is sensitive to changes in the surface albedo, which is to be expected and albedos are probably the most widely used tuning parameters for ice and ice-ocean models. The model also shows some sensitivity to the lateral melt formulation, which is limited and was not shown. Sensitivity to the oceanic nudging timescale and various dynamical parameters is negligible within reasonable ranges for these parameters. For longer simulations a more sophisticated thermodynamics model is needed though, such as Bitz and Lipscomb (1999) or Winton (2000).

Summary and conclusions
In this paper we presented the first comprehensive version of the neXtSIM model, a fully Lagrangian dynamical/thermodynamical model for sea ice. The model is built around the dynamical core previously described in Bouillon and Rampal (2015b). It uses a novel approach to simulate the sea-ice-damaging process and the associated ice cover deformation and mechanical healing.
In order to be able to run simulations for seasonal timescales we have developed and implemented the following numerical and physical components into the model: -local dynamic remeshing accompanied with an efficient and conservative remapping scheme; -a thermodynamics model coupled to a slab ocean; -a healing parameterisation which simulates the restoration of mechanical strength due to refreezing of leads.
In order to evaluate the performance of neXtSIM we used a full Arctic set-up and ran the model for 13 months, starting on 15 September 2008, and using realistic atmospheric forcing. The main evaluation results are as follows: -the model reproduces the local motion of sea ice that is in free drift well; -the model also reproduces the drift of the pack ice well, at local (∼ 10 km) and large (Arctic-wide) spatial scales, and for daily to seasonal timescales; -the model captures the observed spatial multi-fractal scaling of sea ice deformation over 3 orders of magni-tude, from ∼ 10 to ∼ 1000 km, as well as its variability from winter to summer; -the model captures the observed intermittency of sea ice deformation over 2 orders of magnitude, from 1 to ∼ 100 days; -the model produces seasonal cycles of sea ice volume, area, and extent that are all in good agreement with observations.
In conclusion, for scales smaller than a year, neXtSIM performs very well with respect to several important metrics related to sea ice dynamics and thermodynamics. We believe that in its current stage of development, neXtSIM may already be a useful tool for both the scientific and engineering communities. For longer timescales and to study the interactions between sea ice and the ocean, ecosystems, or the atmosphere, more developments are needed, especially on the coupling with other components and the use of a more advanced sea ice thermodynamics model.

Data availability
The atmospheric reanalysis data used in this article are available at http://rda.ucar.edu/datasets/ds631.4/. The ocean reanalysis data used in this article are available at http://marine. copernicus.eu/. The sea ice concentration data sets used in this study are available at http://osisaf.met.no, http://icdc. zmaw.de, and https://nsidc.org/data/search/.