**1. Introduction**

The hydrological processes at the surface-atmosphere interface are a key process in the terrestrial water cycle and impact weather systems [1–3]. Soil moisture change, transpiration, and runoff are the three main components which affect water resources management, agriculture, ecology, and flood prediction [4–7]. Such hydrological processes are becoming increasingly detailed and precise in regional climate models [8–11], especially as recent studies have pointed out that even with inadequate hyper-resolution meteorological and surface data, it is worthwhile to integrate lateral flow dynamics in hyper-resolution land surface models that have a resolution on the order of 100 m at continental scales [12] to better reflect water and energy heterogeneity [13]. Moreover, several numerical models incorporating lateral flow have been developed [14–16].

The Tibetan Plateau (TP) is the highest landform globally, with an average elevation exceeding 4000 m above sea level [17], and is located in the area where the East Asia monsoon and the South Asian monsoon interact. The TP acts as a strong 'dynamic pump' [18], attracting moist air from low latitudes during the warm season. The intense surface radiation and topographic lifting of TP provides a favorable condition for convection [19,20].

**Citation:** Li, G.; Meng, X.; Blyth, E.; Chen, H.; Shu, L.; Li, Z.; Zhao, L.; Ma, Y. Impact of Fully Coupled Hydrology-Atmosphere Processes on Atmosphere Conditions: Investigating the Performance of the WRF-Hydro Model in the Three River Source Region on the Tibetan Plateau, China. *Water* **2021**, *13*, 3409. https://doi.org/10.3390/w13233409

Academic Editor: Aizhong Ye

Received: 17 October 2021 Accepted: 24 November 2021 Published: 2 December 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/).

Summer precipitation is about 60% of annual precipitation [21]. In addition, the thermal forcing convective cloud system in the TP usually moves downstream into eastern China, contributing to heavy rain and severe convective storms in China [22–26]. Precipitation in winter is stored in solid form as snow and ice and released as meltwater when the temperature rises during spring and summer [27]. The huge number of glaciers, snow-packs, lakes, rivers, and a large amount of water storage in TP serves as 'the world's water tower' [20].

The source region of the Three River basins (SRTR), located in the TP, is the source region of the three largest rivers in China, i.e., the Yangtze, Lancang, and Yellow River. The region is suitable for performing the coupled atmospheric-hydrological modeling due to its vital geographic location and essential hydrological functions [28]. Based on meteorological observations and ice core records, the intensified global water cycle [29], especially precipitation and evaporation, might be enhanced in TP because TP is more sensitive to climate change due to high altitudes [30]. Rapid warming in TP [31] could lead to notable changes in the phase and intensity of surface precipitation, which impacts evapotranspiration and runoff generation [32].

Numerical models are powerful tools for studying the hydrology of the SRST, and each of the models has its advantages [33]. There are many models that focus either on the hydrology or on land-atmosphere interactions only. For example, the soil and water assessment tool (SWAT) is a conceptual, distributed parameter model [34], suitable for long-term runoff and pollution transfer simulation in large river basins. Variable infiltration capacity (VIC) is a physically-based macroscale hydrological model, developed to solve water and energy balances [35]. The community land model (CLM) [36] robustly simulates the exchange of water, energy and carbon and nitrogen between the land and the atmosphere. The community Noah land surface model with multi-parameterization options (Noah-MP) [37] is an enhanced version of Noah similar to CLM, but more oriented to applications in regional climate models. These models have been instrumental in the study of the hydrology and atmosphere of SRTR [38–42]. However, due to the design orientation, SWAT and VIC cannot be coupled directly to the atmospheric model. At the same time, the physical processes of lateral terrestrial water flow are absent in the CLM, Noah, and Noah-MP models. As a result, the impacts of the coupled hydrology and atmosphere processes are still less investigated in the SRTR.

