Estimating Ground Snow Load Based on Ground Snow Depth and Climatological Elements for Snow Hazard Assessment in Northeastern China

Extreme snow loads can collapse roofs. This load is calculated based on the ground snow load (that is, the snow water equivalent on the ground). However, snow water equivalent (SWE) measurements are unavailable for most sites, while the ground snow depth is frequently measured and recorded. A new simple practical algorithm was proposed in this study to evaluate the SWE by utilizing ground snow depth, precipitation data, wind speed, and air temperature. For the evaluation, the precipitation was classified as snowfall or rainfall according to the air temperature, the snowfall or rainfall was then corrected for measurement error that is mainly caused by wind-induced undercatch, and the effect of snow water loss was considered. The developed algorithm was applied and validated using data from 57 meteorological stations located in the northeastern region of China. The annual maximum SWE obtained based on the proposed algorithm was compared with that obtained from the actual SWE measurements. The return period values of the annual maximum ground snow load were estimated and compared to those obtained according to the procedure suggested by the Chinese structural design code. The comparison indicated that the use of the proposed algorithm leads to a good estimated SWE or ground snow load. Its use allowed the estimation of the ground snow load for sites without SWE measurement and facilitated snow hazard mapping.


Introduction
Roof snow load for codified structural design is based on the basic ground snow load and a ground-to-roof factor. The basic ground snow load in the Chinese load code (MOHURD 2012) is assigned based on the T-year return period value of the annual maximum ground snow load, where T is equal to 50 years. The basic ground snow load could be estimated based on the historical records of the measured ground snow load, or snow water equivalent (SWE). However, this kind of historical record is available for only a few meteorological stations, while ground snow depth records are more widely available for many meteorological stations. This is because the resources required to measure the SWE are much greater than those needed to measure snow depths (Sturm et al. 2010). Consequently, it is valuable to infer or estimate ground snow load from ground snow depth. This could augment the SWE database and result in a better spatial representation of ground snow load, which can be used to estimate the nominal ground snow load in the context of structural design codes (Hong and Ye 2014;Mo et al. 2016).
There have been several efforts to estimate the ground snow load from the ground snow depth, including Fridley et al. (1994), Jonas et al. (2009), McCreight and Small (2014), Bruland et al. (2015), and Hill et al. (2019), among others. For the estimation, the snow depth and snow load (or Abstract Extreme snow loads can collapse roofs. This load is calculated based on the ground snow load (that is, the snow water equivalent on the ground). However, snow water equivalent (SWE) measurements are unavailable for most sites, while the ground snow depth is frequently measured and recorded. A new simple practical algorithm was proposed in this study to evaluate the SWE by utilizing ground snow depth, precipitation data, wind speed, and air temperature. For the evaluation, the precipitation was classified as snowfall or rainfall according to the air temperature, the snowfall or rainfall was then corrected for measurement error that is mainly caused by wind-induced undercatch, and the effect of snow water loss was considered. The developed algorithm was applied and validated using data from 57 meteorological stations located in the northeastern region of China. The annual maximum SWE obtained based on the proposed algorithm was compared with that obtained from the actual SWE measurements. The return period values of the annual maximum ground snow load were estimated and compared to those obtained according to the procedure suggested by the Chinese structural design code. The comparison indicated that the use of the proposed algorithm leads to www.ijdrs.com www.springer.com/13753 1 3 SWE) are considered to be directly related by the snowpack bulk density. The snowpack bulk density is a function of snow depth, air temperature, air humidity, snow accumulation history, and other climatological variables (Ellingwood and Redfield 1983;Sanpaolesi 1998;Sturm et al. 2010). Consequently, it could be feasible to estimate ground snow loads using these relevant climatological variables and the ground snow depth.
For example, Fridley et al. (1994) considered the measurements from a so-called first-order station (stations at which the SWE is recorded) in the United States and modeled the daily SWE using snow depth and daily average temperature. The proposed model was used to estimate ground snow loads for sites associated with a second-order station, where measurement of the SWE is unavailable. Jonas et al. (2009) developed a snow density model for sites in Switzerland using season, snow depth, site altitude, and site location as input variables. Their analysis represented the relationship between snow density and snow depth by a linear function with model coefficients, which depend on the season, altitude, and site. They indicated that the model variability was of the same order of magnitude as that obtained based on measured snow density at a site. The model developed in Sturm et al. (2010) includes the effects of snow depth, snow aging-that is, day of the year (DOY)-and region. The model parameters were estimated by using a large number of snow measurement records from the United States, Canada, and Switzerland. The developed model allows the SWE to be estimated from the observed snow depth. Hill et al. (2019) regressed the SWE against snow depth, DOY, winter precipitation, and temperature using the measured SWE and the corresponding snow depth from the western United States and British Columbia. In their model, the snow season was separated into the accumulation phase and the ablation phase. An equation was given to evaluate the snow density for each phase. In both phases, the SWE was modeled as a power function of snow depth, DOY, winter precipitation, and temperature. They indicated that their model outperforms some other available models in the literature. A more extensive review on the estimation of the SWE from snow depth is available in Bean et al. (2021), who also proposed an algorithm based on a random forest model and found that it outperforms all existing methods for the U.S. dataset they considered. The random forest model is based on the machine learning technique, which is actually a regression approach (Breiman 2001) and this implies that it is not physically-based.
To develop a physically-based algorithm for evaluating the SWE by using ground snow depth measurements, none of the outlined studies directly used precipitation data as the main input variable to estimate the SWE, although precipitation is an important indicator of the SWE. This may be partly due to the difficulties in collecting the gauge precipitation measurements, determining the fraction of precipitation that falls as snow, and estimating the water loss due to melting and/or evaporation. If the winter precipitation in cold regions represents the snowfall, especially during the accumulation period, the precipitation increment should be equal to the SWE increment. In such a case, if the water loss due to melting and/or evaporation could be estimated, the precipitation data could be directly used to estimate the SWE.
In this study, a new practical algorithm to estimate the SWE is proposed. The algorithm is based on frequently collected records of climatological elements such as precipitation, wind speed, temperature, and ground snow depth. The algorithm uses adjusted precipitation data and an estimation of snow water loss to calculate the SWE. The developed algorithm is applied to estimate the SWE for 57 meteorological stations located in the northeastern region of China. The annual maximum SWE obtained by using the proposed algorithm was compared with that obtained from the actual SWE measurements. Moreover, the application of the proposed algorithm to estimate the SWE for sites without actual SWE measurements was carried out. The obtained annual maximum SWE was used to map the snow hazard for the considered region. The climatological data, as well as the proposed algorithm, are described in the following sections. This is followed by the validation and application of the algorithm in estimating the SWE and mapping the snow hazard.

