Implementation of dynamic snow density within GlobSnow snow water equivalent retrieval methodology

Snow water equivalent (SWE) is an important variable in describing global seasonal snow cover. Traditionally, SWE has been measured manually at snow transects or using observations from weather stations. However, these measurements have a poor spatial coverage, and a good alternative to in-situ measurements is to use spaceborne passive microwave observations, which can provide global coverage at daily timescales. The reliability and accuracy of SWE estimates 10 made using spaceborne microwave radiometer data can be improved by assimilating radiometer observations with weather station snow depth observations as done in the GlobSnow SWE retrieval methodology. However, one possible source of uncertainty in the GlobSnow SWE retrieval approach is the constant snow density used in modelling emission of snow. In this paper, three versions of spatially and temporally varying snow density fields were implemented using snow transect data from Eurasia and Canada and automated snow observations from USA. Snow density fields were used to post-process the baseline 15 GlobSnow v.3.0 SWE product. Decadal snow density information, i.e. fields where snow density for each day of the year was taken as the mean calculated for the corresponding day over ten years, was found to produces best results. Overall, postprocessing GlobSnow SWE retrieval with dynamic snow density information improved overestimation of small SWE values and underestimation of large SWE values, though underestimation of SWE values larger than 175 mm was still significant.


Introduction 20
Snow water equivalent (SWE) is an important property of the seasonal snow cover and estimates of SWE are required in many hydrological and climatological applications, including climate model evaluation (Mudryk et al., 2018) and forecasting freshwater availability. Maximum SWE before the start of spring snow melt is one of the most important snow characteristics for run-off and river discharge forecasts (Barnett et al., 2005;Barry, 2002). 25 Snow depth or SWE can be estimated by interpolating surface snow depth (Dyer and Mote, 2006) or snowfall measurements (Broxton et al., 2016). However, the limited spatial and temporal coverages of the ground-based measurements, especially in northern and alpine regions, limit the quality of estimates (Broxton et al., 2016;Mortimer et al., 2020). An alternative approach for estimating SWE is to use satellite measurements as they can provide global spatial coverage and good temporal resolution.
Passive microwave observations are commonly used as these provide frequent repeat coverage and the influence of atmospheric conditions on the observation is limited. Furthermore, passive microwave radiometer data is also available from 1978 onwards, which allows the analysis of long time series. Many passive microwave radiometer-based approaches for 35 estimating SWE are adopted from an algorithm proposed by Chang and Foster (1987) for estimating snow depth from horizontally polarized Scanning Multichannel Microwave Radiometer (SMMR) measurements. The algorithm is based on the difference in measured brightness temperatures at a frequency insensitive to dry snow, around 19 GHz, and at a frequency sensitive to dry snow, around 37 GHz. The uncertainty of SWE retrievals based on the radiometer measurements alone can be quite high (Mudryk et al., 2015). Retrieval algorithms that use only radiometer data tend to underestimate SWE on deep snow 40 conditions (Derksen et al., 2005) and the performance of these algorithms is even more limited in wet snow conditions (Armstrong and Brodzik, 2001).
To overcome the problems connected to stand-alone passive microwave SWE retrievals, ground-based observation, and satellite radiometer data can be combined as done in the assimilation approach for SWE retrieval introduced by Pulliainen 45 (2006) and complemented by Takala et al., (2011). This assimilation-based approach was used as the baseline method for the Global Snow Monitoring for Climate Research (GlobSnow) initiative of the European Space Agency (ESA). The GlobSnow method has been shown to produce good results when compared to typical stand-alone radiometer algorithms (Mortimer et al., 2020). The GlobSnow version 3.0 (GSv3.0) climate data record with spatial bias correction was used for accurate reconstruction of the northern hemisphere snow mass and its trends for period of the 1979-2018 (Pulliainen et al., 2020). 50 Improving the GlobSnow SWE retrieval methodology will help to further enhance our understanding of the northern hemisphere snow conditions and its changes.
The GlobSnow SWE retrieval utilizes a fixed density of 240 kg m -3 throughout the retrieval regardless of snow depth, location, or time of the year (Takala et al., 2011), which is a known source of uncertainty in the retrieval. The density of snow changes 55 with time and place, and it is greatly affected by surrounding weather conditions. For example, wind breaks down snow crystals, both on the ground and falling from the sky, which allows snow crystals to pack together tightly and increases the density of snow (Jordan et al., 1999). The age of the snow cover also affects its density as the snow slowly packs tighter due to gravity (Maurice and Harold, 1981).

