Brief communication: Evaluation of multiple empirical, density-dependent snow conductivity relationships at East Antarctica

. Nine density-dependent empirical thermal conductivity relationships for firn were compared against data from 10 three Automatic Weather Stations at climatically-different East Antarctica sites (Dome A, Eagle and LGB69). The empirical relationships were validated using a vertical, one-dimensional thermal diffusion model and a phase-change based firn diffusivity estimation method. The best relationships for these East Antarctica sites were identified by comparing the modeled and observed firn temperature at the depth of 1 m and 3m, and from the mean heat conductivities over two depth intervals (1-3m and 3-10m). Among the nine relationships, that proposed by Calonne et al. (2011) appeared to have the best 15 performance. The density and temperature-dependent relationship given in Calonne et al. (2019) do not show clear superiority than other density-dependent relationships. This study provides useful reference for firn thermal conductivity parameterizations in land modelling or snow-air interaction studies on the Antarctica Ice Sheet.


Introduction
In the Earth's climate system, snow cover has two important physical properties, its high albedo and its low thermal 20 conductivity. Both modulate heat exchange between the atmosphere and the surface (Dutra et al., 2010). Heat transport in the near-surface snow layer plays a key role in controlling the upper thermal boundary condition of ice sheets (Ding et al., 2020).
Snow is a porous and inhomogeneous material with thermal conductivity that can be anisotropic and depends on the microstructure of snow: proportion of air and ice, grain shape, grain size, bonds size, etc (Riche and Schneebeli, 2013).
Direct measurements of snow heat conductivity can be made with a needle probe, heated plate and tomographic 3D images 25 (e.g. Sturm et al., 1997;Calonne et al., 2011), all of which require intensive work. Alternative approaches include Fourier analysis methods that can estimate thermal diffusivity and reconstruct snow thermal histories from temperature measurements (Oldroyd et al., 2013), considering that the bulk/apparent heat diffusivity can be more effectively described than the whole physical process of snow metamorphism, as assumed also by needle probe measurement studies (Calonne et al., 2011). Similarly, the spatially averaged thermal diffusivity can be estimated from the changes of amplitude and phase of 30 2 a temperature cycle with depth in the medium (Hurley and Wiltshire, 1993;Oldroyd et al., 2013). The numerical inverse method (optimal control theory) is another possible approach for recovering thermal diffusivity by a least-squared method (Sergienko et al., 2008) or a recursive optimization approach (Oldroyd et al., 2013).
These numerical methods, however, need a relatively large number of temperature measurements, which can be difficult for large scale model studies. Thus, a widely-accepted alternative is to use laboratory-determined empirical 35 relationships to approximate the snow diffusivity and/or conductivity as a function of some typical and easily measured snow parameters such as snow density (e.g., Yen, 1981;Sturm et al., 1997;Calonne et al., 2011).
Density-dependent thermal conductivity relationships are widely used in various model studies. For example, the empirical relationship developed by Jordan (1991) was adopted by the land model CLM, snow model SNTHERM and many land surface energy balance studies, e.g., Wang et al. (2017). Lecomte et al. (2013) used the relationships in Yen (1981) and 40 Sturm (1997) for large scale sea ice-ocean coupling models. Applying the density dependent relationship in Calonne et al. (2011), Hills et al. (2018 investigated the heat transfer characteristics in the Greenland ablation zone. Steger et al. (2017) analyzed the melt water retention in the Greenland ice sheet by adopting the snow density-conductivity relationship given in Anderson (1976). Charalampidis (2016) used the relationship in Sturm (1997) to trace the retained meltwater in the accumulation area of the southern Greenland Ice Sheet. However, none of those relationships have been carefully validated 45 by in-situ data in Antarctica ice sheet.
In this paper, firn temperature data and snow density profile from three sites in East Antarctica were chosen to validate the applicability of these density/conductivity relationships (Table S1). We first describe the meteorological observations. After introducing the method for validating the empirical density-conductivity relationships, we then present the validation results, followed by discussions and conclusions. 50

Site and observational description
Several solar-powered automatic weather stations (AWS) have been deployed along Zhongshan to Dome A traverse route within the cooperative framework between Chinese and Australian Antarctic programs. These include deployments at LGB69 (in January 2002), EAGLE and Dome A (in January 2005). For more than 10 years since then, near-hourly meteorological measurements have been made of air and firn temperature (at several heights/depths), relative humidity, wind, 55 and air pressure. The data from the AWSs is remotely collected and relayed by the ARGOS satellite transmission system. Firn temperatures are measured (using FS23D thermistors in a ratiometric circuit with a resolution of 0.02K) at four depths below the surface. These were 0.1 m, 1 m, 3 m and 10 m when deployed, but they have slowly deepened with time due to snow accumulation. Due to heavy snowfall at LGB69, these data are only available for 2002-2008.
All three sites locate at western side of the Lambert Glacier Basin. LGB69 (70°50'S, 77°04'E; 1854 m a.s.l.) is only 192 60 km away from coast (Figure 1), and has an annual precipitation of 20 cm w.e. per year (~50 cm snowfall), strong wind (~8.5 m/s annual) and a mean annual air temperature of ~-26.10°C. The snow density increases from ~400 to 500 kg m -3 from surface to 10 m depth ( Figure S1). EAGLE (76°25'S; 77°01'E) is a typical "surface glazed" area with hard snow crust 3 because of the effect of drift snow. Its snow accumulation is 10 cm w.e. per year (30 cm snowfall), the snow density increases from ~380 to 550 kg m -3 from surface to 10 m ( Figure S2) depth and the mean annual air temperature is ~-40.80 °C. 65 Dome Argus (80°22′S, 70°22′E; 4093 m a.s.l.) is the highest point of the east Antarctic Ice Sheet. It is also the summit of the ice divide of the Lambert Glacier drainage basin, ~1248 km from the nearest coast, and the surrounding region has a surface slope of only 0.01 % or less. Dome Argus has extremely low surface air temperature (-52.1°C annual mean), specific humidity and snow accumulation rate (around 2 cm w.e. per year) and experiences no surface melt, even at the peak of summer (Ding et al., 2015). The surface snow is very soft here, ranges from ~270 to 450 kg m -3 in the top 10 m layer ( Figure  70 S3). There were no radiation measurements at the site. Therefore, it is nearly impossible to build a complete energy balance model at the snow surface of Dome Argus.

1) Numerical model method.
We validate the heat conductivity by a one-dimensional transient heat diffusion model, where is the firn temperature, K is the heat conductivity, is the heat capacity of snow, ρ s is the density of snow, and z is the depth below the snow surface. The vertical firn density profiles for three sites are shown in the Figures S1-S3. The heat capacity of snow is estimated by assuming snow is a mixture of air and ice,  85 where T d = abs(T model − T obs ), T dm is the mean value of T d , and N is the number of the temperature dataset.
2) Temperature phase change method. In this approach, we approximate the annual temperature cycles as sinusoidal functions (Demetrescu et al., 2007).
where is the firn temperature expressed as a function of depth and time , is the mean annual value of , is the 90 amplitude of the annual firn temperature cycle, is the frequency of the temperature cycle and is the wave phase of the annual cycle.
While it is common to fit a harmonic series (Fourier analysis) rather than a single sine wave to temperature variations, we found that this gave no advantage over Equation 4 for the data at the three sites because the temperature there have non-4 periodic temperature excursions during the "coreless" Antarctic winter. Assuming the snow is horizontally isotropic we can 95 estimate the apparent thermal diffusion, , from the changes of phase at different depths, = 2 ( 2 − 1 ) 2 ( 1 − 2 ) 2

