Probabilistic forecasts of sea ice trajectories in the Arctic: impact of uncertainties in surface wind and ice cohesion

We study the response of the Lagrangian sea ice model neXtSIM to the uncertainty in the sea surface wind and sea ice cohesion. The ice mechanics in neXtSIM is based on a brittle-like rheological framework. The study considers short-term ensemble forecasts of the Arctic sea ice from January to April 2008. Ensembles are generated by perturbing the wind inputs and ice cohesion field both separately and jointly. The resulting uncertainty in the probabilistic forecasts is evaluated statistically based on the analysis of Lagrangian sea ice trajectories as sampled by virtual drifters seeded in the model to cover the Arctic Ocean and using metrics borrowed from the search-and-rescue literature. The comparison among the different ensembles indicates that wind perturbations dominate the forecast uncertainty i.e. the absolute spread of the ensemble, while the inhomogeneities in the ice cohesion field significantly increase the degree of anisotropy in the spread i.e. trajectories drift differently in different directions. We suggest that in order to get a full flavor of uncertainties in a sea ice model with brittle-like rheologies, to predict sea ice drift and trajectories, one should consider using ensemble-based simulations where both wind forcing and sea ice cohesion are perturbed.


Introduction
Sea ice covering the polar oceans is an important component of the Earth System. The dramatic changes of sea ice extent and volume in the Arctic have been regularly reported in the recent decades [1][2][3]. It is therefore crucial to understand the new state and characteristics of the Arctic sea ice cover and how it impacts the regional and global weather and climate [4]. Moreover, reliable sea ice forecasting systems are demanded for both operational and academic purposes [5]. For instance, the thinner sea ice cover offers opportunities in exploiting trans-Arctic shipping routes but its faster dynamics also challenge the safety of operations [6].
One specific and important aim of sea ice models is to represent small-scale dynamics such as the formation of leads and ridges, together with the large-scale drift patterns of big ice plates and small ice floes. In order to achieve this goal, numerical models consider a momentum equation including a specific term involving the internal stresses that accounts for the rheology of the sea ice material. Several rheologies for sea ice have been proposed for being used in continuum models, these include the Viscous-Plastic (VP, [7]), the Elastic-Plastic-Anisotropic (EPA, [8]), the Elasto-Brittle (EB [9]) and the Maxwell-Elasto-Brittle (MEB, [10]). In general, and particularly in view of the different rheologies, evaluating and calibrating models against field/laboratory measurements is necessary to improve their forecasting skill, as well as to identify and quantify the source of uncertainties.
Sea ice motion in the central Arctic is mainly related to the geostrophic winds [11]. However, the uncertainties in the atmospheric reanalysis in the Arctic are higher than that at the mid-latitudes, and observations are insufficient to estimate the statistical characteristics (scale, amplitudes) of the errors. Rabatel et al. [12] investigated the sensitivity of sea ice drift simulated by the Lagrangian sea ice model neXtSIM with the EB rheology to the uncertainties in the surface wind forcing via an ensemble of simulations obtained by applying perturbations to the wind fields. Predefined statistics of ensemble results were assessed and compared to the observed drifting buoy trajectories of the International Arctic Buoy Program (IABP). They concluded that the surface wind forcing accuracy as well as sea ice rheology in regions of highly compact ice cover are important in the probabilistic forecast skill with regard to simulated sea ice trajectories.
Sea ice cohesion, an intrinsic mechanical property of the material, is setting the local resistance of sea ice to pure shear deformation until break up. It is one of the critical parameters used in the brittle class of rheological frameworks (e.g. EB and MEB), which use a Coulomb-type failure envelope as a threshold criterion to weaken the mechanical strength of sea ice. In such rheological framework, this is achieved by increasing of the local, so-called, damage level of the sea ice, e.g. [13]. Because a proper value of cohesion to be used in models should depend on the spatial scale and local inhomogeneities in the ice microstructure [14] that are not explicitly represented in models, this mechanical parameter is commonly assumed constant in time and homogeneous in space. Yet, Bouillon and Rampal [15] showed that the value of cohesion used in neXtSIM significantly affects the properties of the simulated sea ice deformation patterns in both time and space, and therefore should have an impact on the simulation of sea ice trajectories.
Recently, the rheology implemented in neXtSIM has been updated from EB to MEB [13]. The MEB includes extra viscous mechanisms in the stress-strain relation, which can be reduced to the EB if the viscous relaxation time scale of the stresses for undamaged sea ice is taken sufficiently large. neXtSIM-MEB shows remarkable capabilities at reproducing the observed characteristics of sea ice kinematics and dynamics [16], in particular the space-time coupling of sea ice deformation scaling invariances [13] -a property never reproduced by a sea ice model before.
Therefore, using this most recent version of neXtSIM with the MEB rheology, the present study extends the previous work of [12], this time to assess and compare the model response to uncertainties in the surface wind forcing as well as in the sea ice cohesion, and assesses its predictive skill in terms of the forecast of sea ice trajectory forecast. Section 2 presents the methodology, which is based on the analysis of the ensembles of virtual drifter trajectories. Section 3 describes the experiment setup and how we generated the different sets of ensembles used in this study. Section 4 presents the results of the ensemble forecasts. Discussion is given in Section 5. Section 6 gives the conclusions and perspectives for the ensemble forecasting study.

