Modelling the Impact of Vegetation Change on Hydrological Processes in Bayin River Basin, Northwest China

: Vegetation change in arid areas may lead to the redistribution of regional water resources, which can intensify the competition between ecosystems and humans for water resources. This study aimed to accurately model the impact of vegetation change on hydrological processes in an arid endorheic river watershed undergoing revegetation, namely, the middle and lower reaches of the Bayin River basin, China. A LU-SWAT-MODFLOW model was developed by integrating dynamic hydrological response units with a coupled SWAT-MODFLOW model, which can reﬂect actual land cover changes in the basin. The LU-SWAT-MODFLOW model outperformed the original SWAT-MODFLOW model in simulating the impact of human activity as well as the leaf area index, evapotranspiration, and groundwater table depth. After regional revegetation, evapotranspiration and groundwater recharge in different sub-basins increased signiﬁcantly. In addition, the direction and amount of surface-water–groundwater exchange changed considerably in areas where revegetation involved converting low-coverage grassland and bare land to forestland.


Introduction
Vegetation is essential for regional carbon sequestration, soil and water conservation, and climate regulation [1,2]. Arid areas, which account for 40% of the world's land area, are characterised by water shortages and uneven spatiotemporal distributions of water resources [3]. Changes in vegetation and related management practices (e.g., irrigation) in arid areas may lead to the redistribution of regional water resources, which can intensify the competition between ecosystems and humans for water resources [1,4]. In this context, the water demand and water consumption characteristics of vegetation change in arid areas are of particular concern [5,6].
Nowadays, physically based distributed (or semi-distributed) hydrological models can clearly reflect the spatial variability of hydrological processes in a basin, and these models are playing an important role in simulations and predictions of the hydrological cycle in basins [7][8][9]. Notably, SWAT (Soil and Water Assessment Tools) is a typical distributed hydrological model with a strong physical foundation [10]. It is suitable for simulating surface hydrological processes in a complex basin with a variety of soil types, land use types, slopes, and management practices, and it can be used in data-poor regions [11][12][13]. Currently, SWAT is a key component of the USDA-Conservation Effect Assessment Project and the USEPA-Hydrologic and Water Quality System [14]. Nevertheless, SWAT has a weak ability to simulate groundwater processes, thereby limiting its application in arid areas with strong surface-water-groundwater exchange [15][16][17].
The ability of SWAT to simulate groundwater processes can be improved by replacing the groundwater module of SWAT with a well-established groundwater model [15,16].
into the middle and lower reaches, where the human population and industrial and agricultural activities are concentrated, and with frequent surface-water-groundwater exchange. The vegetation in the Bayin River basin has been restored substantially with the implementation of a series of ecological restoration measures, such as the Grain-for-Green program, over the past 20 years. However, irrigation has become essential for such artificial revegetation projects because of the arid climate and the heterogeneity in the spatial distribution of water resources.

SWAT-MODFLOW Model
In this study, the coupling of SWAT with the MODFLOW model was achieved by using the method of Bailey et al. [16]. Given the inconsistency between the two models in terms of computational spatial units, it was necessary to use a GIS platform to unify the spatial resolution before model coupling ( Figure 2). First, specific spatial locations were allocated to the computational units (i.e., HRUs) of SWAT. Second, a mapping relationship was established between the HRUs of SWAT and the computational grid cells of MODFLOW on the same projected coordinate system using a GIS platform [14,15]. Third, the SWAT model was run to simulate the groundwater recharge, evaporation, and extraction with a temporal step of 1 d. Finally, the simulation results were taken as boundary conditions on the corresponding computational grid cells of MODFLOW for groundwater flow modelling [16]. The MODFLOW model was run to simulate the groundwater processes while using the groundwater monitoring data of the basin (provided by Qinghai Provincial Department of water resources) to calibrate and validate the model parameters. Meanwhile, the simulated groundwater table depth from the MODFLOW model was transferred to the computational units of surface water through the abovementioned mapping relationship to impose boundary conditions on the simulation of irrigation groundwater extraction, crop growth, and vegetation transpiration, as well as to test the simulation results [16].

