Snow thermal conductivity and conductive flux in the Central Arctic: Estimates from observations and implications for models

Snow and sea ice are the important components of the Arctic system, regulating the exchange of heat and moisture between the ocean and atmosphere. In order for sea ice to grow, energy must be conducted away from the bottom of the ice for the underlying sea water to cool and freeze. However, snow is a thermal insulator, and its presence on sea ice reduces sea ice growth by slowing energy transport at the surface with an impact that increases as sea ice thickness decreases (e.g., Maykut, 1978; Sturm et al., 2002).

Quantifying the exchanges of energy at the surface is crucial not only for understanding how the Arctic evolves on seasonal to interannual timescales but also for correctly representing the Arctic system in models. In winter over sea ice, net longwave radiation (LWNET), and turbulent heat fluxes are balanced at the surface by heat conduction through the sea ice and snow. This conductive flux represents the loss of energy at the sea ice bottom in order for sea ice to grow. Conduction is proportional to the temperature gradient through ice and snow, which is a function of sea ice and snow thicknesses and the temperature difference between the surface and underlying ocean. Because conduction depends on the surface temperature, it couples atmospheric variability to sea ice growth. As the Arctic continues to warm and sea ice continues to thin, the magnitude and variability of conductive fluxes are likely to change (e.g., Bigdeli et al., 2020; Landrum and Holland, 2022).

Conductive flux also depends on the thermal conductivity of the medium, or how readily energy is transported via conduction. The measurements of snow thermal conductivity range from less than 0.1 W m⁻¹ K⁻¹ for fresh snow to 0.5 W m⁻¹ K⁻¹ or higher for wind slabs from field observations and laboratory measurements (e.g., Sturm et al., 2002; Calonne et al., 2011). This range reflects differences in snow properties, including density, grain structure, and bonding (e.g., Brandt and Warren, 1997). To further complicate matters, these snow properties are known to vary vertically, spatially, and temporally as the snowpack undergoes metamorphism (e.g., Colbeck, 1982). Due to density considerations, all values for snow are less than the thermal conductivity of sea ice, which is around 2 W m⁻¹ K⁻¹, although this too can vary depending on sea ice properties, such as density, porosity, and salinity (e.g., Pringle et al., 2007; Petrich and Eicken, 2010). The spatial variability of snow thermal properties and depth can lead to large variability in conductive flux. At the Surface Heat Budget of the Arctic Ocean project, Sturm et al. (2001) found 4-fold variations in conductive flux across distances of tens of meters due to differences in snow and ice thicknesses and properties.

Despite the large and influential variability of snow properties, the representation of snow on sea ice is relatively simple in most forecast and global climate models as well as reanalyses. Typically, these models treat snow thermal conductivity as a constant for snow on sea ice, commonly 0.3 W m⁻¹ K⁻¹ (e.g., Hunke et al., 2017; Ridley et al., 2018). Not all forecast models even include snow on sea ice (e.g., Bauer et al., 2016; Bengtsson et al., 2017). This exclusion has been related to surface temperature biases, which is important to consider given the common practice of assessing models and forcing models with reanalyses (Batrak and Müller, 2019). The complex modeling of snow is possible (e.g., Liston et al., 2020), but it is typically only done on land because the transient nature of sea ice makes it computationally more expensive to trace snow properties on sea ice that is advected and undergoes strong seasonal evolution (Hunke et al., 2010).

Correctly modeling the variability of conductive fluxes is important for simulating meaningful representations of sea ice thickness and, therefore, sea ice variability. Urrego-Blanco et al. (2016) found that in the Los Alamos Sea Ice Model (CICE), sea ice was more sensitive to the representation of snow than any other parameter, specifically snow thermal conductivity and grain size. Lecomte et al. (2013) tested several parameterizations of snow thermal conductivity in a coupled sea ice-ocean model, and both the mean sea ice state and variability were impacted by the representation of snow thermal conductivity. Wu et al. (1999) decreased snow thermal conductivity from 0.3 to 0.16 W m⁻¹ K⁻¹ in a “zero-layer” coupled atmosphere-sea ice model. This change led to thinner sea ice distributions for both poles and unexpectedly warmer surface temperatures in the Arctic. They explained the warmer surface temperatures as being influenced by changes in sea ice thickness more so than snow thermal conductivity. Blazey et al. (2013) found that increasing snow thermal conductivity in a fully coupled global climate model increased sea ice extent in the Arctic and decreased pan-arctic surface temperatures. Low-level cloud cover also decreased, consistent with the fall sea ice cloud feedback (Kay and Gettelman, 2009), and suggesting that the representation of snow on sea ice is also important for modeling the atmosphere. Yet despite the known importance of snow on sea ice, and in particular its thermal conductivity, better representation is still needed in models (e.g., Batrak and Müller, 2019).

Observations from the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) expedition provide a unique opportunity to study snow as a part of the coupled ice-snow-atmosphere Arctic system. During MOSAiC, the German icebreaker Polarstern was frozen into the sea ice and observations were taken across a regional domain beginning in October 2019 and continuing to the following September (Nicolaus et al., 2022; Shupe et al., 2022). MOSAiC provides some of the most comprehensive observations from the Central Arctic during winter, including coordinated observations of sea ice, snow, and the surface energy budget (SEB).

The goal of this work is to study the variability of snow thermal conductivity on sea ice and conductive flux during the MOSAiC winter. For context, we first present a sensitivity study using a regional coupled forecast model to demonstrate the impact of varying snow thermal conductivity on the SEB. Next, we introduce the observations and present a method for retrieving profiles of snow thermal conductivity and density that are used to calculate conductive flux. Results are then discussed for spatial and temporal variability of snow thermal conductivity and conductive flux through the snowpack. Finally, relationships between temperature, sea ice, and conductive flux from observations are used to further assess the coupled forecast model.

Most modeling studies have focused on the sensitivity of sea ice to changing snow thermal conductivity, but we also expect the adjacent atmosphere to respond since conduction couples the sea ice to the overlying air via the surface skin temperature. Here, we use the Coupled Arctic Forecast System (CAFS; Solomon et al., 2023), a fully coupled regional forecast model, to demonstrate the impact of snow thermal conductivity on the winter SEB over sea ice. This example simply serves as a context for the importance of snow thermal conductivity for modeling the near-surface boundary layer, while more comprehensive simulations and sensitivity studies are warranted but beyond the scope of this study. CAFS uses the following components: CICE version 5 with 3 snow layers and 7 ice layers, WRF ARW version 3.6.1, POP ocean model version 2, and the NCAR CLM land model version 4.5. Global 0.5° GFS forecast fields are used for lateral forcing. Output data are 6-h averages with 10-km resolution over the Arctic basin domain, shown in Figure S1. In addition to a control simulation with a snow thermal conductivity of 0.3 W m⁻¹ K⁻¹, 2 additional simulations were conducted with the values of 0.15 W m⁻¹ K⁻¹ and 0.6 W m⁻¹ K⁻¹ changed over the 16 grid cells poleward of 89.78°N. For each snow thermal conductivity value, a single 10-day simulation is analyzed for January 30 through February 8, 2020. After February 8, the large-scale impact of the varied snow thermal conductivity, for example, shifts in cyclone tracks, makes isolating the impacts on the local SEB difficult. However, the chosen week has both clear and cloudy conditions that uniquely respond to changes in snow thermal conductivity. Average snow depth increases from 17.5 to 20.5 cm during the week, and sea ice thickness is 2.5 m. Only grid cells from the domain where snow thermal conductivity is changed are analyzed.

