Seasonality and timing of sea ice mass balance and heat fluxes in the Arctic transpolar drift during 2019-2020

,


Introduction
The observed decline in Arctic sea ice strongly modifies the exchanges of energy, momentum and matter between the atmosphere and the ocean, thereby significantly affecting the Arctic climate and marine ecosystem.Thermodynamic mass balance of sea ice plays a governing role in these exchange processes, especially on a seasonal scale.The most important processes determining the sea ice thermodynamic mass balance include initial freezing, ice basal growth, and formation of snow ice or superimposed ice at the snow-ice interface, as well as surface, internal, lateral and basal melt (Perovich et al., 2021).Furthermore, snow plays crucial roles in the sea ice mass balance either passively, by slowing down heat release to the atmosphere during winter (Sturm et al., 2002) and buffering sea ice melt in summer (Perovich et al., 2002), or actively, by refreezing at the surface after snow melt (superimposed ice) and flooding (snow ice) (Nicolaus et al., 2003;Merkouriadi et al., 2017).The continued thinning of Arctic sea ice and the widespread occurrence of negative ice freeboard are expected to further enhance the contributions of snow to the sea ice mass balance (Ro ¨sel et al., 2018).
The main challenge to understanding changes in Arctic sea ice mass balance is to accurately quantify sea ice thickness, and particularly its seasonal evolution.Satellite altimeters are capable of determining sea ice thickness on a large scale; however, the accuracy and spatial resolution of the results are limited (Sallila et al., 2019).Measurements based on acoustic or electromagnetic induction soundings collected by submarines, underwater robots, aircrafts, or ships cannot be used to distinguish temporal changes in snow depth and sea ice thickness.In situ manual observations, such as drillings, hot-wire gauges and ablation stakes, are the most accurate methods to quantify snow depth and ice thickness (e.g., Ro ¨sel et al., 2018), but the ability to maintain such measurements on a regular basis is usually very limited due to logistical challenges.A solution for this problem is the use of the autonomous sea ice mass balance buoys (IMBs).The most prominent examples are the IMB and Seasonal IMB designed by the Cold Regions Research and Engineering Laboratory (referred to as CRREL IMB or SIMB) (Richter-Menge et al., 2006;Planck et al., 2019).By combining acoustic sensors with a thermistor chain, these instruments are able to observe simultaneously the changes in snow depth and sea ice thickness as well as the vertical temperature profile.More recently, the snow and ice mass balance array (SIMBA) designed by the Scottish Association for Marine Science Research Services Ltd, Scotland (Jackson et al., 2013) was deployed in significant numbers in the polar regions.SIMBA is a thermistor string-type IMB (Jackson et al., 2013) which measures the environmental temperature (SIMBA-ET) and additionally determines the temperature change around the thermistors after a weak heating applied to each sensor (SIMBA-HT).The combination of SIMBA-ET and SIMBA-HT can yield snow depth and ice thickness, as well as the formation of snow ice or superimposed ice (Provost et al., 2017).Although an individual IMB only collects point measurements, interpretation of such data can still be upscaled, with some limitations that depend mainly on deployment locations and ice types, to represent mesoscale sea ice thermodynamic processes (Perovich and Richter-Menge, 2006).Moreover, clustered IMB observations can further enhance the spatial representativeness (e.g., Provost et al., 2017).
Arctic sea ice circulation is driven mainly by the Beaufort Gyre (BG) and the Transpolar Drift (TPD).In the BG region, the Surface Heat Budget of the Arctic Ocean (SHEBA) campaign in 1997-1998 provided the most comprehensive measurements of sea ice mass balance in the Arctic Ocean to date (Perovich et al., 2003).However, clustered deployments of IMB arrays have been relatively sparse in the TPD region.For example, such an array was deployed further downstream in the TPD during the N-ICE campaign in 2015 (Provost et al., 2017), but those buoys were deployed in the spring and mainly covered the ice melt season.In the last 20 years, sea ice in the regions close to Fram Strait, the main outflow gate of Arctic sea ice, has become thinner (Spreen et al., 2020;Belter et al., 2021), which is related to the influence of Atlantification on the sea ice growth in the TPD upstream region (Polyakov et al., 2017), and the lack of sufficient time for sea ice thickening due to the acceleration of ice advection in the TPD (Spreen et al., 2011;Krumpen et al., 2019;Belter et al., 2021).To assess the contribution of sea ice outflow through Fram Strait to the overall Arctic sea ice loss, obtaining observations of the sea ice thermodynamic processes during the drift with the TPD is of particular importance.
In 2019-2020, the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) expedition was implemented, during which comprehensive, yearround observations of the Arctic sea ice and its interactions with the atmosphere and ocean were obtained (Nicolaus et al., 2022;Rabe et al., 2022;Shupe et al., 2022).The MOSAiC ice camp drifted from the marginal ice zone (MIZ) in the eastern Amundsen Basin in early October 2019 to Fram Strait by early August 2020.In addition to the Central Observatory (CO) installed on an ice floe next to the R/V Polarstern, an extensive Distributed Network (DN) of buoys was deployed in a radius of 40-50 km around the ship (Krumpen and Sokolov, 2020), which incorporated a subarray of IMBs.
In this study, we present results from 23 IMBs to provide a comprehensive characterization of the sea ice mass balance of the wider MOSAiC DN area from October 2019 to July 2020.We deduce the amount, seasonality and timing of thermodynamic growth and decay for the floes with different initial conditions within the MOSAiC DN.In addition, the effects of the snow cover and dynamic processes on the sea ice mass balance are analyzed.Finally, we use auxiliary atmosphere and ocean datasets to determine the main drivers that governed sea ice thermodynamic evolution during the study period.

