2.1 CMIP6 model data

The data of this study consist of CMIP6 climate model data (Table 1) and observational and reanalysis data (Table 2). For CMIP6, we used monthly mean data from those models for which either historical or esm-hist simulations were available for download in August 2021. The historical and esm-hist simulations extend from 1850 to 2014. In historical simulations, the CO₂ concentrations are prescribed, whereas in esm-hist simulations, the models calculate the atmospheric CO₂ concentration interactively based on prescribed CO₂ fluxes (Eyring et al., 2016). In this study, we used three variables from CMIP6 models: SWE (variable “snw”, unit kg m⁻²), surface air temperature (“tas”, unit kelvin), and precipitation (“pr”, unit kg m⁻² s⁻¹). Altogether, there are 17 high-resolution models (100 km) and 21 low-resolution models (250–500 km) in Table 1. We have evaluated the SWE sum over the entire study area for every model but performed the subsequent detailed analysis only for the high-resolution (100 km) models. We decided to leave the models with coarser resolution out of the analysis, as coarser resolution would differ too much from the resolution of the observational datasets, making the comparison more problematic. The number of ensemble members available for the chosen models varies between 1 and 16. For simplicity, we only consider the first member of each model ensemble (r1i1p1f1) in this study. A brief analysis showed that the differences between different ensemble members for the same model were generally smaller compared to intermodel differences. Figure S1 in the Supplement shows all realizations of three different models (CESM2, MPI-ESM1-2-HR, and EC-Earth3), which were chosen in this figure due to a high number of realizations. The figure illustrates that internal variability of each model is smaller than the intermodel variability.

Table 1CMIP6 models used in this study.

Download Print Version | Download XLSX

Table 2Observational and reanalysis datasets used in this study.

Download Print Version | Download XLSX

2.2 Observational and reanalysis data

For SWE reference data, we have used ESA CCI Snow “SnowCCI” (European Space Agency Climate Change Initiative, Snow) data for the non-mountainous regions and MERRA-2 (the Modern-Era Retrospective analysis for Research and Applications, Version 2; Gelaro et al., 2017; GMAO, 2015a), Brown (Brown et al., 2003), and Crocus v7 (Brun et al., 2013) data for the mountainous regions. SnowCCI data, which is the same product as the GlobSnow v3 SWE product (except provided in geographical latitude–longitude grid for easier comparison with climate model data), are based on satellite passive microwave radiometer data and in situ data (Luojus et al., 2021; Pulliainen et al., 2020). The SnowCCI algorithm combines microwave brightness temperature (Tb) data, observed by satellite instruments, with ground-based snow depth measurements from the global network of synoptic weather stations (Luojus et al., 2021). The SWE estimation algorithm is based on the difference in Tb between two frequencies (37 and 19 GHz). The ground beneath the snowpack emits microwaves, which propagate through the snowpack, being partially absorbed during the process. The low-frequency and high-frequency signals attenuate differently as they propagate through the snowpack, which makes the difference in Tb a good indicator for estimating SWE (Cagnati et al., 2004). The attenuation is affected by snow depth, snow grain size, and snow density. The high-frequency signal attenuates more than the low-frequency signal when it propagates through the snowpack, especially for a deep, dense, and large-grained snowpack. Thus, a large difference between high- and low-frequency signals indicates a high SWE (Kelly et al., 2003). The SnowCCI approach combines Tb differences with in situ snow depth observations, which considerably improves SWE estimation relative to a satellite-only retrieval (Pulliainen, 2006; Takala et al., 2011).

A recent bias-correction method combines the original SnowCCI data with extensive ground-based snow course SWE measurements, which significantly reduces the uncertainty of NH SWE estimation (Pulliainen et al., 2020). The method decreases the uncertainty of hemisphere-mean SWE estimation from 33 % to 7.4 %. The bias-corrected SnowCCI data are mapped to a 25 km EASE-Grid and are available from 1979. Mountainous regions, glaciers, and ice sheets are excluded from the data. The original SnowCCI data are available around the year, while bias-corrected SnowCCI data are only available from February to May. Despite limitations in its temporal coverage, we have used the bias-corrected data in this study. We chose to do this as the bias-correction method significantly reduced the uncertainty and made the observational data more accurate, which, in turn, makes the comparison with the models far more meaningful.