Climatological Data and Data Processing
The location of 99 meteorological stations for the considered northeastern region of China, where the required climatological data are available, is shown in Fig. 1a. These stations provided a good spatial coverage for the region. The relevant climatological data for this study-including (accumulated ground) snow depth, (accumulated ground) SWE, precipitation, wind, and temperature-from 1971 to 2018 were acquired from the China Meteorological Administration (CMA). 1 A review of the data indicated that SWE data are available for the 57 meteorological stations shown in Fig. 1b, while the other variables are available for all 99 stations. The considered region is one of the main regions in China associated with heavy snowfall. The mean winter temperature (averaged for December, January, and February) for some areas of this region is as low as −27 °C (Fig. 1c), and the annual mean winter total precipitation could reach 41.5 mm for some locations (Fig. 1d). For such a region, estimating ground snow load is particularly important for structural Int J Disaster Risk Sci design and mitigating risk due to snow load since the ground snow load is used as the basis for assigning the design roof snow load.

Snow Data
The snow data used in this study include the ground snow depth and snow weight per unit area (referred to as SWE hereafter, as it can be calculated to evaluate snow load considering that the water density is a constant) reported for the 57 stations shown in Fig. 1b. The procedure and requirements for measuring the ground snow depth and the SWE in China are specified by the Chinese specifications for surface meteorological observation (CMA 2017a). According to the specification, the ground snow depth is measured daily whenever the ground snow covers more than half of the surrounding ground near the observation field. The SWE is measured every 5 days if the ground snow depth is greater than 5 cm. The measurements are carried out at 8:00 for most cases and at 14:00 or 20:00 h in some special circumstances. For example, if there is no snow on the ground at 8:00 but later snow on the ground reaches the measurement criterion due to new snowfall, a measurement is carried out. Independent of the time when the snow depth is measured, the measured value is assigned to the same measuring day. To partly exclude possible recording and/or measurement errors, records with associated snowpack bulk density (derived from concurrent measurement of the SWE and snow depth) greater than 400 kg/m 3 or lower than 50 kg/ m 3 are eliminated. These adopted bounds are based on the findings of Yang, Wang et al. (1992) and Chen et al. (2011) for the climate in China. Also, erroneous snow Benxi, respectively); c Spatial distribution of mean winter temperature; and d Spatial distribution of annual mean winter total precipitation. depth records (that is, unreasonably high values as compared to the snow depth recorded for the previous and following days) are removed.

