of ‘ Probabilistic forecast using a Lagrangian sea ice model : application for search and rescue operations

I am not really surprised to see that neXtSIM does better than the FD model (especially in winter). In fact, I would like to suggest a different way to present these results. Instead of saying you compare two models, you could say you study the impact of rheology on probabilistic forecasts. The title of your paper could be: Probabilistic forecast using a Lagrangian sea ice model: impact of rheology. This is just a suggestion.


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 the 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 felt 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 2006; Breivik and Allen, 2008;Melsom et al., 2012).
Our main research tool and object of study is the advanced 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 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) and Bouillon and Rampal (2015b)), as well as by the in situ measures 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, yet with different rheologies (e.g.visco-plastic (Coon et al., 1974;Hibler III, 1979), elasto-visco-plastic, (Hunke and Dukowicz, 1997), or Maxwell-elasto-brittle, (Dansereau, 2016)), 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, whereas this behaviour drastically changes when considering a large, continuous and undamaged solid plate.The second reason stands on the fact that surface wind velocity fields provided by atmospheric re-analyses contain large uncertainties in the Arctic.
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.Although these analyses did not use the full complexity of the present version of neXtSIM (in particular they did not include the thermodynamics, nor the re-meshing process), the impact of some mechanical parameters on the ice deformation can still be considered as valid.
This paper is organized as follow: 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 The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-200Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017.CC BY 4.0 License.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 conclusion are drawn in Sect. 5.