As the bias-corrected SnowCCI data are only available for non-mountainous grid cells, we have used an average of the MERRA-2, Brown, and Crocus v7 SWE products for the mountainous regions. MERRA-2 is a NASA (National Aeronautics and Space Administration) atmospheric reanalysis, and it is available from year 1980. The spatial resolution of the data is $\mathrm{0.625}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ (Gelaro et al., 2017). From MERRA-2, we have used the SWE product (Gelaro et al., 2017; GMAO, 2015a). Brown SWE product, in turn, uses a simple snow scheme driven by ERA-Interim reanalysis (Brown et al., 2003). For the third SWE dataset, we have used Crocus version 7 product, which is a physical snow model driven by ERA-Interim reanalysis (Brun et al., 2013). Both MERRA-2 and Crocus v7 tend to slightly overestimate SWE under 150 kg m⁻² and underestimate SWE over 150 kg m⁻² (Mortimer et al., 2020).

For precipitation (P) and temperature (T) reference data, we used GPCC (Global Precipitation Climatology Centre) version 2018 precipitation data (Schneider et al., 2018) and MERRA-2 temperature data (Gelaro et al., 2017; GMAO, 2015b). GPCC is a monthly precipitation product based on data from rain gauge stations, and the data are available on a 0.5^∘ global grid from 1891 to the present (Schneider et al., 2018). The product agrees well with other precipitation products (Behrangi et al., 2016; Sun et al., 2018). All precipitation data are presented here in units of kilograms per square meter per month (kg m⁻² per month), which is equivalent to millimeters per month (mm per month).

For T reference data, we have used the monthly mean 2 m air temperature product, which agrees well with observations in the Arctic (Gelaro et al., 2017; Simmons et al., 2017), and the mean values show very small biases (Bosilovich et al., 2015). MERRA-2 daily temperature tends to have a cool daytime bias and a warm nighttime bias (Bosilovich et al., 2015; Draper et al., 2018). However, this is not a major issue for our study because we use the monthly mean product. In addition, MERRA-2 seems to underestimate global warming trends in the last years of our study period (Gelaro et al., 2017; Simmons et al., 2017). Additionally, we used the ESA CCI Land Cover time series v2.0.7 (ESA, 2017) to study the effect of forest cover on the results.

4.1 SWE in winter

The mean SWE bias in February shows large spatial variability within each model but also varies much between the models (Figs. 3 and S3). The areal-mean model bias varies from about −8 kg m⁻² to 40 kg m⁻². However, the largest negative and positive biases are well concentrated in the same areas in all models. Overall, the models tend to overestimate the SWE in the northern parts of the study area but also in southern Siberia. The negative biases, in turn, occur mostly in the south and especially on the coastal areas of North America. In some models, there are negative biases also in the middle parts of Eurasia. Large biases also occur over Rocky Mountains in every model, but the sign of the bias varies. The CESM2 and NorESM2-MM models show the largest overestimations. For both models, the bias is very high in large regions in northern parts of North America and Eurasia. In these areas the relative bias is typically 150 %–200 %. The BCC-CSM2-MR, E3SM-1-0, and SAM0-UNICON models also show large positive biases, which are concentrated in the same areas; although, the biases are clearly smaller than for the CESM2 and NorESM2-MM models. In other models, the areal-mean biases are closer to 0 kg m⁻²; however, regional differences exist. Overall, the GFDL models are the most consistent with the SWE reference data during February.

All models overestimate NH extratropical precipitation from November to January (Figs. 4 and S4). The largest overestimations are mainly in southern regions and in coastal areas. There are small areas where underestimation occurs, especially in Eurasia around 90^∘ E. In every model, there are large regions where P_cum bias could logically explain the SWE bias (marked with dots), i.e., areas where models overestimate or underestimate both SWE and P_cum. These regions are mostly in the northern parts of the study area, whereas in the south there are more areas where SWE and P_cum discrepancies are more often not consistent with each other.

Figure 3The mean model bias in SWE (ΔSWE, model minus observation) in snow-covered areas in February 1982–2014.

