Impact of rheology on probabilistic forecasts of sea ice trajectories: application for search and rescue operations in the Arctic

. We present a sensitivity analysis and discuss the probabilistic forecast capabilities of the novel sea ice model neXtSIM used in hindcast mode. The study pertains to the response of the model to the uncertainty on winds using probabilistic forecasts of ice trajectories. neXtSIM is a continuous Lagrangian numerical model that uses an elasto-brittle rheology to simulate the ice response to external forces. The sensitivity analysis is based on a Monte Carlo sampling of 12 members. The response of the model to the uncertainties is evaluated in terms of simulated ice drift distances from their initial positions, and from the mean position of the ensemble, over the mid-term forecast horizon of 10 days. The simulated ice drift is decomposed into advective and diffusive parts that are characterised separately both spatially and temporally and compared to what is obtained with a free-drift model, that is, when the ice rheology does not play any role in the modelled physics of the ice. The seasonal variability of the model sensitivity is presented and shows the role of the ice compactness and rheology in the ice drift response at both local and regional scales in the Arctic. Indeed, the ice drift simulated by neXtSIM in summer is close to the one obtained with the free-drift model, while the more compact and solid ice pack shows a signiﬁcantly different mechanical and drift behaviour in winter. For the winter period analysed in this study, we also show that, in contrast to the free-drift model, neXtSIM reproduces the sea ice Lagrangian diffusion regimes as found from observed trajectories. The forecast capability of neXtSIM is also evaluated using a large set of real buoy’s trajectories and compared to the capability of the free-drift model. We found that neXtSIM performs signiﬁcantly better in simulating sea ice drift, both in terms of forecast error and as a tool to assist search and rescue operations, although the sources of uncertainties assumed for the present experiment are not sufﬁcient for complete coverage of the observed IABP positions.

Abstract. We present a sensitivity analysis and discuss the probabilistic forecast capabilities of the novel sea ice model neXtSIM used in hindcast mode. The study pertains to the response of the model to the uncertainty on winds using probabilistic forecasts of ice trajectories. neXtSIM is a continuous Lagrangian numerical model that uses an elasto-brittle rheology to simulate the ice response to external forces. The sensitivity analysis is based on a Monte Carlo sampling of 12 members. The response of the model to the uncertainties is evaluated in terms of simulated ice drift distances from their initial positions, and from the mean position of the ensemble, over the mid-term forecast horizon of 10 days. The simulated ice drift is decomposed into advective and diffusive parts that are characterised separately both spatially and temporally and compared to what is obtained with a free-drift model, that is, when the ice rheology does not play any role in the modelled physics of the ice. The seasonal variability of the model sensitivity is presented and shows the role of the ice compactness and rheology in the ice drift response at both local and regional scales in the Arctic. Indeed, the ice drift simulated by neXtSIM in summer is close to the one obtained with the free-drift model, while the more compact and solid ice pack shows a significantly different mechanical and drift behaviour in winter. For the winter period analysed in this study, we also show that, in contrast to the free-drift model, neXtSIM reproduces the sea ice Lagrangian diffusion regimes as found from observed trajectories. The forecast capability of neXtSIM is also evaluated using a large set of real buoy's trajectories and compared to the capability of the freedrift model. We found that neXtSIM performs significantly better in simulating sea ice drift, both in terms of forecast error and as a tool to assist search and rescue operations, al-though the sources of uncertainties assumed for the present experiment are not sufficient for complete coverage of the observed IABP positions.

