The importance of sub-meter-scale snow roughness on conductive heat flux of Arctic sea ice

Abstract

The conductive heat flux through the snow and ice is a critical component of the mass and energy budgets in the Arctic sea ice system. We use high horizontal resolution (3–15 cm) measurements of snow topography to explore the impacts of sub-meter-scale snow surface roughness on heat flux as simulated by the Finite Element method. Simulating horizontal heat flux in a variable snow cover modestly increases the total simulated heat flux. With horizontal heat flux, as opposed to simple 1D-vertical heat flux modeling, the simulated heat flux is 10% greater than that for uniform snow with the same mean snow thickness for a 31.5 × 21 m region of sea ice (the largest region we studied). Vertical-only (1D) heat flux simulates just a 6% increase for the same region. However, this is highly dependent on observation resolution. Had we measured the snow cover at 1 m horizontal spacing or greater, simulating horizontal heat flux would not have changed the net heat flux from that simulated with vertical-only heat flux. These findings suggest that measuring and modeling snow roughness at sub-meter horizontal scales may be necessary to accurately represent horizontal heat flux on level Arctic sea ice.

1. Motivation and introduction

The decline of Arctic sea ice is a dramatic consequence of and contributor to climate change (e.g., Screen and Simmonds, 2010). However, the observed decline is not fully captured by climate models (Notz and SIMIP Community, 2020). Of the many uncertainties in the sea ice component of climate models, our lack of understanding of how much heat flows through snow and sea ice in winter is one of the greatest uncertainties impacting the amount of Arctic sea ice that models simulate (Urrego-Blanco and others, 2016). The spatial distribution of snow on sea ice is a factor governing this heat flow.

Wind-driven snow redistribution produces a snow cover on Arctic sea ice in which the spatial variability in snow depth is a substantial fraction of the mean, even on level ice. Mallett and others (2022) analyzed 499 snow depth transects (each at least 500 m long) on multi-year ice collected from 1955–1991 and found that the average coefficient of variation, the ratio of the standard deviation to the mean, was 0.417. We tabulated this metric for an additional 24 sites of level, first-year ice (Iacozza and Barber, 1999; Sturm and others, 2002a; Petrich and others, 2012; Webster and others, 2015; Moon and others, 2019) and found an average coefficient of variation of 0.42. Individual sites varied from 0.23 to 0.94. The mean snow depth across these sites was 20.0 cm and the mean standard deviation was 8.5 cm.

Due to fundamental thermodynamics (Fourier, 1822; Sturm and others, 2002b), small-scale (i.e., centimeters to 100 s of meters) spatial variability in snow depth on Arctic sea ice enhances the total amount of conductive heat flux through the snow and ice covers compared to what the conductive heat flux would be for a uniform snow cover with the same mean snow thickness. Conceptually, we can divide this heat flux enhancement into two components based on how we approximate heat flux in models: that which occurs in a ‘vertical-only’ 1D heat flux simulation and that which is simulated with the full 3D heat fluxes–which we term ‘including horizontal’ heat flux. Vertical-only simulation is well known and has been implemented in models (e.g., Abraham and others, 2015). In contrast, just two published works (Sturm and others, 2002b; Popović and others, 2020) have explored including horizontal heat flux simulation on Arctic sea ice, with potentially contradictory findings (discussed below). However, to our knowledge there are no published assessments of the relevant spatial resolution for horizontal heat flux simulation. We address this gap through heat flux simulations using 3–15 cm horizontal resolution measurements of snow surface topography.

Snow is ~10 times more thermally insulating than sea ice (Sturm and others, 2002b). The reason horizontal heat flux is potentially important is that it can reduce the net thermal resistivity of the snow-ice system by transporting heat laterally through the relatively conductive ice to areas of thinner snow (at the cost of a longer transport path length in the ice). Consider a thought experiment: you are at the base of a uniform, 1-m-thick sheet of ice with 20 cm of snow directly above you. You have to choose between heat transport pathways that go directly vertically above you and diagonally through the ice and then through a thinner snow cover. How much thinner does the snow on the diagonal pathway have to be for the diagonal path to be preferable (lower total thermal resistivity) to the vertical pathway? At small horizontal distances, the diagonal distance through the ice is barely greater than the vertical distance. So even small reductions in snow thickness would result in a preferable heat transport pathway. However, as the horizontal distance increases beyond the ice thickness (1 m in this case), the increase in the diagonal distance approaches the increase in the horizontal distance because the cosine of the angle between the diagonal and horizontal approaches one. Thus, much thinner snow, relative to the vertical or shorter diagonal pathways, would be necessary to make longer diagonal pathways preferable. Albeit highly simplified, this thought experiment provides intuition that sub-meter-scale variability in snow depth is likely to be the most important spatial scale of snow depth variability for horizontal heat flux enhancement. Indeed, this thought experiment was the inspiration for the more technical analysis herein.