60
One approach considered for GlobSnow SWE retrieval methodology for estimating temporally and spatially varying snow densities was to use a statistical snow density model presented by (Sturm et al., 2010) which predicts density of snow as a function of the snow depth, day-of-the-year, and snow class (Luojus et al., 2013). However, applying densities obtained using https://doi.org/10.5194/tc-2021-15 Preprint. Discussion started: 22 January 2021 c Author(s) 2021. CC BY 4.0 License. this approach did not improve retrieval skill notably. A different approach for obtaining varying snow density information is to use available snow density data to create snow density fields by applying temporal and spatial interpolation. 65 In this study, dynamic snow density information obtained from ground measurements using interpolation is used to postprocess the GSv3.0 SWE climate data record. Three different versions of the snow densities are implemented for Eurasia for the years 2000 to 2009. Additionally, one version of the dynamic snow densities is also implemented for the whole northern hemisphere for the whole period of GSv3.0 SWE data record, 1979-2018 2 Data and methods

Eurasia
In-situ snow density and SWE measurements were used to obtain dynamic snow density fields and validate the results of SWE retrievals. SWE and density datasets for Eurasia were obtained from Russia (Bulygina et al. 2011) and Finland (Haberkorn,75 2019). These datasets contain snow transect data. Snow transects consist of manual gravimetric snow measurements made at multiple locations along a pre-defined transect several hundreds of meters to several kilometres in length which are averaged together to obtain a single representative (in regard to the spatial resolution of utilized passive microwave radiometers) SWE value for a given transect on a given date.

80
Russia data, from a substantial network of snow transects, has been made available via the All-Russia Research Institute of Hydrometeorological Information-World Data Center (RIHMI-WDC) website. This Russian snow survey dataset contains data from routine snow surveys operated at 515 meteorological station locations and data is available from 1966 to 2019 (Bulygina et al. 2011); data from 1979 to 2018 is used in this study. Routine snow surveys are run through the cold season every ten days or every five days during the intense snowmelt season. The Finnish Environment Institute (SYKE) has a network 85 of about 160 snow survey courses that have been operated from the beginning of the 20th century (Haberkorn, 2019). These 2 to 4 km long snow survey courses that go through different landscapes are visited monthly, and 80 snow depth measurements are made along the snow course through varying landscape about 50 meters apart, eight snow density and SWE measurements are made along each snow course. An aggregate of the SWE measurements is applied to describe the SWE conditions for the snow course for the given sampling date. 90

Pre-processing of Eurasian data
The Russian snow transect dataset was divided into two parts. The division of data was done by finding the nearest neighbours and separating them into different datasets. The first part of the data was used for implementing the dynamic snow density https://doi.org/10.5194/tc-2021-15 Preprint. Discussion started: 22 January 2021 c Author(s) 2021. CC BY 4.0 License. maps and the second part, together with snow transect data from Finland, was used for validating snow density and SWE results. The implementation dataset contains 257 locations and the validation data is formed from data from 625 locations. 95 Implementation and validation snow transect locations are shown in Fig. 1. We divided the Eurasian snow density implementation dataset into three distinct ways to create three different versions of 100 dynamic snow density maps. For the first version, called the multi-decadal version, all data between 1979 and 2018 was used.
The second version, called the decadal version, uses data from 1999 to 2009. The third version, called the annual version, uses data from 2004 to 2008 and produces daily density maps for each year using only data from the year under investigation.
A day-of-the-winter (DOW) value was added to each snow density measurement. DOW values are a modification of the day-105 of-the-year (DOY) values. DOW values start from September 1 and then continue to grow from there until the last day of August. This means that January 1 has a DOW value of 123, and August 30 has a value of 365 or 366. They are used because in the northern hemisphere winter season spans from September until June over the new year.
The Eurasian snow survey data was filtered to remove all negative density observations and all observations larger than 1000 110 kg m -3 . These measurements are most likely erroneous as snow densities typically range between 50 and 550 kg m -3 (Fierz et al., 2009). After filtering, average snow density values were calculated for each DOW that had at least one measurement. Most density measurements have been done systematically on the same DOW from year to year, with few exceptions. This difference in measurement days may cause average densities to fluctuate from one day to another as some density values are 115 not averages but measurements from one specific year. To avoid these fluctuations in multi-decadal and decadal versions, if two consecutive DOWs had measurements, the average density is calculated using data from both days, and this density was assigned the DOW value of the first day. Outlier data points were removed from multi-decadal and decadal datasets after averaged densities were calculated. Outliers were determined to be data points that differ from two previous and two following points by more than 50 kg m -3 . Figure 3 shows how average densities calculated for each DOW and non-consecutive DOWs 120 differ for one snow transect location for the multi-decadal dataset. The depicted snow transect is located in western Russia.   western United States (Serreze et al., 1999). SNOTEL dataset differs from the Canadian and Eurasian datasets, as it consists of automated daily measurements instead of manual snow transect measurements. SNOTEL stations collect data on snowpack SWE, snow depth, precipitation, and air temperature. SWE is measured by a snow pillow filled with an antifreeze solution. 135 Hourly data are available from the snow pillows, but daily measurements were used as they are more robust as hourly data is easily affected by wind and sensor issues.