The hydrologically enhanced version of the weather research and forecasting modeling system (WRF-Hydro) [43] is a high-resolution model that explicitly describes the surface overland flow. And it is a fully coupled atmospheric-hydrological model, since the hydrological processes are added in the WRF model (using Noah-MP as the land surface module). WRF-Hydro contains physical processes including the re-infiltration, surface overland flow and lateral flow of water within the soil layers [11]. The changes in soil moisture due to lateral flow are fed back to the atmosphere. Therefore, the precipitation simulation is different from that of WRF [44–46]. Besides, the accuracy of turbulent heat fluxes increased in Germany using WRF-Hydro [1]. The WRF-Hydro revealed good performance in simulating the diurnal cycle of land surface states and fluxes during the North American monsoon [47]. Coupled WRF-Hydro slightly outperforms the WRF stand-alone concerning sensible and ground heat flux, near-surface mixing ratio and temperature, boundary layer profiles of temperature [1]. The hydrologically enhanced process of WRF-Hydro showed an increased precipitation recycling rate in the inland area of China [48] but decreased slightly in East Africa [49]. In addition, streamflow corresponds with observations at monthly timescales on the south side of the Himalayas [50].

In this paper, the fully coupled WRF-Hydro is chosen to carry out modeling experiments of SRTR. The aim of this work is (1) to validate the applicability of WRF-Hydro in the runoff simulation in the SRTR, (2) in order to investigate the effects of lateral flow on soil moisture and near-surface meteorological variables in a coupled atmospheric-hydrological model, (3) to investigate the effects of coupled atmospheric-hydrological processes versus the uncoupled ones on the boundary layer and regional precipitation.

#### **2. Study Area, Model Description and the Numerical Experiment Design 2. Study Area, Model Description and the Numerical Experiment Design**  *2.1. Study Area*

versus the uncoupled ones on the boundary layer and regional precipitation.

*Water* **2021**, *13*, x FOR PEER REVIEW 3 of 25

## *2.1. Study Area*

The SRTR is located northeast of the TP, with an altitude average of 4510 m and ranging from 2600 m to 6500 m a.s.l. (Figure 1). It is an adjacent area consists of the source region of the Yellow River basin (YRB, 1.22 <sup>×</sup> <sup>10</sup><sup>5</sup> km<sup>2</sup> ), the Yangtze River basin (YARB, 1.37 <sup>×</sup> <sup>10</sup><sup>5</sup> km<sup>2</sup> ), and the Lancang River basin (LRB, 5.3 <sup>×</sup> <sup>10</sup><sup>4</sup> km<sup>2</sup> ). The SRST is mainly covered by grassland [51], and the soil type is loam. The YRB contributes about 35% of the river flow [52]. Glacier cover fractions of YARB and YRB are 0.95% and 0.11% [33], respectively. The proportion of melt contributing to runoff ranges from 3.9 to 6% in four subbasins of the YARB [42]. During 1956–2012, the annual runoff in LRB and YARB increased, while YRB slightly decreased [53]. Under the impact of climate change, frozen ground degradation increased the groundwater discharge rate in winter [54]. Direct snowmelt runoff coefficients are mainly controlled by the air temperature freezing index [55]. The SRTR is located northeast of the TP, with an altitude average of 4510 m and ranging from 2600 m to 6500 m a.s.l. (Figure 1). It is an adjacent area consists of the source region of the Yellow River basin (YRB, 1.22 × 10ହ km2), the Yangtze River basin (YARB, 1.37 × 10ହ km2), and the Lancang River basin (LRB, 5. 3 × 10ସ km2). The SRST is mainly covered by grassland [51], and the soil type is loam. The YRB contributes about 35% of the river flow [52]. Glacier cover fractions of YARB and YRB are 0.95% and 0.11% [33], respectively. The proportion of melt contributing to runoff ranges from 3.9 to 6% in four subbasins of the YARB [42]. During 1956–2012, the annual runoff in LRB and YARB increased, while YRB slightly decreased [53]. Under the impact of climate change, frozen ground degradation increased the groundwater discharge rate in winter[54]. Direct snowmelt runoff coefficients are mainly controlled by the air temperature freezing index [55].