Generalities on the model neXtSIM
In this section, we provide a general description 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, it is not spatial homogeneous, and it 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 center of the mesh elements, and vectors defined at the vertices.The model is using a mechanical framework that has been developed recently (Girard et al. (2009) and Bouillon and Rampal (2015a)), and which is based on the Elasto-Brittle (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 thermodynamic, 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 refreezing, 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 time scales 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; A is the sea ice concentration; 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/source term.The evolution of sea ice velocity comes from the following sea ice momentum equation, integrated over the vertical, The Cryosphere Discuss., https://doi.org/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 with ν being the Poisson's ratio while E(A, d) 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.

Sensitivity analysis
The sensitivity, or uncertainty analysis, are performed in order to understand and to quantify the relative importance of different input sources in the outputs.More specifically, we explore the output space of the model and, after the definition of a region of interest for this output space, we identify which model inputs and which values and uncertainties of the inputs better explain the model results in the chosen region.This type of analysis is an usual important previous step to determine optimal sampling strategies for both probabilistic forecast and ensemble-based data assimilation methods (e.g.Evensen, 2009).

Methodology
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.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 ∀i ∈ {1, . . ., N } , 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), B(t) = N i=1 g i (x 0 , t 0 , t), at the same time t (see the top panel of Fig. 1).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 on the dependence on x 0 and t 0 , to simplify the notation.
Furthermore, we define a 2-dimensional time-dependent orthonormal basis, centred on B(t), and whose axes are the line connecting x 0 to B(t), and its perpendicular.The components/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. 1; they provide informations 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 hands, we compute basic, second-order, statistics.Let consider their means, µ r and µ b , g g g i (t) as our main quantities of interest in the analysis that follows.In particular, we use the standard deviations to compute the ratio that provides a measure of the anisotropy of the ensemble spread of the virtual buoys positions around the barycentre B of the  The wind forcing is taken from the Arctic System Reanalysis (ASR) (Bromwich et al., 2012).This reanalysis product provides wind speeds and direction every 3 hours and at a horizontal resolution of 30 km.For every 3-hourly wind field, we generate spatio-temporal correlated perturbations as described in Evensen (2003), and then add them to the basic wind field from the ASR.This procedure is identical to the one used to produce ensemble runs with the coupled ocean-sea ice model TOPAZ (Sakov et al., 2012;Melsom et al., 2012).The main advantage of the method is that the perturbed wind fields are keeping 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 almost 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 study, we used a decorrelation time-scale of 2 days, a horizontal decorrelation length scale of 250 km, the and the wind speed variance as equal to 1 m 2 s −2 .These values are identical to those used in Sakov et al. (2012) except for a reduced wind speed variance (6 times smaller).
Although the ensemble average of the perturbations is equal by construction to the original wind directions provided by the ASR reanalysis, 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 reanalysis following a method described 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 compare best with the observed ice drift from the OSISAF dataset (Lavergne and Eastwood, 2015).The optimisation is carried out at all times but limited to the region where the ice is in free-drift (see Fig. 2).
Figure 3 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 optimized value for the drag coefficient C a = 5.1 × 10 −3 , lower than the one used in Rampal et al. (2016b) (7.6 × 10 −3 ).
The ocean forcing comes from the TOPAZ4 reanalysis (Sakov et al., 2012).TOPAZ4 is a coupled ocean-sea ice system combined with a state-of-art ensemble Kalman filter data assimilation scheme assimilating both ocean and sea ice observations.In our simulations, we used the 30 m depth currents, 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, respectively from 1 January to 10 May and from 1 July to 20 September, representative of the winter and summer conditions.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), becomes very large and dominant.As a consequence, the ice drift is (on average) much reduced.During the summer period, on the other hand, the ice concentration is lower and the ice pack does not generally reach the coasts, the ice internal stresses are close or equal to zero, and the ice drift closer to a free-drift state (see text below).We remark however that the wind field perturbations are generated using the same procedure aforementioned, for both the winter and the summer, and have thus the same spatial and temporal properties.
We ran a total of 13, in the winter, and 8, in the summer, simulations for successive non-overlapping 10-days long periods.
Limiting the length of the simulations to 10 days ensures that the thermodynamical effect (increase in sea ice concentration or thickness) on the drift can reasonably be considered negligible along the track.The starting positions are spaced by 100 km and cover the domain as displayed in Fig. 4. We ran an ensemble of 12 members, each of them forced by the perturbed wind dataset 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 are thus confident that N = 12 suffices to our purposes.From these ensemble runs we simulated a total of 8000 virtual buoy trajectories over the winter season, and around 3200 trajectories over the summer season.This dataset was used to run the analyses described in Sect. 4 and presented at the 19 th EGU General Assembly (Rabatel et al., 2017).
As already stated, we compared neXtSIM with the so-called free-drift model, hereafter referred to as FD, so that all simulations that follow have been carried out for the two models.neXtSIM, Eq. ( 2) with all terms in 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 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 therefore mimics the drift of a buoy at the surface of the ocean.
Figures 5 and 7 show maps of mean drifting distance and spread (see the definitions of µ b and µ r in Sect.3.1) of the virtual buoys after t = 10 days, averaged over the 13 (winter) and 8 (summer) successive simulations.Similar results are obtained for different time t ∈ [0, 10] days (not shown).The pixels on the maps correspond to boxes of 100 × 100 km centred on the initial positions x 0 where the virtual buoys have been deployed at t 0 .
Figures 6 and 8 are the counterparts of Figs. 5 and 7 and show the average wind speed (left panel) and ice thickness (right panel) for winter (Fig. 6) and summer (Fig. 8) respectively.Note that both figures are relative to neXtSIM, but the free-drift wind speed is identical (same perturbations) and the ice thickness geographical pattern very similar; we have thus omitted to display them to avoid redundancy.
From Fig. 5 and 7 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 completely absent in FD for which the obtained fields appear almost random.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. 6 and 8 for winter and summer respectively).
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.On another hand, in the surrounding Seas (i.e.Beaufort, Bering, Chucky, Kara and Barents Seas from West to East), where the ice is thinner and the winds stronger, the mean of µ r and µ b are larger.Note that in summer these correlations or anti-correlations are even stronger, for example between the means of µ b and the ice thickness (see Figs. 7 and 8).
For FD, and for both winter and summer, the mean values of µ r are less correlated to the thickness field than in neXtSIM, In both winter and summer, the response to wind perturbations is overall lower by 35% in neXtSIM than in FD.This can be attributed to the fact that the ice rheology is taken into account in neXtSIM, thus acting as an additional filter on the momentum transferred from the wind to the ice.In more details, it is interesting to note that the magnitude of the impact of the ice rheology is different whether we look at the drifted distance by advection r or the spread distance by diffusion b and during the winter or the summer.Averaged over the winter, µ r (t) D and µ b (t) D are respectively 35% and 63% lower in neXtSIM than in FD at t = 10 days, whereas over the summer, µ r (t) D and µ b (t) D are respectively 14% and 39% 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 keeps the members closer to each other.During summer, the ice is generally much less packed and the physical/dynamical differences between neXtSIM and FD have a lower impact.
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 concentration over the whole domain in winter is about 0.99 while it drops to 0.83 in the summer.In neXtSIM, this strongly influence 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 buoy's drift in neXtSIM closer to the one 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. In neXtSIM though, the ice rheology is acting in increasing 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 temporal evolution over 10 days of the Arctic-averaged ratio R, Eq. ( 10), that defines the degree of anisotropy of the ensemble spread (1: isotropic; > 1: anisotropic).We observe that R is very close to 1 and relatively constant over time in the FD model.On another 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 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 a behaviour cannot be reproduced by the FD.Note also that R is as large as 2 within day 1 and 2 for neXtSIM in the winter, and it then monotonically decrease for t > 2, but it still remains 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 wind direction away from the origin, but it then starts to break and depart from the barycentre.