Anecdotal experience with snow on Arctic sea indicates that there is considerable depth variability on sub-meter length scales, although there is limited published evidence. From semivariograms (Isaaks and Srivastava, 1989) of snow depth (Iacozza and Barber, 1999; Sturm and others, 2002a; Liston and others, 2018), we estimate that the combination of sub-meter-scale snow depth variability and measurement uncertainty has a standard deviation of 2.5–5 cm. Semivariogram analysis of these data cannot distinguish measurement uncertainty from spatial variability at this scale. Combining the potential for sub-meter-scale snow depth variability with the thought experiment above leads to our hypothesis that it may be important to include sub-meter-scale variability and horizontal heat fluxes when simulating the impacts of snow spatial variability on net conductive heat flux for level, Arctic sea ice.

In this work, we focus solely on the impacts of spatial variability of snow depth. Thus, we will consider only steady-state heat conduction with constant values of thermal conductivity in the snow and ice, and no sources or sinks of heat inside the snow-ice system. These simplifications are commonly made (e.g., Abraham and others, 2015; Popović and others, 2020). The impacts of non-conductive heat transport mechanisms such as vapor diffusion (Calonne and others, 2014), convection (Sturm and Johnson, 1991; Colbeck, 1997), and brine drainage (Niedrauer and Martin, 1979) are beyond the scope of this study but may merit further investigation. For brevity, “conductive heat flux” is shortened to “heat flux” throughout this manuscript. With these simplifications, the temperature distribution in the snow and ice is mathematically represented by the solution to Laplace's Equation (Eqn. (1); Laplace, 1822) and the heat flux at any location within the snow and ice is given by Fourier's Law (Eqn. (2); Fourier, 1822):

(1)$$\nabla \cdot \left(\kappa( x,\; y,\; z) \nabla T( x,\; y,\; z) \right) = 0$$

(2)$$\vec{q}( x,\; y,\; z) = -\kappa( x,\; y,\; z) \nabla T( x,\; y,\; z) $$

where T is the temperature; κ is the thermal conductivity; and $\vec {q}$ is the heat flux. If the heat flux were only vertical, then the steady state conductive heat flux (which we denote q _v) would simplify to (Eqn. (3): Semtner, 1976):

(3)$$q_v( x,\; y) = {T_o - T_a\over ( {h_s( x,\; y) }/{\kappa_s}) + ( {h_i( x,\; y) }/{\kappa_i}) }$$

where T _o and T _a are the temperatures at the ocean-ice and snow-air interfaces respectively; h _s and h _i are the snow depth and ice thickness respectively; and κ _s and κ _i are the snow and ice thermal conductivities respectively. Equation (3) assumes that the snow-air interface temperature is spatially uniform. In theory, spatial variability in conductive heat flux should produce spatial variability in surface temperature that would be balanced by variability in outgoing longwave radiation. However, spatially-resolved measurements of skin temperatures on level ice on the MOSAiC Expedition from a helicopter-borne thermal infrared camera (Thielke and others, 2021) showed minimal spatial variability in skin temperature due to snow depth variability. In this work we will assume the snow-air interface temperature is spatially uniform. Note that if the snow and ice thicknesses were uniform, then all net heat flux would be vertical and uniform.

Vertical-only heat flux simulations are a convenient way to estimate the additional heat flux due to snow spatial variability for two reasons. First, it can be calculated from a collection of independent point measurements of snow and ice thickness (e.g., a transect or stake array). For logistical reasons, most measurements of snow depth and ice thickness come from transects or stake arrays with (sometimes non-uniform) measurements spacing of 1–10 meters (e.g., Hanson, 1965, 1980; Radionov and others, 1997; Sturm and others, 2002a; Perovich and others, 2003; Sturm and others, 2006; Iacozza and Barber, 2010; Petrich and others, 2012; Webster and others, 2015; Liston and others, 2018; Rösel and others, 2018; Moon and others, 2019). However, including the impacts of horizontal heat fluxes in net heat flux estimates would require measuring snow depths and ice thicknesses in a manner that retains their spatial relationships on a dense sampling grid. In fact, we are unaware of any measurement campaigns that have done so at sub-meter horizontal spacing and it may not be possible without disturbing the surface. Second, vertical-only heat flux simulation (Eqn. (3)) can be computed with a single line of Python (Rossum and Drake, 2010) code and trivial computational expense. Estimating three-dimensional heat fluxes (Eqns. (1) and (2)) is much more complicated and computationally expensive.