Instruments and deployment setup
In this study, we used the data collected by two types of buoys: a large number of SIMBAs and one multi-sensor unmanned ice station (UIS).The SIMBA thermistor chain is 4.80 m long and equipped with 241 thermistors (Maxim Integrated DS28EA00) at 0.02-m spacing (Jackson et al., 2013).SIMBA-ET was recorded at a temporal interval of 6 h.Moreover, SIMBA-HT was observed daily after heating cycles of 30s (dT 1 ) and 120s (dT 2 ).
The UIS is a new prototype of IMB assembled by the Polar Research Institute of China, which consists of two separate units (ice and ocean) to measure physical and biogeochemical parameters of the air-snow-sea ice-ocean system (Figure 1).In this study, UIS data were mainly used to characterize physical ice-ocean interactions.For the ice unit, two acoustic sensors were used to measure the relative changes in the positions of the air-snow and ice-water interfaces.Thermistors mounted at 0.03-m spacing along a 4.5-m thermistor chain were used to measure temperature profiles.Air temperature and relative humidity (Vaisala HMP155A), as well as barometric pressure (Vaisala CS106), were measured at 1.5-m height above the initial snow surface.The spectral albedo and transmittance of the ice were measured by optical sensors (not used here).The UIS ocean unit consisted of five conductivity and temperature sensors (RBR duo CT), one conductivity, temperature and depth (pressure) sensor (RBR concerto CTD), and additional sensors measuring selected biogeochemical parameters (not used here).The ocean unit measured upper ocean properties at the depths of about 5-40 m, at the initial depths of 5.4, 10.4, 15.4, 20.4, 25.4, and 40.4 m, which typically exhibit a strong seasonal signal compared with the deeper layer (e.g., Athanase et al., 2019).
In total, 29 SIMBAs and one UIS were deployed during the MOSAiC campaign.In this study, we used the 23 instruments that were deployed at least one month before the onset of sea ice melt in summer 2020, limiting the data to focus on the ice season of 2019-2020 (Figure 2).All buoys used in this study were deployed over initial unponded level ice to ensure the representativeness of the observations with regards to upscaling (e.g., Perovich and Richter-Menge, 2006).
In addition to the CO, the MOSAiC DN included three types of deployment sites: "Large" (L-) sites with an extensive suite of installations; Medium (M-) sites with a handful of different platforms; and Position (P-) nodes typically with only a single GPS buoy (Krumpen and Sokolov, 2020).In the MOSAiC DN area, the sea ice was dominated by first-year ice and second-year ice with diverse initial ice thickness (Krumpen et al., 2020), which provided an excellent opportunity for studying the influence of initial ice conditions on the subsequent sea ice growth processes.Fourteen buoys deployed in the DN setup phase from October 4 to 16, 2019, are used here, including one SIMBA at L1, two SIMBAs at L2, one SIMBA and one UIS at L3 (Figure 1), one SIMBA at each of the 8 M-sites, and an additional SIMBA at one P-site (Figure 2a).Three additional SIMBAs were deployed in the CO (approximately 2 km from the ship) in late October and early November 2019 to complement other measurements of ice stress buoys (IS site) or the coring sites on the first-year (FYI) and second-year ice (SYI).Several buoys were destroyed by floe fracturing or ridge formation after the deployments, but 11 buoys provided valid observations that covered the entire ice growth period (Figure 2b).
The average (± standard deviation) initial ice thickness at the thermistor chain deployment sites was 0.96 ± 0.48 m (n ¼ 13), with minimum and maximum thicknesses of 0.35 and 1.80 m.Comparing to Krumpen et al. (2020), whose on-ice surveys showed an average ice thickness between 0.67 ± 0.54 m (obtained from the 9.6-km survey line) and 1.0 ± 0.81 m (obtained from the 7.9-km survey line) at the L1, L2, L3, and M8 sites, the IMB deployment sites seem representative of the general DN ice conditions.The snow depth at the various sites did not show large differ e n c e s ,w i t ha na v e r a g eo f 0.11 ± 0.04 m (n ¼ 13) at the time of deployments.By December 1, 2019, the buoys were distributed within a distance of approximately 50 km at 85.8 -86.2 Nand 109.5 -116.5 E( Figure 2a).
Six additional SIMBAs were deployed at the CO or the floes (additional P sites) within 2-3 km from the ship between April 4 and 23, 2020, during Leg 3 to supplement the observation sites.By July 1, 2020, the buoy array had drifted over the northwestern Yermak Plateau at 81.53 -82.01 N and 7.7-10.3E( Figure 2a).Most buoys failed or were recovered by the end of July 2020 close to the ice edge in southern Fram Strait.The measurements of the UIS ocean unit lasted until September 28, 2020, because this unit was able to float after the ice melted, and finally failed at 74.09 N( Figure 2b).

Determinations of air-snow-ice-water interfaces
The physical basis for identifying the interfaces among air, snow, sea ice, and water using the ET vertical gradient or the HT rise ratio, defined as the ratio of dT 1 to dT 2 , relies upon the differences in thermal conductivities and thermal capacities for the different materials (Provost et al., 2017).Here, we combined the two methods of the ET vertical gradient and HT rise ratio to reduce the uncertainties because they have respective advantages in identifying different interfaces.As shown in Figure 3,t h e snow or ice (when the snow melts completely) surfaces identified using these two methods are comparable for both thick ice and thin ice.To keep the results consistent, we used the snow or ice surface determined by the ET vertical gradient for the following analysis because of its higher sampling frequency compared to the HT measurement.
During the initial stage after the deployment, however, the freezing front had not reached the ice bottom for the relatively thick ice (Figure 3a).Thus, the ET vertical gradient during this stage was not suitable to identify the ice bottom.Similarly, during the melt season, the appearance of meltwater under the ice resulted in an indistinct icewater interface in the SIMBA-ET profiles (Figure 3a).Thus, we used the HT rise ratio to determine the ice-water interface.The appearance of freshwater under the ice, which originates from snow and ice meltwater, has been termed "under-ice melt pond" (Eicken, 1994).The under-ice meltwater can be refrozen to create a false ice bottom because the thermal diffusion is much faster than salinity diffusion (Notz et al., 2003).However, from both measurements of SIMBA-ET and SIMBA-HT, we are still unable to judge whether the under-ice meltwater was frozen or not.Therefore, as soon as the under-ice melt ponds formed, we used the HT rise ratio to identify the total thickness of the under-ice meltwater plus the false ice bottom without further discrimination.Thus, this thickness was defined as the distance from ice bottom to the lower boundary of under-ice melt pond.Moreover, the HT rise ratio was also used to identify the interface of normal snow and snow ice or superimposed ice when they formed atop the ice layer, as shown in Figure 3b and d.
At the UIS site, the ice bottom evolution was determined using measurements by an underwater acoustic sounder.As the acoustic sensor measurements at the surface was invalid very soon after the deployment, the evolution of the air-snow interface could only be determined using the vertical temperature gradient similar to the SIM-BAs.Overall, the measurement accuracy was 0.1 Cf o r temperature, 0.02 and 0.03 m for the snow or ice surface, and 0.02 and 0.01 m for the ice bottom for the SIMBA and UIS, respectively.Surface melt includes the melt of both sea ice and snow, hence we estimated the mass equivalent surface melt as 1/3h s þh i by assuming that snow has a density of one third of that of sea ice (e.g., Perovich et al., 2014), where h s and h i are the melt thicknesses of snow and ice surface, respectively.Note that, although a surface melt pond can be identified unambiguously from the profiles of the ET or HT rise ratio (Figure 3a-d), the bottom of the surface melt pond was often indistinct, which was possibly related to the high porosity of the sea ice under the pond.Moreover, the melt pond surface is also difficult to determine because of repeated freeze-thaw processes.Thus, the evolution of surface melt ponds was not quantified in this study.Table 1 summarizes the methods used here to identify the interfaces between air, snow, sea ice, and water based on the data recorded by the two types of buoys.Note that no snow ice, superimposed ice or underice melt pond formation was observed in the UIS data.Thus, these processes are only discussed for the SIMBA sites.

