Arctic sea-ice diffusion from observed and simulated Lagrangian trajectories

We characterize sea-ice drift by applying a Lagrangian diffusion analysis to buoy trajectories from the International Arctic Buoy Programme (IABP) dataset and from two different models: the standalone Lagrangian sea-ice model neXtSIM and the Eulerian coupled ice–ocean model used for the TOPAZ reanalysis. By applying the diffusion analysis to the IABP buoy trajectories over the period 1979– 2011, we confirm that sea-ice diffusion follows two distinct regimes (ballistic and Brownian) and we provide accurate values for the diffusivity and integral timescale that could be used in Eulerian or Lagrangian passive tracers models to simulate the transport and diffusion of particles moving with the ice. We discuss how these values are linked to the evolution of the fluctuating displacements variance and how this information could be used to define the size of the search area around the position predicted by the mean drift. By comparing observed and simulated sea-ice trajectories for three consecutive winter seasons (2007–2011), we show how the characteristics of the simulated motion may differ from or agree well with observations. This comparison illustrates the usefulness of first applying a diffusion analysis to evaluate the output of modeling systems that include a sea-ice model before using these in, e.g., oil spill trajectory models or, more generally, to simulate the transport of passive tracers in sea ice.


Introduction
Sea-ice motion can be viewed as a superposition of a mean circulation and turbulent-like fluctuations (Rampal et al., 2009b).Using such decomposition for studying pollutant transport by sea ice was already proposed by Colony and Thorndike (1985), who analyzed sea-ice drift data covering the period 1893-1984 while using arbitrary averaging scales (90 years and 1500 km) to define the mean motion.By using a denser sea-ice drift dataset covering the period 1978-2001 and the theoretical framework introduced by Taylor (1921) for the analysis of turbulent fluids, Rampal et al. (2009b) proposed a methodology to rigorously decompose sea-ice motion into mean and turbulent-like fluctuating parts.The appropriate averaging scales (about 400 km and 5.5 months for winter conditions) were found small enough to clearly separate the interannual variability of the mean circulation from the fluctuating motion due to passing atmospheric perturbations, local oceanic eddies and inertial and tidal motion.This approach, based on the analysis of single particle trajectories, has been widely used to study diffusion properties of Lagrangian drifters in the ocean (see, e.g., Zhang et al., 2001;Poulain and Niiler, 1989) and is now becoming a standard analysis tool for sea-ice dynamics (Lukovich et al., 2011(Lukovich et al., , 2015;;Gabrielski et al., 2015).
Single particle analysis (here referred to as diffusion analysis) is particularly useful for characterizing long-term trajectories as it clearly decomposes the motion into mean (predictable) and fluctuating (unpredictable) parts (Colony and Thorndike, 1985).It also allowed Rampal et al. (2009b) to show that sea-ice diffusion exhibits a clear transition from the so-called ballistic regime to the Brownian regime.This transition is also typical of turbulent fluids and is due to the fast decay of the velocity autocorrelation function.The information coming from the diffusion analysis (mean flow and diffusivity) statistically describe the ensemble of all the potential trajectories that a particle, released at a given location Published by Copernicus Publications on behalf of the European Geosciences Union.P. Rampal et al.: Diffusion in the ice pack: comparison between models and observation at an unknown time during a season, could follow and may then be sufficient to produce a probabilistic forecast of tracer transport.
To simulate tracer transport, one can use continuous or discrete passive tracer models.Continuous models are usually based on the following advection-diffusion equation: where ū is the mean velocity field, K is the corresponding diffusivity and C may describe either the tracer concentration or the probability to find the tracer in a given position after a given time (positional probability) as explained by La-Casce (2008).One can also use discrete passive tracer models for which the displacement dx i in the i direction is defined for independent objects.The simplest approach, known as the "random walk" model or zeroth-order model, is strictly equivalent to the advection-diffusion equation of continuous models and simply defines dx i as where u i is the fluctuating velocity in the i direction and dw i is a Wiener process.The first-order approach is also often used as it can represent the transition between the ballistic and the Brownian diffusion regimes by applying the stochastic term on the evolution of the velocity, leading to the following set of equations: where i is the integral timescale in the i direction.Firstorder approaches can also successfully reproduce the loops often seen in surface drifter trajectories by adding a rotation term similar to the Coriolis term.These approaches have been widely used to study the spread of pollutant and other tracers by eddy turbulence in the ocean and atmosphere (see LaCasce, 2008, for an extensive review).Using such approaches to simulate the spread of sea ice would be inappropriate because sea-ice motion fields are intermittent in time and discontinuous in space, but are valid for reproducing the statistics of individual sea-ice trajectories.The "random walk" model was, for example, used in Colony and Thorndike (1985) with a mean field and statistics on the fluctuations derived by Colony and Thorndike (1984).
Another way to use passive tracer models is to replace the mean field ū in Eqs.(1), ( 2) or (3) by simulated or observed Eulerian velocity fields.This approach is widely used with passive tracer models directly forced by motion fields simulated by an hydrodynamical model (e.g., Nudds et al., 2013), given by a reanalysis (e.g., Gearon et al., 2014) or derived from satellite observations (e.g., Fowler et al., 2004).The diffusion term (or in the discrete models, the Wiener process term) is either neglected or defined such that it accounts for the unresolved part of the fluctuating motion.The unresolved part of the motion could be analyzed with the methodology proposed by Dominicis et al. (2012) for ocean surface tracer modeling, which consists of comparing the characteristics of the fluctuating part of observed trajectories to the ones of trajectories given by a tracer model forced by model output.
Before using one of these approaches one needs to answer a few questions: What are the right averaging scales to define the mean motion field?What are the statistical properties of the fluctuating part of the motion?Is there a transition between different regimes of diffusion?Are the mean and fluctuating parts of the motion correctly reproduced by the forcing field?When not, could this indicate that some processes are missing in the forcing velocity fields?If the fluctuating part is not well reproduced, could it be compensated by adding extra terms in the tracer equation?These important questions are not always answered before running tracer models forced by sea-ice velocity fields and this could strongly impact the validity of the studies based on such results.This is of particular importance in the context of the increasing activities in the Arctic seas (e.g., shipping, fishing, and oil and gas exploration and exploitation), which enhance the risk of pollution in a region where ecosystems are already under threat from the amplified effects of climate change.The extreme conditions (e.g., presence of ice, extreme cold, high winds and the polar night) and the long distance from well-equipped facilities may hamper access to the polluted area for several months (Drozdowski et al., 2011).These conditions also make the detection of the pollution challenging and slow or even stop the natural and artificial degradation processes.In addition, the pollutants are often trapped in or under the sea ice and, therefore, may be transported by the ice over large distances before being released (Rigor and Colony, 1997).To improve the understanding of sea-ice trajectories is therefore crucial for risk assessment and response planning related to pollutant release in Arctic seas.Passive tracer modeling in sea ice also has other applications, for example studying biology and its link with pollutants (e.g., Borgå et al., 2002;Pfirman et al., 1995) or estimating the age of the Arctic sea-ice cover (e.g., Fowler et al., 2004;Hunke, 2014).
In this paper, we demonstrate the utility of applying Lagrangian diffusion analysis to sea-ice trajectories in the context of passive tracer modeling.The analyses presented in this paper are restricted to winter conditions, as this season has been identified as more critical for oil spill recovery operations (Drozdowski et al., 2011).In Sect.2, we apply the same method as in Rampal et al. (2009b) to the International Arctic Buoy Programme (IABP) dataset for the period 1979-2011 (Rigor, 2015).This reference dataset is analyzed to get an overall picture of the characteristic of the mean and fluctuating parts of the motion for the central Arctic domain and to derive the quantities (diffusivity, Lagrangian integral timescale, etc.) needed for tracer models.This section also includes a discussion on how to predict the evolution of the fluctuating displacements variance, which could be used to define the size of the search area around the trajectory predicted by the mean drift.In Sect.3, we follow the methodology proposed by Dominicis et al. (2012), which consists of applying the diffusion analysis to observed and simulated trajectories to evaluate the merits of using these simulated fields for tracer modeling.This evaluation exercise is done for simulations obtained from two different models.Section 4 sums up the main conclusions.

Diffusion analysis on a reference dataset
In this section we present the theoretical framework of the diffusion analysis by showing its application to the IABP dataset from 1979 to 2011 as a reference.Results and theory are then compared and discussed in the context of passive tracer modeling.This is the same analysis as that already done by Rampal et al. (2009b), except that we use the 12hourly buoy positions data instead of the hourly data and extend the analysis period up to 2011.The results we obtained are similar to those reported in Rampal et al. (2009b), but additional evaluation and discussion are presented in Sects.2.4 and 2.5.

Reference dataset
Figure 2 shows all the trajectories analyzed in this section.
We restrict the analysis to the winter period defined as starting on 1 November and ending 15 May.Results for summer can be found in Rampal et al. (2009b).To investigate seaice motion properties in the pack ice, we restrict the IABP dataset to a region located in the center of the Arctic basin (hereafter denoted the central Arctic domain, blue line on Fig. 2).This includes all buoy data north of 70 • N, except between 20 • W and 100 • E where the southern limit is set to 80 • N. In addition only buoy data from more than 100 km off the coast are used.The selected data cover the whole central Arctic basin, but the data coverage in the East Siberian and Laptev seas is sparse (see the maps showing the number of buoys and records in Fig. 2).Sea-ice dynamics in coastal regions are specific with, for example, the presence of land-fast ice and would require a dedicated study.The IABP buoys are mainly deployed over multi-year ice and thus the conclusions from the analyses presented in this study may not be extrapolated for weaker seasonal ice.
The raw IABP buoy positions are sampled irregularly in time with a mean time interval of 1 h, and with errors ranging from 100 to 300 m, depending on the positioning system they used (Thomas, 1999).Before being published, however, the buoy positions are interpolated in time (using a cubic function) to form an homogeneous trajectory dataset, giving for each buoy its position every 12 h.
We manually checked each individual buoy track from the IABP dataset to clean them from unrealistic "jumps" or "spikes" in the trajectories and from obvious errors of the dating system.The unrealistic "jumps" present in buoy trajectories are due either to errors in the positioning system installed in the buoy or to wrong recordings during the deployment or recovery phase of the instruments.A polar stereographic projection is used to change the IABP and virtual buoy positions from geographic to Cartesian (x, y) coordinates.
P. Rampal et al.: Diffusion in the ice pack: comparison between models and observation

Decomposition of the sea-ice motion
The IABP buoy trajectories are geometrically complex, with abrupt changes in direction (Fig. 2).
Therefore, it is helpful to decompose the total sea-ice motion into a mean part, which should be considered homogeneous and stationary, and a fluctuating part, which should contain the unpredictable, local or non-stationary motion.This is done here by following the classical approach used to study Lagrangian particle trajectories (see for example Zhang et al., 2001).This consists of splitting each trajectory into mean and fluctuating parts by using appropriate averaging scales L and T to compute the Lagrangian mean motion at any given location and time.An example for one particular IABP buoy trajectory is shown on Fig. 1.
From the list of positions x i q of a buoy q, one can evaluate its position and velocity at time t i q = t i+1 q + t i q /2 by computing By repeating this for all the available buoys, we build a dataset of 12-hourly velocities and positions, from which we compute the mean and fluctuating part of the motion.
The mean velocity field ūL,T (x, t) is defined for any target position x and time t as ūL,T (x, t) = 1 q,i w i q q,i w i q u i q , where L and T are the spatial and temporal averaging scales and w i q are weight coefficients.In this study we use constant weights, meaning that the mean velocity is simply defined as an average of all the 12hourly velocities available in the dataset that are within a circle of diameter L centered on the target position x and a time window of duration T centered on the target time t.The fluctuating velocities u are then computed at each position x i q and time t i q by subtracting the mean velocity ūL,T ( x i q , t i q ) from the total velocity u i q .Both the mean and fluctuating velocities are then defined for a specific pair of averaging scales L and T .The mean and fluctuating displacements for each buoy (e.g., the one of Fig. 1) are then be defined by integrating in time the mean and fluctuating velocities, respectively.
To verify that the averaging scales L and T are well chosen, one can check that the ensemble average autocorrelation function of the fluctuating velocities rapidly decreases and remains close to 0 for long time intervals τ , as shown in Fig. 3 for L = 1000 km and T = 250 days and the IABP dataset from 1979 to 2011.The ensemble average autocorrelation function χ is defined as where T q max is the duration of each trajectory and C q is the Lagrangian normalized autocorrelation function for each individual trajectory, which is defined as Here u q (t) is the fluctuating velocity of the buoys q at time t and u 2 is the variance of the fluctuating velocities for the whole trajectory, which is defined as

Application of Taylor's diffusion theory
Once the motion is decomposed into mean and fluctuating parts, it is possible to analyze the diffusion properties in the medium by following the theory developed by Taylor (1921) for turbulent fluids.Taylor's diffusion theory is valid for statistically steady and homogeneous turbulent flow without mean flow and whose fluctuating velocity follows a Gaussian distribution.When following a single particle in such conditions, the variance of its fluctuating displacement r 2 (t) = r (t 1 )r (t 1 + t) after a time interval t should be equal to where t 1 is any instant of time in the life time of the particle and the variance of the fluctuating velocity u 2 is constant in time.Note that the subscript q is dropped when dealing with only one particle.
For very long time intervals τ , the autocorrelation vanishes (as shown in Fig. 3) and the integral of C(τ ), is then a constant referred to as the Lagrangian integral timescale.Since we cannot integrate this equation to infinity, the average integral timescale is often computed as where t 0 is the first time χ (τ ) crosses zero (see for instance Poulain and Niiler, 1989;Rampal et al., 2009b).In the example of Fig. 3, t 0 = 11 days, meaning that fluctuating velocities are uncorrelated for larger time intervals, and the integral timescale is equal to 1.71 days.determines the transition between two diffusion regimes.For times much shorter than , the particle is in the ballistic regime and Eq. ( 11) becomes This simply results from the fact that C(τ ) tends to the limiting value unity for small time intervals.For times much longer than , the particle is in the Brownian regime (also called "random walk" regime) and Eq. ( 11) becomes where α is a constant defined as Casce, 2008).We checked that this constant term is very small and can be neglected.This second regime is similar to the one driven by molecular diffusion, i.e., where fluctuating velocities are uncorrelated.
These two regimes are clearly detected in Fig. 4 where is plotted the displacement variance r 2 (t) computed with L = 1000 km and T = 250 days and the IABP dataset from 1979 to 2011, as well as the corresponding solution for the ballistic and Brownian regimes (Eqs.14 and 15).The fluctuating displacements have been computed as the integral in time of the fluctuating velocities for segment periods of 35 days.
Following Lagrangian turbulent theory, diffusivity is defined as In the ballistic regime (with Eq. 14), diffusivity increases with time and may be calculated as The green lines correspond to Eqs. ( 14) and ( 15) and are shown for reference.
In the Brownian regime (with Eq. 15), diffusivity (also called eddy diffusivity in that case) is constant and may be calculated as Note that the values of diffusivity for turbulent fluids are generally much larger than diffusion coefficients linked to molecular diffusion.

