A quasi-objective single buoy approach for understanding Lagrangian coherent structures and sea ice dynamics

. Sea ice drift and deformation, namely sea ice dynamics, play a signiﬁcant role in atmosphere-ice-ocean coupling. Deformation patterns in sea ice can be observed over a wide range of spatial and temporal scales, though high resolution objective quantiﬁcation of these features remains difﬁcult. In an effort to better understand local deformation of sea ice, we adapt the Trajectory Stretching Exponents (TSEs), quasi-objective measures of Lagrangian stretching in continuous media, to sea ice buoy data, and develop a temporal analysis of TSE time series. Our work expands on previous ocean current studies 5 that have shown TSEs provide an approximation of Lagrangian coherent structure diagnostics when only sparse trajectory data is available. As TSEs do not require multiple buoys, we ﬁnd they have an expanded range of use when compared with traditional Eulerian buoy-array deformation metrics, and provide local-stretching information below the length-scales possible when averaging over buoy-arrays. We verify the ability of TSEs to temporally and spatially identify dynamic features for three different sea ice datasets. The ability of TSEs to quantify trajectory stretching is veriﬁed by concurrent ice fracture in buoy 10 neighborhoods ranging from tens to hundreds of kilometers in diameter, as well as the temporal concurrence of signiﬁcant storm events.


Introduction
The Arctic is warming at a much greater rate than the rest of the Earth and sea ice plays an important role in regulating energy exchanges between the atmosphere, cryosphere, and ocean. The recent decline in sea ice extent and the prevalence of younger, 15 thinner ice is well documented (Rothrock et al., 1999;Kwok and Rothrock, 2009;Landy et al., 2022) with thin ice being more susceptible to deformation and fracture. As the ice warms in spring, melt is accelerated around existing fractures due to a reduction in albedo and the presence of more open water. Arctic amplification, the disproportionate warming of the Arctic in a changing global climate, has been partially attributed to the enhanced oceanic heating and ice-albedo feedback caused by diminishing sea ice (Screen and Simmonds, 2010;Dai et al., 2019;Thackeray and Hall, 2019;Jenkins and Dai, 2021). Changes 20 in the Arctic sea ice cover in recent years have also led to modifications of larger atmospheric circulation patterns (Moore et al., 2018), and mid-latitude weather (Siew et al., 2020). It is therefore of great importance to accurately measure sea ice dynamics to understand the state of ice deformation, fracture, and refreezing, as well as the impacts on global climate, infrastructure in Arctic waters, and access for northern communities.
Sea ice deformation is a physical phenomenon that is highly localized in space and time (Oikkonen et al., 2017;Rampal 25 et al., 2019), with local and regional deformation influencing summer sea ice melt rates and future weather. Even short winter storms have been shown to have long-term impacts on sea ice extent and melt (Graham et al., 2019b;Lukovich et al., 2021).
Quantitative methods that identify deformation rates and sea ice flow patterns are thus critical for understanding atmosphereice-ocean exchange processes and grounding sea ice models.
Previous studies have used both displacement grids and sea ice buoy trajectories to quantify observations of sea ice defor-30 mation (e.g., Rampal et al., 2009;Hutchings et al., 2011;Szanyi et al., 2016a;Oikkonen et al., 2017). Over several decades, various approaches have been developed, including both Eulerian diagnostics that focus on ice behavior in a single time slice and Lagrangian diagnostics that assess the evolution of ice flow and the underlying flow patterns (Szanyi et al., 2016a, b).
Perhaps the most common sea ice deformation metrics are the instantaneous values divergence, shear, and total deformation as defined from components of the two-dimensional velocity gradient, ∇v (Leppäranta, 2011): (1) These Eulerian diagnostics (1-3) have been used extensively for both sea ice model validation, and experimental observations.
An alternative approach to studying sea ice dynamics is through methods that robustly identify significant sea ice flow 40 features that persist over time, i.e. Lagrangian coherent structures (LCSs). LCSs are rigorously defined mathematical structures adapted from dynamical systems theory that behave as the underlying skeleton of fluid flow, thus shaping patterns in the evolution of material (Haller, 2015). Of particular relevance to the study of sea ice deformation are hyperbolic LCS. As a material evolves over a particular time window, hyperbolic LCS identify the curves of maximum stretching and compression.
Previous work by Szanyi et al. (2016a, b) studied hyperbolic LCS to determine distinct zones of deformation in Arctic sea ice, 45 and how they varied over time, but in general the extent of their use in sea ice dynamics is currently limited.
In Figure 1 we show an example of a hyperbolic LCS retrieved from ocean surface current data. This visualization reveals why hyperbolic LCS are particularly relevant for studying the deformation of sea ice. In panel 1a, an initial gray square of fluid at time t 0 is shown. The hyperbolic repelling (M R ) and attracting (M A ) manifolds are drawn in blue and red, respectively.
As the fluid evolves from time t 0 to t 2 (panels a-c), fluid particles are attracted to M A with material stretching along that axis, 50 and repel from M R , with compression along that direction. When evaluating sea ice dynamics from an LCS perspective, it follows that separations between distinct dynamic regions may correspond with strong shear zones or separation features, such as cracks and ridges, as previously performed at a pan-Arctic scale by Szanyi et al. (2016a). It is our hypothesis that resolving the time-varying dynamics of these coherent structures may also help identify periods of significant ice modification, such as during influential storms or other breakup events.