Methodology
The Lagrangian trajectory simulations adhere to [12]: ensemble forecasts are conducted with multiple virtual drifters that are seeded in the model at the location and time of interest. The successive positions of the virtual drifters are updated at each time step according to the simulated velocity field, leading to a set of Lagrangian trajectories; the model response to uncertainty sources is estimated statistically using these trajectories. The Euclidean distance within the ensemble of drifters is defined as where ‖•‖ is the L-2 norm operator, indicates the -th ensemble member and N is the ensemble size. # ( $ , $ , ) is the drifter position of the -th ensemble member at time , with $ being the initial position where the drifter is deployed at time $ . ( ) is the barycenter -the mean of the ensemble positions. # ( ) is the distance between the drifter position of the -th ensemble member and the barycenter ( ). The ensemble spread is calculated as the standard deviation of # ( ), denoted as & . We define the position error, ( ), as the difference between ( ) and a reference drifter position ( ) at time as where ( ) could be either from a reference numerical simulation or observations, and ‖ ‖ indicates the distance between the ensemble mean and the reference drifter position. Taking the direction from the origin $ to ( ) as reference, we define the parallel and perpendicular component of ( ) onto the orthonormal basis centered on ( ) , denoted as ∥ ( ) and ( ( ) , respectively, describing the advection and diffusion of the virtual drifters in ensemble predictions. The forecast uncertainty in [12] is described by an anisotropic search ellipse, which is defined by the smallest ellipse encompassing all ensemble drifter positions, centered at ( ) and with its long axis pointed to the origin $ . We still define a region of uncertainty here but use a slightly different definition to account for anisotropy in the way how the ensemble of virtual drifters evolves. We define an anisotropic search ellipse that fits the ensemble of drifter positions using bivariate Gaussian Mixture distribution. The resulting cumulative density function gives the center of the ellipse, and the ellipse size is determined to include 99% of the probability density. Thus, the area of the ellipse and the anisotropy are defined as = and = / , (4) where and are the semi-major and semi-minor axes of the ellipse, respectively. Our definition implies that ellipses are smaller than that defined in [12] and are equal whenever the segment between $ and ( ) overlaps with the long axis.
For the sake of clarity, we denote the spatial average over all buoys of a general quantity as where is the total number of buoys, and # is the position of the i-th buoy. Similarly, we denote the average of over all the periods ( + = 13) as to even out the fluctuations of weekly environmental conditions, where subscript j indicates the j-th period.

The sea ice model neXtSIM
neXtSIM is a full dynamical-thermodynamical Lagrangian sea ice model that has been developed with the aim of better simulating sea ice dynamics and sea ice trajectories in particular. For instance, Rampal et al. [17] and Rampal et al. [18] demonstrated the realism of neXtSIM including the EB rheology with respect to the simulation of large scale sea ice drift and Lagrangian diffusion properties, sea ice cover thickness and extent, as well as the spatial scaling of sea ice deformation against SAR (Synthetic Aperture Radar) observations. The model has recently undergone major upgrades, among which the most relevant for the present study is the change in the rheology from EB to MEB. A detailed description of the model can be found in [13]. Another important development in the model is the enhanced computational efficiency, achieved thanks to the parallelization of the code [19].

