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

: The newly developed WRF-Hydro model is a fully coupled atmospheric and hydrological processes model suitable for studying the intertwined atmospheric hydrological processes. This study utilizes the WRF-Hydro system on the Three-River source region. The Nash-Sutcliffe efﬁciency for the runoff simulation is 0.55 compared against the observed daily discharge amount of three stations. The coupled WRF-Hydro simulations are better than WRF in terms of six ground meteorological elements and turbulent heat ﬂux, compared to the data from 14 meteorological stations located in the plateau residential area and two ﬂux stations located around the lake. Although WRF-Hydro overestimates soil moisture, higher anomaly correlation coefﬁcient scores (0.955 versus 0.941) were achieved. The time series of the basin average demonstrates that the hydrological module of WRF-hydro functions during the unfrozen period. The rainfall intensity and frequency simulated by WRF-Hydro are closer to global precipitation mission (GPM) data, attributed to higher convective available potential energy (CAPE) simulated by WRF-Hydro. The results emphasized the necessity of a fully coupled atmospheric-hydrological model when investigating land-atmosphere interactions on a complex topography and hydrology region. 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.


Introduction
The hydrological processes at the surface-atmosphere interface are a key process in the terrestrial water cycle and impact weather systems [1][2][3]. Soil moisture change, transpiration, and runoff are the three main components which affect water resources management, agriculture, ecology, and flood prediction [4][5][6][7]. Such hydrological processes are becoming increasingly detailed and precise in regional climate models [8][9][10][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][15][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].

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 × 10 5 km 2 ), the Yangtze River basin (YARB, 1.37 × 10 5 km 2 ), and the Lancang River basin (LRB, 5.3 × 10 4 km 2 ). 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 versus the uncoupled ones on the boundary layer and regional precipitation.

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 × 10 km 2 ), the Yangtze River basin (YARB, 1.37 × 10 km 2 ), and the Lancang River basin (LRB, 5. 3 × 10 km 2 ). 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].  Table 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  Table 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 [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.

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].

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.

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 Water 2021, 13, 3409 5 of 23 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.

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].

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 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), 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.

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. 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][82][83][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].

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: where V obs and V mod denote the observation and simulation values. cov(V obs , V mod ) is the covariance of V obs and V mod . σ 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 ith 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: 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.
where y and o 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: 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: 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].

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 −1 , 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.

Comparison with CMA Data
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.
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 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.