For model T biases, the values are mostly positive; however, regional and intermodel variability exists (Figs. 5 and S5). We show ΔT (ΔT_cum divided by 3) instead of ΔT_cum so that the values are more intuitive and easier to interpret. The CESM2 and E3SM-1-0 models generally simulate too warm temperatures, and the largest positive biases are in the northern parts of the study area. The GFDL, MPI-ESM1-2-HR, and MRI-ESM2-0 models simulate too warm temperatures in the north, while the SAM0-UNICON and TaiESM1 models, in turn, simulate too cold temperatures in the north. The BCC-CSM2-MR and EC-Earth3 models are the most consistent with the MERRA-2 data. In every model, there are large areas where the signs of biases for T and SWE are opposite (marked with dots), indicating that biases in T might explain biases in SWE. However, in these areas, the biases are mainly quite small.

The contributions to SWE biases (ΔSWE) due to precipitation biases (P_C), temperature biases (T_C), and other factors (residual R) are quantified using the regression Eqs. (1), (3), and (4). To summarize their relative importance, Figs. 6 and S6 show the areal means of the absolute values of ΔSWE, P_C, T_C, and R. The mean ΔSWE varies from under 30 kg m⁻² in the GFDL-CM4 model to around 50 kg m⁻² in the CESM2 and NorESM2-MM models. In all models, the contribution of P on ΔSWE is clearly larger than the contribution of T. However, the residual is also typically large, indicating that P and T cannot explain the SWE biases alone. This implies that observational uncertainty or model structural factors, such as parameterizations related to the surface energy budget or hydrology, play a considerable part in the observed SWE biases. We studied the terms' statistical significance by using the t test and found that all the terms (P_C, T_C, and R) are significantly different from zero at the 99 % confidence level. The variability in these parameters between the decadal subperiods and the full three-decade analysis period was slight, indicating consistent behavior across time in both models and observations.

Figures 7 and S7 show the spatial distribution of the contributions of P and T and the residual for each model for the entire study period 1982–2014. The spatial distributions of R², β_P, and β_T are shown in the Supplement (Fig. S8). Also, the contributions of P and T and the residual calculated for the shorter time periods are in the Supplement (Figs. S9–S11). Figure 7 shows that, overall, the contribution of P is larger than the contribution of T, as Fig. 6 already indicated. P contributes to ΔSWE especially over northern and coastal regions, with fairly similar patterns for all models considered. The regression coefficient β_P also shows large values (β_P≈1), especially in the northern and middle parts of both continents (Fig. S8), with relatively small intermodel variations. This is consistent with the expectation that in cold regions an increase in precipitation translates into a similar increase in SWE.

Figure 4The mean model bias in P_cum (ΔP_cum, model minus observation) in snow-covered areas from November to January 1982–2014. The dots indicate areas where the models either overestimate both SWE and P or underestimate both SWE and P.

The contribution of T is mostly very weak (Fig. 7); however, for some models, T shows stronger contribution especially over western parts of Eurasia and over northeastern Canada. The regression coefficient β_T is mostly negative or very close to zero (Fig. S8). The negative correlation is strongest in Europe and southern parts of North America and Eurasia. In these regions, the temperature is close to 0^∘ C in winter, which makes temperature an important driver of the SWE. Northern Canada and Siberia, in turn, show areas with positive correlation, meaning that warmer temperatures cause higher SWE. Studies have shown that when the temperature rises due to climate change, the winter precipitation will also increase (Brown and Mote, 2009; Park et al., 2012; Räisänen, 2008). In the south, warming temperatures will shift winter precipitation from snow to rain. In the north, in turn, temperature will stay below 0^∘ C despite the warming, which will lead to an increase in snowfall in the coldest regions of NH and, therefore, to an increase in SWE. This phenomenon is most likely seen here as well; warmer temperatures in the models will increase winter precipitation, resulting in too high SWE in the models. In fact, this exposes a limitation of the regression Eq. (1): it treats ΔT_cum and ΔP_cum as independent variables, which is not fully realistic. When these variables are correlated, their contributions to ΔSWE cannot be fully separated.

The residual shows large spatial and intermodel variability (Fig. 7). Especially for the CESM2 and NorESM2-MM models, the residual shows very large positive values. These large positive residuals are mainly concentrated in the same areas where the models clearly overestimate SWE (Fig. 3). This indicates that, for these models, the large SWE biases in these areas are mainly caused by some other factors than P or T. For other models, the residual shows both positive and negative values across the study area. The possible factors causing the large residual term are discussed in more detail in Sect. 5.2.

