Sea ice inertial oscillations in the Arctic Basin

An original method to quantify the amplitude of inertial motion of oceanic and ice drifters, through the introduction of a non-dimensional parameter M defined from a spectral analysis, is presented. A strong seasonal dependence of the magnitude of sea ice inertial oscillations is revealed, in agreement with the corresponding annual cycles of sea ice extent, concentration, thickness, advection velocity, and deformation rates. The spatial pattern of the magnitude of the sea ice inertial oscillations over the Arctic Basin is also in agreement with the sea ice thickness and concentration patterns. This argues for a strong interaction between the magnitude of inertial motion on one hand, the dissipation of energy through mechanical processes, and the cohesiveness of the cover on the other hand. Finally, a significant multi-annual evolution towards greater magnitudes of inertial oscillations in recent years, in both summer and winter, is reported, thus concomitant with reduced sea ice thickness, concentration and spatial extent.


Introduction
The spectacular evolution of the Arctic sea ice cover over the last few decades is not restricted to the shrinking of the ice extent (Comiso et al., 2008;Stroeve et al., 2008), its thinning (Rothrock et al., 2008;Kwok and Rothrock, 2009), or, consequently, a continued decline of the ice volume (Lindsay et al., 2009).Kinematics is affected as well, and its evolution plays a central role in the changes currently taking place in the Arctic Ocean.As observed from buoy drift data, the sea ice mean speed over the Arctic increased at a rate of 9 % per decade from 1979 to 2007, whereas the mean deforma-tion rate increased by more than 50 % per decade over the same period (Rampal et al., 2009).These two aspects of recent sea ice evolution, i.e. strong decline in terms of ice extent and thickness, and accelerated kinematics, are strongly coupled within the albedo feedback loop.Increasing deformation means increasing fracturing, hence more lead opening and a decreasing albedo (Zhang et al., 2000).As a result, ocean warming, in turn, favours sea ice thinning in summer and delays refreezing in early winter, i.e. strengthens sea ice decline.This thinning should decrease the mechanical strength, therefore allowing even more fracturing, hence larger speed and deformation.A consequence is the acceleration of the export of sea ice through Fram Strait, with a significant impact on sea ice mass balance (Rampal et al., 2009(Rampal et al., , 2011;;Haas et al., 2008), and ice age (Nghiem et al., 2007).Moreover, sea ice mechanical weakening decreases the likelihood of arch formation along Nares Strait, therefore allowing old, thick ice to be exported through this strait (Kwok et al., 2010).
The principle of a strong interaction between the ice state (concentration, thickness), on one hand, and the mechanical behaviour of the cover and its strength on the other hand, although rather intuitive, needs however to be quantified more precisely.Beyond the trivial reduction of strength in proportion to the thinning of an ice plate, one can expect a more complex effect of the compactness of the ice cover on its average strength, e.g. as in granular media (Rajchenbach, 2000;Aranson and Tsimring, 2006).This is the combination of these two effects that we aim to explore here.Ideally, one would like to check whether the sea-ice cover responds differently from year to year to the same mechanical forcing.We here propose to tackle this problem by analysing its response, as deduced from ice drifter trajectories, to the specific inertial forcing.The effect of the Coriolis force on geophysical fluids dynamics has been studied for more than a century.Interestingly, the first studies of oceanic inertial oscillations (Ekman, 1905) were prompted by the observations of Nansen, made during the Fram's journey along the Transpolar drift, that sea ice was moving with a 20 • -40 • angle to the right of the wind direction (Nansen, 1902).Indeed, as the Coriolis force acts perpendicularly to the particle velocity, it induces a deviation of the trajectory to the right in the Northern Hemisphere.This deviation generates inertial oscillations, characterized by a frequency of where φ is the latitude, i.e. close to a semi-diurnal frequency (2 cycles day −1 ) in the Arctic.Within the ice-free ocean, these oscillations are triggered by wind forcing events such as storms or moving fronts (Price, 1983;Gill, 1984;D'Asaro et al., 1995).These events generate inertial oscillations within a few hours, which then decrease progressively with a characteristic time scale of a few days (Park et al., 2009).This damping is essentially due to the internal friction within the Ekman ocean layer, as well as to the radiation of inertial waves toward the thermocline.
On ice-covered oceans, we expect other processes to strengthen the damping of these oscillations, such as friction at the ice/ocean interface or, more importantly, the internal ice stresses resulting from solid mechanical interactions within the ice cover such as fracturing, friction between adjacent floes and shearing of leads, or crushing during convergent deformation and ridge formation, i.e. the "ice internal friction" (Lepparanta, 2004;Colony and Thorndike, 1980).These terms are apparent in the momentum balance of sea ice dynamics: where D/Dt = ∂/∂t + U i , is the Lagrangian time derivative, ρ i the ice density, h i the ice thickness, U i the ice velocity vector expressed in Cartesian coordinates, σ the internal stress tensor, τ a and τ w respectively the wind and oceanic "stresses" (forces per unit ice area).The last term on the left hand side of Eq. ( 2) is the Coriolis force, with f 0 the inertial oscillation frequency, expressed in cycles day −1 , = 2π / 24 rad h −1 = 7 × 10 −5 rad s −1 the Earth rotation velocity and k a unit vector aligned along the South to North Pole axis.
Here, we neglect the contribution of the sea surface tilt to the momentum balance, which is small compared to the other contributions (Steele et al., 1997).In Eq. ( 2), the wind forcing τ a excites the oscillations, whereas the oceanic drag τ w and the internal stress term ∇ •σ damp those oscillations (Colony and Thorndike, 1980).
We expect the response of sea ice to the Coriolis force in the frequency domain to be within two following extremes.A buoy moving in free drift, i.e. fixed on an ice floe that drifts according to wind and ocean currents without any mechanical interaction with other ice floes (the term ∇ •σ is then negligible compared to the others), is expected, in first order approximation, to follow the oceanic fluid parcel and thus to oscillate in a similar way, although another source of damping of the oscillations might come from the friction between the ice bottom and the ocean surface (Lepparanta et al., 2012).In contrast, on a compact ice cover experiencing strong internal stresses, the corresponding contribution ∇ • σ will dominate the other terms in Eq. ( 2) so that the oscillations are immediately damped out and thus become undetectable.
To measure the amplitude of the inertial oscillations of the ice drifters is therefore to estimate the level of mechanical dissipation within the ice cover, and therefore its degree of cohesiveness.As such, its evolution and its spatial pattern can highlight the changes in ice conditions and ice cover mechanical behaviour.Internal ice stress measurements have been performed directly from in situ sensors (Richter-Menge and Elder, 1998;Richter-Menge et al., 2002).The seasonal variations of these local stress measurements have shown an increase of ice-motion induced stresses as the winter season progresses and the cohesiveness/compactness of the ice cover develops (Richter-Menge and Elder, 1998).While of major importance to analyze sea ice mechanical behaviour and rheology (Weiss et al., 2007), these measurements are however limited to the local scale and do not carry any information about a possible large-scale, long-term trend of the mechanical state of the cover.Our approach is thus complementary, as it allows a much better sampling, both in time and space.
The idea to relate the amplitude of the inertial oscillations to the degree of consolidation of the ice cover was already formulated by McPhee (1978) as well as Colony and Thorndike (1980).The latter studied the role of the mechanical behaviour of the ice on the coherency of inertial motion between different buoys, an approach different from what is described below as we analyze trajectories individually.More recently, from buoys trajectories near the Antarctic Peninsula, Geiger and Perovich (2008) observed an increase of inertial motion related to the degradation of the ice pack during the spring breakup.Hutchings et al. ( 2012) also discussed the link between the degree of fragmentation of the ice cover, the mechanical dissipation of energy, and the amplitude of inertial motion.These last authors based their discussion on the analysis of strain-rate records obtained from a dense array of buoys in the Weddell Sea, i.e. they did not analyze individual trajectories.However, although sea ice inertial oscillations have been studied by several others authors (Colony and Thorndike, 1980;McPhee, 1978;Hunkins, 1967;Lammert et al., 2009), including recently on ice strain-rate records (Kwok et al., 2003), this is the first time that a systematic analysis is performed at the Arctic Basin and multidecadal scales.
In this paper, we propose a method to quantify the magnitude of inertial oscillations from Lagrangian (buoys) trajectories (Sect.3).We then apply this methodology to the International Arctic Buoy Program (IABP) dataset covering 30 yr of data in the Arctic Ocean (Sect.4).As shown below, this analysis is in full agreement with the above expectations: inertial oscillations are very weak or absent in a highly cohesive ice cover, such as in winter within the central Arctic, but are well developed in summer at the periphery of the basin, i.e. in regions of less concentrated, loose ice.In addition, a significant strengthening (on average) of these oscillations is observed, suggesting a mechanical weakening of the Arctic sea ice cover.This is confirmed in Gimbert et al. (2012), where a simple ocean-sea ice coupled dynamical model explains these seasonal, geographical and multi-annual variations of inertial oscillation magnitude in term of changes within sea ice internal mechanical properties, through an associating decrease of sea ice internal friction in recent years.

