Integrating Remotely Sensed Leaf Area Index with Biome-BGC to Quantify the Impact of Land Use/Land Cover Change on Water Retention in Beijing

: Maintaining or increasing water retention in ecosystems (WRE) can reduce ﬂoods and increase water resource provision. However, few studies have taken the effect of the spatial information of vegetation structure into consideration when assessing the effects of land use/land cover (LULC) change on WRE. In this study, we integrated the remotely sensed leaf area index (LAI) into the ecosystem process-based Biome-BGC model to analyse the impact of LULC change on the WRE of Beijing between 2000 and 2015. Our results show that the volume of WRE increased by approximately 8.58 million m 3 in 2015 as compared with 2000. The volume of WRE in forests increased by approximately 26.74 million m 3 , while urbanization, cropland expansion and deforestation caused the volume of WRE to decline by 11.96 million m 3 , 5.86 million m 3 and 3.20 million m 3 , respectively. The increased WRE contributed by unchanged forests (14.46 million m 3 ) was much greater than that of new-planted forests (12.28 million m 3 ), but the increase in WRE capacity per unit area in new-planted forests (124.69 ± 14.30 m 3 /ha) was almost tenfold greater than that of unchanged forests (15.60 ± 7.85 m 3 /ha). The greater increase in WRE capacity in increased forests than that of unchanged forests was mostly due to the fact that the higher LAI in unchanged forests induced more evapotranspiration to exhaust more water. Meanwhile, the inverted U-shape relationship that existed between the forest LAI and WRE implied that continued increased LAI in forests probably caused the WRE decline. This study demonstrates that integrating remotely sensed LAI with the Biome-BGC model is feasible for capturing the impact of LULC change with the spatial information of vegetation structure on WRE and reduces uncertainty.


Introduction
Water retention in ecosystems (WRE), also called water conservation, includes the rainwater retained by the canopy, litter and soil of ecosystems to satisfy the water resource demands of the upstream area, but also the downstream area or external areas [1,2]. WRE is a sponge effect that regulates rainfall and releases it more slowly to supply baseflow. WRE is an important ecosystem service, and the magnitude and distribution of WRE are the composite results of climate, terrain, soil, land use/land cover (LULC) and vegetation structure [1,3]. Nowadays, the impacts of LULC changes on WRE are widely reported [4][5][6]. However, few studies have focused on the impact of LULC change on WRE with consideration of vegetation structures. Owing to the acute and frequent LULC changes that Beijing is located in northern China (Figure 1a), and it is also the political, economic, and cultural centre of China. Beijing is in a semi-humid area with an annual mean temperature of 3-15 • C and annual precipitation of 480-640 mm. Mountainous areas account for 62% of the total area of Beijing and the elevation decreases from the northwest to southeast. Forest areas account for approximately 50.97% of the total area (Figure 1b). Cropland and urban areas account for approximately 21.43% and 20.32%, respectively (Figure 1b). The wetland area is the smallest and accounts for only approximately 1.90% (Figure 1b). In past decades, Beijing has experienced serious sandstorms due to vegetation degradation [51]. To alleviate the risk of sandstorms, the Chinese government started the Beijing-Tianjin Sand Source Control Project in 2000 [52]. It is a restoration project which aims to plant more trees and grass to protect soil and reduce dust. In addition, the government implemented the plains afforestation project to improve the quality of the environment and landscape [53,54]. Meanwhile, the urban area has increased quickly due to rapid economic and population growth. All these factors collectively caused dramatic LULC changes and varied vegetation structures.

Study Area
Beijing is located in northern China (Figure 1a), and it is also the political, economic, and cultural centre of China. Beijing is in a semi-humid area with an annual mean temperature of 3-15 °C and annual precipitation of 480-640 mm. Mountainous areas account for 62% of the total area of Beijing and the elevation decreases from the northwest to southeast. Forest areas account for approximately 50.97% of the total area (Figure 1b). Cropland and urban areas account for approximately 21.43% and 20.32%, respectively (Figure 1b). The wetland area is the smallest and accounts for only approximately 1.90% (Figure 1b). In past decades, Beijing has experienced serious sandstorms due to vegetation degradation [51]. To alleviate the risk of sandstorms, the Chinese government started the Beijing-Tianjin Sand Source Control Project in 2000 [52]. It is a restoration project which aims to plant more trees and grass to protect soil and reduce dust. In addition, the government implemented the plains afforestation project to improve the quality of the environment and landscape [53,54]. Meanwhile, the urban area has increased quickly due to rapid economic and population growth. All these factors collectively caused dramatic LULC changes and varied vegetation structures.