SWAT-MODFLOW Model
In this study, the coupling of SWAT with the MODFLOW model was achieved by using the method of Bailey et al. [16]. Given the inconsistency between the two models in terms of computational spatial units, it was necessary to use a GIS platform to unify the spatial resolution before model coupling ( Figure 2). First, specific spatial locations were allocated to the computational units (i.e., HRUs) of SWAT. Second, a mapping relationship was established between the HRUs of SWAT and the computational grid cells of MOD-FLOW on the same projected coordinate system using a GIS platform [14,15]. Third, the SWAT model was run to simulate the groundwater recharge, evaporation, and extraction with a temporal step of 1 d. Finally, the simulation results were taken as boundary conditions on the corresponding computational grid cells of MODFLOW for groundwater flow modelling [16]. The MODFLOW model was run to simulate the groundwater processes while using the groundwater monitoring data of the basin (provided by Qinghai Provincial Department of water resources) to calibrate and validate the model parameters. Meanwhile, the simulated groundwater table depth from the MODFLOW model was transferred to the computational units of surface water through the abovementioned mapping relationship to impose boundary conditions on the simulation of irrigation groundwater extraction, crop growth, and vegetation transpiration, as well as to test the simulation results [16]. Water 2021, 13, x FOR PEER REVIEW 4 of 20

SWAT-MODFLOW Model Coupled with Dynamic HRUs
In view of the inability of the original SWAT model to effectively reflect complete or partial land cover type conversions within the same HRU, this study transformed the HRUs of the original SWAT model to dynamic HRUs to improve the original model. The generation process of dynamic HRUs is illustrated in Figure 3. In contrast to the original HRUs, the generation process of dynamic HRUs involved the defining of spatial units where there were land use/cover changes, i.e., it incorporated the concept of dynamic land use/cover. Such spatial units were combined with soil type and slope data to generate dynamic HRUs such that each had a specific and invariant location, area, and shape with variable attributes. Such dynamic HRUs can more truly reflect the land use/cover changes in the basin. The operational flow chart of the coupled SWAT-MODFLOW model based on dynamic HRUs is shown in Figure 4. First, annual cycle simulations were performed by using corresponding HRUs (land cover). Second, daily cycle simulations were performed in which hydrological processes were simulated by SWAT. The simulation results on each HRU were mapped to the computational grid cells of MODFLOW, where these were taken as boundary conditions for groundwater flow simulations. Third, the simulated groundwater data within the grid cells were mapped to the HRUs of SWAT for subsequent SWAT computations. The above process was conducted in nested loops until the

SWAT-MODFLOW Model Coupled with Dynamic HRUs
In view of the inability of the original SWAT model to effectively reflect complete or partial land cover type conversions within the same HRU, this study transformed the HRUs of the original SWAT model to dynamic HRUs to improve the original model. The generation process of dynamic HRUs is illustrated in Figure 3. In contrast to the original HRUs, the generation process of dynamic HRUs involved the defining of spatial units where there were land use/cover changes, i.e., it incorporated the concept of dynamic land use/cover. Such spatial units were combined with soil type and slope data to generate dynamic HRUs such that each had a specific and invariant location, area, and shape with variable attributes. Such dynamic HRUs can more truly reflect the land use/cover changes in the basin.

SWAT-MODFLOW Model Coupled with Dynamic HRUs
In view of the inability of the original SWAT model to effectively reflect complete or partial land cover type conversions within the same HRU, this study transformed the HRUs of the original SWAT model to dynamic HRUs to improve the original model. The generation process of dynamic HRUs is illustrated in Figure 3. In contrast to the original HRUs, the generation process of dynamic HRUs involved the defining of spatial units where there were land use/cover changes, i.e., it incorporated the concept of dynamic land use/cover. Such spatial units were combined with soil type and slope data to generate dynamic HRUs such that each had a specific and invariant location, area, and shape with variable attributes. Such dynamic HRUs can more truly reflect the land use/cover changes in the basin. The operational flow chart of the coupled SWAT-MODFLOW model based on dynamic HRUs is shown in Figure 4. First, annual cycle simulations were performed by using corresponding HRUs (land cover). Second, daily cycle simulations were performed in which hydrological processes were simulated by SWAT. The simulation results on each HRU were mapped to the computational grid cells of MODFLOW, where these were taken as boundary conditions for groundwater flow simulations. Third, the simulated groundwater data within the grid cells were mapped to the HRUs of SWAT for subsequent SWAT computations. The above process was conducted in nested loops until the The operational flow chart of the coupled SWAT-MODFLOW model based on dynamic HRUs is shown in Figure 4. First, annual cycle simulations were performed by using corresponding HRUs (land cover). Second, daily cycle simulations were performed in which hydrological processes were simulated by SWAT. The simulation results on each HRU were mapped to the computational grid cells of MODFLOW, where these were taken as boundary conditions for groundwater flow simulations. Third, the simulated groundwater data within the grid cells were mapped to the HRUs of SWAT for subsequent SWAT computations. The above process was conducted in nested loops until the end of the simulation. Considering that the dynamic HRU-based SWAT-MODFLOW coupled model could reflect the dynamic changes in land use/cover, the model was referred to as LU-SWAT-MODFLOW.