Ocean and atmosphere forcing fields
The model is used in a stand-alone configuration, driven by ocean and atmospheric reanalysis products. The ocean forcing comes from TOPAZ4, the latest version of a coupled ocean-sea ice data assimilation system covering the North Atlantic and Arctic Oceans [20]. TOPAZ4 is based on the Hybrid Coordinate Ocean Model (HYCOM) and assimilates both ocean and sea ice observations using the ensemble Kalman filter (EnKF) [21]. The ocean forcing provided by TOPAZ4 includes the sea surface height, the current velocity at 30 m depth, the sea surface temperature and salinity, all given as daily mean values in average horizontal resolution of 12.5 km [17].
The atmospheric forcing comes from the Arctic System Reanalysis, (ASR, [22]) at a horizontal resolution of 30 km with 3 hours frequency. The variables used to force the model are 10 m wind velocity, 2 m temperature and mixing ratio, the mean sea level pressure, total precipitation and the fraction of that which is snow, as well as the incoming short-wave and long-wave radiation (see also [23]).
neXtSIM is spatially discretized using the finite element method on a triangular mesh which is adaptive in time, meaning that it is automatically regenerated when and where triangles are too distorted using an efficient remeshing algorithm locally [17]. The nominal mesh resolution (defined as the mean of the squared root of triangles surface areas) used for the experiments is about 7.5 km. The forcing fields are interpolated onto the center of the triangular elements during the model integration.

Simulation setup
We study the sea ice drift and associated Lagrangian sea ice trajectories in the Pan-Arctic Ocean from 1 January to 28 April 2008. During winter the Arctic is mostly covered by sea ice (cf. Figure 1 left and middle panels), and its extent and volume are close to their annual maximum. This provides abundant data of and compact sea ice, thus allowing to study the impact of ice cohesion on the sea ice drift. Initial conditions for sea ice concentration and thickness come from TOPAZ4, whereas the sea ice damage is uniformly set to zero everywhere. The time domain for our simulations is split into 13 successive periods (i.e. 13 start-dates), separated by 9 days. For each period, 10 days long ensemble forecasts are conducted. Over such short time period, the simulated sea ice drift is mainly influenced by the sea ice dynamics while the effect of thermodynamics can be considered negligible. Initial states of those forecasts are provided by a deterministic run, which is done foremost from 1 st January 2007 to 28 th April 2008. The experimental settings are summarized in Table 1. The ensemble runs are carried out in parallel and each 10-day run takes ~0.5-hour wallclock time / ~16 CPU hours on 32 processors of a Lenovo supercomputer. Note that we treat the model uncertainty from the ensemble forecasts that accumulate over time, which is conceptually similar to the forecast errors accumulating during an operational forecast.

Lagrangian trajectories
Virtual drifters are seeded in the model at the initial of each run. We use three sets of drifter trajectories. A first set for which the drifters are initially seeded at the same locations as the IABP 1 (International Arctic Buoy Program) buoys. It is used to evaluate the model skill with respect to the simulation of observed sea ice trajectories (section 4.1). A second set for which the initial positions of the drifters are regularly spaced by 50 km and covering the Arctic Ocean. This set is used to calculate model state statistics (section 4.2). The third and last set is similar to the second except that drifters' initial locations are in this case matching with the OSISAF ice drift dataset 2 (which provides estimated drift vectors with a distance of 62.5 km between them, [24]) and the drifters being redeployed every 2 days in order to be compared in a consistent manner to the OSISAF observations. This later set is used to calibrate the model parameters (section 3.5).
In Figure 1, the right panel shows the trajectories of the IABP buoys during the time period analyzed in this study (red tracks) and the initial layout of the regularly-spaced virtual drifters distant by 50 km (green dots). These virtual drifters are located further than 100 km away from the nearest coast as the same experimental choice was done by e.g. [12,16], in order to prevent the potential bias in the statistics due to the extremely anisotropic drifts present in the vicinity of coastlines. For the same reason, we only account for IABP buoys locating 100 km away from the coast in the following analyses. A brief summary of the ice and wind conditions in the area where the IABP buoys are drifting is given below. The area is mostly covered by sea ice of higher concentration than 95%. The ice thickness increases from 1.7 m in January to 2.2 m in March and then remains as thick until May, following the normal seasonal cycle of ice thickness. The winds are shifting directions and a particular episode of strong winds towards Greenland occurs in mid-March.