Precipitation, Wind, and Temperature Data
Precipitation, wind, and temperature data from 1971 to 2018 were also obtained from CMA for the 99 sites in the study region. Precipitation in China, both solid and liquid, is measured using the Chinese Standard Precipitation Gauge (CSPG), a cylinder of galvanized iron that is 65 cm long and 20 cm in diameter. The gauge is placed at a site without a windshield and at 0.7 m height above the ground surface (Ye et al. 2004;CMA 2017b). According to the Chinese specifications for surface meteorological observations (CMA 2017b), precipitation measurement is carried out at 8:00 and 20:00 h, respectively, for the previous 12 h. Precipitation measured at 8:00 and 20:00 h is denoted as P 20-8 and P 8-20 , respectively, in the following.
It is worth noting that gauge measurements of precipitation have long been recognized as underestimating the actual precipitation due to wind-induced undercatch, wetting losses, evaporation losses, and trace precipitation, among other causes (Larson and Peck 1974;Yang et al. 1995;Goodison et al. 1998;Sevruk et al. 2009;Rasmussen et al. 2012). The underestimation is even more severe for solid precipitation (snow) than for rainfall, as the snowfall trajectory is more easily influenced by the wind. Goodison et al. (1998) indicated that snow measurement using precipitation gauges has systematic losses of up to 100% caused by wind, wetting, and evaporation effects. Consequently, precipitation data should be adjusted or corrected before they can be used for estimating the SWE.
Wind speed is mainly measured by using cup anemometers installed at 10 m height above the ground surface in open country exposure (CMA 2017c). For each station, the recorded data include daily mean wind speed, daily maximum 10-min mean wind speed, and daily maximum 3-s gust mean wind speed. The daily mean wind speed is the average of the 2-min mean wind speed measured at 2:00, 8:00, 14:00, and 20:00 h. Due to the lack of a wind speed record with a higher sampling rate during snowfall events, only daily mean wind speed is used in adjusting the precipitation data in the proposed algorithm to estimate the SWE.
Similarly, the record for the air temperature contains only the daily maximum, minimum, and mean air temperature. Because continuous air temperature measurements are unavailable, only daily mean air temperature is used in the proposed algorithm to estimate the SWE.

Proposed New Algorithm for Estimating the Snow Water Equivalent
Details of the proposed algorithm, including an overview, steps, assumptions, and formulations, are presented and described in this section.

Overview of the Proposed New Algorithm
The mass conservation equation for the snowpack of a unit volume can be written as: where ΔW denotes the change of the SWE, P s is precipitation in the form of snowfall, P rs represents precipitation in the form of rainfall and absorbed by the snowpack, V is mass transport due to wind, E s denotes water loss due to evaporation or sublimation, and M denotes water loss due to snowmelt. As a simplification, a water-holding capacity of 3% (of the snowpack volume) (Zhou et al. 2013) is considered in the following to take into account the effect of water absorption by the snowpack in the case of rainfall. By considering that the wind-induced snow transport reaches an equilibrium state for open country exposure, V in Eq. 1 is neglected in the following. The sum of the water loss due to evaporation or sublimation, E s , and due to snowmelt, M, is treated as a single variable and is assumed to be proportional to the reduction of ground snow depths. Based on these assumptions, Eq. 1 is simplified to: where ΔW P represents the increment of the SWE, which includes precipitation falling as snow, P s , and rain absorbed by the snowpack, P rs ; ΔW L denotes the snow water loss, which is equal to the sum of E s and M. By adopting Eq. 2, the steps included in the proposed algorithm or procedure to estimate the SWE are depicted in the flowchart shown in Fig. 2 and are described in the following.
(1) The snow season is defined as 1 July to 30 June of the next year. Within each considered snow season, precipitation data associated with an observed ground snow depth equal to 0 are discarded, and only those associated with non-zero snow depth are processed further.
(2) Precipitation data associated with non-zero snow depth are treated as snowfall or mixed precipitation. The fractions of snowfall and rainfall are determined accord- ing to the daily temperature and then corrected due to wind-induced undercatch, wetting losses, and trace precipitation.
(3) The daily increment of the SWE is estimated by summing up the snowfall and the rainfall that is absorbed by the snowpack. (4) The daily snow water loss due to melting and/or evaporation is approximated by a decreasing rate multiplying the SWE on the previous day. The decreasing rate of the SWE is assumed to be equal to the decreasing rate of the snow depth after considering a daily densification rate of 5%. The rationale for this step and the use of 5% is explained in Sect. 3.4. (5) The daily SWE is estimated by substituting the results of steps 3 and 4 in Eq. 2.
These evaluation steps are simple but associated with uncertainties as discussed in the following sections.