Model Establishment
Meteorological data required for establishing the SWAT model included the daily monitoring data for precipitation, temperature, relative humidity, wind speed, and solar radiation at the Delingha weather station, which is located at the inlet shown in Figure 1. Agricultural, forestland, and grassland irrigation data were obtained from the Delingha Municipal Water Affairs Bureau ( Table 1). The newest digital elevation model data (Figure 1) of 30 m resolution Shuttle Radar Topography Mission data [28] were used. Soil type data ( Figure 5a) at the scale 1:1,000,000 from China were used [29], and the corresponding soil hydrological attributes were retrieved from the Qinghai Soil Record [30]. Figure 5b presents the sub-basin division map of the study region.

Number of Times of Irrigation Irrigation Duration
Agricultural irrigation 5800 6 March-October Forest irrigation 5400 6 April-November Grassland irrigation 3600 5 April-November

Model Establishment
Meteorological data required for establishing the SWAT model included the daily monitoring data for precipitation, temperature, relative humidity, wind speed, and solar radiation at the Delingha weather station, which is located at the inlet shown in Figure 1. Agricultural, forestland, and grassland irrigation data were obtained from the Delingha Municipal Water Affairs Bureau ( Table 1). The newest digital elevation model data ( Figure 1) of 30 m resolution Shuttle Radar Topography Mission data [28] were used. Soil type data ( Figure 5a) at the scale 1:1,000,000 from China were used [29], and the corresponding soil hydrological attributes were retrieved from the Qinghai Soil Record [30]. Figure 5b presents the sub-basin division map of the study region. Land use type data for the simulated period (2000-2018) were derived from 30 m Landsat images. Remote sensing interpretation marks were created according to spectral features combined with field survey data and relevant geographic maps [31]. Data quality was examined by comparatively analysing field survey patches versus randomly selected patches, and the classification accuracy was determined to be over 90%. Figure 6 presents the regional spatial distribution of land use/cover in 2000 ( Figure 6a) versus 2018 ( Figure 6b). There were six types of land use/cover in 2000, namely, spring wheat, forest, grassland, water, residences, and barren land. There were seven types in 2018, including 'Chinese wolfberry' as a new type in addition to the existing six types. In the study region, Chinese wolfberry was the main tree species used for artificial revegetation. Considering the absence of Chinese wolfberry-related parameters in the built-in land use and vegetation database of SWAT, relevant initial parameters for apple trees in the SWAT model were taken as default parameters for Chinese wolfberry, and these were calibrated using the LAI data to simulate the growth process of Chinese wolfberry. Other relevant parameters of land use/cover types were either set to the default values in the built-in database of the model or obtained by calibration. Revegetation in the study region was mainly characterised by the conversion of farmland to forestland and the conversion of bare land to forestland and grassland.  Land use type data for the simulated period (2000-2018) were derived from 30 m Landsat images. Remote sensing interpretation marks were created according to spectral features combined with field survey data and relevant geographic maps [31]. Data quality was examined by comparatively analysing field survey patches versus randomly selected patches, and the classification accuracy was determined to be over 90%. Figure 6 presents the regional spatial distribution of land use/cover in 2000 ( Figure 6a) versus 2018 ( Figure  6b). There were six types of land use/cover in 2000, namely, spring wheat, forest, grassland, water, residences, and barren land. There were seven types in 2018, including 'Chinese wolfberry' as a new type in addition to the existing six types. In the study region, Chinese wolfberry was the main tree species used for artificial revegetation. Considering the absence of Chinese wolfberry-related parameters in the built-in land use and vegetation database of SWAT, relevant initial parameters for apple trees in the SWAT model were taken as default parameters for Chinese wolfberry, and these were calibrated using the LAI data to simulate the growth process of Chinese wolfberry. Other relevant parameters of land use/cover types were either set to the default values in the built-in database of the model or obtained by calibration. Revegetation in the study region was mainly characterised by the conversion of farmland to forestland and the conversion of bare land to forestland and grassland.  Basin boundaries delineated by the SWAT model were considered as impermeable boundaries to groundwater flow in the MODFLOW model, where the western and eastern outlets of the basin were considered as constant flow boundaries. The river network extracted by the SWAT model was considered to constitute river boundaries in the MOD-FLOW model. The simulated steady-state groundwater head ( Figure 7) was used as the initial head for simulations of transient flow [32]. In addition, the shallow aquifers in the study region were conceptualised as being non-homogeneous and anisotropic according to relevant studies [32], and the groundwater flow was conceptualised as a two-dimensional transient flow. Basin boundaries delineated by the SWAT model were considered as impermeable boundaries to groundwater flow in the MODFLOW model, where the western and eastern outlets of the basin were considered as constant flow boundaries. The river network extracted by the SWAT model was considered to constitute river boundaries in the MOD-FLOW model. The simulated steady-state groundwater head (Figure 7) was used as the initial head for simulations of transient flow [32]. In addition, the shallow aquifers in the study region were conceptualised as being non-homogeneous and anisotropic according to relevant studies [32], and the groundwater flow was conceptualised as a two-dimensional transient flow.
Basin boundaries delineated by the SWAT model were considered as impermeable boundaries to groundwater flow in the MODFLOW model, where the western and eastern outlets of the basin were considered as constant flow boundaries. The river network extracted by the SWAT model was considered to constitute river boundaries in the MOD-FLOW model. The simulated steady-state groundwater head (Figure 7) was used as the initial head for simulations of transient flow [32]. In addition, the shallow aquifers in the study region were conceptualised as being non-homogeneous and anisotropic according to relevant studies [32], and the groundwater flow was conceptualised as a two-dimensional transient flow.