Results
To define the optimal averaging scales, we perform the decomposition of the sea-ice motion and the diffusion analysis for scales ranging from L = 100 to 2000 km by step of 100 km and from T = 50 to 250 days by step of 50 days, and we analyze the deviation (i.e., the root mean square difference) of the obtained displacements variance (black line in Fig. 4) to the theoretical solution for the Brownian regime (from t = 20 to 35 days) given by Eq. ( 15) (upper green line in Fig. 4).The root mean square difference as a function of L and T shows several local minima (not shown).The scales L = 1000 km and T = 250 days are defined as optimal as they correspond to the local minima with the smallest averaging spatial scale.With these optimal averaging scales, the computed fluctuating displacements variance r 2 (t) fits well with Taylor's theory.As in Rampal et al. (2009b), we find a clear transition between the two diffusion regimes (indicated by the green lines in Fig. 4).In the ballistic regime the displacement variance grows with t 2 , whereas in the Brownian regime the displacement variance grows with t.The timescale at which the www.the-cryosphere.net/10/1513/2016/The Cryosphere, 10, 1513-1527, 2016 regime transition occurs corresponds to the integral timescale = 1.71 days.Also the magnitude of the fluctuating displacements variance compares well with the asymptotic values (indicated by the green lines in Fig. 4) predicted by the theory for t (from Eq. 14) and for t (from Eq. 15).We also checked the stationarity of the variance of the fluctuating velocities by comparing u 2 (t) to the mean values u 2 .Having an almost stationary variance was found crucial to having a good match between the computed displacement variance and the asymptotic values predicted by the theory (not shown here).To increase the robustness and statistical significance of the diffusion analysis, we then artificially increase the number of buoy trajectories by splitting each trajectory into 35-day segments starting every 12 h, i.e., every time a new buoy position along-track is available.By doing so, we make sure that the variance of the fluctuating velocities u 2 (t) , where t here goes from 0 to 35 days, is almost constant, with a relative deviation to the mean values u 2 not larger than 10 %.
These evaluation steps ensure that the values given for u 2 , and K (see Table 1) are consistent with the theory and the analyzed data and can then been used with confidence.