Estimations of heat fluxes through the ice
The conductive heat fluxes within the ice column through every 0.10 m below (above) a distance of 0.05 m from the ice surface (bottom) were calculated to analyze the heatcontent regulation efficiency of brine within the ice.The heat flux through the snow and the top 0.05-m ice layer were not estimated here because this layer may involve complex texture compositions, which makes determining the thermal conductivity coefficient difficult (Sturm et al., 2002).Similarly, we defined the lowest ice layer as 0.05-0.15m from the ice bottom to exclude the basal unconsolidated skeletal layer.To calculate the thermal conductivity coefficient of sea ice, we used the ice temperature measured by the thermistor chain and an ice salinity estimated from the ice coring measurements at the CO SYI site.However, because the growth history of sea ice usually has a large variability among different sites, we simplified the observational results by assuming that the salinity for the ice bottom layer (0.05-0.15 m) was 6, and the above layers was 4.Then, the thermal conductivity flux (F c )w a sc a l c ulated using sea ice thermal conductivity coefficients and vertical temperature gradients (Lei et al., 2018).
At the ice bottom, freezing and melting results from a combination of the conductive heat flux (F c ), the equivalent latent heat flux (F l ), the sensible heat flux (F s ), and the oceanic heat flux (F w ).One can estimate the oceanic heat flux (F w ) as the residual term among these heat fluxes (e.g., McPhee and Untersteiner, 1982;Lei et al., 2018).To calculate the equivalent latent heat flux, we assumed aconstant salinity (8) and a constant density (910 kg m -3 )for the basal newly formed ice.The freezing point temperature was time-dependent and determined using the UIS measurement at a water depth of about 5.4 m.

Atmospheric and sea ice climatological data
On the basis of atmospheric reanalysis data and remote sensing products, we compared the sea level air pressure (SLP), near-surface air temperature (2 m, T 2m ), wind speed (10 m, W 10 m ), ice motion speed, ice thickness, and ice concentration obtained along the buoy trajectory from the study year (2019-2020) to the climatology to identify the anomalies of atmospheric and sea ice conditions, which are expected related with sea ice mass balance processes.We calculated the Central Arctic Index (CAI), defined as the difference in SLP between 90 Wand90 Eat84 N (Vihma et al., 2012), to describe the anomalies of atmospheric circulation patterns and characterize the TPD intensity.
Atmospheric reanalysis data, including W 10 m , T 2m ,an d SLP, obtained from the ERA-5 dataset of the European Centre for Medium-Range Weather Forecasts (ECMWF), were used to determine the meteorological anomalies along the buoy trajectory.The daily products of the Sea Ice Motion Vectors Version 4 dataset (Tschudi et al., 2020) and ice concentration data from the Nimbus-7 Scanning Multichannel Microwave Radiometer (SMMR) and its successors (SSM/I and SSMIS) (Fetterer et al., 2017) provided by the National Snow and Ice Data Center (NSIDC) were used to estimate anomalies of ice motion speed and concentration.The weekly merged CryoSat-2-SMOS sea ice thickness product provided by the Alfred Wegener Institute (Ricker et al., 2017) was used to estimate anomalies of ice thickness.All data, except for ice thickness, were interpolated linearly to the position of the MOSAiC CO according to the day of the year.The sea ice thickness product was averaged over an area with a radius of 50 km from the MOSAiC CO.

General atmospheric and sea ice conditions
The reanalyzed SLP, T 2m , and W 10 m and remote sensing product of ice motion speed matched well with the in-situ observations obtained from the UIS (Figure 4a, b, and e) or ship-based measurements (Figure 4d; Knust, 2017), which enhances our confidence in using them to identify anomalies in the study year.There were 60 d with daily SLP lower than the 1979-2020 climatology by one standard deviation during the period from October 9, 2019, to July 18, 2020.This number of days was greater than the 41 ± 18 d obtained in the same period for the 1979-2020 climatology.The low pressure systems also resulted in increases in both air temperature and wind speed (Figure 4b and d).However, along the buoy trajectory, the wind from the northwest was dominant (Figure 4e), which might have brought cold air masses to the buoy location.Thus, the combined effect led to the average T 2m from October 9, 2019, to July 18, 2020 (-15.7 C), being very close to the 1979-2020 climatology (-15.8C), and consequently the comparable cumulative freezing degree days (Figure 4c).In contrast, in the warm season from mid-May to July, the average air temperature was 0.4 C, which was higher than the 1979-2020 climatology by 1.0 C.This higher average was associated with the earliest upward crossing (May 24, 2020) of the 14-day running average air temperature through the melt threshold of -1 C (Rigor et al., 2000) of 1979-2020.These results indicate that the atmospheric conditions in the warm season were conducive to snow and sea ice melting.The average W 10 m from October 9, 2019, to July 18, 2020, was 6.1 m s -1 , and slightly larger than the 1979-2020 climatology (5.7 m s -1 ).More details of meteorological conditions along the MOSAiC trajectory relative to the climatology are given in Rinke et al. (2021).
The average ice motion speed (6.0 ± 4.1 km d -1 ,n¼ 270) from October 9, 2019, to July 5, 2020, was the largest of 1979-2020 (Figure 4e).This result can be related to the increased responses of sea ice motion to atmospheric forcing due to the thinning ice (Spreen et al., 2011) and the enhanced TPD in the study year associated with the strengthened positive Arctic Oscillation in January-March (Krumpen et al., 2021) and CAI in January-June.The average CAI from January to June 2020 was 8.5 hPa, which was the largest in 1979-2020 and over three times the 1979-2020 climatology (2.6 ± 2.4 hPa).The enhanced southward ice drifting during January to June 2020 resulted in the MOSAiC stations being advected faster from the central Arctic Ocean to Fram Strait than in 1979-2011 by about 50 d (Lei et al., 2016).The shorter presence of ice floes within the central Arctic Ocean would lead to a reduction of the time experiencing the relatively cold atmospheric forcing in the north (e.g., Lei et al., 2018).
The data from satellite-based microwave radiometers indicate that the snow depth within the 25-km radius from the MOSAiC CO did not differ much from the long-term mean of 2006-2019 (Krumpen et al., 2021).However, the initial ice thickness of 0.87 m around the MOSAiC CO obtained on October 15, 2019, from the CryoSat-2-SMOS was much smaller than that (1.20 ± 0.45 m) obtained at the same time during years 2011-2019.The CryoSat-2-SMOS thickness increased from 1.03 m in early November 2019 to 2.40 m by early April 2020, with a growth rate (0.0086 m d -1 ) larger than that (0.0052 m d -1 ) obtained by IMB measurements.This difference can be explained by 1) the dynamic thickening due to ridging and rafting that is included in the satellite measurement (e.g., Koo et al., 2021), which is not caught by the IMB measurements, and 2) the potential overestimation of ice growth rate by the CryoSat-2-SMOS data.Through comparing the SIMBA data used in this study and airborne electromagnetic induction sounding within the MOSAiC DN, von Albedyll et al. (2022) identified that ice dynamics increased the mean thickness from fall to early summer by 0.4 m, i.e. 30% of the total change.Over the MOSAiC DN, the observations of ICESat-2 and airborne electromagnetic induction sounding from December 2, 2019, to April 11, 2020, revealed ice growth rates of 0.0078 m d -1 and 0.0073 m d -1 ,respectively (Kooetal.,2021;vonAlbedyll et al., 2022).Both were slightly smaller than that derived from the CryoSat-2-SMOS data in the same period (0.0084 m d -1 ), indicating the potential overestimation of the ice growth rate by the CryoSat-2-SMOS data.

