Evaluation of snow depth and snow cover over the Tibetan Plateau in global reanalyses using in situ and satellite remote sensing observations

The Tibetan Plateau (TP) region, often referred to as the Third Pole, is the world's highest plateau and exerts a considerable influence on regional and global climate. The state of the snowpack over the TP is a major research focus due to its great impact on the headwaters of a dozen major Asian rivers. While many studies have attempted to validate atmospheric reanalyses over the TP area in terms of temperature or precipitation, there have been – remarkably – no studies aimed at systematically comparing the snow depth or snow cover in global reanalyses with satellite and in situ data. Yet, snow in reanalyses provides critical surface information for forecast systems from the medium to sub-seasonal timescales.

Here, snow depth and snow cover from four recent global reanalysis products, namely the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 and ERA-Interim reanalyses, the Japanese 55-year Reanalysis (JRA-55) and the NASA Modern-Era Retrospective analysis for Research and Applications (MERRA-2), are inter-compared over the TP region. The reanalyses are evaluated against a set of 33 in situ station observations, as well as against the Interactive Multisensor Snow and Ice Mapping System (IMS) snow cover and a satellite microwave snow depth dataset. The high temporal correlation coefficient (0.78) between the IMS snow cover and the in situ observations provides confidence in the station data despite the relative paucity of in situ measurement sites and the harsh operating conditions.

While several reanalyses show a systematic overestimation of the snow depth or snow cover, the reanalyses that assimilate local in situ observations or IMS snow cover are better capable of representing the shallow, transient snowpack over the TP region. The latter point is clearly demonstrated by examining the family of reanalyses from the ECMWF, of which only the older ERA-Interim assimilated IMS snow cover at high altitudes, while ERA5 did not consider IMS snow cover for high altitudes. We further tested the sensitivity of the ERA5-Land model in offline experiments, assessing the impact of blown snow sublimation, snow cover to snow depth conversion and, more importantly, excessive snowfall. These results suggest that excessive snowfall might be the primary factor for the large overestimation of snow depth and cover in ERA5 reanalysis. Pending a solution for this common model precipitation bias over the Himalayas and the TP, future snow reanalyses that optimally combine the use of satellite snow cover and in situ snow depth observations in the assimilation and analysis cycles have the potential to improve medium-range to sub-seasonal forecasts for water resources applications.

1 Introduction

Often referred to as the Third Pole, the Tibetan Plateau (TP) region is the world's highest plateau, with an average elevation of 4000 m above sea level. Due to its spatial extent, elevation and geographical position in the midlatitudes, it exerts a considerable influence on regional and global climate. The formation and variability of the Asian summer monsoon, in particular, is affected by the TP through thermal and mechanical effects (Wu et al., 2012, 2015; Xiao and Duan, 2016), with remote impacts both downstream (e.g. Zhang et al., 2004; Xue et al., 2018) and upstream (Lu et al., 2017). In autumn or winter, snow anomalies over the TP have been linked to wave trains extending downstream over East Asia and the Pacific (Liu et al., 2017; Li et al., 2018). Given its importance for climate and given that climate change strongly affects the region (Yang et al., 2014), the cryosphere over the TP is closely monitored in terms of glacier shrinking (Yao et al., 2012; Treichler et al., 2018), lake expansion (Zhang et al., 2017; Treichler et al., 2018) or change in the snowpack (Chen et al., 2017; Wang et al., 2017, 2018). The extent and variability of the snowpack over the TP has been a major focus of investigation because of the role of snow in the surface energy balance and in the hydrological cycle, as well as its potential impact on the large-scale circulation through radiative or thermodynamical feedbacks (Xiao and Duan, 2016; Lin et al., 2016; Henderson et al., 2018). In addition, water storage over the TP area affects the headwaters of a dozen major Asian rivers, affecting the livelihood of 22 % of the world's population (Immerzeel et al., 2009). Several studies have aimed at quantifying the seasonal, interannual and decadal variability of the TP snowpack (e.g. Basang et al., 2017, and references therein). The seasonal snow cover over the TP is unique compared to other midlatitude regions or to higher latitudes because of its location in the 27–40^∘ N latitude belt, its high elevation and its distinct shallow, patchy and short-lived snowpack. In addition, the TP also receives a high amount of solar radiation; hence the shortwave snow albedo effect tends to be strong over that region. The TP is also a challenging region for snow-related research due to the complex terrain and the relative paucity of in situ observation stations in this vast, sparsely populated expanse, especially in its western part.

There is, nevertheless, a substantial number of meteorological stations over the TP region operated by the Chinese Meteorological Administration (CMA), which have been used in many climatological studies (e.g. for recent studies see Basang et al., 2017; Li et al., 2018). Since most of the stations are located in inhabited valleys, below 4000 m in the southeast TP, the representativeness of these in situ data for the TP as a whole is questionable.