Biome-BGC Model
Biome-BGC is an ecosystem process model that simulates biogeochemical and hydrological processes among multiple biomes [55]. This model is driven at a daily time scale with seven daily meteorological values (maximum temperature, minimum temperature, precipitation, vapour pressure deficit (VPD), shortwave radiation and day length). Biome-BGC is a self-equilibrium model that needs spin-up to produce the equilibrium state of the carbon, nitrogen and water states of vegetation, litter and soil pools that are to be saved in its restart files (files that provide initialization parameters used in the normal run). Spin-up is a self-initialization procedure that uses the available climate records to perform repeated simulations until a steady state is reached [56]. Biome-BGC and its modified version have proven to be effective in studies of global or regional water, carbon and nitrogen [56][57][58]. Previous studies have combined remotely sensed data to obtain more accurate phenology information to promote productivity estimation in Biome-BGC [59,60]. Others also combine remotely sensed data or ground-measured data with Biome-BGC to correct soil temperature [61,62], root depth [63,64] and other parameters [65,66] to improve its performance. The hydrological processes in Biome-BGC includes precipitation interception, infiltration, stormwater runoff (SR) and ET. ET is calculated by the Penman-Monteith equation [67] and is separated into canopy evaporation, transpiration and soil evaporation. Canopy evaporation and soil evaporation are calculated from the water

Biome-BGC Model
Biome-BGC is an ecosystem process model that simulates biogeochemical and hydrological processes among multiple biomes [55]. This model is driven at a daily time scale with seven daily meteorological values (maximum temperature, minimum temperature, precipitation, vapour pressure deficit (VPD), shortwave radiation and day length). Biome-BGC is a self-equilibrium model that needs spin-up to produce the equilibrium state of the carbon, nitrogen and water states of vegetation, litter and soil pools that are to be saved in its restart files (files that provide initialization parameters used in the normal run). Spin-up is a self-initialization procedure that uses the available climate records to perform repeated simulations until a steady state is reached [56]. Biome-BGC and its modified version have proven to be effective in studies of global or regional water, carbon and nitrogen [56][57][58]. Previous studies have combined remotely sensed data to obtain more accurate phenology information to promote productivity estimation in Biome-BGC [59,60]. Others also combine remotely sensed data or ground-measured data with Biome-BGC to correct soil temperature [61,62], root depth [63,64] and other parameters [65,66] to improve its performance. The hydrological processes in Biome-BGC includes precipitation interception, infiltration, stormwater runoff (SR) and ET. ET is calculated by the Penman-Monteith equation [67] and is separated into canopy evaporation, transpiration and soil evaporation. Canopy evaporation and soil evaporation are calculated from the water intercepted by the canopy and the water saved in the soil. Transpiration is controlled by conductance under the impact of the meteorological conditions of the minimum air temperature, VPD, solar radiation, and soil water availability [68,69]. SR occurred under two cases: one case is that incoming precipitation induced the soil water content to exceed the saturated soil water capacity, and SR was calculated as the difference between the saturated soil water capacity and the remaining water over capacity. The other case is that while soil water content exceeds field capacity but is lower than saturated soil water capacity, the outflow occurs at a decaying rate with half of the soil water that is over soil water capacity [70]. In addition, litter has been shown to hold large amounts of water and decrease surface runoff [71,72]. Ignoring it may lead to an underestimation of infiltration and an overestimation of SR. Thus, we added a litter hydrological module to Biome-BGC. Details about this are described in the Supporting Information. In addition, to promote its performance in wetland, we enlarged soil depth to increase potential soil water content, which caused evaporation to be restricted only by the availability of radiation rather than water [73]. WRE in urban areas was set to zero owing to these areas being dominated by impervious surfaces.