Finally, two prior studies have explicitly considered horizontal heat flux in Arctic sea ice. Sturm and others (2002b) simulated two-dimensional heat flux in a 40 m long transect of multi-year ice containing ice hummocks and refrozen melt ponds. They found that the spatially variable snow and ice cover resulted in a heat flux 40% greater than uniform snow and ice covers. However, they did not distinguish the impacts of snow spatial variability from ice spatial variability. Popović and others (2020) developed a mathematical model of the snow topography on level, first-year ice based upon the random placement of Gaussian mounds with a fixed aspect ratio. Based on this model, they derived functions relating the mean, variance, and semivariogram range to estimated heat flux for vertical-only and including horizontal cases. For three cases of level, first-year ice near Utqiaġvik, AK, they estimated that simulating vertical-only heat flux for a spatially variable snow cover increased net heat flux by 4–11% compared with a uniform snow cover. However, they estimated that including horizontal heat flux did not result in an increased net heat flux relative to vertical-only heat flux (<1% increase). Because of the smoothness of the Gaussian mounds, it is unclear that their model accurately represents the small-scale surface roughness of snow on sea ice.

We measured the snow surface topography on level, landfast, first-year ice in Elson Lagoon near Utqiaġvik, AK at horizontal resolutions of 3–15 cm. We use this topography to simulate conductive heat flux through snow and ice and estimate the impacts of snow spatial variability on the net heat flux by simulating vertical-only and including horizontal fluxes. For simulated 1-m-thick ice, we find that simulating vertical-only heat flux increases the net heat flux in the largest region we studied (Fig. 1c and Table 1) by 6% compared to a uniform snow cover with the same mean depth. Including horizontal fluxes increases the net flux by 10% (compared to uniform snow) for this region. However, if we had increased our horizontal measurement spacing to 1 m or greater, then including horizontal fluxes has almost no impact compared to the vertical-only simulations (Fig. 2)–supporting our hypothesis that sub-meter-scale spatial variability is important for horizontal heat flux. We discuss the implications of these results and make recommendations for future investigations of heat flux on Arctic sea ice in the Discussion.

Figure 1. Study area in Elson Lagoon, AK on 17 April 2022: (a) photograph, (b) 3D rendering of TLS data (c) assumed snow depth map from TLS data. 3D rendering (b) is from approximately the same viewpoint as the photograph (a)–marked by camera glyph in (c). Camera lens distortion is not modeled. Lighting for 3D scene (b) is from the same elevation and azimuth as sun position in photograph (a) so that shadows generally correspond. The same colormap is used in (b) and (c). The nested regions used in the heat flux analysis are marked in (c). Photo credit for (a): Serina Wesen.

Table 1. Summary statistics and simulated heat fluxes for each region in Figure 1c.

Figure 2. Comparison of how increasing the horizontal measurement spacing (i.e., degrading the resolution) impacts simulated heat flux ratios relative to a uniform snow cover for simulations of vertical-only heat flux and including horizontal fluxes on each region in Figure 1c. Increasing measurement spacing does not impact the vertical-only heat flux ratios (dashed lines), but does reduce the horizontal heat flux ratios (solid).

2. Materials and methods

2.1 Data collection and processing