Advances in satellite remote sensing have however provided invaluable information on the state of this high-elevation snowpack (Basang et al., 2017; Basang, 2018; Yang et al., 2015; Ménégoz et al., 2013), often through blending data from multiple instruments to combine the high spatial resolution of optical data with the all-weather capability of microwave (MW) data. Since 2004, NOAA has provided an operational, daily, multi-sensor snow cover product (the Interactive Multisensor Snow and Ice Mapping System or IMS) at a high resolution (4 km) by combining optical, infrared and MW satellite data, and surface observations (Helfrich et al., 2007). The accuracy of the IMS product compared to station data over the TP has been evaluated by Yang et al. (2015) and Li et al. (2018), who found the pixel matching accuracy (snow or no snow) to be over 90 %. Another widely used satellite optical product is the Moderate Resolution Imaging Spectroradiometer (MODIS) 8 d snow cover (Basang et al., 2017). Joint analysis of MODIS and station data by Basang et al. (2017) indicated that, despite large day-to-day snow cover fluctuations and the non-synchronicity of these two datasets, there is a high temporal correlation (0.77) if the station data are preprocessed in 8 d bins in a similar fashion as the MODIS data, i.e. if the station data are considered snow covered in the 8 d period when snow covered for at least 1 d. Nevertheless, Basang et al. (2017) also showed that the spatial distribution of snow cover is uneven: the snow cover fraction is less than 21 % for 70 % of the TP area, and yet it can be up to 40 % in the eastern part of the TP. They also indicated a discrepancy in the seasonality of the maximum snow cover between station data and MODIS: the former had higher snow cover in the spring, when precipitation increases and falls as snow in the southeast TP. On the contrary, MODIS had a higher winter snow cover, which is more representative of the mean conditions over the plateau. These latter discrepancies are indicative of large spatial differences in the characteristics of precipitation. Snow depth products are also available from satellite remote sensing. Snow depth retrieval from MW observations is difficult in regions with complex topography, hampered by non-heterogeneity of snow grain size and by vegetation cover. Few MW snow depth products have been thoroughly evaluated over the TP, a region characterised by a sparse and rapidly melting snowpack where, given the small snow grain size of fresh snow, the MW retrievals can lead to large errors (Dai et al., 2017). Nevertheless, a long-term MW snow depth product has been developed to account for the special conditions of the TP by Che et al. (2008) and Dai et al. (2017).

Atmospheric reanalyses comprehensively assimilate in situ, balloon-borne, aircraft and satellite observations into a forecast model, and they form an essential tool in meteorological and climate research (e.g. Brönnimann et al., 2018; Fujiwara et al., 2017). Furthermore, reanalyses are sometimes used as initial conditions for seasonal hindcasts or reforecasts. A thorough up-to-date description and intercomparison of atmospheric reanalyses is found in Fujiwara et al. (2017). On the other hand, data assimilation in land surface models is performed separately from the atmospheric data assimilation due to the different nature of observations and methodologies involved (Hersbach et al., 2018). Early land reanalyses were performed as offline simulations without any actual assimilation of observations, but land data assimilation is rapidly evolving. For example, the newest operational seasonal forecasting system 5 (known as SEAS5, Johnson et al., 2019) and the Integrated Forecast System (IFS) used for medium-range prediction at the European Centre for Medium-Range Weather Forecast (ECMWF), both rely on a land surface data assimilation system using a two-dimensional optimal interpolation, which blends IMS snow cover and in situ snow depth observations with the model background (de Rosnay et al., 2014, 2015). This evolution is partly motivated by the renewed interest in tapping into the potential of the land surface, and snow in particular, to improve the prediction skill at the sub-seasonal-to-seasonal timescale (Orsolini et al., 2013). Concerning the TP region more specifically, Lin et al. (2016) assimilated snow cover fraction from MODIS and/or water storage from GRACE (Gravity Recovery And Climate Experiment) in a land model and showed improvement in seasonal atmospheric temperature prediction resulting from the improved land conditions. The impact of realistic land initialisation over the TP region on the Indian summer monsoon (ISM) has also been investigated (Senan et al., 2016; Halder and Dirmeyer, 2018; Rai et al., 2015). In particular, the impact of springtime initialisation on predicting the ISM onset was investigated by Senan et al. (2016) using the coupled ocean–atmosphere seasonal forecasting system 4 of the ECMWF. They found that, in years with anomalously high mean snow depth over the TP, as determined from the ECMWF land reanalyses, the ISM onset was delayed by about 8 d. Half of this delay was attributed to the initialisation of snow over the TP region, highlighting the importance of land initialisation over that region for seasonal forecasting of the ISM onset.

While many studies have attempted to validate atmospheric reanalyses over the TP area in terms of temperature, precipitation or snowfall (Wang and Zeng, 2012; Palazzi et al., 2013; Zou et al., 2014; Viste and Sorteberg, 2015), there have been remarkably few studies aimed at evaluating and comparing the snow depth or snow cover in the reanalyses. Snow in long-term climate reanalyses over the plains of central and northern Russia have been evaluated against station data, showing surprisingly good performance (Wegmann et al., 2017). However, the terrain there is not as complex as over the TP. Previous studies of snow evaluation over the TP have focused purely on the remote sensing snow products (Basang et al., 2017; Yang et al., 2015; Dai et al., 2017). We argue that evaluation of snow reanalyses is of interest per se, since these reanalyses provide initial, critical surface information for the forecast systems and link surface conditions with atmospheric dynamics.