Snow thermal conductivity determines how efficiently energy is transported through the snow to the snow–atmosphere interface. The domain-average conductive flux increases by 7.3 W m⁻², or about 34%, when snow thermal conductivity is doubled to 0.6 W m⁻¹ K⁻¹, and it decreases by 6.6 W m⁻², or about 30%, when snow thermal conductivity is halved to 0.15 W m⁻¹ K⁻¹ (Figure 1a and c). Although conductive flux is linearly proportional to thermal conductivity, the conductive flux is not twice as large after doubling snow thermal conductivity. A larger conductive flux means more energy flows to the surface and the surface temperature is higher, on average by 0.85°C when snow thermal conductivity is doubled from 0.3 to 0.6 W m⁻¹ K⁻¹ (Figure 1b and d). This increase in surface temperature acts as a negative feedback since warmer surface temperatures decrease conduction, mitigating some of the impact of larger snow thermal conductivity. The opposite is also true for halving snow thermal conductivity. The surface temperature decreases by an average of 0.8°C, and the conductive flux decreases by less than half.

The impacts of changing snow thermal conductivity depend on atmospheric conditions. When the surface is warm, the temperature difference between the surface and ocean is smaller, so less conduction occurs. When the surface is coldest, the temperature gradient is relatively large and conduction is greater. This basic relationship is how the atmosphere can impact sea ice growth. When the atmosphere is cloudy, there is more downwelling longwave radiation (LWD) to insulate the surface, such as February 1–2, 2020. Accordingly, the snow thermal conductivity has comparatively less impact on the conductive flux during this time than when the surface is colder. The largest impacts of changing snow thermal conductivity occur when the atmosphere is clear and the surface is coldest, for example, January 30, 2020. Because of the different responses of conduction when the surface is warm or cold, that is, the atmosphere is cloudy or clear, changing snow thermal conductivity also alters the relationship of conduction to LWD (Figure 2a). Note that there are no changes to LWD between simulations. The response, or sensitivity, of conduction to LWD, given by the slope of the best fit line in Figure 2a, increases from 0.21 to 0.38 (an 81% increase) when snow thermal conductivity increases from 0.15 W m⁻¹ K⁻¹ to 0.6 W m⁻¹ K⁻¹. This difference means that the coupling between the surface and atmosphere changes as a result of snow properties.

Because changing snow thermal conductivity alters the surface temperature, upwelling longwave radiation (LWU) and turbulent fluxes also change between simulations. Larger snow thermal conductivity and surface temperature lead to larger LWU, by up to 9 W m⁻², and larger upward sensible heat (SH) flux, by up to 15 W m⁻². Changes to latent heat flux are relatively modest given their already small values (not shown). Again, the sensitivities of LWU and SH to LWD differ between the simulations. As conduction becomes more sensitive to the atmosphere with larger snow thermal conductivity, LWU and SH become less sensitive to LWD (Figure 2b and c). The largest changes in LWU and SH occur when the atmosphere is clear and the surface is colder. With larger snow thermal conductivity, the surface is warmer than in simulations with smaller snow thermal conductivity, so there is more LWU (i.e., Boltzmann’s law). The surface temperatures change less in response to increased LWD when the surface is already warm, so the slope between LWU and LWD decreases from −0.48 to −0.41 when snow thermal conductivity is increased from 0.15 W m⁻¹ K⁻¹ to 0.6 W m⁻¹ K⁻¹ (blue circles to red diamonds in Figure 2b). At the same time, because the surface is warmer with higher snow thermal conductivity, the near surface temperature gradient is reduced during clear and cold periods. Therefore, SH is smaller and its sensitivity to LWD is reduced from −0.29 to −0.19 when changing snow thermal conductivity from 0.15 W m⁻¹ K⁻¹ to 0.6 W m⁻¹ K⁻¹ (Figure 2c). Overall, with higher snow thermal conductivity, conduction becomes more responsive to atmospheric radiative forcing, while SH and LWU become less responsive.

The relationship between LWNET and SH also changes between these simulations, notably when the atmosphere is radiatively clear (LWNET ≤ −20 W m⁻²; Figure 2d). With the lowest snow thermal conductivity (blue circles), LWNET is less negative because the colder surface results in less LWU. SH is more positive because of the larger temperature difference between the atmosphere and surface. Thus, the thermal properties of snow in a coupled model are important not only for representing sea ice but also for the lower atmospheric structure and the interactions among different heat transfer processes.

This experiment presents a practical demonstration for the need to constrain the mean and variability of the thermal properties of snow in coupled sea ice models. This will be the objective of the remainder of this study using observations made across the MOSAiC distributed network (DN).

For the MOSAiC expedition, Polarstern was frozen into an ice floe on October 4, 2019, in the Central Arctic Ocean. We use observations during winter from October 2019 until March 1, 2020, during which time the Polarstern passively drifted with the surrounding sea ice. The Central Observatory (CO) consisted of a suite of instruments on the ice floe within 2 km of the Polarstern. We use broadband LWU and LWD measured at 3 m and 1 m above the surface, respectively, from a station located within a few tens of meters from the 10-m tower. Additional instruments were deployed as part of the DN (Krumpen and Sokolov, 2020) up to tens of kilometers from the CO. At the large sensor sites (L-sites), Atmospheric Surface Flux Stations (ASFSs) were deployed to measure the components of the SEB and meteorological variables (Cox et al., 2023). Broadband LWU and LWD were measured using pyrgeometers mounted at 2 m above the surface, with estimated uncertainty of 2.6 W m⁻² and 1 W m⁻² for LWD and LWU, respectively. Surface skin temperature was derived from longwave fluxes using an assumed surface emissivity of 0.985, which is appropriate for winter Arctic snow conditions. Two conductive heat flux plates were installed within the snowpack within several meters of each ASFS, providing estimates of the conductive flux at a single depth near the sea ice–snow interface. Three-dimensional wind speed was measured by ultrasonic anemometers, while turbulent fluxes were calculated from 10-Hz anemometer data using the eddy covariance technique. SH estimates have an estimated uncertainty of 4 W m⁻². All data used here from the ASFSs and tower have 10-min resolution.

As context for the analyses of snow and ice thermal conductivity and its interaction with the SEB, cloud measurements made from onboard Polarstern are also included. Radar reflectivity is derived from the vertically pointing Ka-band Atmospheric Radiation Measurement (ARM) Zenith Cloud Radar, operated by the U.S. Department of Energy ARM Program (Lindenmaier et al., n.d.).

We use data from 3 seasonal ice mass balance (SIMB) buoys deployed at the DN’s large L-sites (Perovich et al., 2022; Perovich et al., 2023). These SIMBs are spar buoys that included a temperature chain with 192 thermistors at 2-cm spacing and acoustic rangefinders above and below the ice that measure the top and bottom of the snow-ice system, used to determine the ice thickness (Polashenski et al., 2011; Planck et al., 2019). Data were collected at 4-h intervals. We use temperature and sea ice thickness measurements from the SIMB. At installation in October, sea ice thickness at the 3 SIMB sites ranged from 0.30 m to 1.20 m (Table 1), and snow depths are in the range of approximately 0.06–0.12 m. At installation, 0 cm refers to the nominal sea ice–snow interface.