Adjustment for Precipitation Data
The gauge catch deficiencies have been investigated by different researchers, and many articles have been published on this subject (Larson and Peck 1974; Rasmussen et al. 2012;Zhang et al. 2019). The precipitation gauge measurement error increases as the wind speed increases. The error is greater for solid precipitation (snow) than liquid precipitation (rain) because the trajectory of the snowflake is easier to be influenced by the wind (Larson and Peck 1974;Yang et al. 1995;Goodison et al. 1998;Nitu et al. 2018). Consequently, it is necessary to determine whether the precipitation is rainfall or snowfall to account for the possible measurement errors. However, information on precipitation type is not available in the metadata acquired from CMA. To deal with this problem, in this study, the precipitation is classified as rainfall if there is no snow on the ground (ground snow depth is zero). When Correct precipitation on day j using wind data and Eqs. 4-7 Determine increment of snow water W pj using Eq. 8 Determine depth decreasing rate DR dj using Eq. 9 Determine snow water loss W Lj using Eq. 10 Determine SWE on day j using Eq. 2 Last day of snow season?

End
No the ground snow depth measurement is greater than zero, the precipitation is treated as snowfall or mixed precipitation, depending on the climatic conditions. This mixed precipitation is represented by fractions of snowfall and rainfall, where the fractions could be evaluated as suggested by Ye et al. (2004) and Zhang et al. (2019). More specifically, precipitation is classified as snowfall if the daily mean air temperature, T d , is lower than −2 °C, classified as rainfall if T d is higher than +2 °C, and classified as mixed precipitation if T d is between −2 °C and +2 °C. A linear function of temperature is adopted to estimate the fraction of precipitation that has fallen as snow, F s : where F s is equal to 1 when T d is lower than −2 °C and 0 if T d is higher than 2 °C. The fraction of precipitation that has fallen as rain, F r , is equal to (1−F s ).
The physical mechanisms that control the precipitation phase are much more complex than what the aforementioned model represents. Factors that affect the precipitation phase may include stability of the atmosphere, temperature lapse rate, the atmospheric humidity profile, and many others Li et al. 2022). However, these factors are rarely adopted in applications since the implementation of these mechanisms would increase the model complexity and thus introduce more uncertainties . Consequently, the most popular methods to identify the precipitation type in practical applications are those based on near-surface air temperature (Li et al. 2022), which is also considered in this study.
Since this study is focused on the annual maximum values of the SWE, the uncertainty caused by the partition of precipitation may not be very important. This is because the mixed snow-rain events occur only in pre-winter or postwinter periods that are unlikely to affect the annual maximum values of the SWE. Therefore, neglecting the uncertainty in evaluating F s is justified for assessing the snow load for structural design.
To determine the wind-induced error in measuring solid precipitation, the World Meteorological Organization (WMO) Solid Precipitation Measurement Intercomparison was initiated in 1985. Twelve countries, including China, participated and submitted complete data summaries for analysis in the intercomparison (Goodison et al. 1998). The intercomparison results were described in Yang et al. (1991) by using data from Tianshan glaciological station, China, from July 1987 to August 1991. The results indicate that, when compared to the double fence intercomparison reference, the average catch ratio (CR) for rain, rain and snow, wet snow, and dry snow is 95.82%, 89.03%, 83.62%, and 72.62%, respectively (Yang et al. 1991;Goodison et al. 1998). In general, CR depends on the average wind speed during the snow or rain event, U (m/s). According to Yang et al. (1991), it can be estimated using: and, By using the catch ratio, the precipitation is corrected or adjusted according to: for the snow and, for the rain, where P sc and P rc denote corrected snowfall and rainfall, respectively; P g is precipitation measured by the gauge; CR s and CR r are given in Eqs. 4 and 5, respectively, and F s is the fraction of precipitation that represents the snowfall given in Eq. 3. In Eqs. 6 and 7, P w denotes wetting loss that is caused by precipitation that is retained or stuck to the internal wall of the gauge and cannot be measured in precipitation observation. A value of 0.3 mm is taken for P w in this study (Yang et al. 1991;Ye et al. 2004;Zhang et al. 2019). Trace precipitation P t could not be measured by the gauge because of the resolution of the gauge, and a value of 0.1 mm is assigned for it when a trace precipitation event is recorded (Ye et al. 2004;Zhang et al. 2019).
Evaporation loss of precipitation measurements is ignored in this study as it is found that evaporation could be negligible for the Chinese Standard Precipitation Gauge (CSPG), which uses a funnel and a container (Ye et al. 2004;Ren and Li 2007).