We used a Riegl VZ1000 Terrestrial Laser Scanner (TLS) to measure the snow surface topography on an untrodden 200x50 m area of level, landfast, first-year ice in Elson Lagoon, AK (71.349$^\circ$N, 156.526$^\circ$W) on 17 April 2022 as part of the Snow ALbedo eVOlution (SALVO) project. TLS, an application of Light Detection and Ranging (LiDaR), is a standard technique for measuring snow topography (Deems and others, 2013) wherein a tripod-mounted laser scanner maps the surface using a line-of-sight laser emitter and detector. To measure the entire area, we collected TLS data from five scan positions arrayed around the outside of the 200 × 50 m area. We co-located the scan positions with one another within an uncertainty of 3 mm using Riegl 10 cm reflectors mounted on posts frozen into the ice and a plane-matching algorithm in RiSCAN (Riegl's software for TLS acquisition and processing). The ranging precision of the VZ1000 is 5 mm and the beam divergence is 0.3 mrad.

We selected four nested regions in the center of the untrodden measurement area near the central scan position (where our measurement resolution was highest) to simulate heat flux (Fig. 1c and Table 1). For each region, we converted the TLS pointcloud data to gridded data by first constructing a triangulated surface from the pointcloud (Kazhdan and others, 2006) and then sampling this surface at the desired resolution. For the higher resolution regions, the sizes of the regions were limited by our desire to keep the percentage of grid cells without TLS measurement points in them at less than 5% (admittedly an arbitrary cutoff). The size of the largest region was motivated by a desire to exceed the typical structural length scale of 20 m (the distance beyond which the snow depths become uncorrelated; Sturm and others, 2002a; Itkin and others, 2023). Simulating 3D heat flux on domains larger than our largest region was computationally infeasible with our available resources.

2.2 Heat flux simulation

Since our interest is the impact of the snow variability, we will assume that the ice is uniformly 1 m thick with level ocean-ice and ice-snow interfaces (Popović and others, 2020, makes the same assumptions, so these choices facilitate comparing our results). Because ice in Elson Lagoon is first year ice that is protected from substantial ice dynamics, the assumptions of level interfaces is reasonable (e.g., Webster and others, 2014, note also the generally level character of the ice extending to the horizon in Figure 1). We set the same vertical position of the snow-ice interface for all of the regions (as is sensible since they are overlapping) such that all snow depths are greater than 1  cm–consistent with our observation that all ice was snow covered. With these assumptions, the coefficients of variation of snow depth (Table 1) in the larger regions fall within the range typically observed for snow on level, first-year ice.

For each region, we convert the snow and ice volume into a three-dimensional tetrahedral mesh using Gmsh (Geuzaine and Remacle, 2009). We assume constant ocean-ice and snow-air interface temperatures (Dirichlet boundary conditions) and that there is no horizontal heat flux at the lateral borders of the domains (Neumann boundary conditions). We set the snow thermal conductivity at 0.14 W  m⁻¹ K⁻¹ and ice thermal conductivity at 2.0 W  m⁻¹ K⁻¹ (same values used by Popović and others, 2020). For each region, we numerically solve Eqns. (1) and 2 via the Finite Element Method using second-order Lagrange elements (using FEniCS: Kirby, 2004; Kirby and Logg, 2006; Kirby, 2012; Logg and Wells, 2010; Logg and others, 2012a, 2012b, 2012c; Ølgaard and Wells, 2010; Alnaes and others, 2014; Alnæs and others, 2015; Scroggs and others, 2022) assuming a 20 K temperature difference between the ocean-ice and snow-air interfaces. Note that this is a simplified version of Eqns. (1) and (2) where κ(x, y, z) is either κ _s or κ _i in the snow and ice respectively. We numerically integrated the vertical heat flux at the ocean-ice interface to estimate the net mean heat flux. For each region, we also estimated the vertical-only heat flux (for the assumed spatially variable snow cover) and the heat flux assuming uniform snow thickness with Eqn. (3). Finally, to investigate the impacts of the measurement resolution on the results, we down-sampled the gridded snow topographies by factors of 2–20 and computed the heat fluxes.

3. Results

Including horizontal heat flux simulated more total heat flux than that simulated with vertical-only heat flux and the heat flux assuming uniform snow depth (Table 1 and Fig. 2). Increasing the horizontal measurement spacing (i.e., degrading the resolution) does not affect the ratio between the simulated vertical-only heat fluxes and that simulated with a uniform snow cover (the dashed lines are roughly constant in Fig. 2). However, increasing the measurement spacing reduces the impact of including horizontal fluxes in the simulation (the solid lines decrease to the right in Fig. 2). When the horizontal measurement spacing is 1 m or greater, we do not find that including horizontal heat fluxes resulted in a net heat flux greater than the vertical-only simulations. We also examined the sensitivity of our results to ice thickness by simulating region iii with an ice thickness of 2 m (not shown). The overall magnitude of the heat flux decreased, as expected for thicker ice, but the relative effect of including horizontal heat flux was similar. Additionally, we explored the sensitivity of our results to increasing the snow thermal conductivity to 0.3 W m⁻¹ K−1 (Supplemental Material: Table S1 and Fig. S1). This does not change the relative effect of including horizontal heat flux, but does decrease the impact of snow depth variability on heat flux. For example, for region iv at 15 cm resolution, including horizontal heat flux results in 6% greater heat flux than a uniform snow cover, whereas vertical-only heat flux is 3% greater than a uniform snow cover. For comparison, with snow thermal conductivity of 0.14 W m⁻¹ K−1, for region iv at the 15 cm horizontal resolution, including horizontal heat flux results in 10% greater heat flux than a uniform snow cover, whereas vertical-only heat flux is only 6% greater than a uniform snow cover.

4. Discussion and conclusion

The simulation results generally confirm the intuition from the thought experiment in the Motivation: sub-meter-scale snow depth variability increases the net heat flux when our simulations include horizontal fluxes. Furthermore, in this case, the additional net simulated heat fluxes when including horizontal fluxes relative to vertical-only simulations were a similar magnitude, albeit smaller, to the impacts of simulating vertical-only fluxes in a spatially variable snow cover compared to a uniform snow cover. This contradicts the results of Popović and others (2020), although we made the same assumptions about ice thickness and snow and ice thermal conductivities. Popović and others (2020) reported that simulating horizontal fluxes had no impact on the net heat flux beyond vertical-only simulation. We suspect that their model, which represents snow dunes as smooth Gaussian mounds, is not accurately representing sub-meter-scale snow variability. When we downsample our data to horizontal measurement spacings of 1 m or greater, including horizontal heat fluxes no longer simulates additional net heat flux compared to vertical-only simulations (i.e., the same result as Popović and others, 2020). The modeling approaches are also different, Popović and others (2020) derived analytical mathematical relationships between summary statistics (mean, variance, and semivariogram range) and heat flux whereas we conducted finite element simulations with a specific, measured snow topography.

The impacts of degrading the horizontal resolution were qualitatively similar for all four nested regions: degrading the resolution reduced the heat flux simulated when including horizontal fluxes but did not impact the vertical-only simulations. The smaller regions (i and ii) were included to assess whether any novel behavior occurred at extremely high (<10 cm) sampling resolution (e.g., a dramatic increase in the heat flux simulated when including horizontal fluxes), and because technical limitations prevented us from measuring such high resolution topography on larger domains. However, using such small domains introduces quantitative differences in simulated heat fluxes between the regions that are due to the domain size. Most importantly, the variance of snow depth in small regions (e.g., i and ii) tends to be less than that in larger regions (e.g., iii and iv), because snow depth is spatially auto-correlated at distances less than the semivariogram range (typically around 20 m). When snow depth variance is lower (small regions), the impacts of depth variability are smaller (i.e., uniformity is a better approximation). Consider the limiting behavior. If our entire domain were a column of snow and ice with infinitessimal horizontal extent, then all three types of simulated heat flux would be identical. For regions i and ii, sampling resolutions of less than 15 cm (the resolution of region iv) caused insignificant differences in the simulated heat flux (<0.05 W m⁻²). Thus we conclude that the sampling resolution of the larger regions (iii and iv) was sufficient to observe these effects. Since the larger regions are less affected by autocorrelation, we expect them to be more representative of large-scale effects although further work with more computational resources could help extend this analysis to larger domains.

Quantitatively assessing the representativeness of the sub-meter-scale snow roughness on this particular region of level, first-year Arctic sea ice is beyond the scope of this study (and likely requires more data collection of high resolution surface topography). However, the assumed mean, standard deviation, and coefficient of variation of snow depths for our two larger regions (iii and iv) are within the range of typical values on level, first-year Arctic sea ice (Iacozza and Barber, 1999; Sturm and others, 2002a; Petrich and others, 2012; Webster and others, 2015; Moon and others, 2019, the mean and standard deviation are at approximately the 25 percentile, the coefficient of variation is near the median).

Assessing the importance of horizontal heat flux enhancement on large spatial or temporal scales is challenging because there are almost no measurements of sub-meter-scale snow thickness variability. This includes this study, we could only measure snow topography and infer depth based on level-ice assumptions. Collecting TLS scans prior to the first snowfall (as is sometimes done in terrestrial environments: e.g., Hartzell and others, 2015) could address this data gap, although this is logistically challenging on sea ice. Another approach would be targeted studies to see if the fractal-scaling behavior observed in super-meter-scale snow depth (Deems and others, 2006; Trujillo and others, 2007, 2009; Moon and others, 2019) extends to sub-meter-scales on Arctic sea ice. If so, information about sub-meter-scale variability could be inferred from super-meter-scale measurements. Regarding models of the snow and ice cover, our results suggest that if they seek to simulate the impact of horizontal heat fluxes, they should represent sub-meter-scale snow depth variability. How exactly, to simulate or parameterize this variability will require more measurements and process studies of sub-meter-scale variability. Finally, although we consider a single field site and assumed uniform ocean-ice and ice-snow interfaces, the order of magnitude of the heat flux enhancements we simulate may be of use to other researchers. For our largest region (iv), the simulated heat flux is 10% greater including horizontal heat flux than that for uniform snow with the same mean snow thickness. Whereas, simulating vertical-only fluxes increased the net simulated heat flux by 6%.