The aim of this study is to compare snow cover and snow depth over the TP region among a set of modern atmospheric and land reanalyses, and to evaluate them against in situ snow observations and selected satellite-based remote sensing products. A second goal is also to better characterise the snow biases in the reanalyses and to identify their origin. Sensitivity experiments are then carried out with a land surface model to assess the potential relevance of snow processes in driving the identified biases.

2 Data and methods

Maps of our study area with the orography at the resolution of the atmospheric reanalyses are shown in Fig. 1. We make use of three reanalyses produced by ECMWF, namely the latest generation of atmospheric reanalyses (ERA5) (Hersbach et al., 2018), its older counterpart (ERA-Interim, hereafter ERA-I) (Dee et al., 2011), as well a land reanalysis obtained by running the ECMWF land model (HTESSEL) offline, forced by the ERA5 meteorological forcing. Here, we refer to the latter as ERA5L-CTRL since it is not yet the officially released version of ERA5-Land available online since July 2019 (https://climate.copernicus.eu/climate-reanalysis, last access: 20 August 2019). In addition, we also make use of the Japanese 55-year Reanalysis (JRA-55) (Kobayashi et al., 2015) and of the NASA Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) (Gelaro et al., 2017). We did not carry out an exhaustive inter-comparison of all existing global reanalyses, but we chose the aforementioned analyses because they belong to the latest generation (ERA5, MERRA-2 and JRA-55), or else because they incorporate either local snow observations over the TP region (JRA-55) or satellite snow cover observations that encompass the TP (ERA-I and JRA-55). Some general characteristics of these reanalyses are shown in Table 1, and more detailed information about the treatment and assimilation of snow variables is provided in the Appendix. Further information about the reanalysis systems can be found in Fujiwara et al. (2017) and Wright et al. (2018). For the evaluation of reanalyses, we use station data, the multi-sensor IMS snow cover product and satellite MW data, as detailed below.

Figure 1 Topographic maps of the TP region (altitude in metres) at the horizontal resolution of each reanalysis. The locations of the 33 stations are also shown with the coloured circles indicating the altitude range.

Table 1 General characteristics of the reanalyses.

2.1 Station data

We obtained 5 years (2009 to 2013) of in situ daily meteorological observations from 33 stations over the TP from CMA. The observations comprise minimum and maximum temperature, snow depth and total precipitation. The measuring time is 00:00 UTC (08:00 Beijing time). The operational procedures dictate that when more than 50 % of the ground surrounding the weather station is covered by snow, the snow depth (SD) is measured. When converting in situ snow depth to a snow cover fraction (SCF), a 100 % SCF is assumed when snow depth exceeds 1 cm, according to the in situ observing rules of CMA (CMA, 2003). The in situ observations also provide the snow cover days (SCDs), i.e. the number of days when the snow covers more than 50 % of the ground in the sight from the station. The geographical locations of the stations are shown on Fig. 1, along with their altitudes in comparison to the topography resolved by the four atmospheric reanalyses. Most of the stations are located in inhabited valleys, below 4000 m, in the southeastern part and are not representative of the TP region as a whole. There is only one station in the western part of the TP, west of 85^∘ E. In situ observations have several sources of uncertainty. Here, we highlight two sources: (i) the stations might not be fully representative of their local surroundings due to the complex nature of the terrain, and (ii) the quality of the records could be affected by the harsh operating conditions. For example, strong winds limit the instrument ability to record the amounts of falling snow or solid precipitation, a phenomenon called undercatch.

2.2 Interactive Multisensor Snow and Ice Mapping System (IMS)

We use SCF from IMS, a multi-instrument, near-real time daily product covering the Northern Hemisphere with a pixel resolution of 4 km. IMS provides a binary (1/0) snow cover information: either 1 if more than 50 % of the 4 km pixel is covered by snow, or otherwise 0 (snow free). Since, in this paper, we compare IMS data to reanalyses, we aggregate the 4 km product to a 0.25^∘ grid comparable to the highest horizontal resolution among the reanalyses. This is done by counting the number of pixels with a value of 1 in a grid box, assuming they represent 100 % cover – this gives the high estimate. If we assume that a value of 1 represents only 50 % of the 4 km pixel we obtain the low estimate, for which the maximum SCF possible is 50 %. These two estimates provide a range of values, reflecting the uncertainty inherent to aggregating the 4 km binary data, e.g. a value of 1 in a pixel means 50 % to 100 % snow coverage.

The choice of IMS snow cover is motivated by its use in the ECMWF analysis system. A key point is that it is used differently in ERA-I (Drusch et al., 2004) than in ERA5 (de Rosnay et al., 2014, 2015; Hersbach et al., 2018). First, the high-resolution 4 km product is used in ERA5, while the coarser-resolution 24 km product is used in ERA-I. More importantly, IMS data are not used above 1500 m in ERA5, i.e. in high-altitude regions such as the TP, while it was still in use in ERA-I. When comparing the snow cover in IMS with the one in ERA-I or ERA5, one has to recall that these are not independent datasets. Nevertheless, the different usage of IMS data between ERA5 and ERA-I allows us to highlight the importance of snow cover analysis and assimilation over the TP. In that region, in situ data available for the numerical weather prediction community for real-time applications is still scarce.

2.3 Microwave satellite data

A long-term (1978–2010) MW daily snow depth product has been developed to account for the special conditions of the TP by Che et al. (2008) and Dai et al. (2017) at the Cold and Arid Regions Environmental and Engineering Research Institute (CAREERI) of the Chinese Academy of Sciences. We used the 2 years (2009 and 2010) that overlap with our comparison period (2009–2013), the latter dictated by our station data record. Over these 2 years, the MW product is based on measurements from the NASA Advanced Microwave Scanning Radiometer for the Earth Observing System (AMSR-E). The gridded dataset has a horizontal resolution of 0.25^∘.

3 Results

3.1 Snow depth

Figure 2 presents a comparison of daily SD over the 5-year period (2009–2013) between the observations, comprising both in situ station and MW data, and the four reanalyses collocated at the station coordinates. The comparison is an average over the 33 stations. The in situ observations, even in the mean over all available stations, reveal a very thin and rapidly fluctuating snowpack, with episodes of fast build-up and decay, followed by periods void of snow, even in mid-winter. This concurs with previous studies showing that large parts of the TP can remain snow-free even in winter (Basang et al., 2017; Li et al., 2018). A quantitative estimate of the discrepancy between in situ daily observations and the other datasets as a root mean square error (RMSE) is presented in Fig. 4 (left panel). The RMSE is calculated for the 2 years 2009 and 2010 when satellite MW data are available. We also calculated the temporal correlation matrix between the different SD datasets (Table 2), using daily data year-round over the 5-year period, except for the satellite MW data where only the 2 years 2009 and 2010 are used. It is clearly apparent that, with the exception of MERRA-2, the reanalyses show a regular seasonal cycle, with a snowpack that grows nearly steadily during the cold season and culminates in February or March. This is unlike the in situ observations, which show rapidly fluctuating snow increases with little snow accumulation throughout winter. In comparison with in situ observations, MERRA-2 has the best performance among reanalyses for both the RMSE and the temporal correlation, closely followed by JRA-55. The ERA5 reanalysis overestimates the seasonal maximum SD by a factor of 10. Among the ECMWF family of reanalyses, the older ERA-I performs best in terms of RMSE. Nevertheless, the temporal correlation is similar to the one obtained with ERA5, indicating that the newer, higher-resolution ERA5 similarly captures the snow variability despite its large positive bias. The MW satellite data are also overestimating SD, as noted by Yang et al. (2015) or Dai et al. (2017). Like the reanalyses, they show a progressive snow accumulation throughout the cold season. They also poorly captures short-term SD variability (Fig. 2), and its temporal correlation with the in situ data is 0.32, poorer than for the reanalyses, even though calculated over only 2 years. The RMSE error is nevertheless comparable to ERA-I, and the MW data are able to estimate SDs smaller than 5 cm.

Figure 2 Comparison over the 5-year period 2009–2013 of daily SD (in centimetres) between the observations, comprising both in situ station and MW data (over 2009 and 2010), and the four reanalyses collocated at the station coordinates, namely ERA5, ERA-I, MERRA-2 and JRA-55. The average is carried out over the 33 stations shown in Fig. 1.

Table 2 Temporal correlation matrix between in situ observations and reanalyses or satellite data for daily SD (a) and daily SCF (b). The IMS low estimate has been used. Correlations are calculated over the whole 5-year period, except for MW SD data for which the years 2009 and 2010 are used and are averaged over the 33 stations.

It is not surprising that JRA-55 performs well, since it assimilates SDs from a dozen CMA stations over the TP area and snow cover from satellite MW data (see Appendix). Note that the assimilated MW data are not from the same instrument as the one used in the CAREERI MW data. A key factor in the relatively good performance of ERA-I is the fact that it assimilates SCF from IMS in high-altitude regions, hence also over the TP. On the contrary, as mentioned earlier, assimilation of IMS SCF above 1500 m was discontinued in the production of ERA5. Among the differences in the snow treatment between ERA-I and ERA5, the tendency to reduce or remove the snowpack imposed by the IMS observational constraint during assimilation appears to be the most likely key factor in bringing ERA-I SDs significantly closer to the in situ observations.

3.2 Snow cover and snow maps

Figure 3 presents a comparison of daily SCF over the 5-year period (2009–2013) between the in situ station observations and the IMS satellite data and the five reanalyses, collocated at the station coordinates. Again, an average over the 33 stations is presented. The IMS data are represented by the range between the low and high estimates (see Sect. 2). The correlation matrix and the RMSE for SCF are provided in Table 2 and in Fig. 4b respectively, calculated over the 5-year period. There is good agreement between in situ observations and the IMS low estimate: the station-mean year-round correlation (Table 2) over the 5 years (2009–2013) is 0.78. In Li et al. (2018), the correlations between SCF at 55 CMA stations and IMS SCF during the 10 winters between 2000 and 2010 ranged from 0.39 to 0.79, with a 10-winter average of 0.56. The consistency of the temporal correlation with Li et al. (2018) indicates that the satellite snow cover data readily capture the rapidly fluctuating snowy events. It also provides some confidence in the station data despite the harsh operating conditions for in situ measurements, and the spatial degradation applied to IMS (from 4 km to 0.25^∘).

Figure 3 Comparison of daily SCF (between 0 and 1) over the 5-year period 2009–2013 between the observations, comprising both the in situ station and IMS satellite data, with the four reanalyses collocated at the station coordinates, namely ERA5, ERA-I, MERRA-2 and JRA-55. The data have been averaged over the 33 stations. The IMS data are represented by the low and high estimates, resulting from converting the binary 4 km data to a grid comparable to the resolution of the reanalyses.

Figure 4 Estimate of the root mean square error (RMSE) between daily station data and collocated reanalyses or satellite observations for SD (in centimetres) (a) over 2 years when satellite MW data are available (2009 and 2010) and for SCF (between 0 and 1) (b) over the 5-year period (2009–2013).

Table 3 Values of the annual means of monthly snow cover days averaged during the whole 5-year period (the second column) and during the years 2009 to 2013 (the third to seventh columns) for each dataset (a). Temporal correlation matrix between in situ observations and reanalyses for monthly SCD during the whole 5-year period (b). The SCDs are defined as the days when the SCF is greater than 0.5 and are averaged over the 33 stations.

It is clearly apparent that ERA5 again considerably overestimates SCF. JRA-55 is much worse than MERRA-2 for SCF, while their performance for SD was similar. As described in the Appendix, the a posteriori conversion from SD to SCF varies among reanalyses. For JRA-55, a thin 2 cm layer is equivalent to 100 % SCF, and hence JRA SCF is persistently high. On the contrary, in MERRA-2, SCF values are small due to the large SD required to have full snow cover (see Appendix for details). While the in situ station data (e.g. Fig. 3) also display thin layers of a few centimetres, they are very transient. Further validation is provided in terms of the annual-mean SCDs. Table 3 compares SCDs among the datasets, both for the 5-year average and for individual years, based on monthly-mean SCD values. The values of SCDs compare well between in situ observation and IMS. On the other hand, ERA5 and JRA-55 reanalyses largely overestimate SCDs (ERA5 by nearly 3 times), while ERA-I is closer to the observed values. The temporal correlation coefficients for monthly SCDs tend to be higher than that for the daily SCF or SD. Note the values of SCDs in MERRA-2 are small due to the SCF being below the 0.5 threshold used in the SCD definition (see Fig. 3).

Figure 5 Maps of the 5-year mean SD (in centimetres) over the TP region in January for each of the four reanalyses (ERA5, ERA-I, MERRA-2 and JRA-55) with the in situ SD (in centimetres) embedded locally in each map.

Figure 6 Maps of the 2-year mean SD January anomalies (in centimetres) from the satellite MW data for the years 2009 and 2010 over the TP region for each of the four reanalyses (ERA5, ERA-I, MERRA-2 and JRA-55).

Moving to the geographical distribution over the TP region, Fig. 5 shows maps of the 5-year mean SD in January for each of the four reanalyses, with the in situ data embedded in each map. January is chosen to illustrate a mid-winter month when snowfall is relatively common and differences between reanalyses are large. In the southeast TP, where the stations are located, MERRA-2, JRA-55 and ERA-I have SDs comparable to in situ observations, as shown in Fig. 2. In the western station-void part, there are also large differences in the snow depth (up to a factor of 5) among the reanalyses, especially along the arc of the Himalayas and other high mountain ranges. Figure 6 is similar to Fig. 5 but represents January anomalies from the satellite MW data for the combined years 2009 and 2010. In the southeastern part, SDs in MERRA-2, ERA-I and JRA-55 are smaller than the MW data, consistent with Fig. 2. Figure 6 also reveals that ERA5 and JRA-55 have an excessive SD over the high mountains of the western Himalayas compared to satellite MW data, unlike ERA-I with its high-altitude snow cover assimilation. Figure 7 shows maps of the 5-year mean SCF in January for each of the reanalyses, for the two (low and high) estimates from IMS, again with the in situ SCF embedded in the maps. While there is a good agreement between the in situ data and the IMS low estimate, as shown earlier at station locations, the overestimation by ERA5 extends over the eastern and western parts of the TP. Nevertheless, some parts of central Tibet and the Taklamakan desert to the north are snow free in January in both IMS and in most reanalyses. ERA-I agrees much better with the IMS data, while JRA-55 SCF is too high, due to the SD vs. SCF conversion as explained above. On the contrary, the MERRA-2 SCF is exceedingly low and featureless.

Figure 7 Maps of the 5-year mean SCF (between 0 and 1) over the TP region in January for the IMS satellite data (two maps for the high and low estimates) and for each of the four reanalyses (ERA5, ERA-I, MERRA-2 and JRA-55), with the in situ SCF observations (between 0 and 1) embedded locally in each map.

4 Discussion

4.1 Snow cover to snow depth conversion

A first important issue to mention is that the conversion of SD to SCF differs significantly among reanalyses. SCF is perhaps the most important snow-related climate driver for the surface budget, since the shortwave snow albedo feedback is strong on the TP, especially in the spring. The snowpack is also quite thin so that thermodynamical feedbacks linked to snow thickness are less important. A 100 % SCF might correspond to snow depths ranging from 26 to 2 cm, depending on the reanalysis considered. There may also be some dependence on snow density as in the cases of ERA-I and MERRA-2. Hence, while JRA-55 has an excellent performance among reanalyses for SD, the snow cover is exaggerated due to the conversion to 100 % SCF obtained with a small SD threshold value of 2 cm. While MERRA-2 SCF compares best with station data among the reanalyses (see Fig. 4b), being close to the low IMS estimate, it is exceedingly low over the TP overall (Fig. 7) given the comparatively high value (26 kg m⁻²) used for snow mass to get a fully covered snow area (Reichle et al., 2017a).

4.2 Sensitivity in offline ERA5-Land simulations

We carried out simple sensitivity offline experiments with the ERA5-Land model to shed some light on potential explanations for the large biases in the ERA5 reanalysis, which might be partly due to a missing snow removal process or to excessive snowy precipitation. Gordon et al. (2006) reported that blowing snow events are common in the Canadian Prairies and in the Arctic. As snow is carried out in suspension, it may sublimate. In a model study, Pomeroy and Li (2000) predicted that two-thirds of the snow transported from the surface to the atmosphere sublimates in the Prairie regions and one-half sublimates in the Arctic regions. Brun et al. (2013) used Gordon's parameterisation and argued that this was an important process in simulating snow dynamics in northern Eurasia. Processes like the blowing of snow and the sublimation of the blown snow might be important processes in the windy and dry conditions of the TP. To test the hypothesis that sublimation of blowing snow is an important missing process, a simple parameterisation of blown snow sublimation has been introduced, as described in the Appendix. It has been tested for two wind thresholds for blowing snow initiation. Furthermore, the impacts of changing the threshold in the snow depth to snow cover conversion and of reducing snowfall were tested separately. Other processes such as the impacts of blown snow on surface roughness or of snow drifts were not tested. In addition to the ERA5L-CTRL run described in Sect. 2, we hence performed an additional set of five land-only experiments testing (i) the effect of snow sublimation due to blowing snow in BLW and in BLW_L with a lower threshold for blowing snow initiation, (ii) the 100 % SCF threshold, changed from 10 to 5 cm in SCF05 and to 20 cm in SCF20, and (iii) the excessive snowfall by reducing snowfall by 50 % in SNF50. These five experiments with ERA5-Land, as summarised in Table 4, were carried out for the period 1 January 2019 to 31 December 2014, and over the domain 60.125 to 109.875^∘ E and 25.125 to 44.875^∘ N, with a regular latitude and longitude resolution of 0.25^∘. The simulations were driven by the ERA5 meteorological fields interpolated to this regular grid.

Table 4 Details of the offline ERA5-Land experiments.

Table 5 ERA5 and offline ERA5-Land experiments mean annual snow and associated fluxes averaged over the 33 stations. The fluxes are the annual mean accumulation of the following: snow sublimation, snow sublimation due to blowing snow, snowmelt, and snowfall and rainfall intercepted in the snowpack. The last row (OBS) presents the mean snow depth derived from the station data and the mean snow cover from the IMS satellite data (low estimate).

The evaluation of the different budget terms over the 33-station mean snow depth and cover (Table 5) as well as the root mean square errors and the temporal correlation with observations (Table 6) are shown. These results indicate that (i) the offline simulation (CTRL) with the same configuration as ERA5 represents closely the errors found in ERA5 as expected, (ii) the effects of snow sublimation due to blowing snow and of the SCF threshold change are negligible, and (iii) artificially reducing snowfall by 50 % reduces the systematic biases significantly, and it increases the temporal correlation of snow depth as well. The temporal variability of snow cover is captured reasonably well by ERA5 and the other experiments, suggesting that the main dynamic features of snow accumulation and depletion are well represented.

Table 6 Root mean square error (RMSE) and temporal correlation for snow depth (versus stations data) and snow cover (versus IMS satellite) for ERA5 and the different experiments. The time series of the observation and experiments were first averaged over the 33 stations before the score calculations were carried out.

4.3 Excessive precipitation issue

These results suggest that excessive snowfall might be a key factor responsible for the large overestimation of snow depth and cover in ERA5 reanalysis. Excessive snowfall over this region is in fact a common bias among climate and forecast models (Su et al., 2013). In fact, the largest relative climate model bias for precipitation over land, globally, is on the TP (Flato et al., 2013). In order to provide further evidence for the precipitation bias, we present in Fig. 8 the 5-year mean seasonal cycles of minimum temperature, maximum temperature, total precipitation rate (i.e. liquid and snowy precipitation), snow cover and snow depth based on monthly means and averaged over all the 33 stations. Figure 8 reveals that all reanalyses have a cold temperature bias compared to in situ observations, especially in maximum temperature, which is largely consistent with their respective thick snowpack bias. For example, ERA-I is warmer and closer to the in situ observations than ERA5, likely due the latter having a higher snow depth. The good performance of ERA-I (and also of the first-generation MERRA) against in situ station temperature data has been noted by Wang and Zeng (2012). All reanalyses have a precipitation excess (Fig. 8c), except MERRA-2 which uses observed precipitation data, including over the TP region. Incidentally, this may explain the excellent performance of MERRA-2 in terms of mean snow depth over the stations (see Fig. 4a). Palazzi et al. (2013) also reported the precipitation excess in ERA-I, for example. In fact, ERA-I displays a precipitation excess greater than the more recent ERA5, despite having smaller SDs. The improvement in precipitation in ERA5 is likely due to its higher horizontal resolution, allowing better representation of synoptic events. Note that the precipitation seasonal cycle is largely dominated by summer monsoonal precipitation in the southeast TP. While our focus here is on the cold season, when precipitation is much smaller but likely falling as snow, there is still a positive bias in the reanalyses, indicating excessive snowfall. It must be kept in mind that in situ observations of snowfall are associated with a large uncertainty due to undercatch and other technical challenges.

Figure 8 The 5-year mean (2009–2013) monthly-mean seasonal cycle for minimum temperature (^∘C) (a), maximum temperature (^∘C) (b), precipitation rate (in millimetres per day) (c), SD (in centimetres) (d), and SCF (between 0 and 1) (e) averaged over all the 33 stations. The grey line refers to IMS (low estimate) for snow cover or to MW for snow depth.

Hence, a remaining point that must be addressed in future studies is the origin of the seasonal moisture transport to the TP region leading to the excess precipitation in models. While the results of the offline sensitivity tests do not support a role for sublimation due to blowing snow, it could be that the local turbulent surface winds are much higher than in the reanalyses. Hence the importance of the blown snow sublimation cannot be fully discounted. We further note that there is a strong regional seasonality in the precipitation pattern over the TP. The wet season precipitation (May–September) accounts for more than 70 % of the total annual precipitation over the south and southeast TP (where stations are located), a region that falls under the influence of the Asian summer monsoon (Maussion et al., 2014). On the other hand, over the western and northwestern parts of the TP , the wet season precipitation is less than 300 mm. This region is influenced by strong winter and spring westerlies and transient migratory synoptic eddies embedded in them (the so-called westerly disturbances), which provide a significant portion of the annual mean precipitation, especially over the western Himalayas (Tiwari et al., 2017). An effective barrier effect, inhibiting large-scale moisture transport over the region, is not captured by any of the models underpinning the reanalyses used here. The impact of having higher model horizontal resolution to resolve the complex topography has to be considered. Lin et al. (2018) carried out simulations for the summer monsoon season with the Weather Research Forecasting model (WRF) at resolutions of 30, 10 and 2 km and found that the finest resolution (2 km) diminished the water vapour transport to the TP and the precipitation bias. They also noticed a large improvement when model resolution changed from 30 to 10 km. The insufficient sub-scale orography variance and orographic drag seems to be a key factor (Zhou et al., 2019). Based on these results, it appears that the highest resolution in atmospheric reanalyses examined here (near 32 km in ERA5) remains insufficient, as demonstrated by the precipitation excess compared to station data in Fig. 8c.

Finally, we note that the particular relevance of the western TP snowpack for sub-seasonal-to-seasonal prediction, due to the limited snowpack persistence over the eastern and central TP, was already emphasised by Xiao and Duan (2016). Yet, with only one station in our inter-comparison, the western TP snowpack remains the less constrained by in situ data, and further studies on this issue are warranted.

5 Summary

We have shown that several recent, state-of-the-art reanalyses produce an over-extensive snowpack in autumn, winter and spring over parts of the TP. This is at odds with observational studies revealing that snowfall events are very transient, that the snow cover vanishes rapidly on timescales of days and that large parts of the TP can remain snow-free in winter (e.g. Basang et al., 2017; Li et al., 2018). The reanalyses that assimilate local in situ observations or satellite snow cover or precipitation are more capable of representing the shallow, transient snowpack over the TP region. MERRA-2 and JRA-55 have the best performance among reanalyses for snow depth. Considering the family of reanalyses from the ECMWF, we surmise that the underperformance of ERA5, in terms of SD compared to its older counterpart ERA-I, is due to the lack of IMS data assimilation at high altitudes, including over the TP. This issue will be subject of a follow-up publication, using dedicated assimilation experiments. Pending a solution for the common model precipitation bias over the Himalayas and TP, which may require much higher horizontal resolutions than currently used in global reanalyses, improved snow initialisation through better use of observations in the analyses may improve the medium-range to sub-seasonal forecasts.

Appendix A: ERA-Interim

Snow is represented as a single bulk layer on top of the soil column with prognostic temperature, mass, density and albedo. Snow density is constant with depth, increasing exponentially with snow age, and decreasing after snowfall events assuming a constant snowfall density of 100 kg m⁻³. Snow albedo reduces exponentially as snow ages over low vegetation or bare soil, but is constant in time for snow under high vegetation, and is reset to a maximum value of 0.85 after snowfall event.

A snow depth analysis is performed using a Cressman (1959) interpolation with successive corrections. Station observations of snow depth are used; however there are no stations used over the TP. Gridded snow cover from IMS is also assimilated since 2004, using the coarser 24 km resolution product as detailed in Drusch et al. (2004). The snow depth analysis is relaxed toward a climatology when observations are unavailable.

Snow cover fraction is a diagnostic variable derived from the relation SCF $=min(\mathrm{1},\mathrm{SWE}/\mathrm{15})$, where SWE is the snow mass (kg m⁻²). A layer of 15 kg m⁻² represents 100 % snow cover (15 cm depth considering a snow density of 100 kg m⁻³ or 5 cm depth if considering a snow density of 300 kg m⁻³).

Appendix B: ERA5

The snow model is also a bulk single layer as in ERA-I, with the same prognostic variables but with several modifications, in particular a diagnostic of snow liquid water content, a revised SCF formulation and snow density evolution (Dutra et al., 2010). A two-dimensional optimal interpolation snow analysis, including the water equivalent, temperature and density of snow is performed using station observations of snow depth and the gridded SCF product from IMS 4 km resolution product (de Rosnay et al., 2014, 2015). Unlike ERA-I, IMS data are not used above 1500 m, i.e. in high altitude regions, which includes the TP. No station data over the TP are used. Unlike ERA-I, the snow depth analysis is not relaxed toward a climatology when observations are unavailable.

Snow cover fraction is a diagnostic variable derived from the relation SCF $=min(\mathrm{1},\mathrm{SD}/\mathrm{10})$, where SD is the snow depth in centimetres so that a layer of 10 cm represents 100 % snow cover.

Appendix C: ERA5L-CTRL and dedicated experiments

We made several additional offline simulations of the HTESSEL model forced by meteorological data from ERA5. The land model is hence the same as that used in ERA5, but there is no assimilation of snow or land observations in contrast to the ERA5 reanalyses. A simple parameterisation of snow sublimation due to blowing snow, as proposed by Gordon et al. (2006) and used in Brun et al. (2013), is newly implemented in HTESSEL. The sublimation of blowing snow (Q_s, kg m⁻² s⁻¹) is computed as in Gordon et al. (2006):

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}& Q_{\mathrm{s}}=A{\left({\displaystyle \frac{T_{\mathrm{0}}}{T_{\mathrm{a}}}}\right)}^{\mathit{\gamma }}{U}_{\mathrm{t}}{\mathit{\rho }}_{\mathrm{a}}{q}_{\mathrm{sa}}\left(\mathrm{1}-{\mathrm{RH}}_{\mathrm{a}}\right){\left({\displaystyle \frac{U_{\mathrm{a}}}{U_{\mathrm{t}}}}\right)}^B\\ & \mathrm{with}{U}_{\mathrm{a}}>U_{\mathrm{t}},\end{array}\end{array}$$