Air drag coefficient optimization
When using a forced -stand-alone-sea ice modeling system as we do here, the momentum fluxes from the wind and ocean currents to the sea ice need to be calibrated. This can be achieved by tuning the value of the air and/or water drag coefficients denoted . and / , respectively (see Table  2). Following the calibration method presented in [17], the air drag coefficient is optimized by comparing the simulated sea ice drift with observations from the OSISAF dataset [25] only where sea ice is in "free drift" while keeping the default value of the ocean drag coefficient unchanged. We identify the free-drift events when and where the simulated sea ice drift obtained when the sea ice rheology is activated differs by less than 10% of that obtained when it is not. Another criterion we use for the identification is on the sea ice speed, which has to be sufficiently large i.e. between 7 and 19 km/day as used in [17]. The selected drift velocity vectors are compared against the OSISAF drift vectors and for the two components U and V separately. We calculate the correlation coefficient and the Root Mean Square Error (RMSE) for each of the component.
The results obtained from a set of simulations using values of air drag coefficient ranging from 0.003 to 0.008 are presented in Figure 2. We find an optimal value of . about 0.0055, which is in between previously optimized values of 0.0051 and 0.0076 reported by [12] and [17]. We note that in these latter studies, the authors were using the same atmospheric forcing dataset but the neXtSIM-EB model instead of neXtSIM-MEB.

Ensemble experiments by perturbing wind/cohesion sources
The impact of sources of model errors is estimated using an ensemble of simulations with identical settings. Each ensemble member differs from the others by the application of a random perturbation process to two sources of model errors: surface wind and sea ice cohesion. Note that the two have an intrinsic different nature since the former is an external forcing while the latter is a parameter of the sea ice model. However, uncertainties on either or both of them contribute to the forecast errors, and we aim here at quantifying their individual and joint effects.
The wind-perturbation process adopted in this study is identical to [12], which is inherited from the TOPAZ4 data assimilation system [20]. The perturbations are non-divergent time-correlated geostrophic wind fields with a decorrelation time scale of 2 days and a horizontal decorrelation length scale of 250 km. The accumulation of these random perturbations along all model simulations causes the ensemble members to diverge from each other.
We adjust the wind speed variance of the perturbation to 3 m 2 /s 2 to generate sufficient ensemble spread. The perturbation process is applied online to every input wind field. The cohesion being an intrinsic parameter of the sea ice model, its values are perturbed only once at the initial time and then kept constant. Homogeneous scalar cohesion values are used in a first set of sensitivity experiments, then made variable in space in a second set. We thus conduct a series of forecasts using different cohesion values (5.5, 11, 16.5, 22, 27.5, 33, 38.5, 44, 49.5 and 55 kPa), keeping the model configuration almost identical to the deterministic run mentioned above except for the optimal air drag coefficient . = 0.0055 determined above. We note that the cohesion is a scale dependent property of sea ice, and these values are scaled at 7.5 km mesh resolution from laboratory value using the same relationship as in [13]. Figure 4 shows that the spatially averaged errors ⟨‖ ‖⟩ ) , D ∥ E ) , and ⟨ ( ⟩ ) are the lowest when setting the ice cohesion between 20 and 40 kPa, which is consistent with the default -not robustly optimized -cohesion of 25 kPa used in [13] for their neXtSIM-MEB simulations.
To further investigate the effect of sea ice cohesion uncertainties on simulated sea ice trajectories, we conduct ensemble simulations using inhomogeneous initial cohesion fields, so-called perturbed cohesion fields hereafter. At the start of each ensemble simulation, the cohesion value is chosen randomly between 20 and 40 kPa in each element of the model mesh following a uniform distribution and without spatial correlation. Afterwards, the values are kept constant during the simulation except in case of remeshing: the cohesion values in the new elements created by the remeshing process are calculated as the mean of the cohesion of their nearest neighbor elements.   generate the ensembles. Each ensemble contains 20 members, and each member contains 10-day forecasts for the 13 successive periods. Note that in the Monte Carlo techniques, the estimates of ensemble prediction converge by increasing the ensemble size. In our study, the convergence is achieved when the ensemble size exceeds 20 (not shown). Experiment setup is recapped in Table 1. Values of the model parameters used for running the forecasts are given in Table 2.

Results
In this section, we compare the impacts of the perturbation methods. In Section 4.1, results of the experiments are presented by comparing the ensemble spread and bias regarding to the IABP buoys. Section 4.2 presents the results of the analysis of the ensemble statistics of the virtual drifter' trajectories for the whole Arctic region. Note that neXtSIM only tracks the positions of virtual drifters when the concentration in the model grid cell is greater than 15%. In the following analysis, we calculate the statistics of the simulated ensembles only from the drifter trajectories spanning the full forecast period of 10 days.