In this paper, the fully coupled WRF-Hydro is chosen to carry out modeling experiments of SRTR. The aim of this work is (1) to validate the applicability of WRF-Hydro in the runoff simulation in the SRTR, (2) in order to investigate the effects of lateral flow on soil moisture and near-surface meteorological variables in a coupled atmospheric-hydrological model, (3) to investigate the effects of coupled atmospheric-hydrological processes

**Figure 1.** Geographic locations of the observational stations and source region of Three River basins (SRTR) include the Yellow River basin (YRB), the Yangtze River basin (YARB) and the Lancang River basin (LRB). The coordinates of streamflow stations, flux towers and meteorological stations are shown in Table 1. **Figure 1.** Geographic locations of the observational stations and source region of Three River basins (SRTR) include the Yellow River basin (YRB), the Yangtze River basin (YARB) and the Lancang River basin (LRB). The coordinates of streamflow stations, flux towers and meteorological stations are shown in Table 1.

#### *2.2. Model Description 2.2. Model Description*

#### 2.2.1. WRF-ARW and NOAH-MP 2.2.1. WRF-ARW and NOAH-MP

The WRF-ARW model [56] is a time-split non-hydrostatic atmospheric model. It is widely used in the study of land-atmosphere interactions [6], dynamic downscaling The WRF-ARW model [56] is a time-split non-hydrostatic atmospheric model. It is widely used in the study of land-atmosphere interactions [6], dynamic downscaling [57,58], weather and climate research [59,60]. The ability of WRF to resolve strongly nonlinear small-scale phenomena such as stratiform precipitation and convective precipitation makes it suitable for the coupling study of atmospheric and hydrological processes.

Noah-MP [37] is an enhanced version of the Noah land surface model, and introduces a framework for multiple schemes. Due to the enhancement of biophysical and soil freezethaw processes, the Noah-MP has been widely used in the TRSR [41,61,62]. Noah-MP can be selected as the land surface process module of the WRF model.

#### 2.2.2. WRF-Hydro

WRF-Hydro is a modeling framework to facilitate the coupling of WRF and hydrological models. The interactions between hydrological and land surface processes are calculated as follows [43]. Land surface states and fluxes computed by Noah-MP are disaggregated to the high-resolution terrain routing grid. Then the physical hydrologic process is calculated on the routing grid. After that, land surface states and fluxes aggregated from the routing grid are updated to the Noah-MP model grid. Therefore, it can be used as a land surface model offline or fully coupled to the WRF model. The uncoupled WRF-Hydro model is good at spin-up, model calibration, and data assimilation. In contrast, the coupled WRF-Hydro model is used for hydrology-atmosphere coupling research.

There are five major parts in the hydrologic processes of WRF-Hydro [43], namely subsurface flow, overland flow, channel, lake, and a conceptual base flow module. The subsurface flow process calculates subsurface lateral flow [63] and exfiltration from a supersaturated soil column. During overland flow, the infiltration excess and exfiltration calculated in the previous step flow into the river by solving the diffusive wave formula [64]. The channel process simulates the flow in the river network by either fine grid routing using a diffusive wave equation or on a vectorized network of channel reaches by solving the Muskingum or Muskingum-Cunge equation [65]. Besides, WRF-Hydro provides a simple mass balance, a level-pool lake/reservoir routing module, and a conceptual exponential bucket base flow module [43].

#### *2.3. Numerical Experiment Design*

To investigate the performance of WRF-Hydro, uncoupled and coupled simulations were performed over SRST. The first is a WRF-Hydro uncoupled experiment (WRF-H-UP, hereafter), used to find out the calibrated hydrological parameters and test the simulation performance of runoff generation. Secondly, a group of coupled WRF-Hydro (WRF-H, hereafter) and WRF-ARW experiments (WRF-S, hereafter) were conducted to test the effects of hydrological processes on the atmosphere simulation.