Additional snow and sea ice temperatures are obtained from Snow and Ice Mass Balance Arrays (SIMBA), which are thermistor string-type IMBs (Jackson et al., 2013). We use 10 SIMBAs installed during October at the CO and as part of the DN at the L-sites, medium instrumentation sites (M-sites), and at position nodes (P-sites; Lei et al., 2022). SIMBAs consist of a 4.8-m-long thermistor chain with sensors at 2-cm intervals that are sampled at 6-h intervals. SIMBAs were installed vertically through the sea ice and snow, extending from the ocean to the atmosphere. In addition to temperature profiles, SIMBAs measure temperature differences after pulsative heating over 120 s. The ratio of temperature differences between 30 and 120 s after the heating is used to distinguish the sea ice–water interface and calculate sea ice thickness. At installation, sea ice thicknesses ranged from 0.34 m to 1.81 m, and snow depths ranged from 0.06 m to 0.18 m (Table 1). As with the SIMB, 0 cm refers to the sea ice–snow interface at installation. Locations of IMBs (generally referring to all SIMB and SIMBA), the CO, and ASFSs are shown in Figure 3.

To prevent the artificial gradient produced by the calibration offsets between thermistors from propagating into the calculations of thermal conductivity, temperature profiles are smoothed at each time step using a running mean window of 5 levels below 0 cm and 3 levels above. This smoothing is done for both SIMBA and SIMB data used in deriving thermal conductivity and density profiles, discussed next.

To derive thermal conductivity and conductive flux from the temperature profiles measured by IMBs, we solve the one-dimensional heat diffusion equation largely based on the work of Lipscomb (1998). Our derivation further provides the estimates of density profiles. Temperature profiles have been previously combined with forward models and additional constraints to solve for thermal conductivity. Brandt and Warren (1997) and Oldroyd et al. (2013) solved the inverse problem to derive thermal diffusivity in a seasonal snowpack and thermal conductivity at the South Pole. Marchenko et al. (2019) used a forward-model constrained by density measurements to derive the profiles of thermal conductivity on Svalbard.

For this derivation, we assume energy is only transported via vertical conduction through the sea ice-snow column, described by a one-dimensional heat diffusion equation:

where ρ is the density, c_p is the heat capacity, k is the thermal conductivity, z is the vertical spacing, t is the time, and T is the temperature. The left-hand side of this equation represents storage, and the right-hand side is the vertical divergence of conduction. Equation 1 conserves energy and explicitly states that any divergence of the vertical conductive heat flux must manifest itself as a temperature change at that depth. It is appropriate for nonmelting conditions. Equation 1 can be approximated with finite differencing at time step m and level n (n increasing downward, i.e., at the top level n = 0) as

This approximation assumes that density can change with height but does not vary between individual adjacent time steps. Heat capacity is assumed to only depend on temperature (Cuffey and Paterson, 2010). Equation 2 is written for each level between the top of the snowpack (n = 0) to within the sea ice and rearranged in matrix form to solve for k. For each time step, the lowest level N, where Equation 2 is written, is varied from −100 cm to −30 cm, with respect to the nominal sea ice–snow interface from installation, at every 2-cm interval, that is, 36 profiles if sea ice is thicker than 120 cm. At the lowest level N, k_ice is assumed to be 2 W m⁻¹ K⁻¹. We vary the depth of the lowest level because the thermal conductivity at a particular level may not be exactly 2 W m⁻¹ K⁻¹, but averaged over the bulk of the ice, it should be representative; 2 W m⁻¹ K⁻¹ is within the measured range of core sea ice thermal conductivity from Trodahl et al. (2001) and Pringle et al. (2007) for multiyear ice (about 1.8–2.3 W m⁻¹ K⁻¹). Measured thermal conductivity for first year ice is larger, with profile averages of 2.06–2.26 W m⁻¹ K⁻¹ (Pringle et al., 2007) and subsurface averages of 2.06–2.19 W m⁻¹ K⁻¹ (Pringle et al., 2006). If the sea ice is thinner than 120 cm, the lowest depth is set to 20 cm less than the sea ice thickness, so that knowledge of fluxes at the ocean–ice interface is not needed. We also exclude the bottom of the sea ice because it can have different properties than the bulk interior, such as a weak skeletal layer (Weeks, 2010) or brine overturning (Notz et al., 2005).

For each time step and lowest level, initially density is assumed to be constant in the sea ice and snow with values of 900 kg m⁻¹ and 300 kg m⁻³, respectively. The generalized inverse of the thermal conductivity matrix is calculated using its singular-value decomposition (SVD) and including all large singular values with the Python function numpy.linalg.pinv(). If there is no solution to the inverse problem, that is, the determinant is zero, the profile is flagged for removal later in the derivation.

After a first profile of k is derived, density is calculated using Equation 5 from Calonne et al. (2019) for densities below 450 kg m⁻³. For densities above 800 kg m⁻³, Equation 4 from Pringle et al. (2007) is used assuming salinity is 5 ppt. Varying salinity had no measurable impact on snow thermal conductivity. For densities between 450 and 800 kg m⁻³, the following weighted average is used:

where k_Calonne, thermal conductivity from Calonne et al. (2019), and k_Pringle, thermal conductivity from Pringle et al. (2007), are both functions of density. In Equation 3, $ϕ=(1+exp(−2*a*(ρ−450)))−1$⁠, and a = 0.02 m³ kg. Alternative density-conductivity relationships, for example, Yen (1981) and Macfarlane et al. (2023b), lead to differences in snow density estimates ≤30 kg m⁻³ or 10%.

Using the new density profile calculated from k, the generalized inverse of the thermal conductivity matrix is again solved. A density profile is again calculated from the new k profile, and this process is iterated 5 times in total. Profiles typically converge after 1 or 2 iterations, and profiles that do not converge tend to be excluded from the final results because at least 1 profile in the timestep has a determinant equal to zero. The final profiles of thermal conductivity and density for each time step are the averages over the 36 profiles calculated by starting at different depths. For time steps when there is at least 1 profile with no solution for the inverse problem, that is, the determinant is zero, thermal conductivity and density are linearly interpolated temporally over each level. A schematic of these steps can be found in Figure S2.

The same inversion method is used to derive thermal conductivity and density profiles in the bottom part of the sea ice column using an upper boundary condition. Again, the upper boundary is set to k_ice = 2 W m⁻¹ K⁻¹ for each of the multiple calculations with starting depths ranging from −30 to −100 cm. Thermal conductivity profiles are solved from the upper boundary to the first level above the sea ice bottom, and these profiles are appended below the profiles for the ice and snow that were solved using a lower boundary condition. At some time steps, thermal conductivity increases within several centimeters of the ice bottom above realistic values (>3 W m⁻¹ K⁻¹) (Figure S3). These large values could be due to the assumption of constant salinity, the sea ice density-thermal conductivity equation not being representative, or other fluxes (e.g., ocean heat flux) that are not included in Equation 1. Because we are primarily interested in snow and near surface fluxes, further refinement of this method near the sea ice bottom is left for future work.

Finally, conductive flux, F_C, is calculated with unsmoothed temperature profiles:

In the end, each profile of temperature measurements results in profiles of effective thermal conductivity, density, and conductive heat flux, all of which are related to each other via the described relationships.

Key interfaces within the vertical profile, that is, the atmosphere–snow and snow–ice interfaces, can manifest in the measurements in different ways. In particular, due to differences in thermal conductivity and density, temporal and vertical variations in temperature are markedly different for thermistors that are in different media. Additionally, the SIMB systems use sonic rangers to identify the top of the snow, but the SIMBA systems do not have these ranger measurements. To maintain consistent definitions of interface depths across different measurement systems, we identify the top of the snowpack by where the derived thermal conductivity increases beyond 0.6 W m⁻¹ K⁻¹. This threshold generally agrees well with the change in vertical temperature gradient between snow and air. Determining the interface between snow and sea ice is complicated by the presence of refrozen surface scattering layers or transition layers that have density and thermal conductivity between those expected of snow and sea ice. Because the snow–ice interface is not clearly distinguished for many IMBs, we determine the bottom of the snowpack as where the mean vertical profile of thermal conductivity is less than or equal to 0.6 W m⁻¹ K⁻¹, discussed further in Section 4.1.

