Simulation of Vegetation Cover Based on the Theory of Ecohydrological Optimality in the Yongding River Watershed, China

Abstract

During ecological restoration, it is necessary to comprehensively consider the state of vegetation in climate–soil–vegetation systems. The theory of ecohydrological optimality assumes that this state tends to reach long-term dynamic equilibrium between the available water supply of the system and the water demand of vegetation, which is driven by the maximization of productivity. This study aimed to understand the factors that affect the spatial distribution of vegetation and simulate the ideal vegetation coverage (M₀) that a specific climate and soil can maintain under an equilibrium state. The ecohydrological optimality model was applied based on meteorological, soil, and vegetation data during the 2000–2018 growing seasons, and the sensitivity of the simulated results to input data under distinct vegetation and soil conditions was also considered in the Yongding River watershed, China. The results revealed that the average observed vegetation coverage (M) was affected by precipitation characteristic factors, followed by wind speed and relative humidity. The M, as a whole, exhibited horizontal zonal changes from a spatial perspective, with an average value of 0.502, whereas the average M₀ was 0.475. The ecohydrological optimality theory ignores the drought resistance measures evolved by vegetation in high vegetation coverage areas and is applicable to simulate the long-term average vegetation coverage that minimizes water stress and maximizes productivity. The differences between M and M₀ increased from the northwest to the southeast of this area, with a maximum value exceeding 0.3. Meteorological factors were the most sensitive factors of this model, and the M₀ of the steppe was most sensitive to the stem fraction, mean storm depth, and air temperature. Whether soil factors are sensitive depends on soil texture. Overall, the study of the carrying capacity of vegetation in the natural environment contributes to providing new insights into vegetation restoration and the conservation of water resources.

1. Introduction

Vegetation, as the primary component of terrestrial ecosystems, interacts with hydrological processes [1,2]. The spatiotemporal responses of vegetation to water circulation mainly depend on the climate, soil, and physiological characteristics of the vegetation [3,4,5]. Precipitation is linked to vegetation through soil moisture [6]. The uncertainty of both the frequency and amount of precipitation promotes vegetation to optimize productivity [7,8]. Moreover, the soil texture is responsible for the differences in the hydraulic properties of the root zone soil and further impacts the spatial distribution of soil moisture [9]. The water available for use by vegetation in the soil increases from coarse sand, fine sand, sandy soil, loam to clay [10]. Clayey soils generally support shallow-rooted vegetation, and sandy soils favor deep-rooted vegetation [11]. Likewise, the physiological responses of plant species to moisture changes are related to their own structure and strategies for adaptation to drought, such as the crown shape, xylem structures, size and depth of roots, and leaf traits [12,13,14,15]. Diverse vegetation types show differences in water utilization efficiency and stomatal regulation, which are reflected by canopy conductance [16]. For instance, when vegetation is under water stress, reductions in the vegetation height, branch leaf area-to-sapwood area ratio, and leaf size could improve hydraulic efficiency and limit xylem tension [17,18]. The prevalent view is that natural selection and adaptation capability lead to the evolution of vegetation to the optimal state in climate–soil–vegetation systems for long periods of time [19]. Different studies have emphasized this optimal state differently, such as the stomatal optimization model from the perspective of plant physiology or the feasibility of vegetation patterns from the perspective of semiarid watersheds [20,21]. Fewer studies have considered the underlying mechanism of the optimal state from the combined aspects of vegetation physiology and hydrological processes.

This study used the theory of ecohydrological optimality proposed by Eagleson [22] (p. 341), which states that vegetation tends to prefer an average state of zero plant stress under given climate and soil conditions driven by productivity maximization [23]. Based on the precondition that the precipitation process obeys the Poisson distribution, a mathematical and physical expression of the water balance equation was established [24,25,26,27,28,29,30,31,32]. Then, the response of stomata to the available water of vegetation and light is idealized and the potential canopy conductance equation and water balance equation are coupled to formulate this average state of vegetation [23]. This model was improved, and its parameters as well as applicable conditions were subsequently verified by many scholars [33,34,35,36,37,38]. Furthermore, the advancement of satellite measurements as well as the feasibility and universality of ground surveys in China also provide plenty of data for model research [39]. For example, Mo [40] and Cong [41] compared the actual vegetation cover generated by the remotely sensed normalized difference vegetation index (NDVI) with the simulated cover from the model output and verified the applicability of the model in the Horqin Sands in China and the Northeast China Transect. Overall, the ecohydrological optimality model can obtain the ideal vegetation coverage ($M_0$) to maximize productivity and minimize water stress in water-limited or energy-limited ecosystems.

The Yongding River watershed is an important water source conservation area and ecological corridor in the Beijing–Tianjin–Hebei region of China. The upper reaches of the basin are mountainous areas with steep slopes and an uneven distribution of vegetation, and these areas are prone to soil erosion. The lower reaches are plain areas, where the ecosystem is degraded and the river is cut off. China has formulated a series of comprehensive control and ecological restoration measures. The aim of these measures is to restore the local ecosystem to a sustainable natural state. Vegetation cover can represent the growth status of vegetation to a certain extent, and it can be used as a reference index for ecological restoration. Zhang [42] discussed the problem of excessive vegetation restoration in the “Grain-to-Green Program”, which was implemented over the Loess Plateau of China, and considered that the $M_0$ value of the model can be used as a limit for afforestation. Hence, research on the ecohydrological optimality model has practical significance for artificially improving the quality and stability of ecosystems.