Figure 5The mean model bias in T (ΔT, model minus observation) in snow-covered areas from November to January 1982–2014. We show ΔT (ΔT_cum divided by 3) instead of ΔT_cum so that the values are more intuitive and easier to interpret. The dots indicate areas with either cold bias and positive SWE bias or warm bias and negative SWE bias, i.e., the areas where T bias could logically explain the SWE bias.

4.2 Monthly SWE change in spring

Figures 8 and S13 show the mean bias in monthly SWE change (ΔSWE_change) in spring in the snow-covered areas for the whole study period 1982–2014. A figure showing the SWE change rate of multi-model mean, SWE reference data, and bias in multi-model mean is in the Supplement (Fig. S12). Positive ΔSWE_change means that there is less snowmelt in the models compared to the SWE reference data, and negative ΔSWE_change indicates that snow melts faster in models, respectively. The areal-mean ΔSWE_change is mainly negative in every model, which means that snow melts generally faster in the models compared to the SWE reference data. However, there is a large spatial variability in every model, and the intermodel variability is also large. In the CESM2 and NorESM2-MM models, three areas show distinctly positive ΔSWE_change: northeastern Canada, northern Siberia, and eastern Siberia. In all these areas, the SWE bias in February (Fig. 3) was already clearly positive, meaning that these models greatly overestimate SWE in these areas also in spring. The EC-Earth3 model shows clear positive ΔSWE_change in northeastern Canada but also in Eurasia. The area with positive biases in Eurasia is very extensive and differs notably from the other models. The EC-Earth3 model has been found to drastically overestimate the snow cover extent as well (Mudryk et al., 2020). The SAM0-UNICON and TaiESM1 models also show positive values in northern Siberia. The GFDL and MPI-ESM1-2-HR models, in turn, show large areas with negative ΔSWE_change in northern Canada and in eastern parts of Siberia, which differs from the other models. Overall, the model biases in the SWE change rate in spring (Fig. 8) show larger intermodel variations than the corresponding biases in SWE in February (Fig. 3).

As in winter (Fig. 4), models on average overestimate precipitation in spring as well (Figs. 9 and S14). The largest overestimations occur mainly in southern regions and in coastal areas. The regions with mutual biases in P and SWE_change (marked with dots) show large intermodel variability, and they are less extensive than in winter (Fig. 4). This indicates that precipitation is not as important factor in spring as it is in winter.

The mean T bias in spring shows a large spatial variability within each model, but also the variability between models is very clear (Figs. 10 and S15). The BCC-CSM2-MR, E3SM-1-0, MPI-ESM1-2-HR, and MRI-ESM2-0 models show a warm bias in the northern parts of the study area, whereas the SAM0-UNICON and TaiESM1 models show a cold bias in the same area. The EC-Earth3 model, in turn, has a cold bias in eastern Eurasia near 60^∘ N, which clearly differs from the other models. These areas also exhibited less snowmelt than the SWE reference data, which indicates that the bias in snow melt rate can be caused by the bias in T. The sizes and locations of the areas with mutually consistent biases in T and SWE_change (marked with dots) vary greatly between models. Especially in the GFDL-CM4 model, the size of these areas is small, while in most of the models, the dots cover the majority of the study area. Furthermore, in those regions where the biases in T and SWE_change are consistent in spring, the cold or warm temperature biases are typically relatively large, when compared with the corresponding biases in winter (Fig. 5). This indicates that biases in T are a more important driver of biases in SWE in spring than in winter.

Figure 6The areal means of the absolute values of ΔSWE, P_C, T_C, and residual R calculated for the entire study period 1982–2014 (left column, shaded with grey) and for three shorter time periods (1982–1991, 1992–2001, and 2002–2014) for each model in winter. The size of the square indicates the absolute value of ΔSWE of that time period and model, and the color of the square indicates the absolute value of P_C, T_C, and R.

Download

