1. Introduction
Since 10 April 2010, CryoSat-2 has been in a high inclination orbit and has supplied a continuous radar altimetry record of Arctic Ocean surface elevation (Wingham and others, Reference Wingham2006; Laxon and others, Reference Laxon2013). During winter months, the Arctic Ocean surface freezes to form sea ice; a barrier between the cold atmosphere and liquid ocean. Sea ice is distinguishable from space, particularly within longer microwave lengths that enable the continuous monitoring of sea-ice extent (Kern and others, Reference Kern2019). Sea ice is at the mercy of winds and ocean currents (Thorndike and Colony, Reference Thorndike and Colony1982; Heorton and others, Reference Heorton2019) and circulates throughout the Arctic Basin and into neighbouring seas (Kwok, Reference Kwok2009; Babb and others, Reference Babb, Landy, Barber and Galley2019; Wei and others, Reference Wei, Zhang and Wang2019). The circulation of solid sea ice causes it to fracture into sea-ice floes, presenting an opportunity to observe key aspects of its volume (Ricker and others, Reference Ricker, Kauker, Schweiger, Hendricks, Zhang and Paul2021). Sea-ice floes present visible patterns that are tracked to allow for the calculation of sea-ice drift (Liu and others, Reference Liu, Martin and Kwok1997; Kwok and others, Reference Kwok, Schweiger, Rothrock, Pang and Kottmeier1998; Lavergne and others, Reference Lavergne, Eastwood, Teffah, Schyberg and Breivik2010). The cracks between floes at leads reveal the ocean surface in which the sea-ice floes float. CryoSat-2 is able to measure the elevation of both sea ice and leads and by considering the difference between them, the radar freeboard (Peacock and Laxon, Reference Peacock and Laxon2004; Giles and others, Reference Giles, Laxon and Worby2008; Dawson and others, Reference Dawson2022), and sea-ice thickness is estimated (Laxon and others, Reference Laxon, Peacock and Smith2003, Reference Leppäranta2013; Landy and others, Reference Landy2022).
During warmer summer months, sea ice melts, both into the ocean below and within melt ponds upon the surface. Until recently, the presence of melt ponds has hindered the consistent measurement of sea-ice thickness. Melt ponds present a reflecting surface that appears similar to the open ocean in radar altimetry signals, but at the level of the ice floe surface (Kwok and others, Reference Kwok, Cunningham and Armitage2018). Recent advances have allowed for the separation of these two surfaces (Dawson and others, Reference Dawson2022), allowing measurement of radar freeboard using CryoSat-2 within summer months. Radar freeboard, when combined with estimates of snow thickness, allows for the calculation of sea-ice thickness through the balancing of hydrostatic equilibrium. Despite inherent uncertainties in measuring centimetre-sized differences in elevation from an altitude of 719 km (Nab and others, Reference Nab2023), a continuous record of sea-ice thickness is now available that has a good agreement with multiple airborne and in situ observations (Landy and others, Reference Landy2022).
The extent of Arctic sea ice is one of the largest signals of anthropogenic climate change (Notz and Stroeve, Reference Notz and Stroeve2016), with the summer minimum sea-ice extent presenting a record low in 2012 (Parkinson and Comiso, Reference Parkinson and Comiso2013) and the prospect of a summer sea-ice free Arctic ocean now likely (defined as under 1 million km2, Kim and others, Reference Kim, Min, Gillett, Notz and Malinina2023). Despite the clear climatic signal, efforts in forecasting sea-ice extent have had limited success (Stroeve and others, Reference Stroeve, Holland, Meier, Scambos and Serreze2007; Andersson and others, Reference Andersson2021), with theoretical barriers to predictability shown (Bushuk and others, Reference Bushuk, Winton, Bonan, Blanchard-Wrigglesworth and Delworth2020). While the extent of sea ice is a two-dimensional measurement, the controlling processes of sea ice rely on the three-dimensional ice volume and the energy balance of ice growth and melt. This volume and energy balance has been a key aspect of climate models since Thorndike and others (Reference Thorndike, Rothrock, Maykut and Colony1975), with modern sea-ice models considering a complex collection of physical processes including the phase change of sea water, chemical reactions that arise from ocean salinity and the non-Newtonian flow of fracturing sea ice (Heorton and others, Reference Heorton, Feltham and Tsamados2018; Keen and others, Reference Keen2021).
Since early efforts at modelling sea ice as part of global climate, the accurate description of drift patterns has been known to be crucial (Feltham, Reference Feltham2008). Thicker sea ice grows slower than thin ice in the same environmental conditions (Leppäranta, Reference Leppäranta1993), so the redistribution of thin sea ice to colder regions, and vice versa, affects the seasonal growth in sea ice (Kimura and others, Reference Kimura, Nishimura, Tanaka and Yamaguchi2013). Over individual sea-ice floes, the growth of sea ice is more complex, with snow cover (Sturm and others, Reference Sturm, Perovich and Holmgren2002) and ocean processes (Heorton and others, Reference Heorton, Radia and Feltham2017) also contributing (see Nicolaus and others, Reference Nicolaus, Tavri, Tsamados, Wagner, Watkins, Webster and Wendisch2022, for an MOSAiC expedition overview). Over larger length scales similar to the continuum scale used within sea-ice climate models, multiple efforts have been made to assess the volume changes to sea ice. Holland and Kimura (Reference Holland and Kimura2016) used drift and concentration data to reveal regions where dynamics dominate observed changes to sea-ice concentration in both the Arctic and Antarctic. Kwok and Cunningham (Reference Kwok and Cunningham2016) used gridded ice thickness and drift observations to investigate the link between strain invariants and changes to thickness, showing that drift induced deformation accounts for up to 56% of ice volume change in the Lincoln Sea. Ricker and others (Reference Ricker, Kauker, Schweiger, Hendricks, Zhang and Paul2021) similarly separated the dynamic and thermodynamic contributions to sea-ice thickness change but over the wider Arctic up to 82∘ N during winter months. By isolating the thermodynamic changes to sea ice, and by relating this to reanalysis-based atmospheric temperatures, the role of oceanic heat was estimated. Li and others Reference Lin, Pickart, Heorton, Tsamados, Itoh and Kikuchi(2021) used gridded drift and thickness data to investigate changes in Arctic minimum sea-ice volume with a focus upon the Fram Strait. Recently, Anheuser and others (Reference Anheuser, Liu and Key2022) have expanded this approach to include passive microwave brightness estimates of thermodynamic growth rates to constrain Arctic wide winter deformation rates.
Within the Arctic system, there exist separate seas that have contrasting geographies and characteristics. The Beaufort Sea is widely studied due to the presence of an ocean gyre (Giles and others, Reference Giles, Laxon, Ridout, Wingham and Bacon2012; Lin and others, Reference Lin, Pickart, Heorton, Tsamados, Itoh and Kikuchi2023) and high variability in sea-ice circulation (Mallett and others, Reference Mallett2021). Adjacent to the Beaufort Sea ice is the Lincoln Sea directly north of Greenland. This area close to the north pole contains the thickest and oldest sea ice advected by the Beaufort Gyre to the Chukchi Sea and throughout the basin and is heavily compressed into ridges (Kwok and Cunningham, Reference Kwok and Cunningham2016). On the opposite side of the Arctic Basin are the Siberian, Laptev and Kara seas, that in recent years have become increasingly ice-free in summer (Lukovich and others, Reference Lukovich2021), opening the north-east passage. Despite the low summer ice volume in these areas, due to the processes of ice export and divergence, high volumes of ice are formed here each winter (Cornish and others, Reference Cornish, Johnson, Mallett, Dörr, Kostov and Richards2022). On the periphery of the Arctic Ocean are the Greenland and Barents seas and Baffin bay. These regions have seasonal ice cover that has extensive winter marginal ice zones (MIZs), the outer band of the sea-ice cover where open ocean waves interact with sea ice (Brouwer and others, Reference Brouwer2022). Within this region, there are many interactions between sea ice, ocean and atmosphere (Heorton and others, Reference Heorton, Feltham and Hunt2014) and sea-ice growth and melt (Thomas, Reference Thomas2017). Outside the Arctic Basin, there are the seasonally ice covered Labrador and Bering Seas, Sea of Okhotsk and Hudson Bay (Meier and Stewart, Reference Meier and Stewart2023). Within the Basin, there are the ice covered channels of the Canadian Archipelago, where sea ice can persist year round and flow along them (Agnew and others, Reference Agnew, Lambe and Long2008).
As sea ice drifts and circulates through the Arctic, large volumes of sea ice are moved between regions. The Fram Strait is the most widely studied flux route, as it accounts for 90% of sea ice that leaves the Arctic Basin (Sumata and others, Reference Sumata, de Steur, Gerland, Divine and Pavlova2022). While the volume of sea ice that drifts through the Fram strait is less than 1/10th of the seasonal cycle in basin wide sea-ice volume, it has been shown to be correlated to minimum sea-ice extent, showing the link between basin wide drift characteristics and total volume (Williams and others, Reference Williams, Tremblay, Newton and Allard2016; Li and others, Reference Lin, Pickart, Heorton, Tsamados, Itoh and Kikuchi2021). The fluxes between the sub-regions of the Arctic basin have been used to investigate links between winds and sea-ice state. By recording both the areal and volume flux in and out of the Beaufort sea (Mallett and others, Reference Mallett2021) and through the wider Arctic, the link between atmospheric circulation modes such as the Arctic oscillation and summer sea-ice survivability is revealed (Stroeve and others, Reference Stroeve, Maslanik, Serreze, Rigor, Meier and Fowler2011).
The current availability of Arctic wide satellite-based observational data presents an opportunity to close the budget of Arctic sea-ice volume from observational data and reveals the budget components that contribute to the seasonal cycle in volume, over the entire Arctic system and within the constituent regions. Here, we used the year-round sea-ice thickness data from Landy and others (Reference Landy2022) along with sea-ice drift and concentration data (see Section 2). These data were combined using a method similar to previous investigations of sea-ice thickness and drift (Holland and Kwok, Reference Holland and Kwok2012; Holland and Kimura, Reference Holland and Kimura2016; Kwok and Cunningham, Reference Kwok and Cunningham2016; Li and others, Reference Lin, Pickart, Heorton, Tsamados, Itoh and Kikuchi2021; Ricker and others, Reference Ricker, Kauker, Schweiger, Hendricks, Zhang and Paul2021; Anheuser and others, Reference Anheuser, Liu and Key2022), but with an intent of accounting for the entire volume of sea ice throughout the study period of October 2010–May 2022. A key aspect of closing the volume budget is accounting for data coverage and uncertainty. Data handling and coverage are described in Section 2, with the budget analysis and the method developed to propagate data uncertainty described in Section 3 and the appendices. Results are presented for the mean Arctic state (Section 4.1) and inter-annual variability (Section 4.2). The strength of the results in respect to the closure of the volume budget is discussed in detail in Section 4.3. The regional variations and key aspects of the volume budget are described in Section 4.2 and the role of input data and results discussed in Section 5.
2. Data
To consider the changes to Arctic sea-ice volume continuously over the period in which CryoSat-2 thickness data were available (November 2010–April 2022), we required observations of its volume, concentration and drift at a daily time scale. The aim of this study was to access published and validated observationally derived data sources and to obtain results with no modification to the processing algorithms. Daily data were chosen in order to match the temporal resolution of both the available ice concentration and drift data and to enable a consistent calculation of time derivatives. The summer and winter periods were set as 16 April–15 October (summer) and 16 October–15 April the following year (winter), to allow our daily calculations to be consistent with other winter only monthly sea-ice thickness data (Tilling and others, Reference Tilling, Ridout and Shepherd2018). While this choice of period may cause some positive changes to be incorporated within the summer season, making this summation consistent to the true maximum and minimum sea-ice volumes will require a different definition of growth and melt seasons for each region and each year and the careful consideration of the different length of the summation period when performing seasonal comparisons. The code developed for this study is made freely available and its modular form allows for easy adaptation to use alternate data sources and time periods. For regional analysis, the Arctic Basin regions given by the NSIDC were used (Meier and Stewart, Reference Meier and Stewart2023) and combined for regions that present similar characteristics: Beaufort and Chukchi sea, Siberian with the Kara and Laptev seas.
2.1. Sea-ice drift
Selecting sea-ice drift data that gives reliable divergence data remains a challenge. To that effect, we used the sea-ice drift products of Polar Pathfinder (v4, Tschudi, Reference Tschudi, Meier, Stewart, Fowler and Maslanik2019) and OSISAF (OSI-455, Lavergne and Down, Reference Lavergne and Down2023; OSISAF, Reference Parkinson and Comiso2022). These data incorporate floe tracking using satellite images from a variety of sensors, along with buoy drift data for Pathfinder and data filling using a wind driven free drifting sea-ice model. Both drift data sources have been used in order to show results that are consistent between the two. Pathfinder data cover the time series of thickness data used but have issues when calculating divergence due to the assimilation of small-scale buoy drift data (Szanyi and others, Reference Szanyi, Lukovich, Barber and Haller2016). This issue may be the cause of artefacts found in our results when using Pathfinder drift data as commented on in Section 5. The OSISAF data and our results derived from them are all free from these artefacts, but the data are only available up to 31/12/2020. An alternate OSISAF product (OSI-405) has data up to 2023 but has no data for summer months (May–September). Pathfinder is supplied on a 25 km resolution and OSISAF a 75 km resolution. Both datasets have been used for all calculations, with supplemental figures supplied.
The Pathfinder data are created using three combined methods. First, a floe tracking algorithm is applied to AVHRR visible and infrared 50 km images; SMMR, SSM/I, SSMIS and AMSR-E passive microwave 25 km and 12.5 km resolution images. Second, a collection of IABP ice-tethered buoy drift tracks are assimilated. Third, a free drifting ice drift equation is produced using NCAR/NCEP reanalysis 10 m winds (Tschudi, Reference Tschudi, Meier, Stewart, Fowler and Maslanik2019). Ice drift x and y components are supplied along with a single value for the variance in ice drift speed σ 2. We converted this value to a per component value with 
$\sigma/|\bf{u}| = $ 
$\sigma_u/u = \sigma_v/v$.
The OSISAF drift data use a similar but independently created algorithm that is based on only two methods. First is the use of the SSM/I, SSMIS, AMSR-E and AMSR2 passive microwave instruments. Secondly for summer months (June–September), a free drift equation is forced by ERA5 reanalysis 10 m winds to provide drift information where no image-based methods are available (Lavergne and Down, Reference Lavergne and Down2023; OSISAF, Reference OSI SAF2022). Ice drift x and y components (see Fig. 1) are supplied along with a single value of vector error σ converted similarly to the Pathfinder data. For the case of data filling using the free drift model, a constant uncertainty value of 3.9 km per 24 hour period or 0.45 m s−1 is given.
Averaged example ice state data for the budget calculations 2010-2020. Top are the average 15 Oct and 15 Apr sea-ice thicknesses, bottom are the average OSISAF winter ice drift and averaged combined MIZ summer ice loss. Arrows show the average summer and winter drift patterns.
2.2. Sea-ice thickness
In order to close the volume budget over the entire study period, the thickness of Arctic sea ice needs to be known continuously over the whole seasonal cycle. For this, we use the thickness dataset of Landy and others (Reference Landy2022) updated to July 2022, see Fig. 1. These data use the CryoSat-2 SIRAL synthetic aperture radar altimeter, with a novel machine learning algorithm to classify sea-ice floes from the leads that separate them and derive sea-ice radar freeboard, over the full Arctic Basin (Dawson and others, Reference Dawson2022). The freeboards are corrected for elevation biases caused by melt ponds in the summer months (Landy and others, Reference Landy2022). The bias ranges from just a few cm during early melt (May–June) up to 10 cm for the roughest sea ice in mid-summer (July–August). The total freeboard uncertainty is highest, up to 40% of the corrected freeboard, in mid-summer, and constitutes 80–90% of the total ice thickness uncertainty. Freeboard is converted to estimates for sea-ice thickness using snow loading information from the reanalysis system SnowModel-LG (Liston and others, Reference Liston2020). Median ice thickness uncertainty for summer months is estimated by Landy and others (Reference Landy2022) to be 33% of the thickness for first-year ice and 40% for multi-year ice. These data agree well with existing in situ and airborne sea-ice thickness measurements in all seasons. The data are supplied at a 80 km resolution with 15 day (semi-monthly) averages. For daily volume budget calculations, the midpoint of each 15 day window is defined by the original data value, and all other days are linearly interpolated between them. The linear interpolation causes artificially high data covariances, which have been considered when propagating uncertainties (see appendix B).
2.3. Sea-ice concentration
We use the widely used daily ice concentration from the NSIDC NASA Team (Cavalieri and others, Reference Cavalieri, Parkinson, Gloersen and Zwally1996), also derived from the SSMR and SSMIS passive microwave instruments. These data are supplied on a 25 km Polar Stereographic grid. Alternate sea-ice concentration is available from several sources including the NSIDC Bootstrap algorithm, the combined NSIDC CDR product, OSISAF and ESA-CCI. Kern and others (Reference Kern2019) present an in-depth analysis of these data. For this study, our focus was on the volume budget, so our requirements for ice concentration data are for reliability at near 100% concentration and reliability over the time series, both given by the NSIDC NASA Team data. The ice concentration data are used along with the sea-ice thickness data to calculate the sea-ice volume over the ice covered fraction of the data grid. As no information on the uncertainty of the concentration data is given, and as it is beyond the scope of this study to estimate it independently, we use the uncertainty of the thickness data to represent the combined volume. The regions where the ice concentration data accuracy is likely to be most significant are for lower ice concentrations at the sea-ice edge. The ice edge has issues with data coverage that we consider specifically as described in Section 3.1.1.
3. Methods
To evaluate the components of the Arctic sea-ice volume budget, the datasets for thickness, concentration and velocity are regridded and combined. To ease the use of vector data, the native ice drift velocity grid is used with components aligned to the native grid. This gives the calculations using the Pathfinder a 25 km resolution and OSISAF a 75 km resolution. First we calculate the areal unit volume of ice within a grid cell (m3 of ice per m2 of ocean) by multiplying the thickness and concentration to get unit volume 
$V = C\times T$. Taking the time derivative of the unit volume over 2 days, in units of m s−1, we can equate the key components of the volume budget using
\begin{equation}
\frac{\mathrm{d}V}{\mathrm{d}t} = - \mathbf{u}\cdot{\nabla} V - V\nabla\cdot\mathbf{u} + \text{residual}
\end{equation} where the terms, in order, represent: intensification, the static time derivative of volume change; advection, the result of moving ice with a spatially variable volume; divergence, the result of spatial varying velocity and sea ice as a compressible 2-d fluid; and residual, the change in ice growth unaccounted for from dynamical terms. The frame of reference for the following expansion of variables is given by indices 
$i,j$ with spatial distances as 
$\mathbf{x}=(x_i,y_j)$ and velocity components as 
$\mathbf{u} = (u,v)$.
The three terms, intensification, advection and divergence, are calculated numerically using a second-order derivative. For a time index t, the intensification is thus
\begin{equation}
\text{Int}_t=\frac{1}{2 \Delta t}(V_{t+1}-V_{t-1})
\end{equation} with 
$\Delta t$ the number of seconds in a day. For the spatial derivatives required in the dynamics terms, we have for the x component at point i (and with an identical form for y at point j),
\begin{equation}
\begin{aligned}
V_{x}&=\frac{1}{2\Delta x_i}(V_{i+1}-V_{i-1}) \\
u_{x}&=\frac{1}{2\Delta x_i}(u_{i+1}-u_{i-1})
\end{aligned}
\end{equation} where Vx is the x component of the gradient of the volume data, ux is the x component of the divergence of the velocity data, with 
$\Delta x_i$ the local spacing of the x direction of the data grid in meters. The equations for v and y follow similarly. The uncertainty for all components is propagated as described in appendix A.
On a length scale shorter than the native data resolution and implemented data smoothing (see Section 3.1.2), the budget terms can be indicators of observed phenomena, though the methods used in this study do not directly observe these phenomena. Positive divergence within the drift data is often linked to the opening of leads with sea-ice area and the local reduction in sea-ice volume captured by ice concentration data and represented within these data as negative divergence only. Leads, however, can often open and close over daily time scales (Ludwig and others, Reference Ludwig, Spreen, Haas, Istomina, Kauker and Murashkin2019) or quickly freeze (Heorton and others, Reference Heorton, Radia and Feltham2017) and may not be present during the overpass of a satellite sensor (Meier and Stewart, Reference Meier and Stewart2020). Sea-ice concentration datasets also disagree notably at high concentrations (Kern and others, Reference Kern2019).
Negative divergence or compression can cause the closure of leads and the formation of pressure ridges, directly causing thicker sea ice, and thus can be interpreted as being part of the divergence component of the budget. Due to the geometry of ice floes, it is theorized that ridges can occur during wider scale shearing (Wilchinsky and Feltham, Reference Wilchinsky and Feltham2006), which will not be within the divergence term and thus make up a part of the residual. We, therefore, suggest that regions of positive divergence can be associated with a high likelihood of ridging, but not exclusively so, with the residual term in such locations likely to include ridged ice. Anheuser and others (Reference Anheuser, Liu and Key2022) attempt to partition these components, with their deformation and divergence having similar spatial and temporal variability, though differences are within the uncertainty of these terms. If the volume contribution of pressure ridges is to be recorded as part of the residual component in our results, then the thicker ridged ice is required to be included within the thickness data used (see Section 2). However, it is currently not clear how well a single pressure ridge is represented within a CryoSat-2 radar altimeter footprint (Xia and Xie, Reference Xia and Xie2018; Belter and others, Reference Belter2020). It may be possible for future work to include the consideration of both the ice thickness and area separately within Eqn (1) to provide insight of greater sophistication, though this is beyond the scope of this study, and may also require a sub-grid cell distribution and horizontal transport of budget terms.
3.1. Processing steps
The data processing was performed using an object-orientated modular Python script. The coding structure enables flexibility in grid structure, input and output data format and the time scale of output data. This flexibility is particularly important for using a variety of satellite input data (e.g. both Pathfinder and OSISAF drift products) and for the accumulation of uncertainty estimates (see appendix A). The key processing steps and considerations are as follows.
3.1.1. Data loading and masking
Prior to the calculation of budget terms, each data source is analysed to show the available dates. These are recorded with data loading objects to allow for the matching of data from multiple sources and interpolation between dates where required. Data grid objects are defined for input and output formats (including the ice transport gates), and regridding objects are defined to move all data to the ice drift velocity grid. Data and uncertainty values are loaded for each day.
For consistent calculation of spatial derivatives, near land points and the CryoSat-2 pole > 87.5 N are masked in all data. For the calculation of the budget components, data coverage was considered. The total sea-ice extent was taken from the ice concentration data and used as a reference for where all data were required. We found the velocity data were defined over the whole concentration fields, while the thickness data often had gaps. All regions covered by ice concentration data but with no thickness data present are filled with a constant 0.2 m thickness. This thickness was chosen as the missing data are typically from regions near the sea-ice edge. As the first stages of sea-ice growth (the collection thickness of grease ice) are typically of order 0.2 m thick (Heorton and others, Reference Heorton, Radia and Feltham2017), this value was chosen to capture the growth and loss of sea within the MIZ and thus close the total volume budget over the whole Arctic basin over a growth and melt season (see Section 4.2). The volume budget calculations from the regions with filled thickness data are recorded as separate variables in the final data files. Finally, all pixels that contain data while the previous day contained none are recorded as new_ice (along with the opposite case, old_ice) to account for all changes to sea-ice volume. The variable MIZ presented later in this study is the sum of new_ice, old_ice with the intensification from filled ice thickness data cells (int_t). Due to the assumption of filling ice thickness data, the magnitude of int_t is used for the uncertainty of the MIZ variable. While only int_t is included with the MIZ term in this study and is most significant during summer melt (see Figs 1 and S1), the div_t and adv_t terms are key to avoid missing data. As the dynamical terms in Eqn (3) rely on neighbouring grid cells, any gap in unfilled thickness data results in many more missing budget data as discussed in Section 4.1. The extent of filled data has a consistent seasonal cycle (Supplemental Fig. S1). During July and November when ice retreat and advance is at the greatest rate, the filled data can account for up to 10% of the total sea-ice area, though only up to 3% of the total volume (S1b). The majority of the year has less than 1% of the sea-ice volume from filled regions. During July, the MIZ term is at a maximum, with int_t covering the MIZ (Fig. 1). The new_ice and old_ice terms typically cancel each out as the ice edge advances and retreats within a month (S1a).
3.1.2. Data smoothing
After defining data grids and input data formats, the data smoothing options are set. In order to maintain sensible values of divergence and to provide the most reliable budget components for the whole Arctic, a variety of smoothing techniques were experimented with a local grid resolution dependent variable-sized square box convolution kernel used for consistency between different input data sources. To seek information on the total budget of Arctic sea-ice volume, a balance is sought between the detail within the data available and noise, which is a particular issue when taking spatial derivatives, a usage beyond the initial design of these data. To extract information of terms relating to spatial derivatives (described fully in Section 3), smoothing was applied to each of the data as performed by Holland and Kimura (Reference Holland and Kimura2016). The smoothing radius applied was chosen by increasing the value for each input variable until the values for spatial derivatives (Eqn (3)) were typically of the same sign for the majority of neighbouring grid cells, aiming to keep smoothing to a minimum if possible. Sea-ice drift data are smoothed with a 150 km radius, thickness a 100 km radius and concentration a 60 km radius. After smoothing, thickness and concentration data are regridded onto the velocity data grid using a linear interpolator.
3.1.3. Calculations and output
Taking the masked and smoothed data fields for a given day, the advection and divergence terms can be calculated from Eqn (3). For the intensification terms, multiple time steps need to be considered. The three consecutive days of data held within memory are used to calculate time derivatives using Eqn (2). For consistency, daily advection and divergence are given as a 3 day average. Uncertainty estimates are then propagated for all terms using the method given in appendix A. For a given output time period, mean or summed output data are accumulated as required with data history objects. The uncertainty estimates are accumulated using Eqn (7) using the emergent data covariances presented in appendix B. For certain terms with higher covariance between time steps, 12 days of uncertainties are held in memory to be used within Eqn (7). Upon the writing of output data, accumulation parameters and uncertainty values for covariances are reset.
3.2. Comparisons to previous methods
The main budget equation (1) is a direct development of the concentration budget of Holland and Kwok (Reference Holland and Kwok2012) and Holland and Kimura (Reference Holland and Kimura2016). The use of ice thickness in extension to this approach has been performed by Anheuser and others (Reference Anheuser, Liu and Key2022), though only through the exploration of the advection term on weekly timescales, with a focus on independent thermodynamic winter growth estimates and with an estimated uncertainty for the final data using the thickness uncertainty. Ricker and others (Reference Ricker, Kauker, Schweiger, Hendricks, Zhang and Paul2021) also use gridded winter thickness data, but velocity data are only considered at regional boundaries to estimate the net advection and divergence/convergence-related sea-ice volume changes are not directly estimated. The key methodological developments for this study are the use of thickness within advection and divergence components of Eqn (1), the thorough propagation of uncertainty in all data (appendix A), the use of data covering the full year and careful closure of the volume budget (Section 3.1.1).
4. Results
4.1. Budget components
Within the volume budget, the calculations presented in this paper show that, for an average growth and melt season, the residual or thermodynamic component is the largest, with an average 1.76 m of sea-ice grown and 1.53 m of sea ice melted each year (Fig. 2). The intensification (plots a) and residual (plots b) terms are both positive for growth seasons and negative for melt seasons. These values present an imbalance focussed upon in Section 4.2. The advection (plots c) and divergence (plots d) terms have lower average values that show regional variations. Advection, on average, results in 0.7 m of sea-ice loss from the Laptev sea during the growth season and 0.6 m of sea-ice increase in the Beaufort Sea for example. These regions both correspond to regions of higher average drift speed (arrows in all plots). The divergence component leads to a consistent loss of ice within the Beaufort Sea, with other regions typically variable from year to year with the average state reducing to a low mean value. From comparing the intensification maps (a, total change in ice volume) to the calculated residual maps (b, primarily due to thermodynamics), we can see clear spatial differences. The intensification maps are spatially more even than the residual, with the magnitude of changes in the growth season similar to the melt season. The residual has more spatial variation, with higher growth within the Siberian seas, and higher ice loss north of the Fram Strait and immediately north of Greenland and the Canadian Archipelago.
Average volume budget components using OSISAF data over the period 16/10/2010 to 15/10/2020, for intensification (a), residual (b), advection (c) and the divergence (d), for growth seasons (16 Oct–15 Apr) and melt seasons (15 Apr–16 Oct). Each budget component is given as the total volume change per unit area in meters, along with its propagated uncertainty. The budget components and uncertainties are calculated daily and summed to a seasonal value. The average of all seasonal values are presented here. The arrows indicate the average ice drift speed for the season shown. The Pathfinder version of this figure is in the supplemental material S2.
The propagated uncertainty for the intensification, advection and residual components present spatial characteristics seen within the uncertainty of original ice thickness data. This is particularly evident within the growth season (Fig. 2a and c) where rings of differing orbital overlap in the satellite orbits are apparent, though this pattern is not observed in the corresponding data fields that are free from any artefact. This results in a variable number of individual records for the median gridded ice thickness (Landy and others, Reference Landy2022), with both a less defined mean (fewer data), or wider distribution (more data) affecting the uncertainty in the original thickness data. The divergence term, which has a greater dependency upon sea-ice drift, does not show this ring-like pattern in the uncertainty. While the average magnitude of the intensification and residual terms are similar, the uncertainty of the residual is smaller (seen from comparing plots a and b). This result is discussed further in Section 4.3.
We present the in-depth budget components for the region immediately north of the Fram Strait over the year from 01/10/2014 to 01/10/2015 (Fig. 3). This region is selected as an area of high advection and divergence (maps c and d). The map of the budget components on 19/02/2015 shows ice with a thickness between 2 and 4 m drifting southeastwards towards Svalbard (map a). This day has an average drift speed of 0.25 m s−1, one of the days with the highest average drift speed and Fram Strait ice volume flux over the plotted period (plot f). The high drift speed corresponds to a high level of positive advection (influx of relatively thicker ice) and negative divergence (ice spreading out, with lead formation expected). Note here that these budget components are functions of ice thickness, with high magnitudes when the ice is thicker, not just when the ice drift speed is higher. Spatially, the intensification shows a band of increasing ice thickness at the ice edge near the Svalbard coast (Fig. 3c). The budget components show that this increase can be accounted for by the advection of thicker ice into the MIZ. The advection of ice is balanced by a region of negative divergence, showing the spread of ice, and negative advection north of Greenland, where thicker drifting sea ice is replaced by thinner ice. The resulting residual (Supplemental Fig. S3) has increasing ice thickness for the majority of plotted region and the same high temporal variability as the intensification (Fig. 3g).
Volume budget terms using Pathfinder data in the region immediately north of the Fram Strait in detail for the period 01/09/2014–31/09/2015 with maps of budget components (a)–(d) for the day 19/02/2015 as indicated by the vertical line on the time series. The ice export in (i) is calculated as the volume of ice transported normal to the solid black line on the maps. The time series plots for the budget and export components are units of volume change (km3) for the date plotted. Monthly summed data are shown as circles, scaled by 
$1/30$ for comparison. Dashed lines in plots (f)–(i) are for additional OSISAF drift data. A map of the residual component is in Supplemental Figure S3.
The accessed sea-ice thickness data for the region in Fig. 3a is >2 m for the date shown. These maps (Fig. 3a–d) contain the filled ice terms described in Section 3.1.1 that can be seen as the blue cells (< 0.25 m thickness) in Fig. 3a. The average thickness in Fig. 3e considers the filled thin cells and the plotted variance represents the supplied thickness data uncertainty and not the distribution over Fig. 3a. The filled cells are the source of the peak values of advection (Fig. 3c) as this thin ice is advected into a grid cell of thicker ice (dark red) or vice versa (dark blue). While this particular situation is not realistic for the data points next to the coast, this method ensures all ice is accounted for and the budget is closed. Despite the high signal in Fig. 3c, an opposite drift pattern 2 days later balances these peak advection values and the summed thin ice terms have limited total impact in this region. This is seen in Fig. 2 where the advection has no high values near the coast and the summed values for budget components for the filled cells in Supplemental Fig. S1.
4.2. Time series
Figure 4 shows the full time series of budget components for the full Arctic Basin (Central, Beaufort, Chukchi, Siberian, Laptev, Kara and Barents seas, see Fig. 5). The budget is dominated by the residual term (plot a), with intensification and residual accounting for, on average 
$2000\,\mathrm{km}^3$ of sea ice per month over the growth season, covering all of October–April. The melt season is shown to be shorter with larger losses (
$4000\,\mathrm{km}^3$ per month) in sea ice seen particularly in June and July. The dynamical terms (Advection and Divergence, plot b) along with the Net Transport and MIZ (plot c) are smaller, contributing up to 
$500\,\mathrm{km}^3$ of sea ice per month. Advection and divergence have a tendency to oppose each other on a monthly scale, with an average advected ice loss and divergent gain (compressive-like) of sea ice (plot b, mean seasonal cycle).
Full time series of budget components and ice export for the whole Arctic basin (see Figure 5 for region definitions). The units are the total volume change per month, with the average seasonal cycle shown on the far right. Plot (a) is for the intensification and residual, plot (b) for dynamical terms and plot (c) for the ice transport and ice change within the MIZ. The shaded areas represent the propagated uncertainty from the original data. Dashed lines are for terms using OSISAF data, solid for Pathfinder.
Averaged seasonal budget components for key regions of the Arctic sea-ice system using Pathfinder data for the period 15/10/2010–15/04/2022. A version of this figure for OSISAF data in Supplemental Figure S8. The definitions of growth and melt are 6 month split at 15 Oct and 15 Apr. The units represent the volume change due to the listed components per growth season (bold colours) and melt season (pale colours). Note the different scales on the y axes. For the ice transport, each gate is shown along with the total ice transport. The regions are shown in the bottom right. The Arctic Basin region is equivalent to all regions except for the Baffin Bay and the Greenland Sea.
The ice transport is dominated by the Fram Strait, though this term is the smallest of all the dynamical terms with a maximum ice loss of order 
$100\,\mathrm{km}^3$ per month. The MIZ terms show new ice at the sea-ice edge typically in only October and November, while for ice melt similar levels of ice loss last from April through to August. These levels represent two phenomena: first the greater and more accurate ice thickness data coverage during winter months compared to the summer (MIZ terms include the areas not represented by thickness data, a particular issue in early 2011), and secondly that during the growth season the MIZ is only extensive immediately after the sea ice minimum until the whole basin is again recovered by ice, while in the summer MIZ conditions cover a wider band of the sea-ice edge.
Figure 6 compares the total growth and melt seasons for the Arctic Basin. For this period, there is a small negative trend in 15 October sea-ice volume (estimated for the minimum), that is an order of magnitude smaller than the inter-annual variability. The 15 October and 15 April volume anomalies (estimated date for the maximum) cover a similar range (plot a) and the average sea-ice drift speeds are higher during the growth season. The dynamics terms result in sea-ice loss in both melt and growth seasons, which cause the residual to exceed intensification during growth and to be less than the intensification during melt seasons (note that the sign of the melt season residual and intensification are reversed for comparison). While the total minimum and maximum sea-ice volume anomalies are of a similar magnitude to the changing budget terms (
$\pm 2000\,\mathrm{km}^3$), no clear correlations between them exist over the whole basin in this figure. Maximum sea-ice volume in 2013, 2016 and 2018 (plot a) are local minima and correspond to higher average drift speed (plot b) and lower winter intensification and residual terms (plot c). However, the higher drift speed does not correspond directly to particular dynamical states (plots d and e). Drivers for the seasonal changes to sea-ice volume are discussed further in Section 4.2.
Full time series of budget components and ice export for the whole Arctic basin using Pathfinder data (see Figure 5 for region definitions). The units are the total volume change per growth season (solid lines Pathfinder, dot-dashed OSISAF) and melt season (dashed lines Pathfinder, dotted OSISAF). The definitions of growth and melt are 6 month split at 15 Oct and 15 Apr with max/min volume estimated at these dates. Plots (a) and (b) are the input ice thickness and ice drift speed data plot (c) is the intensification and residual with melt season having the opposite sign to ease comparison, plot (d) is dynamical terms and plot (e) for the ice transport and ice change within the MIZ. Plots (d) and (e) have the signs conserved for both seasons. The shaded areas represent the propagated uncertainty from the original data.
4.3. Uncertainty and signal to noise
Following the estimation of uncertainty over the maps of volume budget components, we provide insight into the time and length scales over which the confidence in the calculated budget components is the strongest. We have collected results, with uncertainty estimates accumulated using the method described in appendix A, over three key time scales: daily (for the 2015 growth/melt season), monthly and the separate full growth and melt seasons. To translate the maps of uncertainty to a representation of uncertainty over a region, we use the method described in appendix A using an averaged signal-to-noise ratio (S/N). The averaged S/N and percentage of the data with an S/N ratio over 3 is shown in Table 1. For this table, the budget components and uncertainties are both accumulated over the given time scale, and the S/N ratio calculated from data arrays similar to those shown in Fig. 2. A higher value for S/N indicates data with greater certainty.
Signal-to-noise ratio for each volume budget component and the percentage of each data that has a signal/noise ratio > 3. Italic numbers show values with a ratio < 1.5, bold for > 4. In general, higher values present a greater certainty for each volume component
In general, longer time scales have a higher signal-to-noise ratio, and thus greater certainty and more accurate physical representation. This can be seen in Table 1 where, particularly for intensification and residual, the average S/N is highest for data covering the entire growth season. In the Central Arctic, and many of the seas of the Arctic Basin, over 90% of our data have a signal-to-noise ratio greater than 3, giving a high level of confidence in our calculations for these regions. For the intensification and residual, the certainty in the data increases up to seasonal time scales due to the sign of these terms being consistent, with positive growth values allowing for the signal to grow higher than the uncertainty for longer periods despite the relatively high emergent time covariance between these data (see Table A1 in appendix B). For the advection and divergence terms, the opposite is at play. Typically, these two dynamical terms have highly variable signals from day to day (see the Fram Strait region in Fig. 3 where these terms are positive and negative). While the uncertainty decreases slightly with longer periods (small increase in S/N in Table 1), the variable signal reduces the total budget components to the mean, with the regions with high and variable ice drift (Greenland and Barents Sea, and Baffin Bay), seeing little change in average S/N between monthly and seasonal time scales. The advection terms are in general less certain than the divergence terms. This is due to the advection, to leading order a function of local ice thickness, having a higher temporal covariance than the divergence (see Table A1 and Eqns (4a) and (7)). The summing or reduction to the mean of the dynamical terms has regional characteristics, for example, for growth seasons in the Greenland Sea, where we expect a consistent export of thicker ice, the S/N ratio for advection has its highest value. This contrasts with both advection and divergence in the Barents Sea and Baffin Bay where highly variable sea-ice thickness and drift keep the S/N ratio low over all time scales.
For all cases, the residual values typically have the highest average S/N despite being derived from the more uncertain dynamical terms. This can be accounted for by the nonlinear nature of uncertainty propagation. By propagating the errors in this way, we can show that a proportion of the uncertainty within the intensification calculations, and thus the thickness data itself, can be accounted for by daily advection and divergence, with the underlying thermodynamic changes to sea ice over an entire growth season captured with greater certainty than local daily changes to sea-ice thickness.
4.4. Regional variations
While the dynamical budget components for the whole Arctic Basin tend to cancel out (Fig. 4), when separating the key regions more complex relationships appear. For the Beaufort/Chukchi, Siberian/Laptev/Kara and Central regions, where winter maximum extent covers the majority of the open ocean surface, on average the residual dominates the volume budget (Fig. 5 top row, further information in Supplemental Figs S2–9). The budget terms for these regions have relatively low uncertainty, including for the dynamical terms. For the regions where the open ocean interacts with the sea-ice edge all year round and the ice extent retreats in summer the budget is more complex (bottom row) and the terms have higher uncertainty. For the Greenland Sea ice is advected into the region and away from the coast all year round with continual melt and loss within the MIZ. For all other regions in Fig. 5, the residual has on average positive values for growth and negative values for melt seasons, while the dynamical terms can be of either sign in either season.
By taking the net contribution of growth and melt over a full year of volume budget, the contribution of sea-ice dynamics is revealed (Fig. 7). While the seasonal cycle of atmospheric temperatures balances out much of sea-ice growth and melt, the dynamical redistribution of sea ice depends on the geography of a region and the resultant circulation patterns. Two notable years within the presented record are 2013 (16/10/2012–15/10/2013) and 2019 (16/10/2018–15/10/2019) shown in Fig. 8. 2013 has an increase in minimum sea-ice volume while 2019 has a reduction (see the Basin row in Fig. 7). The 2013 winter has a notable ice circulation in the Beaufort Sea and Transpolar Drift, while 2019 has notable summer sea-ice circulation (compare arrows in Fig. 8 to the mean state in Fig. 2). The 2013 recovery in minimum volume is shown to be due to reduced melt compared to the mean (indicated by the residual), despite low growth everywhere except the Siberian and Laptev seas. 2019 has both higher growth and melt compared to the mean, but the significant factor here is the drift, with ice losses from divergence and advection in the Beaufort/Chukchi and central regions (see anomalies associated with Fig. 8 in Supplemental Fig. S14). While Li and others (Reference Lin, Pickart, Heorton, Tsamados, Itoh and Kikuchi2021) suggest that the Fram Strait was key to the minimum sea-ice volume in 2011, our results show that sea ice lost within the MIZ or not covered by thickness data is likely to be much higher (Fig. 7). 2015 has a greater Fram Strait export, but this is balanced by higher winter residual (Fig. 6).
Total net ice volume budget contributions for each consecutive growth and melt season using Pathfinder data. Each bar for each is the sum from the previous year 16 Oct through to the listed year’s 15 Oct, with the average cycle shown on the far right. Each bar is for a budget component and the lines indicate the change to the minimum volume anomaly.
Sea-ice volume budget components using OSISAF drift data for the 2013 (first two columns) and 2019 (second two columns) seasons shown in Figure 7. The growth (previous 16 Oct–15 Apr) and melt (16 Apr–15 Oct) are displayed separately. The budget component maps are best compared to the mean sate shown in Figure 2, with anomaly plots shown in Supplemental Figure S14. A version of this figure using Pathfinder data is shown in S15.
Focussing further on the Beaufort/Chukchi region as an example, 2013 shows positive total intensification (positive red bar in Fig. 7), the majority of which was due to advected sea ice (positive cyan bar), that was transported into the region (positive black bar), during winter months. There was also positive thermodynamics where growth was greater than melt, and a small fraction of ice was lost to divergence and within the MIZ. This resulted in a recovery of the minimum sea-ice volume within this region. In contrast, 2019 had a slight decrease in minimum volume in the Beaufort/Chukchi region, despite much higher growth than melt, balanced by significant losses to divergence, advection and within the MIZ.
For perfect data coverage in a linear system, the change to minimum sea-ice volume shown in Fig. 7 will be the integration of intensification over that period. However, due to the incomplete data coverage (see Fig. 1d) and as sea-ice drifts and diverges, no regional change in minimum sea-ice volume can be predicted by our intensification data alone, as seen in the 
$R{^2}$ values in Table 2 (the intensification row). The budget terms in Fig. 7 are only presented for times when all original data are present. The MIZ term attempts to account for missing ice volume by tracking all new, old and filled ice thicknesses (see Section 3.1.1 and Supplemental Fig. S1). The MIZ and intensification terms get closer to reconstructing the minimum volume. The remaining difference we account for from the imperfect numerical representation (Eqns (2) and (3)) of the nonlinear movement of sea ice between original and filled regions.
Correlation and multilinear regressions for inter-annual change in 15 Oct sea-ice volume taking the prior year long values up to 15 Oct as predictors as shown in Figure 7. Here B./C. refers to the Beaufort/Chukchi and S./L./K. to the Siberian/Laptev/Kara regions. Displayed 
$R{^2}$ values for each individual component are calculated separately. The combined values at the bottom of the table are the adjusted 
$R{^2}$ from a multilinear regression for the listed components. The boxed terms, while having limited individual correlation, are statistically significant (p < 0.01) for the Adv./Div./Res./Trans./MIZ model
The Beaufort/Chukchi region has the highest correlation between intensification and changes to minimum volume, the Central the least. To explore the contribution of further terms, a multilinear regression is used with the adjusted 
$R{^2}$ displayed to balance the improvements from a greater volume of predictors (Kotz and others, Reference Kotz, Balakrishnan, Read and Vidakovic2006). Adding in first the transport from the separate regions, and then net contribution of the MIZ (which also includes estimates of missing data), the multilinear model has improved correlation (Int./Trans./MIZ row of Table 2). By separating intensification into the budget components (Adv./Div./Res./Trans./MIZ row), all except the Siberian/Laptev/Kara region show improved correlation. While the individual budget components have limited correlation to changes in minimum volume by themselves (displayed R 2 values), they are found to be statistically significant within a multilinear combination, particularly within the Beaufort/Chukchi seas (boxed values that are significant within the Adv./Div./Res./Trans./MIZ model). The results within Table 2 suggest that for the Arctic Basin, better predicting skill and understanding of minimum sea-ice volume and thus extent can be obtained from better representing thin and marginal ice (MIZ terms), while for the Beaufort and Chukchi seas advected sea ice may be of greater importance, as also suggested by Mallett and others (Reference Mallett2021).
5. Discussion
Gridded satellite observation based data for sea-ice thickness, drift and concentration have been combined to provide a continuous record of sources of Arctic sea-ice volume change over the period October 2010–April 2022. For daily data, the derived uncertainty is high (Fig. 3). This is due to typically low intensification signals over short periods with low ice drift and near-constant ice thickness. Despite the uncertainty at daily time scales, our study has provided estimates of changes to sea-ice volume over 6 month growth and melt seasons with relatively low emergent uncertainty. As the time period increases the signal grows faster than noise due to the constant sign of key budget terms over both growth and melt terms (Figs 4 and 5 and Table 1). For dynamical terms over the entire Arctic basin, due to net balancing of the signed divergence and advection, the signal can remain low with respect to the noise with the contributions of these terms over the entire Arctic Basin still uncertain.
The method presented in this study, primarily through the concept of capturing the full budget and general circulation patterns, intentionally neglects shorter, individual sea-ice floe length processes. This has been performed in order to extract spatial derivatives of sea drift that have no unphysically high values and thus allow for a continuous record of the components of the volume budget. Alternate approaches that include finer scale dynamical details have been performed (Hutter and others, Reference Hutter, Zampieri and Losch2019), though these methods are unlikely to cover the whole Arctic basin over a seasonal time scale. It is currently challenging to capture dynamical events with both satellite-derived drift and thickness data (such as shown in Fig. 3) consistently over short (synoptic) time scales. This discrepancy of floe to basin wide dynamics also presents a challenge for in situ validation. Direct validation of the satellite ice thickness data used for this study has been performed by Landy and others (Reference Landy2022) over the growth and melt seasons. However, there are few unambiguous in situ observations available for the different (thermodynamic and dynamic) components of the sea-ice volume budget, so it is not possible here for us to validate the individual budget terms.
In this study, two sources of gridded Arctic sea-ice drift data have been used. The plotted maps of the averaged budget terms (Fig. 2) are for OSISAF data. For the same plot using Pathfinder data (Supplemental Fig. S2), the divergence displays a strong positive and negative signal about the pole hole at 88∘ N (growth season signal, plot S2d and S15), likely due to the numerical incorporation of buoy drift data (Szanyi and others, Reference Szanyi, Lukovich, Barber and Haller2016). This signal is an artefact of the Pathfinder drift data not present when using the OSISAF drift. The artefact makes minor differences to the total summed divergence, seen in Supplemental Figs S6 and S11, as the difference between the Pathfinder and OSISAF based budgets can be accounted for from differing drift speeds. OSISAF data have a faster drift speed than Pathfinder for all regions in growth (3.9 vs 3.2 cm s−1) and melt (3.6 vs 2.2 cm s−1) seasons on average and annually (see Fig. 6 and S9–12), with an increased difference for seasonal seas, seen for the Barents sea in Fig. S9 where the OSISAF ice drift speed (6.7 cm s−1 growth season average) is beyond the estimated uncertainty of the Pathfinder drift speed (3.7 cm s−1 growth season average). The difference in drift speed results in a difference in total divergence (910 vs 
$210\,\mathrm{km}^3$ average basin total yearly ice loss) but limited effect on advection (1420 vs 
$1490\,\mathrm{km}^3$). Despite these differences, the key results shown in Figs 5 and 7 are consistent between data sources (Supplemental Figs S8 and S13 for OSISAF).
Over a full growth season at the Arctic basin scale, dynamical terms balance out, with total sea-ice volume changes predominantly thermodynamic (Figs 4 and 5). The remaining portion of the basin-scale volume budget is accounted for by sea-ice advection out of the Arctic via the Fram Strait. For the melt season, export of sea ice through gates like Fram Strait is lower and volume loss is primarily due to melt, though additionally there are losses within the MIZ. When comparing to Keen and others (Reference Keen2021), converting from values of Gt and removing the pole hole, our average growth residual of 
$10.2\,\pm\,2.8 \times 10^6\,\mathrm{km}^3$ is within the range of modelled thermodynamics from CMIP6 models (8.0–
$13.8\times 10^6\,\mathrm{km}^3$). For melt residual, our values 
$7.3 \pm 2.7 \times 10^6\,\mathrm{km}^3$ are lower than CMIP6 (8.0–9.5
$\,\times 10^6\,\,\mathrm{km}^3$), although the MIZ accounts for 1.0–2.0
$\,\times\,10^6\,\mathrm{km}^3$ of sea ice that will be captured as melt within a climate model. Our dynamics related volume changes 
$1.7\pm0.7 \times 10^6\,\mathrm{km}^3$ are at the low end of the CMIP6 range (2.1–3.7
$\,\times\,10^6\,\mathrm{km}^3$). Our summed regional ice volume budget components are comparable to the mean monthly rates given by Ricker and others (Reference Ricker, Kauker, Schweiger, Hendricks, Zhang and Paul2021): 
$470\,\mathrm{km}^3$ of intensification per month within the Beaufort/Chukchi seas during the growth season, 
$470\,\mathrm{km}^3$ in the Siberian/Laptev/Kara. The sea-ice transport rates are also similar (
$15\,\mathrm{km}^3$, 
$180\,\mathrm{km}^3$ respectively). However their method does not close the volume budget of sea ice, and gridded sea-ice velocity is not considered. The calculations of Anheuser and others (Reference Anheuser, Liu and Key2022) are conceptually similar to those presented here, with gridded sea-ice drift velocity and maps of advection that have very similar structures to ours within the Beaufort Sea. However, these are given as regionally averaged rates of change in thickness that are difficult to compare to our complete volume budgets. The MIZ losses in our data represent both the reduced data coverage in summer months at low ice concentration, as well as the complex physical processes within these regions that cannot be represented by these data. By considering the volume changes of sea ice within the MIZ regions, we can close the full volume budget and reveal the whole seasonal cycle in ice volume.
To fully reveal the role of dynamical terms, we take the net contributions of the volume budget over a whole year (Fig. 7). Here, we define the term over-production as the net yearly summed budget residual, capturing the imbalance of sea-ice growth and melt. By contrasting the total closed volume budget with changes in minimum sea-ice volume, key regional predictors, data gaps and nonlinearities are revealed (Table 2). The uncertainty of these whole season budgets will be dependent upon covariances between regional growth and melt characteristics, which is currently not represented by the technique documented in appendices A and B. This combination of budget terms reveals the differing sea-ice conditions over the regions of the Arctic. The Siberian/Laptev/Kara regions are shown to have a stable over-production of 1000–
$2000\,\mathrm{km}^3$ of sea ice per year, an example of this region as a sea-ice ‘factory’ as described by Cornish and others (Reference Cornish, Johnson, Mallett, Dörr, Kostov and Richards2022). The magnitude of over-production is significant to the total over the entire Basin, thus showing the importance of this region for the survivability of Arctic sea ice (Stroeve and others, Reference Stroeve, Maslanik, Serreze, Rigor, Meier and Fowler2011). Sea ice from the Siberian/Kara/Laptev sea is typically transported to the Central Arctic, though due to the variability of Basin circulation patterns (Mallett and others, Reference Mallett2021), there is a high inter-annual variability for this region. Figure 7 shows that decreasing minimum sea-ice volume can be due to both thermodynamics and dynamics, but recovery years (2013 and 2017) are typified by an over-production of sea ice. For the regions that contain no sea ice in summer (Barents and Greenland seas, Baffin Bay), the seasonal cycle has been stable over the study period 2010–22, though each region shows a different balance. The Greenland Sea has high total advection from Fram Strait, with this sea ice lost to melt and within the MIZ. Baffin Bay is an isolated system (save for a relatively small sea-ice volume transport through the Nares Strait that we cannot resolve with the ice drift data used here (Agnew and others, Reference Agnew, Lambe and Long2008)) that balances growth with MIZ losses.
6. Conclusions
Observational satellite data over the Arctic can now be combined to produce a continuous record of sea-ice volume and drift between November 2010 and April 2022. A budget closure approach can be used to track whether changes to the sea-ice volume come from sea-ice advection, divergence or residual thermodynamic growth and melt. The closed seasonal ice volume budget includes the contributions of individual daily ice dynamic events (Fig. 3) and longer timescale cumulative thermodynamic growth (Section 4.2). These contributions can be combined from separate sea-ice concentration, thickness and draft datasets to produce seasonal results despite high uncertainties at the raw daily time interval (Section 4.3).
By decomposing the spatiotemporal variations of dynamic and thermodynamic processes controlling the sea-ice volume budget, the regional drivers of Arctic sea-ice volume inter-annual variability can be revealed (Section 4.2). Due to the enclosed geography of the Arctic, dynamic sources of sea-ice volume change tend to sum to zero over longer timescales across the full Arctic basin, with pan-Arctic ice volume tendency dominated by thermodynamics. Even if the pan-Arctic sea-ice divergence is anomalously strong for a particular month, it is typically balanced by anomalously large ice advection (Section 4.2). To reveal the importance of sea-ice drift, the regional variations must be considered, along with the net yearly budget residual, defined as over-production. Within particular sub-regions, the characteristic split between dynamical and thermodynamical sources of sea-ice volume change is often consistent between years, for both growth and melt seasons (Section 4.2). The Siberian, Laptev and Kara Seas are shown to consistently over-produce sea ice, with the volume of winter growth exceeding summer melt. Thermodynamic over-production occurs in the Central Arctic too (exceeding the Siberian, Laptev and Kara Seas in 2013) but varies strongly year-to-year, with net ice loss in 2011. There is a clear relationship between over-production and the pan-Arctic sea-ice volume anomaly in the following September (Fig. 7 and Table 2). The ability of the Siberian, Laptev and Kara seas, and Central Arctic to over-produce sea ice will be a key factor in the survivability of Arctic sea ice.
Our observation-based calculations of the sea-ice volume budget present an opportunity for detailed model validation of individual mass-balance components rather than Arctic wide metrics of sea-ice extent and volume (Keen and others, Reference Keen2021). For instance, can models reproduce the enhanced sea ice production over winter 2018–19 in the Beaufort and Chukchi seas associated with anomalously high ice divergence? Or the large inter-annual swings in thermodynamic ice growth in the Central Arctic? To understand how the Arctic sea-ice volume budget, and thus summer sea-ice extent, will evolve in the future, it is crucial that models can reproduce the dynamic and thermodynamic processes controlling inter-annual changes in sea ice currently observed.
Supplementary material
The supplementary material for this article can be found at https://doi.org/.10.1017/aog.2025.3.
Appendix A. Propagation of uncertainty
In order to translate the uncertainty estimates supplied as part of the sea-ice thickness and drift data products to an uncertainty estimate for our final sea-ice volume budget calculations, the uncertainty is propagated using the following formula for a function f of A and B, with variances 
$\sigma _{A}, \sigma_{B}$, and covariance between them of 
$\sigma_{AB} = \rho_{AB}\sigma _{A}\sigma_{B}$ with ρAB the emergent Pearson correlation between A and B (Kotz and others, Reference Kotz, Balakrishnan, Read and Vidakovic2006). The uncertainty of f can thus be expressed using:
\begin{equation}
f=A+B \longrightarrow \sigma _{f}^{2}=\sigma _{A}^{2}+\sigma _{B}^{2}+2\sigma_{AB}
\end{equation}
\begin{equation}
f=A-B \longrightarrow \sigma _{f}^{2}=\sigma _{A}^{2}+\sigma _{B}^{2}-2\sigma_{AB}
\end{equation}
\begin{equation}
f = AB \longrightarrow \sigma _{f}^{2}\approx f^{2}\left[\left({\frac {\sigma _{A}}{A}}\right)^{2}+\left({\frac {\sigma _{B}}{B}}\right)^{2}+2{\frac {\sigma_{AB}}{AB}}\right].
\end{equation}For the case of the budget equations, we have uncertainty estimates for the velocity and thickness data, and we take the uncertainty of the per pixel volume, often called the absolute thickness, as the uncertainty of the thickness measurement. No uncertainty is provided for the ice concentration data. We are thus able to propagate estimates of uncertainty for all variables in the budget equations.
For the derivatives key to the budget equation, 
$dV/dt, \nabla V, \nabla\cdot\mathbf{u} $, the calculations are performed numerically over time and spatially adjacent data points. The variance propagations are performed similarly with:
\begin{equation}
\text{Int}_t=\frac{1}{\Delta t}(V_{t+1}-V_{t-1}) \longrightarrow \sigma_{\text{Int}_{t}}^{2}\Delta t^2=\sigma _{V_{t+1}}^{2}+\sigma _{V_{t-1}}^{2} - 2\sigma_{V_{(t+1)(t-1)}}
\end{equation}
\begin{equation}
V_{x}=\frac{1}{\Delta x_i}(V_{i+1}-V_{i-1}) \longrightarrow \sigma_{V_{x}}^{2}\Delta x_i^2=\sigma _{V_{i+1}}^{2}+\sigma _{V_{i-1}}^{2} - 2\sigma_{V_{(i+1)(i-1)}}
\end{equation}
\begin{equation}
u_{x}=\frac{1}{\Delta x_i}(u_{i+1}-u_{i-1}) \longrightarrow \sigma_{u_{x}}^{2}\Delta x_i^2=\sigma _{u_{i+1}}^{2} + \sigma _{u_{i-1}}^{2} - 2\sigma_{u_{(i+1)(i-1)}}
\end{equation}where 
$\sigma_{V_{(t+1)(t-1)}}$ is dependent upon the emergent correlation 
$\rho_{V_{(t+1)(t-1)}}$, the correlation between sea-ice volume pixels at time points t − 1 and t + 1, and similarly for other covariances. These covariances are discussed in the following appendix. Here, we describe 
$\frac{\partial V}{\partial x}$ as Vx, the x component of the volume gradient and 
$\frac{\partial u}{\partial x}$ as ux as the x component of the velocity divergence. The y and v components are not shown fully here but follow similarly with calculations over the j index.
Thus, the uncertainty for the final dynamic budget terms (ignoring all zero valued emergent covariances) can be calculated with:
\begin{equation}
u\cdot V_{x} \longrightarrow \sigma^2_{u\cdot V_x} = u\cdot V_x\left[\left(\frac{\sigma_{u}}{u}\right)^2 + \left(\frac{\sigma_{V_{x}}}{V_{x}}\right)^{2} \right]
\end{equation}
\begin{equation}
\text{Adv} = u\cdot V_x + v\cdot V_y \longrightarrow \sigma^2_{\text{Adv}} = \sigma_{u\cdot V_{x}}^{2} + \sigma_{u\cdot V_{y}}^{2} +2\sigma_{u\cdot V_{x},u\cdot V_{y}}
\end{equation}
\begin{equation}
\nabla\cdot\mathbf{u} = u_x + v_y \longrightarrow \sigma_{\nabla\cdot\mathbf{u}}^2 = \sigma_{u_x}^2 + \sigma_{v_y}^2 + 2\sigma_{u_x,v_y}
\end{equation}
\begin{equation}
\text{Div} = V\nabla\cdot\mathbf{u} \longrightarrow \sigma_{\text{Div}}^2 = \text{Div}^2\left[ \left(\frac{\sigma_V }{V } \right)^2 + \left(\frac{\sigma_{\nabla\cdot\mathbf{u}}}{{\nabla\cdot\mathbf{u}}} \right)^2 \right]
\end{equation}where Adv. and Div. are the Advection and Divergence components of Eqn (1). Equations (6b) and (6c) having emergent covariances presented in Table A1.
To provide an uncertainty estimate for the residual for each day, we use Eqns (4a) and (4b), considering the correlation between Int. and Adv. (0.08), with no significant correlation between Int. and Div. found. To combine the daily estimates of uncertainty in intensification, advection, divergence and residual to give monthly or seasonal values, covariances between consecutive days were found to be correlated, particularly when the term was dependent on ice volume, that has slow, near constant, rates of change for the 15 day data used here. As the budget code accumulates the mean values from day to day, it was possible to similarly accumulate the uncertainty estimate with the covariance of the previous 
$t_c = 12$ days. The summed Int. uncertainty for ts days (with Adv. and Res. calculated similarly) is thus:
\begin{equation}
\sigma^2_{\text{Int}_{t_s}} = \sum_{t = 1}^{t_s} \left[\sigma^2_{\text{Int}_{t}} + \prod_{m=1}^{t_c} 2\rho_{\text{Int}}^m\sigma_{\text{Int}_{t-m}}\sigma_{\text{Int}_{t}}\right]
\end{equation}with the value of 
$\rho_{\text{Int}}$ a constant parameter taken from the emergent covariance raised to a power m.
The uncertainty propagation gives maps of uncertainty for each budget component. In order to total the volume of ice change over the various regions and time scales presented in this paper, a method of estimating the uncertainty for these sums is required. Adding the σ values either linearly or in quadrature results in values that are dependent on the grid resolution. To remove this dependency, we first calculate the mean signal-to-noise ratio 
${S}_A$ over the spatial region of interest for a given quantity A and its uncertainty σA. We then use this mean ratio to estimate a noise or uncertainty from the final sum values as follows:
\begin{align}
{S}_A & = \frac{1}{n} \sum_{i = 1}^{n} \frac{A_n}{\sigma_{A_n}} \nonumber\\
\sigma_{\Sigma A} &\approx \frac{\Sigma A}{{S}_A}
\end{align}where the index i is over the n grid cells within the map of quantity A. This method allows for the simple estimation of uncertainty independent of grid resolution without the complex task of considering spatial covariances across a 2-d grid. However, it has a drawback in that the magnitude of the uncertainty is proportional to the signal strength 
$\Sigma A$. This is most apparent for values that are close to zero such as the divergence terms in Fig. 4b.
Appendix B. Key covariances within volume budget data
To provide estimates of correlation between key values for the propagation of uncertainty, constant scalar values were sought. For future work, it would be possible to expand this methodology to give maps or seasonal cycles in correlation. While these would be interesting in their own right, the benefits of such complications in estimating uncertainty are unclear and beyond the scope of this study. While seasonal cycles in uncertainty are represented by the original data uncertainties, a benefit of extending the constant values used here to spatially varying values is that any data variation at low latitudes where the orbital overlap is low may be able to emerge from this analysis of covariance. Here, we focus on first finding the order of magnitude of key correlations. The correlations were found by analysing a single growth and melt season from 01/10/2014 to 31/9/2015. Each pair of data were analysed to find the Pearson correlation parameter, with the median value for the year used. This was performed for both raw input data and the smoothed data, with spatial correlations increasing for smoothed data. The values in Table A1 summarize the emergent key correlations found and are then used to parameterize the propagation of uncertainty estimates.
Key covariances for volume budget calculations. Correlations are calculated for each day of data over the period 01/10/2014–31/9/2015, with the median covariance value shown here. Covariances were calculated for the raw input and for smoothed data using both Pathfinder and OSISAF velocity data and grids. The number of decimal places for each term represents the accuracy obtainable using current methods
Four data scenarios are presented in Table A1: Pathfinder and OSISAF drift data and grids, with and without smoothing. For the initial spatial derivatives, 
$\mathrm{d}V/\mathrm{d}x_i$ and 
$\mathrm{d}u/\mathrm{d}x_i$, the application of smoothing increases the covariances for both data sources and thus reduces the uncertainty on these values. For later terms, such as for the addition of divergence components 
${\mathrm{d}u}/{\mathrm{d}x_i}$ and 
${\mathrm{d}v}/{\mathrm{d}x_j}$, the covariances are negative and thus increase the uncertainty compared to a random variable. The covariance becomes greater in magnitude after the application of smoothing. For the case of divergence, this suggests that the smoothed data have more ‘incompressible like’ behaviour than the unsmoothed data. An incompressible flow will have perfectly anti-correlated components. The OSISAF drift and other data on the OSISAF grid all have lesser covariances. This is due to the wider gird spacing of the OSISAF grid compared to the Pathfinder grid.
The key covariances are the initial partial derivatives in Eqn (5a) that are the difference of variables of similar magnitude. If these variables are randomly correlated then the uncertainty will greatly exceed the signal. However, these variables are found to be highly correlated for neighbouring pixels and consecutive time steps, see Table A1, and thus the final term in Eqn (5a) is significant.
Open access