Estimation of Snow Water Increment
A water-holding capacity of 3% (Zhou et al. 2013) is considered in this study. The maximum rainfall absorbed by the snowpack is then estimated by multiplying 0.03 and ground snow depth in mm. The increment of the SWE on day j, ΔW Pj , is then given by: where d j represents the snow depth on day j in mm, and P scj and P rcj are corrected snowfall and corrected rainfall on day j, respectively.

Estimation of Snow Water Loss
Water loss of snowpack is mainly caused by snowmelt, snow evaporation and/or snow sublimation. Male and Gray (1975) indicated that the net amount of evaporation/condensation over a 24-h period during the snowmelt is generally negligible. But accurate estimation of the snowmelt is important for water resources management and flood risk assessment. Thus, the modeling of the snowmelt process has long been studied in hydrology (Anderson 1968;Male and Gray 1975;Kustas et al. 1994;Boudhar et al. 2016). Snowmelt models are mainly based on the energy balance equation of the snowpack (Kondo and Yamazaki 1990; Tarboton and Luce 1996;Zeinivand and De Smedt 2010). A simplified model such as the degree-day model has also been developed and applied (Johannesson et al. 1995;Braithwaite and Zhang 2000).
Unfortunately, the required input data and/or model parameters for the energy balance equation or degreeday model are not reported for ordinary meteorological stations. So rather than use the degree-day model or the energy balance model to estimate the SWE, the proposed algorithm in this study assumed that snow water loss is proportional to the reduction of snow depth.
Note that there is a densification process in the snowpack due to gravity, wind, precipitation, and temperature gradient Sturm and Holmgren 1998). Therefore, a decrease in snow depth does not necessarily imply a reduction in the SWE. Sturm and Holmgren (1998) showed that the average densification rate of snow for regions classified as Maritime, Tundra, and Taiga is 1.31, 0.24, and 0.57 kg/m 3 /d, respectively. The densification rate for freshly fallen snow is relatively higher because the fresh snow is looser and compacts more easily.  indicated that the snow density of fresh snow increases by about 5% on the second day. Considering the typical average density of 70-110 kg/ m 3 for fresh snow, the corresponding densification rate is 3.5-5 kg/m 3 /d. This is close to the 4 kg/m 3 /d given by Wei et al. (2001). However, these values are much lower than the 9.6 kg/m 3 /d given by Chen et al. (2011). The wide range of snow densification rates indicates that snow densification varies significantly and is complicated to predict. As a simplification, an average densification rate of 5% is adopted in this study. This results in a reduction of approximately 5% in snow depth each day with the SWE remaining the same as the previous day.
Densification is assumed to last for 10 days after the occurrence of the snowfall event, and the final density is equal to 1.05 10 (or 1.63) times that of the original density. For a typical fresh snow density of about 70-110 kg/m 3 , the corresponding final density after densification is then estimated to be 115-180 kg/m 3 , which is consistent with the observed values in China (Mo et al. 2016;Mo, Cao et al. 2022). Based on this assumption, the decreasing rate of snow depth on day j, DR dj , is given by the following equation: where d j−1 and d j denote ground snow depth on day j−1 and day j, respectively, and Δd allow,j is the allowance of snow depth reduction on day j that is the summation of potential depth reduction for all snow events within 10 days. Snow water loss on day j, ΔW Lj , is then approximated by applying the decreasing rate to the SWE on the previous day: In addition, the criterion for the limiting values of the density that was applied to the original snow data is also applied to the estimated SWE. That is, the associated snowpack bulk density that is calculated by dividing the derived SWE by the recorded ground snow depth for the same day should not be less than 50 kg/m 3 or greater than 400 kg/m 3 .

Results and Discussion
Based on the proposed procedure, the daily SWE was estimated from the ground snow depth and climatic data for the considered stations shown in Fig. 1b. The estimated daily SWE, as well as the annual maximum SWE, are compared to the measured ones to validate the proposed algorithm. The snow hazard map for northeastern China is then developed using the measured and calculated SWE.