Spatial and temporal properties of the ensemble spread
In Fig. 11, 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 and thus the anisotropy of the ice drift, for neXtSIM exhibit marked spatial correlations that are almost absent in FD.
Another important characterisation of the ensemble spread evolution can be set by looking at the variance of the distance 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's decomposition m = m + m .Figure 12 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 follow two distinct diffusion regimes, one for small time t Γ and one for large time t Γ where Γ is the so-called integral time scale (Taylor, 1921), which is about 1.5 days for sea ice according to Rampal et al. (2009).In neXtSIM, the first regime we found 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 two regimes are reproduced by neXtSIM in our winter simulations.These results are in agreement with the 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), and shows that our experimental setup based on ensemble simulations forced by perturbed winds does not alter the capability of the neXtSIM model to reproduce the properties of sea ice diffusion that was reported recently in Rampal et al. (2016a).One should note that for FD, for both periods, and for neXtSIM during the summer, b i 2 D ∼ t 1.15 for t Γ meaning that the simulated growth of the ensemble is super-diffusive and in disagreement with the observations.

Predictive skills of neXtSIM and of the FD models
We evaluate here how well neXtSIM and FD models are able to forecast real trajectories.As a benchmark, we compare the ensemble runs from each model 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, and are displayed in Fig. 13; the positions of IABP buoys are known every 12 hours.
As a metric for the models skill inter-comparison, 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. 15).The components of e(t) onto the orthonormal basis centred on O (see Sect. 3.1,Fig. 1 and 15), read e (t) and e ⊥ (t). to 28.5 km at day 10 compared to 14.6 km for neXtSIM.
As already deduced from the results in the previous section, the mechanics underlying of the ice drift in neXtSIM and FD get very 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 4 km after 4 days (see the left panel in Fig. 14).The central panel in Fig. 14 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.Nevertheless, the biases in neXtSIM are all as small as one third of the corresponding ones in the FD model.The general positive biases betray a too fast drift in the direction along the ensemble mean drift compared to the observations.Finally, the right panel in Fig. 14 also reveals a bias of the error in the direction across the ensemble mean drift, yet substantially weaker than in the previous case.For FD, e ⊥ still being positive for both periods, corresponds to a drift too far to the right in the observations.This bias should be explained by the fact the Coriolis effect depends on the velocity, which is generally higher in the FD model; the larger the velocity the larger the Coriolis effect.With neXtSIM, these biases are much weaker in both periods.Overall, the performances are best for neXtSIM, especially in winter.
In Hackett et al. (2006); 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 operations time and resources affordable and efficient.
The forecast system should therefore ideally yield a high probability to find 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 to find the object inside the search area, is referred to as the probability of containment, POC, and reads POC is proportional to the ratio of the size of the search area to the square forecast error.A small fore cast 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. 13), 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 ensemble of buoys of 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.2, and illustrated in Fig. 15.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 is increasing with the ensemble spread and contains 100% of the ensemble members.
Short 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 reference.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 a solid, held together by the ice rheology, in particular in high concentration areas, so that 20 The Cryosphere Discuss. 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;Maio et al., 2016).
Figure 16 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 300 km 2 after 2 days and it almost reach 1000 km 2 after 5 days.The search area in FD is about 10 times larger than in neXtSIM in the winter and 4 times larger in the summer.Therefore, even if the forecast errors are smaller in neXtSIM than in FD, its shrunk search areas lead to a smaller POC for neXtSIM than for the FD model (not shown): in practice the probabilistic forecast from neXtSIM is too optimistic, underestimates the uncertainties in the forecast, while the FD forecast overestimates them.

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 (center 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 posed as follows: for an equal area that can 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 km 2 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 possible a straightforward models comparison: the higher the selectivity curve, the better the model ability to locate the searched object.The selectivity curves allow as also for 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} hours, we compute the POC, Eq. ( 13), corresponding to search areas ranging from 1 km 2 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. 17.For a given area, the POCs from neXtSIM are almost always above those from FD except in two cases: in winter at t 0 +1 day and for search areas larger than 200 km 2 , and in summer at t 0 +1 day for search areas smaller than 7 km 2 .As long as the drift is longer than 2 days, the selectivity curves of neXtSIM are systematically above FD.For both periods and both models, all curves exhibit a sigmoid shape with an inflexion point, which position depends on the time horizon (higher POC and larger search areas for longer drift duration).For a 7-days drift in winter and a POC equal to 0.5, the area is smaller than 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 from 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 in winter and 175 km 2 in summer, and the results are shown in Fig. (18).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. 16) respectively in winter and in summer.Figure 18 reveals that, after 2 days in the winter, the POC of neXtSIM is larger of about 0.2 than the POC of FD; most remarkably, such a substantial improved skill is then maintained almost stably up to the last day of simulation ( 10days).During the summer, an 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 3 rd day, the difference between the two models decreases up to vanish completely between day 8 and 9.The fact that most of the superiority of neXtSIM over reveals during winter is, as stated in previous instances, in full agreement with the expectations, given that during the summer the ice mechanics in the two models is similar.

