Reconstruction of a daily gridded snow water equivalent product for the land region above 45° N based on a ridge regression machine learning approach

The snow water equivalent (SWE) is an important parameter of surface hydrological and climate systems, and it has a profound impact on Arctic amplification and climate change. However, there are great differences among existing SWE products. In the land region above 45^∘ N, the existing SWE products are associated with a limited time span and limited spatial coverage, and the spatial resolution is coarse, which greatly limits the application of SWE data in cryosphere change and climate change studies. In this study, utilizing the ridge regression model (RRM) of a machine learning algorithm, we integrated various existing SWE products to generate a spatiotemporally seamless and high-precision RRM SWE product. The results show that it is feasible to utilize a ridge regression model based on a machine learning algorithm to prepare SWE products on a global scale. We evaluated the accuracy of the RRM SWE product using hemispheric-scale snow course (HSSC) observational data and Russian snow survey data. The mean absolute error (MAE), RMSE, R, and R² between the RRM SWE products and observed SWEs are 0.21, 25.37 mm, 0.89, and 0.79, respectively. The accuracy of the RRM SWE dataset is improved by 28 %, 22 %, 37 %, 11 %, and 11 % compared with the original AMSR-E/AMSR2 (SWE), ERA-Interim SWE, Global Land Data Assimilation System (GLDAS) SWE, GlobSnow SWE, and ERA5-Land SWE datasets, respectively, and it has a higher spatial resolution. The RRM SWE product production method does not rely heavily on an independent SWE product; it takes full advantage of each SWE dataset, and it takes into consideration the altitude factor. The MAE ranges from 0.16 for areas within <100 m elevation to 0.29 within the 800–900 m elevation range. The MAE is best in the Russian region and worst in the Canadian region. The RMSE ranges from 4.71 mm for areas within <100 m elevation to 31.14 mm within the >1000 m elevation range. The RMSE is best in the Finland region and worst in the Canadian region. This method has good stability, is extremely suitable for the production of snow datasets with large spatial scales, and can be easily extended to the preparation of other snow datasets. The RRM SWE product is expected to provide more accurate SWE data for the hydrological model and climate model and provide data support for cryosphere change and climate change studies. The RRM SWE product is available from “A Big Earth Data Platform for Three Poles” (https://doi.org/10.11888/Snow.tpdc.271556) (Li et al., 2021).

1 Introduction

The IPCC (Intergovernmental Panel on Climate Change) AR6 (Sixth Assessment Report) notes that the Northern Hemisphere spring snow cover has greatly decreased since 1950, and the feedback effect of the climate system caused by this reduction is extremely large (Masson-Delmotte et al., 2021). In most land areas of the Northern Hemisphere, annual runoff is dominated by snowmelt, and accurately estimating the impacts of such a large amount of snowmelt runoff on ecosystems and human activities is of great significance (Barnett et al., 2005; Bintanja and Andry, 2017; Henderson et al., 2018). Whether through hydrometeorological simulation or global change research, the estimation of the energy budget and mass of snow is very difficult, so a set of highly accurate, long time series snow cover datasets is urgently needed to drive hydrometeorological simulations and land surface process models. Among them, snow water equivalent (SWE) data play an irreplaceable role as an important parameter of the land surface hydrological model and climate model.

At present, there are many forms of SWE data in the world. According to type, these data can be divided into site observational SWE, remote sensing SWE, reanalysis SWE, data assimilation SWE, and model simulation SWE. The remote sensing SWEs are mainly AMSR-E (Kelly, 2009) and AMSR2 (Imaoka et al., 2010; Tedesco and Jeyaratnam, 2019). The reanalysis SWE was mainly based on the ERA-Interim (Dee et al., 2011), MERRA2 (Gelaro et al., 2017), MERRA land (Reichle et al., 2011), and ERA5-Land (Muñoz Sabater, 2019; Balsamo et al., 2015) datasets. The data assimilation SWE mainly includes GlobSnow (Luojus et al., 2021) and the Global Land Data Assimilation System (GLDAS) (Rodell et al., 2004). The site observational SWE mainly includes the GHCN dataset (Menne et al., 2016) and HSSC data (Pulliainen et al., 2020). However, the time ranges of AMSR-E and AMSR-E2 SWE are only from 2003 to the present, which is lacking in terms of time series. Similarly, the GlobSnow SWE dataset is also seriously lacking in time series. Although the reanalysis SWE data have good spatial and temporal continuity and high data integrity, their accuracy is poor, and the mean absolute error (MAE) is 0.65 (Snauffer et al., 2016). The SWE data from stations and meteorological observations cannot meet the needs of hydrometeorological and climate change research. This is mainly because SWE from stations is discontinuous in time series and severely missing. Furthermore, hydrometeorological studies often require spatiotemporally continuous grid data to be derived (Pan et al., 2003). There are great differences among remote sensing SWE, reanalysis SWE data, data assimilation SWE, and observational SWE. For remote sensing SWE, the spatiotemporal characteristics of different passive microwave SWE data differ significantly due to differences in sensors or retrieval algorithms (Mudryk et al., 2015). Data assimilation SWE and reanalysis SWE data also tend to exhibit different spatiotemporal characteristics due to differences in model design, driving data, and assimilation methods (Vuyovich et al., 2014). In summary, although there are a variety of SWE data in the world, the data quality is uncertain.

Previous studies have shown that all kinds of SWE data in the Northern Hemisphere have advantages and disadvantages, and none of these data perform well in all aspects (Mortimer et al., 2020). An effective method was used in a study by Pulliainen et al. (2020), who applied a bias correction to GlobSnow and reanalysis data products based on SWE snow course measurements to obtain improved estimates on annual peak snow mass and SWE in the Northern Hemisphere. Another effective method is to fuse all kinds of SWE data in time and space, integrate the advantages of all kinds of data, and then generate a relatively complete SWE dataset. Many scholars have conducted in-depth studies on SWE data fusion. The main fusion methods can be classified into the following categories: multiproduct direct averaging (Mudryk et al., 2015), linear regression (Snauffer et al., 2016), data assimilation (Pulliainen, 2006), “multiple” collocation (Pan et al., 2015), and machine learning (Snauffer et al., 2018; Xiao et al., 2018; Wang et al., 2020). Studies have shown that even the simplest multisource data average is more accurate than a single SWE product (Snauffer et al., 2018). However, the simple multisource data average cannot highlight the advantages of high-precision data, and it is easily affected by the weight ratio of low-precision data, which reduces the accuracy of fused data (Mudryk et al., 2015). Although the linear regression method can make good use of the actual observational data to correct the original data, it is easy to overfit which causes the overall deviation (Snauffer et al., 2016). The “multiple” collocation method changes the size of the original SWE data before fusion, which easily causes data errors. The data assimilation method is sensitive to the accuracy of input data, and it is difficult to fuse multisource data (Pan et al., 2015). In recent years, machine learning methods have been widely used in data fusion (Santi et al., 2021; Ntokas et al., 2021). Machine learning methods can not only integrate the advantages of multisource data but also make full use of site observational data to train the sample data, which easily generates SWE data products with large spatial scales and long time series (Broxton et al., 2019; Bair et al., 2018).

In summary, based on the existing SWE data products, combining a machine learning algorithm to fuse multisource SWE data is an effective method to prepare SWE products with long time series and large spatial scales and retain the advantages of single SWE data products. The ridge regression model is a biased estimation method specifically designed to address the problem of multicollinear data (Duzan and Shariff, 2015; Saleh et al., 2019). It has good tolerance to “ill-conditioned” data and has a good effect in using SWE data to address the multicollinearity problem (Hoerl and Kennard, 1970b; Guilkey and Murphy, 1975). In this study, we integrated multisource SWE data products of the ridge regression model (RRM) SWE based on the ridge regression model of the machine learning algorithm. We selected ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, AMSR-E/AMSR2 SWE, and ERA5-Land SWE data with relatively complete time series as the original data for the production of the RRM SWE product. The missing parts of the ERA-Interim SWE, AMSR-E/AMSR2 SWE, and GlobSnow SWE data were filled by the spatiotemporal interpolation method. The HSSC dataset (Pulliainen et al., 2020) and Russian snow survey data (Bulygina et al., 2011) were used as training sample data of “true SWE”, and the effect of altitude on the algorithm was also considered. Thus, we prepared a set of spatiotemporal seamless SWE datasets (RRM SWE) covering the land region above 45^∘ N from 1979 to 2019. The spatial coverage of the RRM SWE product covers all land regions north of 45^∘ N.

2 Data and methods

2.1 Research region

The research region of the RRM SWE product is located in the land region north of 45^∘ N (Fig. 1). This region consists of Asia, Europe, and North America. The land region covers Russia, the United States, Canada, Denmark, Norway, Iceland, Sweden, and Finland. This region has a cold climate and a wide area of snow cover.

Figure 1 The DEM and snow survey stations of the research region. The right panel (b) shows the DEM, and the left panel (a) shows the SWE observational stations. HSSC, hemispheric-scale snow course; RSSD, the Russian snow survey station. The spatial range of the RRM SWE product is consistent with that of the DEM.

2.2 Grid SWE data description

In this study, we utilized ERA-Interim SWE data (Dee et al., 2011), GLDAS SWE data (Rodell et al., 2004), GlobSnow SWE data (Luojus et al., 2021), AMSR-E/AMSR2 SWE data (Tedesco and Jeyaratnam, 2019), and ERA5-Land SWE data (Muñoz Sabater, 2019) as the original input datasets for the fusion data (Table 1).

Table 1 Introduction to the SWE data.

GlobSnow is a dataset of global snow cover and SWEs for the Northern Hemisphere released by the European Space Agency (ESA) (http://www.globsnow.info/swe/,last access: 17 February 2022) (Luojus et al., 2021; Pulliainen et al., 2020). The SWE products in this dataset combine the Canadian Meteorological Center (CMC) daily snow depth analysis data (Walker et al., 2011), ground weather site observational data, and satellite microwave radiometer data. We obtained the L3A_daily_SWE product of this dataset. The temporal resolution of the L3A_daily_SWE product is daily, the spatial resolution is 0.25^∘, and the data format is netCDF4.

ERA-Interim is the fourth generation reanalysis data of the European Centre for Medium-Range Weather Forecasts (ECMWF) (Dee et al., 2011). The data provide a global assimilated numerical product of various surface and top atmospheric parameters from January 1979 to the present (https://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/, last access: 17 February 2022). We obtained the SWE dataset with a daily temporal resolution, a spatial resolution of 0.25^∘, and netCDF4 data format. The spatial range of the data is the land region above 45^∘ N.

The Advanced Microwave Scanning Radiometer-Earth Observation System (AMSR-E) is a microwave scanning radiometer on the Aqua satellite of the National Aeronautics and Space Administration (NASA) Earth Observation System (EOS) (Tedesco and Jeyaratnam, 2019). The AMSR-E provides a global daily SWE dataset from 19 June 2002 to 3 October 2011 (https://nsidc.org/data/ae_dysno, last access: 17 February 2022). AMSR2 is a microwave scanning radiometer on the GCOM-W1 satellite launched by the Japan Aerospace Exploration Agency (JAXA) in May 2012. AMSR2 provides a global SWE dataset from 2 July 2012 to the present (https://nsidc.org/data/AU_DySno/versions/1, last access: 17 February 2022). The spatial resolution of the AMSR-E SWE and AMSR2 SWE datasets is 25 km×25 km, the temporal resolution is daily, and the data formats are HDF-EOS and HDF-EOS5, respectively.

The GLDAS is a model used to describe global land information; it contains data, such as global rainfall, water evaporation, surface runoff, underground runoff, soil moisture, surface snow cover distribution, temperature, and heat flow distribution (Rodell et al., 2004). This assimilation system includes data with spatial resolutions of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ and $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ and temporal resolutions of 3 h, 1 d, and 1 month. The GLDAS data are available for download from the Goddard Earth Sciences Data and Information Services Center (GES DISC). We obtain an SWE dataset with a daily temporal resolution, 0.25^∘ spatial resolution, and netCDF4 data format.

ERA5-Land is a reanalysis dataset that provides the evolution of global land parameter data from 1981 onwards (Muñoz Sabater, 2019). The dataset provides eight types of snow parameter data, including snow albedo, snow cover, snow depth, snowfall, the temperature of the snow layer, snowmelt, snow density, and SWE. This dataset provides a global SWE dataset with an hourly spatial resolution, a temporal resolution of $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$, a temporal coverage of January 1981 to the present, and data formats of GRIB (General Regularly-distributed Information in Binary form) and netCDF4.

To maintain consistency in the spatial and temporal resolutions of the fused data, we unified the ERA-Interim SWE data, GLDAS SWE data, GlobSnow SWE data, AMSR-E/AMSR2 SWE data, and ERA5-Land SWE data into a daily temporal resolution, with a spatial resolution of 0.25^∘ and geographic projection of the North Pole Lambert azimuthal equal area.

2.3 Ridge regression machine learning algorithm for preparing the SWE

In this study, we utilize the ridge regression model of a machine learning algorithm to fuse ERA-Interim SWE data (Dee et al., 2011), GLDAS SWE data (Rodell et al., 2004), GlobSnow SWE data (Luojus et al., 2021), AMSR-E/AMSR2 SWE data (Tedesco and Jeyaratnam, 2019), and ERA5-Land SWE data (Muñoz Sabater, 2019) to generate a set of new RRM SWE datasets. The target reference data in this study are the HSSC dataset and Russian snow survey data. The digital elevation model (DEM) was used as an important environmental feature input to the ridge regression model and was included in the model training. The DEM is an auxiliary terrain feature variable in addition to the five SWE prediction feature variables: AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE.

The ridge regression model is a biased estimate regression method for collinear data analysis (Friedman et al., 2010; Hoerl and Kennard, 1970b, a). By abandoning the unbiasedness of the ordinary least squares, this algorithm can obtain the regression method in which the regression coefficient is more practical and reliable at the cost of losing part of the information and reducing the accuracy. The ridge regression model is flexible in the choice of predictor variables and does not require the predictor and target variables to be independent of each other. It can effectively solve the multicollinearity problem of predictor and target variables, as well as reduce the impact of this problem on the training model (Duzan and Shariff, 2015; Saleh et al., 2019). Generally, reanalysis data based on SWE products cannot make the products and models independent of each other; i.e., they are prone to multicollinearity, which leads to distorted model estimation or difficulty in performing accurate estimations. In contrast, the ridge regression model can successfully solve the multicollinearity problem, i.e., the independence of training products and models. In addition, when integrating multiple SWE products, the accuracy of each SWE dataset is likely to differ. A small change in one of the SWE products involved in the training will cause a significant error in the final calculation results, while the ridge regression model has high accuracy and stability for these “ill-conditioned” SWE data. In addition, the main advantage of this model is that SWE products with long time series and large spatial scales are easy to prepare. The principle equation of the ridge regression model is defined as follows:

$$\begin{array}{}\text{(1)}& {\widehat{\mathit{\beta }}}^{\text{ridge}}={\text{argmin}}_{\mathit{\beta }}\left\{\sum _{i=\mathrm{1}}^N{\left(y_i-{\mathit{\beta }}_{\mathrm{0}}-\sum _{j=\mathrm{1}}^p{x}_{ij}{\mathit{\beta }}_j\right)}^{\mathrm{2}}+\mathit{\lambda }\sum _{j=\mathrm{1}}^p{\mathit{\beta }}_j^{\mathrm{2}}\right\},\end{array}$$

where ${\widehat{\mathit{\beta }}}^{\text{ridge}}$ is the extremum solution function of ridge regression; p is the number of gridded SWE product variables involved in training; x_i are the prediction feature variables, which contain two parts: one set contains the main feature variables of the gridded SWE products, and the other part consists of the DEM auxiliary feature variables; y_i is the observed SWE; λ, β, β_j, and β₀ are the parameters to be solved; $\mathrm{1},\mathrm{\dots },N$ is the sample of the training dataset; and $\mathit{\lambda }{\sum }_{j=\mathrm{1}}^p{\mathit{\beta }}_j^{\mathrm{2}}$ is the penalty function term. The total number of samples N in the training dataset is 271 651. The sample sizes of the training dataset, validation dataset, and test dataset are divided according to the ratio of $\mathrm{7}:\mathrm{2}:\mathrm{1}$, where the numbers of training set, validation set and test set samples are 271 651, 77 614, and 38 807, respectively. The model is developed in Python3, and the model framework is based on the “scikit-learn” machine learning library (https://scikit-learn.org/stable/index.html, last access: 17 February 2022). The code is available upon request.

The integration process of the RRM SWE product (Fig. 2) is described as follows:

1. The original ERA-Interim SWE data, GLDAS SWE data, GlobSnow SWE data, AMSR-E/AMSR2 SWE data, ERA5-Land SWE data, DEM data, unified temporal resolution, spatial resolution, projection, spatial range, and unit are preprocessed.

2. The spatiotemporal interpolation method is used to fill in the missing data of AMSR-E/AMSR2 SWE, ERA-Interim SWE, and GlobSnow SWE in space and time. Based on this method, the missing AMSR-E/AMSR2 SWE data at low latitudes and the missing ERA-Interim SWE and GlobSnow SWE data in the time series are added.

3. The SWE data observed at stations from 1979 to 2014 are used as sample training data, and the AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, ERA5-Land SWE data, and DEM data are input into the ridge regression model of a machine learning algorithm for training. During the RRM model training process, we reconstructed the training data to try to extract training samples that are uniformly distributed spatially as much as possible. First, a scan window of 250 km×250 km (10×10 pixels) was created. Then, each gridded SWE data point participating in training is scanned, and the sample numbers in each scan window are counted. Finally, the mean value n of the sample numbers in all scan windows is taken as the number of training samples to be selected in each scan window. For the scan window with sample numbers higher than n, n samples are randomly selected from the scan window. For the scan window with sample numbers lower than n, all samples in the scan window are selected as training samples.

4. When the model was trained, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE were used as the training data between 1979 and 2002 (AMSR-E/AMSR2 SWE data were not available before 2002), and AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE were used as the training data after 2002.

5. Based on the S-fold cross-validation method, the SWE data are continuously trained and validated, and the optimal model and parameters are finally selected and evaluated by the loss function.

6. Based on the trained optimal model, multiple SWE data products are integrated into the time series, missing data are predicted, and a set of spatiotemporally seamless SWE datasets is generated.

7. SWE data observed at stations from 2015 to 2018 are used to evaluate the accuracy of the RRM SWE product.

Figure 2 Flow chart of the RRM SWE data preparation (preparation of spatiotemporal seamless SWE datasets mainly includes three processes: model training, model reasoning, and SWE data preparation).

2.4 Site data and evaluation metrics

2.4.1 Site SWE data for training, validation, and testing

Russian snow survey data (http://aisori.meteo.ru, last access: 17 February 2022) include the average snow depth data and the average snow density data of the station, and the SWE is the product of the measured average snow depth and average snow density (Bulygina et al., 2011). We obtained SWE data from 19 493 stations from 1979 to 2016 from this dataset.

Hemispheric-scale snow course (hereafter referred to as HSSC) observational data are contained in a hemispheric-scale SWE database based on SWE observational datasets from the former Soviet Union/Russia (FSU), Finland, and Canada developed by Pulliainen et al. (2020) (Bronnimann et al., 2018; Brown et al., 2019). This dataset is from the website of the Finnish Meteorological Institute (FMI) (https://www.globsnow.info/swe/archive_v3.0/auxiliary_data/, last access: 17 February 2022). The dataset provides data from 2687 distributed regional snow course observations and contains 343 241 SWE observational data points from 1979 to 2018. The snow courses of the HSSC dataset are transects in which SWE is sampled manually at multiple locations with typical conditions to eliminate uncertainty in the regional-scale spatial variability in SWE due to the influence of snowpack characteristics and land cover type (Pulliainen et al., 2020).

We carefully screened the Russian snow survey data and HSSC data and eliminated some abnormal observational data to ensure the high quality of the training, validation, and test sets. The null and zero values are removed during the HSSC data screening process. The null values, negative numbers, and extreme SWE values greater than 2000 mm are removed during the Russian snow survey data screening process.

2.4.2 Accuracy evaluation method for datasets

Mean absolute error (MAE), root mean square error (RMSE), Pearson's correlation coefficient (R), and coefficient of determination (R²) are used to evaluate the accuracies of AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, ERA5-Land SWE, multisource data-averaged SWE, and the RRM SWE product. The specific equations of accuracy evaluation error are described as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }\text{MAE}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\left|f_i-y_i\right|,\text{(3)}& {\displaystyle }\text{RMSE}& {\displaystyle }={\left[{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(f_i-y_i{)}^{\mathrm{2}}}{n}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}},\text{(4)}& {\displaystyle }R& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}\sum _{i=\mathrm{1}}^n\left({\displaystyle \frac{f_i-\stackrel{\mathrm{‾}}{f}}{{\mathit{\sigma }}_f}}\right)\left({\displaystyle \frac{y_i-\stackrel{\mathrm{‾}}{y}}{{\mathit{\sigma }}_y}}\right),\text{(5)}& {\displaystyle }{R}^{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left(f_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}}},\end{array}$$

where n is the number of samples in the validation dataset, f_i is the SWE dataset product, y_i is the measured SWE at the station, $\stackrel{\mathrm{‾}}{f}$ and $\stackrel{\mathrm{‾}}{y}$ are the averages of SWE products and measured SWEs, respectively, and σ_f and σ_y are the standard deviation of SWE products and measured SWEs, respectively.

To further evaluate the accuracy of the RRM SWE dataset at the spatial scale, we compared it with AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE at different altitude gradients. We also evaluated MAE, RMSE, R, and R² separately for 11 elevation intervals: <100, 100–200, 200–300, 300–400, 400–500, 500–600, 600–700, 700–800, 800–900, 900–1000, and >1000 m. In addition, we evaluated the performances of the RRM SWE product in three representative regions: Russia, Canada, and Finland.

We used the Mann–Kendall trend test (Mann, 1945; Kendall, 1990) method to evaluate the variation trend in the RRM SWE dataset from 1979 to 2019 and analyzed its reliability in terms of time series. Since the AMSR-E/AMSR2 SWE product and the GlobSnow SWE product lack SWE data for Greenland, we removed the Greenland data to maintain consistency in the spatial extent of the comparison data.

3 Results and discussion

3.1 Overall accuracy evaluation of the RRM SWE product

In this study, the accuracies of the RRM SWE, AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE were compared using test datasets from 2015 to 2018. MAE, RMSE, R, and R² were used to reflect the data quality of each SWE product. In addition, we compared the RRM SWE product with the SWE dataset obtained by the multisource data average method.

According to the verification results in Fig. 3 and Table 2, the RRM SWE data have the best overall accuracy, and the MAE, RMSE, R, and R² between the observed SWEs are 0.21, 25.37 mm, 0.89, and 0.79, respectively. The overall accuracy of the GlobSnow SWE and ERA5-Land SWE products is higher than that of other SWE products. The overall deviation of the ERA5-Land SWE products is the smallest except for the RRM SWE data, with MAE and RMSE values of 0.32 and 37.02 mm, respectively. The correlation between the ERA5-Land SWE and observed SWE is the highest except for the RRM SWE data, with R and R² values of 0.84 and 0.71, respectively. Although the overall deviation between the GlobSnow SWE dataset and the measured SWE is small, its correlation with the measured value is low. The overall deviation between the ERA5-Land SWE dataset and the measured SWE is higher than that of the GlobSnow SWE dataset, but its estimation accuracy for the high-value region of the SWE is low. In addition, the overall accuracy of the ERA-Interim SWE dataset and GLDAS SWE dataset is relatively low, but their integrities are higher than those of the GlobSnow SWE dataset and AMSR-E/AMSR2 SWE dataset in terms of temporal and spatial series. The AMSR-E/AMSR2 SWE dataset has a higher estimation accuracy for the low-value SWE region. Moreover, in the land region above 45^∘ N, most of the existing SWE data products with regard to temporal and spatial degrees are missing to various degrees. Obviously, the accuracies of the existing SWE products were uneven as no type of SWE dataset is perfect.

Figure 3 Accuracy comparison of various SWE products. Sector (a) represents the MAE, sector (b) represents the RMSE, sector (c) represents R, and sector (d) represents R². The sector axis represents the size of the error, and the color represents different SWE datasets.

Table 2 Error list for the station data and grid snow water equivalent products.

The verification results also indicate the following ranking orders.

The MAE ranking order is RRM SWE < GlobSnow SWE = ERA5-Land SWE < ERA-Interim SWE < multisource data average SWE < AMSR-E/AMSR2 SWE < GLDAS SWE.

The RMSE ranking order is RRM SWE < ERA5-Land SWE < GlobSnow SWE < ERA-Interim SWE < multisource data average SWE < AMSR-E/AMSR2 SWE < GLDAS SWE.

The R ranking order is RRM SWE > ERA5-Land SWE > GlobSnow SWE > ERA-Interim SWE > GLDAS SWE > multisource data average SWE > AMSR-E/AMSR2 SWE.

The R² ranking order is RRM SWE > ERA5-Land SWE > GlobSnow SWE > ERA-Interim SWE > GLDAS SWE > multisource data average SWE > AMSR-E/AMSR2 SWE.

Compared with the ERA-Interim SWE, AMSR-E/AMSR2 SWE, GLDAS SWE, GlobSnow SWE, ERA5-Land SWE, and multisource data average SWE, the MAE of the RRM SWE and observed SWE is reduced by 0.22, 0.28, 0.37, 0.11, 0.11, and 0.23, respectively. The RMSE of the RRM SWE and observed SWE is reduced by 21.44, 27.02, 39.88, 15.62, 11.65, and 26.63 mm, respectively. The correlation coefficients of the RRM SWE and observed SWE are improved by 0.20, 0.42, 0.37, 0.19, 0.05, and 0.38, respectively. The coefficient of determination of the RRM SWE and observed SWE is improved by 0.31, 0.57, 0.52, 0.30, 0.08, and 0.53, respectively. Although the multisource data average method can improve the accuracy of SWE products to some extent (better than AMSR-E/AMSR2 SWE and GLDAS SWE), the improvement of this method is still very limited. The RRM SWE product has a significant advantage over the multisource data average method, and its accuracy is much higher than that of the simple multisource data average method (Table 2). Based on the above verification results, the accuracy of the RRM SWE is significantly improved; the RRM SWE dataset has higher accuracy than that of any single grid SWE dataset, and it also fills the gap in the original SWE data in terms of spatial and temporal resolutions.

Based on the kernel density estimation method, we analyzed the density distribution of different SWE datasets (Fig. 4). The results show that the RRM SWE dataset is closer to the 1:1 line and has the highest accuracy. The RRM SWE dataset is particularly accurate for SWE estimation in the low-value region, and the test data are concentrated near the 1:1 line in the high-density region (kernel density estimation >0.00015) (Fig. 4). In contrast, the high-density regions of the GLDAS SWE dataset, ERA-Interim SWE dataset, and AMSR-E/AMSR2 SWE dataset deviate significantly from the 1:1 line, resulting in poor accuracy. The AMSR-E/AMSR2 SWE, GLDAS SWE, and GlobSnow SWE are underestimated relative to the SWE measured at the site, among which GLDAS SWE underestimated the observed SWE the most seriously, while ERA5-Land SWE overestimated the observed SWE. Although the accuracies of GlobSnow SWE and ERA5-Land SWE are relatively high, their dispersion degrees are large (the kernel density estimation for most test data is less than 0.0001). Overall, the RRM SWE data have a higher overall estimation accuracy, especially for the low-value area of SWE. For an SWE above 400 mm, the MAE and RMSE of the RRM SWE product and the measured SWE are 0.35 and 43.57 mm, respectively. The estimation accuracy of the RRM SWE product for the high-value range of SWE (SWE>400 mm) is lower than that for the low-value range of SWE (SWE<400 mm) (Fig. 4). The main reason for this is that the training accuracy of the RRM model for the high-value range of SWE is affected by the small number of stations that observe the high-value range of SWE.

Figure 4 Error verification density diagram (a total of 38 807 sample points were used for verification). The color bar represents the value of kernel density estimation. The closer the high-density area is to the 1:1 line, the higher the verification accuracy of the dataset is at most of the measuring stations.

However, in this study, there are still some uncertainties in the ridge regression machine learning algorithm that integrates SWE products. First, this model is strongly dependent on on-site observational data, and the fusion precision of SWE is poor in some areas with sparse observational stations. The fusion accuracy of SWE products will be affected to a certain extent without considering the prior snow cover information. The RRM SWE product is still underestimated in cases of high SWE. Then, in addition to the DEM, meteorological elements, Normalized Difference Vegetation Index (NDVI), land type, and other factors will affect the SWE estimation. Unfortunately, our current RRM presented here does not consider these factors as predictors, which is a limitation of the current RRM SWE product. Finally, in complex terrain with an elevation interval >1000 m, the RRM SWE product performed poorly, with an RMSE of 31.14 mm (Fig. 5), and the integration of SWE products remains challenging (Mortimer et al., 2020).

Figure 5 Comparison of the error between the RRM SWE and AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE at different altitudes (the abscissa represents the altitude gradient, and the ordinate represents different SWE datasets). The color bar indicates the error in each SWE dataset. The closer to red the color is, the higher the accuracy is. MAE: mean absolute error; RMSE: root mean square error; R: Pearson's correlation coefficient; R²: coefficient of determination.

3.2 Accuracy evaluation of the RRM SWE product at different altitudes and regions

The accuracy of each SWE product is not absolute at different altitude gradients based on evaluations of the AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE product accuracies (Fig. 5). The accuracy of a single SWE product is different from its overall accuracy. We consider the influence of altitude in the algorithm and make full use of the accuracy advantage of each SWE data for different altitude gradients.

The above verification results show that the MAE, RMSE, R, and R² between the RRM SWE product and measured SWE perform well at altitude gradients of <100, 100–200, 200–300, 300–400, 400–500, 500–600, 600–700, 700–800, 800–900, 900–1000, and >1000 m (Fig. 5). Overall, the RRM SWE product has the highest accuracy in the elevation intervals of <100, 100–200, 200–300, 400–500, 500–600, 600–700, 700–800, 800–900, and >1000 m. The RRM SWE product itself has the best performance in the elevation interval <100 m. The ERA5-Land product has the best performance in the elevation interval 300–400 m. The GlobSnow product has the best performance in the elevation interval 900–1000 m.

The RRM SWE product has good performance in different regions, and its RMSEs in Russia, Canada, and Finland are 26.39, 29.31, and 25.29 mm, respectively; additionally, the performance of the RRM SWE product in different regions is basically similar (Table 3). The RRM SWE product performs well not only at different altitudes but also in different regions, and it has good stability.

Table 3 Error list for the station data and RRM SWE product in different regions.

3.3 Comparison of spatial distribution patterns between the RRM SWE product and traditional SWE products

A comparison of the spatially distributed annual average SWE distributions is made between the RRM SWE and AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE in 2014, 2015, 2016, and 2017, and their spatial distribution patterns are shown in Fig. 6.

Figure 6 Comparison of the spatial distribution characteristics between the RRM SWE and AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE (the four columns of images represent the comparison results in 2014, 2015, 2016, and 2017).

Overall, the RRM SWE dataset, AMSR-E/AMSR2 SWE dataset, ERA-Interim SWE dataset, GLDAS SWE dataset, GlobSnow SWE dataset, and ERA5-Land SWE dataset have similar spatial distribution patterns in the land region above 45^∘ N, showing a trend of lower SWE at low latitudes and higher SWE at high latitudes. The AMSR-E/AMSR2 SWE dataset covers a limited extent in the land region above 45^∘ N, many data points are missing, and low SWE values exist at low latitudes. In northern Siberia, the ERA-Interim SWE product has a higher SWE, and there are many abnormal, extreme SWE values (SWE>500 mm) in this dataset. In low-latitude regions, such as Alaska, northern Siberia, and the easternmost region of Russia, the SWE of GLDAS SWE products is significantly lower. The GlobSnow SWE product lacks SWE data for Greenland, and this dataset has low SWEs in the regions of Baffin Island, Koryak Mountains, Kamchatka Peninsula, and Alaska. The ERA5-Land SWE products have low SWEs in northeastern Russia, Scandinavia, and northeastern Canada. The RRM SWE dataset is more reasonable for estimating the spatial distribution of SWE in the land region above 45^∘ N, and the data integrity is higher. Moreover, based on the new machine learning algorithm, a variety of SWE data products in different time series are fused, which makes the RRM SWE dataset completely temporally and spatially continuous.

The relative difference between the RRM SWE data and GLDAS SWE data is the highest, and the relative difference is greater than 80 % in most low altitude regions (Fig. 7). The relative difference between the RRM SWE data and the GlobSnow SWE data is relatively small overall, especially in most high-latitude areas where the relative difference is less than 10 % (Fig. 7). Overall, the annual average relative differences in the RRM SWE data and AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE are 37 %, 41 %, 54 %, 25 %, and 29 %, respectively (Fig. 7). Previous studies have shown that the accuracy of the SWE in the Northern Hemisphere estimated by GlobSnow SWE data is higher (Pulliainen et al., 2020), while the spatial distribution pattern of the RRM SWE data is close to the estimation result of GlobSnow SWE. In addition, the single point verification results based on the measured SWE data of meteorological stations in Sect. 3.1 show that the RRM SWE dataset has higher accuracy than the GlobSnow SWE dataset. The RRM SWE dataset has good accuracy.

Figure 7 Temporal and spatial distributions of relative differences (RD%) between the RRM SWE and AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE. Lower-right panel: comparison of annual average relative differences between the RRM SWE and AMSR2 SWE (A), ERA-Interim SWE (B), GLDAS SWE (C), GlobSnow SWE (D), and ERA5-Land SWE (E).

3.4 Comparison of the annual variation tendencies of AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE and the RRM SWE in the land region above 45^∘ N

Based on the Mann–Kendall trend test, we analyzed the changing trend in the region-wide annual average SWE of the AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, ERA5-Land SWE, and RRM SWE in the land region above 45^∘ N from 1979 to 2019.

Based on the Mann–Kendall trend test (see Fig. 8 and Table 4), from 1979 to 2019, the test value of the ERA-Interim region-wide annual average SWE is 1.08, and there is no significant change trend under the significance test level of 0.05. The test value of the GLDAS region-wide annual average SWE was 4.95 and showed a significant increasing trend at the significance test level of 0.05. The test values of the AMSR-E/AMSR2 annual average SWE, GlobSnow annual average SWE, ERA5-Land annual average SWE, and RRM annual average SWE are −3.26, −2.54, −3.43, and −3.00, respectively, and these four SWEs showed a significant decreasing trend at the significance test level of 0.05. Based on the analysis of the RRM SWE product, between 1979 and 2019, the region-wide annual average SWE in the land region above 45^∘ N decreased by 15.1 %. In the Northern Hemisphere, spring snow cover extent has decreased significantly, according to the Fifth Assessment Report (AR5) of the IPCC. Between 1967 and 2010, the spring snow cover extent decreased by an average of 1.6 % per decade, while the June snow cover extent decreased by 11.7 % per decade (Stocker, 2014). Most studies have shown that the annual variation tendency of snow depth and snow cover extent showed a significant decreasing trend in the Northern Hemisphere (Brutel-Vuilmet et al., 2013), which is consistent with the annual variation tendency of the RRM SWE dataset. This dataset can reflect the characteristics of snow cover change in the land region above 45^∘ N in light of climate change and can be used as the driving data for climate models to support climate-change-related research. In addition, this dataset is expected to provide a snow data basis for the study of “Arctic amplification”.

Figure 8 Annual variation tendency of the AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, ERA5-Land SWE, and RRM SWE products from 1979 to 2019 (the dotted line is the trend line calculated based on the Mann–Kendall method).

Table 4 Results of the Mann–Kendall trend test performed for various snow water equivalent products from 1979 to 2019.

^∗ Significance level alpha=0.05.

4 Data availability

The RRM SWE product is available for free download from “A Big Earth Data Platform for Three Poles” (https://doi.org/10.11888/Snow.tpdc.271556, Li et al., 2021). The temporal resolution of the RRM SWE product is daily, and the spatial resolution is 10 km. It spans latitudes of 45–90^∘ N and longitudes of 180^∘ W–180^∘ E. A brief summary and data description document (including data details, spatial range, and usage method) are also provided.

5 Conclusions

In this study, we propose a method to fuse multisource SWE data by a ridge regression model based on machine learning. A new method was utilized to prepare a set of spatiotemporally seamless SWE datasets of the RRM SWE, combined with the original AMSR-E/AMSR2 SWE, ERA-Interim SWE, GLDAS SWE, GlobSnow SWE, and ERA5-Land SWE datasets. In the RRM SWE dataset, the time series of the data is 1979–2019, the temporal resolution is daily, the spatial resolution is 10 km, and the spatial range is the land region above 45^∘ N.

The RRM SWE data product has the best accuracy, especially for the estimation of low SWE. The accuracy ranking of the SWE dataset verified by the test dataset is described as follows: RRM SWE > ERA5-Land SWE > GlobSnow SWE > ERA-Interim SWE > multisource data average SWE > AMSR-E/AMSR2 SWE > GLDAS SWE. The accuracy of the RRM SWE dataset is higher than that of the existing SWE products at most elevation intervals. The RRM SWE product has good performance and stability in different regions. Moreover, the RRM SWE dataset spatiotemporally fills in the missing data of the original SWE dataset.

Compared with traditional fusion methods, machine learning methods have a strong advantage. We find that the simple machine learning algorithm has not only high efficiency but also good accuracy in the preparation of SWE products on a global scale. Without losing the advantages of existing SWE products, this method can also make full use of station observational data to integrate the advantages of various SWE products. The model training process does not rely too much on a specific sample, and this model has a strong generalization ability. In addition, the influence of altitude on the preparation scheme is considered in detail in the model. Compared with the SWE dataset prepared by the traditional method, the spatial resolution is only 25 km, while this new method obtains an SWE dataset with a higher spatial resolution of 10 km.

We propose that the RRM SWE dataset preparation scheme has good continuity and can prepare real-time and high-quality SWE datasets in the land region above 45^∘ N. In addition, the new method proposed in this paper has the advantages of simplicity and high precision in preparing large-scale SWE datasets and can be easily extended to the preparation of other snow datasets. This dataset is an important supplement to the land region above the 45^∘ N SWE database and is expected to provide data support for Arctic cryosphere studies and global climate change studies.