(5)
The conductivity can then be recovered from and the heat capacity of firn.

Results and discussions
In Figures S4, S5 and S6, we show the comparisons of observed and modeled firn temperature using 9 different density-100 conductivity relationships at the Dome A, Eagle and LGB69 AWS station. The deviations of their differences are given in Table 1. At Dome A, the Jor and Sch relationships give us a significant discrepancy between observed and modeled firn temperature (Fig 2). The modeled firn temperature calculated by the Ca1, Lan and Van relationships show a closer agreement with the observed firn temperature. The density and temperature-dependent relationship, Ca2, however, does not appear to have a better performance than its density-dependent version Ca1 at Dome A. 105 In Figure 2 and Table 1, we can see that at Dome A, the Lan relationship gives the best performance at the depth of 1 m and 3 m, followed by the Van and Ca1 relationship. The Lan relationship was derived from in-situ snow conductivity measurements on Filchner Ice Shelf (Lange, 1985). It is the only relationship in this study that is based on in-situ Antarctica firn sample measurements. The Ca1 relationship was derived by analyzing a wide range of different snow samples from a number of different geographical locations with many different snow types (Calonne et al., 2011). The Van relationship is 110 old but is adopted in Cuffey and Paterson (2010), and still shows a nice performance in our model results.
Similarly to the Dome A case, at the Eagle station, the Ca1 relationship out-performs other relationships, followed by the Ca2 and Yen relationship. The Jor, Stu and Lan relationship, however, appear to be not suitable for parameterizing the firn conductivity, compared with other relationships at the Eagle station. At LGB69, however, the Sch and Jor relationship appear to be superior to other relationships, different from the cases at Dome A. The Jor relationship is based on the 115 experimental measurements in Yen (1962). Sch is also an experimental relationship, based on the data given in Mellor (1977). In this case, the Ca2 relationship also gives a smaller temperature difference, compared to the Ca1 result. Note that the same relationship may have a different performance at different depth ranges. For example, for the Lan relationship, the modeled and observed firn temperature shows very good agreement at the depth of 3 m, but has a relatively large discrepancy at the depth of 10 m. 120 We also estimate the spatially averaged annual mean thermal conductivity from the temperature phase shifts between the depth ranges of 0.1 -1 m, 1 -3 m and 3 m -10 m at the Dome A (Fig. S4), Eagle (Fig. S5) and LGB 69 (Fig. S6) Station, and compare them with the mean values corresponding to different density-conductivity relationships (Table S2).
The phase at different levels ( ) (also the phase shift, ) are determined from the least-squared fit to Equation (1), and the thermal diffusivity (conductivity) is then calculated by Equation (5).