IABP buoy dataset
The sea ice drifting buoy dataset, provided by the International Arctic Buoy Program (IABP), consists of approximately 650 buoy tracks recorded from January 1979 to December 2008.These ice drifters are dropped every year at the end of winter, mostly in April, and drift according to the ice motion.Their positions are tracked by GPS receivers or Argos transmitters with a position uncertainty of the order of 100 m and 300 m, respectively.The raw buoy positions are irregularly sampled through time.Thus, in order to get a regular sampling, the buoy positions provided by the IABP result from a cubic interpolation of the raw positions with a 3 h time re-sampling (see the IABP documentation at http://iabp.apl.washington.edu/data.htmlfor further details).This interpolation procedure acts as a low-pass filter, thus most likely reducing the position error around 100 m or below (Lindsay and Stern, 2003).Figure 1 shows all the buoy tracks on a polar stereographic map, where buoy locations are defined as xe 1 + ye 2 , where e 1 and e 2 are two orthogonal unit vectors with e 1 being along the Greenwich meridian and (x, y) = (0, 0) at North Pole.

Oceanic buoy dataset
In order to get an example of the amplitude of inertial oscillations in the absence of sea ice, i.e. when the damping of oscillations is mainly due to the internal friction of the Ekman layer without contribution from internal ice stresses, we consider trajectories of buoys drifting in the North Atlantic, along the northern Spanish coast (Fig. 2).This oceanic buoy dataset, provided by the SHOM (Service Hydrographique et Oceanographique de la Marine), consists of 4 buoy tracks (B239, B241, B242, B245).These buoys were deployed between 6 and 11 December 2006 in the framework of the CONtinental GAScogne (CONGAS) project (Le Cann and Serpette, 2009).The buoy positions were tracked by GPS receivers with a position uncertainty of the order of several tens of meters.The temporal sampling was regular, equal to 1 h.Thus, the buoy positions do not result from any re-sampling or re-interpolation.Buoys B239 and B245 operated during 1 yr, whereas buoys B241 and B242 operated only 1 month.For each latitude (φ)-longitude (ϕ) buoy position, we define the orthogonal base (e 1 , e 2 ) of this coordinate system as x (φ,ϕ) = xe 1 +ye 2 using a Lambert projection centered in the coordinates φ = 43 • and ϕ = 6 • .
3 Quantifying the magnitude of inertial oscillations

Observations
We analyze here several cases of inertial oscillations observed over an ice-free ocean or over sea ice.
Figure 3a shows 1 month of the trajectory of the oceanic buoy B245.Inertial oscillations are noticeable during this time period: being in the Northern Hemisphere, the buoy drift trajectory exhibits cycloids in the clockwise direction.The intermittent character of inertial oscillations is also clearly visible.Periods of strong inertial oscillations, marked by cycloidal loops (red boxes in Fig. 3a), might have been triggered by storms.Following these periods, the oscillations are progressively damped out and the loops nearly disappear (green boxes in Fig. 3a), before a new storm or moving front can trigger new inertial oscillations.We compute the u x ( x, ỹ, t) (along the x-axis) and u y ( x, ỹ, t) (along the y-axis) speeds as follows: with t = 1 h and where the mid-points x, ỹ and t are defined as x = (x(t + t) + x(t))/2, ỹ = (y(t + t) + y(t))/2 and t = ((t + t) + t)/2.The Fourier transform Ûb (ω) of the buoy velocities u x and u y for a selected buoy trajectory, which starts at time t 0 and ends at time t end , is defined as where N is the number of velocity samples along the trajectory and ω = 2πf .This vectorial Fourier transform distinguishes negative and positive frequencies: peaks at f < 0 and f > 0 are associated with clockwise and counter-clockwise displacements, respectively.Figure 3b shows the Fourier spectrum | Ûb (ω)| of the buoy velocities associated with the trajectory plotted in Fig. 3a.The Fourier transform reveals a peak of magnitude 5.2 km day −1 at the inertial oscillation frequency (f 0 ≈ −1.35 cycles day −1 using the average latitude φ = 1/N t φ(t) ≈ 42 • in Eq. 1).This peak corresponds to the component of the buoy motion associated with the Coriolis force.The peak at f = 0 represents the advective component of the buoy motion.Finally, the two peaks observed at f = −2 cycles day −1 and f = 2 cycles day −1 are associated with a semi-diurnal tidal oscillation.Unlike the inertial oscillation, the tidal oscillation does not rotate and the associated peak is therefore observed at positive and negative frequencies.We now consider IABP ice drifters.A buoy moving in free drift, i.e. fixed on an ice floe that drifts according to wind and ocean currents alone, hence without any mechanical interaction with other ice floes, is expected to oscillate in a way similar to e.g. the oceanic B245 buoy of Fig. 3 when subjected to similar forcing conditions, although another source of damping of the oscillations might come from the friction between the ice bottom and the ocean surface.
Figure 4a shows an IABP buoy trajectory within the Fram Strait during the summer period, where ice is highly fragmented and loosely packed.The clockwise-cycloids are clearly visible.We observe bursts in the inertial oscillation intensity (red rectangles) followed by a period of decaying intensity (green rectangle).
The velocity of the buoy is computed using Eq. ( 3), with stereographic coordinates x and y and with t = 3 h.The inertial oscillations are evidenced by a strong peak observed on the velocity Fourier spectrum (Fig. 4b), at the inertial frequency f 0 ≈ −1.97 cycles day −1 (calculated using φ = 80 • in Eq. 1).As the Arctic Basin lies between 70 • and 90 • of latitude, the inertial oscillation frequency varies from −1.88 to −2 cycles day −1 at these latitudes and is thus merged with the semi-diurnal tidal oscillation frequency.The differentiation of these two types of oscillations can be done by looking at the amplitude of the Fourier spectrum with respect to signed frequencies: we assume the spectral peaks associated with the semi-diurnal tidal oscillation to be roughly symmetric at positive and negative frequencies (as a first order approximation, as bathymetry or ocean currents can regionally affect this symmetry), and consider that the excess within the amplitude of the peak at the inertial frequency f 0 is likely to represent the inertial component of the buoy motion.For example, in Fig. 4b, the tidal oscillation generates a small peak at f = 2 cycles day −1 , which means that most of the amplitude of the peak observed at f 0 comes from the inertial oscillation of the buoy.We now consider a buoy tethered to a strongly cohesive multi-year ice pack: the internal stress term of the momentum equation (Eq.2) is then expected to dominate over any external forcing, including the Coriolis force.Consequently, the inertial oscillations that might be generated by storms or moving fronts are immediately damped and thus not visible.Figure 5a shows an example of such a trajectory within the multi-year ice pack at the end of winter.The cycloids seen in Fig. 4a do not appear anymore and the Fourier spectrum plotted in Fig. 5b does not show any peak at the inertial frequency f 0 .Conversely, we can notice that Kwok et al. (2003) found significant oscillations in the same region but in a different time period.This points out the temporal variability in the strength of inertial oscillations, which may depend on both the storm activity and the sea ice state.In the following, a simple method is proposed to quantify the strength of inertial oscillations.Then, performing statistics on this quantity will allow us to interpret space and time variabilities in terms of dissipation of energy through mechanical processes that takes place within sea ice.

Methodology
We now define a parameter M that quantitatively accounts for the time-dependent inertial oscillation magnitude.
The cycloids observed in the trajectories (Sect.3.1) result from the superposition of an advection (f = 0) and a rotation at the inertial frequency (f = f 0 ).In Figs. 3 and 4, the red boxes show parts of the trajectories characterized by cycloidal loops: these loops indicate that the rotation velocity is larger than the advection velocity.In contrast, no full loops are observed within the green boxes, showing that the rotation velocity is lower than the advection velocity.Lepparanta (2004, p. 157) proposed a simple ratio of cycloidal rotation velocity over the advection velocity to qualitatively estimate the strength of inertial oscillations.
We here propose to evaluate in a quantitative way the timedependent oscillation magnitude using the Fourier spectrum at the inertial frequency f 0 .Since inertial oscillations are best quantified by means of a spectral analysis, we here perform such an analysis on a sliding time window.This will enable us to define a time-varying quantity M that measures how inertial oscillations evolve with time as a buoy drifts.For a given buoy location defined by the coordinates (x pst , y pst , t pst ), where "pst" stands for "present", we define a Gaussian window function g pst (t) centered on t pst with a characteristic duration of the order of several inertial time periods: where probe the inertial time scale, since it corresponds to 3 days, and short enough to properly measure rapid time variations in the inertial oscillation magnitude.Thus, 24 data points within which about half the points carry significant weight are used to compute Fourier transform calculations.Gaps (missing data) are present in the IABP buoy records.Their duration varies from one sample every few hours (i.e. a time interval of 6 h between successive positions) to several weeks.Therefore, we introduce in our calculations a selection condition such that no gap of data is allowed for values of t that verify g pst (t) > P , where we set P = 0.01.No M value is thus computed at time t pst if one or more sampling time t verifying this later condition is lacking.Moreover, since missing data even occurring far from t pst (i.e. at values of t such that g pst (t) < P ) can play a non-negligible role on the Fourier transforms, the Gaussian of Eq. ( 5) is truncated at ±1.5 days from t pst .
Under these conditions, we can compute 1.3 × 10 6 M values from the IABP dataset, which represents 83 % of the total number of observations available.
We then compute the velocity W pst (t) as follows: where We compute the normalized Fourier spectrum at time t pst of the velocity time series W pst (t) at the inertial frequency f 0 , where f 0 is made to vary according to the buoy latitude (Eq.1), as where ω 0 = 2πf 0 and t 0 and t end are the starting and ending times of the trajectory.This value is likely to represent the amplitude of the inertial oscillations.We further normalize it in order to define a non-dimensional parameter M that measures the magnitude (still at t pst ) of the inertial oscillation: where is the (current, at time t pst ) mean magnitude of the drift velocity, and the value 1.27 is equal to t end t 0 g pst (t)dt.

Application to the buoy trajectories
Figure 6 shows the values of M computed using the different trajectories plotted in Sect.3.1.Large M values are obtained for the oceanic buoy trajectory plotted in Fig. 3, with an average of 0.56.For the ice drifter of Fig. 4, we also get large M values with an average of 0.8: the M > 1 values are all associated with cycloidal loops observed in Fig. 4a (red rectangles), followed by their decay as characterized by a decrease of M (green rectangles).In contrast, for the ice drifter of Fig. 6c, the M values are much lower: the mean M value over the time period is only equal to 0.13, illustrating the relative absence of inertial loops.
As the raw buoy positions of the IABP ice drifters dataset are irregularly sampled through time, then interpolated (cubic interpolation) and re-sampled at a regular 3 h interval, we investigated the effect of this procedure on the values of M computed for the IABP dataset.To do so, we artificially degraded the oceanic buoy dataset, which samples buoys positions every hour, by randomly removing a given % of raw positions.A cubic interpolation followed by a re-sampling at 3 h is then performed on the degraded dataset, as done for IABP data, and the associated M values (noted M D ) are recomputed (Fig. 7).The cubic interpolation alone, i.e. without degrading the raw data, has almost no influence on M (Fig. 7a).In Fig. 7b, we can see that the deviation of the M D time series from the M time series becomes visible only when more than 50 % of the raw buoy positions are missing.For 50 % of missing data, the error on individual M values is lower than 5 %.More importantly, as we focus in the following on M values averaged over trajectories, the departure of averaged M D values from the average M value for increasing missing data ratios is shown in Fig. 7c.Average M D values do not deviate significantly from the average M value for missing data ratios up to 60 %, which corresponds to an average time sampling of about 02:50 h and a time sampling maximum gap of about 6 h.Thus, for sampling times lower than these threshold values, we can consider that the computed M values are not affected by an irregular sampling through time.While not shown here, similar results have been obtained by considering ice drifters initially regularly sampled at 1 h during the TARA field campaign (Gascard et al., 2008).This indicates the robustness of IABP interpolation procedure coupled to our estimation of the amplitude of inertial oscillations.
These examples demonstrate that the parameter M is appropriate to quantify the inertial oscillation magnitude.Small average values of M (say lower than 0.2) correspond to buoy trajectories within a strongly cohesive ice pack, while large values of M (say greater than 0.6) correspond to buoy trajectories that we can consider to be nearly in free drift condition.

Analysis of 30 yr of IABP data
In this section, we analyse the spatial and temporal patterns of the inertial oscillation magnitude M, and check whether a significant trend can be observed over the Arctic Basin during the 30 yr of the IABP dataset, which would reveal a significant change in ice conditions.As previously explained in Sect.3.1, inertial oscillations are caused by sudden changes  3a, (b) the trajectory of the ice-tethered buoy 1897 plotted in Fig. 4a and (c) the trajectory of the buoy 12825 plotted in Fig. 5a.The red and green rectangles respectively correspond to those shown in Figs. 3 and 4.  in external forces (strong wind), which then decay due to kinetic energy dissipation within the Ekman layer or friction that takes place at the ice/water interface, or due to internal ice stresses.Therefore, more frequent and stronger storms over the years should imply larger M values on average.Cyclonic activity over the Arctic Ocean shows a maximum during summer that could partly explain an annual cycle of the inertial oscillation magnitude (Serreze and Barrett, 2008) (see below).On the other hand, no significant trend in cyclonic activity has been found over the last 50 yr (Serreze and Barrett, 2008).Similarly, the average wind speed over the Arctic Basin, as estimated from the ERA-40 reanalysis dataset, does not show any significant trend over the period 1979-1999 (Rampal et al., 2009).

Seasonal variation
In this section, we compute M for buoys located within the central Arctic Basin, as delimited by the thick black line in Fig. 1, and then investigate intra-annual (monthly) variations by grouping M values into monthly periods, i.e. one average monthly value M for all the M values occurring in January, whatever the calendar year, etc.The central Arctic is here delimited by all buoy positions lying further than 150 km away from the coasts and the Fram Strait.This way, we skip from the analysis buoys possibly stuck on fast ice.
Figure 8 shows the monthly M value averaged over the 30 yr of record.To estimate the associated error bars, we checked using a bootstrap method (Rampal et al., 2009) that the standard deviation M indeed varies as , where N M is the number of M values used to calculate the monthly mean M. Here, N M ranges from a minimum of 6.4 × 10 4 in February to a maximum of 1.1 × 10 5 in May.The monthly mean M value reaches a minimum of 0.168 in May and a maximum of 0.294 in September, and exhibits an obvious annual cycle.We define the summer season by the three months of July, August and September and the winter period by the rest of the year.These results are consistent with other observations: sea ice concentration, sea ice thickness and sea ice deformation also describe an annual cycle (Rampal et al., 2009;Rothrock et al., 2008).The annual cycle of ice thickness within the Arctic Basin has a maximum on 30 April (Rothrock et al., 2008)  , where N M is the number of M values used to calculate the monthly mean M. Here, N M ranges from a minimum of 6.4 × 10 4 in February to a maximum of 1.1 × 10 5 in May.
M value in May.Being thinner and less concentrated in summer, sea ice is less cohesive and more deformable during this time period.This leads, on average, to a lower damping of the inertial oscillations, and therefore larger M values.However, we cannot exclude that a stronger cyclonic activity during summer (Serreze and Barrett, 2008) could reinforce the annual cycle described by the average M values by triggering more oscillations.

Spatial pattern
The question arises as to whether the seasonal changes in the M values evidenced in the previous section can be associated with a spatial pattern.In this section, we select values of M associated to buoys located within the central Arctic Basin and within the Fram Strait (thick black and grey lines in Fig. 1).From theses values, we build a seasonal M dataset by separating the M values computed in winter, i.e. all the M values recorded for all the 30 winter seasons, from the M values computed in summer.The seasonal spatial patterns of M for both seasons are plotted in Fig. 9.The spatially averaged M value, denoted M(X j , Y j ), is computed for each grid point j of coordinates (X j ,Y j ) as where the summation is performed over all the buoy positions (x i , y i ) contained in a circle of radius L = 400 km centered on the coordinates (X j ,Y j ).The weight w ij (x i , y i ) is defined as where d = (x i − X j ) 2 + (y i − Y j ) 2 is the distance buoy position (x i , y i ) and the grid point (X j , Y j ).
The patterns observed in Fig. 9 are, at least qualitatively, consistent with the observed distribution of the sea ice thickness (Kwok and Rothrock, 2009) and concentration (http: //nsidc.org/data/seaice/index.html) (Comiso, 1990(Comiso, , updated 2012)).In summer, large M values are observed at the edge of the basin in the Beaufort, Chukchi and Laptev seas, i.e. in regions where the multi-year ice cover has been progressively disappearing during the last decade (Maslanik et al., 2011).On the contrary, small values of M can be observed along the Canadian coasts, where the average ice thickness is at its maximum (Kwok and Rothrock, 2009).Sea ice behaves more as a strongly cohesive plate in this region, unsensitive to the Coriolis force.This zone corresponds to the multi-year ice still remaining nowadays.Large values of M are also observed south of Fram Strait, which is consistent with the buoy trajectory plotted in Fig. 4a, which we discussed in Sect.3.1.In contrast, the winter pattern does not reveal any particular structure within the basin.The values of M are small over the whole basin, except relatively larger M values computed north of Canadian coasts.These larger M values could be attributed to buoys drifting along the "Circumpolar flaw lead" that consists of a sheared, thus highly fragmented, zone (Lukovich et al., 2011).
To further test the link between the state of the sea ice cover and its cohesiveness, as expressed by the amplitude of the inertial oscillations, we perform a correlation analysis between the M values and the open water concentration 1 − α, where α is the sea ice concentration dataset collected by NSIDC (http://nsidc.org/data/seaice/index.html)(Comiso, 1990(Comiso, , updated 2012)).This dataset has a spatial resolution of 25 km and consists of ice concentration values sampled every two days from 1979 to 1987 (SMMR), and every day since 1987 (SSM/I).For each value of M, we search for the corresponding value of open water concentration as the closest sample in time and space: a 1 − α value is associated with a given M value if we can find, for the same day of record, within one day during the period 1979-1987, a sample that is closer than 25 km.The corresponding correlation coefficient R is equal to 0.245(±0.002).This positive correlation is statistically significant as R is more than 120 times greater than the standard deviation obtained in the null hypothesis at no correlation (numerically computed by randomly reshuffling the M and 1−α values).This is consistent with stronger oscillations characteristic of a less compact sea ice cover.
In order to check whether this global correlation is, or not, only due to the fact that M and 1 − α both describe an annual cycle, we group together these values by months.We then compute the cross-correlation coefficient R separately  10) for each point of a 25 km resolution grid.The graphic representation is a linear interpolation of the gridded M -values.In order to not represent regions with little data, an M value at a given grid point is plotted only if its associated weight is greater than a minimum value we arbitrarily set to 1000.33 Fig. 9. Spatial pattern of inertial oscillation magnitude within the Arctic Basin in summer (left) and in winter (right).These two fields are computed from the seasonal M dataset.A spatially averaged value of M values, denoted M, is computed following Eq.( 10) for each point of a 25 km resolution grid.The graphic representation is a linear interpolation of the gridded M values.In order to not represent regions with little data, an M value at a given grid point is plotted only if its associated weight is greater than a minimum value we arbitrarily set to 1000.
for each month of the year, over the 30 yr of record.While not plotted here, the annual variation of the correlation coefficient is in phase with the annual cycle of M plotted in Fig. 8: the correlation is larger in summer.The correlation coefficient is equal to 0.350(±0.005) in August, which means that the M values are spatially correlated with the associated open water concentration values.In winter, the correlation is much less significant, which can be understood by the fact that, during this period, the sea ice concentrations are equal to 1 almost everywhere within the basin and therefore show very little fluctuations.
Here, all the computed coefficients of correlation remain small.This suggests that other factors than the ice concentration control the M values: we expect both the variability in wind forcing materialized by storm activity and the ice thickness to significantly depress the correlation between M and 1 − α.We will see in Sect.4.3.2 that the correlation between M and 1−α is much improved when setting free of the influence of wind forcing variability, i.e. by considering averaged M and 1 − α values.

Multi-annual variation
We have shown so far that the M parameter is a useful proxy for the cohesiveness of the cover and, consequently, the level of internal stresses.We now analyze whether the M values have changed over the last three decades owing to a change in ice conditions.

M time series
In this section, we use the central Arctic buoy dataset.The temporal sampling of the IABP buoy dataset is heterogeneous, being characterized by variations of data density over the period.Most notably, a larger number of observations are available in the later years, as compared to the early 1980s or during the late 1990s.To circumvent this problem, the evaluation of the multi-annual variation of M is done by equally binning the M values in time.The summer and winter datasets are both separately split into 10 successive bins, all containing the same number of observations.For each bin, an average value of M, associated with an average value of time, is computed.
Figure 10 shows the time series M(t) between January 1979 and December 2008 for the summer and winter seasons.The error in the average M value is computed using a bootstrap method as explained in Rampal et al. (2009).The bootstrap method is performed individually for each bin, and the deviation from the mean M is represented by the error bars.We see on both plots that the most prominent feature is a significant increase of the inertial oscillation amplitude through time.A linear fit gives a positive trend equal to 1.19 (± 0.34) × 10 −5 yr −1 (i.e. a 16.5 % increase per decade) for summer and 5.7 (± 1.9) × 10 −6 yr −1 (i.e. a 11 % increase per decade) for winter.Thus, the increase of inertial oscillations is relatively more marked in summer than in winter.
The use of GPS as the buoy positioning system has become more common since the end of the last century.As noted in Sect.2.1, this system is more accurate than ARGOS positioning.We thus check whether the observed trend on the mean time series M(t) could be a spurious effect caused by reduced noise in recent years.To do so, noised mean time series M N (t) are built by computing M values on buoy trajectories over which Gaussian noise is added.Because the increase in the average M(t) values is particularly marked from the beginning of the 21 century, Gaussian (i.e.16.5 % increase per decade) for summer and 5.7 (± 1.9) × 10 −6 yr −1 (i.e.11 % increase per decade) for winter.
Gaussian distribution with a standard deviation σ .Figure 11 shows realizations of mean time series M N (t) built using σ = 300 m and σ = 1000 m, along with the IABP time series M(t), as in Fig. 10.From the year 2002, the mean time series M N (t) deviate from M(t), showing that noise on the buoy positions has an influence on the parameter M.However, almost systematically, values of M N (t) are larger than M(t), showing that noise most generally acts to increase the M values rather than to decrease them.
This can be explained as follows: noise increases the amplitude of the Fourier spectrum at high frequencies.This has two competing consequences: (i) the norm of the drift velocity W pst slightly increases, which tends to decrease M (see Eq. 9); and (ii) the amplitude at the inertial frequency Ŵpst (f 0 ) is increased, which tends to increase M. The second effect is particularly marked when the inertial motion (so, M) is small, and, most of the time, dominates over the first effect, so that M is increased by noise.Consequently, the evolution observed in Fig. 10 is not the consequence of possible reduced position noise in recent years, as the difference between the periods 1979-2001 and 2002-2008 would have been even larger if the level of noise had remained the same over the entire record.In other words, the observed trend in M(t) is most likely an underestimation of the multi-annual evolution of the inertial oscillation magnitude.The same is also true for the annual cycle (Fig. 7): as the effect of noise is stronger for low M values, winter M values are, in average, more overestimated than summer M values.
In addition, Rampal et al. (2009) reported a significant (9 % per decade) increase of sea ice speeds in the Arctic Basin over the last three decades.Therefore, as increasing advection velocities act to decrease M (Sect.3.2), the multiannual positive trend on the M values reported here most likely underestimates the associated decrease of the mechanical energy dissipation within the ice cover in the later years.

Evolution of spatial patterns
The previous section demonstrated a strong evolution of the average inertial oscillation amplitude, although it is more strongly marked in summer.We here analyze the spatial consequences of the observed multi-annual changes.To do so, we split the whole IABP buoy seasonal dataset into two periods.A Student t-test performed on the winter and the summer M(t)-time series in Fig. 10 shows that the most significant changing point occurs in 5 averaged M values out of the 10.More precisely, splitting the winter time series {M 1 , M 2 , . .., M 10 } into two distributions {M 1 , . .., M 5 } and {M 6 , . .., M 10 }, the probability that these two distributions have the same mean is only 0.92 %.We have verified that this probability is indeed the lowest when splitting the time series just after M 5 .Similarly, we find the same optimal changepoint in summer, with an associated probability that the two distributions share the same mean equal to 0.42 %.This leads to define two distinct periods, with period 1 from 1979 to 2001 and period 2 from 2002 to 2008 (see Fig. 10).We can remark that, according to this periodic splitting, period 1 and period 2 contain the same number of observations, i.e. 5 bins each.As done in Sect.4.2, the spatial patterns of the M values for the summer and winter seasons associated with each period are shown in Fig. 12.As expected, the changes are stronger in summer.They affect most of the Arctic, but are more pronounced on the Siberian side: a drastic increase of the peripheral ice zone area is observed in the later years.Conversely, winter changes are milder north of Canada and Greenland, where the thickest multi-year ice can be found nowadays (Kwok and Rothrock, 2009).
To  10) for each node of a 25 km resolution grid.The smoothing parameter L is equal to 400 km.In order to not represent regions with little data, an M value at a given grid point is plotted only if its associated weight is greater than a minimum value we arbitrarily set to 1000.10) for each node of a 25 km resolution grid.The smoothing parameter L is equal to 400 km.In order to not represent regions with little data, an M value at a given grid point is plotted only if its associated weight is greater than a minimum value we arbitrarily set to 1000.
here, since they correspond to 1 − α values very close to 0 over the whole Arctic Basin.Sea ice concentration values are sampled at buoy positions, as done in section 4.2, and averaged similarly to M values, i.e. following Eq.10.Common features evidenced by M in Fig. 9 (left) and Fig.   10) for each node of a 25 km resolution grid.The smoothing parameter L is equal to 400 km.In order to not represent regions with little data, an open water concentration value at a given grid point is plotted only if its associated weight is greater than a minimum value we arbitrarily set to 1000.
from period 1 to period 2 is also similar to the one observed on M, materialized by a slight migration of the pack zone toward the south, and an increase of the peripheral zone area, materialized by larger open water concentration values.The correlation coefficient computed between M values shown in Fig. 12 and values of 1 − α is equal to 0.57 for period 1 and 0.80 for period 2. These correlation coefficients are much larger than the one obtained in section 4.2.This is explained by the fact that average M and 1 − α values are considered here prior to compute the correlation coefficient, removing like this the influence of wind forcing activity in the variability of M. Thus, we expect the discrepancies that still remain here between M and 1 − α to be mostly due to changes in ice thickness and ice degree of fragmentation.The respective role of ice concentration, ice thickness and ice degree of fragmentation on the inertial oscillation amplitude is thoroughly examined in Gimbert et al. (2012).

