**1. Introduction**

The Tibetan Plateau (TP) is often referred to as "the roof of the world" and "the Asian water tower" because of the considerable impacts it exerts on the regional and global climate and water cycle. As the highest and widest plateau in the world, its complex and mountainous terrain is not only a barrier to the westerly belt at the same latitude but also strengthens the Indian monsoon through strong dynamic and thermodynamic effects, which accelerates the large-scale atmospheric circulation [1,2]. In addition, the TP possesses abundant water resources and is home to the sources of many of Asia's major rivers such as the Yangtze, Yellow, Tarim, Indus, Ganges, and Mekong [3,4]. The TP has attracted continuous attention amongst scientists owing to its sensitivity to climate change and its tendency to warm faster than the global average during the past few decades [5,6]. This enhanced warming has led to significant glacial and snow melt, permafrost degradation, and increasing precipitation [7,8], and these warming-induced changes also have marked effects on the processes of terrestrial evapotranspiration (ET) over the plateau, which accelerates the increase in ET [9].

As a critical component of water cycle, ET connects the land, ocean, and atmosphere; transports 60% of the world's precipitation over land [10]; affects the long-term evolution of vegetation [11,12]; and plays an important role in regulating regional droughts and floods [13]. As for the ET over the TP, research indicates that the mean ET ranges between

**Citation:** Dan, J.; Gao, Y.; Zhang, M. Detecting and Attributing Evapotranspiration Deviations Using Dynamical Downscaling and Convection-Permitting Modeling over the Tibetan Plateau. *Water* **2021**, *13*, 2096. https://doi.org/10.3390/ w13152096

Academic Editor: Gordon Huang

Received: 27 May 2021 Accepted: 27 July 2021 Published: 30 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

250 and 400 mm/year, is lowest in the northwestern TP and highest in the southeastern TP [9,14], and is limited by water in dry areas and by energy in wet areas at the annual scale [15]. Moreover, several studies have focused on the interannual variability of ET in the TP. For example, Zhang et al. estimated the ET in 16 catchments across the plateau and pointed out the rate of increase was 7 mm/(10 years) during 1966–2000 [16]. Li et al. investigated the ET of the Yellow, Yangtze, Qiangtang, and Qaidam basins from 1983 to 2006 and also found the ET showed an upward trend for all seasons [17]. Yin et al. pointed out that the increasing ET trend over the last 30 years is linked with increased precipitation [18].

The terrestrial ET is closely linked with land-surface characteristics (e.g., soil and land-use types). The altitude of the TP is mostly over 4000 m above the sea level, the general tendency is higher in the northwest TP and lower in the southeast TP. Restricted by the terrain height, the land-surface types over the TP are diverse, including grassland, wetland, tundra, snow, and so on [19,20]. Previous studies have found that the multi-year mean ET of low-covered grassland is much lower than that of medium- and high-covered grassland [21,22]; additionally, that of wetland over the TP has been increasing in recent years and may be responsible for wetland degradation [23]. The ET in the permafrost region of the TP has also been found to be affected by the freeze–thaw cycle and presents obvious seasonal changes [24]. Moreover, a recent study indicated that lake evaporation in the south is 1171.9 mm during the ice-free season, which is higher than that of northern TP (1059.7 mm), and the evaporated water amount is about 51.7 km<sup>3</sup> per year when all plateau lakes are included [25].

However, due to the complex topography and harsh natural environmental conditions, observation sites are sparsely distributed, and the representativeness of station observations is limited. Thus, considerable uncertainty exists with respect to the ET over the various land-surface types of the TP.

Given the sparseness of site observations, remote sensing and numerical climate simulations have provided alternative solutions for understanding ET at large scales and causal attributions for the deviations of ET over the TP using the above two ways has been addressed in some studies. Remote sensing usually combined with the traditional methods such as the Penman–Monteith (P-M) method [26], or with the surface energy balance methods, for example, the Surface Energy Balance System (SEBS) [27,28], these methods have a sound physical basis, but the uncertainty in the parameterization of environmental stress factors in the equation usually bring deviations of ET [29–32]. As for the numerical simulations, two kinds of simulation can be employed to model ET: "off-line" and "on-line". The former uses a land-surface model or hydrological model driven by near-surface meteorological variables, which is widely used in ET responses to climate changes. The latter uses a climate model coupled with a land-surface model or hydrological model. Compared to the former, the latter not only involves the ET responses to climate changes but also the ET feedbacks to regional climate. Therefore, it is a more physically reasonable way to simulate the land–atmosphere interaction, and both global and regional climate models can be used in on-line simulations. Although the land-surface model and global climate models are often used for climate studies, the ET over the TP remain poorly modeled, and several studies point out that the deviations of ET over TP were from the deviation of the storage of summer soil water storage, negative bias of precipitation, and the overestimation of downward shortwave radiation flux [15,33]. Some studies also indicated that these imperfections in global climate models (GCMs) may be induced by the coarse resolution and parameterizations for large-scale processes [34,35]. Regional climate models (RCMs) have higher spatial resolution and better depiction of physics parameterizations for meso- and micro-processes compared with GCMs [36–38]. Thus, RCMs are used for dynamic downscaling of the coupled coarse grid GCM output data or the input reanalysis data [39,40], so they provide data with high spatial and temporal resolution, which are to some extent more suitable for heterogeneous regions. Nowadays, this is called dynamical downscaling modeling (DDM) and has been widely applied to study regional climate and

climate change in North America, Europe, Africa, and Asia [41–47]. As for the TP, Gao et al. showed that 30 km DDM can significantly reduce the biases in precipitation and give a more accurate net precipitation (precipitation minus ET) over the TP [48], this is in agreement with Lin et al., who indicated that the finer resolutions can greatly diminish the positive precipitation deviations over the TP, especially from 30 to 10 km [49]. However, because of the limitation of convection parameterization, the highest resolution of DDMs obtained over the TP can only reach 0.25◦ (≈28 km), which is still not high enough for the ET simulations over the TP [50,51].