Comparison of the simulated sea ice trajectories against the IABP dataset
The 12-hourly forecast errors against the IABP buoy trajectories are presented in Figure 4(a). The top panel shows the mean and standard deviation of the spatially averaged errors ⟨‖ ‖⟩ ) as the solid lines and error bars respectively, marked by different colors, one for each experiment. The behavior of the three ensembles relative to each other are generally in agreement across successive forecasts: the horizontal variability range of ⟨‖ ‖⟩ ) in COHESION is generally the smallest, contained within that of WIND and/or JOINT that are similar to each other. This implies the errors are generally larger when winds are perturbed and more variable geographically. One exception occurs in the last period 18 -27 April when COHESION gives the largest errors ⟨‖ ‖⟩ ) . It could be related to the seasonal change of Arctic sea ice. The period from 13 -22 March also stands out as the difference between WIND and JOINT ⟨‖ ‖⟩ ) is larger than usual. This could be a reaction to a strong and horizontally uniform wind event in the area covering the IABP buoys during 11 -18 March, which is observed Figure 4(b) that gives the daily-average wind velocities over the IABP buoys positions from ASR wind dataset. This point is further discussed in Section 5.
The absolute error ⟨‖ ‖⟩ ) is a distance and does not separate the advective and diffusive features of the drift errors. Consistently with [12] the errors are also projected onto the directions parallel and perpendicular to the drift. The spatial averages of these error components, D ∥ E ) and ⟨ ( ⟩ ) , are shown in the middle and bottom panels of Figure 4(a), respectively. The error components ∥ and ( take positive/negative values and follow the relationship ‖ ‖ = P || 3 + ( 3 . A positive/negative ∥ indicates that the virtual drifters move faster/slower than the drifting buoy. And a positive/negative ( indicates that the virtual drifters move to the right/left side of a vector from the initial position to current position of the drifting buoy. Again, the errors from the three ensembles generally follow the same order across successive forecast runs. It also shows that ⟨‖ ‖⟩ ) has its largest contribution from its parallel component D ∥ E ) , even more so when the trajectories are erratic. The quantities in Figure 4 The results from WIND and JOINT are in general agreement with Figure 15 of [12]: the forecast tends to drift too fast and to the right of the observed trajectories on average, although the time series in Figure 4(a) shows the errors differ significantly from one forecast to another and are flow-dependent. The WIND and JOINT ensemble tend to show longer drifts while the COHESION ensemble drifts slightly more to the right of the trajectory. Wind perturbations are constructed to have zero-mean components in the x-and y-directions of the grid, however since the absolute velocity is not a linear function of the x-and y-components, the ensemble average velocity of the perturbed winds is larger than the original unperturbed winds. This may explain the longer drifts in the experiments using wind perturbations. The evolution of the perpendicular error (to the right, then to the center regarding to observations) does not appear like a recurrent feature in Figure 4(a) and is probably not significant.  A complementary indicator of the quality of the ensemble forecasts is the ratio of ensemble spread over the root-mean-square error (RMSE) between the ensemble mean forecast and the observation, i.e., ⟨ & ( )⟩ ) /⟨‖ ( )‖⟩ ) taken as a spatial average. A ratio less than 1 indicates an under estimation of errors in the ensemble forecasts compared to the actual position errors. As shown in Figure 5(a), the ratio in our case is most often below 0.5 but shows strong temporal variations, often in parallel between the three runs, and increases in particular briefly in mid-March. This means that the ensemble spread is generally too small to contain the real buoys, which is consistent with [12] and that the rapid changes of the weather influence the quality of the ensemble forecast in a similar way for all three perturbation setups. The spatially averaged spread ⟨ & ( )⟩ ) is given in Figure 5(b). The spread of WIND and JOINT are almost identical and larger than that of COHESION, excepted for the second week of March. We will come back to this in section 4.2. The phenomena of rapid increase of the spread in 11 -14 March in Figure 5