Is the observed evolution from the IABP buoy dataset representative of the whole Arctic Basin?
In the two previous sections, care was taken to remove the effect of temporal heterogeneities in the IABP buoy dataset.However, the IABP buoy dataset is spatially heterogeneous: changes in the spatial sampling by the buoys could affect the overall trend.In order to check whether the results reported in Fig. 10 truly indicate a significant increase of the M values, we investigate whether this trend could be an artefact of changes in spatial sampling.Following the procedure described in Rampal et al. (2009) for ice velocities, we formulate a null hypothesis: in this hypothesis, there exists no temporal changes in M-maps over the year, but the effect at spatially sampling these maps in different ways at different periods will cause the M-time series to evolve with time.To construct these M-maps, we consider any buoy position associated with a given M value and recorded at a given year and a given season.For such a position, we compute the mean M 0 value for summer, respectively winter, by considering all the M values recorded in summer, respectively winter, and whatever the year, contained within a circle of radius L = 200 km, using Eq. ( 10).This mean value is therefore year-independent, and corresponds to our null hypothesis of no inter-annual changes.Then, we calculate our null-hypothesis M 0 time series as the mean of the summer, respectively winter M values, using the same bins as in Fig. 10.Any inter-annual variation observed in the time series M 0 (t) could only be explained by changes in the spatial sampling rather than by an actual, global trend: for example, a positive trend would be explained if IABP buoy had a tendency to sample regions associated with large M values more often in the later years as compared to earlier years.
Figure 14 shows the mean time series M(t) and M 0 (t) over the whole time period 1979-2008.The trend associated with the null hypothesis is equal to 1.57 (± 0.98) × 10 −6 yr −1 in summer and 1.28 (± 0.51) × 10 −6 yr −1 in winter.This means that changes in spatial sampling are responsible for approximately 10 % of the observed trend in summer and 20 % of the observed trend in winter.We thus conclude that the trend observed in the average M values cannot only be explained by an irregular sampling of the buoys.This increase in M is not an artefact and reveals a genuine increase of the inertial oscillations magnitude over the whole Arctic Basin.