Snow accumulation and surface melt
At the time of the buoy deployments in October 2019, the initial snow depths were 0.10-0.18m.The snow accumulation generally was an intermittent process that was related to synoptic events (grey arrows in Figure 5).The timing of the maximum snow depth differed widely for the various sites, ranging from February 23 to June 1, 2020, because dynamic processes such as snow drifting might be superimposed on the seasonal evolution of the snow cover (e.g., Wagner et al., 2021).From November 2019 to April 2020, the snow depth measured using a magnaprobe in the northern transect of the MOSAiC CO, which was mostly situated on deformed second-year ice (Wagner et al., 2021), was slightly larger than the values obtained from our buoy measurements by about 0.05 m.The difference was likely due to the choice of buoy deployment sites, which usually avoided the deformed sea ice, and hence were more conducive to snow erosion and drifting (Wagner et al., 2021).The drifting snow tended to accumulate on the leeward side of ridges, leading to a relatively larger snow depth obtained from the CO northern transect.
The average snow depth obtained from all available buoys (15) reached a maximum of 0.25 ± 0.08 m by May 2, 2020.After mid-May 2020, the snow cover started to melt when the near-surface air temperature approached -1 C, i.e. the temperature threshold indicating surface melting (e.g., Rigor et al., 2000).The onset of snow melt was earlier than the results from the SHEBA campaign by about two weeks (May 29; Perovich et al., 2003) because the MOSAiC floes experienced a relatively warm melt season.This earlier surface melt onset is likely to have a significant impact on the energy budget during the entire During the storm events in mid-November 2019 and mid-April 2020, the observations of buoy 2019T56 at the CO IS site and of 2019T72 at the M5 site reveal that the ice deformation led to the formation of snow ice because of seawater infiltration over the surface due to a local depression (Figure 6).AttheISsite,thethickness of the consolidated snow ice increased from 0.12 to 0.22 m within 2 d, accounting for 14% of the total ice thickness by November 19, 2019.This layer of snow ice remained stable until June 27, 2020, when melt occurred.The formation of snow ice at this site reduced the heat insulation of snow and consequentially an increased ice growth rate because the snow-to-ice proce s sl edtoani n c re a s ei nbo thd e n s it yan dt h er ma lc o nductivity.From December 1, 2019, to May 31, 2020, the ice basal growth rate was 0.0062 m d -1 at this site, which was slightly larger than that (0.0056 m d -1 )o b t a i n e d from the CO FYI site with a comparable initial ice thickness.At the M5 site, the snow-to-ice process occurred more suddenly, with the thickness of snow ice increasing to 0.18 m within one day, finally accounting for 9% of the total ice thickness on April 18, 2020.Different from t h eI Ss it e ,t h es no w -t o -i c ep r o c e s sa tt h i ss i t ed e c r e a s e d the vertical temperature gradient through the ice cover, leading to a stagnant ice growth period.This effect is attributed to a large fresh snow accumulation after the  As summer approaches, freshwater originating from snow melt might permeate downward and refreeze when encountering the colder ice surface, and subsequently create superimposed ice (Nicolaus et al., 2003).The contribution of superimposed ice to the sea ice mass balance was relatively weak compared to the snow ice at our measurement sites.Only a very thin superimposed ice layer (0.04-0.06 m) was observed at 4 sites, which occurred during late May to mid-June, as shown in Figure 3b and d.The thinness of this layer was because 1) the relatively thin snow cover at the buoy sites could not provide enough meltwater to create superimposed ice, and 2) the rapid snowmeltandicetemperatureriseastheicefloesdrifted to lower latitudes restricted the formation of superimposed ice.

Sea ice mass balance
At 14 available sites, the onset of ice basal growth ranged from October 14 to December 8, 2019, with an average date of November 9.At this time, the ice thickness ranged from 0.36 to 1.72 m, and a 14-d running average air temperature was -19.8 C.These data imply that the onset of ice basal growth had a long lag (about 80 d) relative to the surface freezing onset, which occurred in late August 2019 at the buoy deployment sites, assuming a threshold of -1.0 Co nt h eb a s i so f a 14-d running average reanalyzed T 2m (e.g., Rigor et al., 2000).Among the sites, the thinner ice thickness was significantly related to an earlier onset of ice growth (R 2 ¼ 0.618, P <0.0 1;Figure 7a).The onset of first-year ice growth with the thinnest initial ice thickness of 0.36 m occurred on October 16, 2019, 37 days ahead of that for the thickest multi-year ice with an initial ice thickness of 1.72 m.However, the statistical relationship between the initial ice thickness and timing of maximum ice thickness was insignificant at the level of 0.05.Thus, the significant relationship between the initial ice thickness and ice growth period (Figure 7b) was regulated mainly bytheonsetoficegrowth.Thetotalicegrowthatthe10 sites that monitored the full growth season ranged from 0.64 to 1.38 m and averaged at 0.92 ± 0.28 m.Although t h et h i n n e ri c ew o u l dh a v eal a r g e ro v e r a l lg r o w t hr a t e (Figure 7d) because of the negative conductivity feedback (e.g., Stroeve et al., 2018), its maximum ice thickness still tended to be smaller than the thicker ice (Figure 7c).This tendency implies that the enhanced ice growth caused by the negative conductivity feedback cannot fully compensate the difference in initial thickness.This insufficient sea ice growth was likely because the summer air temperature was abnormally high and the southward advection of sea ice was abnormally fast in the study year, which shortened the ice growth period.
From the time of the annual maximum ice thickness totheonsetoficebasalmelt(23± 13 d), the ice was in a state of thermal equilibrium, with the heat from the ocean being almost equal to the upward conductive heat flux at the ice bottom.However, in this stage, the drainage of meltwater originating from the melting of snow, sea ice at the surface, and internal crystal lattice through the ice cover is still possible (e.g., Perovich et al., 2021).On June 7, 2020, which was the average date of ice basal melt onset (obtained from 14 sites), the bulk average ice temperature was -2.2 ± 0.7 C, suggesting a high permeability (e.g., Golden et al., 1998).The onset of ice basal melt cannot be significantly related to both the initial ice thickness and the maximum ice thickness (Figure 7e and f) because the timing of ice basal melt is regulated mainly by the ocean heat flux, especially in the region close to the MIZ (e.g., Provost et al., 2017;Lei et al., 2018).At 12 sites available for comparison, the onset of ice surface melt (June 20, 2020, ± 9 d) lagged the onset of ice basal melt (June 8, 2020, ± 11 d) by 12 d.This time difference was much shorter than that for the difference in the freezing onsets between the surface and bottom (approximately 80 d with the delayed ice basal freezing onset).Although snow has a buffering effect on the melting of the ice surface, the statistical relationship between snow depth in later spring and the onset of ice surface melt was insignificant at the level of 0.05 because there was no large difference in snow depth among the sites.All the buoys, except for the UIS ocean unit, failed before the sea ice melted completely, which was expected and due mainly to a continued development of melt ponds, the disintegration of the floe, or the enhanced ice lateral melt close to the ice edge.Thus, we defined the deadlines of July 10 or July 20, 2020, to characterize the relative contributions of melt at the surface and at the bottom (Figure 8).From May 15 to July 10, 2020, the total equivalent melt of snow and ice at the surface of 9 available sites was 0.38 ± 0.11 m, which was about three times that (0.12 ± 0.05 m) of ice basal melt (Figure 8a).By July 20, the total surface melt increased to 0.52 ± 0.08 m, about twice that of ice basal melt (Figure 8b).These findings imply that the surface melt dominated in the early melt season.As the sea ice drifted closer to the ice edge, the ice basal melt exhibited an increased relative contribution to the total ice melt.The statistical relationship between mass equivalent ice thickness in later spring (May 15, 2020), calculated as 1/3 of the snow depth plus ice thickness, and the melted thickness of ice by July 10 or 20, 2020, was insignificant at the 0.05 level.This lack of statistical significance suggests that the ice melt process was more complex than that for ice growth, especially in the MIZ.

