Precipitation is a critical input to land surface and hydrology modeling and prediction. Dynamical downscale modeling has added value to representing precipitation, when compared with the performance of coarse-resolution reanalysis and global climate models, over the Tibetan Plateau (TP). Convection-permitting modeling (CPM) may even outperform dynamical downscale models (DDMs). In this study, 4-km CPM results were compared to 28-km DDM results for a snow season (1 October 2013–31 May 2014) over the TP. The CPM- and DDM-simulated precipitation, as well as three merged gridded precipitation datasets, were evaluated against in situ observations below 4800 m. The five precipitation datasets (CPM, DDM, CMFD, COPRPH, and TRMM) showed large differences over the TP with underestimation of TRMM and overestimation of CPM and DDM compared to observations. The most significant difference occurred in the Brahmaputra Grand Canyon. Given the substantial uncertainty in observed precipitation at high mountains, snow cover simulated by a high-resolution land data assimilation system was used to indirectly evaluate the above precipitation data using MODIS observations. Simulated snow-cover fraction was greatly underestimated using all the merged precipitation datasets. However, simulations using the DDM- and CPM-generated precipitation as input outperformed those using other gridded precipitation data, showing lower biases, higher pattern correlations, and closer probability distribution functions than runs driven by the merged precipitation. The findings of this study generally support the assumption that high-resolution CPM-produced precipitation has an added value for use in land surface and hydrology simulations in high-mountain regions without reliable in situ precipitation observations.
Land surface models (LSMs) have undergone significant improvement since they were first developed in the 1980s; however, large uncertainties still exist when applying them to the complex-terrain Tibetan Plateau (TP) (Zhang and Li 2016; Gao et al. 2017). Previous studies have focused mainly on improving the model physics and associated parameters. For instance, the vertically heterogeneous soil and soil organic matter within the topsoil is an important influence on the profile of soil moisture (Yang et al. 2009), reducing the saturated conductivity from the mucilage in the rhizosphere produces better results in topsoil moisture and evapotranspiration (Gao et al. 2015c), and modeling snow albedo reveals considerable impacts on the interaction between the land surface and the atmosphere (Gao et al. 2017). However, atmospheric forcing conditions used to drive LSMs are important uncertainties, which are sometimes more uncertain than parameterization schemes. For example, Gao et al. (2015c) found that different forcing exerts greater impacts on evapotranspiration and surface runoff accumulation because of the precipitation differences. Zhang et al. (2016) claimed that uncertainty in precipitation data exerts more influence on the general behavior of Noah-MP ensemble simulations than the vegetation phenology parameter.
In light of the importance of precipitation in land surface modeling, numerous attempts have been made to achieve reliable gridded precipitation datasets (Sun et al. 2018), such as interpolations based on station records, combining satellite and station observations, and climate modeling. However, most stations in the TP are located in the central and eastern part, and there is no single station situated above 4800 m in the western TP. Furthermore, precipitation records at high altitudes are always somehow problematic due to various limitations, such as precipitation under catchments (Yang et al. 1998; Scaff et al. 2015; Pan et al. 2016), low station density (Maussion et al. 2011; Gao et al. 2015a; Xiao et al. 2016; Li et al. 2017), and the representativeness of the stations (i.e., most of them are in valleys). Thus, substantial discrepancies exist in extrapolations, both spatially and temporally, and this is particularly true for high-altitude mountain regions such as the TP (Gao et al. 2018b).
Attempts have also been made to merge remote sensing and in situ observational records, including the widely used NOAA Climate Prediction Center Merged Analysis of Precipitation (Xie and Arkin 1997), the Tropical Rainfall Measuring Mission (TRMM) (Huffman et al. 2007; Huffman and Bolvin 2014), and the Global Precipitation Climatology Project (Huffman et al. 2001; Adler et al. 2003). However, simulating and assimilating the land surface and hydrology requires subdaily atmospheric variables as inputs rather than coarse temporal resolution (e.g., monthly) datasets. High-frequency precipitation products were generated, such as NOAA’s Climate Prediction Center morphing technique (CMORPH) (Joyce et al. 2004), the Precipitation Estimation from Remotely Sensed Information Using Artificial Neural Networks–Climate Data Record (PERSIANN-CDR) (Ashouri et al. 2015) products, and the Asian Precipitation–Highly Resolved Observational Data Integration Toward Evaluation of Water Resources (APHRODITE) (Yatagai et al. 2009, 2012). Satellite remote sensing provides precipitation information for a broader range of temporal and spatial scales, representing a significant contribution toward mapping gridded rainfall, but thorough verifications are necessary. However, this is hindered by the scarcity of the observation and derivation approaches resulting in large uncertainty in mountainous regions.
Global climate modeling has the potential to produce global precipitation datasets covering the entire TP, but they show substantial biases over the TP owing to the coarse resolution and model physics (Gao et al. 2015a, 2017, 2018a,b). To achieve finer-scale weather and climate information, statistical downscaling (Wilks 1995; Benestad 2004), and dynamic downscaling models (DDMs) (Gao et al. 2011, 2015a, 2018a,b) have been used to provide regional climate information with high spatial and temporal resolutions. The accuracy of the statistical method depends on the sample pool sizes of the observations; hence, accuracy is compromised in data-scarce regions such as the TP region. The dynamic downscaling method uses a regional climate model, in which precipitation over a region of interest is generated based on the atmospheric dynamics and physics of the regional climate model used although the quality of the lateral boundary conditions exert large impact on its performance (Prein et al. 2015; Xu and Yang 2015; Gao et al. 2017). This method has been widely applied for downscaling reanalysis data and global climate simulations to study regional climate and climate change in North America, Europe, Africa, and Asia (Giorgi et al. 1992; Laprise et al. 1998; Kim et al. 2002; Déqué et al. 2005; Duffy et al. 2006; Gao et al. 2006; Ghan et al. 2006; Leung et al. 2006; Plummer et al. 2006; Mearns and Team 2009; Zhang et al. 2009).
The finest resolution of DDMs available over the TP thus far is a quarter degree (Gao et al. 2015a,b, 2017, 2018a,b). They reduce the overestimation of precipitation by around 35% compared to ERA-Interim (Gao et al. 2015a), capture the observed elevation-dependent warming better than the forcing reanalysis (Gao et al. 2015a), and reproduce the observed northwest–southeast contrast in dryness and wetness changes over the TP better than coarse-resolution reanalysis (Gao et al. 2015b). However, substantial overestimations of precipitation still exist, especially for the nonmonsoon season (Gao et al. 2015a).
Recently, convection-permitting modeling (CPM) has produced more realistic precipitation than coarse-resolution DDMs over complex terrain in the North American Rockies (Rasmussen et al. 2011; Prein et al. 2013b; Letcher and Minder 2015; Liu et al. 2016); in European Alps (Prein et al. 2013a; Ban et al. 2014; Coppola et al. 2018) and North American high latitudes (Li et al. 2019). For the TP, Lin et al. (2018) found that finer-resolution simulations (specifically 2 km) reduce the high precipitation bias over the TP. Lundquist et al. (2019) asserts that today’s ability to simulate mountain precipitation is approaching or has even exceeded the ability to characterize mountain precipitation from observations. He et al.’s (2019) precipitation-data comparison study raised a similar argument.
In this study, a CPM spanning eight months of a cold season (from October 2013 to May 2014) was run using the Weather Research and Forecasting (WRF) Model (Skamarock 2004) centered on the TP. The simulated precipitation was evaluated and compared with China Meteorological Administration (CMA) station observations as well as other available merged gridcell datasets. At high elevations where there is no in situ observation, the snow cover simulated by a High-Resolution Land Data Assimilation System (HRLDAS) (Chen et al. 2007) using the CPM precipitation as input was used as an indirect evaluation. The purpose of this work is to answer three key questions:
What are the main differences in current precipitation datasets at high altitudes in the TP among in situ, remote sensing, DDM, and CPM?
Compared to coarser-resolution DDM, does the CPM add value to the simulation of precipitation, in terms of amount, frequency, and intensity, over the TP?
Does the CPM add value to high-resolution simulation of snow cover over the TP?
2. Models, datasets, and methodology
The WRF Model (http://www.wrf-model.org/index.php) is a mesoscale numerical weather prediction system (Skamarock 2004; Powers et al. 2017). It has a nonhydrostatic core and is developed by several research institutes in the United States. It includes many choices for physical parameterizations suitable for a broad spectrum of applications across scales from meters to thousands of kilometers. The dynamic cores in WRF include a fully mass- and scalar-conserving flux from the mass coordinate version. The physics packages include microphysics, cumulus parameterization, planetary boundary layer, shortwave and longwave radiation, and land surface models. In this simulation, the shortwave and longwave radiation were computed by the NCAR’s Community Atmosphere Model shortwave scheme and longwave scheme (Collins et al. 2004), following previous DDMs in our group (Gao et al. 2015a,b). The WRF single-moment 6-class scheme, the Kain–Fritsch (Kain 2004) convection scheme, and the Yonsei University planetary boundary layer scheme (Hong and Pan 1996) were used. The land surface model used here was the Noah LSM four-layer soil temperature and moisture model with frozen soil and snow-cover prediction (Chen and Dudhia 2001). WRF has been used in dynamical downscaling in the TP region to investigate the added value of downscaling when compared with GCMs, and the relative roles of land surface processes and large-scale forcing in dynamic downscaling (Gao et al. 2015a,b, 2017, 2018a). Version 3.8 of WRF was used for this study.
The DDM run was conducted with a 28-km horizontal grid spacing, using 280 × 184 grid cells centered at 37°N, 92°E, which covered nearly the entire Eurasian continent (Fig. S1a in the online supplemental material). The CPM run was conducted with a 4-km horizontal grid spacing (750 × 414 grid cells), which included the entire TP. The topography in the CPM domain is illustrated in Fig. S1b. The DDM and CPM runs share the same parameterizations except that the convection scheme was not used in the CPM run. The vertical levels were set to 27, with the model top at 50 hPa. The whole 2013–14 snow season in a normal climatology year was simulated. The simulation was initialized at 0000 UTC 1 October 2013 and ended at 2300 UTC 31 May 2014. The lateral boundary conditions and SST were updated every 6 h from the ERA-Interim reanalysis dataset (Dee and Uppala 2009; Dee et al. 2011), which has been proven to be the best among the various reanalysis products available for describing the water cycle over the TP (Gao et al. 2014, 2015a,b). Given the huge overestimations over lakes using the nearest SST, observation-adjusted surface air temperature was used as the skin-surface temperature over lakes (Y. Gao et al. 2019, manuscript submitted to Climate Dyn.). The simulation was output in 3-h intervals.
HRLDAS (Chen et al. 2007) was developed based on the Noah land surface model, and uses hourly precipitation and other meteorological variables, in uncoupled mode, to simulate the long-term evolution of the soil state. The goal of HRLDAS is to characterize the soil moisture/temperature and snow variability at small scales over large areas to provide improved initial land and vegetation conditions for the WRF/Noah coupled model. Hence, HRLDAS was configured to duplicate the CPM configuration to ensure consistency in the model grid spacing, physical configuration (e.g., terrain height), soil model, and parameters between the uncoupled soil initialization system and its coupled forecast counterpart. HRLDAS has been applied previously to the TP region to investigate land surface modeling uncertainties (Gao et al. 2015c; Zhang et al. 2016; Li et al. 2018).
To investigate the snow-cover conditions, the HRLDAS run spanned a cold-season period from 1 October 2013 to 31 May 2014. To reach a state of equilibrium, a spinup time of 3–11 years was necessary over the TP depending on physical-scheme options used in land models (Chen et al. 2014; Gao et al. 2015c; Zhang et al. 2016; Jiang et al. 2018). For this, we used China Meteorological Forcing Dataset (CMFD) data from 2012 and reran 5 years for spinup, as in previous studies over the TP. The advantage of HRLDAS is its use of high-resolution land-use and soil texture maps, and the same grid cells as WRF. As a result, it is able to capture finescale heterogeneity at the surface and in the soil. In this study, the Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type product, MCD12Q1, from 2013–14, with a resolution of 500 m, and the Chinese Dataset of Soil Hydraulic Parameters Using the Pedotransfer Functions for Land Surface Modeling is used to specify the soil texture (Shangguan et al. 2014). This dataset is in 0.0083° resolution and was generated based on the 1:1 million China soil map from the second general soil survey and 8595 sites of observations of soil profiles. The monthly mean leaf area index (LAI) for the TP, averaged over 2012–15 from the 0.1° MODIS data, was used.
1) Forcing for HRLDAS (except precipitation)
Seven atmospheric forcing conditions at the desired resolution were used to drive the Noah land surface model: precipitation, downward solar radiation, downward longwave radiation, near-surface temperature, humidity, wind, and surface pressure. Several gridded atmospheric datasets including all these variables and meeting the requirements for land surface modeling over the TP were generated; for instance, the CMFD (Chen et al. 2011; He et al. 2020), the China Meteorological Administration Land Data Assimilation System (CLDAS) (Shi et al. 2011; Jia 2013), and the forcing products of the Beijing Normal University (BNU) (Li et al. 2014; Shangguan et al. 2014). CMFD covers the period of 1979 onward, on a 0.1° × 0.1° grid with a temporal resolution of 3 h. It was generated based on Princeton forcing data (for the period 1979–2008) (Sheffield et al. 2006) and Global Land Data Assimilation System data (from 2009 onward), CMA station data (wind, air temperature, relative humidity, surface pressure, precipitation, and sunshine duration), and TRMM 3B42 satellite precipitation estimates. CLDAS is an hourly dataset with a spatial resolution of 0.0625°. It was generated by combining ECMWF forecast products with regional and national automatic surface meteorological observations (wind, air temperature, humidity, surface pressure, and precipitation) over China. BNU covers the period 1958–2010, on a 5-km grid, with a temporal resolution of 3 h. The spline and residual correction approach was used to merge observational, reanalysis, and remote sensing data for mapping near-surface air temperature, pressure, humidity, precipitation, wind speed, and shortwave and longwave radiation.
Numerous comparisons have been conducted to evaluate these forcing datasets in land surface modeling. Each has its own strengths and weaknesses, and no single product is consistently better than others. For instance, Gao et al. (2015c) compared the CMFD and BNU forcing over the central TP and found BNU/CMFD slightly overestimated the longwave radiation/surface pressure when compared with in situ field observations (Fig. 2 in Gao et al. 2015c). Yang et al. (2017) evaluated the precipitation and shortwave radiation in CMFD and CLDAS, as well as the Global Land Data Assimilation System (GLDAS), over mainland China. They found that the three forcing datasets all captured the key features of the spatial distribution in shortwave radiation, but the estimates from CLDAS and GLDAS were systematically higher than those from CMFD over most of mainland China. Xie et al. (2017) evaluated the CMFD and CLDAS forcing variables in a simulation of blowing snow. They found that the warm biases of CLDAS associated with low temperatures over the TP may have amplified the resistance of snow to drifting, and the humid air likely caused a decrease in snowdrift sublimation. They pointed out that the uncertainties in CLDAS-2 should be paid appropriate attention (Xie et al. 2017). Based on the abovementioned evaluations and characteristics of snow simulation, CMFD was selected as the basis of forcing. However, despite the slightly better performance of CMFD relative to the other datasets, it should be kept in mind that uncertainty using this forcing still exists.
2) Alternative precipitation datasets
Many available precipitation datasets have been evaluated and intercompared for the TP or subregions of the TP. Alazzy et al. (2017) found that CMORPH-CRT and TRMM 3B42 performed better than PERSIANN-CDR and TRMM 3B42RT when reproducing precipitation. For streamflow simulations over the Ganzi River basin, TRMM 3B42 has exhibited superior performance, and even outperformed simulations driven by gauge data during the validation period. Furthermore, the same results were found for the Mekong River basin (Chen et al. 2018) and at the plateau scale (Tong et al. 2014). TRMM 3B42 and CMORPH have demonstrated overall better performance than PERSIANN and the real-time version of TRMM in terms of correlation coefficients and RMSE (Gao and Liu 2013). Wu and Zhai (2013) found CMORPH and TRMM 3B42 to be more reliable over relatively flat terrain than over the complex terrain of the TP. They argued that TRMM 3B42 provides better detail of the regional gamma distributions of daily precipitation amount than CMORPH. Hu et al. (2018) found that all gridded datasets underestimated the observed precipitation over central Asia, especially in mountainous areas. In summary, it is generally agreed that TRMM 3B42 and CMORPH outperform other products for the TP region. Taking the data performance, accessibility, and temporal coverages into account, CMORPH and TRMM 3B42 were chosen as alternative datasets to CMFD.
3) Model evaluation datasets
We used in situ precipitation observations and satellite-derived gridded snow-cover fraction datasets to evaluate the WRF-simulated precipitation and HRLDAS-simulated snow-cover fraction. The in situ observations of precipitation were from the National Meteorological Information Center (http://data.cma.cn/site/index.html) of the CMA. There are 83 stations located on the TP (Fig. 1a), mostly in the central to eastern TP. Nine of them are located around the Brahmaputra Grand Canyon (Fig. 1b).
Five cloud-free snow-cover fraction datasets were evaluated and compared with the CMA in situ snow records after applying a four-step cloud-removal procedure to reduce cloud-cover-day percentage (Jiang et al. 2019). The cloud-free MODIS and Fengyun-3B snow-cover data produced a better annual mean and changes in recent decades when compared against other datasets. Here, the MODIS snow-cover fraction was used to evaluate the HRLDAS simulations since it possesses higher resolution and covers a longer period than Fengyun-3B. MODIS is a 36-channel visible to thermal-infrared sensor of the NASA Earth Observation System (Hall et al. 2002). The daily 500-m snow-cover products (version 6) of the two MODIS sensors (Terra and Aqua) were combined first to remove the cloud for the period 2013–14. Then, the resulting cloud-free MODIS data were used to evaluate the snow-cover fraction simulation. The snow season of October 2013–May 2014 was analyzed.
Comparisons between models and observations or between different resolution gridded datasets are always challenging, particularly for a region with complex terrain like the Tibet. Apart from the surface air temperature, which is elevation dependent, precipitation could not be interpolated as a function of elevation. We first started comparing precipitation datasets at their native grid cells (Figs. S2–S4), but this method may not be fair for coarse-resolution datasets. Chen and Dai (2017) claimed precipitation frequency and intensity are highly sensitive to the spatial and temporal resolutions of the input data and highlight the importance of the need to have similar resolution in comparing between two datasets, or between observations and models, or among different models. Therefore, all the precipitation datasets were interpolated into the coarsest 0.25° grid using the area averaging method that accounts for fractional contributions of the input high-resolution grid points to low-resolution grid point.
Compared to DDM, a more significant improvement in simulated precipitation amount was found in CPM on the native nearest grid cell (Figs. S2 and S3) than 0.25 grid cells (Figs. 2–4), especially in the simulated precipitation frequency. These results indicate 1) CPM shows moderate improvements over DDM even with CPM precipitation regridded to coarse resolution and 2) the largest benefits of CPM over DDM was found on small scales. This is consistent with findings over other mountain regions such as the European Alps and the Rocky Mountains (Prein et al. 2013a, 2015). Since both evaluation methods show consistent improvement in CPM over DDM, comparison on the 0.25° grid cells are used in the following text and the evaluation on native grid cells are presented in supplementary.
To compare with station data, we also used nearest grid point, four-point, and nine-point averages over nearest grids to stations and interpolated to the location of stations, but results show only slight differences among those three interpolation methods. Therefore, the nearest neighbor approach was used in the comparison with observations at 83 CMA stations. The precipitation evaluation was conducted not only on the precipitation amount, but also on precipitation frequency and intensity. Hourly datasets were then accumulated to daily values.
In the HRLDAS simulations, all of the forcing conditions were interpolated into a 4-km resolution, which was the same as that of the CPM grid cells. The MODIS snow-cover fraction was also remapped to the same 4-km HRLDAS grid to facilitate the evaluation of HRLDAS-simulated snow-cover fractions with the MODIS data. The absolute error (Bias), root-mean-square error (RMSE), standard deviation (SD), and pattern correlation (CORR) were used in the precipitation and snow-cover evaluations. The Taylor diagram (Taylor 2001), where SD, CORR, and centered RMSE can be shown on the same plot, was used to judge the performance of the snow-cover fraction simulation driven by the different precipitation datasets. The variations of snow-cover fractions with elevation and the probability distribution functions of the 4-km simulated snow-cover fraction were compared with those from MODIS.
a. Evaluation of precipitation
Figure 2 compares the precipitation amount, frequency, and intensity of the CMFD, CMORPH, TRMM, DDM, and CPM simulations with those from in situ observations obtained from 83 CMA stations. As expected, CMFD had the lowest difference in precipitation amount when compared with the CMA station records, because data from those CMA stations are ingested directly into CMFD. CMORPH ranked as the next best. TRMM significantly underestimated the precipitation; however, the DDM and CPM overestimated it. The differences in the precipitation amount were mainly caused by differences in the precipitation frequency. TRMM underestimated both the frequency and intensity, which led to a substantial underestimation of the precipitation amount. CMORPH and the DDM and CPM all overestimated the precipitation frequency and amount. The underestimation of intensity by CMORPH offset the overestimation of frequency, which resulted in a lower overestimation of precipitation amount when compared with the DDM and CPM results. When compared with the DDM, the CPM overestimated the intensity by more and the frequency by less, resulting in a relatively smaller overestimation of the precipitation amount.
Figure 3 displays the station distribution of differences for the five precipitation datasets compared against the CMA observations from 83 stations. With the lowest TP-average differences, CMFD yielded mild absolute differences of less than 100 mm at all stations (Fig. 3a). CMORPH followed CMFD, showing relatively small differences at most stations. TRMM underestimated the annual precipitation for most stations, especially for stations in the central-eastern TP. The DDM and CPM underestimated the precipitation for the headwaters of the Yangtze and Yellow River basins, but overestimated it in the north and south TP. The largest overestimation was for the southeastern TP. Therefore, it was this substantial overestimation for the southeastern TP that dominated the overall overestimation averaged over the TP (Fig. 2).
The southeastern TP possesses a sharper gradient in terrain elevation height than the northern TP. The Brahmaputra Grand Canyon is a passage through which abundant moisture can be transported to the inner TP. As shown in Fig. 1b, the topography rises steeply from 1000 m up to 6000 m in the Brahmaputra Grand Canyon. The nine CMA stations in this region are all located in valleys with relatively low elevations. The elevation difference between station and the nearest 0.25° grid cells averaged over 83 CMA stations was around 500 m, but increased to over 800 m averaged for nine stations in the Brahmaputra Grand Canyon (Fig. 1c). The differences in precipitation amount, frequency, and intensity for the nine stations in the Brahmaputra Grand Canyon as compared with 83 stations in the entire TP are shown in Fig. 4. As shown, the differences for these nine stations were much higher than the average difference for the 83 stations for the whole TP. CPS produced the largest differences between the 9 stations in the Brahmaputra Grand Canyon and the 83 stations in the entire TP; the 9-station-average bias was 7 times higher than the 83-station average. This increased difference in the Brahmaputra Grand Canyon when compared with elsewhere on the TP appeared not only in the precipitation frequency but also its intensity, except for CMORPH.
The SD of elevations at finer resolution in coarse-resolution grid cells denotes the topographic fluctuation of the given coarse grid. Here, SDs of the elevations of the 4-km simulation within 0.4° grid cells were calculated, and the results showed that SDs range from less than 100 m to greater than 900 m. Figure 5 compares the precipitation among the gridded datasets as well as the in situ observations for the SD of each elevation at intervals of 100 m. The numbers above the bars are the stations numbers at each interval. As shown, most stations were located at intervals of 200–500 m. Observationally, there was a gradual increase along with the SD of elevation. For instance, at intervals less than 100 m, observations showed an average precipitation of 180 mm averaged over two observational records. The station-average precipitation varied slightly, but remained generally consistent at intervals of 100–700 m. The mean and spread of precipitation (ranging between the 12.5 and 87.5 percentiles) increased slightly with the SD of elevation. The precipitation amount of the 87.5 percentile increased faster than the mean. At SD intervals of greater than 900 m, observed precipitation rose to around 500 mm. There was only one station result at the SD interval of greater than 900 m, meaning it was not possible to determine the station spread from the observation. However, this single station recorded a precipitation amount of 500 mm, which was far more than the average precipitation at other intervals. Hence, the observed precipitation increased with the SD of elevation.
CMFD and CMORPH generally reflected the variation of observed precipitation at SD intervals of smaller than 700 m, but greatly underestimated it at intervals greater than 900 m, while TRMM completely underestimated precipitation at every SD interval. The DDM and CPM generally reflected the observation at SD intervals less than 500 m, but overestimated it at intervals greater than 500 m. In particular, huge differences were revealed at intervals greater than 900 m. The DDM produced results that were closest to the observation among the five datasets. The CPM simulated a relatively higher precipitation amount than the DDM. However, the other datasets greatly underestimated the precipitation at the highest SD interval. The significant differences in precipitation between the merged gridded datasets and the WRF simulations hint that there is considerable uncertainty in the simulation of precipitation. The physics-processes-based WRF simulation might be a better representation of the precipitation in areas with large topographic fluctuations when compared with other datasets generated through station interpolation or merging.
It is noteworthy that observations in high-altitude mountain areas have often been underestimated in previous studies. Goodison et al. (1998) reported an underestimation of solid-precipitation measurements because of wind loss. Yang et al. (2005) suggested that the underestimation because of wind could be very large (80%–120%) in winter months and relatively small (10%) in the summer season over high-latitude regions (north of 45°N). In China, Ye et al. (2004) also claimed an underestimation of 6%–62% in the annual mean precipitation averaged over 710 standard meteorological stations. Because observations are scarce and there is tremendous uncertainty in measurements of the steep and high-altitude mountainous region, a few studies have derived the precipitation at high elevations using other variables, such as NDVI (Immerzeel et al. 2015). They found that the high-altitude precipitation was far greater than what was observed at valley stations or estimated by gridded precipitation products. Therefore, in the following analysis, we inversely infer the high-altitude precipitation over the TP by using the snow-cover fraction.
b. Evaluation of the snow-cover fraction simulations
Snow over the TP starts to accumulate in October and melts in April of the following year. Figure 6a shows the snow-cover fraction from the MODIS data in the evaluation period. High snow-cover fractions were located over the Kala Kunlun Mountains in the west, the Himalaya in the southwest, the Hengduan Mountains in the southeast, and the Tanggula Mountains in the central TP. Figure 6b shows that the distribution of the HRLDAS-simulated snow-cover fraction driven by CMFD precipitation was biased. Specifically, there was an underestimation in the central and eastern TP and an overestimation in the western TP. The largest underestimation was detected in the Brahmaputra Grand Canyon, even though the gridded precipitation input showed a high level of agreement with station observations below 4800 m. The same pattern of bias has been found when the surface air temperature is adjusted using the lapse rate of the neighboring 64 grid cells (Y. Jiang et al. 2020, manuscript submitted to J. Geophys. Res. Atmos.), indicating that errors in gridded air-temperature forcing data due to the elevation differences with station altitudes are not the major problem. Thus, we turn our attention to the precipitation data.
Over complex terrains, snow cover exhibits a high spatial heterogeneity (Gerber et al. 2017). On the regional scale, the variability of snow accumulation is mainly determined by the interaction of mountain ranges with the large-scale atmospheric circulation causing orographic precipitation (Stoelinga et al. 2013). To investigate their finescale features, the observed and simulated snow-cover fraction driven by the five precipitation datasets at the great turn of the Brahmaputra Grand Canyon is magnified in Fig. 7. The biases, RMSE, and CORRs between the MODIS snow-cover fraction and the simulations are listed in Table 1. Apparently, the simulation driven by CMFD missed the snow-cover fraction in the north (Fig. 7b). It underestimated the snow-cover fraction in this small region by 26.8%. The simulations driven by CMORPH and TRMM followed the same pattern of snow-cover fraction (Figs. 7c,d), with underestimations of 23.8% and 18.9%, respectively. In contrast, the simulations driven by the DDM and CPM reproduced the most snow-cover fraction distributions (Figs. 7e,f). The RMSE was reduced from the highest value of 39.2% in CMFD to almost half (22.5%) in the CPM-driven run. The DDM slightly overestimated the snow-cover fraction, whereas the CPM slightly underestimated it. The CORR with the MODIS dataset reached 0.79 and 0.80 in the DDM- and CPM-precipitation-driven runs, which was higher than for the runs driven by the three merged precipitation datasets.
The Taylor diagram (Fig. 8) summarizes the performances of the five precipitation datasets in driving the simulation of snow-cover fraction. When compared with the MODIS snow-cover fraction, HRLDAS significantly underestimated the snow cover when driven by CMFD, CMORPH, and TRMM precipitation. The difference in the standard deviation in the CMFD-driven run was the largest, followed by TRMM and CMORPH. In contrast, the simulations using the DDM and CPM precipitation captured the magnitude of the MODIS-observed snow-cover fractions reasonably well. In particular, the CPM-driven run produced the smallest RMSE. The CORR with MODIS was also higher in the DDM- and CPM-driven runs when compared with those driven by the three merged precipitation datasets.
Surface air temperature usually decreases with elevation, and snow-cover fraction increases with elevation height. However, the increase in snow-cover fraction with elevation is nonlinear. Figure 9 presents scatterplots of the snow-cover fraction of MODIS against elevation as well the HRLDAS simulations driven by the five precipitation datasets. Unsurprisingly, MODIS showed an increased snow-cover fraction with elevation height: less than 20% below 2500 m, 0%–60% between 2500 and 3500 m, an evenly distributed 0%–80% between 3500 and 4500 m, and greater than 80% above 5000 m. The median snow-cover fraction (50%) in MODIS expanded from the altitude of 3000–4800 m. Close examination of the median line of snow-cover fraction showed an inverse tangent function with elevation height (Fig. 9a).
The CMFD-, TRMM-, and CMORPH-driven runs showed a snow-cover fraction of almost zero below 3000 m, and a snow-cover fraction below 20% between 3000 and 4000 m. Most grid cells showed a snow-cover fraction below 30% above 4000 m, and only a few grids showed over 80% above 5000 m. On average, there were two major differences when compared against MODIS: 1) three runs missed the high snow-cover fractions below 3500 m and 2) they did not capture the very high snow-cover fraction (>80%) between 3500 and 5000 m (Figs. 9a,c,d). These differences were more pronounced in the CMFD- and TRMM-driven runs (Figs. 9b,c). The CMORPH-driven run, meanwhile, displayed a slightly higher snow-cover fraction between 3500 and 4500 m when compared with the CMFD- and TRMM-driven runs, but the magnitude of the snow-cover fraction was still much lower than that of MODIS (Fig. 9d). The median lines of snow-cover fraction against elevation in the runs driven by the three merged gridded precipitation datasets exhibited a power exponential function, rather than the inverse tangent function in MODIS.
The DDM- and CPM-driven runs, however, displayed better agreement with MODIS than those of the three merged precipitation datasets in terms of the snow-cover fraction distribution against elevation. The median lines of snow-cover fraction against elevation displayed the same inverse tangent function as with MODIS in both simulations. In particular, the CPM-driven run better captured the high snow-cover fractions between 3500 and 5000 m (Fig. 9f). The DDM-driven run slightly underestimated the snow-cover fraction between 2500 and 3500 m, but overestimated it above 3500 m, which formed a relatively narrower width in the snow-cover fraction distribution against elevation when compared with the CPM-driven run and MODIS dataset (Fig. 9e).
In summary, the HRLDAS runs driven by the three merged precipitation datasets substantially underestimated the snow-cover fraction over the TP. They failed to capture the high snow-cover fractions below 3500 m and presented low snow-cover fractions above 4000 m. The DDM showed a MODIS-like snow-cover fraction distribution against elevation, but underestimated the low snow-cover fraction and showed less spread. The CPM-driven run simulated the most MODIS-like snow-cover fraction distribution against elevation among the five runs. It not only better captured the high snow-cover fraction at high elevations, but also better exhibited the spread of snow cover between 3000 and 5000 m.
The advantages and disadvantages of the snow-cover fraction simulations driven by the five precipitation datasets compared to MODIS were also clear from the probability distribution functions of the snow-cover fraction (Fig. 10). MODIS observed the TP to be 43% snow-free, 37% as having a snow-cover fraction of 20%–70%, and 20% with a snow-cover fraction over 80% (snow-cover fractions of less than 15% are treated as snow-free in MODIS). Clearly, substantial spread existed in the snow-free area among the five simulations. The runs driven by the three merged precipitation datasets significantly overestimated the snow-free coverage. The CMFD- and TRMM-driven runs simulated a snow-free coverage of 74% and 70% over the TP, while the CMORPH-driven run simulated an even larger snow-free coverage (86%). They all substantially underestimated the snow-cover fraction of over 20%. The DDM-driven run, however, underestimated the snow-free coverage, at 33% of the area of the TP, while the CPM-driven simulated a 51% snow-free coverage, which was the closest among the five simulations to the MODIS coverage. For the snow-cover fraction of over 20%, the DDM- and CPM-driven runs generally reflected the variation of MODIS, except for a slight overestimation over 90%. DDM captures the PDF for 50%–90% SCF better than CPM, however, overestimates for SCF over 90%. On average, the DDM- and CPM-driven runs better captured the probability distribution functions of the snow-cover fraction, especially the high snow-cover fractions, compared with the runs driven by the three merged precipitation datasets. The CPM-driven run simulated the snow-free coverage even better than the DDM-driven run for the snow free and high SCF areas.
4. Conclusions and discussion
Accurate precipitation is critical for understanding land surface–atmosphere interactions and regional hydrological cycles. However, substantial uncertainties exist in available precipitation datasets in high-altitude mountain regions. On one hand, these regions are far beyond the reach of human accessibility and equipment because of the rugged topography, leading to a scarcity of observations. On the other hand, most of stations are located in valleys, leading to a poor representativeness of the stations. Moreover, significant errors exist in the gauge records due to the windy and cold environment at high altitude or latitude regions. All of these factors hamper our ability to produce high-quality gridded precipitation datasets.
In this study, gridded precipitation obtained from CPM, DDM, and three merged gridded precipitation datasets for 1 October 2013–31 May 2014 were evaluated and compared with available station records over the TP, with a particular focus in the Brahmaputra Grand Canyon with its steep topography. To further evaluate the precipitation at high elevations of above 4800 m where there are no observation stations, the HRLDAS was applied, driven by those five precipitation datasets. The simulated snow-cover fractions were intercompared and evaluated against a cloud-free MODIS snow-cover fraction dataset. The following conclusions were obtained:
The five precipitation datasets exhibited large uncertainties over the TP compared with the CMA observations, with the uncertainty in the precipitation frequency being larger than that in the precipitation intensity. Compared to station observations, TRMM greatly underestimated the precipitation amount; DDM and CPM generally overestimated the precipitation amount, especially in the Brahmaputra Grand Canyon. Note that the overestimation might not be “true” because there is a substantial undercatchment in winter precipitation measurements when using weighing type gauges (Yang et al. 1998; Serreze et al. 2001; Rasmussen et al. 2011; Lundquist et al. 2019).
The HRLDAS runs driven by the observation-based precipitation datasets substantially underestimated the snow-cover fractions in the central and eastern TP, particularly in the Brahmaputra Grand Canyon, although the observational precipitation datasets produced smaller biases than the WRF simulation compared with the CMA station observations. The runs driven by the three merged precipitation datasets substantially overestimated the snow-free area, but greatly underestimated the high snow-cover fraction above 3500 m when compared with the MODIS data.
The simulated snow-cover fraction driven by the DDM and CPM precipitation agreed with the MODIS snow-cover fraction much better than the runs driven by the three merged precipitation datasets. Moreover, the CPM-driven run produced the lowest biases and highest pattern correlations in snow-cover fraction among the five precipitation datasets compared with the MODIS snow-cover fraction. Not only does it capture the high snow-cover fraction above 3500 m better than the DDM, but also the snow-free areal coverage.
The above findings imply a significant underestimation of precipitation over 4800 m in currently available observational precipitation products, and the high-altitude precipitation amount in the Brahmaputra Grand Canyon might be twice what was previously thought. Underestimations in gridded precipitation datasets for high mountains are a common problem. Despite the best efforts, measurements of mountain precipitation are still limited by many factors such as blocked radar signals, spatial details too fine for satellites or sparse gauge networks, and precipitation gauges either capped by snow or missing snow blown over the top of the orifice (Lundquist et al. 2019). Immerzeel et al. (2015) claimed that the amount of precipitation in the upper Indus Basin is far above what was observed at valley stations or estimated by gridded precipitation products. Therefore, one needs to be cautious when evaluating models against precipitation datasets in high mountains.
The added value of WRF CPM simulated precipitation in high mountains is consistent with Lundquist et al. (2019) and He et al. (2019), who suggested that modeled precipitation has crossed a threshold in range-wide accuracy relative to observation-based precipitation datasets in complex terrain. Realizing this and fully taking the advantage of finescale modeling capabilities is also beneficial to improve observational networks and our abilities to understand, model, and forecast our mountain water supplies (Lundquist et al. 2019).
Despite these exciting findings, uncertainties in physics parameterization such as microphysics still remain. Liu et al. (2011) claimed the double-moment schemes are apparently superior to other schemes in snowfall amount simulation. Hughes et al. (2017) overall agree with Liu et al. (2011) but also pointed out the single-moment microphysics tend to better simulates the windward/leeside precipitation contrast than the double-moment ones. Van Weverberg et al. (2014) found that precipitation is similarly simulated in one-moment against two-moment microphysics schemes over Belgium. Morrison et al. (2012) and Rosenfeld et al. (2014) considered that two-moment schemes are highly tunable and many depicted processes are not well understood. Prein et al. (2015) pointed out that sensitivities related to cloud microphysical parameterizations might be even larger in CPM modeling due to a mesh refinement at kilometer scales. Convection-permitting modeling offers the advantage of improving the representation of finescale variations in orography and in land surface characteristics, which is especially beneficial in mountainous regions and in areas with heterogeneous land surfaces. The Noah land model has known biases in simulating snow in complex terrains, which has been somewhat addressed by the addition of multilayer snow in the Noah-MP model (Niu et al. 2011; Pavelsky et al. 2011; Letcher and Minder 2015). Further investigations are necessary to understand impact of these parameterizations and to select optimal CPM physics schemes.
Moreover, this study focused on the assessment of precipitation data impact on land surface modeling, and did not include the impacts of the interactions between precipitation and other forcing variables. Finally, limited by the available computational resources, only a relatively short time period (i.e., 8 months) of CPM simulation was conducted in this study. This is a pilot study using a wide range of datasets for model evaluation in a region that is characterized by observation sparsity and large observational uncertainties, and provides an alternative strategy for evaluating finescale modeled precipitation in the high mountains. Longer simulation period will be beneficial in improving our understanding of the added values of CPM for simulations under current and future climatic conditions. The recent CORDEX Flagship Pilot study, focusing on CPM modeling in the TP (http://www.cordex.org/2018/01/26/endorsed-flagship-pilot-studies/), is expected to yield interesting results.
We appreciate the free access of the Terra and Aqua MODIS snow cover data (MOD10A1 and MYD10A1) and the open access for the atmospheric forcing data (http://www.tpedatabase.cn/portal/MetaDataInfo.jsp?MetaDataId=249369). This research has been supported by the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (2019QZKK010314), the Strategic Priority Research Program of Chinese Academy of Sciences (XDA2006010202), the National Natural Science Foundation of China (91537211), and the NCAR Water System. Y. Gao is involved in the CORDEX-FPS-CPTP.
Denotes content that is immediately available upon publication as open access.