Model Calibration
The  . The dataset is derived by employing a calibration-free nonlinear complementary relationship model with inputs of air and dew-point temperature, wind speed, precipitation, and net radiation from the China Meteorological Forcing Dataset. The dataset is validated in Northwest China and is proved to have good spatial and temporal performance [24]. Furthermore, relevant parameters (conductivity and storage) of the MODFLOW model were calibrated against monthly recorded groundwater table depth at numerous observation wells (Figure 5b) in the basin.
In this study, land cover types in the LU-SWAT-MODFLOW model varied over the years within some of the dynamic HRUs or remained invariant throughout a relatively long period within some of the HRUs. The HRUs in the latter scenario were chosen to calibrate the relevant vegetation growth parameters ( Table 2) against monthly 30 m resolution LAI data [26]. The calibrated parameters were stored in a separate file, which could be visited during SWAT operations to directly tune parameters in HRUs where changes in land cover type were detected. ET-related parameter ( Table 2) calibration at sub-basin scales in the present study was mainly based on the aforementioned remote sensing-derived ET data in accordance with the method of White and Chaubey [33] and Immerzeel and Droogers [25].
calibration-free nonlinear complementary relationship model with inputs of air and dewpoint temperature, wind speed, precipitation, and net radiation from the China Meteorological Forcing Dataset. The dataset is validated in Northwest China and is proved to have good spatial and temporal performance [24]. Furthermore, relevant parameters (conductivity and storage) of the MODFLOW model were calibrated against monthly recorded groundwater table depth at numerous observation wells (Figure 5b) in the basin. In this study, land cover types in the LU-SWAT-MODFLOW model varied over the years within some of the dynamic HRUs or remained invariant throughout a relatively long period within some of the HRUs. The HRUs in the latter scenario were chosen to calibrate the relevant vegetation growth parameters (Table 2) against monthly 30 m resolution LAI data [26]. The calibrated parameters were stored in a separate file, which could be visited during SWAT operations to directly tune parameters in HRUs where changes in land cover type were detected. ET-related parameter ( Table 2) calibration at sub-basin scales in the present study was mainly based on the aforementioned remote sensing-derived ET data in accordance with the method of White and Chaubey [33] and Immerzeel and Droogers [25].

Processes
Parameter Description BLAI Maximum potential leaf area index

Model Performance Metrics
The applicability of SWAT was evaluated in terms of Nash-Sutcliffe Efficiency (NSE), percent bias (PBIAS), and squared correlation coefficient (R 2 ). The NSE can range between −∞ and 1. If the value of the NSE was closer to 1, it was indicative of a better simulation performance and high reliability of the SWAT model. When the NSE was closer to 0.5, the model simulation results were similar to the mean of observations, that is, the model results in general were reliable. A PBIAS between −10% and 10% indicated a good simulation performance of the model. Additionally, larger values of R 2 were indicative of a better simulation performance of the model. The calculation process and significance of the three metrics have been elaborated elsewhere [34]. Model performance during the simulations of the groundwater table depth was evaluated mainly in terms of the absolute error and R 2 in this study.