Under-ice melt ponds
After the onset of ice basal melt, we identified the appearance of under-ice ponds at 10 sites out of the 12 that still had operational IMBs.This percentage suggests that under-ice ponds were more prevalent in the MOSAiC DN than the SHEBA observations in the western Arctic Ocean (about 15%; Perovich et al., 2003).This difference is partly attributed to the facts that the MOSAiC ice camp was located in the MIZ after the ice melt onset, resulting in the higher degree of sea ice melting compared to that for the SHEBA, and that all the MOSAiC buoys were deployed in the level ice, which was conducive to the formation and maintenance of under-ice melt ponds (Perovich et al., 2003).Based on the HT rise ratio, we obtained an onset of under-ice pond formation ranging from June 16 to July 31, 2020, at the 10 sites at 81.8 -82.7 N north of Fram Strait.This process lagged behind the onset of the snow melt by about 30 d because of the time needed to warm the ice column and hence to increase its permeability.The average ice thickness when the snow melted completely at these 10 sites was 1.75 ± 0.13 m, which was significantly smaller than that (2.52 m) of other two sites without an under-ice melt pond.This difference implies that the appearance of under-ice melt ponds is dependent mainly on the permeability of the ice layer.The thinner ice tended to be related to an earlier onset of under-ice melt pond formation, but without a significant relationship between them at the level of 0.05 (Figure 9a).The bulk average ice temperature was -0.9 ± 0.4 C at the onset of under-ice melt pond formation, which promoted brine channel interconnectivity and meltwater discharge (e.g., Golden et al., 1998).The appearance of under-ice melt ponds was often sudden, which was likely related to the abrupt discharge of surface meltwater (e.g., Provost et al., 2019).At the onset of under-ice melt pond formation, the average thickness of the meltwater layer at 10 sites was 0.17 ± 0.10 m.The meltwater thickness increased to the maximum (0.28 to 1.26 m) by July 1 to August 3, 2020, at the 8 sites that continued monitoring through this time.The time was significantly related to the onset date of under-ice melt pond formation (R 2 ¼ 0.680, P < 0.05; Figure 9b).However, the maximum meltwater thickness cannot be related to the onset of formation of under-ice melt pond at a significance level of 0.05 (Figure 9c).The accelerated ice motion and enhanced oceanic dynamics as the ice approached the sea ice edge can promote removing the meltwater layer under the ice through the dissolution or deeper mixing (e.g., Provost et al., 2019).This mechanism weakens the dependence of under-ice pond evolution on the onset of pond formation.Sea ice topography is very important for the development of under-ice melt ponds.For example, the ice bottom sheltered by the ridge keel can act as a trap for meltwater accumulation (Perovich et al., 2003), and the ice surface topography can also play an important role because the meltwater under the ice is likely originating from horizontal permeation (Perovich et al., 2021).However, the influence of the ice topography cannot be assessed quantitatively as no such measurements were conducted at the buoy sites.Here, we used a pair of contrasting examples to describe the effect of surface topography qualitatively.These two sites, with similar ice thicknesses of 1.80 m (L3 SIMBA site) and 1.82 m (CO IS site) when the snow melted completely, had an obviously inconsistent onset (June 16 versus July 28, 2020) and maximum thickness (0.82 m versus 0.28 m) of the under-ice meltwater layer.These inconsistencies are likely due to the difference in local surface topography, possibly also associated with the under-ice roughness, between the sites.The site at L3 had a relatively flat ice surface at the deployment in early October 2019, while the IS site was close to a ridge.Level ice is more prone to surface ponding than deformed ice (Eicken et al., 2004).The total mass equivalent surface melt of the site at L3 was 0.38 m by July 10, 2020, which was more than twice the value at the CO IS site (0.18 m).The development of surface melt pond might increase the meltwater source (e.g., Perovich et al., 2021), thus promoting the evolution of under-ice melt ponds.Moreover, once the keels started to melt at the CO IS site, the meltwater was likely to spill over and redistribute, which was not conducive to the accumulation of meltwater.
After the appearance of under-ice melt ponds, the heat from the ocean becomes partly isolated from the ice by the meltwater layer.However, due to the relatively high temperature of the under-ice meltwater (e.g., Figures 3a,  3c, and 6a), the ice bottom continued to melt gradually.At 4 sites with continuous measurements of under-ice ponds, the ice basal melt rate ranged from 0.0040 to 0.0027 m d -1 in the period of late June to late July 2020.