where A and B are dimensionless constants (A=0.0018, B=3.6) and γ=4, T₀ is the water freezing temperature. The formulation requires the lowest model level (at 10 m height) fields of air temperature (T_a), air density (ρ_a), wind speed (U_a), relative humidity (RH_a) and saturation specific humidity (q_sa). U_t (m s⁻¹) is the threshold for initiation of blowing snow also used by Gordon et al. (2006) following Li and Pomeroy (1997):

$$\begin{array}{}\text{(C2)}& U_{\mathrm{t}}=\mathrm{6.98}+\mathrm{0.0033}{\left(T_{\mathrm{a}}-\mathrm{245.88}\right)}^{\mathrm{2}}.\end{array}$$

In the experiments with reduced wind threshold, the first coefficient (6.98) is replaced by 5. This simple parameterisation only represents the snow mass removal by sublimation due to the blowing of snow, but it does not account for the change in the energy budget. Hence, the energy required to sublimate the snow is not taken from the surface or from the lower atmosphere. While this approximation would not be valid in a coupled land–atmosphere simulation, it seems appropriate in an offline simulation. Other experiments tested the impacts of the threshold in snow cover conversion (see SD parameter in the description of ERA5 above) and of reducing snowfall by half. The experiments are summarised in Table 4.

Appendix D: JRA-55