Comparison of the Calculated and Measured Snow Water Equivalent
Examples of such a comparison are shown in Fig. 3 for four representative sites identified in Fig. 1b, where the data points represent the daily SWE for all available years  for the corresponding station. When plotting the figures, data points that are associated with 0 value in both the calculated and measured SWE are eliminated. Figure 3 shows that the calculated and measured SWE values are in good agreement as they are concentrated near the diagonal line. Their correlation coefficients, r, are greater than 0.83, and the bias ranges from 0.2 to 3.1mm. Differences between the calculated SWE and SWE measurements could also be observed. These differences could be attributed to the following five causes: (1) The difference or gap between the measured SWE and precipitation. In some cases, the gap between recorded precipitation and the SWE is too large to be closed by the correction of the precipitation data. For example, at Site 2, the measured SWE for 30 November and 5 December 1987 is 40 and 51 mm, respectively, which indicates that a total of 11 mm of SWE due to snowfall occurred within the five days. But the corresponding recorded precipitation within these five days was only 0.1 mm and two trace precipitation events, or only 0.6 mm in total after wetting loss correction and trace precipitation correction. This leaves a large difference between the precipitation and the measured SWE.
(2) The criteria for the SWE measurements. The SWE is measured following specified criteria (that is, every 5 days if the ground snow depth is greater than 5 cm; additional measurement is carried out only if the depth increment due to new snowfall exceeds 5 cm). If ground snow depth is less than 5 cm, no SWE is reported even though the actual SWE may not be zero. The absence of SWE measurement, in this case, will certainly cause bias between the recorded and the calculated SWE. (3) The uncertainties caused by the determination of the snow fraction in total precipitation. The variability in the difference between the calculated and the measured SWE could be caused by the uncertainty in the calculated fraction of snowfall (see Eq. 3). Part of this uncertainty is attributed to the fact that the continuous air temperature variation data during a day is unavailable, and only daily maximum, minimum, and average temperatures are reported. For example, a snowfall of 15 mm was recorded on 4 April 1999 for Site 3; the average temperature for that day is 0.6 °C, resulting in a F s of 0.4 according to Eq. 3. The precipitation that fell as snow is estimated to be only 5.4 mm in this case, or 5.8 mm after correction, which is much less than the measured SWE of 15 mm. Such a discrepancy could be reduced or eliminated if the continuous air temperature was recorded, reported, and used to calculate the SWE. (4) The uncertainties caused by the precipitation correction model shown in Eqs. 4-7. The sampling frequency of climatological elements such as wind speed and precipitation is very low. This could affect the calculated (that is, corrected or adjusted) precipitation amount. This problem could only be investigated using data with a higher sampling frequency, which was unavailable for this study. (5) The uncertainties caused by the estimation of snow water loss. The estimation model for snow water loss given by Eqs. 9 and 10 assume that the snow water loss is proportional to the decrease of snow depth upon the consideration of a densification rate of 5%. However, as pointed out in Sect. 3.4, snow densification is a complicated process and varies significantly with time and location. The use of a single 5% value of the densification rate will bring uncertainties to the estimation of snow water loss and further to the estimation of the SWE.
The above observations suggest that a better estimate of the daily SWE may be obtained if more detailed measurements of the climatological elements for a site of interest are reported and used with the proposed algorithm.

