Sea ice diffusion in the Arctic ice pack: a comparison between observed buoy trajectories and the neXtSIM and TOPAZ-CICE sea ice models

Due to the increasing activity in Arctic, sea ice-ocean models are now frequently used to produce operational forecasts, for oil spill trajectory modelling and to assist in offshore operations planning. In this study we propose a method based on a Lagrangian diffusion analysis to evaluate the sea ice drift properties simulated by two sea ice models, TOPAZ-CICE and neXtSIM, used in two different sea ice-ocean systems developed for such applications. We compare their results to the 5 buoy trajectories of the International Arctic buoy Program (IABP) data set and we find that neXtSIM performs better than TOPAZ-CICE in simulating the mean and fluctuating sea ice velocities over the central Arctic in winter. Our analyses indicate that both TOPAZ-CICE and neXtSIM are able to simulate two distinct sea ice diffusion regimes depending on the time scale considered, similarly to what is predicted by the steady and homogenous turbulent flow theory. However, the basin-averaged 10 absolute diffusion computed from the analysis of drifters trajectories simulated with TOPAZ-CICE is almost twice as high as the value estimated from both the corresponding drifters trajectories simulated with neXtSIM and from observed buoy trajectories. Also, the mean Arctic pattern of absolute diffusion obtained from TOPAZ-CICE shows large differences from the one obtained from the observed buoy trajectories, while the neXtSIM results are much more consistent with the results from 15 the buoy trajectories. The information on the mean drift and diffusivity fields provided by our analysis can be used in an advection/diffusion equation or with Lagrangian passive tracers models to study the drift of e.g. pollutants or micro-organisms moving with the ice. More generally, the analysis presented in this paper should be seen as a useful evaluation metric of coupled sea ice-ocean models that aim at being used in operational forecasting platforms, for process and climate studies. 20 1 The Cryosphere Discuss., doi:10.5194/tc-2015-233, 2016 Manuscript under review for journal The Cryosphere Published: 26 January 2016 c Author(s) 2016. CC-BY 3.0 License.


Introduction
The main goal of the present study is to evaluate the sea ice velocity fields provided by two sea ice models that are currently used to produce operational forecast and to simulate pollutant spread in ice-covered areas.A secondary objective is to provide accurate information on specific statistical properties of sea ice motion that could be used to better assess the predictability of the ice trajectories, 25 hence improving the evaluation of their uncertainty.
In order to state the problem, it is interesting to get back to the review of Drozdowski et al. (2011) on oil spill trajectory modelling in the presence of sea ice, in which they made the following statement: "The comprehensive oil spill trajectory modelling system can be expected to provide reliable short term predictions, perhaps for several days to a week.However a spill that occurs in the fall may 30 be inaccessible for 6 or more months.The picture for oil spill trajectory modelling for the long term, say 6 months, is less clear.The multiple pathways for oil movement and the uncertainty in the ocean currents, sea ice drift and the oil-ice interactions means that deterministic model predictions (i.e. a definite prediction of where the oil will go) become unreliable, and that probabilistic predictions (i.e. the model assigns probabilities that the oil will be found at any particular location) are required.The 35 key point is that any ice-ocean model used for the oil spill trajectory modelling needs to provide a credible representation of the large scale spatial variability in the ocean currents and the sea ice drift in addition to the storm driven variability." An interesting research task is to investigate what is the regime of diffusion of a pollutant trapped in sea ice.To do so, one can separate the deterministic from the non-deterministic part of the sea 40 ice motion within the theoretical framework developed for turbulent diffusion in fluids.Indeed, previous studies showed that turbulent theory may be applied to study diffusion properties of sea ice (Lukovich et al., 2011(Lukovich et al., , 2015;;Rampal et al., 2009b;Thorndike, 1986;Colony and Thorndike, 1985).
Following the turbulent theory of Taylor (1921), Thorndike (1986) were able to split sea ice motion into a predictable mean part, ū, and a random fluctuating part, u 0 .The predictable part is forced by 45 the large scale mean motions in the atmosphere and the ocean, while the random part was suggested by Thorndike (1986) to be forced by short term fluctuations in the wind and ocean currents.Later, Rampal et al. (2009b) showed that the fluctuating part of the sea ice motion is also influenced by the mechanical response of the sea ice cover itself, and presented a method based on the turbulent theory of Taylor (1921) to extract this fluctuating part from the total motion.Note that similar methodology 50 has been applied to study diffusion properties from Lagrangian drifters in the ocean (see e.g., Zhang et al., 2001;Poulain and Niiler, 1989).With this method it is possible to study the characteristics of the sea ice motion in a statistical sense, for different conditions and seasons.
In this study we analyse the quality of the sea ice motion fields provided by the TOPAZ-CICE (hereafter referred as TOPAZ) and neXtSIM sea ice models by comparing synthetic buoy trajectories January 2015).For the period 2007 to 2010 we also generate "virtual" buoy trajectories by putting floats in the TOPAZ and neXtSIM models.The floats are initialised at the same time and position as the IABP buoys, and are "killed" when the IABP buoy track stops or when the virtual buoy enters 75 into an area of simulated open water (sea ice concentrations less than 15%).After removing the few non-overlapping periods (i.e. the periods for which the trajectories are not available in the three data sets), three comparable data sets are obtained: i) the observed sea ice trajectories, ii) the trajectories of virtual sea ice floats simulated by TOPAZ, and iii) the trajectories of virtual sea ice floats simulated by neXtSIM.Because the buoy positions of the IABP dataset are sampled every 12 hours, 80 the virtual TOPAZ and neXtSIM float trajectories originally sampled every hour are sub-sampled to obtain 12-hourly positions as well.
We increased the number of buoy trajectories by splitting the buoy trajectories longer than 35 days into as many segments of 35 days as possible, the remaining data being discarded.Since we know that the memory of a piece of ice of its past motion is less than 10 days (Rampal et al., 2009b), 85 each of these 35-days buoy track segments can be assumed to be independent from the others.When removing non-overlapping segments in the three data sets, we are left with 280 individual floats and 3720 35-day segments in each data set for the period 2007 to 2010.
A polar stereographic projection is used to change the IABP and Virtual buoy positions from geographic (lat, lon) coordinates to Cartesian (x, y) coordinates in km from the north pole.Velocities 90 are then calculated along each buoy track as at positions x = (x(t + t) + x(t))/2, ỹ = (y(t + t) + y(t))/2 and time t = t + t/2.The time step, t is 12h for the IABP data and for the virtual buoys.