Although the statistical–dynamical water balance model of Eagleson is well recognized, the theory of ecohydrological optimality is still controversial and needs to be verified in dissimilar regions and vegetation types [43]. In addition, vegetation types were rarely classified in detail and were not comprehensively considered together with soil texture and precipitation characteristics in previous related studies. In this context, the primary objectives of this study were to (1) analyze the spatial distribution characteristics and impact factors of the observed vegetation coverage ($M$) during the growing season from 2000 to 2018; (2) evaluate the applicability of ecohydrological optimality theory in the Yongding River watershed; (3) calculate the simulated $M_0$ value based on local meteorological, soil, and vegetation factors; and (4) investigate the sensitivity of the model to input parameters under varied vegetation and soil conditions. The results of this study provide detailed theoretical guidance and a scientific assessment of the application of ecological restoration and water resource management.

2. Materials and Methods

2.1. Study Area

The Yongding River watershed (Figure 1a), crossing five provinces (Inner Mongolia, Shanxi, Hebei, Beijing, and Tianjin), is located southwest of Beijing (112°–117°45′ E and 38°27′–41°20′ N). The total area of this watershed is 47,016 km², with high elevations in the northwest and low elevations in the southeast, including mountain areas (45,063 km²) and plain areas (1953 km²), with Sanjiadian as the boundary. The average groundwater depth in the plain area is 4–50 m, most of which is more than 20 m (data were developed by the website of Ministry of Water Resources of the People’s Republic of China http://www.mwr.gov.cn/ accessed on 29 January 2019). The area has a temperate continental monsoon climate, with mean annual precipitation from 360 mm to 650 mm, a mean annual temperature of 6.9 °C, and a mean annual runoff of 12.13 mm (averaged from 1963 to 2018). The vegetation growing season in this region is from April to September since the precipitation is concentrated and the temperature is optimum during this period. Major vegetation types include needleleaf forest, broadleaf forest, scrub, steppe, and cultural vegetation, where cultural vegetation refers to the phytocommunity formed by artificial cultivation, including irrigated or nonirrigated crops and deciduous orchards (Figure 1b). The main soil texture types are clay loam and sandy loam (Figure 1c).

2.2. Data Preparation

Digital elevation model (DEM) data were used to divide the watershed boundaries, with a resolution of 30 m, provided by the Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn accessed on 1 December 2009). Then the spatial distributions of vegetation types and soil texture types were obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn accessed on 30 May 2001) with a resolution of 1 km. We obtained NDVI from the MOD13A3 (monthly 1-km) VI product from 2000 to 2018 to analyze the spatial distribution of $M$. The method of converting NDVI data into $M$ was improved by Li [44]:

$$M=\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}$$

(1)

We assumed that the maximum coverage of forest and the minimum coverage of barren land in the study area were 100% and 0%, respectively, corresponding to the maximum and minimum values of NDVI. Due to inevitable image noise, ${\mathrm{NDVI}}_{\max }$ and ${\mathrm{NDVI}}_{\min }$ are the maximum and minimum NDVI values, respectively, in the 95% confidence interval, spanning 19 years. The effective value range of the MOD13A3 product is from −0.2 to 1, while water bodies have NDVI values less than 0. ${\mathrm{NDVI}}_{\max }$ and ${\mathrm{NDVI}}_{\min }$ are 0.71 and 0.16, respectively, in this study.