There are various sources of uncertainty in these derivations of snow thermal properties, ranging from measurement uncertainties to the assumptions within the derivation. Individual thermistor measurements in thermistor strings have an uncertainty of ±0.5°C to −40°C ambient and ±0.25°C to −20°C ambient for SIMB (Planck et al., 2019) and ±0.1°C applicable to the full temperature range for SIMBA (Jackson et al., 2013; Lei et al., 2022). These uncertainties partly manifest as offsets for a given thermistor relative to others in the string. These offsets are of a similar magnitude as the true vertical temperature gradient and its second derivative through the sea ice; hence why we apply the vertical smoothing of temperature profiles when deriving thermal conductivity. Temperature gradients within snow are typically an order of magnitude larger than those in sea ice when the snow is not isothermal, so possible offsets are generally smaller than real temperature differences.

Due to variations in salinity and, to a lesser extent, temperature, the assumption of k_ice = 2 W m⁻¹ K⁻¹ at the lower boundary may not be true at any particular level or time step, and the relative error in this assumption propagates linearly into the derived k values throughout the profile. This assumption impacts the magnitude of derived thermal conductivities in snow but not the relative vertical variability. Typically, derived thermal conductivity profiles for different starting depths at the same time step have the same vertical profile but can be shifted relative to each other (Figure 4a). The standard deviations between profiles derived with different starting depths have median values of 0.021–0.052 W m⁻¹ K⁻¹ in the snowpack (Figure 4c). Varying the level of the lower boundary condition and taking the average at each time step is meant to minimize the impact of this issue, in the sense that if the bulk thermal conductivity of sea ice is 2 W m⁻¹ K⁻¹, the thermal conductivity in the sea ice from the different solved profiles averages to approximately 2 W m⁻¹ K⁻¹.

Subdaily variability of thermal conductivity in the snowpack is shown in Figure 4b. Standard deviations are calculated over each day, that is, across 6 or 4 consecutive time steps for SIMB and SIMBA, respectively, for 6 of the IMBs. Standard deviations range from near zero up to 0.18 W m⁻¹ K⁻¹ with medians around 0.02–0.06 W m⁻¹ K⁻¹ (Figure 4d). This variability is larger from SIMB instruments than SIMBA instruments (Figure S4). We believe this variability is primarily due to different scales between the temporal and vertical resolutions, meaning that there could be significant temperature changes between the time steps that are not adequately sampled by IMBs causing unrealistic variability. For example, several IMBs show temporary changes in thermal conductivity on the order of ±0.7 W m⁻¹ K⁻¹ around November 17, 2019, during a cyclone (Rinke et al., 2021). During this period, ASFS skin temperatures fluctuated by roughly ±10°C within hours. Neither 4 or 6 hourly observations fully capture these temperature fluctuations and their thermal impacts on the snow and sea ice.

It is difficult, maybe impossible, to confidently determine where exactly the snow surface is based on temperature profiles alone (Jackson et al., 2013), although various algorithms try (e.g., Zuo et al., 2018). For this reason, we focus on averages in the snowpack and “near-surface” values rather than values directly at the snow–atmosphere interface. This issue is a significant problem for estimating the surface conductive flux for SEB closure since it can impact the effective surface skin temperature that is consistent with the IMB measurements. Additionally, any potential biases in the IMBs, including those that are temperature dependent, can impact the surface temperature estimated from IMBs compared to that measured by other means. Comparing 7 IMBs closest to ASFS at the L-sites, temperatures at the highest level of resolved profiles (i.e., at the determined snow–atmosphere interface) are typically warmer than ASFS-measured skin temperatures (Figure 5a). Mean differences are +0.7°C–1.2°C; 13% of the data have differences larger than +2°C, and 1% of the data have differences larger than +4°C. Only 7% of data have colder temperatures than ASFS skin temperatures. While horizontal variability in snow-ice properties between the specific sampling sites certainly leads to some of the spread in differences, the statistics of these differences indicate that the IMB temperatures at the implied snow–atmosphere interface are on average warm relative to independent measures of the skin temperature. This apparent warm bias could be due to misidentification of the snow–atmosphere interface as well as a systematic bias in the IMB measurements that appears to be temperature dependent (Figure 5).

The uncertainty of the snow–atmosphere interface creates additional uncertainty for the near-surface conductive flux calculations. To assess this impact, we examine bulk estimates of conductive flux using the temperature difference between the snow–atmosphere interface, T_snow−atm, and the snow–ice interface, T_snow−ice:

where h_eff,snow is the effective snow depth:

with bounds at the snow–ice interface (n = snow-ice) and snow–atmosphere interface (snow-atm). Effective snow depth is the equivalent thickness of snow if it had the same thermal conductivity as a reference value (k_ref = 2 W m⁻¹ K⁻¹). Δz is the spacing between individual thermistors, 2 cm. The snow–ice interface level is defined as where mean thermal conductivity profiles are first less than or equal to 0.6 W m⁻¹ K⁻¹. Equation 5 is calculated twice using ASFS measured skin temperatures and IMB surface temperatures for T_snow-atm. The mean difference between bulk conductive fluxes calculated with these different temperatures is approximately 3 W m⁻², with larger differences at colder temperatures (Figure 5c and d). Thus, while not all details of this apparent IMB measurement bias are clear, this independent assessment suggests that it leads to a negative bias ranging from near 0 W m⁻² at relatively warm temperatures to an average of about 5 W m⁻² at the coldest observed temperatures.

We have quantified various sources of uncertainties in this section; however, it is not a simple task to combine them into a single estimate of uncertainty. Because the individual uncertainty estimates are not independent of each other, they cannot be added in quadrature. The different uncertainties in snow thermal conductivity we have discussed each have similar magnitudes of approximately 0.02–0.06 W m⁻¹ K⁻¹ or approximately 10%–15%. A 10%–15% uncertainty in thermal conductivity would mean an uncertainty in conductive flux around 2–3 W m⁻², on average. The uncertainty due to temperature measurements in the thermal conductivity estimates also contributes to uncertainty in conductive flux. The uncertainty associated with the surface location and SIMB temperature bias is on average 3 W m⁻². Overall, we estimate the conductive flux uncertainty to be around 5 W m⁻² and the thermal conductivity uncertainty to be 0.05 W m⁻¹ K⁻¹.

In this section, we first quantify thermal conductivity in snow and its variability vertically, spatially, and temporally. Next, we assess the conductive flux near the surface as it relates to other, independently measured, surface fluxes. We further determine how environmental factors drive variability in conductive flux through the snowpack, for example, sea ice thickness, snow depth, and atmospheric state. These observation-based relationships between the environment and conductive fluxes are finally used to assess output from CAFS.

Figure 6a and c shows thermal conductivity and density at SIMBA 2019T62, located at the CO near the MOSAiC second year ice coring site. The sea ice grew from 0.80 to 1.48 m over November 1, 2019, to March 1, 2020, and the mean estimated snow depth was 19 ± 2 cm. No flooded sea ice has been identified at this site during this period (Lei et al., 2022). Above 0 cm, the estimated snow–ice interface when the SIMBA was installed, thermal conductivities range from 0.06 to 1.43 W m⁻¹ K⁻¹, with a corresponding density range of 142–721 kg m⁻³. Snow with the lowest density and thermal conductivity (<0.2 W m⁻¹ K⁻¹) appears at the top of the snowpack, associated with time periods where the snowpack increases in depth, for example, January 15–25, 2020, that are consistent with fresh snow (Matrosov et al., 2022). Thermal conductivity near the snow–ice interface increases from less than 0.75 W m⁻¹ K⁻¹ at +2 cm depth to greater than 1.5 W m⁻¹ K⁻¹ below −4 cm depth. Some noise is also apparent in Figure 6a. Notably, on November 17, 2019, there are sudden, temporary changes around 0–10 cm depth. This episode is during a cyclone (Rinke et al., 2021) that caused warmer surface temperatures and potentially changes to temperatures that were poorly sampled by the SIMBA, as discussed in Section 3.3.