Introduction
Large changes in the Arctic sea ice have been observed in recent decades in terms of the ice thickness, extent and drift (e.g. Kwok, 2007;Stroeve et al., 2007;Rampal et al., 2011;Stroeve et al., 2012). These changes, and their underlying driving mechanisms, still need to be fully understood in spite of their being fundamental for building confidence in the forecasting capabilities of current prediction systems. The need for a reliable sea ice prediction platform is particularly important in the modern context of growing economic opportunities with high societal and environmental impacts. For instance, the dramatic decline of sea ice cover in the Arctic is opening new shipping routes, fishing grounds and tourist destinations as well as access to a significant portion of the remaining hydrocarbon resources. Associated with this increasing activity are important risks for pollution of the Arctic environment and to human lives. High-quality predictions of ocean and sea ice in the polar regions are therefore needed in order to measure the risks, to plan future activities and to assist operations in real time.
Current short-term (i.e. within 10 days) sea ice forecasting systems integrate either a stand-alone sea ice model (RI-OPS, Lemieux et al., 2016;Dupont et al., 2015), a coupled ice-ocean model (e.g. ACNFS, Hebert et al., 2015, TOPAZ, Sakov et al., 2012Smith et al., 2015) or more seldom a coupled atmosphere-ice-ocean model (GloSea5, Williams et al., 2015). Seasonal to decadal climate forecasts are more common and include sea ice as part of the Earth system models (see e.g. Carrassi et al., 2016). The sea ice models used in these systems are usually derived from the work of Hibler III (1979), and they treat the sea ice as a continuous medium with a viscous-plastic rheology (Hunke and Dukowicz, 1997;Bouillon et al., 2009). In spite of this development, simple free-drift ice (i.e. in the absence of friction and internal forces) forecasts have remained in use by environment agencies (Grumbine, 1998(Grumbine, , 2003. The forecast skill of these systems based on a free-drift ice has been evaluated in deterministic mode, when a single "best" forecast is provided: despite the lack of realism in the free-drift assumption, the forecast skill of such systems is seen as difficult to beat (Schweiger and Zhang, 2015).
Probabilistic forecasts, widely used in weather forecasting (Molteni et al., 1996;Leutbecher and Palmer, 2008), are still in their infancy in sea ice forecasting. Probabilistic predictions rely on an ensemble of model simulations (i.e. a Monte Carlo simulation) used to describe the forecast uncertainty stemming from errors in the model parameters, initial and boundary conditions, and any external forcing. The resulting cloud of model outputs is used to retrieve statistical information, such as the ensemble mean and its spread (i.e. the standard deviation), which are thus used in place of the deterministic forecast and to estimate the associated uncertainty, respectively. The multiple simultaneous sources of errors usually make the forecast accuracy of the ensemble mean exceed that of the single deterministic prediction (Leith, 1974;Zhu, 2005), although often the spread underestimates the actual forecast error when the sources of error are not all adequately accounted for (Buizza et al., 2005). Monte Carlo techniques are already common practice in different areas (e.g. Dobney et al., 2000;Hackett et al., 2006;Breivik and Allen, 2008;Melsom et al., 2012;Motra et al., 2016;Duraisamy and Iaccarino, 2017) and a common tool for sensitivity analysis.
This study concerns the probabilistic forecast capability of the sea ice model neXtSIM (Rampal et al., 2016b). The work is carried out by performing a Monte Carlo sensitivity analysis of the model with respect to uncertainties in the surface wind velocity. The first goal is to highlight the role of the ice rheology in the ice drift: how do the ensemble mean drift and its standard deviation respond to uncertainties in the wind forcing? To answer this question, we compare the ice drift obtained from neXtSIM to one obtained from a free-drift model. In the second part, we study the skill of the probabilistic forecast using Lagrangian trajectories departing from independent virtual drifting buoys and compare them with real observations without aiming to make it a key objective. We use the conceptual framework of search and rescue operations where a probabilistic forecast is commonly used to draw the search area of the ocean where drifting objects are likely to be found (Hackett et al., 2006;Breivik and Allen, 2008;Melsom et al., 2012). Contrary to these studies, the present simulations are in the context of a "hindcast-forecast", using reanalysed atmospheric forcing fields but assuming that they are affected by errors with the statistical properties that could be expected from a numerical weather forecast in the Arctic. For simplicity, we will use the word "forecast" instead of "hindcast-forecast" throughout the paper.
Our main research tool and object of study is the sea ice model neXtSIM. The model neXtSIM is based on a Lagrangian numerical scheme and on a continuous approach using a newly developed elasto-brittle (EB) ice rheology. This mechanical framework is inspired by the scaling properties of sea ice dynamics revealed by multi-scale statistical analyses of observed sea ice drift and deformation (Marsan et al., 2004;Rampal et al., 2008;Bouillon and Rampal, 2015b) as well as by the in situ measurements of sea ice internal stresses showing that sea ice deformation is accommodated by Coulombic faulting (Weiss et al., 2007;Weiss and Schulson, 2009). For 40 years, a large variety of sea ice models have been developed. Some, like neXtSIM, treat the sea ice as a continuous medium, but with different rheologies (e.g. Coon et al. (1974) and Hibler III (1979) modelled sea ice as an elasto-plastic material; Hunke and Dukowicz (1997) as an elasto-visco-plastic material; and Dansereau (2016) as an Maxwell elasto-brittle material), and are suitable for high ice concentration (> 80 %), while others that treat the ice as a discrete medium (Hopkins et al., 2004;Wilchinsky et al., 2010;Herman, 2011;Rabatel et al., 2015) are more suitable for low ice concentration (< 80 %) such as within the marginal ice zone.
We concentrate here on the impact of the error from the wind field alone. The reasons are twofold. First, the wind is the most influential external force affecting sea ice motion. About 70 % of the variance of the sea ice motion in the central Arctic can be explained by the geostrophic winds (Thorndike and Colony, 1982). However, the sea ice response to winds strongly depends on its degree of damage; sea ice responds in a linear way only when it is fragmented into small floes. Indeed, in this case, the internal forces are negligible and the inertial term is linearly related to the air and water drags, whereas this behaviour drastically changes when considering a large, continuous and undamaged solid plate. The second reason is that surface wind velocity fields provided by atmospheric reanalyses contain large uncertainties in the Arctic due to the limited number of observations. Previous sensitivity analyses of the neXtSIM model have been performed with respect to initial conditions and to some key sea ice mechanical parameters (see Sect. 4 in Bouillon and Rampal, 2015a). These analyses consisted in running the model with different values of the input sources. This allowed the authors to explore and quantify the sensitivity of the ice velocity with respect to the ratio between water and air drag coefficients and of the ice deformation with respect to the compactness parameter value (see Eq. 5), the sea ice cohesion value (see Eq. 10 in Bouillon and Rampal, 2015a), the initial concentration field or the initial thickness field. Al- though these analyses did not use the latest model developments of neXtSIM (in particular they included neither the thermodynamics nor the re-meshing process), the impact of the mechanical parameters on the ice deformation can still be considered as valid. Systematic errors in the mean sea ice drift are evaluated by averaging modelled and observed drift from the OSI-SAF dataset (Lavergne and Eastwood, 2015) over the period between 1 January 2008 and 30 April 2008 and over boxes of 100 × 100 km 2 covering the whole Arctic (see Fig. 1). The largest differences between the observed and simulated mean ice drift are located in the Beaufort, Chukchi, Kara and Barents seas and Fram Strait and in some areas of the East Siberian Sea. In the rest of the domain the error on the mean winter drift is only less than 3 km day −1 , consistently with Rampal et al. (2016b). This paper is organised as follows: Sect. 2 gives a general presentation of the sea ice model neXtSIM with the main equations describing the sea ice dynamical behaviour; Sect. 3 presents the details of the sensitivity analysis based on a Monte Carlo sampling, including the description of the quantities of interest, the construction of the wind perturbations, and the general experimental setup. In the same Sect. 3, we also define the free-drift model that will be used for comparison and benchmark against neXtSIM. Section 4 discusses the results for the ensemble mean, spread and the evaluation of the forecast skills comparing neXtSIM to the free-drift model. Final conclusions are drawn in Sect. 5.

General information on the model neXtSIM
In this section, we provide a general description of neXtSIM. Deliberately, we choose to not go through all model equations here but rather list those that are needed to get an overall understanding of how the model works and that are relevant for the present study. For a more detailed description of the model see Bouillon and Rampal (2015a) and Rampal et al. (2016b).
neXtSIM is a continuous dynamic-thermodynamic sea ice model. It uses a pure Lagrangian advection scheme, meaning that the nodes of the model mesh are moving at each time step according to the simulated ice motion. The model mesh is therefore changing over time, is not spatially homogeneous and can locally become highly distorted, that is, when and where the ice motion field is showing strong spatial gradients. In this case, a local and conservative re-meshing procedure is applied in order to keep the numerical integrity of the model and the spatial resolution of the grid approximatively constant during the simulation. The equations are discretised on a triangular mesh and solved using the classical finite element method, with scalar and tensorial variables defined at the centre of the mesh elements and vectors defined at the vertices. The model uses a mechanical framework that has been developed recently (Girard et al., 2009;Bouillon and Rampal, 2015a) and is based on the EB rheology. The brittle mechanical behaviour of the sea ice is simulated by calculating the local level of damage in each grid cell, a variable which is not considered in classical viscous-plastic sea ice models typically used in the sea ice modelling community. Sea ice thermodynamics, which is parametrised in neXtSIM as in the zero-layer model of Semtner (1976), controls the amount of ice formed or melted at each time step. When a volume of new (and therefore undamaged) ice is formed within a grid cell by thermodynamical freezing, the mechanical strength of the total volume of ice covering that cell is partially restored, and the new damage value is computed as a volume-weighted mean. Note, however, that the damaging process is very fast (i.e. about few minutes) while the mechanical healing process is occurring over much slower timescales of about several weeks. The sea ice variables used in neXtSIM are the following: h and h s are the effective sea ice and snow thickness, respectively (ice and snow volumes per unit area); A is the sea ice concentration (bounded to 1); d is the sea ice damage ranging from 0 (undamaged ice) to 1 (fully damaged); u is the horizontal sea ice velocity vector; and σ is the ice internal stress tensor. The model has two ice thickness categories: ice and open water.
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 the divergence of the horizontal velocity and S φ a thermodynamical sink- 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 parametrised as in Lemieux et al. (2015). The last terms are the Coriolis parameter, f , the upward pointing unit vector, k, the gravity acceleration, g, and the ocean surface elevation, η. The internal stress σ is computed as in Bouillon and Rampal (2015a) and Rampal et al. (2016b). Its evolution equation can be written as where d is the damage and˙ is the deformation rate tensor defined as˙ = 1 2 ∇u + (∇u) T . C can be written as where ν is the Poisson ratio and E(A, d) is the effective elastic stiffness of the ice, which depends on the ice concentration A and the damage d according to where Y is the sea ice elastic modulus (Young's modulus) and α is the so-called compactness parameter. The evolution equation for the damage is written as where d is a damage source term calculated as in Rampal et al. (2016b) (Eq. 8), and S d is thermodynamical sink term which depends on the volume of new and undamaged ice formed over one time step as well as on time (See Rampal et al., 2016b, Sect. 2.3, for more details). The air and oceanic drags, respectively τ a and τ w in Eq.
(2), are written as a force per unit area in the quadratic form using the associated turning angle (Leppäranta, 2011) where . , R θ a , R θ w , u a , u w , ρ a , ρ w , C a and C w are, respectively, the Euclidean norm in R 2 , the rotation matrix through the angle θ a and θ w , the wind velocity, the ocean current, the air density, the water density, the air drag coefficient and the water drag coefficient. The values of the model parameters that are used for the simulations presented in this paper are listed in Table 1. In this study, we perform a sensitivity analysis using a statistical approach based on Monte Carlo sampling of the model inputs. We focus on the response of the model to the uncertainties in the wind velocity field. In particular, we are looking at the response of sea ice drift to wind perturbations representing these uncertainties. Our methodology is based on simulating Lagrangian trajectories of virtual buoys using an ensemble run of the neXtSIM model forced by slightly different (i.e. perturbed) wind forcing (see Sect. 3.2 for more details on the generation of the perturbed winds). The velocity of a given virtual buoy is calculated online, at each time step, as a linear interpolation of the velocities simulated at the nodes of the mesh element containing that buoy (see Lagrangian approach in Sect. 2). Each virtual buoy is associated with an initial position x 0 ∈ D, with D being the initial domain, and a start date t 0 ∈ Y , where Y is the time period of interest of this study (see Sect. 3.2 for more details). A buoy trajectory is denoted g(x 0 , t 0 , t) with t ∈ [t 0 , T ], where T defines the duration of the individual simulations. For each initial position x 0 and start date t 0 , we simulate N trajectories g i i∈{1,...,N } from N model runs, each one corresponding to a different realisation of the wind forcing. If a buoy ends up in an ice-free element, it is then untracked further and its trajectory discarded from the remaining analysis.
For each ensemble member (trajectory), we define the following Euclidean distances The Cryosphere, 12, 935-953, 2018 www.the-cryosphere.net/12/935/2018/ where the quantity r i (t) is the distance of the member position at time t, g i (x 0 , t 0 , t), from its departure origin, x 0 = g i (t = t 0 ). The second quantity, b i (t), represents the distance between the member position at time t and the ensemble mean position (i.e. the barycentre, B(t), of the ensemble), , at the same time t (see the top panel of Fig. 2). We make use here of the convention of using boldface for vectors and matrices and normal face for scalar quantities; hereafter, we drop the explicit mention of the dependence on x 0 and t 0 , to simplify the notation.
Furthermore, we define a two-dimensional timedependent orthonormal basis, centred on B(t), whose axes are two perpendicular lines, one of which connects x 0 to B(t). The coordinates of g i (t) on this basis are hereafter denoted as b i, (t) and b i,⊥ (t), as illustrated in the bottom panel of Fig. 2; they provide information on the spatial and temporal evolution of the ensemble spread and shape and can also be used to look at how the virtual buoy positions are distributed around the ensemble mean over time.
With the individual r i and b i in hand, we compute basic, second-order statistics. Let us consider their means, µ r and µ b , and the standard deviations, σ b and σ b ⊥ , of the components b and b ⊥ , as our main quantities of interest in the analysis that follows. We note that the mean of b i, and b i,⊥ is zero (being barycentric coordinates) and do not appear in the calculation of standard deviations. Throughout the rest of this paper, σ b (t) and σ b ⊥ (t) are only used to compute the ratio which provides a measure of the anisotropy of the ensemble spread of the virtual buoys positions around the barycentre B of the ensemble. It is finally worth observing that the two quantities, r and b, provide complementary information: the former on the advective component of the motion and the latter on its diffusive part. The ensemble mean distance from the starting point, µ r , is a statistical estimate of the distance travelled by an ice parcel according to the ice advection properties of the motion field, while µ b is the (mean) spread relative to the aforementioned distance and accounts for the diffusion properties of the motion; see the top panel of Fig. 2.

Experimental setup
Our domain of study is the region covering the Arctic Ocean. While the coasts are considered as closed boundaries, open boundaries are set at the Fram and Bering straits (see Fig. 3).
The wind forcing is taken from the Arctic System Reanalysis (ASR) (Bromwich et al., 2012). This reanalysis product provides wind speeds and directions at 10 m, every 3 h, at a horizontal resolution of 30 km. No turning angle has been applied (see Table 1). For every 3-hourly wind field, we generate spatio-temporal correlated perturbations as described in Evensen (2003) and then add them to the ASR wind field. This procedure is identical to the one used to produce ensemble runs with the coupled ocean-sea ice model TOPAZ (Melsom et al., 2012) and constitutes the propagation step in the ensemble Kalman filter (Sakov et al., 2012). The method is designed such that the perturbed wind fields keep important physical properties; that is, the wind perturbations are geostrophic (gradients of random perturbations of the sea level pressure) and the wind divergence is kept unchanged. They are built on random stationary Gaussian fields, with a Gaussian spatial covariance function, dimensionalised by the wind error variance and correlated in time. Time series of wind perturbations are assumed to be red noise. For our www.the-cryosphere.net/12/935/2018/ The Cryosphere, 12, 935-953, 2018 Although the ensemble average of the perturbed u and v components of the winds is equal by construction to the original winds provided by the ASR, the wind speed is positively biased. The value of the air drag coefficient (C a in Eq. 7) had previously been optimised in the neXtSIM model when forced by the ASR following the method presented in Rampal et al. (2016b, Sect. 3.2) and set to 7.6 × 10 −3 . We applied the same method here to tune the value of C a so that the simulated ice drift compares best with the observed ice drift from the OSI-SAF dataset (Lavergne and Eastwood, 2015). The optimisation is carried out at all times between 1 January and 30 April 2008 but limited to the region where the ice is in free drift (see Fig. 4), i.e. where the ensemble-averaged simulated ice velocity differs by less than 10 % from the drift simulated by the free-drift (FD) model. Figure 5 shows the comparison, after optimisation of the air drag coefficient, between the observed and simulated ice velocities. As expected for a wind dataset positively biased in magnitude compared to the original one, we found an optimised value for the drag coefficient C a = 5.1 × 10 −3 , lower than the one used in Rampal et al. (2016b) The ocean forcing comes from the TOPAZ4 reanalysis (Sakov et al., 2012). TOPAZ4 is a coupled ocean-sea ice system combined with an ensemble Kalman filter data assimilation scheme assimilating both ocean and sea ice ob- servations. In our simulations, we used the 30 m depth currents, to which we apply a turning angle of 25 • , the surface temperature and salinity, and the sea surface height, all provided as daily means with an average horizontal resolution of 12.5 km, following Rampal et al. (2016b).
Our analysis is based on two periods of the year 2008: from 1 January to 10 May and from 1 July to 20 September, representative of the winter and summer conditions, respectively. We have intentionally studied them separately because winter and summer are characterised by significantly different sea ice mechanical regimes and therefore drift responses. During the winter, the whole Arctic basin is covered by ice and its concentration is close or equal to 100 %; that is, the internal stresses in the ice, and the corresponding ∇ · (σ h) term in Eq. (2), become large and of the same order of magnitude as the wind drag term. As a consequence, the ice drift is (on average) much reduced. During the summer period, in contrast, the ice concentration is lower and the ice pack does not generally reach the coasts, the ice internal stresses are much closer or equal to zero, and the ice drift closer to a free-drift state (see text below). We note, however, that the wind field perturbations are generated using the same aforementioned procedure, for both the winter and the summer, and have thus the same spatial and temporal properties.
We ran a total of 13 simulations in the winter and 8 in the summer during successive, non-overlapping 10-day periods. Limiting the length of the simulations to 10 days ensures that the sea ice state (thickness and concentration) remains as realistic as possible in the free-drift simulation, in which there are no physical limits to the amount of ridging and opening. The starting positions are separated by 100 km and cover the domain displayed in Fig. 3. All ensemble members start from the same initial conditions extracted from a previous -deterministic -neXtSIM simulation by Rampal et al. (2016b)  ice variables: h, h s , A, d, u and σ . We ran an ensemble of 12 members, each of them forced by the perturbed wind generated as explained above. We performed (not shown) a convergence analysis of our results as a function of the ensemble size from N = 3 to N = 20 and observed a convergence from about N = 10 with only minor changes for N ≥ 12, and we are thus confident that N = 12 suffices to our purposes. From these ensemble runs we simulated a total of over 96 000 ( 8000×12) virtual buoy trajectories over the winter season and over 38 000 ( 3200 × 12) trajectories over the summer season. This dataset was used to run the analyses described in Sect. 4 and presented at the 19th EGU General Assembly .
As already stated, we compared neXtSIM with the socalled FD model, so that all simulations that follow have been carried out for the two models. neXtSIM (Eq. (2 with all terms on its right-hand side included) is our reference model. The FD model is equivalent to neXtSIM except that it considers the following simplified version of the momentum equation in which the terms related to the sea ice rheology, the basal stress and the inertial term are neglected: In Eq. (11) the water and air drag forces, the Coriolis force and the gravity force due to the ocean surface tilt are balancing each other. The FD model is therefore analogous to the steady-state drift of an object at the surface. We run the FD model with the same initial conditions as neXtSIM except that d and σ are not used. The drag coefficient is also optimised for the FD model at a value of 3.2 × 10 −3 , which is lower than for neXtSIM, as expected. The optimisation method used for FD is the same as for neXtSIM described above, except that the OSI-SAF drift vectors are used everywhere.

Results
In this section, the notations < .> W and < .> S correspond to winter and summer averages (i.e. over all the 13 and 8 simulation periods of 10 days), respectively. The notations < .> D correspond to the spatial mean over the domain. When considering both spatially and temporally averaged quantities, we use the notations < .> W,D or < .> S,D .  (Fig. 7) and summer (Fig. 9). Note that both figures are relative to neXtSIM, but the freedrift wind speed is identical (same perturbations) and the ice thickness geographical pattern very similar; we have thus not displayed them to avoid redundancy. From Figs. 6 and 8 we see that neXtSIM gives a smoother response to perturbed forcing than the FD model in terms of mean advective drift µ r and mean diffusive spread µ b , in both winter and summer. Indeed, we observe in neXtSIM a clear spatial coherency in both the advection and diffusion of the ice buoys over the domain that is almost absent in FD. We believe that this behaviour is related to the mean ice thickness pattern and, to a lesser extent, to the mean wind speed pattern (see Figs. 7 and 9 for winter and summer, respectively).

Spatial patterns
For neXtSIM, the smallest values for the mean of µ r and µ b averaged over the winter time period are found in the area located north of Greenland and the Canadian Archipelago, which is where the ice is the oldest, thickest (> 4 m) and mechanically the strongest and where the winds are on average weaker as compared to the rest of the Arctic. However, in the surrounding seas (i.e. Beaufort, Bering, Chukchi, Kara and Barents seas from west to east), where the ice is thinner and the winds stronger, the means of µ r and µ b are larger. Note that in summer these correlations or anticorrelations are even stronger, for example between the means of µ b and the ice thickness (see Figs. 8 and 9). For FD, the mean values of µ r are correlated to the wind speed in winter and, to a lesser extent, in summer (left panels in Figs. 7 and 9). It is worth noting that, despite the presence of thick ice in the north of the Canadian Archipelago and low winds, the ice is still advected significantly, as opposed to what is obtained with neXtSIM. Moreover, the spatial pattern of the mean of µ b shows no spatial coherence and resembles the random patterns from the wind perturbations. It is clearly visible in summer, while in winter the sea ice thickness field stays discernible. This may be due to the presence of the ice mass in Eq. (11).
In both winter and summer, the time-averaged response of µ r and µ b to wind perturbations is overall lower in neXtSIM than in FD (except in summer when µ r is 7 % larger in neXtSIM). This can be attributed to the ice rheology being turned on in neXtSIM, thus acting as an additional filter on the momentum transferred from the wind to the ice. In more detail, it is interesting to note that the magnitude of the impact of the ice rheology is different depending on whether we consider the drift distance by advection r or the spread distance by diffusion b and consider the winter or the summer. On average over the winter, µ r (t)  Wind speed (m s −1 ) Ice thickness (m) Figure 9. Summer average of wind speed and ice thickness. Both maps are from the neXtSIM simulations, but a similar thickness field and the exact same wind speed field are obtained for the FD simulations.
at t = 10 days, whereas over the summer µ b (t) D is 21 % lower. This large difference between the two distances, especially in winter, is probably related to the high ice concentration making sea ice harder to break up, and it keeps the members closer to each other. During summer, the ice is generally much less packed and the physical and dynamical differences between neXtSIM and FD have a lower impact. The 7 % larger values of µ r for neXtSIM are likely related to the optimisation of the air drag coefficient that returned a larger coefficient for neXtSIM. As expected, for neXtSIM we observe an increase of µ r (t) of about 51 % and µ b (t) of about 69 % in summer compared to winter. This behaviour differs drastically from the FD for which the values are nearly the same for both periods, and it is presumably related to the decrease in ice concentration due to the summer melting. The averaged sea ice concentra-tion over the whole domain in winter is about 0.99 while it drops to 0.83 in the summer. In neXtSIM, this strongly influences the mechanical behaviour of the sea ice since the effective elastic stiffness E depends non-linearly on the ice concentration (see Eq. 5). Assuming no change in the average level of damage of the ice, a drop by 15 % of the ice concentration between winter and summer implies a reduction of E by 96 %. This reduction of E leads in turn to a significant decrease of the internal stresses within the ice, thus lowering the term ∇ · (σ h) in Eq. (2), which makes the buoys' drift in neXtSIM closer to the ones obtained with the FD model.
The absolute values of µ r and µ b obtained by our analysis reveal that the advection part of the motion is in general larger than the diffusive part, independently of the season under consideration. In FD the ratio γ = µ r (t)/µ b (t) at t = 10 days is about 4.5. In neXtSIM, however, the ice rhe- ology acts to increase this ratio to 7. However, this value presents a strong spatial variability depending on the local thickness and wind speed. Where both are large, γ is large. For example, such areas are observed in the Fram Strait in winter (γ > 10) and in the central Arctic in summer (γ > 12). Where both ice thickness and wind speed are small, γ is small. For example, this is the case around the new Siberian islands in winter (γ < 4) and close to the ice pack edge in summer (γ < 6). Figure 10 shows the probability density function (PDF) of b (t) and b ⊥ (t) at t = 10 days for both neXtSIM and FD and for winter and summer (see Sect. 3.1). The PDFs of b and b ⊥ for the FD case are almost identical, and so we chose to only display one curve (black dashed line). The first aspect to remark from Fig. 10 is that all distributions are unimodal and symmetric, suggesting that the 2-D shape of the ensemble is symmetric around its barycentre B. However, we notice that the ensemble is anisotropic in neXtSIM, i.e. the distributions of b and b ⊥ differ substantially, whereas it is close to isotropic in FD. Figure 11 shows the temporal evolution over 10 days of the Arctic averaged ratio R (Eq. (10), which defines the degree of anisotropy of the ensemble spread (1: isotropic; > 1: anisotropic). We observe on the one hand that R is very close to 1 and relatively constant over time in the FD model. On the other hand, it is systematically larger in neXtSIM, especially in winter, and it also displays a certain short-term variability. Here, again, we encounter the peculiar effect of the  Figure 11. Evolution of the spatial mean of R(t) from t = 1 to t = 10 days for the winter (blue) and summer (red) periods for neXtSIM (solid lines) and FD (dashed lines).

Spatial and temporal properties of the ensemble spread
neXtSIM mechanical response to the external forces, which is to break up and deform along fractures that are dispersing the different members of the ensemble along a preferential direction; such behaviour cannot be reproduced by the FD.
Note also that R is as large as 2 within the first 2 days for neXtSIM in the winter, and then it decreases monotonically for t > 2, still remaining very large (between 1.4 and 1.6 at t = 10 days). This reveals that the ice will first tend to move compactly along the initial fractures (identical for all members at t = 0) away from the origin, but it then starts to break and, after 2 days, the damage pattern becomes significantly different within each member, leading to a more isotropic ice dispersion away from the barycentre. In Fig. 12, we show the maps of the R(t) values computed for each ensemble of trajectories at t = 10 days. These values are represented as coloured squares centred on the starting point x 0 . We observe that highest degree of ensemble anisotropy (R > 1) is found north of Greenland and Canadian Archipelago, where the ice is the thickest and the ice drift and winds the lowest, in overall agreement with the interpretation of the temporal evolution of R for neXtSIM in the winter, provided in relation with Fig. 11. Globally, we observe a high (> 1.5) anisotropy close to the coasts that can be explained by the ice pressure that counteracts sea ice motions towards the coasts (and the associated dispersion as well). In summer, large stretches of the coasts are ice-free and the increase of R is less visible. This is in contrast to the pattern from FD. In absence of internal stresses, the pattern of the anisotropy exhibits no spatial coherence and is similar in both winter and summer periods. Furthermore, as already noticed from Fig. 11, the values obtained for neXtSIM are systematically larger (by about 65 %) than for FD during the winter whereas only 8 % larger during the summer. However, and remarkably, the values of R, and thus the anisotropy of the ice drift, for neXtSIM exhibit marked spatial correlations that are absent from FD.
Another important characterisation of the ensemble spread evolution can be set by looking at the variance of the dis-  tance b between the virtual buoys and the barycentre B over time. The goal is to identify the diffusion characteristics of the ensemble, which can be interpreted in the framework of the turbulent diffusion theory of Taylor (1921). Similar Lagrangian diffusion analysis has been applied to study the regimes of diffusion of surface drifters in the ocean (e.g. Zhang et al., 2001;Poulain and Niiler, 1989) and, more recently, of buoys fixed to the ice cover (e.g. Rampal et al., 2009;Lukovich et al., 2015;Gabrielski et al., 2015;Rampal et al., 2016a). In the analysis performed here, the distance b to the barycentre of the ensemble corresponds to the fluctuating part m of the motion m in the so-called Taylor decomposition m = m + m . Figure 13 shows the temporal evolution of the ensemble average of the distances b i averaged over the Arctic domain D calculated form the buoy's tracks simulated with neXtSIM and FD. We found that the ensemble spread follows two distinct diffusion regimes, one for small time, t , and one for large time, t , where is the so-called integral timescale (Taylor, 1921). In neXtSIM, the first regime we found for winter corresponds to the ballistic regime where b i 2 D ∼ t 2 , and the second to the Brownian regime where b i 2 D ∼ t. These results are in agreement with the wintertime sea ice diffusion regimes revealed by applying the Lagrangian diffusion analysis to the buoy trajectories dataset of the International Arctic Buoy Programme (Rigor, 2002;Rampal et al., 2009) in wintertime and show that our experimental setup based on ensemble simulations forced by perturbed winds does not alter the capability of the neXtSIM model to reproduce these properties, as also shown recently in Rampal et al. (2016a) for the same winter. How-ever, we note that the regime we obtain with neXtSIM for the summer 2008 is super-diffusive, with b i 2 D ∼ t 1.15 for t , and therefore in apparent contradiction with Rampal et al. (2009), who found that sea ice follows a same Brownian regime in both winter and summer when averaging over the period . We suggest that this may actually be the fingerprint of a change in the dynamical behaviour of sea ice in summer that occurred over the most recent years (including 2008), in which the rheology plays a weaker role than it did in the 1980s and 1990s. This is also supported by the results we obtain here with the FD model that neglects the rheology and exhibits super-diffusive regime for 2008, regardless of the season considered.

Predictive skills of neXtSIM and of the FD models
We evaluate here how well the neXtSIM and FD models are able to forecast real trajectories in hindcast mode. As a benchmark, we compare the ensemble runs from both models to 604 (in winter) and 344 (in summer) observed trajectories from the IABP dataset. The simulated trajectories of both neXtSIM and FD are initiated on the same initial positions and at the same time as the IABP buoys (displayed in Fig. 14); the positions of IABP buoys are known every 12 h. It is important to note that most of these buoys were deployed in regions of thick and compact ice, the drift of which is largely influenced by the sea ice rheology. Therefore, we expect the FD model to be less competitive than if the comparison data had been uniformly distributed across the Arctic. As a metric for the models skill intercomparison, we use the linear forecast error vector defined as the distance between the observed IABP buoy position, O(t), and that of the ensemble mean, B(t) (see also Fig. 16). The components of e(t) onto the orthonormal basis centred on O (see Sect. 3.1, Figs. 2 and 16) read e (t) and e ⊥ (t).
We complete this evaluation comparing results from both models with those from a single deterministic forecast in order to verify the advantage of probabilistic forecasts. In this case, we run neXtSIM with parameters found in Rampal et al. (2016b), except the air drag coefficient has been re-tuned to C a = 6.5×10 −3 and unperturbed winds. For this new air drag coefficient, the same optimisation process as the probabilistic case is used against the same observations. Figure 15 shows the average norm of the forecast error, e , and of its components, e (t) and e ⊥ (t), as a function of time, for the experiments with neXtSIM and FD and for both winter and summer. Results reveal that the forecast error is smaller in neXtSIM than FD in both seasons. In winter, the error of the FD model grows almost twice as fast as the error of neXtSIM, up to 26 km at day 10 compared to about 15 km for neXtSIM. As already deduced from the results in the previous section, the mechanics underlying of the ice drift in neXtSIM and FD are similar in the summer, and this is reflected by the two errors being much closer to each other: the difference between the two increases slower, reaching 3 km after 10 days (see the left panel in Fig. 15).
The centre panel in Fig. 15 shows a positive bias of the error in the along-drift component (e ) for both models and both periods except for neXtSIM in winter, which presents a negative bias. The general positive biases betray a too-fast drift in the direction along the ensemble mean drift compared to the observations. Nevertheless, the bias for winter in neXtSIM is 2.5 times smaller than in the FD model, whereas both models perform similarly in the summer.
Finally, the right panel in Fig. 15 also reveals a bias of the error in the direction across the ensemble mean drift, that is substantially weaker than in the previous case. For FD, e ⊥ still being negative for both periods corresponds to a drift too far to the right compared to the observations. This bias to the right could be further reduced by a separate tuning of the turning angle θ w for the FD and neXtSIM models.
Overall, we conclude that the performances are significantly better for neXtSIM in winter but similar in summer, and this would likely remain so even after optimal tuning of the turning angle.
Comparing to a single deterministic neXtSIM forecast, we note that the forecast error is close to the average of the probabilistic run but larger in summer, reaching 34 km at 10 days. The main difference with the probabilistic run is the poorer along-drift component e . Indeed, the error is closer to zero in winter and increases to 15 km in summer.
In Hackett et al. (2006) and Breivik and Allen (2008), Monte Carlo techniques are used to forecast the drift of an object on the ocean surface. They associate the density of trajectories at their end points to a density of probability and use them to define a search area, within which the object is likely to be found. The search area is characterised by a surface centred on the ensemble mean and which size increases with the ensemble spread. The same methodology is followed here for forecasting the location of an object on drifting sea ice. In the context of rescue operations, the search area should be large enough to contain the actual position of the object but not excessively large, so as to keep the rescue operation time period and resources affordable and efficient. The forecast system should therefore ideally yield a high probability of finding the object in the search area, while keeping at the same time the search area as small as possible for the cost efficiency of the rescue procedure.
The probability of finding a drifting object inside the search area is referred to as the probability of containment (POC) and computed by counting the objects falling within the search area divided by the total number of objects. The POC may be interpreted as the ratio of the size of the search area to the square forecast error. Thus, a small forecast error compared to the search area leads to a strong POC; conversely, a small search area (ensemble spread) compared to the forecast error leads to a poor POC. In order to evaluate the probabilistic forecast capabilities of neXtSIM and the more classical FD model, the context of a search and rescue operation is adopted. We assume that an IABP buoy has been lost for 10 days: its initial position, x 0 (see Fig. 14), is assumed to be its last known position. The search area is then defined as the smallest ellipse centred on the ensemble mean position, B(t), encompassing all simulated members of the ensemble at time t. The main axes of the ellipse, a and a ⊥ , are aligned, respectively, with the parallel and perpendicular directions from the initial position, as defined in Sect. 3.1 and illustrated in Fig. 16. Similarly to Eq. (10) an anisotropy ratio R = a /a ⊥ can be defined: R can be large due to the sea ice rheology. A search area defined in this way increases with the ensemble spread and contains 100 % of the ensemble members.
Due to a lack of related literature for search and rescue in sea ice, we consider the values of open-ocean search areas and POC found in Breivik and Allen (2008) and Melsom et al. (2012) as references. These are, respectively, of the order of 1000 km 2 and 0.5 after 2 days of drift in the North Atlantic. We do not expect, however, a direct correspondence of these values to those of this section. First, the sea ice is solid and held together by the ice rheology, in particular in high concentration areas, so that the ensemble spread is expected to be smaller than in the open ocean. Second, the currents in the North Atlantic are generally stronger than in the Arctic Ocean. Finally, the search areas may be more complex than just an ellipse; it may well be a set of disjoint areas, each one with an associated different POC (e.g. Abi-Zeid and Frost, 2005;Breivik and Allen, 2008;Guitouni and Masri, 2014;Di Maio et al., 2016). Figure 17 shows the evolution of the ellipse areas, averaged over all IABP buoys. The increase is nearly linear for both model configurations and seasons. After 2 days of drift in neXtSIM, the area does not exceed 100 km 2 in summer and not even half as much in winter. The area is larger in  FD, and there is very little difference from winter to summer. The area for the FD is around 200 km 2 after 2 days and it reaches 500 km 2 after 5 days. The search area in FD is about 7 times larger than in neXtSIM in the winter and 2.5 times larger in the summer. Therefore, even if the forecast errors are smaller in neXtSIM than in FD, its shrunken search areas lead to a smaller POC for neXtSIM than for the FD model (not shown).
On Fig. 18, we show the spread-error relationship for both periods and both models. The curves represent the spatial mean of the forecast error. Overall, the curves are above the black line, indicating that the forecast error is larger than the spread for both models. The probabilistic forecasts from both model during both period are therefore too optimistic: they underestimate the uncertainties of their forecast. However, it is interesting to note two properties of neXtSIM. First, for spreads larger than 4 km, the forecast error from neXtSIM becomes independent from the spread, unlike FD, the errors of which grow monotonically. Second, for large spreads (greater than 3.5 km in winter and 6 km in summer) the curves from neXtSIM are consistently below those from FD and getting closer to the spread. Contrarily to the previous results, the FD and neXtSIM models behave very differently in the summer.
The small values of the spread correspond to shorter forecast lead times (see Fig. 17) and these are the times when the neXtSIM model is still heavily influenced by its initial conditions of damage, as previously noted on the anisotropy ratio (Fig. 11). As the damage is irrelevant to the FD model, the initial error grows slower initially, but keeps growing while the rheology maintains the errors closer to the spread in neXtSIM.
It should be no surprise that the two models underestimate the errors since this is a common behaviour of probabilistic forecast systems, but the differences of the shape of the spread-error relationships indicate that the two models un- derestimate the errors for different reasons: the neXtSIM ensemble is lacking spread during the initial times of the forecast but the asymptotic convergence of the spread to the errors tends to be the cause of the constitution of the initial ensemble.
If one had considered the more linear spread-error relationship in the FD model alone, it would have been tempting to increase the variance of wind perturbation errors until a perfect match of the spread to the errors was obtained, but this would have over-tuned the variance of the wind and masked that the FD model suffers from unresolved physics.

Relevance for search and rescue operations
Whether a prediction model is too optimistic or too pessimistic may be equally problematic in view of search and rescue operations. In practice, the resources available for search and rescue operations are limited and only a given  area can be covered, although the shape of the area (centre and eccentricity in the case of an ellipse) may not influence the cost significantly. Thus, rather than looking at the size of the search area as estimated from the ensemble model prediction, the search-and-rescue operation can be posed as follows: for a given area to be searched, which model forecast gives the ellipse that is most likely to contain the object?
The ensemble forecast provides the expected position, B(t), and the anisotropy, R(t), of the ellipse as defined previously, but the ellipse area is left free to grow homothetically from 1 to 3000 km 2 . The POC increases then accordingly as the observed buoy position is more and more likely to fall within the ellipse. The dependency between the search area and the associated POC defines the so-called selectivity curve, which makes a straightforward model comparison possible: the higher the selectivity curve, the better the model's ability to locate the searched object. The selectivity curves also allow an immediate evaluation of the rate at which predictive skill is lost as a function of time.
For each time t 0 + t, with t ∈ {12, 24, 36, 48, . . ., 10 × 24} h, we compute the POC corresponding to search areas ranging from 1 to 3000 km 2 for both models and seasons. Results from neXtSIM (solid lines) and FD (dashed lines) are shown at t 0 + 1, 2, 3 and 7 days in Fig. 19. In winter, for a given area, the POCs from neXtSIM are almost always above those from FD except in two cases: at t 0 + 1 day for search areas larger than 100 km 2 and at t 0 + 2 days for search areas larger than 500 km 2 . If we neglect the anisotropy for these cases (i.e. consider circular search areas), the POCs from neXtSIM become larger than FD. This indicates that the strong anisotropy in neXtSIM is more a disadvantage for small time horizon and large search areas in this experiment. Otherwise for smaller areas, larger time horizons or in summer, considering circular or ellipsoidal search areas makes no difference (not shown).
As long as the drift is longer than 3 days, the selectivity curves of neXtSIM are systematically above FD. However, in summer, for any time horizon and any POC, the results are very similar with a faint advantage to neXtSIM. When comparing to POCs of ellipses centred on forecasts from a deterministic neXtSIM run (not shown), the results are identical in winter and poorer with the deterministic run in the summer. For both periods and both models, all curves exhibit a sigmoid shape with an inflexion point, the position of which depends on the time horizon (higher POC and larger search areas for longer drift duration). For a 7-day drift in winter and a POC equal to 0.5, the area is around 300 km 2 with neXtSIM, while it reaches 1000 km 2 in FD. In the summer, a larger area is necessary to obtain the same POC for both models. For a given search area, the gap between the POCs from neXtSIM and FD seems independent of the drift duration in summer, whereas in winter it increases with the time prediction horizon. It is interesting to note the lowermost value of the POC for small areas in winter, which remains above 0.1 for neXtSIM. This could be a consequence of the capability of neXtSIM to simulate immobile ice, while the FD ice is always in motion with the winds and currents.
How do the different models perform for different forecast time (i.e. drift duration)? To answer this question, we study the time evolution of the difference between the neXtSIM and FD POCs: when this difference is positive (negative), neXtSIM (FD) is outperforming FD (neXtSIM). The POC for both models is evaluated for a fixed search area -a vertical section across the selectivity curves -equal to 50 km 2 www.the-cryosphere.net/12/935/2018/ The Cryosphere, 12, 935-953, 2018 in winter and 175 km 2 in summer, and the results are shown in Fig. 20. The chosen values of the search areas, 50 and 175 km 2 , correspond to the mean ellipse areas based on the ensemble spread from neXtSIM after 3 days, averaged over the IABP dataset (see Fig. 17), respectively, in winter and in summer. Figure 20 reveals that, after 2 days in the winter, the POC of neXtSIM is larger by about 0.2 than the POC of FD; most remarkably, such a substantially improved skill is then maintained almost stably up to the last day of simulation (10 days). During the summer, the POC of neXtSIM is also generally higher than the one of FD, but the difference is half of the one observed in the winter. Furthermore, after the third day, the difference between the two models decreases to vanish completely between days 8 and 9. The fact that most of the superiority of neXtSIM is found during winter is logical and should be no surprise given that during in the summer the ice mechanics of the two models are similar.
The negative values for lead time shorter than 1 day in winter are again likely caused by the initialisation of the neXtSIM ensemble and another reason to constrain its initial anisotropy to observations.

Discussions and conclusions
The ensemble model sensitivity experiment carried out with neXtSIM and with an FD model reveals the prominent role of the rheology, which marks the key difference between the two models. On average over the whole Arctic neXtSIM is less sensitive to the wind perturbations than the FD, although large seasonal and regional differences are observed. This is exemplified by the imprint of the ice thickness field in the ensemble spread from neXtSIM and the much smaller sensitivity of neXtSIM in winter than summer, in contrast to the FD model (Figs. 6 and 8). Both aspects point clearly to the role of the rheology, which accounts for the ice thickness and compactness. This behaviour should be expected to hold also for other sea ice rheologies than the EB.
The two models have been tuned on different observations of ice, seen as in free drift by each model, so that the different performances originate by the differences in the resolved model physics at their best performance. The diffusion regimes of neXtSIM and FD are very different in winter (Fig. 13): the offset between the curves indicating differences of sensitivity and the slopes indicating different rates of increase and thus sea ice diffusivity. The expected differences between summer and winter are only represented when the rheology is turned on.
Due to the dispersive properties of the sea ice, the shape of the ensemble of simulated buoys positions is generally anisotropic. Such anisotropy is a signature of the underlying mechanism that drives the dispersion of the members, which is the shear deformation of the ice cover along active faults/fractures in the ice. This mechanism is missing in the absence of rheology (like in the FD model) and represents a clear strength in principle for the EB rheology in neXtSIM, although with the present ensemble initialisation it did not prove to be a practical advantage. Other rheological models, such as the elasto-viscous plastic model, also present some degree of anisotropy (Bertino et al., 2015), although the two models have not been compared in the same conditions.
The performance of the two models differs significantly when forecasting the trajectories of IABP buoys. The ensemble mean position errors are larger in the summer (5 km after 1 day and 12.5 km after 3 days of drift for neXtSIM, about 16 % below the FD results) and consistent with the values reported by Schweiger and Zhang (2015) (RMSEs of 6.3 and 14 km, respectively, but using different time periods). The corresponding errors are smaller in winter, especially for neXtSIM (25 % smaller than FD), and down to 4 km for a 1-day drift and 7.5 km for a 3-day drift. These values seem competitive compared to the year-round average RMSE of 5.1 km per day in the TOPAZ4 reanalysis (Xie et al., 2017), even though the ice drift measurements are assimilated in TOPAZ4 (Sakov et al., 2012). The RMSEs of the free-drift model in Grumbine (2003) also seem to be higher than 5 km per day.
The model sensitivity to wind perturbations has been evaluated, yielding (for 10 days of drift) a spread from 5 to 10 km, for winter and summer, respectively, but this is smaller than the corresponding errors (15 km from the barycentre to the observations in Fig. 15). Still, since the diffusion regime is respected (at least in the winter), we are confident that the spread simulated by the model is physically consistent. However, other methods for perturbing the winds should be tested to remove the super-diffusive behaviour in summer.
To further improve the spread-error relationship, alternative sources of errors should be considered such as model initial conditions and forcings (ice thickness, concentrations, damage, ocean currents). Since the errors are increasing faster in the first days of the simulations, the more likely source of local and short-term errors lies in the position and orientation of the sea ice fracture network, which is left unconstrained in any of the experiments presented here.
Although we would expect an increase of the ensemble spread if the ice thickness, concentrations and ocean currents had been taken into account in the ensemble initialisation, we do not believe it would lead to a much larger spread, especially in the winter. We suggest instead that, in the perspective of efficient sea ice forecasting, major efforts should be directed toward assimilating the observed fractures (as of satellite images). The assimilation of fracture (as objects rather than quantitative observations) represents a priori a challenging avenue in terms of data assimilation, which traditionally deals with quantitative scalar or vector observations; however, we envision that the damage variable in neXtSIM, showing localised features, can be constrained quantitatively to deformation rates as derived from observed high-resolution ice motions and serve as "object assimilation".
The Cryosphere, 12, 935-953, 2018 www.the-cryosphere.net/12/935/2018/ In spite of the biases, the selectivity curves indicate that a probabilistic forecast using neXtSIM is largely more skilful than the traditional free-drift model, and it has the larger potential for practical use in search and rescue operations on sea ice. Since the Arctic is not easily accessible, forecast horizons of 5 to 10 days are probably the most relevant for logistical reasons. On those timescales, the differences of POC shown in Fig. 20 indicate that the free-drift model gives a poorer information in winter because of the biases in the central forecast location and the lack of anisotropy, while in the summer the use of an EB rheology is only marginally advantageous. The comparison of deterministic versus probabilistic forecast gives, as expected, an advantage to the average of the probabilistic forecast, although it is rather small and surprisingly more important in the summer, although the model non-linearities are stronger in the winter.
The physical consistency of the ensemble sensitivities is a necessary condition to the success of ensemble-based data assimilation methods (Evensen, 2009), which constitutes one follow-up research direction the authors are currently considering. Combining the modelling and physical novelty of neXtSIM with modern observations of the Arctic is seen as a major asset for forecast and reanalysis applications.
Besides the potential use of observations of fractures, as mentioned above, which is indeed another unique advantage of models such as neXtSIM, ice drift data are also crucial. Observations of ice drift are still seldom used for data assimilation, and when it is the case, the success is limited by the lack of sensitivity of the sea ice model (see e.g. Sakov et al., 2012). Nevertheless, the main fundamental issue related to the use of data assimilation, and particularly ensemble-based methods, is related to the nature of the Lagrangian mesh of neXtSIM, which also includes the possibility of re-meshing (Rampal et al., 2016b). This feature, while essential to the skill of the model in describing the mechanics of the sea ice with great details, represents a challenge in developing compatible data assimilation schemes, as the dimension of the state space can change over time when these re-meshing occur. This problem has recently attracted attention in the data assimilation research community (see e.g. Bonan et al., 2017;Guider et al., 2017;Carrassi et al., 2017) and it is also a main area of ongoing investigation of the authors, following the present study.
Author contributions. The sensitivity analysis was implemented and performed by MR. The results were analysed by MR, PR, AC and LB. The manuscript was written by MR and PR, then reviewed and improved with inputs from all authors.