Assessment of model sensitivity
The sensitivity of forecast is assessed using the statistics of ensemble trajectories. We look at the uncertainties in the ensemble forecasts (WIND, COHESION, and JOINT) using the simulated drifters' trajectories which have been seeded evenly (50 km distance from each other) over the whole Arctic domain at the beginning of each forecast.
At any time and location, we calculate the area of the search ellipse that encompasses the drifter positions from all ensemble members as in Eq. (4). The temporal evolution of this quantity indicates the diffusion properties of the virtual drifters that are initialized at the same location. Figure 6 presents the spatially averaged ellipse area over time, ⟨ ( )⟩ ) . The mean and standard deviation of ⟨ ( )⟩ ) are indicated as the solid lines and the error bars respectively; colors showing different experiments. ⟨ ( )⟩ ) is further averaged over the 13 forecast periods as ⟨⟨ ( )⟩ ) ⟩ , , and displayed in the enclosed panel. In this analysis, data exceeding three times the median absolute drift are removed as outliers. Figure 6 shows that the ellipse area increases monotonically with time as expected. The estimates from WIND and JOINT are close to each other and are on average 4 times larger than in COHESION. Besides, the lower typical areas (one standard deviation below average) from WIND and JOINT are generally larger than the higher areas (one standard deviation above) from COHESION. This signifies that dominant forecast uncertainties are due to the perturbation of the wind in most locations, which will be further confirmed below. There are however exceptions to this rule as the second week of March the COHESION ellipses grow faster than in previous weeks and the larger ellipses can become larger than the smaller ellipses from both WIND and JOINT. The area of the WIND and JOINT ellipses also tend to diminish in April and overlap more often with COHESION than earlier. We calculate the temporal average of ( ) of the ellipses evolved from the same starting positions (i.e., green dots in Figure 1(right)), ⟨ ( )⟩ , . Figure 7 shows the spatial distribution of ⟨ ( )⟩ , at forecast horizons 3, 5, and 7 days for each experiment. ⟨ ( )⟩ , are presented at the starting positions. Throughout the temporal domain, the statistics of ellipse areas from WIND and JOINT are in agreement, with a growth rate of ellipse area of about 89 km 2 /day, while the ellipse area of COHESION is much smaller with a growth rate of about 22 km 2 /day. Note that the larger ellipse areas, in all experiments, are found near the ice edge in the Chukchi/Beaufort Sea to the west, as well as the Nansen Basin, Barents Sea and Kara Sea to the east. Strong and narrow currents in those regions may contribute to increase sea ice diffusion. The color scale is capped at 1000 km 2 to highlight regional differences, thus very large spreads saturate at the ice edge in the Fram Strait, related to strong local currents.
We further compare the forecast uncertainty generated by the wind perturbation and cohesion perturbation alone. ( ) is defined as the ratio of the intersection area between the wind-ellipse and the cohesion-ellipse over the area of the cohesion-ellipse at time t. Hence, = 1 indicates that the wind-ellipse includes completely the cohesion-ellipse, while = 0 indicates the two ellipses are entirely disjoint. Figure 8 shows that averaged over all virtual drifters is higher than 0.9 throughout the studied period, indicating that the ellipse due to cohesion perturbations is mostly included in wind-driven ellipses. This is the case over the entire Arctic with the exception of an approximately 100 km×100 km region near the coast between north of Greenland and the Canadian Archipelago (not shown): we argue that this may be related to the specific role of sea ice cohesion in areas of very thick ice. Further, we present the spatial average of the anisotropy against time ⟨ ( )⟩ ) in Figure 9. Results show that ⟨ 4567)85" ⟩ ) is nearly twice ⟨ 98": ⟩ ) but D ;58", E ) is only slightly larger than ⟨ 98": ⟩ ) , for all , so that time is dropped from the notation for clarity. This signifies that the impact of the cohesion perturbations on the sea ice shear deformation is predominant and substantial. The same quantity averaged also over all periods ⟨⟨ ( )⟩ ) ⟩ , , is shown in the enclosed panel. The spatialtemporal averaged anisotropy drops quickly after the first two days in all experiments, but ⟨⟨ 4567)85" ⟩ ) ⟩ , decreases at a slower rate than ⟨⟨ 98": ⟩ ) ⟩ , and RD ;58", E ) S , , which is consistent with Figure 11 in Rabatel et al. [12] obtained for their experiment with wind perturbations. The quantitative differences may however stem from the different definitions of the ellipses adopted here (see Section 2) and not only from the differences in the rheology. The decreasing behavior may be explained that ice floes (tracked by the virtual drifters in our ensemble runs) tend to move along the same initial fractures from = 0. With diverse ice-broken pattern due to the applied perturbation, significant different ice fractures could be developed among the ensemble members in the first 2 days. It leads to a more isotropic dispersion (smaller ) of ensemble drifters with respect to their barycenter [12]. Figure 10 shows the spatial distribution of the anisotropy averaged over the 13 periods at 7-th day from = 0, ⟨ (day = 7)⟩ , , for all the experiments. Similar pattern is observed in other days and is not shown. The maximum of the color scale is kept at 3 to better visualize the most central part of the Arctic basin, however the maximum anisotropy does reach 13 near the coasts. In all experiments, the anisotropy is lower in the central Arctic and higher near the Canadian-Greenland and Siberian coasts. The results of the WIND and JOINT experiments are also quantitatively in agreement, showing lower anisotropy than that of COHESION. This is due to the dominance of isotropic wind perturbations, which masks the anisotropy caused by the cohesion perturbations. Moreover, the distribution of anisotropy is opposite to the distribution of ellipse area shown in Figure 7. This pinpoints areas where the dispersion of ensemble drifter follows both small and narrow shapes. This may be related to the thick ice or to the presence of landfast sea ice in these regions of the Arctic. Both result in fewer ice fracturing, opening of leads, and therefore significantly slower sea ice drift and diffusion.