Conclusions
From the IABP ice drifter trajectories recorded between 1979 and 2008, we analyzed the average amplitude of inertial oscillations over the Arctic sea ice cover.To do so, we defined a non-dimensional parameter M that quantifies the magnitude of inertial oscillations relative to advection motion along the trajectories.For a sea ice cover of low concentration, constituted from a loose floe field moving nearly in free drift, inertial motions and M values are large.On the reverse, a highly cohesive sea ice pack, characterized by strong internal stresses, is expected to be associated with low M values.
From appropriate averaging of this 30-yr dataset at different space or time scales, we have shown the following: (i) The M values describe an annual cycle with a minimum reached in May and a maximum in September, in a qualitative agreement with the corresponding annual cycles of sea ice extent (Comiso et al., 2008), concentration, thickness (Rothrock et al., 2008;Kwok and Rothrock, 2009), advection velocity and deformation rates (Rampal et al., 2009).(iii) A significant increase of average values of M is observed from 1979 to 2008.This increase, although more marked in summer, is observed in both seasons and is associated with the reduction of the thick, multi-year ice zone in recent years.
From the expected link between the magnitude of the oscillations and the degree of consolidation of the ice cover (McPhee, 1978;Colony and Thorndike, 1980;Geiger and Perovich, 2008), we believe that point (iii) is a signature of the mechanical weakening of the Arctic sea ice cover in recent years.However, one may argue from the momentum balance of sea ice (see Eq. ( 2)), that this strengthening of inertial motion might simply and more directly result from the observed thinning of the cover, i.e. a reduction of ice mass per unit area.In addition, this evolution could, to some extent, be the result of a modification of vertical penetration of turbulent momentum within the ocean boundary layer.In Gimbert et al. (2012), we have shown, from a simple ocean-sea ice coupled dynamical model, that these two explanations cannot fully account for the evolution of inertial motion in the Arctic, which actually reveals a genuine mechanical weakening, through an associate decrease of the sea ice internal friction magnitude, of the cover at the basin scale.Such mechanical weakening has (will have) strong consequences in terms of ice drifting speeds, deformation rates, export (Rampal et al., 2009) and therefore on mass balance (Rampal et al., 2011).