Heat fluxes through and below the ice cover
The phase change between brine and ice crystals within the sea ice column is expected to regulate the heat flux through the ice column, in turn weakening the seasonal responses at the change in temperature in lower ice layers or of ice growth to the atmospheric forcing.Thereby, with ice thickening and an increase in air temperature, the heat flux in the top layer had a relatively larger seasonal attenuation compared to the bottom (Figure 10a and b).Our observations show that the ratio between them decreased significantly from 21.5 ± 26.8 in October 2019 (n ¼ 20) to 0.5 ± 0.6 by May 2020 (n ¼ 20).In June, the direction of conductive heat flux in the top and bottom layers was opposite, with heat conducted into the ice interior, which promoted the internal melting of the ice.At the measurement sites, the thinner snow-covered sea ice had a larger bulk average conductive heat flux (Figure 10c) and consequentially a larger ice growth rate, which constitutes the so-called negative feedback mechanism of ice thicknessconductive heat flux.The difference between the sites decreased seasonally with the relatively faster thickening of thinner ice (Figure 10c and d).Thus, the thinning of sea ice over the pan Arctic Ocean would promote a heat loss from the ice-ocean system to the atmosphere, especially during the autumn and winter, when the ice growth rate is relatively large.
For the heat balance at the ice bottom, we first analyzed the dataset obtained by the UIS to characterize the physical coupling mechanism between the sea ice and the upper ocean layer (Figure 11).At this site, the ice growth lasted from November 13, 2019, to May 25, 2020, with the ice growth rate being regulated mainly by the conductive heat flux similar to other sites (Figure 11a-c).The ice basal growth led to the formation of dense water in the oceanic mixed layer, which contributed to eroding vertical density gradients during November and December 2019, as recorded by the UIS ocean unit (Figure 11d  and e).The strong low-pressure synoptic systems were likely to promote the upward entrainment of warm water from the lower layer (Figure 11d and e).The rapid seasonal heat loss from the mixed layer and the convection associated with brine drainage induced the supercooling that occurred at the topmost mixed layer, with a temperature lower than the local freezing point (e.g., Hoppmann et al., 2020;Katlein et al., 2020).Although the onset of ice basal growth at the UIS site occurred in mid-November 2019, the supercooling at the top of the mixed layer appeared already in late October.This temporal gap was likely because the UIS site had relatively thick ice, while the thinner ice on the other part of the floe may have started to grow earlier.For example, ice basal growth at the SIMBA site on the same floe started earlier on October 14, 2019, due to the thinner initial ice thickness (0.50 m).
During October and November 2019, the layer with supercooled temperature oscillated between 10 and 20 m, with temporary shoaling associated with the upward mixing of lower warm water.The supercooled water occupied a significant part of the oceanic mixed layer at this stage.According to observations by the CT sensor packages of the UIS ocean unit, the supercooled layer steadily thickened from about 14 to 21 m during early and mid-December, which was likely related to the accelerated heat loss.It remained in a relatively stable state until mid-February 2020, when the layer started to steadily shoal and disappeared by mid-April 2020.Accordingly, the negative departure from the freezing point interpolated to a depth of 10 m averaged at 0.007 ± 0.002 C (n ¼1464) from October 10 to December 10, 2019, which increased to 0.010 ± 0.006 C(n¼1464) during December 11, 2019, to February 10, 2020, and decreased again to 0.007 ± 0.002 C(n¼1080) during February 11 to March 28, 2020.The supercooled water may lead to the formation of platelet ice within the water column or attaching to the ice bottom.From late December to mid-March, widely distributed platelet ice was observed under the MOSAiC CO floe by a remotely operated underwater vehicle (ROV) (Katlein et al., 2020), indicating that the supercooled water appearing under the ice at the L3 site was not an isolated exception within the MOSAiC DN area.However, from our thermistor chain measurements, we could not identify the appearance of platelet ice.The distribution of platelet ice was likely mostly patchy, which makes identifying by observations of individual thermistor chains difficult.
From mid-February 2020 onward, the oceanic mixed layer became much saltier as the floe approached and in the average and standard deviation of conductive heat flux through the entire ice cover, as well as the average thermal equivalent ice thickness obtained from all available sites.Note that the thermal equivalent ice thickness, different from the mass equivalent ice thickness, was estimated by assuming that the snow has a thermal insulation capacity (inversely proportional to the thermal conductivity coefficient) 6.7 times that of sea ice (Sturm et al., 2002).
crossed the Gakkel Ridge (Figure 11f), also making the transition of hydrologic regimes from the Amundsen Basin to the Nansen Basin.In the Nansen Basin, the weaker halocline allows vertical ventilation of the Atlantic Water up to the surface (Korhonen et al., 2013;Athanase et al., 2019).Thus, the disappearance of the supercooled state of the oceanic mixed layer after April 18, 2020, was due mainly to the spatial changes in the hydrologic regime.From June 4 to July 17, 2020, the UIS drifted across the northwest Yermak Plateau, over the shallowest bathymetry of 573 m along the trajectory.During this period, the oceanic mixed layer became warmer, which was likely due to topographically enhanced upward mixing (e.g., Fer et al., 2015) and potential tidal mixing because this region is a tidal hotspot with elliptic tidal currents (e.g., Padman et al., 1992;Rabe et al., 2022).From July 29, 2020, onwards, the floe drifted onto the northeast continental shelf of Greenland, with the upper ocean warming further.
Although the UIS and one SIMBA (2019T70) were deployed on the same floe with a distance about 100 m between them, the estimated oceanic heat fluxes revealed some deviations at times (Figure 11f).These deviations were likely related to local changes in the ice basal topography and the typical uncertainty (1-2 W m -2 )o ft h e "residual estimation" method (Lei et al., 2018).Previous studies have shown that the oceanic heat flux may change greatly at the scale of tens of meters due to the variability of the ice bottom topography through changing the friction between ocean and ice (Wettlaufer, 1991).
Second, we also estimated the ensemble average of oceanic heat flux using the data from all available buoys to reduce the uncertainty of the "residual estimation" method.The oceanic heat fluxes obtained from various buoys have a generally consistent seasonal variation (Figure 12a).The ensemble average oceanic heat flux (Figure 12b) decreased from 4.4 ± 0.7 W m -2 in mid-October 2019 to the annual minimum of 1.6 ± 1.0 W m -2 by early December.It remained at a low level until the end of April 2020.The overall average oceanic heat flux was 2.8 ± 1.1 W m -2 from December to April, which was very close to the typical winter value (2.1 ± 2.3 W m -2 ) in the Arctic Eurasian Basin (Peterson et al., 2017).This nonnegligible value suggests that the supercooling at the topmost mixed layer did not completely isolate the oceanic heat flux from the lower layer.From May 2020 onwards, the heat flux increased rapidly to 10.0 ± 2.6 W m -2 by mid-June 2020.In Fram Strait and on the northeast shelf of Greenland during July to September 2020, the relatively large heat content in part of the oceanic mixed layer (Figure 11f) suggests a higher oceanic heat flux.

Seasonal and spatial aspects
Based on data from 8 buoys with long-lasting operation from October or early November 2019 to July 2020 (Figure 13), the average onset of ice basal growth was November 9, 2019; the bulk ice temperature reached a minimum (-11.0C) by March 10, 2020; from May 16, 2020 onwards, the bulk ice temperature was warmer than -5.0 C, possibly increasing the ice permeability (Golden et al., 1998); the ice thickness reached a maximum (1.90 m) by May 20, 2020; the ice basal and surface melt commenced on June 8 and 20, 2020, respectively; and the under-ice ponds appeared on June 26, 2020.The occurrence of these events obtained from the above 8 selected buoys was very close to that from all available samples (generally 5-7 more), with the largest difference of 7 d for the onset of under-ice melt pond formation.Thus, the statistical results obtained for specific variables using the data from all available buoys given in the above sections are generally robust.The relatively large difference for the onset of under-ice pond formation was because it is a complex process depending on many different factors such as ice surface and bottom topographies, permeability and texture of the ice column, and the amount of available meltwater.
The results given in this study not only characterize the seasonality of the sea ice mass balance, but also include the superposition of a spatial pattern, especially during the accelerated regime transition from the central Arctic Ocean to the MIZ north of Fram Strait.When the ice floes drifted to the MIZ by mid-July 2020, they still had a thickness greater than 1 m (Figures 3 and 13).Subsequently, oceanic heat and dynamic breakup played a critical role in the final melting of the sea ice.After entering Fram Strait, the southward drift of the sea ice was accelerated.This acceleration led to the ice absorbing more solar radiation, and experiencing higher air temperature and more frequent rainfall events at the south latitudes, which further promoted the melt and ultimately resulted in disintegration of the floes.
In 1979-2011 (Lei et al., 2016), the average transport time of ice floes from the central Arctic Ocean to Fram Strait would have taken nearly two months longer compared to the MOSAiC floes.In these earlier years, the growth period of the sea ice would have been prolonged and the time to start melting would have been delayed.As a result, the possibility of sea ice surviving through summer would increase.As a comparative example, one SIMBA was deployed by the Chinese National Arctic Research Expedition (CHINARE) north of the Laptev Sea (85.1 N, 147.4 E) on September 3, 2012 (Lei et al., 2018), which was close to the initial MOSAiC deployment location (85.1 N, 132.9 E) on October 4, 2019.Ten months after the initial deployment, the MOSAiC DN had drifted to Fram Strait by early August 2020; however, the CHI-NARE buoy was still in the region close to the North Pole (88.9 N, 59.7 W).In the 2012-2013 ice season, the ice growth at the CHINARE site ceased by late June, which was about one month later than that obtained from the MOSAiC buoys.However, the total ice growth (0.31 m) was only about one third of the MOSAiC buoys (0.92 ± 0.28 m) because of the thinner initial thickness at the MOSAiC sites than the CHINARE site (0.96 ± 0.48 m versus 2.40 m).The earlier onset of ice basal melt observed at the MOSAiC sites also resulted in a relatively strong summer ice melt.The total equivalent-ice melt from surface and bottom obtained from the MOSAiC buoys by mid-Summer (July 20, 2020) was 0.77 ± 0.18 m, which was even slightly larger than that obtained from the CHINARE buoy for the entire melt period in 2013 (0.60 m).The enhanced melting was also a potential contribution to the widespread emergence of under-ice melt ponds at the MOSAiC buoy sites.
This difference between two ice seasons can ultimately be related to the difference in the atmospheric circulation pattern.During the corresponding 10-month period from October to August, the CAI index during the 2012-2013 ice season was 2.0 hPa lower than the 1979-2020 climatology, while that in 2019-2020 was 5.0 hPa higher than the climatology.The relatively high CAI in 2019-2020 was likely a consequence of large-scale low-pressure anomalies prevailing around the Barents-Kara-Laptev Sea region between January and March.Therefore, the atmospheric circulation of 2019-2020 was more conducive to the enhanced outflow of sea ice via the TPD.In contrast, the CHINARE buoy had drifted to Fram Strait on December 25, 2013, having spent about six months longer in the central Arctic Ocean compared to the MOSAiC buoys.