Recently, the common calculation errors caused by the use of convection parameterization schemes were mentioned by some studies [52–56]. Meanwhile, Prein et al. and Ou et al. also pointed out that the convection parameterization schemes start to break down as deep convection starts to be resolved, explicitly when grid resolution becomes smaller than 10 km, and a so-called "gray zone" existed when grid spacing was between 10 and 4 km, at these scales, the individual convective cells cannot be resolved and may lead to insufficiently resolved deep convection [54,57]. Thus, convection-permitting modeling (CPM), which not only removes the convection parameterization scheme, but also greatly improves the data resolution and the representation of complex orography, is a more advanced method than DDM and has become a better choice for researchers studying the TP. Previous studies found that CPM usually gives more accurate results over the TP, especially for precipitation [49,58]. However, detailed comparison and evaluation of CPM-simulated ET over TP is rare. Therefore, in this study, we investigated (1) the performance of "on-line" simulations in ET over TP compared to remote sensing and "off-line" simulations; (2) whether dynamical downscaling (DDM and CPM) performs better in terms of the heterogeneity of ET—for instance, over different land-surface types—over the TP than global models; and (3) what causes the deviation of ET simulated by CPM and DDM over the TP, and whether it varies with the seasons.

The paper is structured as follows: Section 2 briefly introduces the Weather Research and Forecasting (WRF) model and model configurations, datasets, and methods. Section 3 evaluates the datasets mentioned above in terms of annual and seasonal means, as well as dominant land-use types, and then explores the factors contributing to the deviation in ET. Finally, the main conclusions and some further discussions are presented in Sections 4 and 5, respectively.

#### **2. Model, Datasets, and Methods**

## *2.1. Model and Configurations*