Fig. 1 . 26 Fig. 1 .
Fig. 1.Map of the Arctic basin showing the buoy trajectories of the IABP dataset.The positions are sampled every 3 h from January 1979 to December 2008 and plotted following a stereographic projection centered on the North Pole.The Laptev sea is poorly covered by this dataset.The thick solid black and grey lines define the Central Arctic basin and the Fram strait region, respectively.figure Fig. 3. (a) 1 month sample of the trajectory of buoy B245 (25 May to 25 June 2007).The trajectory is plotted in latitude-longitude coordinates.Its beginning is marked by the red circle.The red and green boxes outline periods with strong and low cycloidal loop activity, respectively.(b) Fourier spectrum of the buoy velocity.The velocity is computed following Eq.(3) with t = 1 h.

Fig. 4 .
Fig. 4. (a) 1 month sample of the trajectory of IABP buoy 1897 (30 August to 30 September 1987).The location of the buoy is indicated by the red box in the inset of Fig. 4a.The red boxes show periods when cycloidal loops are best observed, whereas the amplitude of the loops is much lower in the green boxes.(b) Fourier spectrum of the buoy velocity.The velocity is computed following Eq.(3) with t = 3 h.

Fig. 5 .
Fig. 5. (a) 1 month sample of the trajectory of buoy 12825 (20 April to 21 May 1992).The location of the buoy within the Arctic Basin is indicated by the red box in the inset of Fig. 5a.(b) Fourier spectrum of the buoy velocity.The velocity is computed following Eq.(3) with t = 3 h.