Comparison of the Calculated and Measured Annual Maximum Snow Water Equivalent
This section presents the comparison of the annual maximum SWE extracted from the calculated and measured SWE values. The comparison is focused on the probability distribution and the T-year return period value of the annual maximum SWE. For each of the 57 meteorological stations depicted in Fig. 1b, a probability distribution fitting was carried out for the extracted annual maximum SWE. For the distribution fitting, the Gumbel distribution and the lognormal distribution were considered as the candidate probability distribution. The consideration of the Gumbel distribution is justified since it is recommended by the Chinese structural design code (MOHURD 2012) and used for the structural design code calibration for other jurisdictions (Bartlett et al. 2003;BSI 2003;NRCC 2015;Hong et al. 2021;Mo, Cao et al. 2022). The lognormal distribution is considered because it is preferred for the annual maximum ground snow depth or SWE (Ellingwood and Redfield 1983;ASCE 2016;Mo et al. 2016;Mo, Ye et al. 2022).
The Gumbel distribution can be written as (Coles 2001): where α and μ are the scale and location parameters, respectively. The lognormal distribution can be written as: where m lnx and σ lnx are the distribution parameters that represent the mean and standard deviation of ln(x), respectively, and Ф( ) denotes the standard normal distribution function.
To identify the preferable probability model for the annual maximum SWE, the Akaike Information Criterion with a small sample size (AICc) (Burnham and Anderson 2003) is considered: where k is the number of model parameters, n is the sample size, and L is the maximum likelihood of the considered model. The preferred model among the candidate models is the one that leads to the lowest AICc value. The distribution fitting results indicate that the Gumbel distribution and lognormal distribution are preferred by 18 and 39 stations, respectively, if the measured SWE is considered, and by 9 and 48 stations, respectively, if the estimated SWE is considered. It was observed that in many cases, the differences in AICc for the considered two models are not very large. The preference for the lognormal distribution is consistent with the observation made by Mo et al. (2016), which indicated that the lognormal distribution is preferable for the annual maximum SWE or annual maximum ground snow depth from about two-thirds of the sites considered in their study.
Since there is no clear spatial trend on the preferred distribution of the SWE, the preferred lognormal distribution is considered in the following. Moreover, since the use of the method of L-moments (Hosking and Wallis 1997) is recommended in the literature, it is used to fit the annual maximum SWE. The empirical and the fitted distributions for the four representative sites identified in Fig. 1b are presented in Fig. 4. Figure 4 indicates that the empirical and fitted probability distributions of the annual maximum SWE extracted from the calculated SWE match well with those of the annual maximum extracted from the measured SWE, especially for Sites 1 and 4. The differences between the empirical distributions based on the calculated and measured data are appreciable. As has been pointed out in the previous section, the differences are attributed to various reasons and could be reduced if more detailed measurements of the climatological elements for a site of interest are reported and used with the proposed algorithm..
To better quantify the impact of using the SWE calculated using the proposed algorithm in estimating the snow hazard, the 50-year return period value of the annual maximum SWE, s 50 , for the considered four sites in Fig. 1b, was calculated based on both the calculated and measured SWE.
The calculated values are shown in Table 1. Note that s 50 shown in the table is presented in kPa, which is converted from mm by considering a water density of 1,000 kg/m 3 (13) AICc = 2k − 2 ln L + 2k 2 + 2k n − k − 1 and a gravitational acceleration of 9.8 m/s 2 . The comparison indicated that the relative difference-that is, (s 50 based on the calculated SWE-s 50 based on the measured SWE) / (s 50 based on the measured SWE)-is site-dependent; the absolute value of the relative error is less than about 17% for the four representative sites. The analysis that is carried out for the four representative sites is repeated for each of the 57 meteorological stations depicted in Fig. 1b. The obtained results are shown in Fig. 5a. The figure indicates that the proposed algorithm estimates s 50 fairly well for locations with s 50 less than about 0.8 kPa, but underestimates the s 50 for locations where s 50 is higher than about 0.8 kPa. More specifically, a histogram of the estimation error, Err, (Fig. 5b) shows that the distribution of the error is left-skewed, implying that the proposed algorithm tends to underestimate the s 50 in general. A correction to the estimated s 50 is desired to eliminate such a trend.
For the purpose of correcting the estimated s 50 , attempts are made to plot Err versus the 30-year normal of winter precipitation (December, January, February), normal of maximum winter wind speed, and normal of mean coldest month  temperature. Unfortunately, no trend is observed. However, when Err is plotted versus s 50 estimated from the measured SWE (Fig. 6a), an apparent decreasing trend for s 50 greater than about 0.6 kPa, and an increasing trend for s 50 less than 0.6 kPa could be observed. Although the reasons for these trends are not clear in this study-and whether the trends are applicable to other regions requires further investigation when snow data for other regions become available-the trends provide a basis for the correction of the estimated s 50 from a practical point of view, at least for the region of interest.
Since the proposed algorithm is designed for stations where SWE measurements are not available, the correction is made according to 50-year ground snow depth, d 50 , instead of s 50 . The variation of Err with d 50 is depicted in Fig. 6b. Comparing Fig. 6a and b shows that the trends of Err in Fig. 6b are quite similar to those in Fig. 6a. Regression analysis is carried out to model Err and it is found that a polynomial of order 2 is adequate to model the error: where μ E is the modeled value of Err (see the red line in Fig. 6b), ε is the residual of the regression, and Err and d 50 are in the unit of kPa and m, respectively. The regression residual ε is normally distributed (Fig. 7b). The corrected s 50 , s 50,cor , could then be expressed as: where s 50,cal is the s 50 estimated from the calculated SWE. Similar to Fig. 5a, the corrected estimation of s 50 is compared to s 50 estimated from the measured SWE in Fig. 7a. The figure shows that the corrected estimation of s 50 follows the diagonal line fairly well, implying that a good estimation  of s 50 was obtained. It also indicates that the correction improves the estimation effectively when compared to Fig. 5a. The estimation error of the corrected s 50 is the same as the residual shown in Fig. 7b, which could be assumed to be normally distributed.
To see the ground snow load in the context of the structural design code making, it is noted that the Chinese load code for the design of building structures (MOHURD 2012) recommends that s 50 could be evaluated based on the average snowpack density and the 50-year annual maximum ground snow depth when SWE measurements are not available: where ρ avg is the applicable average snowpack density for the considered site, which is equal to 150 kg/m 3 for northeastern China, g = 9.8 m/s 2 is the gravitational acceleration, and d 50 is the 50-year return period value of the annual maximum snow depth, which is calculated by assuming that the annual maximum ground snow depth is Gumbel distributed (as recommended in the code). The obtained s 50 by using Eq. 16 is also included in Fig. 7a for the considered sites. The figure indicates that the error of the estimated s 50 based on the SWE evaluated using the proposed algorithm, in general, is smaller than that suggested by the code procedure. This suggests that for sites without measured SWE, the use of the proposed algorithm to evaluate the SWE can result in a more accurate estimate of s 50 .