The input data and parameters (which are listed and explained in

Table 1. The symbols, description, and values of meteorological data and vegetation data of five vegetation types.

Symbols		Description	Units	Values
				Needleleaf Forest	Broadleaf Forest	Scrub	Steppe	Cultural Vegetation
Meteorological data	$m_h$	mean storm depth during growing season	mm	13.90	13.41	14.25	12.93	13.63
	$m_{\upsilon }$	number of independent storms per growing season	time−1	28.41	28.87	28.49	29.30	28.73
	$m_{t_r}$	mean storm duration during growing season	days	0.54	0.53	0.54	0.53	0.53
	$m_{t_b}$	mean time between storms during growing season	days	5.91	5.81	5.89	5.72	5.85
	$T_{\tau }$	mean daily air temperature during growing season	°C	19.64	18.17	19.30	16.37	18.29
	$\mu$	mean daily wind speed at 2 m height during growing season	m/s	2.14	2.20	2.13	2.36	2.29
	$RH$	mean daily relative humidity during growing season	%	56.38	55.21	56.75	55.41	55.78
	$N$	mean sunshine duration per day during growing season	h	7.79	7.94	7.73	8.03	7.87
	$P_d$	precipitation during nongrowing season	mm	57.09	56.81	58.92	55.95	56.44
	$E_{ps\tau }$	mean potential rate evaporation during growing season	mm/day	4.14	4.14	4.12	4.12	4.14
	$E_{psd}$	mean potential rate evaporation during nongrowing season	mm/day	1.35	1.30	1.34	1.27	1.30
	$\Delta$	slope of saturation vapor pressure curve	kPa/°C	0.15	0.14	0.14	0.13	0.14
	${\gamma }_0$	surface psychrometric constant	kPa/°C	0.06	0.06	0.06	0.06	0.06
	$m_{\tau }$	length of the growing season	days	183	183	183	183	183
	$m_d$	length of the nongrowing season	days	182	182	182	182	182
Vegetation data	$M_{\tau }$	mean vegetation coverage during growing season	dimensionless	0.78	0.59	0.69	0.50	0.46
	$M_d$	mean vegetation coverage during nongrowing season	dimensionless	0.52	0.39	0.51	0.35	0.31
	$L_t$	mean foliage area index during growing season	dimensionless	1.82	1.29	1.47	0.92	0.86
	$\beta$	cosine of leaf angle	dimensionless	0.45 a	0.31 a	0.45 a	0.73 a	0.73 a
	$hs/h$	stem fraction, stem height/tree height	dimensionless	0.40 b	0.40 b	0.20 b	0.00 b	0.15 b
	${\eta }_0$	stomatal leaf area/illuminated leaf area	dimensionless	2.50 c	1.00 c	1.00 c	1.00 c	1.00 c
	$m$	exponent relating shear stress on foliage to horizontal wind velocity	dimensionless	0.50 d	0.50 d	0.50 d	0.59 d	0.59 d
	$n$	number of sides of each foliage element producing surface resistance to wind	dimensionless	2.00 d	2.00 d	2.00 d	2.00 d	2.00 d

Note: The symbols and description of meteorological and vegetation data are from Eagleson [23]. “^a”, “^b”, “^c”, and “^d” represent the empirical value of the parameter, while the other values are the annual average values of the observed data from 2000 to 2018. “^a” from Rauner [46] and Ripley [47], “^b” from Cong [41], Rauner [46], and Zhang [42], “^c” from Waring [48], and “^d” from Amiro [49]. The above empirical value of cultural vegetation was obtained without considering seeding and harvesting.

and

Table 2. The empirical values of the soil parameters of five soil texture types.

Symbols	Description	Units	Values
			Loamy Clay	Clay Loam	Sandy Clay Loam	Loam	Sandy Loam
$n_e$	effective soil porosity	dimensionless	0.45	0.45	0.45	0.35	0.25
$d$	diffusivity index of soil	dimensionless	4.30	4.30	4.30	2.90	2.30
$\Psi \left(1\right)$	saturated matrix potential of soil	mm	900	900	900	450	250
${\rm K}\left(1\right)$	effective saturated hydraulic conductivity of soil	mm/day	29.40	29.40	29.40	29.40	2940
$m_s$	soil pore-size distribution index	dimensionless	0.44	0.44	0.44	1.20	3.30

Note: The symbols, description, and values of the soil parameters are from Eagleson [23].

) required to run the model include vegetation data, meteorological data, and soil data, among which the empirical parameters are derived from previous studies and the leaf area index (LAI) was obtained from the MCD15A2H (8-day 500-m) product. Multiple meteorological datasets were required to calculate the precipitation characteristic factors ($m_h$, $m_{\upsilon }$, $m_{t_r}$, and $m_{t_b}$) and the potential evaporation using the Penman–Monteith equation [45]. Daily precipitation, average air temperature, maximum air temperature, minimum air temperature, sunshine duration, wind speed, and relative humidity data from 2000 to 2018 were acquired from 66 meteorological stations (Figure 1a) across the study area by the China Meteorological Data Service Center (http://data.cma.cn accessed on 5 October 2019). The data on precipitation characteristics and potential evaporation were interpolated to the whole study area with ordinary kriging in ArcGIS software, and the spatial resolution was set to 1 km. In addition, the observed evapotranspiration (ET) data were used to verify the theory of ecohydrological optimality and were from the MOD16A2 (8-day 500-m) product.

The above remote-sensing data products are collection 6, which was obtained from NASA EARTHDATA (https://earthdata.nasa.gov/ accessed on 1 April 2020) and includes MOD13A3, MCD15A2H, and MOD16A2. These products were defined to the WGS84 geographic coordinate system and Albers projection coordinate system and resampled to 1 km by the MODIS Reprojection Tool (MRT). Briefly, simulated $M_0$ was calculated in each grid cell with an area of 1 km × 1 km, and a variety of parameters were set in light of the actual vegetation and soil texture types (Figure 1b,c). It was assumed that there was no energy or water exchange between grids. After removing urban and built-up lands, water bodies, wetlands, and barren areas, a total of 47,690 grids were calculated.

2.3. The Ecohydrological Optimality Model

On the basis of the statistical–dynamic water balance approach, the ecohydrological optimality theory was developed by Eagleson. Later, the water balance equation and potential canopy conductance equation were simultaneously solved, and thus, $M_0$ in equilibrium was obtained.

2.3.1. Water Balance Equation

The statistical–dynamical water balance model and the analytical solutions of its components were suggested by Eagleson with the approach of event-based derived distribution [23]:

$$P_{\tau }+m_{\mathsf{\tau }}w=m_{\upsilon }E\left[E_r\right]+m_{\upsilon }E\left[R_{sj}\right]+E\left[E_{T\tau }\right]+m_{\mathsf{\tau }}v+\Delta S,$$

(2)

where $E\left[\right]$ refers to the expected value of $\left[\right]$ and the subscript $\tau$ indicates the growing season. The first term to the left of the equal sign is precipitation ($P_{\tau }$, mm); the second term to the left is the capillary rise from the water table, where $m_{\mathsf{\tau }}$ is the length of the growing season, days, and $w$ is the rate of capillary rise from the water table to the surface, mm day⁻¹. The first term to the right is the surface retention capacity, where $m_{\upsilon }$ is the number of independent storms per growing season, time⁻¹, and $E_r$ is the surface retention depth of each independent storm, mm. The second term is the surface runoff, where $R_{sj}$ is the surface runoff of each storm, mm. The third term is ET, where $E_{T\tau }$ is the total evapotranspiration, mm. The fourth term is the percolation of soil water, where $v$ is the rate of percolation to the water table, mm day⁻¹; and the last term is the seasonal change in soil moisture storage from the beginning to the end of the growing season ($\Delta S$, mm). All quantities in Equation (2) are the annual average in the growing season that has the same length in this study. The detailed calculation process refers to Appendix B.

In Equation (2), $E\left[E_{T\tau }\right]$ consists of vegetation transpiration (represented by $m_{\upsilon }{E}_{ps\tau }{m}_{t_b^{\prime }}{M}_{\tau }{k}_v^{*}{\beta }_v$) and bare soil evaporation (represented by $m_{\upsilon }{E}_{ps\tau }{m}_{t_b^{″}}\left(1-M_{\tau }\right){\beta }_s$). Therefore,

$$E\left[E_{T\tau }\right]=m_{\upsilon }{E}_{ps\tau }\left[m_{t_b^{\prime }}{M}_{\tau }{k}_v^{*}{\beta }_v+m_{t_b^{″}}\left(1-M_{\tau }\right){\beta }_s\right],$$

(3)

where $k_v^{*}$ is the potential canopy conductance, dimensionless, which is identically equal to $E_{pv}/E_{ps}$, where $E_{pv}$ is the long-term average potential transpiration rate, mm day⁻¹, and $E_{ps}$ is the long-term average rate of the potential evaporation from the bare soil, mm day⁻¹; ${\beta }_v$ is the canopy transpiration efficiency, ${\beta }_v\equiv \frac{E_v}{E_{pv}}$, dimensionless, where $E_v$ is the long-term average transpiration rate, mm day⁻¹; ${\beta }_s$ is the bare soil evaporation efficiency, dimensionless; $m_{t_b^{\prime }}$ is the average interstorm time available for transpiration, days; and $m_{t_b^{″}}$ is the average interstorm time available for bare soil evaporation, days.

A state of maximum vegetation productivity, the critical moisture state, was proposed by Eagleson for the given species and the environment in terms of the water balance equation. Vegetation is assumed to be under initial stress at the end of the average growing season, and thus, the stomates are assumed to be fully open throughout the average interstorm period. We assume that when stomates are fully open, ${\beta }_v$ = 1. The maximum canopy moisture flux is

$$M_{\tau }{k}_v^{*}{E}_{ps\tau }=\frac{V_e}{m_{t_b^{\prime }}},$$

(4)

where the left-hand side of Equation (4) is the transpiration rate without water stress (the rate of vegetation water demand), and the right-hand side is the rate of available water supply for a specific environment. $M_{\tau }$ and $k_v^{*}$ are two vegetation state variables; $V_e$ is the water available for vegetation during the average interstorm period, mm, which is calculated from the water balance equation (the detailed calculation process can be found in Appendix B Equation (A1)). If the long-term average vegetation water demand is greater than the long-term average available water supply, the long-term average of $M_{\tau }$ will decrease to maintain the water balance.

2.3.2. Potential Canopy Conductance Equation

We assume that the bare soil is saturated after the storm. Potential evaporation from a saturated surface for the given atmospheric conditions was proposed by Penman [45]:

$$\lambda E_{ps}=\frac{\Delta R_n+\rho c_p\left[e_s-e_a\right]/r_a}{\Delta +{\gamma }_0},$$

(5)

while the vegetation transpiration is expressed by

$$E_v=\frac{\Delta R_n+\rho c_p\left[e_s-e_a\right]/r_a}{\Delta +{\gamma }_0\left(1+\frac{r_c}{r_a}\right)}$$

(6)

In this study, the stomates were assumed to always be fully open; hence $E_v=E_{pv}$. The potential canopy conductance is

$$k_v^{*}\equiv \frac{E_{pv}}{E_{ps}}=\frac{1+\frac{\Delta }{{\gamma }_0}}{1+\frac{\Delta }{{\gamma }_0}+\frac{r_c}{r_a}}=\frac{1+\frac{\Delta }{{\gamma }_0}}{1+\frac{\Delta }{{\gamma }_0}+\left(1-M\right){\left(\frac{r_c}{r_a}\right)}_{M\to 0}+M{\left(\frac{r_c}{r_a}\right)}_{M=1}},$$

(7)

where $\lambda$ is the latent heat of vaporization, MJ kg⁻¹;$R_n$ is the net solar radiation, MJ m⁻² day⁻¹; $\rho$ is the mean air density at constant pressure, kg m⁻³; $c_p$ is the specific heat at constant pressure, MJ kg⁻¹ °C⁻¹; $e_s$ is the saturation vapor pressure, kPa; $e_a$ is the actual vapor pressure, kPa; ${\left(\frac{r_c}{r_a}\right)}_{M\to 0}$ is the resistance ratio for open canopies, dimensionless; ${\left(\frac{r_c}{r_a}\right)}_{M=1}$ is the resistance ratio for closed canopies, dimensionless; $r_a$ is the lumped atmospheric resistance over 2 m above the canopy top, s mm⁻¹; and $r_c$ is the lumped resistance to flow through the canopy, s mm⁻¹.

A state of maximum vegetation productivity, the optimum foliage state, was additionally suggested by Eagleson (for the given species and the environment). To use light reasonably, the appropriate fixed values of the cosine of the leaf angle ($\beta$) and the foliage area index ($L_t$) were produced by natural selection; thus, the resistance ratio has its minimum when canopies are closed, and as a result, the $k_v^{*}$ value is at a maximum (the detailed calculation process can be found in Appendix B Equation (A3)) [50].

2.3.3. Ideal Vegetation Coverage

The $M_0$ required in this research is the vegetation coverage corresponding to the state of maximum productivity in the natural environment during the growing season. This state, also known as the optimum canopy state, must reach the critical moisture state and the optimum foliage state at the same time. Equations (4) and (7) are drawn in Figure 2. If

$$\frac{V_e}{m_{t_b^{\prime }}{E}_{ps\tau }}\ge \frac{1+\frac{\Delta }{{\gamma }_0}}{1+\frac{\Delta }{{\gamma }_0}+{\left(\frac{r_c}{r_a}\right)}_{M=1}},$$

(8)

the canopy is light-limited,$k_v^{*}=\frac{1+\frac{\Delta }{{\gamma }_0}}{1+\frac{\Delta }{{\gamma }_0}+{\left(\frac{r_c}{r_a}\right)}_{M=1}}$ and its $M_0$ = 1.

In contrast, if

$$\frac{V_e}{m_{t_b^{\prime }}{E}_{ps\tau }}<\frac{1+\Delta /{\gamma }_0}{1+\frac{\Delta }{{\gamma }_0}+{\left(\frac{r_c}{r_a}\right)}_{M=1}},$$

(9)

the canopy will be water limited, and $M_0$ is obtained by the simultaneous execution of Equations (4) and (7). At this time, the water supply of the system and the water demand of vegetation reach a dynamic equilibrium.

2.4. Multiple Linear Stepwise Regression Analysis

After checking for correlations between the independent variables, the bidirectional elimination stepwise method was utilized in the multiple linear regression analysis to screen variables and solve the problem of multicollinearity among variables. The bidirectional elimination stepwise method combines forward selection and backward elimination while considering the independent variables one by one [51]. The variance test and the minimum principle of the Akaike information criterion (AIC) were used whether selecting or eliminating variables. Each variable is standardized to obtain standardized partial regression coefficients and thus reflects the linear influence of each independent variable on the dependent variable when other independent variables are fixed.

2.5. Sensitivity Analysis

The relative sensitivity analysis method was used in this study [52,53].

$$C=\frac{{{\displaystyle \sum }}_{i=0}^{z-1}\left|\frac{\frac{Y_{i+1}-Y_i}{Y_0}}{\left(F_{i+1}-F_i\right)}\right|}{z},$$

(10)

where $C$ is the sensitivity coefficient; $Y_i$ is the output value of the ith operation of the model;$Y_0$ is the output value of the model that inputs the initial parameter value;$F_i$ is the percentage change in the parameter value of the ith operation of the model relative to the initial parameter value; and $z$ is the number of model runs.

3. Results

3.1. Spatial Distribution and Impact Factors of Observed Vegetation Coverage

The theory of ecohydrological optimality focuses on the long-term average status of the growing season rather than dynamic trends. From the perspective of the spatial distribution, the $M$ value showed horizontal zonal changes, with an average of 0.502 (Figure 3b). The average $M$ values of needleleaf forest, broadleaf forest, scrub, steppe, and cultural vegetation were 0.784 ± 0.176, 0.593 ± 0.241, 0.691 ± 0.212, 0.501 ± 0.171, and 0.456 ± 0.139, respectively, in the 19-year growing season. There was a gradual decrease in $M$ from the eastern edge of the mountainous area to the central section. The plain area is mainly planted with cultural vegetation, which has low vegetation cover.

The correlation of eight meteorological variables was analyzed and then the highly correlated variables were manually removed (here, high correlation means that the correlation coefficient is higher than 0.80 and the results of correlation can be found in the “Supplementary Materials”). We found that the number of independent storms ($m_{\upsilon }$), the mean time between storms ($m_{t_b}$), and the mean daily air temperature ($T_{\tau }$) were highly correlated except steppe. $m_{\upsilon }$ and $T_{\tau }$ were thus excluded due to our primary focus being the response of vegetation to available water during the average interstorm period. Next, the bidirectional elimination multiple linear stepwise regression was used to analyze the influence of meteorological variables on the average growing season $M$ for the five vegetation types from 2000 to 2018 (

Table 3. Standardized partial regression coefficients of multiple linear stepwise regression between the observed vegetation coverage ($M$ ) value and meteorological variables for the five vegetation types.

Meteorological Variables	Vegetation Types
	Needleleaf Forest	Broadleaf Forest	Scrub	Steppe	Cultural Vegetation
$m_h$	0.5515	0.5256	0.2060	0.1305	0.7294
$m_{\upsilon }$	-	-	-	-	-
$m_{t_r}$	0.1907	0.2778	0.02631	0.05210	0.04273
$m_{t_b}$	−0.1144	−0.1214	−0.05810	−0.1769	−0.2335
$T_{\tau }$	-	-	-	−0.2071	-
$\mu$	-	−0.1674	−0.1687	−0.2834	−0.06677
$RH$	-	-	0.3956	0.5016	-
$N$	-	−0.08726	-	-	-

Note: “-” indicates that the variable corresponding to the regression coefficient was eliminated based on the bidirectional elimination stepwise method. All the selected variables and regression equations passed the significance test at a significance level of 0.001. The adjusted coefficients of determination (adjR²) of the regression equations are 0.2864, 0.5907, 0.5389, 0.3376, and 0.3419 for the five vegetation types, respectively.

). The standardized partial regression coefficient results suggest that the $M$ values of all the vegetation types were significantly affected by the mean storm depth ($m_h$), the mean storm duration ($m_{t_r}$), and $m_{t_b}$, of which $m_h$ is the most significant variable. For the needleleaf forest, the $M$ value was only significantly affected by the precipitation characteristic factors. For the other vegetation types, the mean daily wind speed ($\mu$) and the mean daily relative humidity ($RH$) are also important to the $M$ value. It is worth noting that the $M$ value of the steppe was negatively correlated with $T_{\tau }$. Generally, the associations between the precipitation characteristic factors and the $M$ value are stronger than other variables and the $M$ value on a long-term average time scale.

3.2. Verification of Ecohydrological Optimality Theory

To verify the theory of ecohydrological optimality, the $k_v^{*}$ value corresponding to the actual vegetation coverage is obtained by inputting the $M$ data into Equation (7). Then, the $k_v^{*}$ and $M$ values are substituted into Equation (4) and combined with the water balance equation to obtain the simulated ET. Figure 4a,b presents similar spatial distributions of the average simulated ET data and average observed ET data during the growing season from 2000 to 2018. The simulation results and the observation results are drawn in the scatter plot and fitted by a linear function. As illustrated in Figure 4d, the Nash–Sutcliffe efficiency coefficient (NSE) and the coefficient of determination (R²) demonstrate that the ecohydrological optimality theory is applicable to this watershed in general, with the NSE value being 0.316 and the R² value being 0.770.

Interestingly, the slope and intercept of the linear function are 1.300 ± 0.003 and −43.084 ± 0.729, respectively. The simulated ET is slightly higher than the observed ET. With successive increases in the observed ET, the degree of deviation of the simulated ET increases. The root mean square error (RMSE) of the scatter plot is 58.280, which indicates that there are marked differences between the simulated ET and the observed ET in some areas. The spatial distribution of these differences is provided in Figure 4c. Compared with the vegetation type distribution (Figure 1b), it can be seen that vegetation types with large differences between the simulated ET and the observed ET are forest (including needleleaf forest and broadleaf forest) and scrub, with high vegetation coverage.

3.3. Simulated Vegetation Coverage Based on the Ecohydrological Optimality Model

The simulated $M_0$ under the natural environment that a specific climate and soil can maintain was calculated based upon the meteorological, soil, and vegetation data of the watershed for 19 years, with an average of 0.475 (Figure 3a, 2000–2018). The average $M_0$ values of needleleaf forest, broadleaf forest, scrub, steppe, and cultural vegetation in the growing season were 0.433 ± 0.066, 0.395 ± 0.089, 0.498 ± 0.090, 0.548 ± 0.116, and 0.458 ± 0.100, respectively, from 2000 to 2018. The simulated values were highest in the north and east of the study area and lowest in the central and southwest regions, where the simulated values were roughly similar to the observed values. $M$ exhibited a wide range, from 0.058 to 1, while $M_0$ exhibited a range from 0.212 to 0.949. The differences between the $M$ and $M_0$ values increase from the northwest to the southeast of this area, with a maximum value exceeding 0.3 (Figure 3c).

3.4. Sensitivity of the Ideal Vegetation Coverage to the Input Data of the Model

The input parameters of this model were divided into physical parameters ($m_h$, $m_v$, $m_{t_r}$, $m_{t_b}$, $T_{\tau }$, $\mu$, $RH$, $N$, $P_d$, $E_{psd}$, $L_t$, $M_d$) and empirical parameters ($\beta$, $h_s/h$, ${\eta }_0$, $m$, $n$, $n_e$, $d$, $\psi \left(1\right)$,$K\left(1\right)$, $m_s$). The relative sensitivity analysis method facilitated the screening of insensitive parameters in the model. Therefore, the efficiency of the model operation and validation will be improved. When the absolute value of $C$ was less than 0.05, the parameters were considered insensitive and are not shown in Figure 5. The $M_0$ was not sensitive to soil conditions if the whole watershed was simulated as sandy loam soil (Figure 5a). However, the results of the model were sensitive to all parameters if the soil texture was loam (Figure 5b). When the soil texture was simulated as loamy clay, clay loam, or sandy clay loam, the outputs were sensitive to only effective saturated hydraulic conductivity of soil ($K\left(1\right)$) among all soil parameters (Figure 5c). In all situations, $L_t$ was an insensitive parameter, which implied that it may be displaced by an empirical value in future research. Remarkably, the reason the $C$ value of the stem fraction ($h_s/h$) of the steppe could not be calculated by this method is that the initial value of this parameter was 0.00. Taken together, these results indicated that the accuracy of meteorological data primarily affected the model results.

The relative sensitivity analysis method can only qualitatively analyze the input parameters and produce the relative sensitivity ranking of the model parameters. To compare the sensitivities of $M_0$ at disparate vegetation and soil texture types to the same parameter, further quantitative analysis is needed in accordance with the actual range of the parameters. The observed range of parameters $\beta$ and $h_s/h$ came from previous studies [41,42,46]. As shown in Figure 6 and Figure 7, the variation range of $M_0$ is quite narrow when $\beta$ and $h_s/h$ change in terms of their actual range. $M_0$ increases with increasing $\beta$ or decreasing $h_s/h$. In addition, the change in parameter $h_s/h$ changes only the water demand curve, which has nothing to do with soil properties. The results illustrated that the sensitivity of parameter $h_s/h$ is affected by only vegetation type, and the sensitivity of parameter $\beta$ is mainly affected by soil texture. $M_0$ is more sensitive to $\beta$ when the soil texture of the study area is simulated as sandy loam instead of other types. The $M_0$ of the steppe is more sensitive to $h_s/h$ than that of the other types.

To shed more light on the response of $M_0$ to changes in meteorological factors, it is necessary to determine the threshold values of these factors over the years. Obviously, the $m_h$, $m_{t_b}$, and $T_{\tau }$ values have a great influence on the $M_0$ value (Figure 5). $M_0$ increases as $m_h$ increases and $m_{t_b}$ and $T_{\tau }$ decrease (Figure 8). When the parameter $m_h$ changes in the same range, the variation value of $M_0$ for sandy loam soil is almost twice that for other soil types. Meanwhile, the $M_0$ of steppe changes more than that of scrub and broadleaf forest. When the parameter $m_{t_b}$ changes in the same range, the fluctuation of $M_0$ for scrub is slightly larger than that for steppe and broadleaf forest. The response ranges of dissimilar vegetation types to temperature sorted by size are steppe, scrub, and broadleaf forest. In summary, the response of $M_0$ to meteorological factors fluctuates most intensively if the soil type is sandy loam. $m_h$ and $T_{\tau }$ strongly affected the $M_0$ of the steppe, and the change in $m_{t_b}$ had the greatest influence on the $M_0$ of scrub.

4. Discussion

In this study, we analyzed the relationship between meteorological factors and the vegetation distribution pattern, evaluated the applicability of the ecohydrological optimality theory in the Yongding River watershed, and calculated the $M_0$ value in this area based on the local climatic–soil–vegetation conditions. Our observations showed that $M$ had the greatest correlation with the precipitation characteristic factors. The annual average $M$ value in the growing season generally showed horizontal zonal changes. Compared with the simulated $M_0$ values, there were several differences between $M_0$ and $M$ in some areas. To use the model more conveniently, we analyzed the sensitivity of the input parameters as well. Collectively, our results indicate that the ecohydrological optimality model is suitable to simulate $M_0$ in this watershed, and these results could provide reference indicators for ecological restoration.

4.1. Causes for the Distribution Pattern of Vegetation

A large number of studies have shown that vegetation coverage is affected by human activities, precipitation, temperature, and soil moisture conditions [54,55,56]. Analysis of $M$ and meteorological factors enables us to speculate the causes of the distribution pattern of vegetation. Hydrothermal conditions directly affect the spatial distribution of $M$ [57]. On the one hand, increases in the air temperature and radiation will increase evaporation rates, resulting in drier soils, which impede the growth of vegetation [58]. On the other hand, increases in the air temperature and radiation promote photosynthesis, which accelerates the increase in vegetation cover [59]. The steppe was negatively correlated with $T_{\tau }$, possibly because increases in $T_{\tau }$ affect the shallow root community in areas such as the steppe to absorb water from the surface of the soil. The heat conditions of the study area are sufficient for steppe growth, which could not offset the negative effect of soil drying [60]. Secondly, the distinct effects of meteorological factors on the vegetation distribution pattern are associated with the vegetation type [61,62]. The mean sunshine duration per day ($N$) and $\mu$ significantly affected the broadleaf forest, but not the needleleaf forest. Changes in the needleleaf forest are more dependent on precipitation than ET due to their smaller leaf area and thicker leaf blade, which prevents heat accumulation and water evaporation on the leaf surface [63]. In addition, the variability in the precipitation characteristic factors, especially $m_h$ and $m_{t_b}$, changes the spatiotemporal availability of soil water uptake by plants and has a significant impact on vegetation [64]. Zhang [65] found that the number of rainy days, heavy-rainfall events, and consecutive dry days are three important indicators affecting the NDVI in the growing season across the Sahel. Although the multiple stepwise regression quantified the linear impact of the meteorological variables on $M$, further mathematical and physical deduction is needed for mechanism analysis.

Geographical factors indirectly affect the spatial distribution of $M$ [61,66]. Eastern and southeastern China face the Pacific Ocean, and the amounts of precipitation and humidity decrease from southeast to northwest. Nevertheless, a part of the Tai-hang Mountains is located on the southeast edge of the mountainous area of the Yongding River watershed. The high-altitude mountains influence the distribution of not only hydrothermal conditions but also soil moisture and nutrient conditions [67]. As a consequence, the $M$ of the southeast slope is higher than that of the northwest slope of the Tai-hang Mountains.

Furthermore, the average $M$ value of the steppe and cultural vegetation was lower than that of the other vegetation types in the growing season. This difference may be related to human activities and the drought escape mechanism of vegetation [68]. Because of the undeveloped xylem, herbaceous plants have a low water storage capacity and are sensitive to climate change [69]. Some plants will complete their life cycle ahead of water deficits [70]. Some types of cultural vegetation have the ability to avoid drought by delaying ripening [71]. In addition, this difference is also related to the calculation principle of NDVI. The chlorophyll content and the photosynthetic capacity of the forest were greater than those of the steppe, reflecting the presence of less visible red light, more near infrared light, and a higher NDVI.

4.2. Applicability of the Ecohydrological Optimality Theory

We introduced the MOD16A2 products from NASA to validate the theory. The MOD16A2 data have been verified by the observed latent heat flux from 232 global watersheds and 46 field-based eddy covariance flux towers, and the results, the observed ET, should be considered robust [72]. The annual average $M$ value (Figure 3b) in the growing season is the value obtained under actual local water supply conditions, which is the performance of actual water consumption of vegetation. Similarly, the simulated ET was calculated under the same conditions based on the ecohydrological optimality theory. Consequently, the results presented here can be used to demonstrate the applicability of this theory for cultural vegetation and other types of vegetation (Figure 4). Our findings are consistent with those of Zhang [42]. Conversely, the water supply condition calculated by the water balance equation is under a natural environment without interference (Figure 3a). The simulated $M_0$ obtained by the model is an ideal value. The correlation between it and the $M$ may be interfered with by human beings.

One unanticipated finding was that the average simulated ET calculated by substituting $M$ in the model was higher than the average observed ET in the region with forest and scrub during the growing season for 19 years. This result implies that the vegetation did not reduce $M$ to avoid water stress but did reduce the transpiration amount in these areas on the scale of 19 years. In other words, this result contradicts the hypothesis that transpiration always occurs at the rate without stress until the vegetation is under initial stress just at the end of the average growing season. Kerkhoff [43] proposed that this seems to be a critical oversight that can cause the transpiration simulated by models to exceed the amount of water available to vegetation. This phenomenon may be related to the simulation process in which the actual situation is replaced by the long-term average value. In a short period of time (within one drought episode), for both trees and herbaceous plants, physiological regulation, such as adjusting the stomata to cope with the situation that the soil moisture is not enough to maintain the maximum transpiration rate in the later interstorm period, will occur in vegetation [73,74]. However, for long time scales, the selection pressure of productivity and the competitive pressure for water allow the types of vegetation that close the stomata much later to survive [75]. Diversified vegetation has different feedback on the frequency of stomatal closure. Most herbaceous plants cannot tolerate premature stomatal closure or frequent stomatal closure, and thus a low vegetation coverage is selected to reduce the demand for water. On the contrary, some woody plants can tolerate or avoid intermittent short-term perturbations of water stress and thereby prevent reductions in vegetation coverage [76,77]. Specifically, longleaf pine responds to extreme drought by better controlling its stomata to limit water loss, and transpiration recovers rapidly during the next precipitation event [78]. The leaves of Robinia pseudoacacia increase the elasticity of the cell wall and veins and the stomatal density to adapt to long-term water stress [15]. Obviously, this long-term, frequent resistance of woody plants to water stress was not considered by Eagleson. Therefore, the $M_0$ value calculated by ecohydrological optimality theory may be lower than the $M$ value in some areas of forest and scrub, mainly in areas with a high vegetation coverage.

In the long run, the optimal vegetation coverage proposed in this paper is not the maximum vegetation coverage under the same water supply conditions, but the vegetation coverage obtained by minimizing water stress and maximizing productivity. Although trees can tolerate a certain frequency of hydraulic failure, such events may cause embolism in the xylem or roots [79]. The mechanism of drought tolerance or drought avoidance of trees will gradually degenerate as the number of droughts increases [80]. Some research has pointed out that decreases in hydraulic function and photosynthesis for long durations may increase the risk of death after many years [81]. In summary, ecohydrological optimality theory proposes a kind of vegetation coverage that can maintain ecosystem stability.

4.3. Comparison between Simulated and Observed Vegetation Coverage

Our results indicated that the growth of vegetation in this watershed under a natural environment was limited by water, not light, because $M_0$ < 1. The comparison of $M_0$ and $M$ contributes to the more effective development of ecological restoration measures. The simulated $M_0$ was calculated according to the local hydrological cycle, without considering the artificial water supply and consumption. Although built-up lands with intense human activity were not the focus of this research, artificial irrigation, planting, or fertilization of both forest land and farmland occur locally. This condition is one of the reasons $M$ is larger than the simulated $M_0$ in some areas. In addition, we found that the simulated values along the southeastern edge of the study area were not as high as the observed values. This difference may be caused by the shortage in the number and uneven distribution of meteorological stations, which fails to simulate the actual differences in the hydrothermal distribution due to the impacts of topography. Interpolation of meteorological stations also means that the simulated values are not as detailed as the observed values for the description of the dynamic distribution of vegetation. In general, combined with Figure 4c, when $M$ is obviously larger than $M_0$ and the observed ET is similar to the simulated ET, it implies that the vegetation has an extra water supply or artificial cultivation and is not under water stress. However, if $M$ is clearly larger than $M_0$ and the simulated ET is much larger than the observed ET, regardless of whether or not there is extra water supply, the vegetation will suffer from water stress and evolve drought resistance measures.

In the practical process of ecological restoration, we should pay more attention to the quality of vegetation rather than pursue high vegetation coverage. Vegetation restoration should be carried out as the water capacity permits. Excessive vegetation restoration not only keeps vegetation under water stress but also threatens water resources. As the limit of vegetation restoration, $M_0$ should be classified and discussed in this region. Most needleleaf forests and broadleaf forests in this region have reached $M_0$. At this time, we should strengthen the management and protection of existing trees and adopt measures such as forest tending, pruning, and thinning to ensure the growth quality of vegetation [80]. For cultural vegetation, we should use water-saving irrigation technologies and adjust the structure of agricultural planting. For the steppe, we should use fencing and strengthen management measures to ensure natural growth [82].

4.4. Sensitivity of the Ideal Vegetation Coverage to Meteorological, Soil, and Vegetation Factors

Diverse vegetation types and soil texture should be discussed separately in sensitivity analyses because the vegetation and soil parameters are set to different values derived from the actual situation of each grid cell when the model runs. The soil texture mainly affects surface runoff, which results in the distinction of the water available for vegetation. Vegetation type primarily affects the water demand capacity of vegetation. The results demonstrated that the ideal coverage of vegetation planted in sandy loam soil is highly sensitive to changes in meteorological and vegetation factors. This result could be attributed to the runoff calculated from soil properties being extremely small and more water being supplied to the vegetation in this case; thus, soil properties have little effect on the model. The variation in $\beta$ changed both the water supply curve and the water demand curve. The less surface runoff there is, the more important the surface retention capacity determined by $\beta$ is. The higher the value of $\beta$ is, the lower the value of the water demand curve.

These findings suggest that it is better to measure the soil parameters on site if the soil texture is loam or utilize weather stations with high precision to measure the meteorological factors if the soil texture is sandy loam to improve the accuracy of simulations. In addition, scrub vegetation with larger $\beta$ and smaller $h_s/h$ values should be considered in place of forests that are under water stress to increase the vegetation coverage or save water resources.

5. Conclusions

On a long-term average timescale, the optimal coverage of vegetation in the Yongding River watershed is limited by water rather than light. In this study, the spatial patterns and impact factors of vegetation were explored at the watershed scale. Meanwhile, we verified the applicability of the ecohydrological optimality theory in this region and calculated the $M_0$ value in terms of the local meteorological, soil, and vegetation factors from 2000 to 2018. Finally, the sensitivity of $M_0$ to the input parameters of the model under disparate conditions of vegetation and soil was also considered. Additionally, a large number of studies with long-term observations of precipitation, runoff, transpiration, soil moisture, and root depth of vegetation are needed in future work. The main findings obtained from this research are as follows:

• From the spatial point of view, the holistic $M$ exhibits horizontal zonal changes that basically correspond to topographic factors, which gradually decrease from the eastern edge of the mountainous area to the central section. Temperature increases limited the water needed by the steppe. The $m_h$, $m_{t_r}$, and $m_{t_b}$ values had the highest correlation with the $M$ value, as they controlled the spatial and temporal distribution of the water available for vegetation.
• The ecohydrological optimality model simulates a kind of vegetation coverage that maximizes the productivity of biomass in the absence of water stress rather than the maximum vegetation coverage, as this model cannot completely simulate the adaptive strategies of vegetation to water stress. The vegetation coverage proposed by ecohydrological optimality theory can reduce the risk of embolism in the xylem or roots or vegetation death in the years following drought events, which is significant for maintaining the stability of ecosystems.
• The $M_0$ should be used as the limit of vegetation restoration, regardless of whether the vegetation evolves drought resistance measures. For forests that have reached $M_0$, we should strengthen the management and improve the quality of existing vegetation. For cultural vegetation that has not reached the $M_0$ value, we should use water-saving irrigation technologies and adjust the structure of agricultural planting. For the steppe, we should use fencing and strengthen management measures to ensure natural growth.
• In this model, the accuracy of meteorological data is essential. The $M_0$ value was the most sensitive to the meteorological factors and $\beta$ when the soil texture of the study area was simulated as sandy loam instead of the other types. The $M_0$ of the steppe is most sensitive to $h_s/h$, $m_h$, and $T_{\tau }$.