Baseline SWE retrieval
The European Space Agency (ESA) GlobSnow and succeeding projects have produced a family of daily satellite-based SWE 145 climate data records spanning over 40 years. The most recent GSv3.0 data record is based on methodology introduced in (Pulliainen, 2006;Takala et al., 2011) and the latest version is presented in detail in (Luojus et al., 2021). The retrieval algorithm combines satellite-based passive microwave measurements with ground-based synoptic weather station snow depth observations by Bayesian non-linear iterative assimilation. The GlobSnow approach uses two vertically polarized brightness temperature observations at 19 and 37 GHz and a scene brightness temperature model (the HUT snow emission model (Pulliainen et al., 1999)). First, an effective snow grain size is estimated for grid cells that coincide with weather station snow depth observations. The snow grain sizes are used to construct a kriging interpolated background map of the effective grain size, including an estimate of the effective grain size error. This spatially continuous map of grain size is then used as an input for HUT model inversion to provide an estimate of SWE. The 155 daily weather station snow depth measurements are also used to form a continuous background field of snow depth independently from passive microwave measurements. The interpolated snow depth field is fused with space-borne brightness temperature observations, using the scene brightness temperature model in a Bayesian approach that weights all information sources with their estimated variances, to provide the final SWE estimates. A constant value of snow density is used (240 kg m -3 ). 160 The retrieval method does not produce SWE estimates for mountainous areas, glaciers, or Greenland. The data record is based on data from the SMMR aboard NIMBUS-7, and SSM/I and SSMIS sensors onboard DMSP 5D F-series satellites, and synoptic weather station snow depth data from the northern hemisphere.

Creation of snow density fields 165
Two main steps for generating snow density fields are temporal and spatial interpolation. Snow density measurements are made usually every 10 or 15 days and thus, there are many days without snow density observations. Simple linear interpolation was used to obtain estimates of snow density values for the days lacking observations. The length of the snow season depends on the year and place but to get similar series for each station, and for each version, 170 the average of three first existing densities was added to DOW 30 if a station did not have measurements from DOW 30 or before this day. Similarly, if the station did not have any measurement after DOW 280, the average of the last three densities was calculated and added as the density for DOW 280. After this procedure, all stations have density values from DOW 30 to DOW 280, and interpolation could be performed for this period. Figure 4 shows the results of interpolation for all three versions for one snow transect location. The interpolation of the multi-decadal dataset produces smoothest results and yearly data shows 175 the most fluctuation.

180
We used ordinary kriging interpolation for interpolating snow density values for areas without observations. Kriging interpolation is a spatial interpolation method that predicts values for the location with no measurements based on the spatial autocorrelation of measured values (Goovaerts, 1997). This means that two closely located points are more likely to have similar values than two points further afield. The main advantages of kriging interpolation compared to nongeostatistical interpolation methods are that kriging interpolation provides variance of predicted values and that the spatial smoothing is 185 defined through variogram. The model of spatial variability can be expressed as (Høst, 1999): where ( ) denotes the predicted value at some location, ( ) is the deterministic function describing the trend component of ( ), and ( ) is the stochastic locally varying but spatially dependent residuals.