#### 2.3.1. Model Calibration

A WRF-Hydro simulation is set up in offline mode (WRF-H-UP) for parameter calibration and runoff simulation. Since it is still a challenge to reproduce daily runoff using fully coupled models [48], only runoff from the uncoupled simulation is analyzed in the study. The domain of WRF-H-UP corresponds to the area of the inner domain (d02) used by the subsequent coupling run, shown in Figure 2. The simulations of WRF-ARW drive the WRF-H-UP to verify the forecasting ability of WRF-Hydro on the runoff discharge. The runoff simulations are produced after the parameters have been calibrated. Because WRF-Hydro includes many parameterized nonlinear physical processes, the default parameters given by WRF-Hydro are only valid over a small region; calibration of related model parameters is often required to use in the new domain [66].

The parameters are calibrated manually in this study. Most of the parameters are insensitive during the calibration, except those that control base flow (maximum depth, Zmax). By increasing the Zmax from 10 to 50, WRF-Hydro can calibrate errors caused by overestimated precipitation by storing water in conceptual groundwater buckets. Other parameters involved in the calibration process include the control of water movement in the soil, such as reference infiltration factor (REFKDT), soil evaporation exponent (FX-EXP\_DATA), deep drainage coefficient (SLOPE), saturated soil hydraulic conductivity (SATDK), saturated hydraulic conductivity coefficient of lateral flow (LKSATFAC), porosity (MAXSMC), and filed capacity (REFSMC). Parameters are relevant to overland flow, that is, the roughness coefficient (SFC\_ROUGH), and Manning's roughness coefficient (MannN). Parameters required for lake routing include coefficient (WeirC), weir length (WeriL), orifice coefficient (OrificeC), orifice area (OrificeA). Groundwater parameters include the coefficient of baseflow (Coeff) and exponent (Expon) of the bucket model. However, adjusting

these parameters had little effect in the runoff simulations, so all default values were used in coupled simulations except for Zmax, which was calibrated from 10 to 50.

#### 2.3.2. Coupled Simulations

The WRF-ARW (WRF-S, hereafter) model is configured standalone on a two-way nested domain of 20.4 km cell resolution, covering east Asia and SRTR, respectively (Figure 2). The simulation has 33 vertical layers up to 50 hPa and is driven by the European Centre for Medium-Range Weather Forecasts Reanalysis (ERA) Interim, which has a resolution of 0.75 degrees. The land surface static physiographic input is generated by WRF Preprocessing Tools, in which the MODIS IGBP 21-category data is used to interpolate land use categories. Soil data is from a comprehensive 30 arc-second resolution grided soil characteristics data of China [67]. Grell-Devenyi [68] is selected as Cumulus parameterization options in domain1, and the convection-permitting module is used for domain2. The free-drainage approach [8,69] is selected as the lower boundary condition. A previous study shows that parameter calibration cannot resolve the deficiency of groundwater table-based parameterizations in simulating runoff [62]. Other selected parameterization schemes in this study are shown in Table S1 (in Supplementary Materials). The hindcast simulation period is from 1 June 2018 to 30 November 2018, with the first 45 days for warm-up and the rest for comparison. The reason for setting such a long warm-up time is that 1.5 spin-up days are needed for precipitation, but a much more extended period is necessary for discharge due to the influence of soil moisture [70]. *Water* **2021**, *13*, x FOR PEER REVIEW 6 of 25

**Figure 2.** The terrain elevation map of the WRF nested domains. WRF-Hydro has the same settings as WRF, except the inner domain (d02) is coupled with a routing subgrid in WRF-Hydro. **Figure 2.** The terrain elevation map of the WRF nested domains. WRF-Hydro has the same settings as WRF, except the inner domain (d02) is coupled with a routing subgrid in WRF-Hydro.