SIMBA 2019T63 at L2 shows a different vertical structure than SIMBA 2019T62. Sea ice grew from 1.12 to 1.56 m at L2 during the considered period, and snow depth was 10 ± 2 cm from November through January and increased to 19 cm after February 1, 2020. SIMBA 2019T63 (Figure 6b and d) has a lower range of snow thermal conductivity and density compared to SIMBA 2019T62, 0.10–0.55 W m⁻¹ K⁻¹ and 180–503 kg m⁻³ above 0 cm. Excluding time periods with noise due to cyclones (e.g., November 17, 2019), thermal conductivity varies less vertically within the snowpack than at SIMBA 2019T62, although low thermal conductivity values also appear near the top of the snow. SIMBA 2019T63 additionally shows a more gradual transition between snow and sea ice, at least in terms of density and thermal conductivity values. Mean thermal conductivity increases from 0.31 W m⁻¹ K⁻¹ at 0 cm to 1.5 W m⁻¹ K⁻¹ at −20 cm, whereas at SIMBA 2019T62, thermal conductivity increases the same amount in about half the thickness.

Different vertical structures are apparent between IMBs when comparing mean vertical profiles of thermal conductivity (Figure 7a and b). SIMB L3 has the largest thermal conductivity throughout the snowpack (approximately 0.45 W m⁻¹ K⁻¹). Several SIMBAs (2019T58, 2019T62, 2019T67, 2019T63, and 2019T65) all have minima in thermal conductivity within the middle of the snowpack with lowest values just above 0.2 W m⁻¹ K⁻¹. 2019T62 (solid blue) varies the most within the snowpack, from approximately 0.4 W m⁻¹ K⁻¹ at 2 cm to approximately 0.2 W m⁻¹ K⁻¹ at 12 cm. Three SIMBAs (2019T47, 2019T64, and 2019T62) show larger thermal conductivity in the topmost layer. These values could be due to wind slab formation or the result of how the top of the snowpack is determined using a threshold of k_eff = 0.6 W m⁻¹ K⁻¹.

Thermal conductivity profiles from all IMBs have some “transition layer” at least 6 cm thick between the snow and the sea ice, where mean thermal conductivity values are between 0.7 W m⁻¹ K⁻¹ and 1.5 W m⁻¹ K⁻¹. Ten of the IMBs have transition layers at least 8 cm thick, and 3 of those 10 have a transition layer that is 14–16 cm thick. IMBs installed in first year ice (SIMBAs 2019T47 and 2019T70 and SIMB L3) have thinner transition layers (6 or 8 cm), while the thickest transition layers are all observed in second year ice. Past observations found sea ice cores with lower densities and therefore lower expected thermal conductivities, near the top of sea ice due to granular ice, transition layers, and/or the high desalination above the sea level (e.g., Eicken et al., 1995; Huang et al., 2013). Surface scattering layers were also observed at MOSAiC (Macfarlane et al., 2023a), and these transition layers could be remnants of the surface scattering layer from the prior melt season for the second year ice. Due to vapor diffusion and other processes, some sort of transition layer is expected between the snow and ice, and we have no reason to believe the observed transition layers are an artifact of the measurements or processing. Moving forward, we focus on thermal conductivity and conductive fluxes within the snowpack only, defined based on where the mean vertical profile of thermal conductivity is less than or equal to 0.6 W m⁻¹ K⁻¹ (Figure 7a and b).

Spatial differences in snow thermal conductivity are typically larger than temporal variability at a given location. Figure 7c shows the 2-week running means of median thermal conductivity above the snow–ice interface. For most IMBs, snow thermal conductivity varies in time by approximately 0.05–0.1 W m⁻¹ K⁻¹. Local minima, such as around November 25, 2019, for IMBs at L2 (orange solid and dashed lines), are due to higher derived snow top estimates with lower thermal conductivities consistent with snow accumulation. When the upper snow is redistributed, compressed, or undergoes other transformations, the median snow thermal conductivity increases again. Over the full winter, 7 IMBs show increasing trends, consistent with Macfarlane et al. (2023b). Three of the IMBs show overall decreases in snow thermal conductivity, and the remaining 3 have no clear trend and/or large variability.

Median snow thermal conductivity values can vary as much between sites as within a given site. In Figure 8, values are included for levels above the sea ice–snow interface (Figure 7) that are also less than 0.8 W m⁻¹ K⁻¹. This threshold is used to exclude time periods with synoptic events that temporarily and unrealistically increase the apparent thermal conductivity. How the sea ice–snow interface is defined does impact the magnitude of the following statistics but not the relative relationships between IMBs or changes over time. Over winter, median snow thermal conductivity ranges from 0.33 W m⁻¹ K⁻¹ (SIMBA 2019T65 at L2) to 0.47 W m⁻¹ K⁻¹ (SIMB at L3). At L1 (red) median, thermal conductivities are 0.39 and 0.41 W m⁻¹ K⁻¹. The SIMB at L2 has a similar median as IMBs at L1, 0.39 W m⁻¹ K⁻¹ but is significantly higher than the SIMBAs at L2 that have the medians of 0.34 and 0.33 W m⁻¹ K⁻¹ (orange). The SIMB at L3 median is also larger than the SIMBA at L3, 0.47 compared to 0.41 W m⁻¹ K (yellow). At a given L-site, the individual systems are typically within tens of meters of each other. SIMBAs at 3 of the 4 M-sites (light blue) have medians that have a similar range, 0.38–0.43 W m⁻¹ K⁻¹, despite being kilometers apart. The largest interquartile ranges are from SIMB measurements. Over all sites, the median is 0.39 W m⁻¹ K⁻¹ with an interquartile range of 0.34–0.44 W m⁻¹ K⁻¹ (black).

In this section, we focus on conductive flux and its temporal and spatial variability during winter and among sites. First, a case study is used to highlight the relationships between SEB terms, their variability, and the difficulty of trying to close the SEB. Next, we investigate variability due to sea ice thickness, snow depth, and temperature and quantify how differences in snow impact the coupling between conduction and the atmosphere. Finally, we use these relationships to reassess conductive flux in CAFS.

A case study from November 18 to December 13, 2019, shows the temporal variability of conductive flux and its relationships to other surface fluxes (Figure 9). The first half of this period is dominated by clouds, with some intermittent cloud-free times around November 25, 28, 30, and December 1–2, 2019 (Figure 9a). The atmosphere is mostly cloud-free after December 9. The changes in LWNET (blue) and skin temperature (red) follow cloud cover: LWNET and skin temperature are most negative during cloud-free time periods, for example, November 27 and December 2, 2019 (Figure 9a and c). When the atmosphere transitions to being cloudy, LWNET increases closer to 0 W m⁻² and the surface warms by approximately 10°C. SH (orange in Figure 9c) frequently follows changes in skin temperature, although high wind speed can also increase it. SH is directed into the surface during cloud-free and cold periods (December 1–2) and out of the surface during cloudy and warm periods (December 3–5), both on the order of tens of W m⁻². LH (green in Figure 9c) is less than 5 W m⁻² except for December 3–5.