Fig. 6 .
Fig.6.M values for (a) the trajectory of the oceanic buoy 255 plotted in Fig.3a, (b) the trajectory of the ice-tethered buoy 1897 plotted in Fig.4a and (c) the trajectory of the buoy 12825 plotted in Fig.5a.The red and green rectangles respectively correspond to those shown in Figs.3 and 4.

Fig. 7 .
Fig. 7. (a) M values for the trajectory of the oceanic buoy 255 plotted in Fig. 3 (raw data: blue line, same as Fig. 6a) and the trajectory obtained after a cubic interpolation and a 3 h resampling of the raw positions (red line).(b) M D time series computed from the degraded datasets of missing data ratio varying between 0 and 60 % (the red line is the same as in a).(c) M D values averaged over the whole trajectory (1 month) as a function of the missing data ratio.20 realisations have been done at each given value of missing data ratio.The grey dashed line shows the average raw M value, equal to 0.56.

Fig. 8 .
Fig. 8. Average annual cycle in the M values.The average is performed considering the central Arctic Basin buoy dataset.The error bars show the deviation M from the average M value estimated using the central limit theorem and computed as M = M √ N M

Fig. 9 .
Fig.9.Spatial pattern of inertial oscillation magnitude within the Arctic basin in summer (left) and in winter (right).These two fields are computed from the seasonal M dataset.A spatially averaged value of M values, denoted M , is computed following Eq.(10) for each point of a 25 km resolution grid.The graphic representation is a linear interpolation of the gridded M -values.In order to not represent regions with little data, an M value at a given grid point is plotted only if its associated weight is greater than a minimum value we arbitrarily set to 1000.

