1. Introduction
Sea ice dynamics, including drift and deformation, are crucial in the polar climate system, significantly influencing mass and heat exchanges between the atmosphere, ice and ocean. These processes directly contribute to Arctic sea ice loss through enhanced ice outflow to the North Atlantic Ocean. This outflow accounts for about 10% of the total Arctic sea ice area (Kwok, Reference Kwok2009), and is associated with the Atlantic Meridional Overturning Circulation (Serreze and others, Reference Serreze2006).
Sea ice kinematics is primarily driven by momentum transfer from wind and currents (Leppäranta, Reference Leppäranta2011), modified by ice conditions such as concentration, thickness and consolidation (Heil and others, Reference Heil2008, Reference Heil, Massom, Allison and Worby2011; Oikkonen and others, Reference Oikkonen, Haapala, Lensu, Karvonen and Itkin2017). The decline in sea ice extent, coupled with the prevalence of younger and thinner ice during the recent decades in the Arctic Ocean (Kacimi and Kwok, Reference Kacimi and Kwok2022), has made ice cover increasingly susceptible to deformation and fracture during drifting. Intensified surface wind stress associated with cyclone events can also dramatically alter the ice kinematic regime, promoting fragmentation of ice floes, especially in the marginal ice zone (MIZ; Itkin and others, Reference Itkin2017). Measurements by buoy arrays deployed in the central Arctic Ocean indicate that the temporal evolution of sea ice motion and deformation was closely correlated with the cyclonic strength and tracks. Sudden changes in wind forcing during extreme synoptic events can enhance the inertial component of ice motion (Lei and others, Reference Lei, Gui, Heil, Hutchings and Ding2020, Reference Lei2021), prolonging the activities of ice kinematics after the cyclone passes (Haller and others, Reference Haller, Bruemmer and Mueller2014). Ice inertial motion is damped by the friction at the ice-ocean interface and internal ice stress, particularly in compacted and consolidated ice fields. Thus, the strength of ice inertial oscillation can, in turn, serve as an index to measure the consolidation of the ice field (Lammert and others, Reference Lammert, Bruemmer and Kaleschke2009). In Arctic peripheral seas with shallow waters, tidal forces significantly influence ice kinematics, allowing ice to oscillate (Artana and others, Reference Artana, Provost, Ferrari, Bricaud, Poli and Park2022) in similar frequency bands to the inertial oscillation. Thus, analysing the ice kinematics at specific timescales using the complex Fourier transform or wavelet analysis is essential for quantifying the respective contributions of forces to the ice motion (e.g. Haller and others, Reference Haller, Bruemmer and Mueller2014).
This study focuses on the Arctic Transpolar Drift (TPD), the region associated with the Arctic sea ice outflow (Babb and others, Reference Babb2023). The sea ice conditions in the TPD region exhibit significant spatial heterogeneity (Krumpen and others, Reference Krumpen2019). Sea ice at the downstream end of the TPD has the widest boundary with the open water, likely involves some shallow water, and has experienced more rapid sea ice retreat in winter than other Arctic sectors over the past decades (Gerland and others, Reference Gerland2023). Therefore, the spatial heterogeneity of sea ice dynamics from the pack ice zone (PIZ) in the north to the MIZ in the south across the TPD is expected to be more pronounced than in other Arctic regions. Furthermore, sea ice in this region is increasingly susceptible to cyclones because of both the thinning of sea ice and increased cyclone activity from the North Atlantic Ocean (e.g. Zhang and others, Reference Zhang2023). The interactions between the atmosphere, sea ice and ocean in this region are incredibly active and complex, especially during extreme cyclones and/or in the MIZ. During the Norwegian Young Sea ICE (N-ICE 2015) expedition in the downstream region of the TPD, younger and thinner Arctic sea ice (Itkin and others, Reference Itkin2017) and winter storms (Graham and others, Reference Graham2019) resulted in unexpectedly enhanced ice deformation, significantly affecting the local sea ice mass balance through creating new ice over the leads (Itkin and others, Reference Itkin2018). However, the lack of synchronous observations makes it unclear what the differences are in sea ice kinematics and deformation in the downstream TPD affected by cyclonic activities compared to those in the highly consolidated ice pack in the central Arctic Ocean beyond the cyclonic influence.
As Arctic sea ice decreases, there is a significant northward shift of the MIZ (Strong and others, Reference Strong, Foster, Cherkaev, Eisenman and Golden2017). Insufficient observations also limit our understanding of the dynamic processes affecting sea ice in the MIZ, particularly regarding the application of free-drift theory (Brunette and others, Reference Brunette, Tremblay and Newton2022). Additionally, two MIZ definitions, based on the sea ice concentration (SIC) and SIC variability thresholds (Vichi, Reference Vichi2022), reveal distinctions, especially in the seasons with high ice mobility. Thus, there is still no clear consensus on defining MIZ from the perspective of sea ice dynamics.
The traditional approach to determining ice deformation intensity is to calculate the strain rate tensor from a triangular network formed by buoy arrays. The deformation rates calculated by this method are sensitive to the length scale and geometry of the buoy array. However, constructing a sustained and uniformly distributed buoy array is challenging in the TPD region because of sea ice drift patterns and strong ice dynamics. As an alternative option, trajectory-stretching exponents (TSEs) estimated from the trajectory of a single buoy can be used to characterise ice deformation during drifting (Haller and others, Reference Haller, Bruemmer and Mueller2014, Reference Haller, Aksamit and Encinas-Bartos2021). Compared with traditional deformation metrics during N-ICE 2015 and MOSAiC drift, local TSE maxima correlated well with spikes in total deformation rate by full spatial derivatives, and had greater sensitivity to sea ice deformation than the common polygonal approach (Aksamit and others, Reference Aksamit, Scharien, Hutchings and Lukovich2023). This method allows comparison of ice deformation rates derived from sparse buoys, capturing the heterogeneity at large scales. Moreover, Lagrangian observations often have mixed spatial variations with seasonal changes in sea ice kinematics, particularly for ice drifting along the TPD from high-latitude PIZ over deep waters to low-latitude MIZ partially over shallow shelves. This limitation emphasises the importance of synchronous observations in regions with different sea ice conditions over the TPD.
On 21 January 2022, an Arctic cyclone with record low sea level pressure (SLP) formed in East Greenland and tracked northward over the Barents and Kara seas, resulting in the record weekly reduction of regional sea ice area in the Barents Sea in 1979–2022 (Blanchard-Wrigglesworth and others, Reference Blanchard-Wrigglesworth, Webster, Boisvert, Parker and Horvat2022). Coincidentally, an array of ice-drifting buoys was deployed in the Arctic TPD region in late summer 2021 by the Chinese National Arctic Research Expedition (CHINARE-2021). As these buoys drifted, their trajectories spanned regions from the PIZ in the central Arctic Ocean to the ice edge north of Svalbard, capturing the impact of this extreme cyclone in the MIZ during late January 2022. The observations of buoys close to the centre of this cyclone provided an excellent opportunity to clarify the response mechanisms of sea ice kinematics to cyclones over the MIZ. Combining the data measured by the buoys deployed in other regions of the Arctic Ocean during the same period, we can further characterise the regional differences in ice kinematics with various ice conditions and degrees of cyclonic impact.
Thereby, this study primarily aims to address the following questions using data collected from the ice season of 2021–2022 as a case study:
1. What are the differences in ice kinematics between the MIZ and PIZ in the TPD region during the melting and freezing seasons?
2. How do strong storms impact sea ice kinematics in the downstream TPD?
3. What changes occur in the sea ice kinematics under intensified atmospheric and oceanic forcing in the MIZ of the downstream TPD?
Addressing these questions is particularly important for further understanding the spatiotemporal changes and influencing factors of sea ice kinematics in the Arctic TPD region.
2. Data and methods
2.1. Buoy observations
During CHINARE-2021, we deployed nine ice-tethered drift buoys (Ice Tracker, Pacific Gyre IT, USA) and one Snow and Ice Mass Balance Array (SIMBA, the Scottish Association for Marine Science Research Services Ltd, Scotland) in the TPD region in August 2021. In addition, data from 14 buoys deployed over the Arctic Ocean in the same year, obtained from the International Arctic Buoy Program (IABP), were used as a supplement and comparison. We used the position data from these 24 buoys (Fig. 1) to describe temporal and spatial variations in ice kinematics between August 2021 and February 2022, the transition period from melting to freezing.
(a) Buoy operation periods, with red line representing the main array, and orange and blue lines denoting two buoy arrays of CA and CB. (b) Buoy drift trajectories are overlaid on the ocean bathymetry (available from the International Bathymetric Chart of the Arctic Ocean, Version 4.0), with a coloured line denoting the drift date. The thin red and blue lines represent the sea ice edges (SIC = 15%) on 1 August 2021 and 28 February 2022, respectively. Coloured dots are buoys in the main (red), CA (orange) and CB (blue) arrays.
According to the initial deployment area, this study divided these 24 buoys into three groups. Twelve buoys, including 10 buoys deployed by CHINARE-2021 and 2 IABP buoys, constituted the main buoy array. The main array was deployed on ice floes in the downstream part of the TPD close to the MIZ north of the Kara Sea during summer 2021. They drifted in the TPD, through the transition to ice growth and involvement in the PIZ, and finally into the MIZ again, close to the ice edge north of Svalbard. The trajectory of the main buoy array was mainly located in the Atlantic sector of the Arctic Ocean, 80°–87°N and 30°–140°E, which is a corridor for both Arctic sea ice outflow and northward tracks of cyclones. By the end of February 2022, two-thirds of the main buoys drifted into the MIZ in the Fram Strait or north of Svalbard.
The other 12 IABP buoys were divided into two groups of CA and CB, representing a seasonal-ice dominated zone in the upstream TPD and a multiyear ice (MYI) dominated zone in the central Arctic Ocean, respectively (e.g. Babb and others, Reference Babb2023). The CA array of seven buoys was used to compare the difference between sea ice kinematics at the TPD upstream and downstream regions. The data collected by five CB buoys can be used to create a benchmark representing the ice kinematics in the PIZ region throughout the entire study period. Comparison against this benchmark allows us to more clearly identify the impacts of the changes in ice conditions and extreme storms on the kinematics of ice floes with the main buoy array.
2.2. Other auxiliary data
Meteorological Data, including SLP, 2-m air temperature and 10-m wind speeds in 1979–2022, obtained from the ERA5 global reanalysis of the European Centre for Medium-Range Weather Forecasts (ECMWF; Hersbach and others, Reference Hersbach2023), was used to construct a climatology and identify anomalies of atmospheric forcing in the study period along the buoy trajectories. The daily products of sea ice motion and SIC derived from the Nimbus-7 Scanning Multichannel Microwave Radiometer (SMMR) and its successors (SSM/I and SSMIS), provided by the National Snow and Ice Data Center (NSIDC; Tschudi and others, Reference Tschudi, Meier and Stewart2019), were used to estimate anomalies of ice speed and concentration. The merged CryoSat-2-SMOS sea ice thickness product from 2010 to 2022, provided by the European Space Agency (2023), was used to depict ice thickness anomalies.
The lead area fraction along the trajectories in the study year was obtained from the ice lead product provided by Hoffman and others (Reference Hoffman, Ackerman, Liu, Key and McConnell2022). This lead fraction was used to analyse the impacts of sea ice advection and deformation on the formation and evolution of leads. Ocean bathymetry data with a resolution of 200 m, taken from the International Bathymetric Chart of the Arctic Ocean (IBCAO) Version 4.0 (Jakobsson and others, Reference Jakobsson2020), was used to determine the bathymetry along the buoy trajectories.
2.3. Storm events and cyclone tracking
Storm events and their impacts are determined primarily by periods of strong wind speed and secondarily by decreases in SLP (Cohen and others, Reference Cohen, Hudson, Walden, Graham and Granskog2017). This study defines storm events as periods with 24-hour moving average wind speeds exceeding 8 m · s−1 and lasting longer than 6 h, without breaks longer than 6 h. Among the storms, those with the maximum SLP decrease rate exceeding 1.0 hPa per hour are categorised as ‘Severe’ (Se) storms and others are classified as ‘Moderate’ (Mo) storms, as shown in Table 1.
Summary of the moderate (Mo) and severe (Se) storm events involving the main buoy array
Arctic storms were mainly caused by the intrusions of cyclones from the south. We derived cyclone tracks from hourly ERA5 SLP fields using the Lagrangian detection and tracking algorithm following Crawford and others (Reference Crawford, Schreiber, Sommer, Serreze, Stroeve and Barber2021). As cyclone impacts on sea ice dynamics are generally confined within a 400 km radius from the low-pressure centre (Peng, Reference Peng2019), we mainly focus on cyclones with the centre situated less than 400 km from the main buoy array. There were three events (Se1, Se2 and Se3) that have been identified as severe storms, with the minimal central SLP lower than 980 hPa, approaching the main buoy array on 1 November, 3 December 2021 and 21 January 2022, respectively. In particular, the Se3 storm was associated with the record-breaking Arctic cyclone, which had the 1979–2022 record low central SLP of 935 hPa occurring on 21 January 2022 (Blanchard-Wrigglesworth and others, Reference Blanchard-Wrigglesworth, Webster, Boisvert, Parker and Horvat2022). This cyclone formed in the Greenland Sea and gradually approached the main buoy array as it deepened. During this event, a maximum wind speed of 24.1 m · s−1 was observed at the buoy locations. At that time, the cyclone centre was about 120 km from the geometric centre of the main buoy array. This cyclone also caused a northward retreat of the sea ice edge by about 50 km in the northern Barents Sea due to the enhanced northward wind crossing the ice edge, and a decrease of average SIC at the main buoy array from 100% to less than 80% within one week from 21 to 27 January 2022. That is to say, this synoptic process rapidly changed the sea ice regime from PIZ to MIZ at the main buoy array, following the MIZ definition by Strong and Rigor (Reference Strong and Rigor2013). It likely caused the failure of one buoy on 28 January 2022 at 78.9°N, 1.0°W, about 10.3 km from the ice edge. Note that most winter extratropical cyclones propagating towards the Arctic Ocean have the potential to promote the increase in ice concentration because the potential ice divergence can increase ice freezing in the newly formed leads (Schreiber and Serreze, Reference Schreiber and Serreze2020). The pronounced decreases in ice concentrations along the main buoy trajectories observed by this study were mainly because (1) the southward advection of ice floes towards the ice edge would reduce the surrounding ice concentration around the trajectories, and (2) the warm seawater close to the ice edge would prevent the leads from refreezing.
By comparing the kinematic behaviours of sea ice between the main buoy array and those of the other arrays during or after the storm events, we can further reveal the influence of storms on the ice motion and deformation.
2.4. Analysis of sea ice kinematics and deformation
Before calculating ice drift speed, position data measured by the buoys at either 0.5- or 1-h intervals were interpolated to an hourly regular interval (τ). Generally, Lagrangian dispersion statistics were used to quantify dynamical regimes in the ice cover (e.g. Lukovich, Reference Lukovich, Hutchings and Barber2015). For an ice buoy in an ensemble buoy array, the absolute dispersion is defined as
\begin{equation}{A^2} = {\left| {{x_i}\left( t \right) - {x_i}\left( 0 \right) - \lt {x_i}\left( t \right) - {x_i}\left( 0 \right)} \gt \right|^2},\end{equation} where x i is the zonal or meridional position of the ith buoy in the array, as a function of the elapsed time 
$t$. The angle brackets denote ensemble averaging of the array. The absolute dispersion of buoy’s trajectories has strong scale effects and depends on the properties of the ice field. That is to say, the flow dynamics of ice field can be characterised by the scaling exponent 
$\alpha $:
\begin{equation}
A^2 \sim t^\alpha
\end{equation} where 
$\alpha $ < 1 corresponding to a sub-diffusive dynamical regime which captures oscillatory behaviour in the Lagrangian velocity field, 
$\alpha $ = 1 to a diffusive regime characteristic of random motion, and 
$\alpha $ > 1 to a super-diffusive regime including subcategories where 
$\alpha $ = 2 for a ballistic regime, 
$\alpha $ = 5/3 for an elliptic regime and 
$\alpha $ = 5/4 for a hyperbolic regime. In this study, the scaling exponent 
$\alpha $ was utilized to divide the study period into various stages according to the difference in ice kinematic regime.
Furthermore, we use multiple metrics to comprehensively analyse the spatiotemporal changes in sea ice kinematic characteristics and the associated influencing factors.
On the temporal scales of a day or longer, wind is the primary force for the balance of sea ice momentum. We assumed that there is an approximately linear relationship between ice velocity (U) and 10-m wind (
${U_{10}}$), following Thorndike and Colony (Reference Thorndike and Colony1982):
\begin{equation}U = \,k{e^{ - i\theta }}{U_{10}} + \varepsilon ,\end{equation} where 
$U = u + iv$ is the ice drift vector, 
$k$ is the wind factor, 
$\theta $ is the turning angle between wind and ice, 
${U_{10}} = {u_{10}} + i{v_{10}}$ is the wind vector, and 
$\varepsilon $ is the residual. The residual 
$\varepsilon $ can weigh the other contributions for the ice motion, resulting from ocean current, internal ice stresses, Coriolis force and other forcing, as well as error in the ERA5 wind data.
On the sub-daily temporal scale, inertial oscillations and tidal influence were enhanced, but not always synchronously. In the high latitudes of the Northern Hemisphere, semi-diurnal ice oscillations induced by inertial and tidal forcing can be easily confused due to their close frequencies (∼ 2 cycles d−1). However, inertial oscillations are clockwise (CW). Conversely, semi-diurnal tidal forcing would enlarge amplitudes around 2 cycles per day for both CW and counterclockwise (CCW) rotational spectra. This discrepancy allowed us to distinguish the contributions of inertial and tidal forcing by the complex Fourier spectrum. In this study, we used the inertial motion index (IMI) (Gimbert and others, Reference Gimbert, Marsan, Weiss, Jourdain and Barnier2012) and the positive-phase amplitude (PHA) to estimate the inertial oscillation and semi-diurnal tidal components of ice velocity. Additionally, a 3-d Gaussian window was applied to analyse the storm-associated temporal variations of these parameters. For a given buoy location at the present time of 
${t_{pst}}$ (
${x_{pst}}$, 
${y_{pst}}$), the IMI is derived from the Fourier spectrum of ice velocity at the negative inertial frequency 
$ - {f_0}$:
\begin{equation}
IMI = \frac{\left| \hat{W}_{pst}(-f_0) \right|}{\int_{t_0}^{t_{end}} g_{pst}(t) \, dt} \times \frac{4}{\pi \bar{W}_{pst}}
\end{equation} where 
${f_0}$ varies according to the buoy latitude, 
${g_{pst}}\left( t \right) = exp\left( {\frac{{ - {{\left( {t - {t_{pst}}} \right)}^2}}}{{2{{\left( {n{T_0}} \right)}^2}}}} \right)$ is the Gaussian window function, 
${\hat W_{pst}}\left( { - {f_0}} \right)$ is normalised Fourier spectrum at 
${t_{pst}}$, 
${\bar W_{pst}}$ is the mean magnitude of the drift velocity within the Gaussian window. Similarly, the PHA derived from the Fourier spectrum of ice velocity at the semi-diurnal tidal frequency of 
${f_0} = \,2$ cycles per day.
Buoy position data were used to estimate sea ice deformation. However, the sea ice deformation rate definition from the two-dimensional velocity gradient is not applicable for the buoys used in this study. The locations of any three given buoys are unable to form consistent triangles with interior angles >15°, which is the threshold to ensure the reliability of the estimations (Hutchings and others, Reference Hutchings, Heil, Steer and Hibler2012). The distances between any given buoy pairs also have a significant difference, ranging from 10 to 500 km, even for the main buoy array. Meanwhile, such differences will undergo drastic changes over different arrays and through the study period. Thus, it limits the possibility of identifying spatiotemporal variations in ice field deformation based on changes in the relative position of buoys (e.g. Lukovich, Reference Lukovich, Geiger and Barber2017), as well as comparing differences in the regions where various arrays are located. Since the TSE of the buoy, corresponding to the Lagrangian stretching of ice flow, is robust against the initial length scale and has nothing to do with the geometric configuration of the buoy array (e.g. Aksamit and others, Reference Aksamit, Scharien, Hutchings and Lukovich2023). Thus, we adapt TSEs to characterise the ice deformation rate and identify the differences among various buoys and ice regions. Note that, the TSE cannot differentiate the respective contributions from divergence and shear components to the stretching of the ice flow. Instead, it denotes the stretching in the tangent direction of the ice trajectory and is estimated as
\begin{equation}TSE_{{t_0}}^{{t_N}}\left( {{x_0}} \right) = \frac{1}{{{t_N} - {t_0}}}log\frac{{\left| {\dot x\left( {{t_N}} \right)} \right|}}{{\left| {\dot x\left( {{t_0}} \right)} \right|}},\end{equation} and the 
$\overline{\text{TSE}}$ gives a cumulative average of stretching over the period from 
${t_0}$ to 
${t_N}$.
Note that, the choice of integration time defines the average stretching extent can help us focus on and preserve signals of short-term forcing, especially intraday variations. While constructing time series of 
$\overline{\text{TSE}}$, we can also identify the changes caused by the forcing at synoptic (a couple of days) and seasonal scales. Through sensitivity analysis, our main conclusions are not sensitive to integration time, ranging from 1 to 5 d.
Lagrangian observations of ice kinematics with long time series include the components of temporal and spatial changes. Specifically, sea ice motion and deformation rate are regulated by both seasonal and spatial variations in ice conditions and external forcings. Thus, we further use the wavelet analysis on the time series of TSE. Here, the Morlet wavelet is used to take the Continuous Wavelet Transform (Grinsted and others, Reference Grinsted, Moore and Jevrejeva2004), which can be considered a consecutive series of band-pass filters applied to the time series of TSE. Wavelet analysis provides a robust framework for quantifying the temporal changes in the intensity of TSE signals in varying periods. Consequently, it allows for the identification of the predominant forcings, with various frequency characteristics, that drive sea ice deformation throughout the study period.
3. Results
3.1. General atmospheric and ice conditions during the buoy operation
As illustrated in Fig. 2, a relatively moderate westward wind in September and October 2021 increased the tortuosity of the buoy trajectories of the main array. Subsequently, a large-scale low-pressure system entered the Kara Sea from the European sector in November 2021, strengthening the TPD. This is an example of a Nordic sea cyclone intruding into the high Arctic, following an eastern pathway along the Arctic peripheral seas (Kenigson and Timmermans, Reference Kenigson and Timmermans2021). Meanwhile, the Central Arctic Index (CAI), defined as the SLP difference between 90°E, 80°N and 90°W, 80°N (Vihma and others, Reference Vihma, Tisler and Uotila2012), was strongly positive at that time (>20 hPa, Fig. S1), favouring sea ice outflow to the Fram Strait. In December 2021, a stronger high-pressure anomaly over Svalbard produced a northward meridional wind anomaly over the Fram Strait, which hindered Arctic sea ice outflow. In January and February 2022, under frequent cyclonic activity, low-pressure anomalies persisted in the Barents Sea, driving the main buoy array towards the Fram Strait until it approached the ice edge.
Monthly mean sea level air pressure (shading) and 10-m wind (arrows) anomalies from September 2021 to February 2022, with respect to 1979–2021. The cyan lines denote the complete drift trajectories of the main buoy array through the study period, and the red lines denote the buoy trajectories during the respective month.
For the main buoy array, we determined the MIZ using the monthly SIC variability (σ) > 0.11, according to Soleymani and Scott (Reference Soleymani and Scott2023b). Compared to the SIC threshold-based definition of 15–80% for the MIZ, the SIC variability-based definition (Vichi, Reference Vichi2022) contains information on the ice moveability, which causes the ice advection and thus changes in SIC. To consolidate the rationality of the analysis, we also used the temporal scaling laws (
$\alpha $) from absolute dispersion (Eqs 1–2) of the main array to verify the changes in dynamic regime (Fig. 3b). As a result, we divided the whole study period into three stages based on monthly SIC variability and dynamic regimes determined by the changes in 
$\alpha $ (Fig. 3a).
(a) The monthly SIC variability (standard deviation) at the main array. (b) Temporal scaling laws (α) from total (black), zonal (red) and meridional (blue) absolute dispersions of the main array. (c) Bulk-averaged ice temperature obtained from SIMBA measurement. (d) Sea ice thickness (SIT) along the trajectories of the main array. (e) Distance from the main buoys to the ice edge. (f) Ocean bathymetry along the trajectories of the main buoys. The thick line and shading in panels d–f represent the average and standard deviation, and the thin grey lines represent the data derived from individual buoys.
In Stage 1 from August to September 2021, the main array displayed sub-diffusive (
$\alpha $ < 1) zonal and meridional dispersions in MIZ. The relatively high ice temperature (Fig. 3c) also revealed that the sea ice had not yet sufficiently consolidated. In Stage 2 from October 2021 to January 2022, with the onset of ice cooling (Fig. 3c), the main array drifted in PIZ far away from the sea ice edge (Fig. 3e). The total dispersion revealed the sub-diffusive (
$\alpha $<1) dynamic regime during three severe storms, diffusive behaviour (
$\alpha \sim$1) in mid-November 2021 and late December 2021, and elliptic behaviour (
$\alpha \sim$5/3) in January 2022. In Stage 3 in February 2022, after the severe storm Se3, the main array in the MIZ of downstream TPD revealed the dynamic regimes of super-diffusion, dominant ballistic, and enhanced meridional dispersion (
$\alpha $>2). Because most buoys in the main array rapidly drifted towards the Fram Strait and approached the ice edge (<100 km) in February 2022 (Fig. 3e), with 8 buoys drifting over the Yermak Plateau north of Svalbard (Fig. 3f). Thus, both the SIC variability and the changes in ice dynamic regime determined by the 
$\alpha $ confirm that our division of the study period is well founded, at least from the perspective of sea ice kinematics. During the last stage, enhanced oceanic heat flux contributed to the ice melting (Fig. 3d) and tide forcing on the ice floes is expected as the main array drifted over the shallow water.
As shown in Fig. 4, the SLP, 10-m wind speed, and 2-m air temperature along the main buoy trajectories from 10 August 2021 to 1 March 2022 were close to the 1979–2021 climatology, aside from episodic anomalous events associated with cyclones. The low-pressure systems usually corresponded with increased wind speeds and air temperatures. However, the average 10-m wind speed over the entire study period of 2021–2022 (6.4 m ⋅ s−1) was slightly lower than the 1979–2021 climatology (6.9 m ⋅ s−1, Fig. 4b). This indicates that cyclonic activities did not enhance wind forcing throughout the study period. However, the average sea ice speed of 0.10 ± 0.07 m · s−1, derived from the daily remotely sensed ice motion product, was even slightly larger than the 43-year average of 0.09 ± 0.06 m · s−1 (Fig. 4c) due to the enhanced ice moveability.
Atmospheric and sea ice conditions along the trajectories of the main buoys in the study year of 2021–2022 (red) and the average (blue) with standard deviation (shadow) of the climatology: (a) sea level air pressure, (b) 10-m wind speed, (c) sea ice speed, (d) 2-m air temperature, (e) sea ice thickness, and (f) sea ice concentration. The multiyear statistics of ice thickness were derived from 2010 to 2021, and other parameters were from 1979 to 2021.
Satellite-based data shows that the ice thickness along the main buoy trajectories in 2021–2022 was within one standard deviation of the 21-year average. The onsets of ice growth in early November, the onset of melt in early February, and the seasonal evolution of ice thickness were also almost consistent between the study year and the climatology.
3.2. Sea ice kinematics
As shown in Fig. 5a, during the study period, the average drift speed of the main array was 0.14 m · s−1, about twice that of the CB array (0.07 m · s−1) at higher latitudes. Ice speed increased, particularly when the buoy array encountered cyclones or drifted close to the ice edge north of Svalbard in Stage 3. On a daily scale, sea ice drift can be considered a linear response to wind, with a high correlation coefficient of 0.74 ± 0.08 (P < 0.01). The average ice-wind speed ratio (IWSR) and turning angle of the main array were 2.2% and 20.0°, respectively, significantly higher than those of the CA array (1.8% and 16.2°) and the CB array (1.5% and 16.1°) (Fig. 5c and 5d), suggesting a significant wind-driven ice motion in the downstream TPD compared to other regions. During storm events, the response of sea ice to wind forcing aligned more closely with theoretical values (Fig. 5b). The average bias during storm events was 0.02, about half that in the 10 d preceding the storms, indicating that the response of ice motion to wind stress was significantly enhanced during storm events.
Time series of (a) observed daily drift speed for three arrays, (b) residual 
$\varepsilon $ of drift speed estimation from Equation 3, (c) ice-wind speed ratio, (d) turning angle between ice and wind, and the probability of IWSR for all arrays over (e) the entire study period, and in the (f) Stage 1, (g) Stage 2 and (h) Stage 3. The top colour bar represents three stages of buoy operation. The colour shade in panels a–d denotes the storm events, with moderate events in green and severe events in cyan. The horizontal black lines in panel c denote the wind factor, which is used to estimate drift speed. The numbers in panels e–h denote the median values of various buoy arrays.
In Stage 1, the IWSR for the main array ranged from 2.6% to 3.2% (Fig. 5c), exceeding the average of the study period (2.2%) and the CB array in PIZ (2.3%). Note that for the main array in this stage, with SIC variability >0.11, the SIC was close to 100%. This indicated that the SIC threshold-based definition of MIZ may underestimate the extent of MIZ at the end of the melt season from a perspective of ice kinematics. The SIC variability-based definition can better capture the changes in ice kinematic characteristics. In Stage 2, due to the strengthened consolidation of the ice pack, the IWSR distribution narrowed relative to Stage 1 for all three arrays. The average IWSR of the main array (2.1%, except buoy-12 drifted to the MIZ earlier) was comparable to that of the CA array (2.0%) at the same latitude, but larger than that of the CB array (1.5%) at higher latitudes. The relatively low IWSR implies suppressed momentum transfer from wind to ice pack due to intensified internal ice stress. In Stage 3, as the main array drifted over the Yermak Plateau and approached the ice edge, the IWSR of the main array increased to 3.3% and was about triple that of the other two arrays. The higher IWSR confirms the weakened internal ice stress due to the main array drifting in the MIZ with reduced SIC and SIT after the Se3 storm event. The wider distribution of IWSR is likely related to the fact that the oscillating tide forcing potentially increased or reduced ice speed, with a high uncertainty, as the main array drifted over the shallow waters.
The probability distribution of IWSR further quantifies the seasonal and spatial variations of ice kinematics. For all buoys, the median IWSR was higher in Stage 1 (2.7%) than in Stage 2 (1.8%) and Stage 3 (1.8%) (Fig. 5e). However, for the main array, the median of IWSR increased evidently in Stage 3 by 3.3%, relative to those in the two preceding stages (2.7% in Stage 1 and 2.0% in Stage 2, respectively). Thus, the enhanced meridional gradient in ice speed in the downstream TPD likely contributed to greater seasonal changes in sea ice kinematics here compared to the central Arctic Ocean.
As shown in Fig. 6, the internal motion index (IMI) revealed a significant seasonal variation, with the average IMI from all buoys decreasing from 0.215 ± 0.023 in Stage 1 to 0.056 ± 0.005 in Stage 2. Notably, the IMI rapidly reduced by 60% during early October 2021, when the onset of ice freezing occurred. The monthly IMI was significantly correlated with the 2-m air temperature, with correlation coefficients of 0.91, 0.88 and 0.87 (P < 0.01) for the main, CA, and CB arrays, respectively. Because the dropped air temperature would strengthen ice pack consolidation and increase the internal ice stress, which can regulate the ice inertial motion (e.g. Oikkonen and others, Reference Oikkonen, Haapala, Lensu, Karvonen and Itkin2017). The situation deviating from the above relationships was observed for the main array in Stage 3. Despite the typical low air temperature (<−15°C) in Stage 3, the average IMI of the main buoys was about 4–6 times that of the CA and CB array. It is likely due to the enhanced ice pack fragmentation, the reduction of SIC and SIT as the main array drifted closer to the ice edge, and the transition of the ice dynamic regime driven by the extreme storm.
Time series of (a) inertial motion index (IMI) and (b) positive phase Amplitudes (PHA) at semi-diurnal frequency obtained from the 3-d Gaussian window. Panel (c) shows a box-and-whisker plot of the coefficient of determination (R2), indicating the influence of wind forcing on IMI variability, with stars representing the average coefficient and open circles denoting outliers. Power spectral density (PSD) clockwise (CW) and counterclockwise (CCW) rotation components for the wind or ice velocities for the main array in (d) Stage 1, (e) Stage 2, (f) Stage 3. Solid lines and shading show the average and standard deviation. Vertical dashed lines denote the frequency of the local peak of PSD at about 2 cycles per day.
At the short-term scale of 3 d, there was a significant negative correlation between the IMI and the 10-m wind speed, with a correlation coefficient of −0.43 ± 0.11 (P < 0.01). In Stage 1, the coefficient of determination by the wind speed was significantly higher, especially for the main array (R2 = 0.71 ± 0.05) (Fig. 6c). This is because the high wind speeds could promote the mobility of sea ice, which suppresses ice inertial oscillations through enhancing the interactions between ices, especially for the Stage 1, when the sea ice has not fully consolidated. In Stage 2, inertial oscillations were only observed following extreme storm events. The enhanced internal ice stress within a consolidated ice cover suppressed sea ice inertial oscillation, resulting in a reduced coefficient of determination between wind speed and IMI for all arrays.
In Stage 3, robust signals were observed in both the CW and CCW components of ice velocity for the main array within the semi-diurnal and diurnal frequency bands, which can be ignored for wind velocity (Fig. 6d–6f). The PHA increased to approximately twice that observed in the deeper basin during earlier stages or for the other buoy arrays. Thus, during this stage, the sea ice motion at the sub-daily scale of the main array was primarily driven by tidal forcing. The enhanced tidal oscillation was associated with eight of the 12 buoys in the main array that drifted over the shallow water of Yermak Plateau. Additionally, a significant PHA peak was also observed in the CB buoy array in early February 2022 (Fig. 6b), likely due to the temporary influence of tides over shallow bathymetry, with two buoys drifting over Cape Morris-Jesup north of Greenland (Fig. 1).
3.3. Sea ice deformation
Generally, storm events cause significant stretching of the buoy trajectory and 
$\overline{\text{TSE}}$ achieved local maximum during the storm events (Fig. 7a, b). TSE and 
$\overline{\text{TSE}}$ also have the most significant variance during peak stretching events, especially after the storm events, suggesting the variance of ice deformation between buoys would be amplified during enhanced dynamic events. The duration of 
$\overline{\text{TSE}}$ peaks of CA and CB arrays were relatively short compared to those obtained from the main array, indicating that the ice deformation events exhibited higher intermittence for the more compacted ice pack.
(a) Average TSE time series for all buoys (black line) and the main array (red line). (b) Time series of 3 d 
$\overline{\text{TSE}}$ for the main array and their standard deviation (red line and shade), as well as those obtained from the arrays of CA and CB. (c) The lead fraction was obtained from the local region of 50 km from each buoy of the main array. The probability density of 3-d 
$\overline{\text{TSE}}$ for all arrays (d) over the study period and in the (e) Stage 1, (f) Stage 2 and (g) Stage 3. Storm events involving the main array are shown in panels (a)–(c), with dark green for severe events and light green for moderate events.
During Stage 1, the ice pack deformation was relatively active. For the main array, the median of 
$\overline{\text{TSE}}$ was 5.7 d−1, which was higher than that (3.6 d−1) obtained over the whole study period and that (3.8 d−1) obtained by the CB array in the central Arctic Ocean during the same stage (Fig. 7e). In Stage 2, the median of 
$\overline{\text{TSE}}$ for the main array reduced to 2.3 d−1, which was close to that (2.6 d−1) of the CA array but still much larger than that (1.4 d−1) of the CB array (Fig. 7f). This indicated that the ice freezing and consolidation from autumn to winter constrain the ice deformation. This transition was more pronounced in regions (main array and CA array) where first-year or second-year ice dominated than in regions (CB array) where multiyear ice dominated. In Stage 3, once the main array drifted into MIZ over the Yermak Plateau, the enhanced tidal oscillation and warmer seawater prevented sea ice from refreezing after the fragmentation caused by severe storms in late January. 
$\overline{\text{TSE}}$ peaked after the storm event on 18 February 2022, which also corresponds to ice floes drifting into the MIZ and the ice melt onset (Fig. 4a). After this, a sharp decline was observed in 
$\overline{\text{TSE}}$ when the ice began to disintegrate at the ice edge. Thus, the value corresponding to the peak of 
$\overline{\text{TSE}}$ probability distribution in the main array increased moderately to 3.4 d−1 (Fig. 7g), and was larger than twice that from the CB array of 1.6 d−1. Furthermore, it is worth noting that this parameter obtained for the CA array also increased to 3.4 d−1, although other kinematic parameters, such as IWSR and IMI, still maintained at relatively low levels in this stage (Figs. 5c and 6a). This was coincident with a reduced CAI index from 11 hPa in January to 4 hPa in February 2022, resulting in the weakening of previously abnormal CCW atmospheric and sea ice circulations in the Beaufort Gyre region, as shown in Fig. 1. The transition of ice advection regime caused by the change in atmospheric circulation pattern is likely to enhance sea ice deformation and promote lead formation in this region (e.g. Qu and others, Reference Qu, Lei, Liu and Li2024).
During the main array drift, the lead fraction was correlated with the 
$\overline{\text{TSE}}$ on the time scale of 3 d, with a correlation coefficient of 0.56 ± 0.11 (P < 0.01). Specifically, widespread lead formation was observed, associated with the extreme storms Se1 and Se3 (Fig. 7c). However, such a situation was not observed during or after the Se2 event. This is likely because the dominant southerly wind caused by this event opposed the mean ice advection (Fig. 5d), reducing ice divergence and lead formation. Therefore, the relative orientation between the cyclone path and the advection direction of the ice pack plays a decisive role in the formation of leads. Here, we also note that the remote sensing product of a 1-km resolution cannot be used to identify small-scale leads, limiting the evaluation of the impacts of moderate storms on lead formation.
Given that different forcing mechanisms on ice kinematics exhibit distinct timescale characteristics, we employ wavelet analysis for the TSE time series to identify the dominant forcing governing sea ice deformation at various timescales and different stages. During the entire study period, the most significant power of TSE spectra remained at the synoptic scale with peaks at about 10 d, which was likely due to the momentum transferred from the enhanced wind forcing caused by the extreme storms (Fig. 8). Especially, in late January 2022, a power peak at the periods of 8–14 d was observed in all three arrays, attributed to the influence of the Se3 event. The impact of this super cyclone during 21–27 January on sea ice deformation spanned almost the entire TPD, despite its centre being located close to the main array (Fig. S2c). The occasionally enhanced power of TSE spectra at the 2-d temporal scales is likely due to the wind force, as the corresponding power at the same timescale was also observed in the wind speed wavelet spectra.
Changes in the TSE wavelet power spectrum of (a) main, (b) CA and (c) CB arrays. The shade of white indicates the cone of influence, and the black contours indicate the 95% significance level. The horizontal dotted line indicates two notable periods of 0.5 and 10 d.
The impact of semi-diurnal inertial oscillations throughout the entire study period or tidal forcing in Stage 3 on sea ice deformation was always more pronounced in the region with the main array than in those with the other two arrays. In Stage 1, the TSE wavelet power spectrum shows larger semi-diurnal power for the main array than for the CB array, which emphasises the influence of inertial oscillation on the ice deformation for the ice field with relatively unconsolidated characteristics. In Stage 2, the wavelet power for the main array was significantly weakened as the main array drifted in PIZ. In Stage 3, the TSE wavelet power intensity of the main array was enhanced again due to strong tidal forcing as they drifted over shallow waters.
Calculating the average wind speed and TSE wavelet power spectrum at the corresponding time scales, we focus on the regime of sea ice deformation at significant temporal scales, including sub-daily scale and synoptic scale (7–14 d).
On the sub-daily scale corresponding to inertial or tidal frequency (Fig. 9a), the peaks of TSE average variance were almost an order of magnitude larger than those of wind speed, especially in Stage 1. This indicates that sea ice localised stretching was primarily caused by inertial oscillation. In Stage 2, the peaks of TSE average variance were only observed after the storm events and correspond well with the IMIs (Section 3.3), suggesting that sea ice inertial motion stimulated by severe storms forced ice stretching. In Stage 3, the peaks of TSE average variance did not strengthen at the inertial frequency, while there were stronger inertial oscillations in ice motion (Fig. 6c). This is likely due to the enhanced and complex ocean forcing, including the accelerated southward average current in the strait region (e.g. Lei and others, Reference Lei, Heil, Wang, Zhang, Li and Li2016) and tidal currents over shallow bathymetry (e.g. Fer and others, Reference Fer, Müller and Peterson2015).
Changes in average variance of TSE and wind speed at (a) sub-daily scale and (b) synoptic scale (7–14 d) of the main buoy array.
On the synoptic scale of 7–14 d (Fig. 9 b), wind plays a dominant role on sea ice deformation (r = 0.91, P < 0.001). During the entire drift, there were three significant power peaks in wind speed spectra on this scale, corresponding to three storm events of Mo1, Se1, and Se3, respectively. Weak peaks were observed during the Se2 event because the wind direction opposite the main array drift enhanced the ice compaction. In addition, during the Mo1 event, the correlation of TSE with the wind speed was insignificant at the 0.05 level, which indicates this synoptic event did not significantly regulate ice deformation. This is because inertial oscillation dominated the ice deformation at that time, which likely led to the decoupling of sea ice kinematics with wind force. During the Se1 and Se3 events, the TSE peaks averaged 7–14 d, coinciding with the peaks of wind speed but with delays of about 3 d for Se1. Unlike the previous events, during the strongest storm event, Se3, the response of TSE to wind forcing was enhanced and timely without any delay.
4. Discussion
Since the CB array drifted in the central Arctic Ocean during the entire study period, its kinematics represent PIZ conditions and serve as a benchmark to identify changes or anomalies in the ice kinematic parameters in other regions. As shown in Figure 10, higher ice speed and IWSR were observed in the main and CA arrays, being approximately 1.5–2.5 times those in the CB array. This is likely because the main and CA arrays were located at lower latitudes, where the sea ice consolidation was relatively weaker. After the Se3 events, the differences in both ice speed and IWSR between the main array and CB array significantly increased to about 20 and 10 times, respectively. The transition of ice conditions at the main array from PIZ to MIZ could explain the amplified differences in its ice kinematics compared to those of the CB array in early February 2022.
The ratios for four ice kinematic or deformation parameters of the main and CA arrays relative to the CB array: (a) ice speed, (b) IWSR, (c) IMI and (d) 
$\overline{\text{TSE}}$. Blue lines represent the benchmark obtained from the CB array. The colour shade in panels (a)–(d) denotes the storm events, with moderate events in green and severe events in cyan.
Focusing on the impacts of synoptic events, enhanced differences in sea ice kinematics between arrays were observed during moderate storms than during the three severe storms. Specifically, the ice speed of the main array was 2.8 times that of the CB array during moderate storms and was 2.3 times that of the CB array for 5 d after these storms. The corresponding magnifications for severe storms are 2.2 and 1.9 times, respectively. Similarly, the IWSR of the main array was 1.7 times that of the CB array during moderate storms, which was also slightly larger than that (1.5 times) during severe storms. The seemingly counterintuitive results highlight the localised impacts of moderate storm events. Strong storm events, identified based on the atmospheric conditions at the main array, were associated with strong and larger-scale cyclones. For instance, the impact of the Se3 event extended across the Arctic basin, influencing the sea ice dynamics beyond the main array, reaching the regions with other buoy arrays. In contrast, the moderate storms were more likely to be experienced only by the main array, with a more limited extent of impact on ice kinematics. On average, the 
$\overline{\text{TSE}}$ of the main array was 1.8 (1.4) times that of the CB array during moderate (severe) storms. However, during the 5-d period after the storms, the 
$\overline{\text{TSE}}$ ratio of the main array to the CB array for the severe storms (1.6 times) was larger than that for the moderate storms (1.1 times). This means that strong storms have a more sustained impact on sea ice deformation relative to moderate storms.
Since this work presents a case study based on observation data collected during one ice season in 2021–2022, the representativeness of the sea kinematics during the study year is worth discussing. Using the NSIDC sea ice motion product and ERA5 reanalysis wind speed along the buoy trajectory, we compared the IWSR obtained from 2021 to 2022 with those from the past 43 years. To construct multiyear averages, we define two periods of the last 10 years, between 2011 and 2022, and the full 43 years, as the sea ice has experienced significant changes in the last decade.
As shown in Fig. 11, the IWSR probability distributions reveal differences between the study year and the past 10 or 43 years. Although the wind speed and its probability distribution didn’t reveal significant changes (Fig. S3), the median and proportion of high IWSR in all stages and regions were larger in the past 10 years than those in the past 43 years, revealing the enhanced response of ice kinematics to the wind forcing in the recent years. Furthermore, the IWSR in 2021–2022 had higher values than those in the past 10 years, even for the CB array in the central Arctic Ocean. This indicates that there is also a slight enhancement of the ice response in the study year, even compared to that of the past 10 years. Especially for the main array in Stage 3, the IWSR in the study year was higher (Fig. 11c) even though fewer storms occurred (Fig. S3c), with relatively high values of >0.02 occupying 38% of the time, which was 1.3 and 1.7 times those of past 10 and 43 years (29% and 22%), respectively. This was likely caused by the enhanced fragmentation and mobility of the ice pack after the super cyclone in this stage.
Comparison of ice-wind speed ratio probability density between 2021 and 2022 and the past 10 (2011–2022) or 43 years (1979–2022). The IWSR is calculated using the daily products of NSIDC sea ice motion and ERA5 reanalysis wind speed for the study year and previous years for consistency.
Compared to the IWSR obtained from the past 10 or 43 years, we can further conclude that the sea ice kinematics of study year were jointly affected by the seasonal evolution and large-scale spatial differences in ice conditions (similar to past decades) and abnormal sea ice conditions and extreme synoptic events occurring in the study year.
5. Conclusions
In this study, we employed a buoy array to investigate the dynamic characteristics of sea ice during TPD from the MIZ north of the Kara Sea, through transition to PIZ, and finally to the sea ice edge north of Svalbard and Fram Strait. In addition, the data collected by the other two arrays drifted in the upstream TPD and the central Arctic Ocean, were used to reveal the disparities in ice kinematic characteristics across the PIZ and MIZ. Our analysis demonstrated how atmospheric and oceanic forcing and local ice conditions influenced the seasonal evolution and spatial heterogeneity of sea ice motion and deformation.
A comprehensive suite of metrics was used to describe ice kinematics. Along buoy trajectories, linear regression of ice velocity against wind speed, the ratio of their absolute values, rotational spectroscopy analysis and complex Fourier transform in a sliding Gaussian window of ice velocity were combined. The comparisons of these ice kinematic parameters illustrate the response mechanisms of sea ice motion to wind and other forces, as well as the changes in the ice pack itself. Additionally, the tangential stretching exponent of the buoy trajectory and its wavelet analysis were obtained to describe and compare ice deformation across large scales.
The results reveal that in the downstream TPD, sea ice advection, temporal changes and spatial differences in sea ice, atmospheric and oceanic conditions, and synoptic events collectively shape the seasonal variations in sea ice kinematics and deformation. From mid-August 2021 to late February of the following year, the ice kinematic characteristics proceeded through three stages: a relatively active stage at the end of ice melt season through late September in the MIZ north of the Kara Sea, a freezing stage lasting until late January in the PIZ, and a stage with rapidly enhanced dynamics close to the ice edge in February. The IWSR, IMI and TSE obtained from the main buoy array reduced from 2.7%, 0.22 and 5.7 d−1 in Stage 1 to 2.0%, 0.07 and 2.3 d−1 in Stage 2 and increased again to 3.3%, 0.08 and 3.4 d−1 in Stage 3. In Stage 3, the ice kinematics in the region close to the ice edge north of Svalbard and Fram Strait differed from those in the central Arctic Ocean, with the amplified IWSR and IMI, which were about 3.5 and 2 times that in the central Arctic Ocean.
Extreme cyclones amplify local sea ice deformation by 2–3 times, with effects lasting for 3–5 d after the events and probably driving the transition of the ice dynamic regime. We also observed that although the extreme storm event occurred in the Barents Sea in late January 2021, the impact of this synoptic event on sea ice kinematics extended to the North Pole, as the synoptic system spanned almost the entire Atlantic sector of the Arctic Ocean. Moreover, the shallow bathymetry of the Yermak Plateau north of Svalbard enhanced the tidal-induced oscillations of sea ice motion. These factors contributed to the significantly different seasonality in sea ice kinematics in the downstream TPD compared to the upstream or the central Arctic Ocean.
However, a case study using observations obtained from one year is still insufficient. Continuous buoy monitoring, combined with remote sensing products of other geophysical parameters of sea ice, is necessary to complete our understanding of the kinematic properties of sea ice and its responses to climate and sea ice changes in this region. Especially in the MIZ close to the ice edge of the north Barents Sea or Fram Strait, the understanding of the response mechanism of the ice kinematics to severe storm events is often derived from short-term observations, for instance, this study and observations of N-ICE 2015 (Itkin and others, Reference Itkin2017). Whether there are seasonal differences in this response mechanism and how it relates to seasonal changes in local ice conditions and the relative position of cyclone paths are still deficiencies in knowledge that deserve further investigation. How the kinematics and deformation of sea ice affected by cyclones can be fed back to local changes in sea ice and waters beneath the ice, from perspectives of physics, ecology, or biogeochemical cycles, is also a scientific issue that deserves attention and requires interdisciplinary field investigations and numerical simulations to seek the solutions to deconstruct the issue.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.10099.
Acknowledgements
This work was financially supported by the National Natural Science Foundation of China (Grant Nos. 42325604, 42206259 and 42476250), the Ministry of Industry and Information Technology of China (Grant No. CBG2N21-2-1), and the Program of Shanghai Academic/Technology Research Leader (Grant No. 22XD1403600).
Open access