**3. Validation Data and Analysis Metrics**  *3.1. Validation Data*  For model validation, in situ data were collected from three hydrological discharge To investigate the influence of the WRF-Hydro extension, the coupled WRF-Hydro model (WRF-H, hereafter) is set using the same parameters, driven by the same forcing data in the same hindcast period with the WRF-Standalone model (WRF-S, hereafter),

stations, two eddy covariance stations, and fourteen meteorological stations (Figure 1).

a.m. each day. Data periods are from July to November 2018. The Northwest Institute of Eco-Environment and Resources, Chinese Academy of Science in SRTR, set up two eddycovariance systems to record the sensible and latent heat fluxes every 30 minutes. Meteorological station data were supplied by the China Meteorological Administration (CMA), and the observations were made at 0:00, 6:00, 12:00, and 18:00 UTC, and data were pro-

cessed as daily averages.

except the inner domain of WRF-H is coupled with a sub-grid at a 400-m resolution to compute overland and river flow. The refined routing grids in the WRF-H extension are created by WRF-Hydro GIS Pre-Processing Tools, Version5.1, with the aggregation set to 10. The input elevation data is from the Shuttle Radar Topography Mission (SRTM), which has a spatial resolution of 90 m [71]. The simulations start from 1 June using the same initial and boundary conditions, and the data before 15 July are used for spin up.

#### **3. Validation Data and Analysis Metrics**

#### *3.1. Validation Data*

For model validation, in situ data were collected from three hydrological discharge stations, two eddy covariance stations, and fourteen meteorological stations (Figure 1). The primary information of the observation data is shown in Table 1, and the geographical coordinates of the stations are shown in Table 2. The runoff data was observed at 12:00 a.m. each day. Data periods are from July to November 2018. The Northwest Institute of Eco-Environment and Resources, Chinese Academy of Science in SRTR, set up two eddy-covariance systems to record the sensible and latent heat fluxes every 30 min. Meteorological station data were supplied by the China Meteorological Administration (CMA), and the observations were made at 0:00, 6:00, 12:00, and 18:00 UTC, and data were processed as daily averages.

**Table 1.** Information of validation data, including China Meteorological Administration observation (CMA), Hydrologic data, Eddy covariance data, Global the Precipitation Mission (GPM) level\_3 IMERG final product, the Soil Moisture Active Passive (SMAP) mission level-4 product, Global Land Evaporation Amsterdam Model (GLEAM), and China Meteorological Forcing Dataset (CMFD).


In addition, the grided dataset in Table 1 is used for validation. The GPM level\_3 IMERG final product [72] is the rainfall estimates combining data from all passive-microwave instruments in the GPM constellation [73]. The Soil Moisture Active Passive (SMAP) Mission is a satellite-based soil moisture product [74]. SMAP mission level-4 product is a modeled product. Based on brightness temperature observations (SMAP L1C\_TB), it assimilates SMAP L1C\_TB into the NASA catchment land surface model using a spatially distributed ensemble Kalman filter [75,76], producing estimations of surface (0–5 cm) and root zone (0–100 cm) soil moisture estimation. GLEAM (Global Land Evaporation Amsterdam Model) is a set of algorithms [77] providing evapotranspiration and estimating its components separately. The dataset of version 3.3b [78] is mainly based on satellite data. The China Meteorological Forcing Dataset (CMFD) [79] is a gridded near-surface meteorological dataset, which is made by fusion ground-based observations with several gridded datasets from remote sensing and reanalysis. The CMFD dataset was only used for verification in this study.


**Table 2.** The coordinates of hydrographic stations (red pentagons in Figure 1), turbulent heat flux stations (pink pentagrams in Figure 1), and China Meteorological Administration ground observation (CMA) stations (black triangles in Figure 1).