Comparison with CMA Data
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 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.   Table 2.    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% in RMSE and 8.31% in CC than WRF-S.

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 −2 for the latent heat flux and reduces it by 23.77 Wm −2 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.  Figure 5 shows the CC and RMSE distribution of soil moisture, evaporation, and air temperature at 2 m against the remote sensing or reanalysis data from 15 July to 30 November. The simulations and validation data are processed as daily averages to avoid the daily cycle dominating the evaluation measure. And bilinear interpolation in space was performed for the simulations to match with the validation data.
When comparing the soil moisture, WRF-H exhibits a lower correlation and a larger RMSE than WRF-S in most areas, especially in the west part of SRTR (Figure 5a-d). Due to the wetter soil moisture simulation, more evapotranspiration was simulated by WRF-H. Result in a slightly low CC and high RMSE simulated by WRF-H than WRF-S (Figure 5e-h). But the WRF-H exhibits advantages in simulation temperature at 2 m, with lower RMSE at the same CC. From the spatial distribution of soil moisture and evaporation, the simulations are similar in the eastern part of the Three Rivers. In contrast, the simulation of WRF-H is worse than WRF-S in the central and western regions.
After spatial validation, WRF-Hydro improved the simulation of temperature and deteriorated the humidity and evapotranspiration. The WRF-Hydro does not show the advantages in space scale as it does at the point scale. It is probably because the CMA stations were built in the more inhabitable areas of the plateau, such as in the valleys. These places are usually crossed by rivers and have a high soil moisture content. And these wetter areas cannot be resolved by GLEAM data at 0.25 • resolution but can be identified by the 4 km WRF-Hydro model. The anomaly correlation is convenient to detect similarities in the patterns of departures. And it avoids the influence of bias in the grid-based validations as they are essentially a modeled product. Figure 6a illustrates that WRF-H exhibits a higher anomaly correlation coefficient (0.955 versus 0.941). They reflect the spatial distribution of soil moisture improvements due to lateral soil water flow simulated by WRF-Hydro. Models that do not incorporate this process may lack the ability to reproduce anomaly patterns of soil moisture. Figure 6b,c show that the anomaly correlation coefficient scores achieved by WRF-H and WRF-S were close, and scores will be lower in winter. Soil water lateral flow has little effect on evapotranspiration and 2 m air temperature patterns. Figure 5. Evaluation of soil moisture, evapotranspiration, and air temperature at 2 m over the SRTR. (a) Spatial distribution of correlation coefficient (CC) between satellite retrievals or reanalysis data and soil moisture simulations from WRF−H during the study period; (b) Same as (a), but for spatial distribution of RMSE; (c,d) Same as (a,b), but simulations are from WRF−S; (e−h) Same as (a−d), but for evapotranspiration; (i−l) Same as (a−d), but for air temperature at 2 m. The dots in (a) and (c) denote the correlation coefficient and are significantly above the 95% confidence level, and the correlation coefficient in each grid point of (e,g,i,k) passed the significance test.
The anomaly correlation is convenient to detect similarities in the patterns of departures. And it avoids the influence of bias in the grid-based validations as they are essentially a modeled product. Figure 6a illustrates that WRF-H exhibits a higher anomaly correlation coefficient (0.955 versus 0.941). They reflect the spatial distribution of soil moisture improvements due to lateral soil water flow simulated by WRF-Hydro. Models that do not incorporate this process may lack the ability to reproduce anomaly patterns of soil moisture. Figure 6b,c show that the anomaly correlation coefficient scores achieved by WRF-H and WRF-S were close, and scores will be lower in winter. Soil water lateral flow has little effect on evapotranspiration and 2 m air temperature patterns.    Figure 7 compares the total precipitation, average soil moisture, latent heat and sensible heat fluxes simulated by WRF-H and WRF-S. The spatial distribution of the two sets of the models is similar (Figure 7a,b). The difference (Figure 7c) is prominent in local areas but small when averaged across the region, consistent with other study [43]. The characteristics of the soil moisture, in contrast, are opposite to the spatial distribution of precipitation. WRF-H spatially overestimates soil moisture throughout the SRTR and, more significantly, in the northern part of the SRTR (Figure 7f). The result is that WRF-S simulates soil moisture as wet in the southeast and dry in the northwest (Figure 7e), while WRF-H will hamper such a spatial distribution feature (Figure 7d). Figure 7 compares the total precipitation, average soil moisture, latent heat and sensible heat fluxes simulated by WRF-H and WRF-S. The spatial distribution of the two sets of the models is similar (Figure 7a,b). The difference (Figure 7c) is prominent in local areas but small when averaged across the region, consistent with other study [43]. The characteristics of the soil moisture, in contrast, are opposite to the spatial distribution of precipitation. WRF-H spatially overestimates soil moisture throughout the SRTR and, more significantly, in the northern part of the SRTR (Figure 7f). The result is that WRF-S simulates soil moisture as wet in the southeast and dry in the northwest (Figure 7e), while WRF-H will hamper such a spatial distribution feature (Figure 7d). The spatial distribution of latent heat has the same characteristics as soil moisture (Figure 7g-i), while the sensible heat flux is the opposite (Figure 7j-l). Changes in the distribution of sensible and latent heat fluxes affect the boundary layer development and influence the precipitation structure. Overall, the pattern of spatial changes in soil moisture, sensible and latent heat flux is quite similar. The temporal variation will be analyzed next. Figure 8 demonstrates simulated daily variable averages on the TRSR from 15 July to 30 November. Although the precipitations from both models match over the entire analysis period, the soil moisture simulated by WRF-H is more humid than WRF-S during the The spatial distribution of latent heat has the same characteristics as soil moisture (Figure 7g-i), while the sensible heat flux is the opposite (Figure 7j-l). Changes in the distribution of sensible and latent heat fluxes affect the boundary layer development and influence the precipitation structure. Overall, the pattern of spatial changes in soil moisture, sensible and latent heat flux is quite similar. The temporal variation will be analyzed next. Figure 8 demonstrates simulated daily variable averages on the TRSR from 15 July to 30 November. Although the precipitations from both models match over the entire analysis period, the soil moisture simulated by WRF-H is more humid than WRF-S during the analysis period. Latent and sensible heat flux declines with seasonal changes, and even negative values emerge in sensible heat flux. The heat fluxes in the Tibetan Plateau are more susceptible to the freeze-thaw process than the high-latitude frozen soil regions [95]. The latent heat flux simulated by WRF-H is greater than WRF-S before the soil freezes, but the sensible heat flux is more petite than WRF-S. After the soil is frozen, WRF-H coincides with the curve of the WRF-S simulation. Due to the energy balance, the surface skin temperature simulation also showed differences before freezing. negative values emerge in sensible heat flux. The heat fluxes in the Tibetan Plateau are more susceptible to the freeze-thaw process than the high-latitude frozen soil regions [95]. The latent heat flux simulated by WRF-H is greater than WRF-S before the soil freezes, but the sensible heat flux is more petite than WRF-S. After the soil is frozen, WRF-H coincides with the curve of the WRF-S simulation. Due to the energy balance, the surface skin temperature simulation also showed differences before freezing. The albedo map shows that the two models' albedo is relatively stable in summer and dramatic in autumn and winter. The steady increase of albedo in summer is caused by the decrease in soil moisture and phenology. The change in albedo due to inconsistent precipitation simulation time is small-the accumulation and sublimation of snowfall The albedo map shows that the two models' albedo is relatively stable in summer and dramatic in autumn and winter. The steady increase of albedo in summer is caused by the decrease in soil moisture and phenology. The change in albedo due to inconsistent precipitation simulation time is small-the accumulation and sublimation of snowfall cause the fluctuated albedo. In September, due to differences in surface temperature simulations, the albedo of the WRF-H is higher than the WRF-S, then the albedo of the two models coincide as the temperature difference decreases.