190
The spatial autocorrelation is modelled with a semivariogram (Goovaerts, 1997): with assumption of intrinsic stationarity (variance on the right-hand side is only dependent on the vector difference d) and assumption that the process is isotropic (the autocorrelation is only dependent on the distance between observations). The empirical semivariogram can be estimated from the observations as follows (O'Sullivan and Unwin, 2010): 195 where ( ) and ( + ) are sampled data pairs at distance and ̃ is the distance between observations. In this study, exponential function is used for the variogram: ( ) = 1 * exp( * 2 ) + 0 . 200 Snow density values were predicted for each pixel in 25 km Equal Area Scalable Grid (EASE-Grid version1, to match GSv3 processor grid) between latitudes from 42° N to 80° N, and for longitudes from 20° E to 180° E. The snow density values were estimated using coordinates for centres of each pixel. Density maps were made for all three versions of snow densities for each day starting from DOW 30 and ending at DOW 280. Pixels that are not on land are assigned a value of -1, and land areas outside the area of interpolation are designated a constant density of 240 kg m -3 (effectively retaining the original retrieval 205 methodology for those regions where snow density could not be reliably established).

Usage of snow density information and validation
The derived snow density information is used to post-process baseline GSv3.0 SWE retrieval. Post processing of SWE retrieval means that SWE values obtained are scaled with the ratio of dynamic and constant snow density: where has values of 240 kg m -3 . Scaling is performed for each pixel within the area for which dynamic densities are available, for the regions outside dynamic snow density information, the constant density consideration is retained.
The obtained snow densities and post-processed SWE datasets were validated using independent validation data. Validation locations were separated from the data used for generating snow density fields to ensure independent cross-validation. Root-215 mean-squared (RMS) error, bias, correlation coefficients and mean absolute error (MAE) are the four statistical measures used for assessing the performance.  Figure 5 shows the comparison of estimated and observed densities for multi-decadal, decadal, and annual snow densities.
Multi-decadal and decal version of snow densities exhibit similar behaviour with the multi-decadal version having slightly smaller RMSE and larger correlation coefficient but larger bias than the decadal version. MAE is equal for these two versions. 230 The annual version differs from the two other versions more and it estimates snow densities below 200 kg m -3 and above 300 kg m -3 better than the other versions. The multi-decadal and decadal versions of snow densities are produced from data that are averaged from a large number of measurements. This averaging means that the highest and lowest measurements have a diminished effect on estimated snow densities. Annual densities have a wider range of densities present in the data used for deriving these dynamic densities compared to the range of densities used for deriving the other two sets of the density maps.

Post-processed SWE
The three different density sets were used to post-process the baseline GSv3.0 SWE data between 2000 and 2009. Obtained SWE values were compared to validation SWE measurements of the Eurasian dataset. Validation was performed for SWE values up to 500 mm and SWE values up to 150 mm, as the bulk of the observations are below this value. Table 2 summarizes the results of the SWE validation (for different snow density realizations). Figure 6 shows mean errors for the baseline, multi-245 decadal, decadal and annual versions.
Post-processing with the multi-decadal and the decadal densities produces similar results and post-processing with these two versions improves the baseline product. These two versions have smaller RMSE, MAE, and bias than the baseline retrieval.   considerably. Estimates of smaller SWE values are similar to estimates of the other two post-processed versions. This can be caused by the annual density dataset having a larger range of densities than the other two density datasets. If SWE estimation has been close to correct, but the density used in the retrieval is far from the estimated density, post-processing causes SWE estimation to change significantly. A wider range of densities causes more significant changes in the post-processing.

Northern hemisphere, 1979-2018 270
Based on the results obtained for the years 2000-2009, the decadal version of snow densities was extended to cover the whole northern hemisphere and the whole period of the baseline retrieval. Four separate sets of density maps were made, each covering one decade starting from the 1980s and ending in the 2010s. The density maps were made using the methods explained in chapter 2, with a few exceptions : 1) the spatial interpolation was performed for latitudes from 35° N to 80° N, and for longitudes from -180° W to 180° E to cover the same area as the baseline retrieval and 2) the North American dataset was 275 filtered before adding DOW values. The filtering was done in a way that locations that were within the same EASE-grid cell were combined and the average value of measurements was calculated for each day. The behaviour of decadal densities over a longer period is very similar to their behaviour over a shorter period. The RMSE for decadal densities in Eurasia between 1979-2019 was 50.1 kg m -3 and the correlation coefficient 0.694.