Conclusions and outlook
The seasonality and timing of thermodynamic mass balance and heat fluxes of sea ice in the Arctic Transpolar Drift were investigated using observations from an array of 23 IMBs spanning a distance of approximately 50 km as part of the MOSAiC drift experiment in 2019-2020.This IMB dataset represents one of the most comprehensive observational records of sea ice growth and decay in the TPD.In the early setup stage of October and early November 2019, initial ice thicknesses ranged between 0.35 m and 1.80 m at the buoy sites.The seasonality of thermodynamic sea ice mass balance can be divided into four stages (Figure 13) :o n s e to fi c eb a s a lf r e e z i n gi nm i d -October-November (I), a rapidly-growing stage in December-March (II), a slowly-growing stage in April-May (III), and a melting stage from June onwards (IV).In stage I, the average date of basal freezing onset was November 9, 2019 (± 18 d), which depended mainly on the initial ice thickness.During stage II, the thinner ice exhibited a relatively large ice growth rate caused by the negative conductivity feedback.Most snow accumulation also occurred during this stage, and was related to intermittent synoptic events.Supercooling in the topmost oceanic mixed layer, observed during stages I and II, was caused by heat loss from the mixed layer and the convection due to brine drainage.During stage III, the increases in ice thickness, air temperature, and oceanic heat flux resulted in a decrease of the ice growth rate.The average ice thickness increase was 0.92 ± 0.28 m for the entire ice growth season, which on average lasted 193 ± 26 d and ended by May 17, 2020 (±11 d).During the melting stage (IV), the oceanic heat flux increased rapidly to 10.1 ± 2.6 Wm -2 by the end of June.This increase was due mainly to an increased interaction between the ice pack and the warmer upper ocean, caused mainly by a bathymetryinduced upward mixing, especially over the Yermak Plateau.The surface melt dominated in the early melt season.As the sea ice drifted closer to the ice edge, the ice basal melt increased due to the increased oceanic heat flux.From June 26, 2020, onwards, the sea ice concentration dropped intermittently below 95%.I c ef l o e sb e g a nt o melt from the base, surface, and edges.Snow and ice meltwater frequently appeared under the thinning ice to form under-ice melt ponds during this stage.
The seasonal signal of the atmospheric forcing was stronger than its spatial differences with respect to driving the sea ice thermodynamic mass balance.For example, near-surface air temperature, shortwave radiation and longwave radiation all showed strong seasonal signals during the MOSAiC observational period (Rinke et al., 2021).However, the spatial change in the oceanic forcing was comparable to its seasonal signal.For example, both the transition of hydrologic regimes from the Amundsen Basin to the Nansen Basin in March-April 2020 and the advection of the MOSAiC floes over the Yermak Plateau and Fram Strait during June-July 2020 promoted rapid changes in the oceanic mixed layer (Rabe et al., 2022), and consequently regulated the heat transfer from the ocean to the ice.Unusually high air pressure gradients perpendicular to the TPD, especially during January-June 2019, resulted in faster advection of sea ice from the central Arctic to Fram Strait in the study period compared to the climatology.Subsequently, the sea ice experienced a relatively warm summer in the region close to Fram Strait compared to the climatology.Thus, the observed atmospheric circulation pattern led the sea ice to enter the melting period earlier under the combined forcing of atmosphere and ocean in the region close to or within Fram Strait.When using the MOSAiC sea ice time series data for numerical modelling studies, both the representativeness and peculiarity of the atmospheric and oceanic forcing experienced by the MOSAiC floes relative to the climatology will be important to keep in mind.
The total sea ice mass balance is generally governed by both thermodynamic (e.g., freezing and decay) and dynamic (e.g., advection and deformation) processes.
While this study focused on the thermodynamics of sea ice derived from IMB data, a comparison to the airborne electromagnetic induction sounding measurements across the MOSAiC DN showed that the thermodynamic thickening of the sea ice accounts for 70% of the total thickening from fall 2019 to early summer 2020 (von Albedyll et al., 2022).However, the thermodynamic and dynamic processes of sea ice growth are not isolated from each other.For example, storm events brought rises of air temperature, snowfall, and enhanced drifting and blowing snow and ice deformation.Ice deformation led to local surface depression, resulting in the formations of slush and snow ice.The acceleration of sea ice motion, associated with synoptic events, led to the upward mixing of lower warm water and affected the growth of sea ice.
To obtain a full picture of the sea ice mass balance still requires the use a combination of multi-source observations, e.g.repeated ROV sounding, thickness stake, ice coring, snow pits, transect measurements, and airborne observations of electromagnetic induction.The comprehensive measurements during MOSAiC constitute excellent data sets to identify the variations in sea ice composition and the distribution of ice thickness on multiple scales.Combined with the observations of the whole MOSAiC buoy array or Synthetic Aperture Radar to describe the ice deformation, we can identify in future analyses the contributions of various components to the total sea ice mass balance dynamically and thermodynamically.The seasonality and timing of sea ice thermodynamic mass balance derived from this study can support many future MOSAiC multi-disciplinary studies, in particular, to explain the seasonal changes in energy and matter exchanges between atmosphere, sea ice and ocean, the salt and heat budgets of the upper ocean, and the biological communities within and under the ice.