55
Hyperbolic LCS are typically approximated by ridges of the finite time Lyapunov exponent (FTLE) (e.g., Szanyi et al., 2016a). Calculation of FTLE fields, however, rely on advecting grids of particles in time-resolved velocity fields and quantifying spatial gradients of the flow map, as discussed in detail in Appendix A1 and by Haller (2015). Obtaining gridded sea ice velocity data that can be used to calculate rate-of-strain invariants (1-3) and FTLE fields is a significant hurdle. This is currently only feasible at large scales via motion tracking algorithms and remotely sensed data. Several of these gridded velocity 60 products are also known to be hindered by artificial velocity discontinuities caused when assimilating different data streams (Bouillon and Rampal, 2015;Szanyi et al., 2016b). These limitations prevent calculation of reliable Eulerian and Lagrangian diagnostics, especially at high (sub-daily) temporal resolution. In contrast, buoy GPS trajectories are often sampled at hourly or sub-hourly timescales and can provide a wealth of true sea ice motion data. Lagrangian and Eulerian diagnostics from these data sources can help address some of the limitations in gridded velocity studies.
at statistics for all suitable triangles (e.g. no small angles) in an array . For an N -sided polygon, with vertices indexed in counterclockwise orientation such that (x N +1 , y N +1 ) = (x 1 , y 1 ), ∇v can be approximated as follows: where A is the area of the polygon.
Complications with the Green's theorem method include sensitivities to user choice of array, accounting for GPS signal-to-80 noise ratios, and several physical inaccuracies arising from discretization. As arrays of buoys passively follow ice floes, they typically do not maintain a shape that allows accurate calculation of the spatial velocity gradients from drift trajectories. Some of these issues have been previously discussed by Lindsay and Stern (2003); Hutchings et al. (2012) and Dierking et al. (2020), including the polygon array selection method. Because of the necessary discretization of a continuous boundary, integral approximations by finite sums in the Green's theorem approach also generate errors. Lindsay and Stern (2003) previously 85 addressed boundary representation errors for nonlinear velocity fields, but only errors from discontinuities (leads) that intersect polygon boundaries were studied. Uncertainty also stems from the trapezoid rule approximation of the contour integral (4).
For a known continuous velocity field (e.g. without cracks), the upper bound of integral approximations can be quantified (Atkinson, 1989), though an explicit impact of this error cannot be found in the sea ice literature. As is detailed in Appendix A2, the trapezoid rule error can cause equilateral triads to indicate both divergent and convergent conditions in a steady, 90 divergence-free, and continuous (fracture-free) flow, depending on array orientation.
In light of these complications, alternative methods that do not rely on buoy geometry or orientation, such as single trajectory metrics, could be largely beneficial, even if only to identify and separate dynamically active regions. One primary issue with single-trajectory metrics is that they are not frame-indifferent. That is, they depend on the reference frame of the observer. One may argue that it is sufficient to require all computations to be performed in the same reference frame, but one of the main 95 axioms of mechanics is that the material response (e.g. ridging, fracturing) of a material is independent of the observer (Gurtin, 1981). One intuitive interpretation of this axiom is a fracture formed in ice is independent of whether the ice is viewed from an airplane, or from the ice surface. Frame-indifference is, therefore, a foundational benchmark that diagnostics must meet in order to identify coherent features in any deforming continuum. Invariants of the rate-of-strain tensor are frame-indifferent, as is FTLE, but velocity and velocity gradients are not.