Regional Vegetation Change
The main types of vegetation change in the study region from 2002 to 2018 are summarised in Table 3. Specifically, the main types of vegetation change pertained to the conversion of low-coverage grassland, farmland, and bare land to forestland, in which the converted areas amounted to 2,877,518 and 321 hm 2 , respectively. Correspondingly, the irrigation water volume of each revegetation plot underwent dramatic changes. In addition, the conversion of farmland to forestland mainly occurred in the irrigated area of the Gahai Lake in the south-eastern part of the basin, while the conversion of low-coverage grassland and bare land to forestland mainly occurred in the irrigated area of Delingha, which is situated in the north-western part of the basin ( Figure 6). Table 3. Main types of revegetation in the study region.

Change in the Annual Irrigation Rate (m 3 /hm 2 ·a)
Conversion of low-coverage grassland to forestland 2877 0→5400 Conversion of farmland to forestland 518 5800→5400 Conversion of bare land to forestland 321 0→5400 Figure 9 shows the HRUs generated by the original SWAT-MODFLOW model versus those from the LU-SWAT-MODFLOW model. The original SWAT-MODFLOW model generated 1304 HRUs, and the LU-SWAT-MODFLOW model generated 2978 HRUs. The higher number of HRUs in the new model was due to the land use/cover data of LU-SWAT-MODFLOW being a superposition of years-long data and thereby covering a higher number of patches. The higher number of HRUs also implied that the operation and parameter tuning of the LU-SWAT-MODFLOW model would be more complicated [35]. those from the LU-SWAT-MODFLOW model. The original SWAT-MODFLOW m generated 1304 HRUs, and the LU-SWAT-MODFLOW model generated 2978 HRUs higher number of HRUs in the new model was due to the land use/cover data of SWAT-MODFLOW being a superposition of years-long data and thereby coveri higher number of patches. The higher number of HRUs also implied that the opera and parameter tuning of the LU-SWAT-MODFLOW model would be more complic [35].    (Table 4). This indicates that the LU-SWAT-MODFLOW model was more accurate than the original SWAT-MODFLOW model in simulations of the monthly LAI after both models were calibrated and validated.  (Table 4). This indicates that the LU-SWAT-MODFLOW model was more accurate than the original SWAT-MODFLOW model in simulations of the monthly LAI after both models were calibrated and validated.

Comparison of the ET Simulation Results
During the calibration and validation period, the performance metrics of the original SWAT-MODFLOW model were NSE > 0.65, PBIAS of −20%-20%, and R 2 > 0.63 in simulations of the monthly mean ET for each sub-basin, while the counterparts for the LU-SWAT-MODFLOW model were NSE > 0.72, PBIAS of −20%-20%, and R 2 > 0.73 ( Figure 11). This indicates that the LU-SWAT-MODFLOW model was more accurate than the original SWAT-MODFLOW in simulations of the monthly mean ET for most of the sub-basins. Figure 12 shows the multi-year mean of the simulated ET from the original SWAT-MODFLOW model (Figure 12a) versus that of the LU-SWAT-MODFLOW model (Figure 12b) in comparison with the multi-year mean of remote sensing-derived ET (Figure 12c). The multi-year mean of remote sensing-derived ET exhibited a spatial distribution pattern of high values in the north-eastern mountains and low values in the southwestern plains, and such a distribution pattern existed for both calibrated models. Table Depth Groundwater table depth data were scarce within the study region. Observation wells 1, 2, and 3 only provided monthly data for 2009-2011, and observation well 4 only provided monthly data for 2013-2015; meanwhile, observation well 5 only provided monthly data for 2014-2015. Thus, the observed groundwater table depth of wells 1, 2, and 3 were used to calibrate the two models, and the rest were used to validate the models. Linear regression results of the simulated groundwater table depth on observed groundwater table depth were compared between the original SWAT-MODFLOW model and the LU-SWAT-MODFLOW model ( Figure 13). Both models performed well in simulating the changes in the groundwater table depth of the study region, with an R 2 > 0.95 and absolute error within 0.5 m. In addition, the simulation performance of LU-SWAT-MODFLOW was slightly better than that of SWAT-MODFLOW, which was likely attributed to the detailed consideration of the spatiotemporal changes in irrigation and land cover by the former model versus the latter model. original SWAT-MODFLOW in simulations of the monthly mean ET for most of the subbasins. Figure 12 shows the multi-year mean of the simulated ET from the original SWAT-MODFLOW model (Figure 12a) versus that of the LU-SWAT-MODFLOW model ( Figure  12b) in comparison with the multi-year mean of remote sensing-derived ET (Figure 12c). The multi-year mean of remote sensing-derived ET exhibited a spatial distribution pattern of high values in the north-eastern mountains and low values in the southwestern plains, and such a distribution pattern existed for both calibrated models.   Table Depth Groundwater table depth data were scarce within the study region. Observation