GPM is reliable for rainfall and cold season solid precipitation estimates [80], making it suitable for evaluating precipitation in SRTR [81–84]. Previous assessments evaluating the SMAP against two soil moisture observation networks on TP showed that the SMAP retrievals could well capture the amplitude and temporal variation of the soil moisture [85]. When atmospheric water balances, evapotranspiration was taken as the evapotranspiration (ET) baseline, GLEAM matches low ET estimates, and errors are amplified when ET is higher [86]. Previous studies showed that the surface air temperature downscaled from CMFD is more accurate than ERA-interim data over the TP [87]. This dataset is used for spatial validation after compensating for the deficiencies of the sparse observation of CMA weather stations in western China and systematic bias of reanalysis/remote sensing datasets [48].

#### *3.2. Skill Metrics*

The Pearson correlation coefficient (CC, range −1 to +1) and root mean square error (RMSE, range 0 to ∞) evaluate the model simulation effectiveness against the in situ observations or remote sensing data. In addition to RMSE, the Nash-Sutcliff Efficiency (NSE, range −∞ to 1) is used in runoff assessment. The formula of the criterion at a given site or grid point is:

$$\text{CC} = \frac{\text{cov}(V\_{\text{obs}}, V\_{\text{mod}})}{\sigma\_{\text{obs}}\sigma\_{\text{mod}}} \tag{1}$$

$$RMSE = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left(V\_{mod}^{i} - V\_{obs}^{i}\right)^{2}} \tag{2}$$

$$NSE = 1 - \frac{\sum\_{i=1}^{n} \left(V\_{mod}^{i} - V\_{obs}^{i}\right)^{2}}{\sum\_{i=1}^{n} \left(V\_{mod}^{i} - \mu\_{obs}\right)^{2}} \tag{3}$$

where *Vobs* and *Vmod* denote the observation and simulation values. *cov*(*Vobs*, *Vmod*) is the covariance of *Vobs* and *Vmod*. *σobs* and *σmod* are the standard deviation of the observation and simulation. *µobs* is the average value of the observation. *V i obs* and *V i mod* are the observation and simulation of the *i*th time step. In addition, Taylor diagrams providing summary statistics of CC, RMSE and standard deviation (SD) were used for model performance analysis [88].

The *t*-statistic tests the statistical significance of the linear correlation coefficient:

$$t = \frac{r\sqrt{n-2}}{\sqrt{1-r^2}}\tag{4}$$

where *r* is the linear correlation coefficient, *n* is the sample size. The significance level in this study is set at 95%.

Uncentered anomaly correlation coefficient (a number range 0 to 1) [89,90] is used in the verification of spatial fields.

$$\text{ACC} = \frac{\sum\_{m=1}^{M} y\_m' o\_m'}{\sqrt{\sum\_{m=1}^{M} \left(y\_m'\right)^2 \sum\_{m=1}^{M} \left(o\_m'\right)^2}} \tag{5}$$

where *y* 0 and *o* 0 are the departure of simulation and verifying data from a reference state at each grid point (*m*).

The final metric is a quantitative method to study the precipitation frequency-intensity structure proposed by [91]. In this method, two quantitative parameters are obtained by fitting a double exponential function with the statistical results of the frequency of different rainfall intensities. The equation is:

$$\text{Fr}(I) + 1 = \exp\left[\exp\left(\alpha - \frac{1}{\beta}I\right)\right] \tag{6}$$

where Fr(*I*) is the frequency when the precipitation intensity is in *I* (mm) categories. α and *β* are the parameters we want. For the left and right sides of Equation (6), take two natural logarithms at the same time to get:

$$\ln(\ln(\text{Fr}(I) + 1)) = \alpha - \frac{1}{\beta}I \tag{7}$$