100
As an alternative to Eulerian contour integral diagnostics, we propose using novel single-buoy stretching diagnostics, the trajectory stretching exponents (TSEs) (Haller et al., 2021 correspond to proximity to hyperbolic LCS, which are regions of significant attraction, repulsion, and shear. We can then use single buoys to find spatial and temporal domains when sea ice is behaving like it is near an underlying hyperbolic LCS. There are limited frame-indifferent Lagrangian alternatives for coherent structure identification with sparse buoy data to 110 compare to. One notable exception is single-point (squared) relative dispersion (Haller and Yuan, 2000), also known as twoparticle dispersion, which monitors the separation between a pair of particles. Similar variations have been applied in various contexts by Rampal et al. (2009), and Lukovich et al. (2014, 2015, 2017. These approaches are not evaluated here as one cannot generate both spatial maps and time series in a manner comparable to TSEs or (1-3). This is because relative dispersion relies on initially nearby buoys whose motions decorrelate at longer time scales. Relative dispersion diagnostics may show 115 spatially coherent structures for a short window of time, but as the initially close buoys become decorrelated, one must reselect new buoy pairs, thus preventing a single meaningful time-series of local conditions as is possible with TSEs or (1-3). In this way, methods based on relative dispersion are more often used for understanding scaling relationships of sea ice motion.
In the present research, we show that TSEs localize regions and periods of significant stretching and compression, some of which is not evident with conventional buoy metrics. We independently verify the TSE stretching and compression features 120 with two approaches. While TSEs are not a fracture-specific diagnostic, concurrent lead formation in remotely-sensed data provides a frame-indifferent confirmation of material sea ice response. We also compare our event detections with synoptic storm analysis from a high-resolution sea deformation experiment. This approach provides coupling between TSEs and the passage of storms previously verified to have significant sea ice influence. We discuss several logistical, and computational advantages with the single-buoy diagnostics, as well as future insights possible with TSE applications. We also compare our 125 quasi-objective Lagrangian ice deformation analysis with standard Eulerian array-based rate-of-strain metrics.

Analytical Methods
For a continuously differentiable velocity field v(x, t), a particle's trajectory x(t) is governed by the ordinary differential equationẋ(t) = v(x(t), t), where t represents time. Consider a material curve γ(t; s) ⊂ U ⊂ R 2 , parameterized by the scalar parameter s ∈ R at time t that has evolved from an initial curve γ(t 0 ; s). One can quantify the stretching of vectors tangent 130 to this material curve using the equation of variations. In a steady flow (with no time dependence), particle trajectories are themselves material curves, and the Lagrangian velocity vector v(t) evolves as v(t) = ∇v(x(t))v.
(5) Haller et al. (2021) showed that the average material stretching over the time interval t ∈ [t 0 , t N ] of the Lagrangian velocity vector, v 0 = v(x 0 ), can be written as For our situation, equation (6) gives an objective measure of sea ice stretching from only a single buoy velocity v(x(t)).
This degree of stretching also only depends on initial and final conditions. That is, interim cycles of stretching and compression can cancel out. In sea ice, cycles of stretching and compression may lead to significant ridging and fracturing that one would not see in a standard fluid flow. To quantify the cumulative impact of repeated stretching, relaxation, or compression over a 140 given time window, we can also utilize the averaged hyperbolicity strength with a strictly positive integrand, Equation 7 adds up all hyperbolic (stretching and compression) action, giving a measure of cumulative changes, but does not differentiate stretching or compression. For slowly-varying (steady) flows, λ t N t0 approximates (equals) the finite-time Lyapunov exponent associated with the initial vectorẋ 0 at the initial position x 0 , an objective Lagrangian measure of deformation over 145 time (Ott and Yorke, 2008;Haller et al., 2021).
Remark 1: In an unsteady flow, the right hand side of Eq (5) would also include a δ t v(x(t), t) term. For unsteady flows, the magnitude of the approximation error of tangential stretching in Eq (6) strongly depends on the slowly-varying nature of the flow, such as when Eulerian acceleration is much smaller than Lagrangian acceleration. The pointwise slowly-varying nature of sea ice velocities are assessed in Appendix A3 using Polar Pathfinder Daily 25 km EASE-Grid Sea Ice Motion Vectors (Tschudi 150 et al., 2018), similar to the analysis conducted for ocean currents by Haller et al. (2022). Our assessment of the slowly-varying assumption does not depend on a specific integration time as it is a point by point comparison of Eulerian and Lagrangian acceleration. It does, however, depend on the temporal and spatial resolution of the underlying Pathfinder velocity field. A finer resolution velocity product could reveal sharper velocity gradients at smaller scales, such as across fractures, and further support the slowly-varying hypothesis. A shorter sampling period may also reveal stronger temporal gradients surrounding 155 fracture events. However, this cannot be investigated until higher resolution and artifact-free velocity fields become available, and is a topic of future research.
Remark 2: Tangential stretching is rigorously defined when the trajectory is traveling in a continuous, incompressible medium. While it is true that sea ice velocity fields are not continuously differentiable, we can assume they are piecewise continuously differentiable (between fractures), and thus trajectory stretching exponents are well-defined in large continuously deforming ice covers. In a steady incompressible two-dimensional flow, a negative stretching exponent, λ t N t0 , equals the growth exponent of a vector normal to the material line formed by the trajectory x(t) (Figure 2) (Haller et al., 2021). For unsteady compressible flows, such as sea ice, λ t N t0 andλ t N t0 are proxies for the material change in a plane normal to the trajectory, and may indicate stretching, ridging, or fracture. Negative λ t N t0 values indicate compressive processes whereas positive values indicate stretching along the trajectory ( Figure 2). The degree to which this differs from the growth exponent of the plane or-165 thogonal to the trajectory depends on slowly-varying and compressibility conditions, as well as the rheology of the ice. These considerations parallel those for standard rate-of-strain metrics when dealing with a fractured ice cover, varying rheology, or regions with sea ice concentration less than 100%. For example, one is no longer solely quantifying the divergence of sea ice material when the ice is fractured and leads are present, rather the divergence of the mixed water-ice continuum. With these considerations in mind, we focus primarily on mid-winter and early spring ice dynamics to minimize extensive fracturing of 170 the ice cover, but briefly discuss dynamics during ice disintegration in Section 4.2. We further suggest due consideration when calculating stretching from ice buoys as ice concentration decreases and their trajectories become more representative of ocean dynamics than ice cover dynamics.
For real discrete observational trajectory data with an initial position x 0 = x(0), we approximate the integrals in eqs (6) and (7) and define the trajectory stretching exponents (TSEs) used for the remainder of this research: As with λ, TSE is positive for stretching, and negative for compression along a trajectory. TSE does not allow for this cancellation as the summand is strictly positive, and gives a measure of cumulative changes, but does not differentiate stretching or compression, as with λ. TSE and TSE have been shown to accurately separate dynamically distinct regions (eddies and fronts) 180 in sparsely-sampled open ocean flows (Haller et al., 2021), but have not yet been studied in a sea ice context.
In contrast to eqs. (2-3), TSEs do not differentiate between contributions from divergence or shear to the stretching of a material. Rather, they quantify stretching in the direction of vectors tangent to the buoy trajectory. We discuss the mathematical relationship between (2-3), the closely-related Lagrangian averaged divergence (Szanyi et al., 2016a), and TSEs in Appendix 1. TSEs are calculated using only buoy speed and do not require projection to orthogonal velocity components as in Green's 185 theorem approximations from arrays. Speed can be easily calculated using geodesics between GPS locations, which prevents any inconsistencies of results due to map projections. Furthermore, TSEs are parameter-free with integration time being the only user-chosen value. The choice of integration time defines the average stretching extent, but does not prevent identification of events at timescales shorter or longer than that. See Section 4.1 for one such example of hourly timescale events being identified by a multi-day integration window.

190
The sampling period may influence the estimated speed of each buoy (Lei et al., 2021) and should be mentioned when performing analysis with TSEs. Once raw data has been collected, the user may then resample the data at a higher or lower resolution. Changing this sampling period is akin to changing the resolution of a Riemann sum that is approximating an integral.
That is precisely what eqs. (8 & 9) are with respect to eqs. (6 & 7). The choice of interpolation method would have analogous effects.

195
While Eulerian strain-based metrics (2-3) provide an imperfect comparison to a Lagrangian Coherent Structure-based approach to ice deformation, they are commonly used in sea ice dynamics research, and generate long time series when polygons are suitably chosen (e.g., Hutchings et al., 2012;Itkin et al., 2017). TSEs are Lagrangian diagnostics and are calculated over a forward-looking window from a starting time, t 0 (eqs. 8 & 9). By incrementally increasing t 0 , we can also generate time series of TSE and TSE for each step in a buoy trajectory. In this way, large TSE and TSE time series values should slightly precede 200 significant storms, sea ice stretching, or breakup events as the summation occurs from t 0 forward in time. This also means that TSE and TSE series will be slightly out of phase with concurrent Eulerian polygon-based divergence, shear, and deformation series. Indeed, we find this phase difference is on the order of the length of the integration window when maximizing the cross-correlation between TSEs and div and D (Eq. 1-4, not pictured). We refer the reader to Appendix A1 for a thorough numerical comparison of Eulerian and Lagrangian diagnostics in a simple geophysical model.

205
Results here focus on the ability of TSEs to accurately characterize spatially and temporally localized coherent stretching during short-lived mid-winter storms. This choice is motivated by the availability of recent high-resolution experiments and local deformation being an open problem in sea-ice dynamics (Oikkonen et al., 2017). As with other Lagrangian methods, the integration time is user-defined and reflects the duration over which we are seeking significant and coherent stretching. For example, previous LCS analysis using 6-month FTLE fields identified three distinct dynamic regions at pan-Arctic scales: the Fram strait, Beaufort gyre, and Northwind Ridge (Szanyi et al., 2016b). Similarly, TSEs are applicable to both short or long integration times (t N − t 0 ), and reveal time-averaged dynamics at those specified scales. We refer the reader to the Appendix A1 for a comparison of TSEs when varying integration times.
The choice of TSE integration window is only limited by the accuracy of the GPS trajectory data. Shorter integration times likely affected more significantly by small discrepancies in GPS locations. The instantaneous limit of the TSEs are also 215 well-defined and provide approximations of hyperbolic (or parabolic) objective Eulerian coherent structures, which are the instantaneous limits of hyperbolic (or parabolic) LCSs (Serra and Haller, 2016;Nolan et al., 2020).

Datasets
Our results verify that TSEs can identify significant ice dynamics events using three different data sources at different spatial scales: synoptic storm analysis from a targeted experiment, linear kinematic feature (LKF) formation in subsequent SAR frames 220 on the order of tens of kilometers, and sea ice brightness temperature measurements that show landfast and free-ice fractures at scales of hundreds of kilometers. We exhibit the utility of TSE and TSE spatio-temporal analysis for two experimental sea ice buoy data sets, the N-ICE2015 (Itkin et al., 2015 and MOSAiC expeditions (Krumpen et al., 2020), and data from the International Arctic Buoy Program (IABP) (International Arctic Buoy Programme, 2022). sampling frequency using a linear interpolant, and we follow this convention for our N-ICE2015 analysis. After interpolation, subsequent buoy speeds that exceeded 5 km/day were removed and buoy positions were resampled using a linear interpolant. Itkin et al. (2017) calculated divergence, shear, and deformation using a tessellation of buoy triads, removing some of the user-dependence, but potentially introducing triads with inappropriate geometry for the contour integral approach. In particular, the winter deployment of buoys for the N-ICE2015 campaign was hindered by logistical challenges as the researchers deployed buoys in the polar night by snow machines and on skis. This resulted in an initially quasi-linear buoy array geometry, with many

245
The Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC), was the largest multidisciplinary Arctic expedition to date, spanning the winter of 2019-2020 (Krumpen et al., 2020). The study centered around the research icebreaker Polarstern (Knust, 2017;Nicolaus et al., 2022;Rabe et al., 2022) which was moored to an ice floe for an entire year. MOSAiC provides an unprecedented look at winter ice dynamics with 213 unique buoy deployments and complementary atmosphere, ocean, ecology and biogeochemistry data. 250 We focus here on the paths of 101 buoys deployed within 40 km of the Polarstern. This public data set was documented by (Bliss et al., 2022). The half-hourly buoy track data was cleaned following (Hutchings et al., 2012). Triads were also handpicked from the MOSAiC buoys with data spanning October 2019 to June 2020, and is the focus of a forthcoming publication. The arrays were selected to maintain reasonable shapes (no small angles, area greater than 1km 2 ) from the beginning to the end of the time series and the public data was resampled to uniform 6-hourly intervals. Handpicking triads, however, does require 255 user discretion. Buoy tracks were resampled to match the triad sampling rate. The arrays used are shown in Figure 4. A deeper comparison and refinement of geometrically suitable arrays in the MOSAiC data is a current topic of research. The method we use here is in line with previous work (Hutchings et al., 2011(Hutchings et al., , 2012.  Though MOSAiC analysis is now just beginning, there is currently limited mid-winter buoy data available for the Arctic.
MOSAiC and N-ICE2015 provide the best testbeds for evaluating the efficacy of TSEs to identify significant stretching events 260 at high spatial resolution in coherent ice floes.

IABP
The international arctic buoy program (IABP, https://iabp.apl.uw.edu/data.html) is a network of drifting buoys deployed in the Arctic ocean by a collection of nations since the 1970's. In comparison with the previous two experiments, IABP buoy density is significantly lower, but the length of the IABP buoy record and its considerable spatial coverage allows for analysis 265 of ice dynamics at much greater scales (see Figure 5). IABP data have been a useful ground-truth of sea ice motion for several decades, and continue to provide invaluable information for both model verification and understanding the complexities of ice dynamics (Rampal et al., 2009;Bouchat and Tremblay, 2020). IABP buoy tracks are available in both raw and cleaned formats.
We focus our analysis on the publicly available 3-hourly L2 IABP buoy data. We downsampled the data a 6-hour sampling to maintain consistency with the MOSAiC data using a linear interpolant. We focus our analysis on trajectories during 2017 270 ( Figure 5), with particular interest in the well-studied Beaufort sea.

N-ICE2015
For our first evaluation of TSEs, we utilize previous storm analysis from the N-ICE2015 expedition. Cohen et al. (2017) studied winter storms during during this sea ice cruise, and found they typically occurred with a duration between 3 and 10 days. Using 275 a buoy-array approach (Eq. 1-4), Itkin et al. (2017) found that sea ice deformation events could be roughly coupled with the passage of these storms. As we wish to quantify coherent ice deformation over a given time window, we choose a three-day integration period for stretching exponent calculations so as to capture the storm-scale Lagrangian stretching without averaging over too wide of a time window. In this way, TSE provides a before-and-after stretching measure, while TSE calculates the accumulation of trajectory stretching and contraction during the three-day window.

MOSAiC
Our second example utilizes MOSAiC buoy data from November 2019 to April 2020. We again choose a 3-day integration window in the absence of an alternative integration scale. Should an exhaustive storm study, such as that conducted for N-ICE2015, be performed, revisiting this analysis may be insightful. While the effect of integration scale does have some effect We plot all TSE time series available (101) and again find there is impressive coherence in stretching exhibited by the entire cluster. Local values of TSE and TSE also have the largest variance during peak stretching events. This allows us to analyze both the temporal localization of sea ice fracture by TSE, as well as potentially relate spatial differences of trajectory stretching to nearby LKFs. In Figure 7b and 7d, blue dots indicate the instantaneous value from each buoy triad, and the black lines show 310 the mean over all triangles.
Qualitatively, prior to the January 31 event there is a reasonable agreement between TSEs and triad-based events, when accounting for the three-day integration window. For the entire period, one finds that all strong divergence and deformation events are accounted for as extrema (> 1d −1 ) in the TSEs, though often with a different relative magnitude. In mid-November 2019, large TSE magnitudes generate a spike in TSE that is also evident as spikes in div and D. Local TSE maxima in mid  and late December also correlate well with spikes in div and D. Mean div provides a suggestion of whether divergence or convergence was more prevalent amongst the triads at a given time, but the significant scatter of both positive and negative values obscures such behavior at finer spatial scales. In contrast, there is much less ambiguity of the characteristic ice behavior for the entire cluster in the TSE series suggesting the cluster remained part of a largely coherent structure, that stretched and compressed together, with variance between buoys appearing during dynamic events.

320
After January 31, several relatively large peaks in TSE appear as minor oscillations in div or D, or are not present. This indicates local dynamic stretching and potential fracture events may not be evident with the contour integral metrics. One such example is shown for the peak on April 17 highlighted by the vertical black dashed line. In div and D plots, we also show the 3-day integration window in green. This example was chosen due to the concurrent availability of pre-and post-storm SAR data during the period after the major January 31 event, as well as a relatively minor signature in div and D. The large  In the 3-day window following the Apr 17 TSE and TSE peak, the mean buoy divergence oscillated around zero ( Figure   7b), with the magnitude staying below 0.1d −1 . This is approximately 1% of peak values of mean divergence, suggesting a relatively insignificant period of divergence. This is in contrast to TSE and TSE on April 17 which sits at approximately 50% of their total peak values, suggesting a relatively localized motion with a larger contribution to ice dynamics at the same time.
The contribution to total deformation from shear was slightly larger, with a peak mean deformation of 0.46d −1 at 17:00, April 17. The choice of triangles dictates which LKF are spanned during this dynamic event, and may explain why there was a weak deformation signature from the triads. As well, the lack of objectivity or trapezoid rule errors when using Green's theorem approximations with triads, as detailed in Appendix A2, may also contribute. In this scenario, the stretching and relaxation measured by TSE presents a clear correlation with material deformation of the ice and suggests TSEs may provide ice behavior insight during times when Green's theorem methods are not possible, such as when there are too few buoys or 340 they are their orientation is inappropriate, and when array-based approaches have underestimated dynamic behavior.
As sea ice melts in spring, ice concentration drops and the ice pack loses integrity, TSEs and rate-of-strain diagnostics are no longer strictly representing the deformation of the ice cover. In Figure 8, we look at mean TSE and mean D during different stages of the melt period identified by (Lei et al., 2022). In late spring TSE values increase as the ice becomes more dynamic.
Once the ice floe enters warmer upper ocean water over the Yermak plateau, basal melting increases and TSE oscillations trend 345 upwards. At the same time, oscillations in D also begin to increase. By June 25, snow and snow-ice had completely melted from the sea ice cover, further enhancing the surface melt rates. After this, we obtain peak TSE values for the year, followed by a sharp drop. It is at the end of this sharp decline in stretching that Krumpen et al. (2021) note the ice floe begins to disintegrate.
We hypothesize that TSE spikes precisely because this onset of disintegration enters the integration window.
At the timing of this drop in TSE, and onset in disintegration, D values also obtain large values never seen before in the

355
In this section, we focus specifically on behavior in the Beaufort sea from October 2016 to October 2017. During this winter, thinner than usual sea ice may have caused the collapse of a typical high-pressure system over the Beaufort and anomalous surface winds (Moore et al., 2018). We again choose a 3-day integration window for our study to maintain comparability with the previous examples, though other timescales are equally applicable.
In Figure 9a, we show the three-day TSE for each buoy in the Beaufort during that period in grey. The maximum distance

Conclusions
We find that TSEs successfully identify significant local material deformation tangent to individual sea ice buoy trajectories.
Periods of strong local stretching are representative of changing ice dynamics in a neighborhood to 10-1000 km. That is, though TSEs are local in nature, their values provide insight into much larger coherent sea ice structures, and distant fracture events.

395
This local-regional connection provides a valuable avenue for researching sea ice dynamics in sparsely sampled regions, and understanding the changing rheology of sea ice. In contrast, conventional polygon-based approaches provide an estimate of the spatial-average of divergence, shear, or deformation over an ice volume, with errors increasing as area decreases, number of buoys decreases or the array becomes skewed.
Approaching sea ice dynamics through quasi-objective stretching, we were able to capture coherent deformation events 400 in concentrated and sparsely-sampled buoy experiments. Spatial and temporal signatures of stretching with TSE are well correlated with formation of nearby leads and changing ice transport patterns, as well as with the concurrence of well studied sea-ice-impacting storms. For two high resolution mid-winter buoy deployments, we find that the TSE had greater sensitivity to sea ice deformation than the common polygonal approach, and potentially avoided a false positive identification.
Specifically, large TSE values coincided with major storms in the N-ICE2015 experiment, and did not identify a mysterious 405 storm-free dynamic ice-event that was described by Itkin et al. (2017). For the first half of the MOSAiC experiment, we found a good qualitative agreement between polygon-based metrics and the TSEs, though this correlation was much weaker following a large midwinter storm. During this latter half, we verified TSEs were able to identify an influential stretching event that was much less evident with Green's theorem methods. Lastly, TSEs were able to spatially and temporally isolate major fracture events during Beaufort sea ice breakup using IABP data at a much lower spatial resolution. A buildup of stress in the Beaufort 410 sea was detailed by increasing local TSE maxima, until an initial separation formed at the land-fast ice during the seasonal TSE peak. During subsequent diminishing peaks the ice fractured and further detached from land-fast regions, creating a coherent mass of rotating ice that was delimited spatially with TSE. Increasing amplitude TSE events preceded spring breakup events in both high-resolution MOSAiC and sparsely-sampled IABP data. This suggests predictive abilities of TSE may possibly be developed with further investigation.

415
The single-buoy quasi-objective trajectory stretching exponents (TSEs) identify dynamic sea ice events that are potentially significant in terms of understanding spatially and temporally varying sea ice deformation. As sea ice dynamics plays an important role in atmosphere-ice-ocean exchange processes, we find the further event-detection sensitivities possible with TSEs are a valuable complement to common, polygon-based divergence, shear, and deformation approximations. We find TSE usage provides some distinct advantages and potential improvements for future research: 420 1) TSEs identify deformation events from a single buoy and are thus unaffected by buoy array geometry, or the number of buoys deployed. As buoys passively follow ice floes, we can still obtain insights into ice dynamics as arrays become heavily skewed and non-uniform. TSE maintains integrity in this situation, whereas arrays of buoys become aligned and unsuitable for strain estimates due to shearing in the ice pack. This expands the domain of potential analysis, and reduces logistical burdens of Arctic and Antarctic expeditions as stretching dynamics can be studied with TSE from linear, randomly distributed, or 425 organized arrays. This is particularly valuable for the Antarctic where there are significantly fewer buoys deployed than in the Arctic.
2) As TSEs are quasi-objective metrics in continua, their values provide a useful proxy for frame-indifferent (objective) measures of stretching, with a quantitative difference depending on the compressibility of the ice and slowly-varying nature of the flow. In a heavily fractured ice context, the degree to which these single-trajectory metrics approximates along-trajectory 430 ice stretching varies, but this work suggests positive correlations with remotely-sensed dynamic fracture and breakup events.
Further numerical and observational experiments can help improve these correlations.
3) TSE values do not need to be separated or averaged based on length scales (e.g. Itkin et al., 2017), and showed great success in identifying fracture regions at a wide range of spatial scales.
4) TSE calculations are mathematically simple as TSE is calculated using only buoy speed and does not require projection 435 to orthogonal velocity components. Speed can be calculated using geodesics between GPS locations, which prevents any inconsistencies of results due to choice of map projection or projection distortions, or needing to perform differential calculus in an elliptic geometry. Furthermore, TSE is parameter-free with integration time being the only user-chosen value.
Calculating TSE from single buoy trajectories requires significantly less effort than SAR-based approaches, is not subject to the same errors (e.g. Bouchat and Tremblay, 2020), and also supplements data in the pole hole during winter (e.g. Krumpen 440 et al., 2021). If we can further verify the slowly-varying nature of sea ice at sub-diurnal time scales, TSE may also fill in temporal gaps due to the sparse monitoring of satellites. Calculating TSE fields from more general feature trajectories, such as from X-Band and HF radar or SAR datasets, would be a straight-forward application of Eq (8) & (9) and can also enrich the analysis of sea ice dynamics in existing datasets.
The ability of TSEs to overcome some of the short-comings of other buoy approaches may also provide an additional source 445 of deformation information, such as is necessary to constrain and improve sea-ice models (Bouchat and Tremblay, 2020). To obtain rate-of-strain invariants for sea ice deformation, it is still necessary to use a high-density buoy array. Such an array also reveals gradients of trajectory stretching and further enhances precise stretching localization with TSEs. A modeling study of compressibility, slowly-varying impacts, and sampling frequency would be useful for TSE applications in model development, as would a further comparison with divergence, deformation, and shear to aid the broader community in the 450 physical interpretation of TSE signals. This would help interpret TSEs as complementary sources of stretching information when buoy arrays are aligned linearly due to strong shear, have edges that span large distances encompassing multiple fractures or coherent structures, or have areas below tolerable error thresholds. In remote polar regions where data is still difficult to obtain, but the changing climate has an outsized impact, TSEs provide new physics-based insights into ice dynamics while only requiring single GPS tracks.

455
Our findings suggest that spatio-temporal analysis of local and regional sea ice dynamics should be revisited with a comparison between single-trajectory and spatial velocity gradients. Our comparisons suggest potentially significant dynamic events in terms of understanding atmosphere-ice-ocean exchange processes may go undetected using conventional polygon-based methods, or in the case of N-ICE data, may be inaccurately detected. This advantage stems from the TSEs' ability to identify stretching in a small neighborhood of individual buoys and connect local changes with broader ice behavior. A deeper under-460 standing of TSEs connection to more broadly used ice dynamics metrics will help researchers understand how TSE can inform ice responses to various forcings.

Appendix A: Introduction
The first section of the Appendix develops the rigorous connection of trajectory stretching exponents to more common deformation diagnostics through the rate-of-strain tensor and Cauchy-Green strain tensor. The second section of the Appendix 465 provides a simple analytic example of buoy motion to illuminate Green's theorem estimation errors for a buoy-array diagnostic. This example shows how the Green's theorem approximation fails to provide self-consistent divergence quantification (in a divergence-free flow), depending on equilateral triad orientation. The last section of the Appendix provides an assessment of the assumption that sea ice is a slowly varying flow to satisfy criteria for the quasi-objectivity of TSEs.
We can express the common Eulerian sea ice deformation metrics, divergence and shear, as functions of the eigenvalues 475 (principal stresses) of (1): To quantify the deformation of a material over time, one considers the flow map that takes an initial position x 0 to its current position x(t; t 0 , x 0 ) following the trajectory defined by the different equatioṅ From the gradient of the flow map, we define the right Cauchy-Green strain tensor as Invariants of C are commonly used to study the deformation of material in continuum mechanics (Truesdell and Noll, 2004).

485
In an incompressible two-dimensional flow, the rate of length change for an infinitesimal material element vector l based at the point x is Likewise, consider a curve γ parameterized by the dummy-time variable s that is evolving under the flow map, F t t0 (x 0 ). This could represent the evolution of a line of dye in a fluid, or a physical transect along sea ice. A vector ξ 0 tangent to γ at x 0 490 evolves as ξ t = ∇F t t0 ξ 0 . It then follows that In the case that l or ξ 0 are eigenvectors of their respective tensors, (A7) and (A8) equal the respective eigenvalues λ i , which represent the stretching of an infinitesimally small sphere into an ellipse along its principal axes Haller (2015). When ξ 0 is the eigenvector associated with the the larger Cauchy-Green eigenvalue, λ 2 , one easily obtains the widely-used finite-time

495
Lyapunov exponent, Otherwise, one can define the averaged stretching exponent also known as the finite-time Lyapunov exponent associated with the initial vector ξ 0 . Due to the exponential growth of fluid particle separation and chaotic nature of fluid flows, one often utilizes stretching exponents in this context instead of eigenvalues 500 of deformation tensors, though they are clearly closely related.
Lagrangian averaged versions of shear and divergence can also be computed as respectively. LAD was previously used by Szanyi et al. (2016a) to differentiate between different contribution to FTLE at pan-Arctic scales.

505
To calculate the rate-of-strain tensor, once must calculate spatial derivatives of the velocity field. Similarly, the right Cauchy-Green strain tensor requires the spatially and temporally resolved velocity field information to accurately calculate the gradient of the flow map. This sort of spatially and temporally resolved velocity data is unavailable from sparse buoy trajectories. TSEs are designed to complement deformation metrics in sparse data settings.
Consider a trajectory x(t). One can calculate the Lagrangian velocity vector along the trajectory 510 v(t) = v(x(t)) =ẋ(t), but we do not have any information about nearby velocities. Haller et al. (2021) show that, in a reference frame where the velocity field is steady (does not change with time), x(t) is a material curve, and v(t) evolves as a tangent vector to the trajectory. Thus the stretching of v(t) is a measure of material stretching along the trajectory and we can closely approximate the averaged stretching exponent in that frame without the need of spatial derivatives. That is, in a steady flow (no time dependence), TSE measures the correct stretching of Lagrangian velocity vectors as they materially evolve in that reference 515 frame, as determined by the gradient of the flow map (A8).
The instantaneous limit of TSE corresponds with the analogous Eulerian rate of length change (A7) for the choice of vector l = v(t 0 ). In the following examples we highlight the ability of TSE and TSE to delineate distinct coherent flow features without the use of spatial derivatives and compare their behavior with other diagnostics. These models are idealized flows without many of the complication from sea ice motion. That is, they represent steady flow incompressible flow fields of a 520 continuous medium and are thus ideal for comparing these diagnostics.

A1.2 Examples
For our first example, we compare Eulerian and Lagrangian diagnostics for the widely studied model of geophysical fluid flow (see, e.g., Rypina et al., 2007), the steady Bickley flow. The stream function of the flow is given by where we use the geophysically based parameters U 0 = 62.66e−6, L = 1.77, c = U 0 /2, A = 0.1, r 0 = 6.371, k = 6/r 0 . As the velocity field is derived from the stream function, the velocity field is divergence-free. The flow consists of a central westward jet, two eastward jets above and below, and clockwise vortices. Figure A1 shows streamlines and velocity vectors for this flow, as well as divergence, shear, and vorticity fields calculated on a 500×150 grid. Note that the div = 0 by design, and that the vortex cores align when comparing streamlines and vorticity contours. Shear (Fig A1c), shows a shear-free core, and high shear 530 regions adjacent to the jet and vorticies, though contours of shear do not agree with the topology suggested by the streamlines or vorticity. This indicates that shear is not a good indicator of the edges of all coherent structures. Figure A2 details TSE and TSE fields for a 500×150 grid of initial conditions, calculated from a trajectory integration time of 90 days. We also calculate LS (Fig. A2c) for each trajectory x(s). For brevity we omit LAD as the flow is divergence free.
Contours of TSE or TSE organize particles with similar trajectory stretching and reveal coherent structures of the flow.

535
Note that both TSE or TSE can identify the edges of the central jet and eddies as shown in the streamlines of this steady flow. TSE also reveals where cumulative trajectory-tangent stretching is greatest, near the edges of the eddies. TSE provides intricate stretching and compression information surrounding the vortex cores and hyperbolic saddles, suggesting a complex stretch-relax cycle for material in the flow. LS also does a much better job than Eulerian shear to extract the boundaries of each feature, and also highlights the stronger shear near the edges of the central jet.
These underlying structures, and their subsequent influence on sea ice deformation is exactly the kind of information that TSEs provide for sea ice dynamics analysis. As can be seen, mathematically there is limited relationship between the Eulerian metrics (eqs. A2-A4), and the stretching exponents (A9), particularly because ξ is not necessarily aligned with the eigenvectors of S. Qualitatively, one can see there is no relation between div and TSEs for this simple example, and very limited similarity with shear, except when evaluating a Lagrangian average (LS). 545 We also examine the influence of integration time on TSEs for this steady flow. In Figure A3, we increase the integration time of particle trajectories from 1 day to 30 days and plot TSE and TSE contours. The eddies, saddles between eddies, and central jet are evident for all integration times, but the clarity of boundaries increases as integration time increases. Note this improvement is not always possible with time-dependent flows as features may change. In more complex turbulent flows, such clarity at higher integration times can also be hindered by trajectories traveling near multiple uncorrelated coherent structures.

550
For our second example we use there is in fact no reliable relationship to be inferred between TSEs and shear with a simple divergent flow It is easy to analytically calculate that shear and vorticity are constantly zero, and divergence is uniformly equal to 2. This is 555 reflected in the contour plots of Figure A4. Stretching is also constant along trajectories in this flow due to the linear acceleration away from the origin. In this way, TSE and TSE reveal the same uniform deformation structures as the divergence field ( Figure   A5). As shear is everywhere zero, LS is also zero, and reveals the same uniform structure as the TSEs, divergence, and LAD.

A1.3 Conclusion
The previous mathematical exposition and examples show there is no simple relationship between TSE and TSE with div, 560 shr, or vorticity. This is in part because TSE and TSE are Lagrangian metrics, they measure material deformation over time, whereas div, shr, and vorticity are Eulerian, relying only on instantaneous rates of change at one point in time. TSE and TSE do not approximate div, shr, or vorticity.
As such, TSE and TSE fields are most successfully used for identifying dynamic flow features, such as the edge of a gyre, or deformation across a jet. This stems from an integration of temporal changes in material stretching along a trajectory. As value of div G is sensitive to the triad's rotation on the unit circle. Figure S2a shows the range of div G calculated as the triad is rotated. Depending on how the vertices are oriented, this approach will suggest either a divergent or convergent flow, with the correct value only appearing for two specific orientations.
In contrast to errors that have been discussed and quantified before (e.g., Lindsay and Stern, 2003;Hutchings et al., 2012;Dierking et al., 2020), this issue arises from the approximation of the partial derivative contour integrals by the trapezoid rule.

580
For a generic real-valued twice differentiable function f (x) defined on the interval [a, b] ⊂ R, the difference between the true integral and trapezoid rule is bounded by (Atkinson, 1989): That is, the potential error of each partial derivative scales with the cube of the distance between buoys and the second derivative of the associated velocity component along that interval. If we increase the number of uniformly distributed buoys on the unit 585 circle in this velocity field, we find the average divergence approximation inside the polygons quickly converges to the true value ( Figure S2b). In fact, using a four sided polygon can significantly reduce the error for this simple flow, as an equal number of buoys can be positioned on either side of the axis dominant flow structure ( Figure S2c). This improvement can be expected for other flows as well, though more complex flow structures would require a greater number of vertices to avoid these same errors. Analogous errors are evident in Green's theorem approximations of shear and total deformation as well, 590 reminding researchers to consider the number of buoys and methods used when quantifying sea ice dynamics.

A3 Slowly-varying condition of sea ice
TSE is a quasi-objective metric, indicating it approximates an objective metric under a given condition. For the work presented here, this requires that the velocity field of interest be slowly varying. That is, |v t (x(t), t)| |a(t)|, where a(t) is the Lagrangian velocity. In order to calculate the temporal derivative on the LHS, we rely on daily gridded sea ice velocity data from 595 Tschudi et al. (2018) as a best estimate for the flow. While one can calculate a(t) from trajectory data, discrete calculations of v t (x(t), t) rely on repeat measurements of velocity at the same location. In this way, we require gridded velocity data to validate the slowly-varying assumption.
The use of Pathfinder sea ice displacement grids comes with its own errors, such as the spatial gradient artifacts previously identified by Szanyi et al. (2016a). Furthermore, the low temporal resolution of this data smooths the short-time variability of 600 surrounding instantaneous shear and fracture events. Prior to the development of a more robust gridded dataset of ice velocities, this remains the best available option.
The probability function in figure A7 shows that for the majority of the Arctic domain, the slowly varying condition is satisfied pointwise along sea ice trajectories. This probability is representative of mid-winter ice conditions, and may vary in the summer months. As this is a pointwise ratio, and not an integrated comparison, the slowly-varying condition does not 605 depend on a specific time window, the resolution of the underlying Pathfinder velocity field. Video supplement. The supplemental videos provide a NASA Worldview animation, and TSE buoy animation of the Beaufort Sea ice breakup detailed in Section 3.3. The animations can be viewed at https://youtu.be/WHCsOaL4Nks and https://youtu.be/l2tOJSnTfSY .   to Eulerian shear measurements. Lagrangian averaged divergence is omitted as it is everywhere zero. Figure A3. Comparison of TSE and TSE contours for increase integration times from 1 to 30 days (L to R) for a subsection of the Bickley jet.
The eddies, saddles between eddies, and central jet are evident for all integration times, but the clarity of boundaries increases as integration time increases. Note this improvement is not always possible with time-dependent flows as features may change. Figure A4. a) Streamlines and velocity vectors for the divergent flow (A13). b) The flow is designed to be uniformly divergent divergence value of 2. c) Instantaneous values of shear no shear. d) Vorticity fields are also uniformly zero. Figure A5. a-b) TSE and TSE show uniform stretching due to the uniform linear acceleration away from the origin. This matches the structure revealed by the uniform divergence field. c) Lagangrian averaged divergence is also constant, thus presenting the same pattern as TSE and TSE. LS is omitted as it is constantly zero in this shear-free flow.