Mapping the Snow Hazard Using the Measured and Calculated Snow Water Equivalent
To illustrate the application of the proposed algorithm in snow hazard mapping, first, the calculated s 50 -by using the (16) s 50 = avg gd 50 measured SWE for the 57 considered meteorological stations depicted in Fig. 1b-is presented in Fig. 8a to map the snow hazard. Second, the calculated s 50 for the same stations but considering the SWE evaluated based on the proposed algorithm is presented in Fig. 8b. Third, to improve the spatial coverage, s 50 -for all stations shown in Fig. 1a-was estimated by using the SWE evaluated based on the proposed algorithm. The obtained s 50 was then mapped and is shown in Fig. 8c. Fourth, the snow hazard was mapped as shown in Fig. 8d by using s 50 , which was estimated based on the measured SWE if the measurements are available from a station; otherwise it was estimated based on the SWE calculated using the proposed algorithm. Finally, s 50 -calculated according to the structural design code-suggested practice (see Eq. 16)-was mapped in Fig. 8e. Four major observations could be drawn from Fig. 8: (1) The spatial trends and the magnitudes of the mapped s 50 estimated based on the measured SWE (Fig. 8a) and based on the calculated SWE by applying the proposed algorithm (Fig. 8b) are consistent except the outlier located in the central-east of the region, where the s 50 is underestimated significantly by the proposed algorithm.
(2) The use of s 50 by considering all stations shown in Fig. 1a improves the granularity of the mapped hazard (by comparing Fig. 8b and c), indicating the usefulness of applying the proposed algorithm in augmenting the SWE data. (3) The use of s 50 estimated based on the measured SWE and the calculated SWE, when the measured SWE is unavailable, leads to a snow hazard map with refined resolution. A comparison of Fig. 8a and d indicates that the consideration of extra stations could reduce the influence of outliers on the hazard map. (4) Comparison of Fig. 8a-e indicates that the use of the code-recommended approach results in the mapped hazard that differs from those obtained based on the measured SWE or calculated SWE. It indicates that the code-recommended approach (MOHURD 2012) should be reassessed.

Conclusions
For snow hazard-prone regions, the estimation of extreme ground snow loads for locations where only ground snow depths are available is of great importance. Given the ground snow depth alone, this estimation is challenging since the snowpack bulk density is uncertain and with significant temporal and spatial variability. By considering the relationship between the snow water equivalent (SWE) and precipitation, a new algorithm or procedure is proposed to estimate the SWE (or ground snow load). Using the proposed algorithm, the ground snow depth records, and the temperature, precipitation, and wind speed records, the SWE (or ground snow load) history was calculated and compared with the recorded SWE values for 57 locations in northeastern China. Also, an analysis of the 50-year return period value of the ground snow load was carried out based on the measured SWE, the calculated SWE using the proposed algorithm, and the procedure recommended by the Chinese load code. Snow hazard maps were developed based on the estimated s 50 . Three main conclusions can be drawn from the numerical results: (1) In general, the SWE obtained by using the proposed algorithm agrees well with the measured SWE.
(2) The spatial trends and the magnitudes of the mapped s 50 estimated based on the measured SWE and based on the calculated SWE by applying the proposed algorithm are consistent, indicating the usefulness of applying the proposed algorithm in augmenting the SWE data.
(3) The error analysis and comparison of the mapped snow hazards indicate that the code-recommended approach should be reassessed by the committee that reviews the snow load for the load code (MOHURD 2012) since it leads to the 50-year return period value of the snow load deviating significantly from that estimated based on the measured SWE.  Fig. 1a; d Using s 50 estimated based on the calculated and measured SWE (when available); and e Using s 50 estimated based on code-suggested procedure (see Eq. 16).