Discussion
In this study, we revisited the results from our precursor work [12] using the Lagrangian sea-ice model neXtSIM, this time with the updated Maxwell-Elasto-Brittle rheology of [10] and evaluated the impact of the assumed main intrinsic (ice cohesion) and extrinsic (surface winds) source of uncertainties on sea ice trajectory forecasts. Compared to the previous EB rheology, the MEB introduces an extra viscous relaxation of stresses mechanism in the damaged sea ice. It has an unprecedented capacity to reproduce the main characteristics of sea ice deformation [26], and therefore a positive impact on the overall sea ice drift properties and the simulation of ice trajectories.
The study domain is the Arctic Basin and the time period is from 1 st January to 28 th April 2008. This 'wintertime' period is characterized by the sea ice cover extending from coast-to-coast across the Arctic Basin. It is thus convenient to exhibit the role of cohesion on ice damage, an intrinsic parameter of the brittle-like MEB rheology, which plays a major role on the mechanical behavior of sea ice when the latter is packed. We conduct three types of ensemble forecasts. The first one is generated by perturbing the wind, the second by perturbing the ice cohesion, and the third one by perturbing jointly the wind and the cohesion. The sea ice drift is assessed through the analysis of Lagrangian sea ice trajectories simulated by the model, obtained by seeding and tracking the virtual drifters all over the Arctic region.
The results from the three ensemble experiments are first compared to the trajectories of the IABP buoys. The results of WIND generally agree with an earlier model version [12]. Both the errors and the uncertainties of dynamic sea ice drift are overall larger with the perturbation of winds than ice cohesion, according to the time-averaged forecast errors and spread (Figures 4, 5 and 6). The spread underestimates the errors in all cases, but even more so if only cohesion is perturbed. The ice cohesion is known to play an important role for ice drift in a complete ice cover. However, its effect reduces when ice cover is fragmented, e.g., broken due to the wind, in which case the uncertainty of ice drift is more sensitive to the uncertainty of winds. This is exemplified by the strong and uniform wind event from 11 to 18 March in the area where the IABP buoys were present. The ensemble spread increases rapidly from 11 -14 March in Figure 5(b), and unusually quick in the case of COHESION.
The effect is extended to the forecasts from 13 to 22 March that the spreads of WIND and JOINT increases quickly in the first few days. In contrast, the spread of COHESION grows slower because the ice cohesion is re-initialized at the beginning of each forecasting period, thus the wind effect on ice dynamic is limited. It explains the peak of ratio of spread over the forecast error in WIND and JOINT on 15 March since the position errors between the simulated drifters and real buoys are small. Because the model response is more sensitive to perturbing the wind than perturbing the cohesion, under the wind event, the difference is more obvious that the virtual drifters in WIND drift farther than that of COHESION regarding to the real buoys in Figure 4(a).
Further, we studied the sensitivity of the simulated diffusion and dispersion of drifter trajectories -characterized by the surface area and eccentricity of the ellipse encompassing the ensemble. We found that the spatial and temporal characteristics of the diffusion and dispersion in the WIND and COHESION experiment are different, with the WIND ellipses being much larger ( Figures 6 and 7), and therefore the JOINT experiment mostly resembled the WIND experiment. In addition, the ellipses from the COHESION experiment are generally contained into the corresponding ellipses from WIND experiment ( Figure 8). This confirms that the uncertainties in probabilistic forecasts arise mainly from the wind perturbations. Note however that the results are obtained for given choice of the wind and cohesion perturbations amplitude and could be sensitive to their value, although the order of the results and our conclusions would be unlikely to be reversed under a different choice.
Nevertheless, we show that there are a few times and locations where the cohesion is an important source of error and moreover its uncertainties systematically increases the eccentricity of the ellipses, i.e. the degree of anisotropy of the ensemble dispersion. In other words, the individual members tend to disperse along preferred directions (see Figure 9 and 10). This phenomenon can be explained by the ice deformations starting from locally low-cohesion regions leading to cracks in the ice field with inhomogeneous cohesion. Under further internal/external forces, neighboring cracks are connected to form longer linear fractures. As a result, drifters in these forecasts tend to move along the linear features with least resistance from ice internal stress, leading to higher anisotropy in the ensemble as compared to forecasts initialized with a homogeneous cohesion field. Overall, the anisotropy is on average lower in the JOINT experiment than in COHESION due to the applied isotropic wind perturbations. It is worth noting that increasing the ensemble size reduces but does not remove the anisotropy. In summary, we believe that the anisotropy is inherited from the intrinsic mechanical behavior of the sea ice, which is characterized by the formation of fractures, geometrical features that are themselves -by nature-anisotropic.
The reduction of anisotropy with time, observed in the small panel in Figure 9 could be explained as follows. All members of the ensemble initially share identical sea ice properties, including the sea ice damage. During the first few steps of the forecast run, all the drifters starting from the same position in an ensemble tend to drift in the same direction, along the preexisting linear kinematic features, although at different velocities. Thus, the ellipse generated from the ensemble drifter's positions at early times is highly anisotropic. As the simulation progresses, the ice damage level eventually increases under the effect of further break up events, leading to the formation of new -additional-linear kinematic features in other directions, also different from one ensemble member to another. The drift of the drifters in an ensemble becomes gradually more isotropic over time.