Conductive flux from 7 IMBs that are within tens of meters of the ASFS are plotted in Figure 9d (pink). As previously discussed, the precise location of the top of the snowpack is not known from the IMBs and is difficult to extrapolate from the ASFS snow depth measurements due to strong spatial variability in snow depth. Therefore to be certain the flux is within the snow, we plot conductive flux from 2 levels (4 cm) below the estimated snow top from the thermal conductivity derivation. Conductive flux from the 2 flux plates installed just above the snow–ice interface at each ASFS is also shown in black. Comparing conductive flux estimates from IMBs and flux plates is complicated because of the constant height of the flux plates relative to the sea ice–snow interface and the spatial heterogeneity of the snow. As snow accumulates over the flux plates, the measured conductive flux is reduced, while a storage term should increase. Relating the change in conductive flux to snow depth is complicated, though, because there was only 1 snow depth measurement at each site, while the 2 flux plates can vary by a factor of two. Despite these caveats, conductive flux estimates from IMBs and flux plates in Figure 9d have similar magnitudes and temporal variability, for example, November 18–24. By the end of winter, conductive flux measurements from the different flux plates have smaller magnitudes and temporal variability compared to the IMB estimates (not shown).

The conductive flux responds to changes in the atmosphere, both in terms of magnitude and spatial variability, since it partly depends on the surface temperature. Higher surface temperatures, coincident with cloud occurrence, reduce the temperature gradient through the snowpack and therefore reduce the conductive flux. For example, during the transition from cloudy to cloud-free conditions around November 20, the conductive flux from IMBs increases from near 0 W m⁻² to approximately 15–40 W m⁻². Variability among flux plates and IMBs also tends to be larger when the surface is cold (November 26–December 1) and smallest when the surface is warm and the atmosphere is cloudy (November 19 and 20). The flux plates at ASFS 50 differ by upward of 40 W m⁻² during colder time periods (November 26–December 1) despite being within meters of each other. The differences between flux plates that are meters apart are often larger than differences between IMBs that are kilometers apart, consistent with Sturm et al. (2001). This local variability also demonstrates the modulating effects of sea ice and snow on conduction, although it complicates comparisons between instruments.

SEB closure from surface-based observations, such as those employed here, is a well-recognized challenge (i.e., Foken, 2008). There are multiple reasons why the energy budget from these observations may not fully close, including different scales of spatial variability for each term, different footprints of influence for each measurement, different time resolutions for the measurements, and different locations and snow-ice properties between IMBs and ASFSs. Unsurprisingly, individual IMBs rarely balance the net atmospheric flux (Figure 9d). Before November 23, conductive flux estimates from IMBs are within a few W m⁻² of the net atmospheric flux, but the conductive flux is generally smaller than the net atmospheric flux. The significant heterogeneity of snow makes it difficult to compare point measurements of conductive flux, either from IMBs or flux plates, to the footprint of the atmospheric fluxes, which range from a few meters to a few hundred meters depending on parameters and conditions. Therefore, we also compare averages. The pink shading in Figure 9e shows the range of conductive flux from IMBs within 6-h windows, and the solid pink line gives the mean. The cyan shading and line are ranges and mean from ASFS net atmospheric flux, multiplied by −1, during the same 6-h windows; 6-h averages are also calculated from the flux plates, shown in black, to show what another instrument measures even though we know a storage term should be included for them to balance the SEB. At the beginning and end of the time period, the range of conduction overlaps the net atmospheric flux, but the mean conductive flux is still biased low by 10 W m⁻² at most. Means from the IMBs and flux plates agree fairly well during these times. Over November 23 through November 30, mean conductive flux is approximately 20–30 W m⁻² smaller than the net sum of atmospheric fluxes. This time period has largest difference between IMBs and ASFS fluxes during the winter. For most of the winter at least, the ranges from IMBs and ASFS net atmospheric flux estimates overlap.

Overall, the residual from 6-h mean IMB conduction and ASFS net atmospheric flux has a mean of −6 W m⁻² and standard deviation of 7 W m⁻² over the entire winter (Figure 10a). It is unsurprising that there is a bias considering the imprecise knowledge from IMBs, including uncertainty in the thermistor measurements, the uncertainty associated with the location of the surface, and the general aggregate uncertainty in deriving thermal conductivity. Additionally, the atmospheric fluxes carry significant uncertainty themselves (Cox et al., 2023), which, when combined in quadrature, result in an aggregated uncertainty for the net atmospheric flux of approximately 5 W m⁻². The residuals tend to be larger during radiatively clear time periods, with large negative LWNET, and during periods of high wind (Figure 10b and c). Over the entire winter, the average residual during radiatively clear times (LWNET < −20 W m⁻²) is −8 W m⁻² but only −3 W m⁻² during radiatively cloudy times (LWNET > −20 W m⁻²). Although there is uncertainty in IMB and ASFS measurements, the conductive flux being too small during cold time periods would be consistent with the potential temperature-dependent bias in thermistors discussed in Section 3.3. Some disagreement between IMBs and ASFS could be due to spatial variability, which is generally greater during cold time periods, as discussed more in the following sections. The association between larger wind speeds and residuals could mean that energy is transported at the surface by means other than vertical conduction and/or that the atmospheric turbulent heat flux estimate has a larger uncertainty. However, further investigation into closing the SEB is beyond the scope of this work. While the average residual is nonzero, the fact that it is within the combined uncertainties from atmospheric and conductive fluxes provide some confidence that the IMB estimates are physically reasonable.

Given the large variability, both spatially and temporally, seen in Figure 9d, we now turn to exploring what environmental variables drive that variability. Because snow is a stronger thermal insulator than sea ice but physically thinner, we look at the relationships between conduction and sea ice and snow thicknesses as well as effective depth (h_eff). Effective depth is the equivalent thickness of snow and sea ice scaled by a reference value (k_ref = 2 W m⁻¹ K⁻¹):

where k_n are the thermal conductivities over each 2-cm layer Δz that are summed over the portion of the column for which thermal conductivity is derived above −30 cm. The upper bound is the snow–atmosphere interface “snow–atm.” Depths are relative to the 0 cm defined at IMB installation. The right-hand term represents the depth of sea ice below the snowpack and possible transition layer (h_ice) using the average thermal conductivity for the bulk sea ice at and below −30 cm (⁠$k¯ice$⁠), generally close to 2 W m⁻¹ K⁻¹. Effective snow depth is calculated in the same manner as effective depth but only for levels within the snowpack (Equation 6).

Effective depth is useful for considering the thickness impacts of snow compared to sea ice. For example, 10 cm of snow with k_snow = 0.3 W m⁻¹ K⁻¹ is equivalent to 67 cm of sea ice with k_ice = 2.0 W m⁻¹ K⁻¹, in terms of how readily each medium conducts energy. We further use effective snow depth on its own to compare the sensitivity of conduction to sea ice and snow depths. In Figure 11, snowpack-average conductive fluxes are plotted against snow depth, effective snow depth, sea ice thickness, and total effective depth, and colored by surface temperature. Best fit lines are calculated for 3 temperature ranges shown in the legends of each subplot; 27% of the data falls within the coldest temperature range, 17% is in the warmest temperature range, and 55% is in the middle temperature range. Using temperature ranges that evenly split the percentage of data in each temperature range results in the same basic relationships but different magnitudes.