Discussion
The first outcome of the diffusion analysis is to provide a simple and rigorous method to separate the mean circulation from the fluctuating motion, which can then be analyzed separately.
A second outcome is to quantify the diffusion properties of sea ice that can then be compared to the diffusion properties of passive tracers in the ocean.We note that the integral timescale (about 1.71 days) and the diffusivity K (1.17 × 10 3 m 2 s −1 ) are of the same order of magnitude as those found for ocean drifters (e.g., Poulain and Niiler, 1989;Zhang et al., 2001).
A third outcome of this analysis is to give a way to estimate the evolution of the fluctuating displacement r .The magnitudes of diffusivity and integral timescale can be used to evaluate the fluctuating displacements variance (and its standard deviation) for any time t with Eqs. ( 14) and ( 15).
The standard deviation of the fluctuating displacement is a crucial piece of information for the planning of a recovery operation in the case of drifting oil or some other pollutant that is trapped in or attached to the ice, as it gives an estimate of how the size of the search area around the predicted mean drift should increase with time, in a statistical sense.This is illustrated in Fig. 5, which shows the evolution of the norm of the fluctuating displacement for every thousandth segment retrieved from the IABP trajectories for the winter periods from 1979 to 2011.This norm indicates the distance between a given buoy and the trajectory predicted by the mean drift.About 67, 96 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.
For forecasts longer than a few days, typically only the mean ice drift can be trusted.Then, the long-term average standard deviation provided here could be used to define the size of the search area around the position predicted by the mean drift.For example, the search area could be defined as a circular region with a radius equal to 3 standard deviations of the fluctuating displacement.The search radius would then be about 84 km after 5 days (corresponding to a surface area of 22 200 km 2 ) and about 222 km after 30 days (corresponding to a surface area of 154 900 km 2 ).More examples are given in Table 2.
The diffusion analysis may also be used to predict longterm (typically seasonal) sea-ice trajectories based on continuous or discrete tracer models.The average mean velocity needed by the tracer model may be defined from observations taken over the last few months, whereas the term reflecting the effects of the unpredictable fluctuations may be defined by the values of diffusivity and integral timescale derived from the diffusion analysis.
Compared to the analysis of the same data performed by Rampal et al. (2009b), several improvements should be highlighted.First of all, the mean velocity field is defined with a simpler weighted average which does not depend on the rank or the distance of the observations to the target point.Second, the criterion to define the optimal averaging scales does not depend on an arbitrary criterion of convergence but rather on the minimization of the deviation from Taylor's diffusion Table 2.This table gives an estimate of the search radii and areas corresponding to 1, 2 and 3 standard deviations, respectively, and for time horizons ranging from 1 to 30 days.These numbers are averaged over the whole domain and period analyzed in Sect. 2 and should be reevaluated for specific applications, periods and domains of interest, for example by using model outputs having passed the evaluation test proposed in Sect.3. Note that we checked that about 67, 96 and 99.6 % of the fluctuating displacements are smaller than 1, 2 and 3 standard deviations.theory.Finally, we checked that the results presented here fit with Taylor's diffusion theory.In Rampal et al. (2009b), the fluctuating displacement are underestimated by a factor 100.This mistake, which the authors (Rampal et al., 2009b) are aware of, simply comes from a wrong conversion factor and has no impact on the rest of their study, but it is worth mentioning here, especially with regard to the present discussion.Note that the value of diffusivity given in Rampal et al. (2009b) is twice as low as the one found here and is not consistent with the other results presented in their study.We do not know the reason for this inconsistency.
There are, however, some limitations when using the mean motion and mean diffusivity to force passive tracer models: the averaging smooths out local mean circulation features such as coastal currents; the method is not well suited for studying dispersion as it assumes no spatial correlation; the diffusivity could differ spatially and be not well represented by the basin-wide mean value; the diffusivity could be affected by the long-term trend in the mean speed identified in Rampal et al. (2009a) and may therefore not be well represented by a single mean value derived from these decades.
These issues do not occur when the tracer models are directly forced with sea-ice velocity fields representing correctly both the mean and fluctuating parts of the motion field, as well as their gradient at all scales.The diffusion analysis presented here can also help to assess the representation of the mean and fluctuating parts.In the following section such an assessment is carried out on the results of the two sea-ice(-ocean) modeling platforms, neXtSIM and TOPAZ.
3 Diffusion analysis on observed and simulated sea-ice trajectories Model output or reanalyses are often used to directly force passive tracer models (e.g., Nudds et al., 2013;Gearon et al., 2014).In that case, it is important to check whether the simulated trajectories adequately represent the mean and fluctuating parts of the sea-ice motion before pursuing an analysis of these trajectories.In some cases, a specific term is added to the tracer model (either via a diffusive term or a random perturbation) to represent the effect of the unresolved physics on the tracer evolution.When applied in the ocean, the stochastic part represents molecular and turbulent diffusive processes that are not included in the velocity fields simulated by the ocean model.However, as sea ice, especially compact pack ice, does not behave as a turbulent fluid, the stochastic term to be used for sea ice should not be taken to represent the same underlying physical processes as in the ocean, and it therefore may have a different form or be scaled with a different coefficient.In this section, we follow an approach that is frequently used for oceanic drifters (e.g., Dominicis et al., 2012).This approach consists of applying the P. Rampal et al.: Diffusion in the ice pack: comparison between models and observation diffusion analysis to observed and simulated trajectories to determine how the mean and fluctuating motion are represented and how unresolved physics could be taken into account.