The land surface model in JRA-55 is the Simple Biosphere model (SiB) (Sellers et al., 1986), which represents snow mass on the ground and its evolution by updating several parameters and calculations. A separate optimal interpolation-based snow depth analysis is performed once per day (at 18:00 UTC). The first-guess background state combines the land-surface analysis and Special Sensor Microwave/Imager (SSM/I) and the Special Sensor Microwave Imager Sounder (SSMIS) snow cover satellite observations. The analysis also considers in situ observations of snow depth over China from CMA archives, including a dozen stations over the TP (Onogi et al., 2007; Kobayashi et al., 2015).

Snow cover fraction is a diagnostic variable derived from the relation SCF $=min(\mathrm{1},(\mathrm{SD})/\mathrm{2})$, where SD is the snow depth in centimetres so that a layer of 2 cm represents 100 % snow cover.

Appendix E: MERRA-2

The evolution of snow mass, depth and heat content is simulated within each catchment using a three-layer snow model (Reichle, et al., 2017b). Density in each layer evolves via parameterised representations of compaction due to snow overburden, and melting and refreezing (Stieglitz et al., 2001). Snow is redistributed among layers as necessary to keep the surface layer shallow enough to respond to diurnal variability. The albedo of snow-covered land surface depends on snow density and vegetation type.

No snow data assimilation was performed. However, MERRA-2 performs an online correction for precipitation using two different gridded precipitation datasets from NOAA Climate Prediction Center, namely the Unified Gauge-Based Analysis of Global Daily Precipitation (CPCU) product and the Merged Analysis of Precipitation (CMAP) product (Reichle et al., 2017a, b).

Snow cover fraction is a diagnostic variable derived from the relation SCF $=min(\mathrm{1},\mathrm{SWE}/\mathrm{26})$, where SWE is the snow mass in kilograms per square metre. A layer of 26 kg m⁻² represents 100 % snow cover (i.e. a 26 cm snow depth considering a density of 100 kg m⁻³).