Impacts of Vegetation Change on Hydrological Processes
The case study area is located in a water consumption area of an inland river and almost never generates runoff. Thus, here we focus on the analysis of the vegetation change impacts on ET and groundwater processes.

Impacts on ET
The LU-SWAT-MODFLOW model was run in the following two scenarios to accurately analyse the impacts of revegetation and the related extensive irrigation on ET: (1) revegetation was assumed absent while considering the actual changes in other types of land use/cover; and (2) the actual changes in land use/cover were considered, including those pertinent to revegetation (irrigation) and other types of land use/cover. Figure 14a

Impacts of Vegetation Change on Hydrological Processes
The case study area is located in a water consumption area of an inland river and almost never generates runoff. Thus, here we focus on the analysis of the vegetation change impacts on ET and groundwater processes.

Impacts on ET
The LU-SWAT-MODFLOW model was run in the following two scenarios to accurately analyse the impacts of revegetation and the related extensive irrigation on ET: (1) revegetation was assumed absent while considering the actual changes in other types of land use/cover; and (2) the actual changes in land use/cover were considered, including those pertinent to revegetation (irrigation) and other types of land use/cover. Figure 14a shows the simulated monthly ET in the revegetation-absent scenario versus the revegetation-present scenario from 2002 to 2018. The results indicate that revegetation and related irrigation did not change the trend of monthly ET in the basin, in which the monthly ET in the revegetation-present scenario was only 1.5 mm higher than that in the revegetation-absent scenario for most months. Similarly, the trend of annual ET was almost the same in both scenarios. In 2004 and later years, ET showed weakly higher values in the revegetation-present scenario than in the revegetation-absent scenario (Figure 14b). Figure 13c illustrates the difference in the multi-year mean ET between the two scenarios in each sub-basin. Such a difference was greater than 10 mm in sub-basins 4, 12, 13, 14, 26, and 33, that is, the ET increase was most obvious in these sub-basins. Comprehensive comparisons of the land use/cover map ( Figure 6) with the LAI map for the study region during the study period further confirmed that relatively obvious revegetation had been achieved in these sub-basins.
Water 2021, 13, x FOR PEER REVIEW 15 o shows the simulated monthly ET in the revegetation-absent scenario versus the revege tion-present scenario from 2002 to 2018. The results indicate that revegetation and rela irrigation did not change the trend of monthly ET in the basin, in which the monthly in the revegetation-present scenario was only 1.5 mm higher than that in the revegetati absent scenario for most months. Similarly, the trend of annual ET was almost the sa in both scenarios. In 2004 and later years, ET showed weakly higher values in the reve tation-present scenario than in the revegetation-absent scenario (Figure 14b). Figure  illustrates the difference in the multi-year mean ET between the two scenarios in each s basin. Such a difference was greater than 10 mm in sub-basins 4, 12, 13, 14, 26, and 33, t is, the ET increase was most obvious in these sub-basins. Comprehensive comparisons the land use/cover map ( Figure 6) with the LAI map for the study region during the stu period further confirmed that relatively obvious revegetation had been achieved in th sub-basins.