Figures 11 and S16 summarize the areal means of the absolute values of ΔSWE_change, the contribution of P (P_C), the contribution of T (T_C), and the contribution of other factors (the residual R). The contributions of P and T are quite similar in magnitude in almost all the models. None of the variables shows a large dependence on the period considered. The EC-Earth3 model stands out from the other models, as T_C is larger than in the other models. The large contribution of T in ΔSWE_change was already seen in Figs. 8 and 10, as the spatial distribution of the biases were similar. This also suggests that the overestimated snow cover extent in the EC-Earth3 model (Mudryk et al., 2020) might be at least partly due to the cold bias. In all models, the residual is larger than P_C or T_C. This indicates that, overall, the biases in snow melt rate in spring are dominated by other factors than T or P. These factors will be further discussed in Sect. 5.2. We also studied the terms' statistical significance by using the t test. Even though the contributions of T and P are considerably smaller than the residual term, all the terms (P_C, T_C, and R) are significantly different from zero at the 99 % confidence level.

Figures 12 and S17 show the spatial distribution of the contributions of P and T and the residual for each model for the entire study period 1982–2014. The spatial distributions of R², β_P, and β_T are displayed in the Supplement (Fig. S18). Also, the contributions of P and T and residual calculated for the shorter time periods are in the Supplement (Figs. S19–S21). P contributes to ΔSWE_change mostly in Alaska and northern Siberia, but the effect varies between models. Furthermore, even though T showed clear warm and cold biases in many areas (Fig. 10), the contribution of T is mostly quite weak, because of the small regression coefficient β_T (Fig. S6). However, exceptions exist; in particular, the EC-Earth3 and EC-Earth3-Veg models stand out, as in Eurasia there is a large area where a negative bias in T (Fig. 10) contributes substantially to a positive bias in SWE_change (Fig. 8). Also, for the CESM2 and GFDL-ESM4 models, T shows a stronger contribution over northern parts of North America. Overall, however, the contributions of both P and T are small compared to the residual R, which is consistent with Fig. 11. This indicates that other factors than T or P are the dominant drivers for the SWE_change discrepancies. The residual is mostly negative in all the models, which means that snow would melt too fast in the models if T and P were simulated correctly. This issue is discussed further in Sect. 5.2.

5.1 The contribution of P and T to the SWE bias

In winter, the models mostly overestimate SWE (Fig. 3), but spatial and intermodel variability exists. The overestimations are mostly concentrated in the same areas in all models, but the magnitude varies greatly between the models. The models also overestimate precipitation in winter (Fig. 4), suggesting that overestimated SWE is caused by the overestimated P. The regression coefficient β_P also shows very clear correlation between ΔSWE and ΔP_cum (Fig. S8). Therefore, P clearly contributes to the SWE bias, whereas the contribution of T is substantially smaller (Fig. 7). This result is consistent with the expectation that precipitation is needed to initiate the snow cover. In other words, even too cold temperatures cannot cause too high SWE without precipitation. Also, studies have shown that both the P bias (Zamani et al., 2020) and the SWE bias (Mudryk et al., 2020) have decreased from CMIP5 to CMIP6. This indicates that there is a connection between these variables, which is consistent with the results from this study.

Figure 10The mean model bias in T (ΔT, model minus observation) in snow-covered areas in spring 1982–2014. The dots indicate areas with either cold bias and positive SWE bias or warm bias and negative SWE bias, i.e., the areas where T bias could logically explain the SWE bias.

The link between ΔT_cum and ΔSWE is strongest in the warmest regions of the study area (Figs. 7 and S8). Especially in the coldest regions, where T is well below 0^∘ C, variations in T do not significantly affect the amount of snow on the ground. In regions where T is closer to 0^∘ C in winter, T plays a more significant role and has a greater impact on SWE. This physically reasonable behavior suggests that climate models might be able to simulate SWE trends in the warming climate correctly, even if SWE itself is not reproduced accurately.

In spring, ΔSWE_change and ΔT show quite similar patterns in many models (Figs. 8 and 10), which indicates that biases in T affect the biases in SWE_change. This result is to be expected because the increase in T is the main factor that causes snow to melt in spring. A relationship between temperature and snow cover in spring has also been observed in CMIP5 models (Brutel-Vuilmet et al., 2013; Mudryk et al., 2017). The CMIP5 models have been found to underestimate the observed trend towards a reduced snow cover extent due to an underestimation of the spring warming trend (Brutel-Vuilmet et al., 2013).