That is to say, the double logarithm of the occurrence frequency of a certain precipitation intensity is a linear function of precipitation intensity. Previous studies have shown that the precipitation distribution of meteorological station observations, satellite remote sensing, and numerical simulation obey this law well [92,93]. In this function, α and β are related to weak and intense precipitation occurrence frequency, respectively. Numerical simulations are associated with an increased frequency of light precipitation but decreased frequency in heavy precipitation, resulting in a larger α and smaller β in simulation than observation. [92].

#### **4. Results**

#### *4.1. Streamflow Simulation*

The performance of the runoff simulation was evaluated using the parameters obtained from the calibration procedure in Section 3.2, i.e., Zmax = 50. The NSE of the source region of the Yangtze, Lancang and Yellow rivers are 0.12, 0.19 and 0.36, respectively, based on daily observations from hydrological stations (Figure 3). The effect of setting Zmax to 50 is to take the overestimation of precipitation in the WRF-S and store it in the conceptual groundwater bucket to keep the simulated streamflow close to the observations. Such results were consistent with the study in the eastern Alps [45], where after parameter-calibrated that the bias in precipitation does not deteriorate the reproduction of runoff discharges.

As the three river sources are interconnected, considering them as a whole reduces errors in runoff simulation caused by the errors in the precipitation fall zone. This precipitation fall zone error is generated by the GCM and introduced into the flow production

simulation system by the precipitation as driving information. The simulated and observed total flows of the three rivers are obtained by adding up the flows of Zhimenda, Changdu and Tangnag. Then the Nash-Sutcliff Efficiency (NSE) and root mean square error (RMSE) between simulation and observation of runoff in SRTR is 0.55 and 324.2 m s−<sup>1</sup> , respectively (Figure 3d). This considerable improvement in NSE suggests that the overestimation of precipitation by the GCM and the error in the precipitation fall zone are sources of error in flow production.

Nevertheless, The WRF-Hydro produces high temporal resolution estimates of yield flow with acceptable results. And it illustrates that the bias from meteorology is essentially addressed through calibration. The accuracy of model simulations on the daily scale will be further improved as the driving information error decreases due to more observations. *Water* **2021**, *13*, x FOR PEER REVIEW 10 of 25

**Figure 3.** Time series of observed and simulated (uncoupled WRF-Hydro) discharge of (**a**) Zhimenda station; (**b**) Changdu station; (**c**) Tangnag station; (**d**) amount of three stations. **Figure 3.** Time series of observed and simulated (uncoupled WRF-Hydro) discharge of (**a**) Zhimenda station; (**b**) Changdu station; (**c**) Tangnag station; (**d**) amount of three stations.

#### *4.2. Comparison with In-Situ Observations 4.2. Comparison with In-Situ Observations*

#### 4.2.1. Comparison with CMA Data 4.2.1. Comparison with CMA Data

than between models.

As the first step toward evaluating the coupled WRF-Hydro modeling system's performance, we present comparisons between observed and simulated variables at fourteen individual CMA stations in Figure1. The compared variables are listed in Table 1. Simulations in the WRF-S/WRF-H model grid nearest to the CMA stations are converted from hourly data to daily for comparison with observations. As the first step toward evaluating the coupled WRF-Hydro modeling system's performance, we present comparisons between observed and simulated variables at fourteen individual CMA stations in Figure 1. The compared variables are listed in Table 1. Simulations in the WRF-S/WRF-H model grid nearest to the CMA stations are converted from hourly data to daily for comparison with observations.

Taylor diagrams for different variables are shown in Figure 4. As thermodynamic variables, both the air temperature and surface pressure are simulated with a slight standard deviation. These two variables and surface skin temperature strongly correlate with observations due to daily and seasonal variation. The simulated wind speed has shown the worst correlation among all variables. The positions of the stations on the Taylor map are essentially the same on the humidity map as they are on the precipitation map. The precipitation RMSE between WRF-H/WRF-S and observations is mainly derived from the Taylor diagrams for different variables are shown in Figure 4. As thermodynamic variables, both the air temperature and surface pressure are simulated with a slight standard deviation. These two variables and surface skin temperature strongly correlate with observations due to daily and seasonal variation. The simulated wind speed has shown the worst correlation among all variables. The positions of the stations on the Taylor map are essentially the same on the humidity map as they are on the precipitation map. The precipitation RMSE between WRF-H/WRF-S and observations is mainly derived