Impacts on Groundwater Recharge
Groundwater is the most important water resource in the arid endorheic river wat shed. Changes in groundwater recharge may affect the groundwater storage and furt impact the ecological environment. Figure 15 shows the monthly (Figure 15a) and yea (Figure 15b) groundwater recharge in the entire study area. After revegetation, groundwater recharge increased by approximately 1.27 mm on average per month a 14.02 mm on average per year. Fan et al. [36], Yang and Lu [37], and Qubaja et al. [ showed that canopy interception and root water absorption would lead to reduction soil water, surface runoff, and groundwater recharge in woodland. However, here, hough considerable areas of low-coverage grassland, farmland, and bare land were c verted to forestland, the groundwater recharge with revegetation was evidently hig than that without revegetation. We reported the yearly average groundwater recha after revegetation in the entire study area (Figure 14c). The groundwater recharge in irrigation district where the revegetation was applied was the highest (>14.51 m 3 /da that is, the irrigation for the recovered vegetation strongly affected the groundwater charge.

Impacts on Groundwater Recharge
Groundwater is the most important water resource in the arid endorheic river watershed. Changes in groundwater recharge may affect the groundwater storage and further impact the ecological environment. Figure 15 shows the monthly (Figure 15a) and yearly (Figure 15b) groundwater recharge in the entire study area. After revegetation, the groundwater recharge increased by approximately 1.27 mm on average per month and 14.02 mm on average per year. Fan et al. [36], Yang and Lu [37], and Qubaja et al. [30] showed that canopy interception and root water absorption would lead to reduction of soil water, surface runoff, and groundwater recharge in woodland. However, here, although considerable areas of low-coverage grassland, farmland, and bare land were converted to forestland, the groundwater recharge with revegetation was evidently higher than that without revegetation. We reported the yearly average groundwater recharge after revegetation in the entire study area (Figure 14c). The groundwater recharge in the irrigation district where the revegetation was applied was the highest (>14.51 m 3 /day); that is, the irrigation for the recovered vegetation strongly affected the groundwater recharge.

Impacts on Surface Water and Groundwater Exchange
There was frequent surface-water-groundwater exchange in the study region, which dominated the regional hydrological processes. We analysed the surface water and groundwater exchange affected by revegetation. Figure 16 shows the amount of groundwater recharge and discharge in the revegetation-absent scenario minus the value in the revegetation-present scenario. Specifically, river reach I was in the upper study region, where groundwater was recharged by river water. River reach II was situated in the lower study region, where significant amounts of groundwater were discharged to the river. In river reach III, both surface water recharge to groundwater and groundwater discharge to surface water were present. River reach III was situated in the irrigated area of Delingha, where it was greatly affected by agricultural, forestland, and grassland irrigation, which led to a relatively complex pattern of surface-water-groundwater exchange. The area where river reach III was situated was also the main revegetation area of the study region. Comparisons of Figure 16 revealed that the direction of surface-water-groundwater exchange was reversed in six grid cells, which was attributed to the changes in the irrigation volume within these grid cells after revegetation.

Impacts on Surface Water and Groundwater Exchange
There was frequent surface-water-groundwater exchange in the study region, which dominated the regional hydrological processes. We analysed the surface water and groundwater exchange affected by revegetation. Figure 16 shows the amount of groundwater recharge and discharge in the revegetation-absent scenario minus the value in the revegetation-present scenario. Specifically, river reach I was in the upper study region, where groundwater was recharged by river water. River reach II was situated in the lower study region, where significant amounts of groundwater were discharged to the river. In river reach III, both surface water recharge to groundwater and groundwater discharge to surface water were present. River reach III was situated in the irrigated area of Delingha, where it was greatly affected by agricultural, forestland, and grassland irrigation, which led to a relatively complex pattern of surface-water-groundwater exchange. The area where river reach III was situated was also the main revegetation area of the study region. Comparisons of Figure 16 revealed that the direction of surface-water-groundwater exchange was reversed in six grid cells, which was attributed to the changes in the irrigation volume within these grid cells after revegetation.

Discussion
The LU-SWAT-MODFLOW model was more accurate than the original SWAT-MODFLOW model in simulating the monthly LAI after both models were calibrated and validated. LAI plays a key role in SWAT for estimating ET, canopy interception, and

Discussion
The LU-SWAT-MODFLOW model was more accurate than the original SWAT-MODFLOW model in simulating the monthly LAI after both models were calibrated and validated. LAI plays a key role in SWAT for estimating ET, canopy interception, and biomass accumulation [35]. The enhanced modelling of LAI could improve the performance of the SWAT model in eco-hydrological processes [26,38]. However, accurate simulation of LAI relies on many parameters which are difficult to calibrate. Generally, parameters of SWAT-MODFLOW are calibrated with observed data in watershed outlets or sub-basins [23,27]. Only few studies have calibrated the parameters at the HRU level because of its difficulty and complexity. In this study, the remote-sensed monthly LAI data were used to calibrate the SWAT-MODFLOW and LU-SWAT-MODFLOW at the HRU level using the SWAT-CUP software (https://swat.tamu.edu/software/swat-cup/) (accessed on 1 October 2021) with a satisfactory result. This suggests that the model calibration at the HRU level is possible and effective if the related observation data exist.
The LU-SWAT-MODFLOW model was more accurate than the original SWAT-MODFLOW in simulating the monthly mean ET for most sub-basins. Ma et al. [26] reported that canopy interception and soil water content would be seriously affected by LAI in SWAT, which would further affect ET. Therefore, the enhancements of LU-SWAT-MODFLOW in modelling the monthly ET can be attributed to the more accurate simulation of the LAI.
Revegetation projects have been conducted in both the Gahai Lake irrigated area and the irrigated area of Delingha, but the revegetation had a relatively high impact on the direction and amount of surface-water-groundwater exchange in the latter area; in the former area, there was an almost negligible impact. This discrepancy was attributed to the fact that revegetation in the Gahai Lake irrigated area was mainly characterised by the conversion of farmland to forestland, and the irrigation volume did not differ significantly between the two land cover types [38]. In contrast, the irrigated area of Delingha was dominated by the conversion of low-coverage grassland to bare land and forestland, and the former two land cover types required no irrigation, while the latter land cover type required a large irrigation volume.
This study is subjected to some limitations. On the one hand, we used land use/cover map to analyse the revegetation process in our study area. In fact, plant density, age, and growth status were not considered because of the limitations in the SWAT model. Moreover, these factors may affect the eco-hydrological processes in such an arid area [26]. On the other hand, meteorological data were scarce in both original SWAT-MODFLOW and LU-SWAT-MODFLOW models. This may impact the model performance in formulating the water budget [39][40][41]. Nonetheless, these limitations should be addressed in future studies by using and analysing different datasets.

Conclusions
This study was carried out in the middle and lower reaches of the Bayin River basin in the north-eastern part of the Qaidam Basin, China, where there is frequent surface-watergroundwater interaction and evident vegetation change. A LU-SWAT-MODFLOW model was developed by integrating a coupled SWAT-MODFLOW model with dynamic HRUs in view of their ability to reflect the actual land cover changes in the basin. The impacts of revegetation and related irrigation on the main hydrological processes in the basin were more accurately simulated and analysed by the LU-SWAT-MODFLOW model than by the original SWAT-MODFLOW model.
The LU-SWAT-MODFLOW model generated dynamic HRUs by pre-defining spatial units where land use/cover changes occurred during the simulated period, thereby overcoming the inability of the original SWAT model to effectively reflect the complete or partial land cover type conversion within the same HRU. This new model outperformed the original SWAT-MODFLOW model in simulating the LAI. The LAI is an important parameter of SWAT as it affects a series of processes, such as ET and infiltration; therefore, accurate simulations of the LAI are a key to accurate hydrological simulations. Moreover, the LU-SWAT-MODFLOW model outperformed the original SWAT-MODFLOW model in simulating the ET and groundwater table depth of the basin.
The LU-SWAT-MODFLOW model was run in two different scenarios, one with revegetation and the other without it, to assess the impacts of revegetation and related irrigation on the main hydrological processes in the study region. The results showed that after regional revegetation, ET in the different sub-basins increased by approximately 1.5 mm per month and by 6 mm per year. After revegetation, the groundwater recharge increased by approximately 1.27 mm on average per month and 14.02 mm on average per year. Irrigation for the recovered vegetation strongly affected the groundwater recharge. Meanwhile, the direction and amount of surface-water-groundwater exchange underwent evident changes in areas where revegetation was characterised by the conversion of low-coverage grassland and bare land to forestland. In areas where revegetation was characterised by the conversion of farmland to forestland, the irrigation volume was not greatly altered; thus, this transition had a weak impact on the direction and amount of surface-water-groundwater exchange. Changes in the direction and amount of surface-water-groundwater exchange may lead to a series of ecological and environmental issues. To avoid problems in the future, water-saving irrigation techniques should be advocated when conducting revegetation in arid inland river basins. In addition, our findings indicate that it would be advantageous to preferentially apply revegetation measures that promote the conversion of farmland to forestland/grassland provided that they do not adversely affect regional economic development.

Data Availability Statement:
The processed data required to reproduce these findings cannot be shared at this time as the data is also a part of an ongoing study.

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