280
These new decadal snow density maps were used to post-process the whole GSv3.0 baseline dataset and these results were compared to validation datasets from Eurasia and North America. Again, density values up to 500 mm and 150 mm were considered. Table 3 summarizes the results of this evaluation. Table 3 shows results for the whole northern hemisphere and separately for Eurasia and North America. The RMSE was 43.34 mm and 42.49 mm for the baseline and post-processing with decadal densities, respectively. Figure 8 shows scatter plots for the whole northern hemisphere for the baseline and the post-processed datasets. As expected, post-290 processing improves overestimation of SWE values between 10 and 100 mm. When focusing on shallow to moderate conditions, SWE below 150 mm, the accuracy is considerably better with RMSE for the baseline being 31.87 mm and for the post-processed dataset, the error is 30.87 mm. This underlines the challenge for the satellite-based GlobSnow approach to estimate deep snow conditions as noted in (Mortimer et al., 2020).

295
Figures 9 and 10 show density scatter plots for the baseline and the post-processed datasets for Eurasia and North America datasets, respectively. The post-processed Eurasian dataset shows similar improvements as the overall northern hemisphere dataset, overestimation of SWE values between 10 and 100 mm has been mitigated. North American datasets shows also this behaviour but not as clearly as the Eurasian dataset. Figure 7 shows SWE maps for February 6, 2011, for GSv3.0 retrieval and post-processed retrieval. Maps show how values of SWE are clearly lower (darker colours) in the post-processed map in 300 Eurasia. In North America, this difference is not as clear. North American reference measurements consist of both snow transect and point measurements which may affect the results. SNOTEL data can also be unrepresentative of areas surrounding the stations as these stations are often located at a higher elevation and in areas that accumulate deeper snowpacks than the surrounding areas (Molotch and Bales, 2006).     snow density. Similarly, the post-processing helps to improve the overestimation of small SWE values. The constant snow density used in the retrieval tends to be too large in the early winter when the snow is fresh, and the density is at its lowest. 335 Post-processing with dynamic densities can be used to improve existing datasets and the post-processing procedure is straightforward to perform. However, implementing dynamic snow densities directly into the retrieval could possibly produce even more notable improvements and is a key research topic to be studied in the future. Snow density is one of the input parameters of the HUT snow emission model, determining the absorption coefficient in snow, refraction, and transmissivity 340 at the air-snow interface and transmissivity at the snow-ground interface through modelled permittivity of the snow layer. The HUT model is used for determining effective snow grain size over weather station locations as well as for obtaining the final SWE estimates using numeric model inversion. The final snow grain size (and its variance) at each location is the average grain size of six nearest stations. If the true snow density between stations changes significantly the variance of the estimated snow grain sizes increases. This in turn potentially reduces the weight of radiometer measurements on the final SWE 345 estimation, as well as the accuracy of the individual grain size estimates. Similarly, applying a wrong density value in the final retrieval step potentially deteriorates the accuracy of the HUT model and thus retrieval skill.
Linear interpolation was used for temporal interpolation to obtain density measurements for days without measurements.
However, snow density measurements contain some errors and these errors affect the interpolation. Thus, alternative methods 350 for determining behaviour of snow density throughout the snow season could be evaluated in future investigations.

Conclusions
In this study spatially and temporally changing snow density fields were implemented using snow density measurements from Eurasia and North America. The dynamic snow density fields were used to post-process GlobSnow version 3.0 SWE retrievals.
Post-processing was found to improve overestimation of small SWE values between 10-100 mm and underestimation of large 355 SWE values. This indicates that the constant density used in the baseline retrieval is too large for early winter and too small for late winter. The overall results indicate a clear path forward to improve the overall GlobSnow SWE retrieval methodology by application of dynamic snow density in the post-processing scheme.