Conclusions and perspectives
The present study demonstrates that using the neXtSIM stand-alone model and applying wind perturbations identical to those in the TOPAZ4 HYCOM forecasting system [20] have a significant general impact on sea ice drift, and lead to a large spread of sea ice trajectories. Wind forcing should be considered as a primary source of uncertainty, not only for the sea ice drift, but also for sea ice concentration, thickness, damage and snow thickness on ice. Although to a lesser extent, uncertainty in ice cohesion also occasionally impacts the forecast skill by enhancing the preferential directions of the local drift as a result of the formation of fractures in the sea ice. Surprisingly, the addition of cohesion as a new source of uncertainty did not improve the ensemble forecasting skills, in particular not in the wake of the strong wind event of the 11 th -14 th March that made the model more sensitive to cohesion. There are two possible explanations to this: either the anisotropy brought by perturbations of cohesion did not align with the actual anisotropic features in the field and deteriorated on average the drift forecasts or the benefits of the perturbations were masked by the variable quality of the model forcing (winds and currents may have been more inaccurate during the passage of Arctic storms, thus annihilating the benefits of higher ensemble spread).
Future efforts should therefore focus on the assimilation of sea ice drift observations, which should reduce position errors for all tracers and help aligning the simulated with the observed anisotropic features. Further, errors in modeled sea ice thermodynamics will in the longer term affect the same tracers. Fortunately, observations of sea ice concentration are readily available (e.g., ESA SICCI [27]) and observations of sea ice thickness as well, although in the winter months only (e.g., CS2SMOS [28]). For an efficient multivariate assimilation of these observations, the sources of errors in model thermodynamics should be included in the above framework with perturbations of e.g., surface air temperature, snow precipitation and cloud cover. This will allow the assimilation of observed variables to update the unobserved ones, similarly to TOPAZ4 [29]. The present framework for ensemble simulations would then provide -in the data assimilation jargon -the ensemble forecast, while the analysis can be obtained by applying a generic ensemble data assimilation package.
Our future work is to implement the ensemble Kalman filter (EnKF, [21]) in the current neXtSIM stand-alone model. The EnKF has been successfully used in several geophysical systems to improve the forecasts [30]. One of the main challenges is to use EnKF on the Lagrangian framework used in neXtSIM, which is an adaptive triangular non-conservative mesh -i.e. whose dimension can also change in time. Aydoğdu et al. [31] developed a modification of the EnKF for moving nonconservative mesh models in one-dimension, and we are currently studying the extension of their methodology to two-dimensions and its application to neXtSIM. Results in this study suggest that, for the ensemble in the EnKF to maintain the spread of trajectories, it should be preferably generated by perturbing the surface wind forcing but should not exclude cohesion errors.