5
Clearly, we can see in Table S2 that there is an increasing trend of conductivity with depth. Similar to the case in Table   1, the Ca1 and Ca2 relationship give close values than the conductivity values recovered from the phase change method at 3 different depth intervals, in consistent with the comparison results in performances. This is possibly a result of our assumption that the vertical density profile is kept constant in time, and that the heat capacity of firn is a linear relationship of the capacity of air and ice (Eqn 2). In addition, only one relationship is a funciton of firn temperature (Ca2). As temperature is an important parameter in affecting firn heat conductivity (Calonne et al., 2019), considering only density in other eight relationships may introduce model uncertainties in our evaluation. But since we consider a long time span of observation and model years (7 years for Dome A, 8 years for Eagle and 3 years for 135 LGB69), the overall deviation of modeled and observed temperature should be accountable in quantifying the performance of different density-dependent conductivity relationships.

Conclusions
In this study, we apply two methods to validate nine different density-conductivity relationships: 1) by applying a 1D vertical heat diffusivity model, we compare the modeled firn temperature at the depth of 1m and 3m with observations; 2) 140 we compare the mean empirical snow conductivity at three depth intervals (0.1-1m, 1-3m and 3-10m) according to the phase change derived temperature variations.
It is found that some empirical density relationships have generally good model performance and agree well with phase change recovered conductivity, but they show diverse behaviors at different depth levels. Based on these two methods, we find that the relationship proposed by Calonne et al. (2011) (Ca1) generally has an overall best performance. The Jordan 145 (1991) relationship (used in snow models like CLM and SNTHERM), however, does not present a very good model results for Dome A, Eagle and LGB69 Stations. All in all, no conductivity-density relationship is optimal at all sites and the performance of each varies with depth.
The 3 AWS sites in in the paper cover a large range of elevation and distance from coast. We thus argue that our findings can shed some lights on firn thermal studies (e.g., the applicability of different firn density-conductivity 150 relationships) in Antarctica.

Supplement
The supplement related to this article is available online at: ****** 155 LGB69 (e, f).