uncertainties in the model, forcing the heavily overestimated precipitation in ERA-Interim over the TP [94]. In general, the simulation performance varies more between observations

from the uncertainties in the model, forcing the heavily overestimated precipitation in ERA-Interim over the TP [94]. In general, the simulation performance varies more between observations than between models. *Water* **2021**, *13*, x FOR PEER REVIEW 11 of 25

**Figure 4.** Taylor diagram for simulations compared with 14 CMA observational daily data for; (**a**) precipitation; (**b**) air temperature at 2 m (T2); (**c**) surface skin temperature (SST); (**d**) surface pressure; (**e**) relative humidity at 2 m (RH); (**f**) wind **Figure 4.** Taylor diagram for simulations compared with 14 CMA observational daily data for; (**a**) precipitation; (**b**) air temperature at 2 m (T2); (**c**) surface skin temperature (SST); (**d**) surface pressure; (**e**) relative humidity at 2 m (RH); (**f**) wind speed at 10 m (V10). The number near the dot in the figure is the ID of the stations listed in Table 2.

speed at 10 m (V10). The number near the dot in the figure is the ID of the stations listed in Table 2. Table 3 shows the average value of each variable for 14 stations. By comparison, all six variables from WRF-H have a smaller RMSE than WRF-S, and five variables have a better correlation with observations in WRF-H. WRF-H made an improvement of 13.48% Table 3 shows the average value of each variable for 14 stations. By comparison, all six variables from WRF-H have a smaller RMSE than WRF-S, and five variables have a better correlation with observations in WRF-H. WRF-H made an improvement of 13.48% in RMSE and 8.31% in CC than WRF-S.

> Relative humidity (%) 15.7 15.7 0.6 0.59 14 14 Wind speed (msିଵ) 1.0 1.1 0.29 0.28 12 12

in RMSE and 8.31% in CC than WRF-S. **Table 3.** Averaged Root-Mean-square Error (RMSE) and Correlation Coefficient (CC) between the measurements collected from the 14 China Meteorological Administration (CMA) stations and the simulated near-surface and land surface meteorological variables for the period from 15 July 2009 **Table 3.** Averaged Root-Mean-square Error (RMSE) and Correlation Coefficient (CC) between the measurements collected from the 14 China Meteorological Administration (CMA) stations and the simulated near-surface and land surface meteorological variables for the period from 15 July 2009 to 30 November 2018. The last column is the number of sites that passed the significance test across the 14 observatories.


## 4.2.2. Comparison with Turbulent Heat Flux Observation

We also evaluated the WRF-H/WRF-S simulations against turbulent flux measurements shown in Figure 1. Table 4 shows the mean RMS error and CC in all comparison periods. All correlation coefficients in the table pass the significance test with a confidence level of 95%. The result at Elinghu shows that WRF-H simulates a larger CC than the WRF-S in both latent and sensible heat flux. For RMSE, the WRF-H and WRF-S are superior to each other in the simulation of sensible heat flux and latent heat flux. The simulation of WRF-H increases RMSE by 15.91 Wm−<sup>2</sup> for the latent heat flux and reduces it by 23.77 Wm−<sup>2</sup> for sensible heat flux. In contrast to the Elinghu station, the WRF-H simulates a smaller RMSE than WRF-S in both latent and sensible heat flux but a larger CC for latent heat simulations and smaller CC for sensible heat simulations.

**Table 4.** Root-Mean-Square Error (RMSE) and correlation coefficient (CC) computed between WRF-Hydro/WRF and observation from 15 July to 30 November 2018, where the correlation coefficient passed the significance test.

