*2.1. Study Area*

The Bayin River, which is the fourth largest river in the Qaidam Basin, is situated in the north-western region of China (Figure 1). The basin is in an arid area with annual precipitation of approximately 200 mm. The main land cover types of the area are grassland, shrubland, barren land, and farmland. The Bayin River flows out of the mountains 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. 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 artifi‐ cial revegetation projects because of the arid climate and the heterogeneity in the spatial distribution of water resources.

into the middle and lower reaches, where the human population and industrial and agri‐ cultural activities are concentrated, and with frequent surface‐water–groundwater ex‐ change. The vegetation in the Bayin River basin has been restored substantially with the

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

**Figure 1.** Study area. **Figure 1.** Study area.

tion results [16].

#### *2.2. SWAT‐MODFLOW Model 2.2. 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 relation‐ ship 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 extrac‐ tion 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 pro‐ cesses 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 waterthrough the abovementioned map‐ ping relationship to impose boundary conditions on the simulation of irrigation ground‐ water extraction, crop growth, and vegetation transpiration, as well as to test the simula‐ 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].

**Figure 2.** Integration of Soil and Water Assessment Tool (SWAT) and MODFLOW computing units. **Figure 2.** Integration of Soil and Water Assessment Tool (SWAT) and MODFLOW computing units.

#### *2.3. SWAT‐MODFLOW Model Coupled with Dynamic HRUs 2.3. SWAT-MODFLOW Model Coupled with Dynamic HRUs 2.3. SWAT‐MODFLOW Model Coupled with Dynamic HRUs* In view of the inability of the original SWAT model to effectively reflect complete or

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. In view of the inability of the original SWAT model to effectively reflect completeor 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. 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.

**Figure 3.** Production of dynamic hydrological response units (HRUs). **Figure 3.** Production of dynamic hydrological response units (HRUs). **Figure 3.** Production of dynamic hydrological response units (HRUs).

The operational flow chart of the coupled SWAT‐MODFLOW model based on dy‐ namic HRUs is shown in Figure 4. First, annual cycle simulations were performed by us‐ ing 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 subse‐ quent SWAT computations. The above process was conducted in nested loops until the The operational flow chart of the coupled SWAT‐MODFLOW model based on dy‐ namic HRUs is shown in Figure 4. First, annual cycle simulations were performed by us‐ ing 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 subse‐ quent 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. pled model could reflect the dynamic changes in land use/cover, the model was referred to as LU‐SWAT‐MODFLOW.

end of the simulation. Considering that the dynamic HRU‐based SWAT‐MODFLOW cou‐

*Water* **2021**, *13*, x FOR PEER REVIEW 5 of 20

**Figure 4.** Flowchart of the SWAT‐MODFLOW coupled model based on dynamic HRUs. **Figure 4.** Flowchart of the SWAT-MODFLOW coupled model based on dynamic HRUs.

#### *2.4. Model Establishment 2.4. 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 (Fig‐ ure 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. 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.

**Table 1.** Irrigation volume in the study region.


Forest irrigation 5400 6 April–November

Grassland irrigation 3600 5 April–November 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. From 2000 to 2018, the years when evident vegetation changes (restoration) occurred in the study region were 2005, 2008, 2015, and 2018. Accordingly, the land use data for 2000, 2005, 2008, 2015, and 2018 were used to generate dynamic HRUs. In contrast, land cover types in the other years were only weakly altered and therefore ignored. *Water* **2021**, *13*, x FOR PEER REVIEW 6 of 20

**Figure 5.** (**a**) Soil types in the study area and (**b**) sub‐basins input to the SWAT model. **Figure 5.** (**a**) Soil types in the study area and (**b**) sub-basins input to the SWAT model.

forestland and grassland. From 2000 to 2018, the years when evident vegetation changes (restoration) occurred in the study region were 2005, 2008, 2015, and 2018. Accordingly, **Figure 6.** Revegetation status of the study region in (**a**) 2000 and (**b**) 2018. **Figure 6.** Revegetation status of the study region in (**a**) 2000 and (**b**) 2018.