IABP
Originally, the raw IABP buoy positions were sampled irregularly in time with a mean time interval of 1 hour, and with errors ranging from 100 to 300 m depending on the position system (Thomas, 1999).Before being delivered to the scientific community, the buoy positions were interpolated (using a cubic function) to form an homogeneous trajectory data set giving for each buoy its position 100 every 12 hours.
To investigate sea ice motion properties in the pack ice, we restricted the IABP dataset to a region located in the centre of the Arctic basin (hereafter denoted the Central Arctic domain), i.e. above 70 N, at least at 100 km from the nearest coast and from the section of the 80 th parallel going from Greenland to Severnaya Zemlya.Sea ice dynamics in coastal regions are specific with for 105 example the presence of land-fast ice and would require a dedicated study.We also restrict the analysis to the winter period defined as starting the 1 st of November and ending the 15 th of May.
Each individual buoy track from the IABP dataset has been checked manually to clean them from unrealistic "jumps" or "spikes" in the trajectories, or when the dating system was giving obviously wrong times.These unrealistic "jumps" are due to errors in the positioning system embarked on 110 the buoy or in the recordings during the deployment or recovering of the instruments.The period 2007-2010 has been selected for its relatively good coverage, with more than 40 buoys recording their positions simultaneously every day.Buoy trajectories shorter than 30 days have been removed.2.2 Float tracks generated with the TOPAZ model 120 The forecast system TOPAZ is the official Arctic Ocean forecast platform of the European Copernicus Marine Environment Monitoring service (http://marine.copernicus.eu).The ocean part of TOPAZ uses HYCOM version 2.2, with 28 vertical layers divided into isopycnal layers in the stratified inte- rior of the open ocean and z-coordinates in the unstratified surface mixed layer (Bleck, 2002).The ocean model is coupled to a one thickness category sea ice model based on CICE, the Los Alamos 125 Sea Ice Model, version 4 (Hunke and Lipscomb, 2010).The part that solves the sea ice dynamics is built around a standard EVP rheology (Hunke and Dukowicz, 1997) while the one solving sea ice thermodynamics is described in Drange and Simonsen (1996).The sea ice strength is set to 27 500 N and the advection of sea ice scalar variables is calculated using a 3 rd order WENO scheme (Jiang and Shu, 1995) with a 2 nd order Runge-Kutta time discretisation.The sea ice-ocean model is run 130 here in free-run mode (i.e., no data assimilation is applied).The model is initialised from a restart file taken from the free-run simulation described in Sakov et al. (2012).The TOPAZ coupled sea ice-ocean system provides sea ice variables, such as sea ice velocities, concentration, and thickness, at a high temporal resolution (i.e.every hour).A more detailed description of the TOPAZ system may be found in Sakov et al. (2012).The float tracking with TOPAZ is performed off-line.Hourly sea ice motion simulated by the model is used to force the off-line float tracking system that moves the floats in the quasi-homogeneous TOPAZ Arctic grid.We use that grid instead of a regular longitude/latitude grid in order to avoid singularity errors at and around the Pole.The advantage of the off-line float tracking system is that it allows us to easily perform experiments for a given time period and location and with a reasonably 145 large number of floats.We checked that for the time scale and spatial resolution considered here, this off-line tracking method is giving similar results to using an online tracking system, while remaining computationally efficient.
The float-tracking system initialises floats at a prescribed start position (longitude and latitude) and starting time (year, month, day and hour), then moves the floats with a simple Eulerian method 150 and kills the floats at a prescribed end time (year, month, day and hour) so that it matches the end of the corresponding real IABP buoy's trajectories.The Eulerian sea ice velocities given by the TOPAZ model are interpolated to the position of the virtual Lagrangian floats every hour.Floats drifting out of the simulated sea ice cover, i.e.where sea ice concentration is below 15% in the present case, are not tracked further.The virtual buoy tracks are stored as longitude and latitude positions.

Float tracks generated with the neXtSIM model
neXtSIM is a fully-Lagrangian thermodynamic-dynamic sea ice model, using an adaptive finite element mesh and a mechanical framework based on the elasto-brittle rheology (Rampal et al., 2015).The neXtSIM model is able to simulate the observed evolution of the sea ice volume, extent and area over an annual time scale (Rampal et al., 2015).However, the simulated drift and deformation fields have only been extensively evaluated for the winter season in Rampal et al. (2015) and Bouil- AMSR-E passive microwave sensor (Kaleschke et al., 2001;Spreen et al., 2008, data obtained from the Integrated Climate Date Center, University of Hamburg, Germany, http://icdc.zmaw.de)and the 175 TOPAZ ice thickness, within the area reported by AMSR-E as being covered with ice.The modelled ice thickness of the TOPAZ model is known to be biased (too low) and so we increase the initial thickness uniformly so that the total volume is the same as that given by the PIOMAS model on September 15 th 2007, 2008 and 2009 (Zhang and Rothrock, 2003).The good performance of PIOMAS in simulating Arctic sea ice volume as compared to available observations is reported in 180 Schweiger et al. (2011).The temperature and salinity of the slab ocean model are initialised with temperature and salinity from TOPAZ.
The model is forced with the ocean state of the TOPAZ4 reanalysis (see Sakov et al., 2012).The oceanic forcing variables are sea surface height, velocity at 30 m depth, and sea surface temperature and salinity (for more details, see Rampal et al. (2015)).The atmospheric state comes from the Arctic Interim is a high resolution atmospheric reanalysis (30 km) known to reproduce particularly well the near-surface wind fields in the Arctic region (Bromwich et al., 2015).
The frictional drag parameter for the atmosphere-ice boundary layer is closely linked to the ap-190 plied atmospheric forcing field and directly impacts sea ice motion.The common approach to optimise this parameter is to compare the simulated and observed sea ice drift over the whole Arctic basin.However, this approach leads to interdependence between the optimisation of the mechanical parameters and the optimisation of the air drag coefficient.Instead, we developed a tuning methodology that has the advantage of differentiating the choice of the air drag coefficient from the mechanical 195 parameters.This method and its results are presented in Rampal et al. (2015).We here use the values for the air drag parameter (c a = 0.0076) found to be optimal when using the ASR wind forcing and the same value of water drag parameter as in TOPAZ (c w = 0.0055) The float tracking with neXtSIM is performed at run time.The main reason for doing this is that the Lagrangian advection used in the neXtSIM model offers some additional challenges to a 200 post-processing approach over the Eulerian advection used in TOPAZ.Using this approach, the tracking system is initialised with floats that have the same positions and starting dates as the IABP buoys.These positions then change as the underlying mesh moves and are interpolated onto the new mesh after each remeshing step.(See Rampal et al. (2015) for more details on the remeshing procedure.)The virtual drifters are removed when the corresponding IABP buoy tracks stop or when 205 the simulated local ice concentration falls below 15%.The procedure of adding and removing buoys is therefore the same as the one followed for the aforementioned off-line float tracking system used with TOPAZ.

Methods
The method used to analyse the Lagrangian particles/floats obtained from IABP, TOPAZ and neXtSIM 210 is similar to the one presented in Rampal et al. (2009b) and is based on the theory proposed by Taylor (1921) to study turbulent fluid.For the decomposition of the motion into a mean and a fluctuating part, we follow the classical approach used to study Lagrangian trajectories (see for example Zhang et al., 2001).Taylor (1921) 215 Following the theory developed for turbulent fluids by Taylor (1921), diffusion describes in a statistical sense how an individual particle moves apart from any of its previous position.Taylor's diffusion theory is valid for steady and homogeneous turbulent flow without mean flow and for which the fluctuating velocity follows a Gaussian distribution.When following a single particle in such conditions, the variance of its fluctuating displacement hr 02 (t)i should in theory evolves as

Diffusion theory of
where hu 02 i, the variance of the fluctuating velocity, is constant in time and C(⌧ ), the Lagrangian normalized autocorrelation function, is defined as where T max is the duration of the particle trajectory (here 35 days) and u 0 (t) is its fluctuating velocity 225 at time t defined following the method described in Section 3.2.
For very long time intervals ⌧ , the autocorrelation vanishes and the integral of is then a constant.This constant is referred as the Lagrangian integral time scale and determines the transition between two diffusion regimes.For time much shorter than , we are in the "ballistic" 230 regime and Equation 3 becomes (this simply comes from the fact that C(⌧ ) tends to the limiting value unity for small t).For time much longer than , we are in the "Brownian" regime (also called "random walk" regime) and Equation 3becomes where ↵ is a constant defined as , 2008).The second regime is similar to the one driven by molecular diffusion (i.e., where fluctuating velocities are also uncorrelated) but with much larger diffusion coefficients.
Following Lagrangian turbulent theory, diffusivity (noted K) is defined as follows: In the "ballistic" regime (with Equation 6), diffusivity increases with time and may be calculated as In the "Brownian" regime (with Equation 7), diffusivity (also called absolute diffusivity in that case) is constant and may be calculated as In the present study we will use the absolute diffusivity to compare the different data sets.

Decomposition of the sea ice motion
As shown in the previous section, one may look at statistical properties of fluctuating displacements to study diffusion properties.Therefore, we need to decompose the total sea ice motion into a mean 250 and a fluctuating part.We here shortly describe how this decomposition is performed and how the fluctuating part of the motion is then described.A more complete and detailed presentation of the methodology, in particular of how to choose the appropriate temporal and spatial averaging scales, can be found in Rampal et al. (2009b).

Mean velocity calculation 255
From the list of positions x i q of a buoy q given with a regular time interval t, one can evaluate its position and velocity at time t = t i q + 0.5 t by computing: By doing the same for all the available buoys, one can build a data set recording all the velocities 260 computed at a time scale t.
From that data set a mean velocity field ūL,T (x,t) can be defined for a given spatial and temporal scale, L and T , by averaging all the buoy velocities u k recorded at a distance less than L/2 from x and within the time window [t T/2; t + T/2].The averaging operator used here is where w k are the weight coefficients defined as The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -233, 2016 Manuscript under review for journal The Cryosphere Published: 26 January 2016 c Author(s) 2016.CC-BY 3.0 License.
This mean velocity is then evaluated at each recorded position x q (t) and subtracted from the recorded velocities to define the fluctuating velocities by computing 270 As seen in Equation 15, the fluctuating velocities depend on the value of the averaging scales, L and T .This dependance and the fact that the autocorrelation function has to remain close to 0 for long time interval ⌧ provide an elegant method to determine the appropriate value for L and T .First we compute for each buoy trajectory q and pair of averaging scales L and T , the normalized autocorrelation function, C L,T q , using Equation 4. Secondly, we calculate an ensemble average 275 autocorrelation function as where the influence of each buoy is weighted by the time length of its trajectory, T q max .Finally, we select the averaging scales for which L,T (⌧ ) remains close to 0 for large ⌧ .An example of L,T (⌧ ) function obtained with L = 400 km and T = 165 days is presented in Figure 2. The results of this 280 analysis are detailed and discussed in Section 4.
An equivalent method to the one described above is to define the appropriate spatial and temporal averaging scales as the lowest scales for which L,T remains quasi constant (i.e., less than 1% change).Since we cannot integrate Equation 5 to infinity, the average integral time scale L,T for the selected period is computed as where t 0 is the first time L,T (⌧ ) crosses zero (see for instance Poulain and Niiler, 1989;Rampal et al., 2009b).In the example of Figure 2, t 0 = 6 days and the integration time scale for L = 400 km and T = 165 days is equal to 1.5 days.Note that the absolute diffusivity (K = 1.0 ⇥ 10 3 m 2 s 1 ) for the example of Figure 2 is also given for information but will be discussed in Section 4. We verified 290 by plotting L,T (not shown) that a plateau is reached for averaging scales larger than L = 400 km and T = 165 days, which is then the appropriate averaging scale for the example shown here.

Results
The result of the diffusion analysis applied to the whole IABP data set from 1979 to 2011 is used as a reference to get an overall picture for the Central Arctic region.In addition, we present a comparison 295 between observed and simulated float trajectories for the winter seasons 2007 to 2010.We restrict our analysis here to winter seasons because of the too sparse buoy trajectories available in the IABP dataset during summer, and because the performance of the neXtSIM model has been extensively assessed so far for winter only (Rampal et al., 2015).
Figure 3 shows the maps with all the buoy trajectories for each winter season from the IABP 300 data set, the TOPAZ model and the neXtSIM model.From this figure, we can already see that the trajectories from the neXtSIM model are more similar to observations than the ones from the TOPAZ model.We clearly see that, in general, the TOPAZ trajectories are longer than their corresponding ones in the IABP data set.The short IABP trajectories north of the Canadian Arctic Archipelago (clearly seen in 2007/2008 winter) are the result of the significantly thicker and more ridged sea ice, 305 which leads to nearly zero speed in that region.This dynamical sea ice response, i.e. the expected drift behaviour of a solid thick plate of ice surrounded by closed boundaries and stressed by external forces, is captured by neXtSIM but not by TOPAZ.To better describe the differences between the simulated and observed trajectories, we analyse separately the mean drift and the fluctuating part of the motion by applying the so-called Taylor's decomposition to the sea ice velocities.

Decomposition of the sea ice motion
To compute the mean part of the motion, Rampal et al. (2009b) determined that the appropriate averaging scales to be used for the entire IABP data set are L = 400 km and T = 5.5 months (165 days) for winter, and L = 200 km and T = 2.5 months for summer.These spatial and temporal averaging scales are defined as the smallest values for which the integral time L,T remains quasi constant.By 315 performing a similar analysis as in Rampal et al. (2009b) for the IABP data but only from 2007 to 2010, we found that the same averaging scales should be used and so we did to compute the mean part of the motion at each location of buoy along their track in the present study.
Applying the averaging scales L = 400 km and T = 165 days with the procedure described in section 3.2.1,we split the velocity field into a mean, ū, and a fluctuating part, u 0 .The fluctuating sea 320 ice displacement is then derived from the fluctuating velocities as with r 0 x (t = 0) = 0 and r 0 y (t = 0) = 0 and where t is the time step in the data set.Figure 4 shows an example of the partition of the displacement of an IABP buoy into its mean and fluctuation trajectories.With this partitioning method, the mean motion can be considered ho- The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -233, 2016 Manuscript under review for journal The Cryosphere Published: 26 January 2016 c Author(s) 2016.CC-BY 3.0 License.fluctuating velocities should be symmetrically distributed around zero.This is the case in our results (not shown), meaning that one can directly look at the speed (i.e., norm of the velocity) without losing information.The PDFs of the fluctuating speeds are plotted in Figure 7 with the Gaussian and exponential fits indicated for reference.We clearly see the fluctuating speeds of the IABP buoy 365 data follow an exponential distribution with a mean equal to 6.9 cm s 1 .It is important to note here that the data follow an exponential distribution instead of a Gaussian distribution, as expected in fully developed turbulence (Batchelor, 1960;Frisch, 1995) and observed in different turbulent fluids (Van Atta and Chen, 1970;Zhang et al., 2001).This means that the sea ice fluctuating speed can be much larger than a standard deviation away from the (zero) mean.However, such non-Gaussian 370 distributions for fluctuating speeds are also observed for oceanic surface currents during energetic events associated with large organised structures such as jets and vortices.Such a signature for sea ice may indicate that sea ice dynamics are dominated by the passage of large perturbations over the Arctic, whereas less energetic features have less impact on sea ice motion.This selective sensitivity to energetic events may be related to the intrinsic properties of solids associated with threshold 375 mechanics.This seems to be supported by the fact that for weaker seasonal sea ice, the observed fluctuating velocities rather follow Gaussian statistics (Lukovich et al., 2011).
The fluctuating speeds from TOPAZ are too high on average (by about 30%) with a mean value equal to 8.97 cm s 1 , and their statistics do not follow an exponential distribution, hence missing the highest values of observed fluctuating speed.The fluctuating speeds from neXtSIM are slightly too 380 low (by about 10%) with a mean value equal to 6.14 cm s 1 and follow an exponential distribution within the range 0 to 30 cm s 1 .The neXtSIM model also misses the highest values of observed sea ice fluctuating speed.The differences between the performances of the two models might come from the different rheologies, but also from the differences in the initial conditions and external forcings that are used.For example the too thin ice pack of TOPAZ at the onset of the freezing season may 385 contribute to the overestimation of the winter drift (Rampal et al., 2011;Ólason and Notz, 2014) and could be more prone to have fluctuating speeds following a Gaussian distribution as discussed by Lukovich et al. (2011).The forcings are also not the same since the atmospheric reanalysis ASRinterim, used here to force neXtSIM, has higher spatial and temporal resolution than ERA-interim, used here to force TOPAZ.Another difference to be notified is the fact that the results of TOPAZ 390 presented here come from a free run with no data assimilation whereas the results of neXtSIM are produced with the ocean forcing of the TOPAZ reanalysis including data assimilation.To find and explain the origins of the models performances is out of scope of the present paper and would require a dedicated analysis.However, as the forcing fields and the rheology used in TOPAZ are widely used in the sea ice modelling community, we expect the performances shown here by the TOPAZ system 395 as a reasonable reference of most of the available coupled ocean-sea ice forecasting platforms.
Using the fluctuating part, u 0 of the velocity we calculate the fluctuating displacement components r 0 x and r 0 y , see Equation 18, and the norm r 0 of the fluctuating displacement.By analysing By doing so, we make sure that the variance of the fluctuating velocities hu 02 (t)i, where t here goes from 0 to 35 days, is almost constant.We verified that for each data set the deviation of hu 02 (t)i to the mean values hu 02 i is maximum of about 10%. 410 The standard deviation of the fluctuating displacement may be a crucial piece of information for the planning of a recovery operation in a case of drifting oil or pollutant that is trapped or attached to the ice, as it gives an estimate of how the size of the searching area around the predicted mean drift should increase through time, in a statistical sense.If an operator can only trust the mean drift, which is the case for forecast longer than a few days, and if the size of the polluted ice is so small that it can 415 be considered as a single particle (i.e. one can assume no relative dispersion of the pollutant), the searching area could be defined as a circular region with a radius depending on the standard deviation of the fluctuating displacement.We verified that about 68.9%, 95.9% and 99.6% of the fluctuating displacements are smaller than 1, 2 and 3 standard deviations, which means that the fluctuating displacement distribution is in the Gaussian attraction basin.Another way of representing the same 420 information is to draw the relative displacement as if seen from the mean drift (see bottom panel of Figure 8).In this example, the searching radius (defined here as equal to 3 standard deviations) should be about 87 km after 5 days (corresponding to a surface area of 24000 km 2 ), and about 206 km after 30 days (corresponding to a surface area of 134000 km 2 ).As the mean speed and deformation in the Central Arctic are increasing (Rampal et al., 2009a), we may expect higher values for recent 425 and coming years.
Another way of analysing the results is to look at the variance of the fluctuating displacement, hr 02 (t)i, as a function of time (see Figure 9).This variance is sometimes called "dispersion" or "spread" around the trajectory given by the mean drift.In this paper, we rather use the term "single particle dispersion" or "absolute dispersion" to mark the difference with the relative dispersion, 430 which indicates how the distance between two particles evolves through time.In these log-log plots the two diffusion regimes are indicated by the dashed lines showing the slope of the initial "ballistic" regime where the displacement grows with t 2 , hr 02 (t)i ⇠ t 2 , and the later "Brownian" regime where the displacement grows with t, hr 02 (t)i ⇠ t, as in the Taylor (1921) turbulent diffusion theory.The time scale at which the regime transition occurs corresponds to the integral time scale, which is 435 observations and those obtained from TOPAZ and neXtSIM, and found 0.7 and 0.85, respectively.
One should note that the diffusivity fields could, for example, be used to evaluate the fluctuating displacement variance for t >> by using the relationship hr 02 (t)i = 2Kt, which in turn could be used to estimate the searching radius in case of pollutant spill in ice-covered waters as r s = 3 p hr 02 (t)i.
55 simulated by the models to the corresponding buoy trajectories of the International Arctic buoy Program (IABP) data set.The statistical properties of these trajectories are analysed to evaluate how ing Buoy Pressure, Temperature, Position, and Interpolated Ice Velocity, Version 1. subset C. Boulder, Colorado USA.NSIDC: National Snow and Ice Data Center.http://dx.doi.org/10.7265/N53X84K7.

Figure 1
Figure 1 shows all the trajectories analysed in this study.The IABP data set from 1979 to 2011 covers almost the whole Central Arctic domain except the East Siberian and Laptev Seas.The dif-115

135
The applied atmospheric forcing fields are the 6-hourly 10-meter wind velocities from the ERA interim reanalysis (ERAi) distributed at 80 km spatial resolution (http://www.ecmwf.int/en/research/The Cryosphere Discuss., doi:10.5194/tc-2015-233,2016 Manuscript under review for journal The Cryosphere Published: 26 January 2016 c Author(s) 2016.CC-BY 3.0 License.climate-reanalysis/era-interim, ECMWF (2011)).The frictional drag parameters for the atmosphereice stress (c a = 0.0016) and for the ice-ocean stress (c w = 0.0055) are those currently used and optimised for the TOPAZ operational platform.140 155 The domain used for the simulations with the neXtSIM model is defined from TOPAZ (version 4) coastlines and open boundaries.It covers the Arctic and North-Atlantic Oceans, extending from an 160 open boundary at 43 N in the North-Atlantic to an open boundary in the Bering Strait.The mean resolution of the finite element mesh used by neXtSIM is about 10 km.Thermodynamic growth and melt of the ice are based on Semtner (1976) and the ice model is coupled to a slab ocean model.The temperature and salinity in the slab ocean model are constantly nudged towards values of the TOPAZ reanalysis in order to simulate oceanic heat and salt transfer.A complete description of this 165stand-alone version of the neXtSIM sea ice model as well as his coupling to the slab ocean model may be found inRampal et al. (2015).
170lon andRampal (2015).The three simulations used for this study are therefore starting on September 15 th and finishing on May 15 th 2007, 2008 and 2009, and our analysis is restricted to winter time, i.e. from November to May.The model is initialised with the ice concentration derived from the The Cryosphere Discuss., doi:10.5194/tc-2015-233,2016 Manuscript under review for journal The Cryosphere Published: 26 January 2016 c Author(s) 2016.CC-BY 3.0 License.

The
Cryosphere Discuss., doi:10.5194/tc-2015-233,2016   Manuscript under review for journal The Cryosphere Published: 26 January 2016 c Author(s) 2016.CC-BY 3.0 License. the characteristics of the norm we quantify how the distance between any particular trajectory and the mean trajectory evolves through time.To increase the robustness and statistical significance of 400 our analysis, we artificially increase the number of buoy trajectories by splitting each trajectory into 35-day segments starting every 12 hours, i.e. every time a new buoy position along-track is available.

405Figure 8
Figure 8 shows the evolution of the norm of the fluctuating displacement for every tenth segments retrieved from the IABP trajectories for the winter periods from 1979 to 2011.The dashed red line indicates the distance corresponding to 3 standard deviations of the fluctuating displacement.The standard deviation of the fluctuating displacement is the square root of hr 02 (t)i, which is defined as