Fig. 10 .
Fig. 10.Time series M(t) from January 1979 to December 2008 computed using (a) the summer dataset and (b) the winter dataset.For each season, the dataset is equally split into 10 bins over which the average is computed.One bin corresponds to approximately 3.1 × 10 4 M values in summer and 8.0 × 10 4 M values in winter.The horizontal lines associated with each data point indicate the time period associated with each bin.Vertical lines are error bars computed from a bootstrap method.Bold lines are linear fits: the trends are 1.19 (± 0.34) × 10 −5 yr −1(i.e.16.5 % increase per decade) for summer and 5.7 (± 1.9) × 10 −6 yr −1 (i.e.11 % increase per decade) for winter.

Fig. 11 .Fig. 12 .
Fig. 11.Influence of noise, prescribed on the positions of buoys, on the monthly mean time series M(t) for (a) the summer season and (b) the winter season.Gaussian noise, with a constant standard deviation σ , has been added to all buoy positions recorded since year 2002.Results for noised mean time series M N (t) (in yellow and green) for σ = 300 m and σ = 1000 m are compared with the IABP mean time series M(t) as in Fig. 10. 36

Fig. 12 .
Fig. 12. Spatial repartition of M within the Arctic Basin in summer and winter computed using the seasonal dataset of the period 1979-2001 and the seasonal dataset of the period 2002-2008.An average mean value of M, denoted M, is computed following Eq.(10) for each node of a 25 km resolution grid.The smoothing parameter L is equal to 400 km.In order to not represent regions with little data, an M value at a given grid point is plotted only if its associated weight is greater than a minimum value we arbitrarily set to 1000. 12

Fig. 13 .
Fig. 13.Spatial repartition of open water concentration within the Arctic Basin in summer computed using the seasonal dataset of the period 1979-2001 and the seasonal dataset of the period 2002-2008.An average mean value of open water concentration, denoted 1 − α, is computed following Eq.(10) for each node of a 25 km resolution grid.The smoothing parameter L is equal to 400 km.In order to not represent regions with little data, an open water concentration value at a given grid point is plotted only if its associated weight is greater than a minimum value we arbitrarily set to 1000.

Fig. 14 .
Fig. 14.Mean time series M(t) (a) in summer (red) and (b) in winter (blue), as in Fig. 10 and M 0 (t) (black lines), resulting from the null hypothesis.The best linear fits are also plotted, showing that changes in spatial sampling only account for a small fraction of the observed trend.

(
ii) The spatial pattern of M over the Arctic Basin is in agreement with the sea ice thickness and concentration patterns.Low M values are observed in western Arctic, whereas large values are observed within a peripheral zone (Beaufort Sea, eastern Arctic) and south of Fram Strait.