the land use data for 2000, 2005, 2008, 2015, and 2018 were used to generate dynamic HRUs. In contrast, land cover types in the other years were only weakly altered and there‐ fore ignored. Basin boundaries delineated by the SWAT model were considered as impermeable boundaries to groundwaterflow in the MODFLOW model, where the western and eastern outlets of the basin were considered as constant flow boundaries. The river network ex‐ tracted 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 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

to relevant studies [32], and the groundwater flow was conceptualised as a two‐dimen‐

The years 2000–2001 were used as the model warm‐up period, 2002–2011 as the model calibration period, and 2012–2018 as the model validation period. Vegetation growth parameters in the SWAT model were calibrated at HRU scales against the monthly 30 m resolution LAI data of 2002–2018 provided by the National Tibetan Plateau/Third Pole Environment Data Centre [29]. The dataset is the fusion of MODIS LAI and observed

ET simulated by the LU‐SWAT‐MODFLOW model was calibrated at sub‐basin scales against the ET data at a 0.1° × 0.1° resolution (Figure 8) provided by the National Tibetan Plateau/Third Pole Environment Data Centre [29]. The dataset is derived by employing a

**Figure 7.** Initial groundwater head for MODFLOW.

LAI in Qilian mountainous area (including the Qaidam Basin).

sional transient flow.

*2.5. Model Calibration*

*Water* **2021**, *13*, x FOR PEER REVIEW 7 of 20

(**a**) (**b**) **Figure 6.** Revegetation status of the study region in (**a**) 2000 and (**b**) 2018.

relevant studies [32], and the groundwater flow was conceptualised as a two-dimensional transient flow. to relevant studies [32], and the groundwater flow was conceptualised as a two‐dimen‐ sional transient flow.

Basin boundaries delineated by the SWAT model were considered as impermeable boundaries to groundwaterflow in the MODFLOW model, where the western and eastern outlets of the basin were considered as constant flow boundaries. The river network ex‐ tracted 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

**Figure 7.** Initial groundwater head for MODFLOW. **Figure 7.** Initial groundwater head for MODFLOW.

#### *2.5. Model Calibration 2.5. Model Calibration*

The years 2000–2001 were used as the model warm‐up period, 2002–2011 as the model calibration period, and 2012–2018 as the model validation period. Vegetation growth parameters in the SWAT model were calibrated at HRU scales against the monthly 30 m resolution LAI data of 2002–2018 provided by the National Tibetan Plateau/Third Pole Environment Data Centre [29]. The dataset is the fusion of MODIS LAI and observed LAI in Qilian mountainous area (including the Qaidam Basin). The years 2000–2001 were used as the model warm-up period, 2002–2011 as the model calibration period, and 2012–2018 as the model validation period. Vegetation growth parameters in the SWAT model were calibrated at HRU scales against the monthly 30 m resolution LAI data of 2002–2018 provided by the National Tibetan Plateau/Third Pole Environment Data Centre [29]. The dataset is the fusion of MODIS LAI and observed LAI in Qilian mountainous area (including the Qaidam Basin).

ET simulated by the LU‐SWAT‐MODFLOW model was calibrated at sub‐basin scales against the ET data at a 0.1° × 0.1° resolution (Figure 8) provided by the National Tibetan Plateau/Third Pole Environment Data Centre [29]. The dataset is derived by employing a ET simulated by the LU-SWAT-MODFLOW model was calibrated at sub-basin scales against the ET data at a 0.1◦ × 0.1◦ resolution (Figure 8) provided by the National Tibetan Plateau/Third Pole Environment Data Centre [29]. 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 dew‐ point temperature, wind speed, precipitation, and net radiation from the China Meteoro‐ logical Forcing Dataset. The dataset is validated in Northwest China and is proved to have good spatial and temporal performance [24]. Furthermore, relevant parameters (conduc‐ tivity and storage) of the MODFLOW model were calibrated against monthly recorded

groundwater table depth at numerous observation wells (Figure 5b) in the basin.

**Figure 8.** Locations of ET data points. **Figure 8.** Locations of ET data points.