Datasets
The acquired data used to drive the Biome-BGC model are listed in Table 1, and explicit details about the data are introduced in the following parts. The meteorological data were obtained from the National Meteorological Centre (http: //www.data.cma.cn/ (accessed on 20 July 2021)) and were based on 31 site observations ( Figure S3a). They included the daily maximum temperature, daily minimum temperature, and daily precipitation. The inverse distance weighted interpolation method provided by ArcGIS 10.2 was used to interpolate these meteorological variables for the whole of Beijing. Meanwhile, we also assessed the confidence of the interpolation by comparing the interpolated gridded precipitation with observed precipitation at the daily scale. The comparison showed a highly consistent relationship between them, implying the reliability of our interpolated gridded precipitation (R 2 = 0.98, Figure S3b). Next, elevation, latitude and these variables were input into MT-CLIM software to produce the daily VPD, average daytime temperature, daily shortwave radiation, and day length data [74]. Finally, all of these meteorological variables were prepared as the input meteorological variable for the BIOME-BGC model.

LULC and LAI
LULC data were derived from the China Ecosystem Assessment and Eco-security Database (http://www.ecosystem.csdb.cn (accessed on 18 July 2021)). The image classifications were derived from Landsat TM/ETM+ and OLI images using the object-oriented multi-scale segmentation and decision tree procedures. The algorithm for hierarchical classification was employed to effectively control classification quality. The independent validation indicated that the overall accuracy reached 86% [75]. The LULCs were categorized into forest, grassland, wetland, cropland and urban (Table S1). The LAI data are an 8-day product with a spatial resolution of about 500 m, representing the MCD15A2H V6 data derived from Terra MODIS, and can be downloaded from the website (https://earthexplorer.usgs.gov/ (accessed on 22 July 2021)).

Elevation and Soil Data
The elevation data used in this study were the SRTM DEM data. These data are produced by NASA with a spatial resolution of 90 m and can be downloaded from the USGS website (https://earthexplorer.usgs.gov/ (accessed on 22 July 2021)). The soil data included the sand (%), silt (%), and clay (%) contents, which were obtained from the Harmonized World Soil Database v 1.2 with a spatial resolution of about 30 arc seconds. These data can be downloaded from the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/ (accessed on 22 July 2021)).

Framework of Integrating RS Time Series in Biome-BGC to Estimate WRE
In the original Biome-BGC model, the LAI was estimated from the model itself by multiplying leaf C (kg/m 2 C) by specific leaf area (SLA) (m 2 /kg C) [76]. Leaf C is the total dry leaf mass C per unit area, while SLA is the leaf area per unit of dry leaf mass C [76,77]. However, the LAI estimated by Biome-BGC does not consider the vegetation structure of varied LULC change trajectories and disturbances (fire, pests and diseases) [55]. Thus, we used the remotely sensed LAI time series derived from MOD15A2H V6 as a surrogate. First, a modified temporal spatial filter was used to fill the gaps and replace lower-quality records with higher-quality records according to the quality control and filled value information [78]. Then, a Savitzky-Golay filter with a time window of 7 was used to smooth the LAI time series [79][80][81]. Finally, we interpolated the smoothed 8-day LAI time series into a daily LAI of one year (total 365 days) to meet the input needs of Biome-BGC.
Next, Equation (1) was used to calculate WRE based on water balance theory, and all required variables were the output variables from Biome-BGC ( Figure 2): where PRE k is the precipitation on day k; ET k is the total evapotranspiration on day k; SR k is the stormwater runoff on day k; CE k is the canopy evaporation on day k; T k is the transpiration on day k; SE k is the soil evaporation on day k; LE k is the litter evaporation on day k; SW k is the soil water content on day k; SW sat is the saturated soil water capacity; and SW f c is the field soil water capacity. To test the confidence of the simulated WRE, we compared our simulated ET with Multiple datasets were used as inputs to drive the Biome-BGC. Meanwhile, these datasets were varied in format, and spatial and temporal resolution. Thus, they were integrated to meet the needs of model input. First, all datasets were converted into raster grids via WGS84/Albers Equal Area Conic projection. Then, they were regrouped and resampled into a 1 km × 1 km grid for matching. Next, the spin-up run was used to produce the initial values of the water-state variables, which were saved in the restart files. Then, the normal run was used to model the hydrological process and export the hydrological variables. Finally, these output hydrological variables were used to calculate WRE based on Equations (1) and (2). In addition, we calibrated the eco-physiological parameters with local field data of prior studies in China, which would improve the performance of the Biome-BGC model in hydrological processes simulation. All eco-physiological parameters were obtained from the published literature and are given in Table S2.

Evaluating the Effect of Land Cover Changes on WRE 2.5.1. Model Performance Evaluation
To test the confidence of the simulated WRE, we compared our simulated ET with the ET results from the MOD16A2 product derived from Terra MODIS (https://earthexplorer. usgs.gov/ (accessed on 25 July 2021)) and the results from the CLSM (Catchment Land Surface Model) [82], VIC (Variable Infiltration Capacity Model) [83] and NOAH [84] derived from the Global Land Data Assimilation System (GLDAS) (https://ldas.gsfc.nasa. gov/gldas (accessed on 25 July 2021)) ( Figure S2a). Meanwhile, we also compared the simulated SR we estimated and that estimated based on the methods of Gong et al. [1] and Hamel et al. [17]. The comparisons showed a good consistency among them ( Figure S2b). In addition, we also compared the WRE capacity we estimated with other estimations that were conducted in the same climate region. All comparisons showed good consistency, indicating the feasibility of our results ( Figure S2c).

Data Analysis
(1) Quantifying LULC changes and their impact on WRE Here, we created the LULC change map by overlaying the LULC map of 2015 onto the LULC map of 2000 for each pixel cell to track the LULC change trajectories ( Figure 3). Meanwhile, we summarized LULC conversion categories, as shown in Table 2, for analysis of the results. To reveal the difference in vegetation structure across LULCs undergoing different LULC changes, we collect the annual maximum LAI of all pixels in each LULC. Then, we mapped the frequency distribution of annual maximum LAI in each LULC. Here, we used the maximum LAI rather than the mean or median LAI to highlight the size gap of LAI across LULCs with varied LULC change trajectories. In addition, we also added the median value line to better describe the difference in vegetation structure across different LULCs.
To quantify the impact of LULC change on WRE, we used the WRE map from 2015 minus the WRE map from 2000 to generate a map of WRE differences. Then, we used the LULC change map produced above to overlap the WRE difference map to compute the volume of WRE and WRE capacity across each LULC.
(2) Determining factors that affect WRE based on random forest (RF) algorithm The RF model is an ensemble learning method for regression based on constructing a multitude of decision trees to infer relationships between explanatory and response variables [85]. It has been widely used to address the influential predictors in large predictor sets in ecology-related studies [86][87][88]. Here, RF was used to identify the relative importance of different environmental factors in determining WRE. WRE obtained from each pixel was used as the target variable, while the maximum LAI, clay content (%), litter content (g/m 2 ), annual total precipitation (mm) and slope calculated from DEM ( • ) were obtained from each pixel and were used as input variables. All these input variables were chosen by referring to prior studies that demonstrate which factors connect tightly with WRE [89,90].
Next, the 'feature importance' produced by RF was used to quantify the relative importance of different factors to WRE.
(3) Using loess regression to fit a smooth curve between LAI and WRE/ET/SR Loess regression is a regression method that creates a smooth line through a scatter plot, displaying the relationship between variables and helping readers to foresee trends [91]. Compared with other regression methods, loess regression does not require the specification of a function to fit a model to all data in the samples. Here, we used loess regression to fit a smooth curve between LAI and the WRE-related variables to better understand the dynamic response of LAI to WRE/ET/SR.  (2) Determining factors that affect WRE based on random forest (RF) algorithm The RF model is an ensemble learning method for regression based on constructing a multitude of decision trees to infer relationships between explanatory and response variables [85]. It has been widely used to address the influential predictors in large predictor sets in ecology-related studies [86][87][88]. Here, RF was used to identify the relative importance of different environmental factors in determining WRE. WRE obtained from each pixel was used as the target variable, while the maximum LAI, clay content (%), litter content (g/m 2 ), annual total precipitation (mm) and slope calculated from DEM (°) were obtained from each pixel and were used as input variables. All these input variables were chosen by referring to prior studies that demonstrate which factors connect tightly with WRE [89,90]. Next, the 'feature importance' produced by RF was used to quantify the relative importance of different factors to WRE.
(3) Using loess regression to fit a smooth curve between LAI and WRE/ET/SR Loess regression is a regression method that creates a smooth line through a scatter plot, displaying the relationship between variables and helping readers to foresee trends [91]. Compared with other regression methods, loess regression does not require the specification of a function to fit a model to all data in the samples. Here, we used loess regression to fit a smooth curve between LAI and the WRE-related variables to better understand the dynamic response of LAI to WRE/ET/SR.    Table S3). Likewise, increased urban areas w mostly converted from cropland (990.35 km 2 ). Meanwhile, cropland expansion and banization were the main causes inducing the loss of wetlands (Figure 4b, Table S3). Furthermore, the density distribution of the LAI was varied in LULCs with differ LULC change trajectories. For forests, the LAI of unchanged forests was mainly distr uted in the high-level range (4-6) ( Figure 5). In contrast, the LAI of new forests conver from other LULCs mainly occurred in the low-level range (0-3) ( Figure 5). Meanwhile, LAI in grassland and cropland with different LULC change trajectories mainly occurr in the low-level range (0-3) ( Figure 5).   Furthermore, the density distribution of the LAI was varied in LULCs with different LULC change trajectories. For forests, the LAI of unchanged forests was mainly distributed in the high-level range (4-6) ( Figure 5). In contrast, the LAI of new forests converted from other LULCs mainly occurred in the low-level range (0-3) ( Figure 5). Meanwhile, the LAI in grassland and cropland with different LULC change trajectories mainly occurred in the low-level range (0-3) ( Figure 5). (Figure 4b), respectively, while wetlands and croplands declined by 34.19% (161.89 km 2 ) and 27.02% (1301.66 km 2 ), respectively (Figure 4b). Dramatic transfer in and transfer out occurred in different LULCs. Increased forests were mostly converted from cropland (364.76 km 2 ), while forest loss was mostly owing to cropland expansion (52.83 km 2 ) and urbanization (65.95 km 2 ) (Figure 4b, Table S3). Likewise, increased urban areas were mostly converted from cropland (990.35 km 2 ). Meanwhile, cropland expansion and urbanization were the main causes inducing the loss of wetlands (Figure 4b, Table S3). Furthermore, the density distribution of the LAI was varied in LULCs with different LULC change trajectories. For forests, the LAI of unchanged forests was mainly distributed in the high-level range (4-6) ( Figure 5). In contrast, the LAI of new forests converted from other LULCs mainly occurred in the low-level range (0-3) ( Figure 5). Meanwhile, the LAI in grassland and cropland with different LULC change trajectories mainly occurred in the low-level range (0-3) ( Figure 5).

Effect of Varied LAI on WRE across Different LULCs
To understand the importance of different environmental factors in determin WRE, these factors were sorted based on RF. Under the condition of considering the motely sensed LAI in Biome-BGC, Table 3 shows that the LAI was the most import factor (RF score was 2.12). Clay content (%) (RF score was 0.97) and litter ranked as second (RF score was 0.65) and third most important factors, respectively. In addition, we compared the response of ET, SR and WRE with LAI dynamics different LULCs. CE increased with increasing LAI but with varied slopes in all LUL (Figure 7). Transpiration also increased with the rising LAI except in forests, in wh transpiration gradually became saturated when the LAI was over 4 (Figure 7). Converse SE decreased with increasing LAI with varied slopes. LE decreased with increasing L but gradually became saturated (Figure 7). SR also decreased with increasing LAI w varied slopes in different LULCs. Specifically, SR decreased with an LAI of less than 4 was depleted after that in forests. The response of WRE to the LAI varied across differ LULCs. In forests, the WRE increased with increasing LAI until 4 but decreased after LAI exceeded 4. However, WRE increased slowly with increasing LAI at various rates other LULCs (Figure 7).

Effect of Varied LAI on WRE across Different LULCs
To understand the importance of different environmental factors in determining WRE, these factors were sorted based on RF. Under the condition of considering the remotely sensed LAI in Biome-BGC, Table 3 shows that the LAI was the most important factor (RF score was 2.12). Clay content (%) (RF score was 0.97) and litter ranked as the second (RF score was 0.65) and third most important factors, respectively. LAI: annual maximum leaf area index; CLAY: clay content (%); LITTER; litter content (g/m 2 ); PRCP: annual total precipitation (mm); SLOPE: slope calculated from DEM ( • ).
In addition, we compared the response of ET, SR and WRE with LAI dynamics in different LULCs. CE increased with increasing LAI but with varied slopes in all LULCs (Figure 7). Transpiration also increased with the rising LAI except in forests, in which transpiration gradually became saturated when the LAI was over 4 (Figure 7). Conversely, SE decreased with increasing LAI with varied slopes. LE decreased with increasing LAI but gradually became saturated (Figure 7). SR also decreased with increasing LAI with varied slopes in different LULCs. Specifically, SR decreased with an LAI of less than 4 but was depleted after that in forests. The response of WRE to the LAI varied across different LULCs. In forests, the WRE increased with increasing LAI until 4 but decreased after the LAI exceeded 4. However, WRE increased slowly with increasing LAI at various rates in other LULCs (Figure 7).

Contribution of Forests to WRE
The results show that forests are important in guaranteeing regional water resource security; they use less than 30% of the total area but contribute over half of the total WRE (Table S3, Figure S5). This was mostly related to a higher WRE capacity in forests than other LULCs. Compared with forests, the simpler species structure and composition in grassland and cropland caused a lower LAI ( Figure 5). The varied LAI among different LULCs is probably the main cause explaining the divergence of WRE capacity (Table 3, Figure 5). Meanwhile, the increased volume of WRE contributed by unchanged forests (14.46 million m 3 ) was much greater than that of new-planted forests (12.28 million m 3 ), but the increase in WRE capacity per unit area in new-planted forests (124.69 ± 14.30 m 3 /ha) was much greater than that of unchanged forests (15.60 ± 7.85 m 3 /ha). Unchanged forests hold higher LAI and grow faster to cause more ET, which was probably the main

Contribution of Forests to WRE
The results show that forests are important in guaranteeing regional water resource security; they use less than 30% of the total area but contribute over half of the total WRE (Table S3, Figure S5). This was mostly related to a higher WRE capacity in forests than other LULCs. Compared with forests, the simpler species structure and composition in grassland and cropland caused a lower LAI ( Figure 5). The varied LAI among different LULCs is probably the main cause explaining the divergence of WRE capacity (Table 3, Figure 5). Meanwhile, the increased volume of WRE contributed by unchanged forests (14.46 million m 3 ) was much greater than that of new-planted forests (12.28 million m 3 ), but the increase in WRE capacity per unit area in new-planted forests (124.69 ± 14.30 m 3 /ha) was much greater than that of unchanged forests (15.60 ± 7.85 m 3 /ha). Unchanged forests hold higher LAI and grow faster to cause more ET, which was probably the main factor to explain the lesser WRE capacity improvement compared with new-planted forests ( Figures 7, S4 and S8). Thus, afforestation would probably quickly increase the volume of WRE in an area.
Although the WRE capacity in most new-planted forests was lower than that of unchanged forests (Figure 7), WRE capacity continued to improve with vegetation growth. New-planted trees occur mostly due to restoration efforts, which are usually of a single species. These artificially planted trees are dominated by poplar and pine and are far from mature [92][93][94]. Young stand age, homogenized species and sparse understory vegetation commonly result in a low LAI [88]. New-planted forests will continue to grow with time, probably for decades. The continued growth of forests inducing increased LAI leads to increased WRE. Nonetheless, the inverted U-shape relationship between LAI and WRE in forests implies that growth did not always promote the increase in WRE (Figures 7 and S8). Although complex canopy structure is more efficient to intercept rainwater to reduce SR, it becomes saturated after forest growth enters a certain stage (Figure 7). ET continues to increase with rising LAI, which in turn causes a decrease in WRE (Figure 7). Prior studies have already demonstrated that the increase in LAI causes more ET after restoration, which leads to a decline in soil water [95][96][97]. Furthermore, the soil water consumption was gradually aggravated along with the increasing forest age, in particular, after 30-35 years of plantation [98][99][100]. Excessive consumption of soil water was the most important reason for the wilt or mortality of forests [101].
Debate exists as to whether forests increase or reduce water resources over a long period [102]. Although previous studies have consistently concluded that increasing forest coverage would decrease SR to reduce total runoff [103][104][105], contradictions have mainly focused on the relationship between forests and baseflow. Some studies proposed that forest decreases WRE, which in turn reduces baseflow [106], while others concluded that forest would increase WRE to promote baseflow [107][108][109]. Indeed, the relationship between forests and WRE is the result of the composite effect of climate, terrain, vegetation and soil [110][111][112]. The relationship between forests and WRE will change as the forest age increases [2,113]. Moreover, the impact of forest age on WRE will be aggravated in rainwater-limited areas such as Beijing [13,114]. Thus, to reduce the possible risk of declining WRE, slow-growth tree species with higher water use efficiency should be given priority over exotic species during restoration [115,116].

Remotely Sensed LAI Role in Biome-BGC
Integrating remotely sensed LAI into Biome-BGC is an efficient way to describe the impact of LULC change with varied vegetation structures on WRE. LAI estimated from Biome-BGC was controlled by climate, terrain, soil and vegetation types but did not consider the impact of LULC change trajectories [55,117]. Thus, it was difficult to describe the spatial information of the vegetation structure owing to varied LULC change trajectories, causing LAI estimated from Biome-BGC to be distant from reality ( Figures S6 and S7). In contrast, the remotely sensed LAI can capture the long-term dynamics in vegetation structure, and especially in capturing the difference in vegetation structure of unchanged old forests and increased forests ( Figure S6). The fact that the importance score of remotely sensed LAI was higher than modelled LAI in determining WRE implied that integrating remotely sensed LAI into Biome-BGC is helpful to describe the constraint of vegetation on hydrological processes (Table 3, Figures S6 and S7). Prior studies also demonstrated that remotely sensed LAI was efficient to describe canopy structure information, which has been used in hydrological variable inversion [48,118,119]. Remotely sensed LAI provides longtime-series, multi-temporal, high-resolution and stable images. Although field surveys and experiments probably obtain more direct and confident LAI, shortages of time and funding restrict the application of field surveys to large areas. Conversely, indirectly measured remotely sensed LAI offers a comprehensive understanding of the landscape at a regional scale. Additionally, the addition of a litter water module could help perfect the water cycle process. Litter increases water interception and infiltration [71,120], and it also reduces surface runoff and soil erosion, especially in old-growth forests [71].
However, the limitation of remotely sensed LAI is that it is difficult to capture accurate information on understory vegetation [121,122]. The understory vegetation would induce higher ET and intercept more precipitation to reduce the SR. Thus, our results probably underestimate the ET and overestimate the SR. It is still unclear whether these opposite effects on WRE would be counteracted totally. In most cases, the size and biomass of understory vegetation are much lower than those of canopy vegetation [123][124][125]. The effect of understory vegetation on WRE might be limited [126][127][128]. Some studies have implied that Lidar has a stronger canopy penetration ability and the potential to detect information on understory vegetation [129]. However, high acquisition costs and low availability restrict its application, especially in developing countries [130]. In future, we will attempt to build the relationship between canopy structure and understory vegetation [131] based on measured data. Meanwhile, we will modify EPM to include processes that consider the understory vegetation structure and its impact on hydrological processes.

Conclusions
Previous studies that focused on the impact of LULC change on WRE neglected the effect of varied vegetation structures. This neglect has hampered a more comprehensive acknowledgment of the impact of LULC changes on WRE. In this study, we integrated remotely sensed LAI into the ecosystem process-based Biome-BGC to analyse the impact of LULC changes on WRE in Beijing between 2000 and 2015, including consideration of the varied vegetation structures. We found that the volume of WRE in 2015 increased by approximately 8.58 million m 3 compared with 2000. Unchanged forests and new-planted forests contributed approximately 14.46 million m 3 and 12.28 million m 3 , respectively. In contrast, urbanization, cropland expansion and deforestation caused WRE to decline by 11.96 million m 3 , 5.86 million m 3 and 3.20 million m 3 , respectively. Although the increased WRE contributed by unchanged forests (14.46 million m 3 ) was much more than that of new-planted forests (12.28 million m 3 ), the increase in WRE capacity per unit area in newplanted forests (124.69 ± 14.30 m 3 /ha) was much greater than that of unchanged forests (15.60 ± 7.85 m 3 /ha). Higher LAI in unchanged forests induced more ET, which was the main cause leading to the lower WRE capacity per unit area. Our results suggest that the effect of vegetation structure on WRE cannot be ignored. Integrating the remotely sensed LAI with Biome-BGC was efficient in analysing the effect of varied vegetation structures when quantifying the impact of LULC change on WRE. Furthermore, our approach can be used to better understand the comprehensive results of the impact of LULC change on WRE and guide land use policies to maintain or increase regionally available water resources.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

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