Even though ΔSWE_change and ΔT are quite consistent with each other, the contribution of T is not very strong (Fig. 12), because the regression coefficient β_T is small (Fig. S18). The only model showing clear contribution of T is EC-Earth3, as there is a large area in eastern Eurasia where positive values of ΔSWE_change co-occur with a cold bias. The spatial distribution of the bias in SWE_change in the EC-Earth3 model differs considerably from other models (Fig. 8). The model has been found to drastically overestimate snow cover extent though not so much the SWE (Mudryk et. al. 2020). This suggests that while the positive bias in SWE_change (i.e., underestimated snowmelt rate in spring) is related to the cold bias, the cold bias might be maintained by a too high surface albedo resulting from overestimated snow cover extent.

Figure 11The areal means of absolute values of ΔSWE_change, P_C, T_C, and residual R calculated for the entire study period 1982–2014 (left column, shaded with grey) and for three shorter time periods (1982–1991, 1992–2001, and 2002–2014) for each model in spring. The size of the square indicates the mean absolute value of ΔSWE_change of that time period and model, and the color of the square indicates the mean absolute value of P_C, T_C, and R.

Download

Several factors can weaken the regression coefficient β_T. The analysis shows that in spring the relationship between ΔT and ΔSWE_change is strongest in western parts of Eurasia and in southern parts of North America (Fig. S18). These are regions with the earliest snowmelt onset. In these areas, T is the dominant factor causing snow to melt throughout the spring season. In the northernmost parts of the study area, the melt season begins later, so that early spring still belongs to the snow accumulation season. As a result, P may still be the dominant factor influencing ΔSWE_change in early spring, and T will become a more significant factor later in the spring. This can reduce the correlation and affect R², β_P, and β_T in the northern parts of the study area. Additionally, if SWE in winter is simulated incorrectly, this can affect the melt rates in spring, as there can be too little or too much snow that can melt.

5.2 The residual term

In both winter and spring, the residual term R of the regressions is also significant (Figs. 6, 7, 11, and 12). This means that the model biases in SWE in winter and SWE_change in spring cannot be entirely explained by the biases in P and T, but other factors also contribute to these biases. These factors may include inaccuracy in the model parameterizations related to snow and surface energy budget but also inaccuracy in the observational datasets.

The residual term R is particularly pronounced for the SWE change rate in spring, when it is typically larger than either the contribution of P or T. Interestingly, the residual is mostly negative (Fig. 12). The negative residual means that if P and T were simulated correctly in the models, snow would melt too fast in spring. While understanding the origins of this bias would be worth a separate study, a previous study with ECHAM5 (Räisänen et al., 2014) is of interest here. ECHAM5 is a predecessor of the atmospheric part of MPI-ESM1-2-HR, for which the residual R in Fig. 12 is especially strongly negative. This is consistent with the finding that in ECHAM5 snow melted generally too fast in spring, despite a cold bias in T (Räisänen et al., 2014). A major factor for this was that T was not calculated separately for snow-covered and snow-free parts of the grid cell. Because of that, T was not able to rise above 0^∘ C if there was snow left in the grid cell, and, therefore, a too large fraction of the available energy was used in melting the snow (Räisänen et al., 2014).

Figure 12Spatial distribution of the P contribution, the T contribution, and the residual for each model in spring 1982–2014.

The parameterization of the surface albedo is another factor that influences the snowmelt rate. In spring, the surface albedo feedback is strong, and if the albedo is not correctly simulated in the models, it can cause large biases in SWE in spring (Thackeray et al., 2018). CMIP5 models have had challenges in simulating correctly the surface albedo on snow-covered land (Thackeray et al., 2015). The intermodel variability in surface albedo feedback is large in CMIP5 models (Qu and Hall, 2014; Thackeray et al., 2018), and the large variability has been found to persist in CMIP6 models (Thackeray et al., 2021), suggesting that there are still large uncertainties in simulating the surface albedo. Specifically, regarding the discussion on ECHAM5 and MPI-ESM1-2-HR above, the only region in which MPI-ESM1-2-HR displays a positive residual in Fig. 12 is southeastern Siberia. In this very region, ECHAM5 featured delayed snowmelt, related to overly high albedo in the presence of vegetation over snow (Räisänen et al., 2014). While the specific mechanisms leading to too fast snowmelt might differ in different models, the example of ECHAM5 highlights the importance of the treatment of the surface energy budget in the presence of snow.