The Weather Research and Forecasting (WRF) model is a next-generation mesoscale numerical weather prediction system (http://www.wrf-model.org/index.php, accessed on 20 June 2021) [59], which was initially developed in the late 1990s and is still being improved today by scientists all over the world. WRF includes several dynamic cores, such as a fully mass- and scalar-conserving flux from the mass coordinate version, as well as many different physical parameterizations (e.g., microphysics, cumulus parameterization, planetary boundary layer, shortwave and longwave radiation, and land-surface models) [51,60,61].

The DDM run (outer domain) covered the entire Eurasian continent with a 28 km horizontal grid spacing, while the CPM run (inner domain) covered the entire TP with a 4 km horizontal grid spacing. The topography of the TP is shown in Figure 1a,b. Clearly, the 4 km resolution topography shows far more detail in terms of mountain tops and valleys, especially in the southern TP.

**Figure 1.** (**a**,**b**) The distribution of topography (unit: m), and (**c**,**d**) the dominant land-use categories over the Tibetan Plateau (TP). Panels (**a**,**c**) are from DDM (dynamical downscaling modeling), and panels (**b**,**d**) are from CPM (convectionpermitting modeling). (**e**) The land-use distribution (unit: %) of ten dominant categories. In (**c**–**e**), the numbers 1–10 stand for cropland, grassland, shrubland, mixed grassland/shrubland, forest, water bodies, wetland, barren or sparsely vege-**Figure 1.** (**a**,**b**) The distribution of topography (unit: m), and (**c**,**d**) the dominant land-use categories over the Tibetan Plateau (TP). Panels (**a**,**c**) are from DDM (dynamical downscaling modeling), and panels (**b**,**d**) are from CPM (convection-permitting modeling). (**e**) The land-use distribution (unit: %) of ten dominant categories. In (**c**–**e**), the numbers 1–10 stand for cropland, grassland, shrubland, mixed grassland/shrubland, forest, water bodies, wetland, barren or sparsely vegetated, tundra, snow or ice, respectively.

tated, tundra, snow or ice, respectively. Following our group's previous study on DDM and CPM [51], the NCAR (National Center for Atmospheric Research) CAM (Community Atmosphere Model) radiation scheme [62], the WSM6 (WRF Single-Moment 6-class microphysics scheme) [63], Yonsei University planetary boundary layer scheme [64], the Kain–Fritsch convection scheme [65] (used in the DDM run but not in the CPM run), and the Noah LSM (land-surface model) four-layer soil temperature and moisture model [66] were employed in this study too. The specific model configurations are shown in Table 1. We selected 2014 as the research period because the precipitation in this year was close to the climatological mean, which can help avoid the influence of certain extreme conditions and better test the ability of CPM and DDM to simulate ET under the climate mean state. The lateral boundary conditions and sea surface temperature (SST) were provided by the ERA-Interim reanalysis Following our group's previous study on DDM and CPM [51], the NCAR (National Center for Atmospheric Research) CAM (Community Atmosphere Model) radiation scheme [62], the WSM6 (WRF Single-Moment 6-class microphysics scheme) [63], Yonsei University planetary boundary layer scheme [64], the Kain–Fritsch convection scheme [65] (used in the DDM run but not in the CPM run), and the Noah LSM (land-surface model) four-layer soil temperature and moisture model [66] were employed in this study too. The specific model configurations are shown in Table 1. We selected 2014 as the research period because the precipitation in this year was close to the climatological mean, which can help avoid the influence of certain extreme conditions and better test the ability of CPM and DDM to simulate ET under the climate mean state. The lateral boundary conditions and sea surface temperature (SST) were provided by the ERA-Interim reanalysis dataset at 6 h intervals. For both DDM and CPM, the simulations were separated into two stages: the first was initialized at 0000 UTC 1 October 2013, ended at 2300 UTC 31 May 2014, and

archived in 3 h intervals; and the second was initialized at 0000 UTC 1 June 2014, ended at 2300 UTC 31 December 2014, and archived in 1 h intervals.

**Table 1.** The Weather Research and Forecasting (WRF) model configurations for dynamical downscaling modeling (DDM) and convection-permitting modeling (CPM).


#### *2.2. Datasets*

Besides the DDM and CPM simulations, a gridded remote sensing dataset (EB), a dataset from off-line simulations (GLDAS), and two large-scale reanalysis datasets (ERA-Interim and ERA5) were adopted for comparison.

The EB dataset (http://data.tpdc.ac.cn, accessed on 20 June 2020) was derived from satellite data and a surface energy balance method for the global land area. It includes daily and monthly datasets from 2000 to 2017, with a 5 km horizontal resolution. This dataset has been proven to be robust across a variety of land-cover types and performs well in providing spatial and temporal information on the water cycle and land–atmosphere interactions for the Chinese landmass [67]. Therefore, EB was used as the reference ET in this study.

GLDAS [68] simulates the ET from four land-surface models driven by a series of conventional and satellite-derived observations. GLDAS products have a spatial resolution of 0.25◦ × 0.25◦ and a temporal resolution of 3 h. For a long time, GLDAS has been the only gridded dataset for ET available over the TP [60,69] and has been found to perform well in terms of surface air temperature and precipitation forcing [70].

ERA-Interim and ERA5 are commonly used reanalysis datasets released by the ECMWF (the European Center for Medium-Range Weather Forecasts). ERA-Interim [71] covers the period 1979 to 2019 and has a horizontal resolution of 0.7◦ . It was produced with a 12 h 4D-Var (four-dimensional variational) data assimilation scheme and an IFS (Integrated Forecast System) release Cy31r2 forecast model. Previous studies have shown that it is the outstanding performer among all reanalysis datasets in describing the water cycle over TP [48,69].

ERA5 [72] is the successor to ERA-Interim, with a temporal resolution of 1 h and horizontal resolution of 0.1◦ . With a more sophisticated hybrid incremental 4D-Var system [73] and advanced forecast model (IFS release Cy41r2), ERA5 covers a longer period from 1950 to the present day and outputs more meteorological elements. Moreover, its radiative transfer model, land-surface model, and snow data assimilation are all different from those of ERA-Interim. These improvements may have advantages in exploring regions with complex terrain, such as the TP.

#### *2.3. Methods*

To intercompare the six datasets with different temporal and spatial resolutions, all of the datasets were first monthly averaged and then interpolated into the same 0.25◦ grid through the local area averaging method. The comparison was conducted in terms of the annual and seasonal means as well as seasonal variabilities. According to the typical seasonal variation of precipitation over the TP [2], the monthly data were further averaged to the monsoon season (May–September) and the non-monsoon season (October–April).

Since the calculation of ET is directly land-use dependent, we compared the dominant land-use types in the DDM and CPM simulations. There are 24 land-use types for the USGS (United States Geological Survey) categories in the WRF model (http://www2.mmm. ucar.edu/wrf/users/docs/user\_guide\_V3.8/AWUsersGuideV3.8.pdf, accessed on 20 June 2021) [74]. However, not all of them exist over the TP. To present the categories clearly, we reclassified these land-use types over the TP into 10 categories (cropland, grassland, shrubland, mixed shrubland/grassland, forest, water bodies, wetland, barren or sparsely vegetated, tundra, and snow or ice) (provided in Table 2). The distributions of these 10 categories are shown in Figure 1c–e. Considering that the EB data do not include the underlying surface of water bodies, all comparisons in this study do not consider this surface type. The distributions according to DDM and CPM are quite similar, with very slight differences in percentages for each category (Figure 1c–e). Therefore, the distribution of the 10 dominant categories of CPM can be used in the comparison among the 0.25◦ grid cells via the nearest neighbor algorithm. To investigate the contributions from each land-use category to the TP average, the ET in the same category was averaged and intercompared among the six datasets.

**Table 2.** Reclassified 10-category land-use categories.


In addition, to better judge whether the ET deviation from the other five datasets over the dominate land-use categories of TP is significant compared to EB, we also conducted a 99% significance test with the Monte-Carlo method [75], which does not require the data to be normally distributed. The Monte-Carlo method regards the two datasets with sample sizes of *n<sup>a</sup>* and *n<sup>b</sup>* as a method of randomly extracting *n<sup>a</sup>* samples from the total (*n<sup>a</sup>* + *nb*) samples, it calculates the absolute value, *D<sup>i</sup>* , of the sample mean difference from m resampled group as follows:

$$D\_i = |\overline{\mathfrak{X}\_{ai}} - \overline{\mathfrak{X}\_{bi}}| \tag{1}$$

where *xai*, *xbi* represent the averaged value of the two new datasets in each group; *i* = 1, 2, 3, . . . , *m*; *m* is re-sampling times, and *m* = 10,000 in this study. Then, the obtained *D<sup>i</sup>* is sorted from smallest to largest and compared with the original sequence mean difference, *Dori*. If the *Dori* satisfies the following condition:

$$D\_{ori} < D\_{\frac{\mu}{2}} \text{ or } D\_{ori} > D\_{\frac{1-\mu}{2}} \tag{2}$$

where *α* = 0.01 in this study, it is considered that the samples have passed the 99% Monte-Carlo significance test, and there is a significant difference between the two compared datasets.

Many factors play a key role in the simulation of ET, especially precipitation [76–78]. In addition, radiation provides the energy source for ET. However, the benchmark data EB is a satellite-derived dataset and has no other matching variables. Thus, it cannot be used as

an attribution reference criterion. In order to explore the factors influencing the deviations of simulated ET, first, EB was used as the standard to find the nearest and farthest datasets to it, and then the nearest one was used as the reference in the comparison and the farthest to explore the source of deviations in attribution analysis (in Section 3.3).

Then, to explore the influence of precipitation on deviations in ET, the deviations of ET (dET) of the selected datasets were separated into those from precipitation (dP) and all other types of deviation (dET-dP), which were obtained from

$$\mathbf{dA}\_{i} = A\_{si} - A\_{refi} \tag{3}$$

where *dA<sup>i</sup>* is an abstract symbol, which can be replaced by deviations from ET (mm/d), precipitation (mm/d), net radiation flux (w m−<sup>2</sup> ), or soil temperature (K) in this equation; *Asi* is the corresponding variable (ET, precipitation, net radiation flux, or soil temperature) of the dataset furthest from the ET value of the EB data (mm/d); *Are f i* is the corresponding variable (ET, precipitation, net radiation flux, or soil temperature) of the dataset closest to the ET value of the EB data (mm/d); *i* is the land-use type range from 1 to 10 without 6 (the water bodies).

To further explore the contributions from radiation and the land-surface model in the non-monsoon season, the dET-dP was further analyzed by comparing deviations from surface net radiation (dR) and soil temperature (dST), which can be calculated using Equation (1). Considering the different units, a relative contribution was defined by normalized deviations, and defined as

$$\mathbf{dB}\_{r\bar{j}} = \frac{|dB\_{\bar{j}}|}{|dB\_{\text{max}}|} \times 100\% \tag{4}$$

where d*B<sup>j</sup>* is an abstract symbol, which can be replaced by the deviations from net radiation flux (w m−<sup>2</sup> ) or soil temperature (K) in this equation; *dBrj* is the relative deviation (%) of the corresponding variable (net radiation flux or soil temperature); *dBmax* is the maximum deviation of the corresponding variable [net radiation flux (w m−<sup>2</sup> ) or soil temperature (K)] of all nine land-use types (without the water bodies); *j* is land-use type 9 (the tundra) or 10 (snow or ice).

The relative contributions from surface net radiation (dRr) and the land-surface model (dSTr - dRr) are represented by the magnitudes of the normalized deviations.

#### **3. Results**

#### *3.1. Annual and Seasonal Mean ET*

Figure 2 presents the annual mean distribution of ET in EB and its differences with those of the GLDAS, ERA-Interim, ERA5, DDM, and CPM datasets. EB exhibits a decrease in gradient from the southeastern TP to the northwestern TP. The daily ET reaches 2 mm d −1 in the southeastern TP and is almost down to zero in the northwestern TP (Figure 2a). The spatial distribution is generally in line with the pattern of precipitation over the TP [50,51]. Compared to EB, GLDAS shows a very similar annual mean spatial distribution (Figure 2b). ERA-Interim and ERA5 both significantly overestimate the ET almost for the entire TP (Figure 2c,d). As for the DDM and CPM, the simulations are very similar, both of them producing lower ET than EB in the central and southern regions of the plateau, with CPM being relatively lower than DDM (Figure 2e,f).

**Figure 2.** Distribution of the daily terrestrial evapotranspiration (ET) in 2014 in the (**a**) EB data (units: mm/d) and (**b**–**f**) the **Figure 2.** Distribution of the daily terrestrial evapotranspiration (ET) in 2014 in the (**a**) EB data (units: mm/d) and (**b**–**f**) the deviation from the EB distribution in the (**b**) GLDAS, (**c**) ERA-Interim, (**d**) ERA5, (**e**) DDM and (**f**) CPM data.

deviation from the EB distribution in the (**b**) GLDAS, (**c**) ERA-Interim, (**d**) ERA5, (**e**) DDM and (**f**) CPM data. Figure 3 shows the annual and seasonal mean ET for the six datasets in 2014 averaged over the TP. Averaged over the TP, EB estimates around 0.79 mm of moisture was evaporated per day in 2014. Broken down, this figure is 1.12 mm d−1 for the monsoon season, which accounts for most of the annual ET, and 0.52 mm d−1 for the non-monsoon season. GLDAS shows obvious deviation in both monsoon season and non-monsoon season in comparison with EB. ERA-Interim and ERA5 have the largest annual deviation, which was mainly from the monsoon season. The DDM and CPM are the closest datasets to EB Figure 3 shows the annual and seasonal mean ET for the six datasets in 2014 averaged over the TP. Averaged over the TP, EB estimates around 0.79 mm of moisture was evaporated per day in 2014. Broken down, this figure is 1.12 mm d−<sup>1</sup> for the monsoon season, which accounts for most of the annual ET, and 0.52 mm d−<sup>1</sup> for the non-monsoon season. GLDAS shows obvious deviation in both monsoon season and non-monsoon season in comparison with EB. ERA-Interim and ERA5 have the largest annual deviation, which was mainly from the monsoon season. The DDM and CPM are the closest datasets to EB in the monsoon season (Figure 3 and Table 3), while CPM is the only one that underestimates the ET.

in the monsoon season (Figure 3 and Table 3), while CPM is the only one that underesti-

mates the ET.

**Figure 3.** The annual and monsoon/non-monsoon seasonal mean ET in 2014 from the EB, GLDAS, ERA-Interim, ERA5, DDM, and CPM data averaged over the TP (units: mm/d). **Figure 3.** The annual and monsoon/non-monsoon seasonal mean ET in 2014 from the EB, GLDAS, ERA-Interim, ERA5, DDM, and CPM data averaged over the TP (units: mm/d).

**Table 3.** The mean ET from the EB, GLDAS, ERA-Interim, ERA5, DDM, and CPM data, and the standard deviation (SD), the root-mean-square error (RMSE), and spatial correlation coefficients (CORR) of GLDAS, ERA-Interim, ERA5, DDM, and CPM compared to EB. **Table 3.** The mean ET from the EB, GLDAS, ERA-Interim, ERA5, DDM, and CPM data, and the standard deviation (SD), the root-mean-square error (RMSE), and spatial correlation coefficients (CORR) of GLDAS, ERA-Interim, ERA5, DDM, and CPM compared to EB.


CORR Ann --- 0.94 0.95 0.91 0.95 0.93 Monsoon --- 0.94 0.92 0.91 0.95 0.94 Non-Monsoon --- 0.90 0.95 0.87 0.90 0.84 The seasonal mean ET of EB shows a great decrease in the non-monsoon season when compared to that in the monsoon season (Figure 3), which is mainly due to the large decrease in ET over the southern and central TP in the non-monsoon season in comparison with the monsoon season (Figure 4a,b). As for the GLDAS, it is relatively larger in the southeastern TP in the monsoon season and underestimated over the central and southern TP in the non-monsoon season (Figure 4c,d), which results in a balance in the annual mean. The overestimation of ERA-Interim and ERA5 exists almost across the whole TP but is largest in the northwestern TP in the monsoon season (Figure 4e,g). This is unsurprising given the known positive precipitation bias of ERA-Interim and ERA5 over the TP [79]. As for the non-monsoon season, ERA-Interim performs best compared to EB (Figure The seasonal mean ET of EB shows a great decrease in the non-monsoon season when compared to that in the monsoon season (Figure 3), which is mainly due to the large decrease in ET over the southern and central TP in the non-monsoon season in comparison with the monsoon season (Figure 4a,b). As for the GLDAS, it is relatively larger in the southeastern TP in the monsoon season and underestimated over the central and southern TP in the non-monsoon season (Figure 4c,d), which results in a balance in the annual mean. The overestimation of ERA-Interim and ERA5 exists almost across the whole TP but is largest in the northwestern TP in the monsoon season (Figure 4e,g). This is unsurprising given the known positive precipitation bias of ERA-Interim and ERA5 over the TP [79]. As for the non-monsoon season, ERA-Interim performs best compared to EB (Figure 4f and Table 3), while ERA5 shows an obvious underestimation especially for the southeastern TP (Figure 4h). DDM and CPM are the closest datasets to EB in the monsoon season (Figure 4i,k); however, in the non-monsoon season, DDM is significantly lower than EB for most of the TP but especially in the southeast (Figure 4j), while CPM ranks as the next best in comparison (Figure 4l and Table 3).

soon

4f and Table 3), while ERA5 shows an obvious underestimation especially for the southeastern TP (Figure 4h). DDM and CPM are the closest datasets to EB in the monsoon season (Figure 4i,k); however, in the non-monsoon season, DDM is significantly lower than EB for most of the TP but especially in the southeast (Figure 4j), while CPM ranks as the

next best in comparison (Figure 4l and Table 3).

**Figure 4.** Distributions of the daily ET from the (**a**,**b**) EB data (units: mm/d), and (**c**–**l**) the differences of the (**c**,**d**) GLDAS, (**e**,**f**) ERA-Interim, (**g**,**h**) ERA5, (**i**,**j**) DDM, and (**k**,**l**) CPM data compared to EB, averaged over the monsoon (left-hand panels) and non-monsoon (right-hand panels) sea-**Figure 4.** Distributions of the daily ET from the (**a**,**b**) EB data (units: mm/d), and (**c**–**l**) the differences of the (**c**,**d**) GLDAS, (**e**,**f**) ERA-Interim, (**g**,**h**) ERA5, (**i**,**j**) DDM, and (**k**,**l**) CPM data compared to EB, averaged over the monsoon (left-hand panels) and non-monsoon (right-hand panels) seasons.

#### sons. *3.2. ET over Dominant Land-Use Categories*

shown in Table 4.

*3.2. ET over Dominant Land-Use Categories*  Figure 5 presents box-and-whisker plots of the ET of nine land-use types based on the six datasets. The ET values averaged over the monsoon and non-monsoon seasons are given in Figures 6 and 7, respectively. In addition, we also used the Monte-Carlo test to Figure 5 presents box-and-whisker plots of the ET of nine land-use types based on the six datasets. The ET values averaged over the monsoon and non-monsoon seasons are given in Figures 6 and 7, respectively. In addition, we also used the Monte-Carlo test to detect whether the simulated ET is significantly different from EB, and the results are shown in Table 4.

detect whether the simulated ET is significantly different from EB, and the results are

**Figure 5.** Box-and-whisker plots of ET of nine land-use categories based on six datasets averaged in 2014: (**a**) cropland; (**b**) grassland; (**c**) shrubland; (**d**) mixed grassland/shrubland; (**e**) forest; (**f**) wetland; (**g**) barren or sparsely vegetated; (**h**) tundra; (**i**) snow or ice (units: mm/d). The top and bottom of each box are the 20th and 80th percentiles in grid cells with the same land-use type. The lines in each box represent the median; red dots represent the mean. **Figure 5.** Box-and-whisker plots of ET of nine land-use categories based on six datasets averaged in 2014: (**a**) cropland; (**b**) grassland; (**c**) shrubland; (**d**) mixed grassland/shrubland; (**e**) forest; (**f**) wetland; (**g**) barren or sparsely vegetated; (**h**) tundra; (**i**) snow or ice (units: mm/d). The top and bottom of each box are the 20th and 80th percentiles in grid cells with the same land-use type. The lines in each box represent the median; red dots represent the mean. *Water* **2021**, *13*, x FOR PEER REVIEW 12 of 21

**Figure 6.** Box-and-whisker plots of ET of nine land-use categories based on six datasets for the monsoon season in 2014: (**a**) cropland; (**b**) grassland; (**c**) shrubland; (**d**) mixed grassland/shrubland; (**e**) forest; (**f**) wetland; (**g**) barren or sparsely vegetated; (**h**) tundra; (**i**) snow or ice (units: mm/d). The top and bottom of each box are the 20th and 80th percentiles in grid cells with the same land-use type. The lines in each box represent the median; red dots represent the mean. **Figure 6.** Box-and-whisker plots of ET of nine land-use categories based on six datasets for the monsoon season in 2014: (**a**) cropland; (**b**) grassland; (**c**) shrubland; (**d**) mixed grassland/shrubland; (**e**) forest; (**f**) wetland; (**g**) barren or sparsely vegetated; (**h**) tundra; (**i**) snow or ice (units: mm/d). The top and bottom of each box are the 20th and 80th percentiles in grid cells with the same land-use type. The lines in each box represent the median; red dots represent the mean.

**Figure 7.** Box-and-whisker plots of ET of nine land-use categories based on six datasets for the non-monsoon season in 2014: (**a**) cropland; (**b**) grassland; (**c**) shrubland; (**d**) mixed grassland/shrubland; (**e**) forest; (**f**) wetland; (**g**) barren or sparsely vegetated; (**h**) tundra; (**i**) snow or ice (units: mm/d). The top and bottom of each box are the 20th and 80th percentiles in grid cells with the same land-use type. The lines in each box represent the median; red dots represent the mean.

**Figure 6.** Box-and-whisker plots of ET of nine land-use categories based on six datasets for the monsoon season in 2014:

grid cells with the same land-use type. The lines in each box represent the median; red dots represent the mean.

**Figure 7.** Box-and-whisker plots of ET of nine land-use categories based on six datasets for the non-monsoon season in 2014: (**a**) cropland; (**b**) grassland; (**c**) shrubland; (**d**) mixed grassland/shrubland; (**e**) forest; (**f**) wetland; (**g**) barren or sparsely vegetated; (**h**) tundra; (**i**) snow or ice (units: mm/d). The top and bottom of each box are the 20th and 80th percentiles in grid cells with the same land-use type. The lines in each box represent the median; red dots represent the mean. **Figure 7.** Box-and-whisker plots of ET of nine land-use categories based on six datasets for the non-monsoon season in 2014: (**a**) cropland; (**b**) grassland; (**c**) shrubland; (**d**) mixed grassland/shrubland; (**e**) forest; (**f**) wetland; (**g**) barren or sparsely vegetated; (**h**) tundra; (**i**) snow or ice (units: mm/d). The top and bottom of each box are the 20th and 80th percentiles in grid cells with the same land-use type. The lines in each box represent the median; red dots represent the mean.

Results from Figure 5 indicate that the ET of EB is generally larger over the vegetated ground of the TP. Cropland evaporates the most, with an annual mean value of 1.6 mm d−<sup>1</sup> , while barren or sparsely vegetated ground only has 0.3 mm d−<sup>1</sup> . The tundra and snow or ice surfaces present relatively high values for the annual mean. GLDAS shows relatively high consistency with EB in its annual mean over most land-use types, except grassland and tundra. Moreover, although the ET of GLDAS is also seen to be lower than that of EB over the snow/ice land surface, it is not significant (Table 4). The annual mean ET of ERA-Interim and ERA5 are relatively higher than those of EB, especially for the shrubland, mixed grassland/shrubland, forest, and wetland, while DDM and CPM produce lower annual mean ET over vegetated land.

Figure 6 presents the ET value of nine land-use types in the monsoon season. The ET of the EB dataset in the monsoon season is in line with its annual mean situation, the vegetated land evaporates generally more than the barren land. However, the performance of GLDAS in terms of the seasonal mean is unsatisfactory. In the monsoon season, clear overestimation exists for shrubland and grassland, which cover the largest area of the TP (Figure 1). The values provided by ERA-Interim and ERA5 are generally higher than those of other datasets, and this phenomenon exists and is significant for nearly all surfaces in the monsoon season (Table 4). The most obviously overestimated land-surface type of the reanalysis during the monsoon season is mixed grassland/shrubland, with the deviations almost reaching 1.1 mm d−<sup>1</sup> . Moreover, for DDM and CPM, they both perform much better than other datasets and are closer to the EB data over most surfaces in the monsoon season.


**Table 4.** The Monte-Carlo significance test of ET at confidence level of 99% over nine land-use categories of GLDAS, ERA-Interim, ERA5, DDM, and CPM compared to EB. \* Indicates that the value has passed the 99% significance test, which means the ET is significantly different from that of EB over the corresponding land-surface type.

> As for the non-monsoon season shown in Figure 7, the tundra and snow or ice surfaces of EB still present relatively high values. Apart from the monsoon season, GLDAS shows large underestimations over cropland in the non-monsoon season (Figure 7 and Table 4) and needs further study as well. ERA-Interim agrees well with EB for nearly all dominant land surfaces in comparison to the non-monsoon season, while the ET in ERA5 is generally lower than that of ERA-Interim, except for wetland and barren or sparsely vegetated land. The lower annual mean ET over vegetated land produced by DDM and CPM (shown in Figure 5) is mainly due to underestimation in the non-monsoon season, especially for DDM. As shown, DDM simulates considerably lower ET over cropland, grassland, shrubland, mixed shrubland/grassland, and forest, while CPM underestimates the ET over cropland, grassland as well as shrubland, which ultimately leads to the underestimation over the whole of the TP. Moreover, the underestimations for tundra and snow/ice in our simulations are significant, as Table 4 shows, and still need to be improved in future simulations.

#### *3.3. Attribution of ET Deviations*

The DDM- and CPM-simulated ET deviations compared to ERA-Interim and their counterparts from precipitation and other factors are presented in Figure 8. The solid lines represent a constant dET. Most points are located along the solid lines, which demonstrates that the ET deviation is equivalent for most land-cover categories, either in the monsoon or

effects.

non-monsoon season. The dashed lines combined with the solid lines are used to separate the ET deviation. Points above or below the angle between the dashed and solid lines indicate that the ET deviation is more a result of the precipitation deviation. Points to the left or right indicate that the ET deviation is more a result of the other effects. Meanwhile, the points outside the red dashed box indicate that the deviations from precipitation and other effects are larger than one SD, which means the deviations are significant and cannot be ignored. *Water* **2021**, *13*, x FOR PEER REVIEW 15 of 21

**Figure 8.** Deviations (monsoon: ERA-Interim minus DDM or CPM; non-monsoon: DDM or CPM minus ERA-Interim) in ET (units: mm/d) from precipitation and other factors for each land-use category in the (**a**,**b**) monsoon and (**c**,**d**) nonmonsoon seasons. The ordinate is the ET deviation from precipitation, the abscissa is the ET deviation from the other effects. The black solid line indicates 1:1 correspondence of the first and third quadrants, and the gray dashed line indicates 1:1 correspondence of the second and fourth quadrants. The red dashed box indicates the SD from precipitation and other **Figure 8.** Deviations (monsoon: ERA-Interim minus DDM or CPM; non-monsoon: DDM or CPM minus ERA-Interim) in ET (units: mm/d) from precipitation and other factors for each land-use category in the (**a**,**b**) monsoon and (**c**,**d**) non-monsoon seasons. The ordinate is the ET deviation from precipitation, the abscissa is the ET deviation from the other effects. The black solid line indicates 1:1 correspondence of the first and third quadrants, and the gray dashed line indicates 1:1 correspondence of the second and fourth quadrants. The red dashed box indicates the SD from precipitation and other effects.

In the monsoon season, DDM and CPM agree well with EB. However, ERA-Interim and ERA5 show a large difference. Considering that the boundary conditions and SST in our simulations came from ERA-Interim, exploring its deviation from the two analog values can also lead to a better understanding of where the WRF mode has improved. Thus, DDM and CPM were selected as the reference data to explore the contributing factors to the ET deviations of ERA-Interim in the monsoon season. As shown in Figure 8a,b, all land-use categories in ERA-Interim, except snow or ice, show clear positive deviations in In the monsoon season, DDM and CPM agree well with EB. However, ERA-Interim and ERA5 show a large difference. Considering that the boundary conditions and SST in our simulations came from ERA-Interim, exploring its deviation from the two analog values can also lead to a better understanding of where the WRF mode has improved. Thus, DDM and CPM were selected as the reference data to explore the contributing factors to the ET deviations of ERA-Interim in the monsoon season. As shown in Figure 8a,b, all land-use categories in ERA-Interim, except snow or ice, show clear positive deviations in precipitation in comparison with DDM, which indicates that the overestimation of ET

precipitation in comparison with DDM, which indicates that the overestimation of ET in ERA-Interim mainly comes from its large precipitation bias, especially for grassland. Im-

mances of DDM and CPM are not as good as those of the monsoon season. ERA-Interim shows the best performance and was, therefore, selected as the reference dataset in the non-monsoon season. The deviation from other factors dominates over most land surfaces of the TP, except cropland, grassland, and forest. As shown in Figure 8c,d, although most points are inside the red dashed box, the points for tundra and snow/ice show significant in ERA-Interim mainly comes from its large precipitation bias, especially for grassland. Importantly, the two WRF simulations greatly improve the result.

In the non-monsoon season, however, according to the previous analysis, the performances of DDM and CPM are not as good as those of the monsoon season. ERA-Interim shows the best performance and was, therefore, selected as the reference dataset in the non-monsoon season. The deviation from other factors dominates over most land surfaces of the TP, except cropland, grassland, and forest. As shown in Figure 8c,d, although most points are inside the red dashed box, the points for tundra and snow/ice show significant deviation. The other-factor deviations of the snow or ice surface type almost reach <sup>−</sup>3.0 mm d−<sup>1</sup> in CPM. This is in line with the former analysis that precipitation in the monsoon season dominates the ET differences, while other effects such as surface radiation or deviations in the land-surface model play a relatively important role in the non-monsoon season. *Water* **2021**, *13*, x FOR PEER REVIEW 16 of 21 deviation. The other-factor deviations of the snow or ice surface type almost reach −3.0 mm d−1 in CPM. This is in line with the former analysis that precipitation in the monsoon season dominates the ET differences, while other effects such as surface radiation or deviations in the land-surface model play a relatively important role in the non-monsoon season.

As known, in the non-monsoon season, radiation is the dominant forcing that controls the energy balance of melting. To further explore the ET deviations over the tundra and snow/ice surfaces in the non-monsoon season, the relative contributions from surface net radiation deviation and land-surface process deviations were calculated and intercompared (Figure 9). All points are located above the 1:1 line, regardless of DDM or CPM, indicating that the lower ET of these two surfaces is mainly because of the underestimation of radiation, which may cause stronger freezing in the model. Specifically, radiation deviation of the tundra and snow/ice surfaces accounts for 90% and 62%, respectively, while the deviation from the land-surface processes is only 10% and 28%, in DDM. In CPM, the radiation deviation accounts for an even larger proportion, exceeding 90%. As known, in the non-monsoon season, radiation is the dominant forcing that controls the energy balance of melting. To further explore the ET deviations over the tundra and snow/ice surfaces in the non-monsoon season, the relative contributions from surface net radiation deviation and land-surface process deviations were calculated and intercompared (Figure 9). All points are located above the 1:1 line, regardless of DDM or CPM, indicating that the lower ET of these two surfaces is mainly because of the underestimation of radiation, which may cause stronger freezing in the model. Specifically, radiation deviation of the tundra and snow/ice surfaces accounts for 90% and 62%, respectively, while the deviation from the land-surface processes is only 10% and 28%, in DDM. In CPM, the radiation deviation accounts for an even larger proportion, exceeding 90%.

**Figure 9.** Relative contributions (normalized deviations between DDM (dynamical downscaling modeling) or CPM (convection-permitting modeling) and ERA-Interim) from surface net radiance deviation and land-surface processes deviation (unit: %) for tundra and snow or ice land-use categories in the non-monsoon season. The ordinate is the deviation from radiance; the abscissa is deviation from the land-surface processes. The black solid line indicates 1:1 correspondence of the first and third quadrants. **Figure 9.** Relative contributions (normalized deviations between DDM (dynamical downscaling modeling) or CPM (convection-permitting modeling) and ERA-Interim) from surface net radiance deviation and land-surface processes deviation (unit: %) for tundra and snow or ice land-use categories in the non-monsoon season. The ordinate is the deviation from radiance; the abscissa is deviation from the land-surface processes. The black solid line indicates 1:1 correspondence of the first and third quadrants.

#### **4. Discussion 4. Discussion**

The ET from remote sensing, off-line and on-line simulations was intercompared over the TP, and some significant deviations were found. There are, however, some drawbacks in this study. The ET from remote sensing, off-line and on-line simulations was intercompared over the TP, and some significant deviations were found. There are, however, some drawbacks in this study.

Site observations are generally considered to be the direct and most accurate measurement of ET, some atmospheric field experiments have gradually been carried out over the TP since 1970 [80–82], and several observation stations have been established [83].

#### *4.1. Uncertainties of Usage Data*

Site observations are generally considered to be the direct and most accurate measurement of ET, some atmospheric field experiments have gradually been carried out over the TP since 1970 [80–82], and several observation stations have been established [83]. However, due to the limited measurement space of one instrument, which can only measure a limited area within a few hundred to several kilometers [84], the existing observation sites on the TP are too few and insufficient to represent the ET over the whole TP.

Therefore, we selected the high-resolution remote sensing dataset EB, which performs well in land–atmosphere interactions over China, as the benchmark in this study. EB is based on the SEBS method, which requires parameterization of excessive resistance and is sensitive to the errors of air temperature (Ta) and land-surface temperature (Ts) [29]. Unlike the clear-sky conditions, the Ts used for ET calculation in EB can only be obtained by interpolation on cloudy days, and there is always error in interpolation, leading to uncertainty in the EB data. Moreover, many previous studies have noted that the uncertainty of lake parameterization and its impacts on regional climate are non-negligible [50,85,86] in current models; however, the EB data exclude the surfaces covered by lakes, meaning these areas were excluded in the present study, this is also a limitation of the EB dataset. Thus, multiple datasets from remote sensing can be considered for use in the following studies to reduce the uncertainties caused by a single dataset such as EB.

The other five datasets used in the study are all simulated from different models. GLDAS is the "off-line" simulation, which has wide applications in research. The data quality of GLDAS is mainly affected by the driving variables such as temperature, relative moisture, and wind speed [87]. Meanwhile, the lack of feedback mechanism also leads to the lack of the land–atmosphere interaction processes, which increases the uncertainties of GLDAS data as well.

The reanalysis datasets ERA-Interim and ERA5 are objective datasets obtained by conducting quality control on a series of observation data, including site observations, satellites, radar, soundings, aircraft, and ship, and then calculated through the numerical weather forecast data assimilation system. Therefore, these many sources of data bring to the reanalysis not only the advantages of a long-time scale and high spatial resolution but also the problem of poor data continuity [70,79]. The inherent uncertainties also come from the forecast model and data assimilation [88].

As for the "on-line" simulated DDM and CPM datasets, the atmospheric variables such as precipitation are also from simulation rather than being simply input as driving data. The selection of parameters in the model and the uncertainties of the parameterization scheme will influence the accuracy of ET, precipitation, and other variables simulated by DDM and CPM.

In addition, limited by the available computational resources, only a one-year simulation was conducted in this study. Uncovering whether the above results are of climatological significance would require further analysis with multi-year data.

#### *4.2. Parameterizations of the Land-Surface Process*

The analysis of the attribution of the ET deviation was only separated into precipitation and radiation deviation at present, with consideration of the water balance and energy balance. As is known, ET consists of soil evaporation, canopy transpiration, and canopyintercepted water evaporation. Further in-depth study and analysis regarding the impacts of land-surface processes on ET deviations need to be carried out.

Land-surface processes can be solved by parameterization in models; however, previous studies indicated that despite the improvements offered by current simulation capabilities, there are still numerous uncertainties in the parameterization of land-surface processes. The DDM and CPM simulations underestimated the ET over vegetated ground, and the inaccurate representation of canopy light use, interception loss, and root water uptake processes resulting in underestimation of plant transpiration might be the reason for this [89,90]. Furthermore, the difference in the months when the coldest soil temperature

appears between the reanalysis data and the WRF simulations (not shown) also indicates an insufficient consideration of snow parameterization, such as the higher snow albedo feedback mechanism, which results in lower incoming radiation and ultimately a cold deviation of the soil temperature [19,51].

Therefore, further investigations are necessary to improve the simulation effect of models and to obtain a more comprehensive and accurate understanding of the ET and water cycle process over the TP.

#### **5. Conclusions**

In this study, the ET from satellite merged data (EB), an off-line run (GLDAS), two global climate datasets (ERA5 and ERA-Interim) and two WRF dynamical downscaling simulations (DDM and CPM) were intercompared for a whole year, 2014, which is a normal year for precipitation, over the TP in terms of the annual mean, seasonal variation, and land-use types. Factors contributing to the deviations compared to EB were investigated in order to provide some insight into subsequent model improvements and land–atmosphere interaction simulations. The major conclusions are as follows:

Compared with the EB data, GLDAS generally reproduces the annual mean magnitude and spatial distribution, but seasonal variations are poorly presented. ERA-Interim and ERA5 generally overestimate ET, mainly in the monsoon season. ERA-Interim better reproduces the ET in the non-monsoon season. DDM and CPM perform well in the monsoon season but worse in the non-monsoon season.

ET is underestimated over tundra and snow or ice by all five datasets. GLDAS, ERA5, DDM, and CPM underestimate the ET over cropland, grassland, and shrubland in the non-monsoon season. GLDAS, ERA-Interim, and ERA5 overestimate the ET over grassland and shrubland in the monsoon season. ERA-Interim and ERA5 also overestimate the ET over the other four land-use types in the monsoon season.

Atmospheric forcing plays a dominant role in the simulation of ET. The considerable overestimation of precipitation dominates the ET deviation for the whole year, which is responsible for the particularly high results of ERA-Interim and ERA5 in the monsoon season. In the non-monsoon season, surface net radiation plays a secondary role after precipitation in the simulation of ET. The land-surface model exerts less impact than surface net radiation.

**Author Contributions:** Conceptualization, methodology: J.D. and Y.G.; validation: J.D., Y.G. and M.Z.; formal analysis, investigation: J.D.; data curation: Y.G.; writing—original draft preparation: J.D.; writing—review and editing: J.D., Y.G. and M.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was jointly supported by the Second Tibetan Plateau Scientific Expedition and Research Program (grant 2019QZKK010314) and the Strategic Priority Research Program of the Chinese Academy of Sciences (grant XDA2006010202).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors thank Xuelong Chen for his valuable suggestions on this research. The authors thank the National Supercomputer Center in Tianjin for computational resources. The authors appreciate the National Tibetan Plateau Data Center for free access to the EB dataset (http://data.tpdc.ac.cn, 20 June 2021), the National Aeronautics and Space Administration (NASA) for free access to the GLDAS dataset, and the European Center for Medium-Range Weather Forecasts (ECMWF) for free access to the ERA5 and ERA-Interim datasets.

**Conflicts of Interest:** The authors declare no conflict of interest.