Thicker snowpack and warmer surface temperature reduce the variability of conductive flux. Differences in snow depth, both physical and effective, become less important at warmer temperatures because the temperature gradient is reduced. From linear regressions (solid lines in Figure 11a–c), the sensitivity of conductive flux to physical snow depth decreases from −0.81 W m⁻²/cm for the coldest temperatures (−40°C to −28°C) to −0.47 W m⁻²/cm or the warmest temperatures (−17°C to −4°C). The sensitivity to effective snow depth shows a similar decreases, by over a factor of two, from −0.11 to −0.04 W m⁻²/cm. Put another way, differences in snowpack thickness result in more spatial variability in conductive flux when the surface is coldest. In addition to linear regressions, exponential fits are shown with cyan dashed lines in Figure 11a–f. Exponential fits may be more representative of the relationship between conduction and snow depth than a linear regression, given that the analytical solution of the heat equation using Fourier series is exponential. Exponential fits are also more physically intuitive since conduction will not become increasingly negative for thicker snowpacks. This reasoning appears true for the middle and warmest temperature ranges; however, the exponential fit for the coldest temperature range still appears linear, reiterating the point that conductive fluxes in the coldest temperature range have higher sensitivity, that is, a steeper curve, to effective snow depth than the warmest temperatures. The exponential fits for the medium and warmest temperature ranges approach asymptotes around 5–10 W m⁻², whereas in the case of infinitely thick snow, we might expect the conductive flux to approach 0 W m⁻². The nonzero asymptote could be due to the relatively thin sea ice below that still allows some conduction to occur. We can also consider the inverse relationship: how snow thickness influences the sensitivity of conduction to surface temperature. As snow depth increases, the range of conductive fluxes decreases (Figure 11a). Similarly, for effective snow depth less than 100 cm, the standard deviation of conductive flux is 10.4 W m⁻², whereas for effective snow depth thicker than 200 cm, the standard deviation is only 3.7 W m⁻², less than half the variability (Figure 11d–f).

In contrast to the decreasing sensitivity of conduction to snow depth with temperature, conductive flux in snow is about equally sensitive to sea ice thickness in each temperature range, −0.09 to −0.11 W m⁻²/cm (Figure 11g–i). This consistency could be due to the fact that sea ice grows as the surface temperature overall decreases throughout the winter. For the coldest temperature range, conductive flux is about as sensitive to sea ice thickness as effective snow depth, but at warmer surface temperatures, conductive flux is more sensitive to sea ice thickness. This sensitivity could drive the variability between IMBs during warm periods in Figure 9d. Because conduction is still sensitive to sea ice thickness when the surface temperature is warm, the sensitivity of conduction to effective depth is not quite as varied between temperature ranges as the sensitivity to effective snow depth. Conduction is less than half as sensitive to effective depth over the warmest temperatures (−0.03 W m⁻²/cm) as compared to the colder temperatures (−0.08 W m⁻²/cm; Figure 11j–l).

Conductive flux is coupled to atmospheric forcing through the surface temperature, which is driven by the balance of longwave radiation, turbulent heat fluxes, and conduction. The slope from the line of best fit from regressing conductive flux against LWD (Figure 12a) represents the sensitivity of conductive flux to winter radiative forcing. For each additional W m⁻² of LWD, conductive flux decreases by 0.22 W m⁻². For comparison, an additional W m⁻² of LWD increases LWU by approximately 0.5 W m⁻² (not shown). The sensitivity of conduction to LWD is actually the same regardless of effective snow depth, although the range is larger for thinner effective snow depths than thicker effective snow depths.

Previous work has pointed to the bimodal winter radiative states of the Arctic sea ice system, which are largely defined by the presence or absence of clouds (e.g., Shupe and Intrieri, 2004; Stramler et al., 2011). This bimodality is evident in the distribution of LWNET (Figure 12b), with −25 W m⁻² serving as the approximate threshold between radiatively clear and cloudy conditions. The joint histogram between LWNET and conductive fluxes also has distinct clear and cloudy clusters (Figure 12d). Despite the distribution of conductive fluxes appearing unimodal (Figure 12e), mean conductive fluxes in clear versus cloudy states are statistically different: 25 W m⁻² (10 W m⁻² standard deviation) for radiatively clear conditions and 14 W m⁻² (10 W m⁻² standard deviation) for radiatively cloudy times. The difference in conduction between clear and cloudy states has a slight dependence on effective snow depth. The difference between radiatively clear and cloudy states is slightly larger for thinner snowpacks than thicker snowpacks. For thinner effective snow depths (<92 cm), the means are 29 W m⁻² versus 16 W m⁻² for radiatively clear versus cloudy. The means drop to 20 W m⁻² versus 11 W m⁻² for thicker effective snow depths (>92 cm). These differences imply that thicker snow modestly dampens the coupling between the atmosphere and sea ice.

Given these relationships between conductive flux and the environment, we now return to CAFS to reassess how it compares to the observations. CAFS data are analyzed from a simulation over the MOSAiC year with output from grid cells nearest the Polarstern and ASFSs. For this run, the first 48 h of each daily forecast are used. Fluxes are positive toward the surface. For comparisons, near-surface conductive flux estimates are taken from IMBs, and data from CAFS and observations are averaged to daily timescales.

The relationships between sea ice thickness, snow depth, and conductive flux differ between CAFS and observations. In CAFS, sea ice thickness ranges from 161 to 254 cm over the winter (blue in Figure 13b). These values are generally larger with less variability over the winter compared to observations that vary from 50 to 208 cm (orange in Figure 13b). Snow depths in CAFS are within the observed range (4–29 cm) but have a narrower spread (7–20 cm). Conductive flux in CAFS depends less on both sea ice thickness and effective snow depth than in observations (Figure 13a). In observations, the thicker the sea ice or effective snow is, the less conduction. However, in CAFS even the thickest sea ice (approximately 250 cm) experiences the full range of conductive flux (about −10 to 60 W m⁻² in Figure 13d). Sea ice and snow do not appear to modulate conductive flux in CAFS in the same way that they do in observations.

Differences between observed and modeled conductive fluxes are also tied to atmospheric state. In CAFS, daily average conductive fluxes span a larger range (−10.9 to 62.4 W m⁻²) than observations (−3.7 to 48.0 W m⁻²). The larger conductive fluxes in CAFS are associated with colder surface temperatures. In Figure 13g, CAFS has more occurrences of surface temperatures below −30°C compared to surface temperatures from IMBs and skin temperatures from ASFS. These colder surface temperatures occur during cloud-free periods, which CAFS has slightly more occurrences of compared to observations (Figure 13f). This difference suggests that CAFS simulates an optically thinner atmosphere during clear-sky conditions. CAFS also has more instances of conduction going downward from the surface, that is, negative conductive fluxes, than in the observations; 4% of conductive fluxes in CAFS are less than 0 W m⁻² compared to 1% from IMBs and less than 1% from the flux plates. In observations, downward conduction at the surface typically only occurred during cyclones when the atmosphere was warmer than the surface. If the surface is more often colder than the atmosphere, negative conduction is more likely to occur. These connections between conduction, surface temperature, and the atmosphere underscore the importance and complications of modeling the Arctic as a coupled system.