In addition to issues related to the models, uncertainties in the observational data can also affect the residual values. With the bias-correction method, SWE data are more accurate than before, but the uncertainty in hemisphere-mean values is still 7.4 %. In mountainous regions, the complex terrain poses a challenge in the SWE estimation, but averaging over multiple products can improve the accuracy (Mortimer et al., 2020). There are also uncertainties associated with the GPCC and MERRA-2 datasets that can cause errors in assessing the model biases; for example, MERRA-2 underestimates global warming trends in the last years of our study period compared to other reanalyses (Gelaro et al., 2017; Simmons et al., 2017). Snow cover in spring is especially sensitive to warming (Hernández-Henríquez et al., 2015), and, therefore, the uncertainties in MERRA-2 can affect the results especially in spring.

The uncertainties associated with the GPCC precipitation product can be divided into two categories: the systematic measuring error and stochastic sampling error (Schneider et al., 2014). The systematic error results from evaporation and from droplets and snowflakes being drifted across the gauge funnel. This error almost always causes underestimation of precipitation (Schneider et al., 2014). The sampling error, in turn, is caused by the sparse network density of the in situ stations, and the error increases when the density of the network decreases. A correction method taking into account the weather conditions has been implemented in the GPCC to improve the P estimate (Schneider et al., 2014). In general, estimating precipitation in high latitudes is a challenge. The estimates have improved, but especially regional discrepancies still exist between different precipitation products (Behrangi et al., 2016). This causes uncertainties in assessing the model bias and means that the choice of the precipitation product may influence the results.

We additionally studied the dependency between the residual term and fractional forest cover (Fig. S22), and we did this separately for the entire study area and for non-mountainous regions. In the entire study area, the spread of the residual term is quite large when the fraction of forest cover is small in both winter and spring. When looking at only the non-mountainous areas, the spread decreases compared to the entire study area, especially in spring. This illustrates that the residual term is particularly large in mountainous areas, indicating that in mountainous areas, T or P can only explain a small fraction of SWE bias, but other factors contribute SWE bias considerably. Also, in mountainous areas, the intermodel variability of the residual term is substantial. The complex terrain and large SWE variability make SWE estimates challenging. Also, the original resolution of CMIP6 models can be too coarse for accurately describing SWE in mountainous terrain; biases in snow cover and snow depth have been found to decrease considerably with downscaled regional climate models in mountainous regions (Matiu et al., 2019). Resampling the datasets to a common grid can cause uncertainties in mountainous areas due to the high topographic variability within a grid cell.

For most models, the residual term in non-mountainous areas in winter depends little on the fractional forest cover, indicating that the residual arises from other sources. For the CESM2 and NorESM2-MM models, however, the residual decreases when forest cover fraction increases, and this dependency is strikingly similar for these two models. Since CESM2 and NorESM2-MM both employ the CLM5 land surface model (Lawrence et al., 2019), this hints that the large positive SWE residual in these models originates from CLM5. Also, a unique feature within CESM2 that may cause high residual values in the coldest regions is that the model allows for a very large maximum SWE to enable the simulation of firn production over ice sheets (van Kampenhout et al., 2017). Even though our study does not cover ice sheets, this feature may cause positive SWE bias in the coldest regions of the study area (e.g., in the northernmost regions of Canada).

In spring in non-mountainous areas, the residual of SWE_change depends quasi-systematically on the forest cover fraction. While the level of the residual terms varies between the models, the values increase with increasing forest cover fraction in all models. The residual varies mostly between −25 and −12 kg m⁻² per month for open areas and between −13 and −2 kg m⁻² per month when the forest cover fraction is close to 100 %. This leaves room for at least two possible interpretations. On one hand, the larger intermodel spread of the residual over open areas suggests that problems in snow albedo parameterization (e.g., how to account for snow metamorphism, Colbeck, 1982) could be at least partly responsible for the large residual term. Due to the shading effect of the canopy (Malle et al., 2019), snow albedo treatment differences in the models should manifest more strongly over open tundra regions and less so over dense forest cover. On the other hand, the quasi-systematic dependence of the residual on the forest cover fraction could also point to some subtle issues in the SWE estimates of the SnowCCI dataset in spring.

Overall, the reasons behind the residual term are complex, and more detailed model-specific investigations for all participating CMIP6 models are beyond the scope of this study. Further analysis would be required in the future to fully understand the factors behind the residual term.