Observed and simulated trajectories datasets
In this section we compare observed trajectories from the IABP dataset to trajectories of virtual buoys (here called "floats") whose motion is forced by sea-ice fields coming from two different model setups.Due to limited available computational time, this analysis is restricted to three consecutive winters.The period 2007-2010 has been selected for its relatively good data coverage, with more than 40 IABP buoys recording their positions simultaneously every day.The float simulations are initialized at the same time and position as the IABP buoys (280 individual floats).The positions of each float are sampled every 12 h at the same time as the IABP buoys and stop when the IABP buoy track stops or when the float enters into an area of simulated open water (sea-ice concentrations less than 15 %).By doing so, three comparable datasets with the exact same number of positions are obtained: (i) the observed sea-ice trajectories, already discussed in Sect.2.1; (ii) the trajectories of virtual sea-ice floats forced by a free run of the TOPAZ sea-ice-ocean data assimilation system; and (iii) the trajectories of virtual sea-ice floats simulated by the neXtSIM sea-ice model.

TOPAZ trajectories dataset
TOPAZ is a coupled sea-ice-ocean data assimilation system (Sakov et al., 2012) used in the operational Arctic Ocean forecast platform of the European Copernicus Marine Environment Monitoring service (http://marine.copernicus.eu).It has also been used to build a 23-year reanalysis , also distributed by the Copernicus marine service.The ocean part of TOPAZ uses the HYCOM model version 2.2, with 28 vertical layers, whereas the sea-ice part uses a singlethickness-category sea-ice model whose thermodynamics are described in Drange and Simonsen (1996) and dynamics are built around a standard EVP rheology (Hunke and Dukowicz, 1997) as it was implemented in the CICE sea-ice model (the Los Alamos sea-ice model) version 4 (Hunke and Lipscomb, 2010).
The sea-ice-ocean model of TOPAZ (hereafter called the TOPAZ model) is used here in free-run mode (i.e., no data assimilation is applied) in the same configuration as in Sakov et al. (2012).The model grid covers the Arctic and North Atlantic oceans with a mean resolution of approximately 12 km over the Arctic.The three TOPAZ simulations analyzed in the following start on 15 September and finish on 15 May for three consecutive winters from 2007 to 2010 with initial conditions extracted from the free-run simulation described in Sakov et al. (2012).The applied atmospheric forcing fields are the 10 m wind velocity, the 2 m temperature and mixing ratio, mean sea level pressure, total precipitation and the fraction of that which is snow, and the incoming short-wave and long-wave radiation from the ERA interim reanalysis (ERAi) distributed at 80 km spatial resolution and 6 h temporal resolution (http://www.ecmwf.int/en/research/climate-reanalysis/era-interim, ECMWF, 2011).
The TOPAZ model in free-run mode has been evaluated in Sakov et al. (2012) and was found to overestimate seaice drift by about 3 km day −1 compared to buoy data.To try to solve this issue, the frictional drag parameters for the atmosphere-ice stress has been reduced to c a = 0.0016 in the TOPAZ operational platform (not yet documented).The same value is also used here.The mean sea-ice thickness in the free run is generally underestimated and has reduced horizontal gradients.The sea-ice thickness is slightly better in the TOPAZ reanalysis, which applies assimilation of sea surface temperature, in situ profiles, sea-ice concentration and sea-ice drift, but the total volume is still too low (Sakov et al., 2015, http://marine.copernicus.eu/documents/QUID/CMEMS-ARC-QUID-002-003.pdf).
The float tracking with the TOPAZ model is performed offline by using the hourly mean sea-ice velocity fields simulated by the model.The float-tracking system moves the floats with a simple Eulerian method.The virtual floats move in the quasi-homogeneous TOPAZ Arctic grid in order to avoid singularity errors at and around the North Pole that would arise with a regular longitude / latitude grid.The seaice velocities given by the TOPAZ model are interpolated with a bilinear method to the position of the virtual Lagrangian floats every hour.We checked that, for the timescale and spatial resolution considered here, this tracking method gives similar results to computing the float position during runtime with the advantage of remaining computationally efficient.

neXtSIM trajectories dataset
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., 2016).Thermodynamic growth and melt of the ice are based on the zero-layer model of Semtner (1976) and the ice model is coupled to a slab ocean model, whose variables are the slab ocean temperature and salinity.While still being under development, the model is already used in an experimental sea-ice forecast platform covering the Kara Sea (https://www.nersc.no/data/nextsim-f).
The configuration used here is the same as the one presented and evaluated in Rampal et al. (2016) neXtSIM simulations used here start on 15 September and finish on 15 May, for three consecutive winters from 2007 to 2010.The model is initialized with the ice concentration derived from the 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 TOPAZ reanalysis ice thickness, within the area covered with ice.As the modeled ice volume of the TOPAZ reanalysis is known to be too low (Sakov et al., 2015), we increased the initial thickness uniformly so that the total volume is the same as that given by the PIOMAS model (Zhang and Rothrock, 2003)  ).The ASR-Interim is a high-resolution atmospheric reanalysis (30 km with output every 4 h) known to reproduce particularly well the nearsurface wind fields in the Arctic region (Bromwich et al., 2016).
The neXtSIM model is able to simulate correctly the observed evolution of the sea-ice volume, extent and area for the freezing season (from September to May) but simulates a too-rapid melt from May onwards (Rampal et al., 2016).This limitation does not impact the results of this study as we only analyze simulated drift in winter with simulations restarted every September.
For the winter season 2007-2008, the simulated drift fields have been extensively evaluated in Rampal et al. (2016) against satellite derived products, showing a high correlation (higher than 0.85), only a negligible bias in the 3-day drift and a good representation of the mean circulation.
The float tracking with neXtSIM is performed at runtime.The main reason for doing this is that the Lagrangian advection used in the neXtSIM model offers some additional challenges to a post-processing approach using Eulerian fields due to the remeshing technique applied (see Rampal et al. (2016) for more details on the remeshing procedure).

Results
Figure 6 shows maps with all the buoy / float trajectories for each winter season from the IABP, TOPAZ and neXtSIM ice trajectories datasets.The three winter subsets from 2007 to 2010 are much sparser and cover a smaller portion of the Arctic Ocean than the reference dataset of 1979-2011 analyzed in Sect.2.1.As it is not practical to directly compare the simulated and observed trajectories, we analyze separately the mean drift and the fluctuating part of the motion by applying the decomposition method presented in Sect.2.2.
The mean velocity fields for each winter season and for the three datasets are shown in Fig. 7.The two main features of the Arctic-wide mean circulation are the Beaufort Gyre and the Transpolar Drift.The lowest drift speeds are observed along the Canadian Arctic Archipelago, where the ice is significantly thicker and more ridged.The strength and the extent of the Beaufort Gyre, as well as the strength of the Transpolar Drift, vary from one year to the next.This interannual variability is well represented by both TOPAZ and neXtSIM.The quality of the mean drift simulated by the two models is constant in neither time nor space.In 2007/2008, the simulation with TOPAZ largely misses the low drift speed along the Canadian Arctic Archipelago, when the simulation with neXtSIM correctly reproduces it.For the other winters, it is less obvious to distinguish clear differences in the quality of the simulated mean drift fields.
The statistical distribution of the mean velocity gives valuable information and can be used to evaluate the simulated mean drift more objectively.Figure 8 shows the probability density function of the mean speed Ū = √ ū2 + v2 as computed from the IABP buoy data and from the TOPAZ and neXtSIM virtual buoy data for the period 2007-2010, as well as the Gaussian and exponential fits.None of the three distributions fit well with the Gaussian or exponential distributions on the whole range of values.The mean speed distribution obtained with the TOPAZ dataset seems to have not enough values within the range of 0 to 1.5 cm s −1 and too many values within the range of 1.5 to 4 cm s −1 .The mean speed distribution obtained from the neXtSIM dataset is close to the observed distribution within the range of 0 to 4 cm s −1 but differs for the larger values.The mean value of the observed mean speed is equal to 1.95 cm s −1 , when the one of the TOPAZ dataset is 2.89 cm s −1 , which is about 48 % larger than the observations, and the one from the neXtSIM dataset is 1.65 cm s −1 , which is about 15 % lower than the observations.
When removing the mean part of the velocity field we are left with the fluctuating velocity field u (x).If the mean part is removed correctly (according to the Taylor's theory), each component of the fluctuating velocity should be symmetrically distributed around zero.This is the case in our results (not shown), meaning that one can directly use the norm of the fluctuating velocity instead of the components separately, without losing information.The probability density functions (PDFs) of the fluctuating speeds are plotted in Fig. 9 with the Gaussian and exponential fits indicated for reference.
The fluctuating speeds of the IABP buoys clearly follow an exponential distribution with a mean equal to 6.98 cm s −1 .The fact that the data follow an exponential instead of a Gaussian distribution means that the sea-ice fluctuating speed can be much larger than a standard deviation away from the (zero) mean.Such non-Gaussian distributions for fluctuating speeds are not expected for fully developed turbulence (Batchelor, 1960;Frisch, 1995) but have been observed for oceanic surface currents during energetic events associated with large organized structures such as jets and vortices.Such a signature for multi-year 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 mechanics (Rampal et al., 2009b).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 mean value of the fluctuating speeds from the TOPAZ setup is too high by about 30 % (8.97 cm s −1 instead of 6.98 cm s −1 ).Also their statistics do not follow an exponential distribution as in the observations but are closer to the Gaussian fit centered on the value 10 cm s −1 , resulting in an underrepresentation of both the very high speed (larger than 30 cm s −1 ) and very low speed (smaller than 5 cm s −1 ).The mean value of the fluctuating speeds from neXtSIM is too low by about 10 % (6.2 cm s −1 ).The fluctuating speeds follow an exponential distribution but only in the velocity range 0 to 30 cm s −1 .The high speed values are missed by the two models.
Figure 10 shows the evolution of the fluctuating displacements variance for the observed and simulated trajectories.constantly overestimated with TOPAZ, by about 40 % in the ballistic regime and almost 50 % in the Brownian regime.This is consistent with the fact that the integral timescale and fluctuating velocity variance are overestimated by 15 and 30 %, respectively (see Table 1).In the Brownian regime, these overestimation are combined resulting in the overestimation of the diffusivity by about 50 %.With neXtSIM the fluctuating displacements are underestimated in the ballistic regime but correctly reproduced for the Brownian regime.This is consistent with the underestimation of the fluctuating velocities variance by about 23 %, which is balanced in the Brownian regime by the overestimation of the integral timescale by about 25 %, leading to a diffusivity almost equal to the one coming from the observations for the same period.
Figure 11 shows the regional distribution of the diffusivity fields.The results obtained with the TOPAZ setup have a correlation coefficient against the diffusivity map obtained from observations of about 0.71 but generally overestimate the diffusivity.The results obtained with the neXtSIM setup have a correlation coefficient against diffusivity map obtained from observations of about 0.85 and adequately represent the magnitude of the diffusivity.
The diffusivity field computed from the IABP buoys is not uniform and seems to be related to the spatial distribution of sea-ice thickness shown in Fig. 12, with rather low diffusivity values along the Canadian Arctic Archipelago (i.e., as small as about 0.5 × 10 3 m 2 s −1 ) and larger values in the Beaufort and East Siberian seas (1.5-2.5 × 10 3 m 2 s −1 ).The spatial distribution of the diffusivity obtained with the TOPAZ and neXtSIM setups also correlate well with the simulated seaice thickness pattern (shown in Fig. 12).

Discussion
The goal of the present analysis is not to compare the model systems themselves but to illustrate how the simulated motion fields differ from observations and what would be the impact of using such model outputs to force passive tracers models to study, for example, trajectories of pollutant trajectories in sea ice.The differences between the simulated and observed motion may be due to many factors, ranging from the internal characteristics of the sea-ice models (their rheology, drag parameterization, etc.) to external causes, such as the initial conditions, atmospheric forcing and impact of the ocean.To distinguish the effects of each factor would require one to run the same model with different initial conditions, forcings and set of parameters or to run different models in the same configuration (initial conditions, forcings, parameters, etc.).Other diagnostics than the diffusion analysis would also be necessary.For example, the effect of the rheology would be better analyzed by a dispersion analysis (double particle diffusion) as in Rampal et al. (2008) as it is directly related to sea-ice deformation.Nevertheless, even when the present analysis cannot clearly distinguish the sources of the differences between the simulated and observed trajectories, it provides pertinent information on the quality of the simulated trajectories.
The TOPAZ model reproduces the very basic characteristics of the Arctic sea-ice mean circulation, with interannually varying Beaufort Gyre and Transpolar Drift.However, the averaged mean drift is overestimated by about 48 %.This overestimation partly comes from missing the mean drift speed smaller than 5 cm s −1 that are, in the observations, localized north of the Canadian Arctic Archipelago (CAA) in a region covered by thick multi-year ice.
In the simulations analyzed here, TOPAZ misses both the low (larger than 10 cm s −1 ) and high (larger than 30 cm s −1 ) values of fluctuating speed.As the low values largely dominate the observed distribution, missing those values leads to an overestimation of the fluctuating velocities mean and variance.When combined to the overestimation of the integral timescale, it leads to a large overestimation (by about 50 %) of the long-term fluctuating displacement and absolute diffusivity.
The neXtSIM model also reproduces the interannually varying Beaufort Gyre and Transpolar Drift.The statistical distribution of the mean circulation is close to the one obtained from the observations for the range of 0 to 4 cm s −1 , but the largest values of the mean circulation (mean drift larger than 4 cm s −1 ) are underestimated, resulting in an underestimation of about 15 % of the averaged mean drift.
In the simulations analyzed here, neXtSIM represents well the statistical distribution of the fluctuating speed until 30 cm s −1 but misses the higher values, leading to an underestimation of the fluctuating velocities mean and variance.When combined to the overestimation of the integral timescale, this leads to long-term fluctuating variance and diffusivity having the same magnitude as in the observations.
The two model setups used here have a common deficiency at representing the fluctuating speed higher than 30 cm s −1 .This likely comes from the missing local and rapidly varying winds in the atmospheric reanalyses (ASRinterim and ERA-interim) used here to force the models.This lack of variability in the forcings may also explain the overestimation of the integral timescale seen in the two setups.These too long integral timescales may also partly come from missing physics or too weak coupling (e.g., too low frequency in the coupling and forcing) in the model setups that may lead to a misrepresentation of the inertial oscillations and impact 12-hourly drift statistics.
It is common practice to add an extra diffusive term to the tracer evolution equation, as discussed earlier.In the case of the TOPAZ setup, adding such an extra term to the tracer evolution equation would not help, as the model already overestimates the fluctuation both in the ballistic and Brownian regimes.Adding such a term when using the neXtSIM setup presented here may improve the evolution of the fluctuating displacement in the ballistic regime.Adding a random term would increase the fluctuating velocity variance but decrease the integral timescale.This may allow us to maintain the good performance in reproducing the long-term displacements and diffusivity fields, without impacting the long-term mean drift.

Summary and conclusions
In the first part of the paper (Sect.2), we analyze IABP buoys trajectories for the winter periods between 1979 and 2011 and we estimate the values of the integral timescale (about 1.7 days), the 12-hourly fluctuating velocities variance (59 km 2 day −2 ) and the diffusivity (1.17 × 10 3 m 2 s −1 ).We additionally verify that the computed displacement variance is consistent with Taylor's diffusion theory.
This information can be used in the context of pollutant tracking to evaluate the proper size for the search area around the long-term trajectory predicted by the mean drift.If one defines the search area as a circular region with a radius equal to 3 standard deviations of the fluctuating displacement, which would include the polluted area with 99.6 % confidence, we find that on average the search radius should be about 84 km after 5 days (corresponding to a surface area of about 22 000 km 2 ) and about 222 km after 30 days (corresponding to a surface area of 155 000 km 2 ).Before using these estimates, one should be reminded that they are given for the whole Arctic basin and for the whole period 1979-2011 and may then differ from estimates computed for specific regions and time periods.Also as the IABP buoys are mainly deployed on multi-year ice, these estimates may not be valid for seasonal ice and for future applications.
The estimates of the mean drift field, diffusivity and integral timescale computed here could also be used within a passive tracer model (either with an advection-diffusion equation or a Lagrangian stochastic approach) to estimate the probability for a particle to be in a given position after a given time.The limitations of that approach would be the excessive smoothing of local mean circulation patterns (e.g., coastal currents), the inability to represent dispersion and the potential misrepresentation of the spatial and temporal distribution of the diffusivity values.
In the second part of the paper (Sect.3), we analyze trajectories of virtual buoys whose motion is forced by simulated sea-ice velocity fields.This approach eliminates the limitations of using tracer models forced by mean fields but relies on the good representation of the sea-ice drift by the models.To illustrate how one could evaluate sea-ice model output before using it for trajectory modeling, we applied the diffusion analysis to three similar datasets, one from the IABP data, one from the TOPAZ model and one from the neXtSIM model, and we compared the numbers obtained from simulated and observed trajectories.
The mean velocities in the simulations using TOPAZ are on average 50 % too high and generally miss the very low mean drift located near the Canadian Arctic Archipelago.The long-term displacement variance and absolute diffusivity are also overestimated by about 50 %.Using the output of this TOPAZ setup for tracer studies would produce too long trajectories and too large displacement variance, potentially affecting the conclusions of such studies.
The mean velocities in the simulations using neXtSIM are on average 15 % too low; they reproduce well the very low mean drift located near the Canadian Arctic Archipelago but miss the highest values of the mean motion.The long-term displacement variance and absolute diffusivity fit well the observations.Tracer studies based on such results could be trusted except for the ballistic regime (first few days), where the simulated displacement variance is too weak.
Using the outputs of the simulations made with neXtSIM would give better sea-ice trajectories than using the outputs of the TOPAZ simulations analyzed here.However, whether this difference mainly originates from (a) the different initial conditions, forcing and ocean or (b) the sea-ice model itself (different rheologies and thermodynamics) cannot be clearly answered.As a follow up of this study, it would be interesting to investigate the causes of the missing high values of fluctuating velocities and the overestimation of the integral timescale, first by examining the impact of the atmospheric forcing resolution and second by checking how the inertial/tidal oscillations are represented by the two modeling platforms used here.To better asses the quality of the simulated sea-ice dynamics, it would be interesting to also perform a dispersion analysis as in Rampal et al. (2008) or to specifically study sea-ice deformation as in Bouillon and Rampal (2015) with data from models and observations.Finally, it is worth noting that the representation of the mean sea-ice circulation depends on many processes (mean circulation of the atmosphere and ocean, spatial and temporal variation of the ocean-ice and air-ice drag coefficients as a function of the ice age/type, representation of the ocean-ice and atmosphere-ice boundary layers (McPhee, 2012), etc.), the respective roles of which would also need to be further explored.

Figure 1 .
Figure1.Example of a 35-day long trajectory from a buoy of the IABP dataset (thin red line) partitioned into a mean (thick black line) and fluctuating (thin black line) parts using the method described in Sect.2.2 and the averaging scales L = 400 km and T = 150 days.The starting point of the three trajectories is the same and arbitrarily set to be the origin of the axes.

Figure 2 .
Figure 2. Buoy tracks from the IABP dataset for the winter periods 1979-2011 (left panel) and the corresponding number of buoys (middle panel) and records (right panel).The central Arctic domain is delineated by a blue line.

LFigure 3 .
Figure 3. Ensemble average autocorrelation function of the fluctuating velocities computed from the IABP dataset for winter seasons 1979 to 2011 with three different averaging scales.The first zerocrossing time t 0 is given for indication.

Figure 4 .
Figure 4. Ensemble mean of the variance of the fluctuating displacements < r 2 > for the winter seasons 1979-2011 for IABP.The green lines correspond to Eqs. (14) and (15) and are shown for reference.

Figure 5 .
Figure 5.Time evolution of the norm of the fluctuating displacement r for every thousandth 35-day segments extracted from the IABP buoys tracks for the winters 1979-2011 (left).The solid lines in color indicate 1, 2 and 3 standard deviations of the fluctuating displacement, respectively.Illustration of search area after 35 days (right) estimated statistically from the IABP buoys dataset for the period 1979-2011.

Figure 10 .
Figure10.Ensemble mean of the variance of the fluctuating displacement < r 2 > for the winter seasons 2007-2011 for IABP, TOPAZ and neXtSIM.The dashed lines indicate the theoretical slopes for the two diffusion regimes: the ballistic regime ( r 2 (t) ∼ t 2 ), which has a slope equal to 2, and the Brownian regime ( r 2 (t) ∼ t), which has a slope equal to 1.

Figure 11 .
Figure 11.Diffusivity fields obtained from the analysis of the IABP buoys trajectories (left), TOPAZ float trajectories (middle) and neXtSIM float trajectories (right) for the winters 2007-2010.The diffusivity is averaged over boxes of 400 by 400 km.

Figure 12 .
Figure 12.Mean ice thickness in the central Arctic obtained from ICESat satellite observations (Kwok et al., 2009) and the two models.ICESat results are only available in February and March of 2008, while the model results have been averaged over the winters of 2008, 2009 and 2010.

Table 1 .
Taylor (1921)ows the total number of floats (Nrf), the calculated integral timescale ( ), the variance u 2 and the calculated diffusivity K for the different dataset (IABP, TOPAZ and neXtSIM) and time periods used in this study.All these Lagrangian statistics were computed following the diffusion theory ofTaylor (1921)and using L = 1000 km and T = 250 days as averaging scales to calculate the Lagrangian mean velocities.