Our thermal conductivity values in snow are within the range of past observations and other measurements from MOSAiC, albeit on the high end. Including all measurements leads to a median snow thermal conductivity of 0.39 W m⁻¹ K⁻¹, which is higher than the value used by large scale models, although it is sensitive to how the sea ice–snow interface is defined. Macfarlane et al. (2023b) measured a median snow thermal conductivity of 0.26 ± 0.05 W m⁻¹ K⁻¹ using a finite element method and microCT scans of snow samples at MOSAiC. They defined snow as samples with density less than 550 kg m⁻³. From our estimates, none of the 13 IMBs have median values in snow less than 0.3 W m⁻¹ K⁻¹, and our range extends above the maximum of 0.5 W m⁻¹ K⁻¹ in Macfarlane et al. (2023b). Our larger estimates could, in part, be due to differences in defining snow and processes are included in the “effective” definition of thermal conductivity. Macfarlane et al. (2023b) consider processes that contribute to vertical conduction, while our methods assume any measured temperature changes are due to vertical conduction. Perovich et al. (2023) also derived snow thermal conductivity from the SIMBs at MOSAiC around 0.39–0.48 W m⁻¹ K⁻¹ by assuming flux continuity at the sea ice–snow interface. Our median estimates from the same SIMB instruments are similar; however, including results from additional IMBs increases the spread of thermal conductivity. Future work should be done to synthesize the different thermal conductivity estimates from MOSAiC.

A key point of this work is that the representation of snow thermal conductivity is important for modeling the SEB and the lower atmosphere. In agreement with past work (Wu et al., 1999; Blazey et al., 2013), changing snow thermal conductivity in a coupled model also changes the surface temperature and balance of surface fluxes. While we are not the first to acknowledge this, the representation of snow has received far less consideration in atmospheric modeling as compared to sea ice modeling. To wit, not all forecast models even include snow on sea ice (e.g., Bauer et al., 2016; Bengtsson et al., 2017). Given the observed variability of snow thermal conductivity, future studies should test the extent to which spatial and temporal variability of snow thermal conductivity is important in coupled models.

Although conductive flux estimates are often within measurement uncertainties for closing the SEB, future work can improve upon SEB closure. The estimates of conductive flux at the snow–atmosphere interface from thermistor string measurements are uncertain due to both the challenge in identifying the surface location as well as potential temperature-dependent biases in the thermistors. If the thermistors have larger positive biases at colder temperatures (as suggested by Figure 5), the temperature gradients during these time periods would be biased low. Underestimated temperature gradients would cause underestimated conductive flux, consistent with our conductive flux estimates being typically smaller than the net atmospheric flux. Future studies could further assess, quantify, and correct this potential bias in the temperature profile data.

We have estimated effective snow thermal conductivity during October–February over the Central Arctic sea ice during the MOSAiC expedition using temperature profiles from sea ice mass balance buoys. Our method provides vertical profiles of thermal conductivity, density, and conductive flux that are continuous in time. Median values from 13 IMBs range from 0.33 to 0.47 W m⁻¹ K⁻¹, with spatial variability greater than the estimated 0.05 W m⁻¹ K⁻¹ uncertainty. The median from all available data during October through February is 0.39 W m⁻¹ K⁻¹, and the interquartile range is 0.34–0.44 W m⁻¹ K⁻¹. The temporal and vertical variability of snow thermal conductivity and density are consistent with environmental changes, such as new snowfall.

Estimated thermal conductivity was used to calculate conductive flux through the snowpack. Conductive fluxes are larger and more spatially variable during cold and radiatively clear periods as compared to warmer and cloudier times. Typically, the derived conductive flux is smaller than the observationally estimated net atmospheric flux, leading to a mean residual of −6 W m⁻² over the winter, with larger residuals during cold, cloud-free periods and diminished residuals during warm, cloudy periods. The mean residual is on the same order of magnitude as uncertainties in the conductive flux and the net atmospheric flux (each approximately 5 W m⁻²), but there can still be considerable residuals at any specific time.

Using a coupled regional forecast model, we found that varying snow thermal conductivity impacts the balance of all surface heat fluxes. Importantly, the snow thermal conductivity modulates the interplay between atmospheric radiation and turbulent heat fluxes. Detailed and time-variable thermal conductivity and conductive flux measurements can be coupled with the time-varying estimates of atmospheric surface flux terms to better explore the processes that help to drive the surface temperature and the overall SEB. In both observations and models, coupling the atmosphere and surface is essential for linking atmospheric variability to the thermal state of the snow and sea ice as well as to sea ice growth in winter. Future work should further explore these processes for model assessment and improvement.

Coupled Arctic Forecast System output is being archived on PANGAEA. Atmospheric Surface Flux Stations and Met City measurements are available from https://doi.org/10.18739/A2PV6B83F, https://doi.org/10.18739/A2XD0R00S, https://doi.org/10.18739/A25X25F0P, and https://doi.org/10.18739/A2FF3M18K. SIMB data are available at https://doi.org/10.18739/A2SB3X05G. SIMBA data are available from PANGAEA at the following: https://doi.org/10.1594/PANGAEA.940393 (2019T58), https://doi.org/10.1594/PANGAEA.940387 (2019T47), https://doi.org/10.1594/PANGAEA.940231 (2019T62), https://doi.org/10.1594/PANGAEA.940593 (2019T63), https://doi.org/10.1594/PANGAEA.940617 (2019T64), https://doi.org/10.1594/PANGAEA.940634 (2019T65), https://doi.org/10.1594/PANGAEA.938128 (2019T67), https://doi.org/10.1594/PANGAEA.940650 (2019T68), https://doi.org/10.1594/PANGAEA.940659 (2019T70), and https://doi.org/10.1594/PANGAEA.940668 (2019T72). Derived snow and ice thermal properties are available from https://doi.org/10.18739/A2CF9J80P.

The supplemental files for this article can be found as follows:

Supplemental material is available in sledd_elementa_supplemental_20231207.pdf.

Data sets used in this article were produced as part of the International Multidisciplinary drifting Observatory for the Study of the Arctic Climate (MOSAiC) with the tag MOSAiC20192020 and the Project_ID: AWI_PS122_00. Some data used here were collected by the Atmospheric Radiation Measurement User Facility, a DOE Office of Science user facility operated by the Biological and Environmental Research program. The authors thank all people involved in the expedition of the Research Vessel Polarstern (Alfred-Wegener-Instut, Helmholtz-Zentrum für Polar-und Meeresforschung, 2017) during MOSAiC in 2019–2020 as listed in Nixdorf et al. (2021). They additionally thank Darren Jackson for feedback on the manuscript.

This work was supported by DOE Atmospheric System Research Program (DE-SC0021341), the National Science Foundation (OPP-1724551), NOAA’s Physical Sciences Laboratory (PSL; NOAA Cooperative Agreement NA22OAR4320151), and by NOAA’s Global Ocean Monitoring and Observing Program (GOMO; FundRef https://doi.org/10.13039/100018302). RL is funded by the Program of Shanghai Academic/Technology Research Leader (22XD1403600). DP was supported by NSF-OPP 1724424, NSF–OPP–2034919, and NSF–OPP–2138785.

The authors declare that they have no conflict of interest.

Contributed to conception and design: AS, MDS, AS, CJC.

Contributed to acquisition of data: MDS, AS, CJC, DP, RL.

Contributed to analysis and interpretation of data: AS, MDS, AS, CJC, DP.

Drafted and/or revised the article: All authors.

Approved the submitted version for publication: All authors.

How to cite this article: Sledd, A, Shupe, MD, Solomon, A, Cox, CJ, Perovich, D, Lei, R. 2024. Snow thermal conductivity and conductive flux in the Central Arctic: Estimates from observations and implications for models. Elementa: Science of the Anthropocene 12(1). DOI: https://doi.org/10.1525/elementa.2023.00086

Domain Editor-in-Chief: Detlev Helmig, Boulder AIR LLC, Boulder, CO, USA

Guest Editor: Zoe Courville, US Army Engineer Research and Development Center Cold Regions Research and Engineering Laboratory, Hanover, NH, USA

Knowledge Domain: Atmospheric Science

Part of an Elementa Special Feature: The Multidisciplinary Drifting Observatory for the Study of Arctic Climate (MOSAiC)