Discussions and Conclusions
The ensemble sensitivity experiment carried out with neXtSIM and with a 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 35% less sensitivity to the wind perturbations than the FD, albeit 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, contrarily to the FD model (Fig. 5 and 7).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 for other sea ice rheologies than the elasto-brittle.
The two models have been tuned on a common dataset of observations of ice in free drift, so that the different performances originate solely by the differences in the resolved model physics.The diffusion regimes of neXtSIM and FD are very different (Fig. 12): 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 and advantage of the elasto-brittle rheology in neXtSIM.
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 drift for neXtSIM, about 20% below the FD results), and consistent with the values reported by Schweiger and Zhang (2015) (RMS errors of 6.3 km and 14 km respectively, but using different time periods).The corresponding errors are smaller in winter, especially for neXtSIM (31% smaller than FD) and down to 4 km for a 1-day drift and 7.5 km for a 3-days drift.These values seem competitive compared to the year-round average RMS error 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 model sensitivity to winds has been evaluated, yielding (for 10-days 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. 14).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.Alternative sources of biases must be called such as, for example, other model inputs (thickness, con- 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 initialization, yet we do not believe it would lead to a much larger spread, especially in the winter.Still it is the wind forcing being the key player in the spread evolution.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 localized features, can be constrained quantitatively to deformation rates as derived from observed high-resolution ice motions and serve as "object assimilation".
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 time scale, the differences of POC shown in Fig. 18 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 a elasto-brittle rheology is only marginally advantageous.
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 of the 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 most the sea ice models (see, e.g.Sakov et al., 2012).
Nevertheless, the main fundamental issue related to the use of data assimilation, and particularly ensemble-based, methods, stands in the nature of the Lagrangian mesh of neXtSIM, which also include 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., 2016;Guider et al., 2017) and it is also a main area of on-going investigation of the authors, following the present study.