Figure 1 .
Figure 1.Schematic diagram of the unmanned ice station (UIS).The UIS was deployed at site L3 of the MOSAiC Distributed Network (see Figure 2), which included sensors for conductivity, temperature, and depth (CTD), as well as others not used in this study: optical sensors and those for dissolved oxygen (DO) and chlorophyll a (Chl a).

Figure 2 .
Figure 2. Trajectories and operational period of all buoys in this study.(a) Drift trajectories of the IMBs (thin blue lines) and the UIS ocean unit (colored line denoting the time), with two insets denoting the buoy positions on December 1, 2019, and July 1, 2020; (b) the operation periods of the buoys, with solid and dashed lines denoting the buoys with measurements that fully or partially covered the ice growth period, and red, brown, blue and green lines denoting the buoys deployed at the CO, L, M, and P sites, respectively.

Figure 3 .
Figure 3. Identifications of air-snow-ice-water interfaces.The colored plots show the changes in profiles of (a) temperature and (b) HT rise ratio obtained by SIMBA 2019T63 deployed in relatively thick ice at the L2 site and (c) temperature and (d) HT rise ratio by 2019T70 deployed in relatively thin ice at the L3 site.Also shown are the snow or ice surface (dark grey solid line), the ice bottom (purple solid line), and the lower boundary of the under-ice melt pond (white solid line), as well as the initial snow-ice interface at the 0-m mark (grey dashed line).The black rectangles in (b) and (d) highlight the layer of superimposed ice.The white areas in (a) and (b) are data gaps.

Figure 4 .
Figure 4. Atmospheric and sea ice conditions relative to the climatology.(a) Daily air pressure obtained from the UIS measurements and atmospheric reanalysis data in 2019-2020 and the average and standard deviation of 1979-2020; (b) the same as panel a, but for air temperature; (c) the cumulative freezing degree days derived from the reanalyzed T 2m of 2019-2020 and 1979-2020; (d) the same as panel a, but for wind speed, and the measured data were obtained from the ship-based weather station, which was unavailable from May 29 until June 13, 2020, when the ship temporarily left the ice camp; (e) directional distribution of the ship-based wind velocity in October-July, January-July, and March-June; (f) the same as panel a, but for ice motion speed; and (g) changes in ice thickness obtained from the CryoSat-2-SMOS dataset in 2019-2020 and the average and standard deviation of 2011-2020, as well as that derived from 11 IMBs with continuous measurements from November 2019 to April 2020.

Figure 5 .
Figure 5. Changes in the air-snow or air-ice interfaces, and near-surface air temperature.Variations in air-snow or air-ice interfaces at the sites with snow-ice formation (brown and purple lines) and other available sites (grey lines), as well as near-surface air temperature measured by the UIS.The red line and blue shade denote the average and standard deviation of snow depth obtained from 8 buoys operating from October 2019 to at least the onset of ice surface melt.Also shown are the temporal range of maximal snow depth (blue double arrow) and the average (blue arrow), the temporal range of ice surface melt onset (green double arrow) and the average (green arrow), as well as four distinct storm events (grey arrows).The zero line denotes the initial snow-ice interface.
f o r m a t i o no fs n o wi c e ,w i t ha b o u t0 .3 0mo fn e ws n o w overlaying the snow-ice layer by April 20, 2020.

Figure 6 .
Figure 6.IMB measurements at two sites with snow ice formation.Profiles of (a) temperature and (b) HT rise ratio obtained by the buoy 2019T56 at the CO IS site and (c) temperature and (d) HT rise ratio by 2019T72 at the M5 site.The blue solid line denotes the interface between snow ice and normal snow.Also shown are the snow or ice surface (dark grey solid line), the ice bottom (purple solid line), and the initial snow-ice interface at the 0-m mark (grey dashed line).

Figure 7 .
Figure 7. Linear regressions between the timing of ice mass balance and ice thickness.Statistical relations between initial ice thickness and (a) onset date of ice basal growth, (b) ice growth period, (c) maximum ice thickness, (d) ice growth rate, (e) onset date of ice basal melt, and (f) between maximum ice thickness and onset date of ice basal melt.The dashed blue lines denote a linear regression.No significant relation was identified for the paired variables in (e) and (f).The red circles denote the measurements that cover the entire ice growth season; the black circles, only cover the two concerned variables.

Figure 8 .
Figure 8. Basal and surface equivalent ice melt.The melted thickness of ice by (a) July 10 and (b) July 20, 2020.Also shown is the total mass equivalent ice thickness on May 15, 2020 (red line).

Figure 9 .
Figure 9. Statistics of timing of under-ice melt pond formation.Relations between (a) the onset of under-ice pond formation and the ice thickness at the time when the snow melted completely, (b) the time with maximum under-ice meltwater thickness and pond onset, and (c) the maximum under-ice meltwater thickness and pond onset.The dashed line in (b) denotes the linear regression but no significant relation was identified for the paired variables in (a) and (c).

Figure 10 .
Figure 10.Changes in conductive heat flux and thermal equivalent ice thickness.Monthly conductive heat flux through (a) the top ice layer (0.05-0.15 m from the initial ice surface), (b) the basal ice layer (0.05-0.15 m from the ice bottom), and (c) the entire ice cover as a function of the thermal equivalent ice thickness.Also shown are (d) changesin the average and standard deviation of conductive heat flux through the entire ice cover, as well as the average thermal equivalent ice thickness obtained from all available sites.Note that the thermal equivalent ice thickness, different from the mass equivalent ice thickness, was estimated by assuming that the snow has a thermal insulation capacity (inversely proportional to the thermal conductivity coefficient) 6.7 times that of sea ice(Sturm et al., 2002).

Figure 11 .
Figure 11.Interaction mechanism between sea ice and the oceanic mixed layer at the L3 site.Changes in (a) near-surface (1.5 m) air temperature and pressure, (b) ice temperature profiles, (c) conductive heat flux (F c ) through the ice (positive upward), (d) seawater salinity at depths of 6 to 40 m, (e) the departure of temperature from local freezing temperature, and (f) oceanic heat content at layers of 6-35, 6-15, 15-35 m obtained by the UIS.The lines shown in (b) and (c) denote the interfaces between air, snow, ice and water as Figure 3.The dark grey areas in (d) and (e) denote data gaps.The blue line in (e) indicates the lower boundary of the supercooling layer.The red line in (f) denotes the bathymetry obtained from the Version 4.0 International Bathymetric Chart of the Arctic Ocean.

Figure 12 .
Figure 12.Seasonal changes in oceanic heat flux.(a) Changes in the 15-day average oceanic heat flux obtained from all sites, and (b) their average and standard deviation (n corresponding to the buoy number shown in the above panel).The dark blue color in (a) indicates data gaps.

Figure 13 .
Figure13.Time series of sea ice mass balance.Time series of averages (solid line) and standard deviations (shade) of snow depth (dark blue), ice bottom (light blue), and bulk average ice temperature (red) obtained from 8 long-lasting buoys.The ice bottom is only plotted after all these buoys were deployed because of the large difference in ice thickness between the sites.The vertical purple lines denote the timing of sea ice mass balance obtained from the 8 selected buoys; the vertical yellow lines, from all available buoys.The horizontal black arrows denote the four stages for the entire ice season.

Table 1 .
Identification methods of the interfaces between air, snow, sea ice and water