Time Variation
WRF-H is more consistent with the WRF-S in winter because of the weakening of overland flow and subsurface lateral flow. The decrease of soil hydraulic conductivity of frozen soil reduced the subsurface flow. Furthermore, most winter precipitation falls to the ground in snowfall and dissipates through sublimation, rarely infiltrating into the soil layer. The snowmelt will contribute to discharges mainly during the rainy and peak flow periods [70]. Figure 9 shows the fitting results of GPM precipitation, WRF-H, and WRF-S. The circle points represent the natural double pairs with different rainfall intensities, and the straight lines represent the fitted curves. The simulations of WRF-H and WRF-S are similar, but WRF-H is slightly closer to the GPM. According to Equation (7), the fitted α is 2.39, 2.74 and 2.8, and the fitted β is 26.63, 9.9 and 9.12 for GPM, WRF-H and WRF-S, respectively. The difference between simulation and GPM, WRF-H is reduced by 15% compared to WRF-S in a and β is reduced by 5%. The fitting results show that soil water lateral flow in WRF-H makes the frequency-intensity structure closer to the observation. Since convectivepermitting is used in both models, the boundary layer that triggers convection needs to be further evaluated. Figure 9 shows the fitting results of GPM precipitation, WRF-H, cle points represent the natural double pairs with different rainfall straight lines represent the fitted curves. The simulations of WRF-H a lar, but WRF-H is slightly closer to the GPM. According to Equation (7 2.74 and 2.8, and the fitted β is 26.63, 9.9 and 9.12 for GPM, WRF-H tively. The difference between simulation and GPM, WRF-H is reduce to WRF-S in a and β is reduced by 5%. The fitting results show that so in WRF-H makes the frequency-intensity structure closer to the obs vective-permitting is used in both models, the boundary layer that needs to be further evaluated.

Boundary Layer Variables
In Figure 10a, WRF-H shows a higher convective available potent most of the unfreezing period, and also that coupling hydrology proc tional energy more readily available in SRTR. The amplitude of co (CIN) did not change in the unfreezing or freezing period between t WRF-H exhibits a slightly lower lifting condensation level (LCL) in m period, indicating a low cloud base. The level of free convection (LFC) in the unfrozen season but becomes stable when it uplifts to the top o frozen season. The improvement of precipitation structure is mainly crease of CAPE caused by wetter soil.

Boundary Layer Variables
In Figure 10a, WRF-H shows a higher convective available potential energy (CAPE)in most of the unfreezing period, and also that coupling hydrology processes make convectional energy more readily available in SRTR. The amplitude of convective inhibition (CIN) did not change in the unfreezing or freezing period between the two simulations. WRF-H exhibits a slightly lower lifting condensation level (LCL) in most of the unfrozen period, indicating a low cloud base. The level of free convection (LFC) fluctuates violently in the unfrozen season but becomes stable when it uplifts to the top of model layers in the frozen season. The improvement of precipitation structure is mainly attributed to the increase of CAPE caused by wetter soil.

Discussion
We evaluated the runoff simulation skills of uncoupled WRF-Hydro and the effect of lateral flow of soil water on the atmosphere by an enhanced fully coupled WRF-Hydro model in SRTR. This is one of the earliest studies using WRF-Hydro in the SRTR, which is located on the Tibetan Plateau and is an important watershed for ecological conservation and water management in China. The results reveal the following main conclusions: 1. The reproduction of the daily discharge amount of the three stations has a Nash-Sutcliffe model efficiency generally above 0.55, demonstrating the potential of WRF-Hydro for hydrological forecasting in SRTR. 2. In the coupled experiment, WRF-Hydro made an improvement of 13.48% in RMSE and 8.31% in CC above WRF-ARW when compared against CMA data, and an improvement of 6.6% in RMSE and 1% in CC when compared against turbulent heat fluxes observation. WRF-Hydro tends to overestimate the soil moisture of the western part of the SRTR, but less so in the eastern part and areas with wetter soil moisture. This difference occurs mainly during the period when the soil is not frozen. Although the root mean square error of WRF-Hydro is larger than that of WRF-ARW compared against SMAP and GLEAM, WRF-Hydro scores higher in the anomaly

Discussion
We evaluated the runoff simulation skills of uncoupled WRF-Hydro and the effect of lateral flow of soil water on the atmosphere by an enhanced fully coupled WRF-Hydro model in SRTR. This is one of the earliest studies using WRF-Hydro in the SRTR, which is located on the Tibetan Plateau and is an important watershed for ecological conservation and water management in China. The results reveal the following main conclusions: 1.
The reproduction of the daily discharge amount of the three stations has a Nash-Sutcliffe model efficiency generally above 0.55, demonstrating the potential of WRF-Hydro for hydrological forecasting in SRTR.

2.
In the coupled experiment, WRF-Hydro made an improvement of 13.48% in RMSE and 8.31% in CC above WRF-ARW when compared against CMA data, and an improvement of 6.6% in RMSE and 1% in CC when compared against turbulent heat fluxes observation. WRF-Hydro tends to overestimate the soil moisture of the western part of the SRTR, but less so in the eastern part and areas with wetter soil moisture. This difference occurs mainly during the period when the soil is not frozen. Although the root mean square error of WRF-Hydro is larger than that of WRF-ARW compared against SMAP and GLEAM, WRF-Hydro scores higher in the anomaly correlation coefficient. These findings show the significance of lateral flow in soil moisture simulation. 3.
The coupled WRF-Hydro results in an increase in latent heat flux, a decrease in sensible heat flux, and a decrease in soil surface temperature due to the moist soil. The change in turbulent heat flux gives the WRF-Hydro simulation an enormous CAPE and easier convection, reducing precipitation intensity-frequency errors.
The NSE of runoff in this study is not high, perhaps because the agility of WRF-Hydro might be unnecessarily constrained by its complex process [96]. This is probably the reason that NSE is low in some watersheds in similar studies [97][98][99][100]. The calibrated model parameters (Zmax = 50) in this study are to offset the impact of continuously excessive precipitation on runoff production. The extra precipitation stored in the conceptual groundwater bucket during the simulation period will deteriorate runoff simulation in the next period, for the discharge of baseflow will increase the runoff in winter. If the simulation time is increased and the accumulated precipitation overestimation is large enough, then the WRF-Hydro cannot obtain an accurate runoff through calibration and the NSE will decrease as a result. It suggests that the output of the GCM is not recommended as a driver for multi-year runoff hindcasting in SRTR. Accurate precipitation-driven data, such as CMFD, is necessary for such simulation. However, using Zmax = 50 does not worsen the coupled simulation, as the parameter does not affect soil moisture and cannot further influence weather processes.
The lateral flow from WRF-Hydro leads to wetter soils in the SRTR, similar to the semiarid environment in the USA [101]. For areas where WRF-Hydro tends to overestimate the soil moisture, it is recommended to calibrate the parameters affecting soil moisture simulation, such as the reference infiltration factor (REFKDT). This study did not calibrate soil moisture because the calibrated parameters could artificially lead to better simulation results for one variable than the other.
The limitations of the study are mainly threefold. First, the overestimation of GCM precipitation in SRTR [102] resulted in its inability to be used as a force to simulate multiyear runoff. Secondly, although the soil texture dataset [67] has the highest resolution in China right now, the small number of sampling points on SRTR leads to a larger uncertainty in soil texture than in other regions in China. Thirdly, the CMA station was established in a habitable place on the plateau. These areas tend to be low-lying, with moist soils and rivers passing through them. This means that the stations are located where the lateral flow of soil water flows in, not out. More observations are needed where lateral soil water flows out.
The low but acceptable NSE of runoff simulation in the study shows the potential of WRF-hydro in the hydrological simulation of the SRTR. Meanwhile, the high anomaly correlation coefficient scores achieved by WRF-Hydro in soil moisture simulations, as well as the closer to observed precipitation intensity-frequency structure, suggest that the WRF-Hydro module is worthy of being incorporated in convective-scale simulations. Reducing the precipitation overestimation of the GCM in SRTR is urgently needed for the future application of the coupled WRF-Hydro model in SRST.

Conclusions
In the present study, we set up two sets of experiments to investigate the prospects of WRF-Hydro in the hydrological forecasting and its impacts on land-atmosphere interaction of SRTR. Results reveal the WRF-Hydro has shown potential in runoff prediction in SRTR. Coupled WRF-Hydro with soil water lateral flow increases wet soil bias in the western part of the SRTR, but improves the soil moisture anomaly pattern. The coupled model also enhances CAPE, and produces a more reliable precipitation intensity-frequency structure. Overall, our results illustrate the effect of WRF-Hydro in the coupled hydrologyatmosphere simulation system. GCM with less overestimation of precipitation in SRTR to drive coupled WRF-Hydro is desirable for future work.

Institutional Review Board Statement: Not applicable.
Data Availability Statement: The numerical model simulations upon which this study is based are too extensive to archive or transfer. Instead, we provided all the information needed to replicate the simulations; we used model version WRF-Hydro v5.0.3. The model code, compilation script, initial and boundary condition files, WRF-Hydro routing grid files, and the settings (namelist) are available from the authors.