Figure 1 .
Figure 1.From a 12 bouquet of simulated trajectories of a virtual buoy drifting during 10 days of which only two of them, denoted i and j, are drawn, we represent the distances r, b (top) and the coordinates b , b ⊥ (bottom) for the virtual buoy i and j at time t. and the standard deviations, σ b and σ b ⊥ , of the components b and b ⊥ ,

Figure 2 .
Figure 2. Spatial distribution of the number of occurrences of free drift events between 1 January 2008 and 30 April 2008.These are the instances used for the optimisation of the air drag coefficient.

3. 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 4).

TheFigure 3 .
Figure 3. Scatter plots for the two components, x (left) and y (center), of the simulated (neXtSIM, x-axis) and observed drift (OSISAF dataset, y-axis) after the air drag optimisation procedure.The cumulative distribution of the ice velocity errors is shown in the rightmost panel.

Figure 4 .
Figure 4. Maps showing the Arctic domain considered for this study.The red lines are open boundaries, while the black coastlines are closed boundaries.The starting points of the ensemble trajectories simulated with the neXtSIM and FD models are represented by the blue crosses.

Figure 5 .Figure 6 .Figure 7 .Figure 8 .
Figure 5. Mean over the winter period of µr(t) and µ b (t) at t = 10 days.The calculated values are represented by coloured squares centred on the starting points x0 shown in Fig. 4.

Figure 9 Figure 9 .
Figure9shows 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 (seeSect.3.1).The PDFs of b and b ⊥ for the FD case are almost identical, thus we chose to only display one curve (black dashed line).The first aspect to remark from Fig.9is that all distributions are uni-modal and symmetric, suggesting that the 2D-shape of the ensemble is symmetric around its barycentre B. However, we notice that the ensemble is anisotropic in neXtSIM, that is, the distributions of b and b ⊥ differ substantially, whereas it is close to isotropic in FD.

(R > 1 )Figure 10 .
Figure 10.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).

Figure 11 .
Figure 11.Means over the winter (top panels) and summer (bottom panels) of R(t) at t = 10 days.The values are represented by coloured squares centred on the starting points x0.The spatio-temporal mean values (i.e.domain averaged and seasonally averaged) are calculated for each model and displayed above each panels for reference.

10 Figure 14 Figure 12 .Figure 14 .
Figure14shows the average module 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

TheFigure 15 .
Figure15.Illustration of the forecast error and the anisotropic search area.Blue dots represents the position of one member, while the barycentre of the ensemble (its mean) is B(t).The observation O(t) is in green and the forecast error is defined in Eq. (12).See text for definitions of the search ellipse and anisotropy ratio.

Figure 16 .
Figure16.Time evolution of the averaged ellipse areas for neXtSIM (solid lines) and FD (dashed lines) in winter (blue) and in summer (red).

10 23Figure 18 .
Figure18.Time evolution of the POC difference between neXtSIM and FD for a search area equal to 50 km 2 in winter (blue) and equal to 175 km 2 in summer (red).

The
Cryosphere Discuss., https://doi.org/10.5194/tc-2017-200Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017.CC BY 4.0 License.centrations, 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 not constrained at all in these experiments.
∈ 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 ], and 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 }