Application of SWAT Model with a Modified Groundwater Module to the Semi-Arid Hailiutu River Catchment, Northwest China

Abstract

In the original soil and water assessment Tool (SWAT) model (SWAT-O), the contributions of shallow aquifers and deep aquifers to streamflow are simulated using the linear reservoir method. The movement of groundwater was limited in the hydrological response unit which is a minimum calculation unit in the SWAT. However, this computational method may not be suitable for the areas where a groundwater system is complicated, and the river is predominately recharged by groundwater. In this paper, we proposed an enhanced groundwater module which divides shallow aquifers into upper and lower aquifers, integrates all the deep aquifers of a sub-basin into a regional aquifer, and simulates interactive water amount between lower aquifer and deep aquifer using water depth difference. The modified groundwater module was introduced to the original SWAT model, hereby referred to as SWAT-MG. The SWAT-MG and SWAT-O models were applied to the Hailiutu River catchment, which is a semi-arid wind sandy grass shoal catchment. Results showed that both models underestimated streamflow in peak flow, while the simulated streamflow of SWAT-MG was closer the observed values than that of SWAT-O. Three evaluation criteria (NSE, RSR, PBIAS) were applied to evaluate the performance of the models and the results showed that SWAT-MG had a better performance than SWAT-O. The baseflow index of Hailiutu River which was calculated by the results of SWAT-MG was 96.78%, which means the streamflow is predominately recharged by groundwater, and this conforms to the actual situation of Hailiutu River catchment. This indicates that a SWAT model with a modified groundwater module could better represent the groundwater flow behavior in the study area.

1. Introduction

The soil and water assessment tool (SWAT), developed by the United States Department of Agricultural Research Service, is widely used to simulate the water cycle and transportation of agricultural chemicals and sediment transportation [1,2,3,4,5]. In addition, it is used to quantify the relationship between climate variability and hydrological regime, and anthropogenic influence on a hydrological regime [6,7,8,9,10]. It can also be applied to predict the hydrological change impact of land use, land management practices, and climate change [11,12,13,14].

Despite its widespread application, the performance of the SWAT model varies significantly in different cases. Consequently, some researchers have modified the SWAT model to satisfy their specific requirements. Qi et al. [15] proposed an energy balance module for SWAT to improve the simulation accuracy for snowmelt. Chen et al. [16] modified the SWAT auto-irrigation functions to improve the simulation accuracy of agricultural irrigation management. Tsuchiya et al. [17] proposed an enhanced paddy simulation module to improve the simulation accuracy of paddy water balance. In the original SWAT, a hydrologic response unit (HRU) is a minimum calculation unit, and there is no exchange of water between different HRUs. In each HRU, two-layers aquifer model, shallow/unconfined aquifer and deep aquifer/confined aquifer, are used to represent the aquifer system, and a linear reservoir model to simulate groundwater flow. The water of a shallow aquifer could move to the deep aquifer, while the reverse process is not allowed. Using this computational method to simulate groundwater flow may lead to inaccurate results for the areas where a groundwater system is complicated, and the river is predominately recharged by groundwater. Shirmohammadi et al. [18] observed that the SWAT underestimates groundwater flow, especially during wet periods, in a piedmont physiographic region. Amatya et al. [19] found that the SWAT was unable to get accurate monthly outflows simulation results when it was applied to simulate streamflow in a karst watershed where runoff was dominated by groundwater. Cheng et al. [20] observed that the SWAT model was not effective in simulating the base flow in the Kuye River basin, a basin with coarse sand in the middle of the Yellow River watershed with a semi-arid climate. Similarly, Pfannerstill et al. [21] found out that the SWAT model performed poorly for low flow periods in North German lowlands. Wu et al. [22] also pointed out that the SWAT model underestimates groundwater flow in a dry year. Several studies have indicated that outflow from groundwater storage is not linearly proportional to storage [23,24]. Therefore, some researchers have modified the groundwater module of the SWAT model to improve the simulation accuracy. Wang et al. [25] modified the linear reservoir model, which is used to calculate groundwater flow, to a nonlinear reservoir model in a Karst area. Similarly, Gan and Luo [26] used a nonlinear reservoir model to calculate base flow in glacier and snowmelt-dominated basins. Some researchers used multi linear reservoirs to reflect the nonlinearity of outflow. Pfannerstill et al. [21] divided shallow groundwater storage into fast storage and slow storage in a lowland catchment.

Using the modified groundwater module could achieve better results in reproducing streamflow processes, especially during low flow periods. In their study area, in the dry season, the observed discharges decreased to a very low value, even zero. Nonetheless, that is not the real situation of the streamflow process during the dry season in some catchments where baseflow is the major component of streamflow. In these catchments, a relatively steady discharge is typically maintained even in the dry season. The groundwater system of these catchments is generally more complex than other catchments. In order to obtain a reasonable representation of groundwater flow, the SWAT model coupled with the modular three-dimensional finite-difference groundwater flow (ModFlow) model is effective [27,28,29,30,31]. This approach also requires more detailed input data and extensive computational time when obtaining better groundwater simulation accuracy [31].

In this paper, a modified groundwater module, which divides a shallow aquifer into upper and lower aquifers and integrates all the deep aquifers of a sub-basin into a regional aquifer to simulate groundwater flow, was incorporated to the original SWAT model. This modified SWAT model was applied in a semi-arid wind-sandy grass shoal catchment—the Hailiutu River catchment in northwest China—to evaluate its performance. In this catchment, groundwater has been found to be the primary source of recharge for streamflow [20,21]. The original and modified model results were compared. The results indicated whether the modification of the groundwater module of the SWAT model can improve the simulation of the streamflow which is steadily recharged by groundwater during the dry season (October to May) in the semi-arid wind-sandy grass shoal catchment.

2. Materials and Methods

2.1. Study Area

The Hailiutu River catchment (38°02′–38°50′, 108°37′–109°14′), located in the middle section of the Yellow River basin in northwest of China, covers an area of approximately 2500 km² (Figure 1). This catchment is a typical wind-sandy grass shoal area and a section of Maowusu sandy land. The catchment has a relatively flat topography, with altitude ranging from 1011 m to 1475 m, and the main landform is sand dunes. The catchment belongs to a semi-arid climate, the annual average precipitation is about 380 mm, annual mean temperature is about 8.5 °C, and annual average sunshine time is about 2780 h. The percentage of urban lands in the Hailiutu River basin is small and most are scattered in the basin. The population density of the study area is approximately 30 people/km². In the study area, farmers practice agriculture substantially, with maize as the primary crop [32,33]. Most agricultural lands are located in the river valley and in the Bulang River sub-catchment. Studies have shown that the underground of the study area is rich in mineral resources such as coal, oil, and natural gas [34,35]. Therefore, in recent years, the government has continuously increased its investment in the construction of base-building in this region for the country’s economic growth. For environmental protection, the government implemented some projects such as ‘Three North Forest Shelterbelts’ and took the policy of ’Closing Sandy land and Forbidding Herding’ to increase shrub and grass land by decreasing bare land in the study area [36]. The Hailiutu River, which is the major river in the catchment, is one of the tributaries of the Wuding River which is a major tributary of the Yellow River. The streamflow of the Hailiutu River is primarily recharged by groundwater [33,34,35,36]. Yang et al. [32] found that the hydrological processes are dominated by the evapotranspiration and direct infiltration of precipitation. Yang et al. [32] also observed with isotopic separation method that up to 96.4% of the total streamflow was composed of baseflow in the Hailiutu River catchment. Zhou et al. [33] observed that only about 5% of the Hailiutu River catchment has a surface slope steeper than 5°and the aquifer could get substantial recharge from precipitation infiltration due to permeable Aeolian sand formation. The flat land surface and low drainage density (0.18 km/km² from Zhou et al.) assist the Hailiutu River in maintaining a relatively steady discharge even in the dry season (October to May).

The geological formation in Hailiutu River catchment can be divided into four strata: (1) the Holocene Aeolian sand with a thickness of about 0 to 50 m; (2) the upper Pleistocene Shalawusu sandstone formation of Quaternary age with a thickness of 5 to 90 m and mainly consists of sand of the lacustrine origin; (3) the Luohe sandstone of Cretaceous with a thickness of 180 to 330 m and overlain by the Shalawusu formation; and (4) the bedrock, which consists of impermeable Jurassic sedimentary with about 250 m [32,37,38]. Consequently, the complicated geological formation leads to a complex groundwater cycle in the Hailiutu River catchment (Figure 2b).

2.2. Data Collection

In this paper, the required data include digital elevation model (DEM), land use data, soil data, meteorological data, daily precipitation and streamflow data (

Table 1. Summarized information on the data used for this study.

Data	Description of Data	Data Sources
Topographic	30 × 30 m resolution digital elevation model (DEM)	Geospatial Data Cloud of China
Soil map/layer	1 × 1 km resolution map; soil layer attributes for each soil layer	Environmental and Ecological Science Data Center for West China
Land use	30 × 30 m resolution map	Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences
Daily meteorological data	Daily wind speed, minimum and maximum temperature and relative humidity from 1970 to 1985	China Meteorological Sharing Service System.
Daily rainfall and streamflow	Daily precipitation and streamflow data in/around watershed	Yellow River Conservancy Commission

). The DEM data, with a spatial resolution of 30 × 30 m, were obtained from the Geospatial Data Cloud of China (Figure 1). Land use data for 1980 s (1:100,000) were provided by Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (Figure 1). The Hailiutu River catchment consists of five major land use types, including agricultural land, shrub and grass land, water body, urban land, and bare land. The soil data were obtained from the Harmonized Word Soil Database (HWSD) supplied by the Environmental and Ecological Science Data Center for West China. Meteorological data from 1970 to 1985, were acquired from the China Meteorological Sharing Service System. These data included maximum and minimum temperature, wind speed and relative humidity. The observed daily precipitation and streamflow of the Hailiutu River were provided by the Yellow River Conservancy Commission. It was observed that the time series of streamflow of the Hailiutu River catchment showed an abrupt recession around the year 1985. This was due to increased irrigation water and river diversion water [36,39]. Therefore, in this paper, the period prior to 1986 (1970–1985) was selected to assess the SWAT model.

Using the terrain processing module of Arc-Hydro Tools software with ArcGIS interface, based on D8 algorithm and DEM data, the basin and sub-basin boundary of study area were delineated with a threshold area of 3000 ha. In the watershed delineation and subdivision processes, the actual river system was applied to force the generated streams of the SWAT model to follow the actual stream. The land use data, slope data, soil data were applied to delineate the HRU. As a result, Hailiutu River catchment was divided into 10 sub-basins and 138 HRUs in the SWAT model.

2.3. Introduction of SWAT Model

2.3.1. Original SWAT Model

In the SWAT model, the source of streamflow was mainly divided into four components, namely surface flow (Q_sf), tile flow ($Q_{tf}$), lateral flow ($Q_{lf}$) and groundwater flow ($Q_{gw}$), as illustrated in Figure 3. The groundwater flow could be sub-divided into shallow groundwater flow and deep groundwater flow which come from shallow and deep aquifers, respectively.

Simulation of $Q_{sf}$, $Q_{tf}$ and $Q_{lf}$ are represented in detail in the theoretical documents of the SWAT model [40] and thus will not be described repeatedly in this paper. The original SWAT model partitions underground saturated aquifer into shallow aquifer and deep aquifer. The water balance of shallow aquifer can be expressed as

$$S_{sh}^i=S_{sh}^{i-1}+W_{rg,sh}^i-Q_{gw,sh}^i-W_{rp}^i-W_{pump,sh}^i$$

(1)

where, $S_{sh}^i$ and $S_{sh}^{i-1}$ refer to the amount of water stored in the shallow aquifer on day $i$ and $i-1$, respectively, $Q_{gw,sh}^i$ is the groundwater flow from shallow aquifer on day i, $W_{rg,sh}^i$ represents the amount of recharge entering the shallow aquifer on day $i$, $W_{rp}^i$ represents the amount of water entering the soil zone for evaporation and transpiration on day i, and $W_{pump,sh}^i$ is the amount of water removed from the shallow aquifer by pumping on day i. The linear reservoir model expressed as Equation (2), which assumes that there is no re-evaporation and no pumping, is used to calculate groundwater flow.

$$Q_{gw,sh}^i={\alpha }_{gw,sh}{S}_{sh}^i$$

(2)

where, ${\alpha }_{gw,sh}$ is the groundwater recession constant of shallow aquifer. The small groundwater recession constant contributes to limit the groundwater flow into the stream. The groundwater from shallow aquifer contributed to the stream on day i is expressed as:

$$Q_{gw,sh}^i=Q_{gw,sh}^{i-1}{e}^{-{\alpha }_{gw,sh}}+W_{rg,sh}^i(1-e^{-{\alpha }_{gw,sh}})$$

(3)

The amount of recharge entering the shallow aquifer on day i is calculated as:

$$W_{rg,sh}=(1-{\beta }_{dp})W_{rg}^i$$

(4)

$$W_{rg}^i=(1-e^{-\frac{1}{{\delta }_{gw,sh}}})W_{seep}+e^{-\frac{1}{{\delta }_{gw,sh}}}{W}_{rg}^{i-1}$$

(5)

$$W_{rg,deep}={\beta }_{dp}{W}_{rg}^i$$

(6)

where $W_{rg}^i$ and $W_{rg}^{i-1}$ are the amount of recharge entering the aquifers (shallow and deep) on day $i$ and $i-1$, respectively, ${\delta }_{gw,sh}$ is a coefficient that reflects the delay time of recharge for shallow aquifer, and $W_{seep}$ is the total amount of water exiting the bottom of the soil profile on day $i$, and ${\beta }_{dp}$ is a coefficient that reflects the amount of percolation water moving into the deep aquifer.

Similar to the Equation (3), the groundwater from a deep aquifer contributed to the stream on day i is calculated as:

$$Q_{gw,deep}^i=Q_{gw,deep}^{i-1}{e}^{-{\alpha }_{gw,deep}}+W_{rg,deep}^i(1-e^{-{\alpha }_{gw,depp}})$$

(7)

It is worth mentioning that the SWAT model uses threshold value ($Q_{trd,sh}$) to decide whether groundwater contributes to main channel, and if shallow aquifer storage falls below $Q_{trd,sh}$, $Q_{gw,sh}^i$ is set to zero. Similarly, revap (the process that draws water from aquifer for evaporation and transpiration in the theoretical documentation of SWAT [40] is allowed to occur only if the amount of water stored in the shallow aquifer exceeds a threshold value ($r_{trd,sh}$).

2.3.2. Modified SWAT Model

In the original SWAT, based on independent HRUs, two groundwater aquifers were adopted to calculate groundwater flow. In this way, it is difficult to obtain a better simulation results in some areas where a groundwater system is complicated, and the stream is predominately recharged by groundwater. Using more active aquifers in the SWAT model to calculate groundwater flow can provide a more realistic representation of the hydrological cycle. However, more aquifers would add complexity to the model and make the calibration of parameters more difficult. According to the geological formation of the study area (see Section 2.1), similar with Pfannerstill et al. [21], we sub-divided the shallow aquifer into two active aquifers in the groundwater module of SWAT model, i.e., upper aquifer and lower aquifer, and the deep aquifer was remained in the modified SWAT to simulate the groundwater flow from confined aquifer.

The groundwater flow from upper and lower aquifers are calculated as following (modified from Equation (3)):

$$Q_{gw,u}^i=Q_{gw,u}^{i-1}{e}^{-{\alpha }_{gw,u}}+W_{rg,u}^i(1-e^{-{\alpha }_{gw,u}})$$

(8)

$$Q_{gw,l}^i=Q_{gw,l}^{i-1}{e}^{-{\alpha }_{gw,l}}+W_{rg,l}^i(1-e^{-{\alpha }_{gw,l}})$$

(9)

where, $Q_{gw,u}^i$ and $Q_{gw,l}^i$ represent the groundwater flow from upper aquifer and lower aquifer into the stream on day i, respectively; ${\alpha }_{gw,u}$ and ${\alpha }_{gw,l}$ mean the groundwater recession constant of upper aquifer and lower aquifer, respectively; $W_{rg,u}^i$ and $W_{rg,l}^i$ represent the amount of recharge entering the upper aquifer and lower aquifer on day $i$, respectively.

In the Hailiutu River basin, there is no regional aquitard between the unconfined aquifer and confined aquifer. Therefore, the two aquifers have some hydraulic connection. In the original SWAT, a fixed positive value was adopted to reflect this interaction (see more in Equation (4)). The water of the shallow aquifer could recharge the deep aquifer; however, the water of the deep aquifer recharge to shallow aquifer was prohibited. In order to better reflect the groundwater cycle, this process was improved in the modified SWAT. Considering that both the upper and lower aquifer belong to the shallow aquifer, the interaction of these two aquifers was simulated according to the method of the original SWAT.

The recharge to the upper aquifer on a given day is calculated as:

$$W_{rg,u}^i=(1-{\beta }_l)W_{rg}^i$$

(10)

$$W_{seep,u}^i={\beta }_l{W}_{rg}^i$$

(11)

where, ${\beta }_l$ is a coefficient that reflects the amount of percolation water that moves into the lower aquifer, $W_{seep,u}^i$ is the leakage of water from upper aquifer on day i. Similar to Equation (5), the recharge to lower aquifer on a given day is expressed as

$$W_{rg,ul}^i=(1-e^{-\frac{1}{{\delta }_{gw,l}}})W_{seep,u}^i+e^{-\frac{1}{{\delta }_{gw,l}}}{W}_{rg,ul}^{i-1}$$

(12)

where $W_{rg,ul}^i$ and $W_{rg,ul}^{i-1}$ mean the amount of recharge from the upper aquifer on day $i$ and $i-1$, respectively. ${\delta }_{gw,m}$ is a coefficient that reflects the delay time of recharge for lower aquifer.

The interaction of lower aquifer and deep aquifer was simulated according to the method of ModFlow:

$${\beta }_d=\frac{1}{\frac{Th{k}_l/2}{K_l}+\frac{Th{k}_d/2}{K_d}}$$

(13)

where, Thk means the thickness of aquifer (m), K means vertical hydraulic conductivity (m/d), subscript $l$ and $d$ mean lower aquifer (unconfined aquifer) and deep aquifer (confined aquifer), respectively. Some researchers adopt ModFlow to simulate the groundwater cycle in the Hailiutu River basin. The hydraulic conductivity of unconfined aquifer was calibrated to be 2–12 m/d. The hydraulic conductivity of confined aquifer was calibrated to be 0.25 m/d. The thickness of the confined aquifer was set to be an average value (250 m). Therefore, the value of ${\beta }_d$ could be estimated as 0.002 d⁻¹ in this study.

The specific yield of unconfined aquifer and specific storage of confined aquifer were assumed as one. In this case, the amount of water of interaction of these two aquifers could be calculated as:

$$W_{rg,deep}={\beta }_d(h_l-h_d)\Delta t$$

(14)

where, $h_l$ and $h_d$ represent the water depth of lower aquifer and deep aquifer (mm), $\Delta t$ is calculation step (d).

The groundwater flow system could be divided into local groundwater flow system, intermediate groundwater flow system and regional groundwater flow system [32]. The local groundwater flow system is mainly affected by topography, it originates from local high topographic area and ends in the lower topographic area. The recharge and discharge area of intermediate and regional groundwater flow are defined by a series of topographic highs and lows. In this case, the upper aquifer and lower aquifer, which are based on an independent HRU could be considered as a local groundwater flow system. The deep aquifers of HRUs in each sub-basin were considered as a regional groundwater system in the modified SWAT. Studies show that both surface water and groundwater flow from northwest to southeast of the Hailiutu River basin (Figure 2a) [41,42]. Therefore, the interaction of each regional groundwater in the modified SWAT could be simulated according to the route of runoff concentration (Figure 4).

Step 1. Calculate the average depth ($S_{d,sub}$) of sub-basin and add the amount of regional aquifer recharge from upstream to $S_{d,sub}$.

$$S_{d,sub}=\frac{{\displaystyle \sum _{i=1}^N{S}_{d,hru,i}{A}_{hru,i}}}{A_{sub}}+\frac{Q_{rchrg,sub}}{A_{sub}}$$

(15)

where, $N$ is the number of HRU in the sub-basin, $S_{d,hru,i}$ and $A_{hru,i}$ mean depth and area of deep aquifer of i-th HRU, respectively, $A_{sub}$ is the area of sub-basin, $Q_{rchrg,sub}$ is amount of regional aquifer recharge water from upstream.

Step 2. Determine whether to trigger a ’recharge to river’ event. If no, then go to Step 3. If yes, then calculate the amount of recharge water based on linear reservoir model:

$$Q_{rvr}^{}={\alpha }_{gw,sub}\cdot (S_{d,sub}-S_{min,sub})A_{sub}$$

(16)

where, ${\alpha }_{gw,sub}$ is the recession constant of regional aquifer. $S_{min,sub}$ is a threshold value to decide whether regional groundwater contributes to the river.

Step 3. Calculate amount of flow out water for regional groundwater cycle:

$$Q_{out}^{}={\alpha }_{gw,out,sub}{S}_{d,sub}{A}_{sub}$$

(17)

where, ${\alpha }_{gw,out,sub}$ is outflow coefficient of regional aquifer for regional groundwater cycle. Meanwhile, in order to reflect the travel time of regional groundwater between two sub-basins. Refer to Equation (5), the amount of water actually participated in the regional groundwater cycle on i day is calculate as:

$$\{\begin{array}{l}{Q}_{out2,sub}^i=(1-e^{-\frac{1}{{\delta }_{gw,sub}}})Q_{out,sub}^i+e^{-\frac{1}{{\delta }_{gw,sub}}}{Q}_{out2,sub}^i\\ {\delta }_{gw,sub}={\theta }_{delay,sub}Di{s}_{sub}\end{array}$$

(18)

where, ${\delta }_{gw,sub}$ is a coefficient that reflects the delay time of recharge from one sub-basin to another sub-basin, $Di{s}_{sub}$ is the distance between two interrelated sub-basins, ${\theta }_{delay,sub}$ is an dimensionless positive value.

Step 4. Based on the amount of recharge to river and flow out for regional groundwater cycle, update the depth of deep aquifer in each HRU in the sub-basin.

2.4. Evaluation Criteria

In this paper, three evaluation criteria which proposed by MORIAS [43] were applied to evaluate the performance of SWAT. These evaluation criteria are NSE, RSR and PBIAS and calculated by Equations (19)–(21):

$$NSE=1-\frac{{\displaystyle \sum _{t=1}^n(}{O}_t-S_t{)}^2}{{\displaystyle \sum _{t=1}^n(}{O}_t-O_m{)}^2}$$

(19)

where, O and S represent the observed data and simulated data, respectively; subscript m represents average value of data.

$$RSR=\frac{\sqrt{{\displaystyle \sum _{t=1}^n(}{O}_t-S_t{)}^2}}{\sqrt{{\displaystyle \sum _{t=1}^n(}{O}_t-O_m{)}^2}}$$

(20)

$$PBIAS=\frac{{\displaystyle \sum _{t=1}^n{O}_t-S_t}}{{\displaystyle \sum _{t=1}^n{O}_t}}\times 100$$

(21)

3. Results and Discussion

For model evaluation, the whole dataset was divided into calibration period (1970–1980) and validation period (1981–1985). A warm-up phase from 1970 to 1973 was selected to achieve the steady state for the SWAT model. Studies have shown that the contribution of confined aquifer to streamflow only occurred in the downstream of Hailiutu River basin [42,44]. In the original SWAT, there is no threshold for outflow of the deep aquifer. In order to reflect the actual situation of the Hailiutu River basin, the groundwater recession value of deep aquifer was set as a minimum value, except for the sub-basin, which was located in the downstream of study area. The SUFI-2 method and daily observed streamflow data were used to calibrate both the original SWAT model (SWAT-O) and modified SWAT model (SWAT-MG). The calibration results of parameters for original and modified models are presented in

Table 2. The parameters’ values for the calibration of the original and modified SWAT models.

parameters	Initial Range		Calibrated Range				Optimal Value
	Up	Low	SWAT-O		SWAT-MG		SWAT-O	SWAT-MG
			Up	Low	Up	Low
$r\_CN2\mathbb{F}$	−0.65	−0.10	−0.53	−0.51	−0.55	−0.52	−0.51	−0.54
$v\_esco$	0.35	0.90	0.73	0.74	0.55	0.6	0.73	0.56
$a\_awc$	−0.035	−0.025	−0.025	−0.024	−0.027	−0.026	−0.025	−0.027
$v\_lat\_ttime$	0.50	5.00	2.92	3.24	2.52	2.70	3.04	2.59
$r\_sol\_k$	−0.80	−0.20	−0.31	−0.28	−0.71	−0.68	−0.29	−0.69
$v\_r_{trd,sh}$	30.0	70.0	34.5	36.7	72.2	74.4	35.7	73.7
$v\_g{w}_{revap}$	0.04	0.20	0.16	0.17	0.16	0.18	0.17	0.17
$\mathrm{v}\_{\alpha }_{gw,sh}$	0.0001	0.01	0.0006	0.0008	/	/	0.0007	/
$v\_{\alpha }_{gw,deep}$	0.0001	0.01	0.0079	0.0083	/	/	0.0081	/
$\mathrm{v}\_{\delta }_{gw,sh}^{}$	250	600	322	356	/	/	354	/
$\mathrm{v}\_Q_{trd,sh}$	10.0	30.0	11.8	12.5	/	/	12.2	/
$\mathrm{v}\_{\beta }_{gw,sh}$	0.30	0.60	0.38	0.42	/	/	0.40	/
$v\_{\alpha }_{gw,up}$	0.300	0.700	/	/	0.461	0.492	/	0.466
$v\_{\alpha }_{gw,low}$	0.0010	0.0100	/	/	0.0022	0.0025	/	0.0023
$v\_{\alpha }_{gw,out,sub}$	0.001	0.01	/	/	0.0007	0.0008	/	0.0007
$v\_{\alpha }_{gw,sub}$	0.0001	0.01	/	/	0.0012	0.0013	/	0.0013
$v\_{\delta }_{gw,up}^{}$	0.50	2.00	/	/	1.25	1.62	/	1.51
$v\_{\delta }_{gw,low}$	150	500	/	/	502	541	/	534
$v\_{\theta }_{gw,sub,delay}$	1.00	10.00	/	/	2.65	3.22	/	2.91
$v\_Q_{trd,up}$	10.00	35.00	/	/	13.78	17.36	/	14.28
$v\_{\beta }_{gw,up}$	0.950	0.985	/	/	0.978	0.982	/	0.979

Note: r means to multiply by original value, a means to add or subtract original value and v means to replace original value.

.

Compared the calibrated range of sol_k of two SWAT version, the sol_k values of SWAT-MG are smaller than SWAT-O. In this case, compared with SWAT-O, the lateral flow of SWAT-MG decreased and will allow more water to recharge groundwater. In the SWAT-O, a fixed value was applied to reflect the recharge of deep aquifer and the calibrated range was 0.38–0.42. In this way, shallow aquifer and deep aquifer would receive recharge at the same time, the difference is solely the amount of recharge water. In the SWAT-MG, the interaction of shallow aquifer and deep aquifer was simulated by the water depth difference of these two aquifers. In this case, the shallow aquifer will receive recharge earlier than deep aquifer. Consequently, the simulated water depth of shallow aquifer in the SWAT-MG will be more than that which in SWAT-O. Under such a circumstance, the calibrated groundwater evaporation threshold value ($r_{trd,sh}$) in SWAT-MG would more than that in SWAT-O.

Groundwater recession constant and delay time are important parameters which affect the recharge and outflow of groundwater aquifer. In the study area, the calibrated groundwater recession constant of deep aquifer (${\alpha }_{gw,deep}$) has a bigger value, however, the outflow of deep aquifer was only allowed in the downstream of the study area. Therefore, shallow aquifer was the main recharge source of streamflow. In the SWAT-O, the calibrated range for groundwater recession constant (${\alpha }_{gw,sh}$) and delay time (${\delta }_{gw,sh}$) of shallow aquifer are 0.0006–0.0008 and 322–356, respectively. In this case, the groundwater could drain into the stream for a long period while the responding time to recharge events is long. In the SWAT-MG, based on geological formation of the study area, the shallow aquifer was sub-divided into upper aquifer and lower aquifer. The calibrated range of groundwater recession constant (${\alpha }_{gw,up}$) and groundwater delay time (${\delta }_{gw,up}$) of upper aquifer were 0.461–0.492 and 1.25–1.62, respectively. These values indicate upper aquifer will quickly receive substantial recharge and drain groundwater into the stream during rainfall events. The calibrated range of groundwater recession constant (${\alpha }_{gw,low}$) and groundwater delay time (${\delta }_{gw,low}$) of lower aquifer were 0.0022–0.0025 and 502–541, respectively. Under such a circumstance, the groundwater from lower aquifer could continuously recharge stream with negligible influence from the rainfall events. The calibrated range of groundwater recession constant of regional aquifer (${\alpha }_{gw,sub}$) is similar to the calibrated range of ${\alpha }_{gw,low}$. That means the outflow characteristic of regional aquifer is similar to the lower aquifer in the study area.

Based on the calibrated range of parameters of Table 2, 1000 group parameters were generated by the Latin hypercube sampling method for SWAT-O and SWAT-MG, respectively. These group parameters were applied to drive two SWAT versions to investigate the difference of performance in SWAT-O and SWAT-MG. Figure 5 illustrated the comparison of 1000 group evaluation criteria values between the SWAT-O and SWAT-MG during the calibration and validation period, respectively. As shown in the Figure 5, in both SWAT versions, the distribution of PBIAS is more even than the other two evaluation criteria. For the SWAT-O, the performance had a significant difference between the calibration and validation period. For instance, the median of NSE was about 0.47 in the calibration period while it decreased to about −0.62 in the validation period. Different from SWAT-O, SWAT-MG has a similar performance during calibration and validation period. For streamflow simulation, the model performance can be identified as satisfactory if NSE > 0.5, RSR < 0.7 and −25 < PBIAS < 25 [43]. That means the performance of SWAT-MG could be considered as satisfactory in both calibration and validation period. In general, SWAT-MG has a better performance than SWAT-O on all evaluation criteria.

Combined with the performance of calibration and validation period, the best 40 simulated streamflow of SWAT-O and SWAT-MG were selected to calculate the variation range of flow duration curve (FDC). For comparison, the FDC of observed streamflow was also presented (Figure 6). In both calibration and validation period, SWAT-O and SWAT-MG have similar performance in the FDC segment from 0% to 10%, both SWAT versions tend to underestimate the peak flow. Compared with SWAT-O, SWAT-MG produces less underestimation. In the FDC segment from 10% to 60%, both SWAT versions have similar performance in the calibration period, however, the SWAT-MG performs better than SWAT-O in the validation period. In the FDC segment from 60% to 100%, SWAT-MG could better reproduce the variation of observed FDC while SWAT-O tended to overestimate the streamflow in the calibration period. In the validation period, SWAT-MG mainly overestimates the streamflow while SWAT-O significantly underestimates the streamflow. Studies shown that, in the 1980s, the agricultural land has increased and amount of irrigation water also increased [33,39]. These variations were difficult to exactly reflect in the SWAT model. Consequently, SWAT-MG tended to overestimate the observed streamflow in the validation period.

Based on the aforementioned analysis, a group of parameters that could obtain the best performance was selected as the optimal parameters for SWAT-O and SWAT-MG (Table 2). The best simulated streamflow of each SWAT version could be obtained by using optimal parameters to run the model. Figure 7 shows the scattered comparison of the simulated streamflow of two SWAT versions with observed streamflow. In the study period, the average annual streamflow was about 2.73 m³/s. The daily streamflow, which was more than 4 m³/s generally occurred in the rainy season and this value located at about 5% of FDC (Figure 6). Therefore, the daily streamflow could be considered as peak flow if it is more than 4 m³/s in the study area. It is indicated from Figure 7, although two SWAT versions underestimated streamflow during peak flow, simulation results show that SWAT-MG produced less underestimation, which means SWAT-MG is better than SWAT-O. For the non-peak flow, the simulated points of two SWAT versions were distributed relatively uniform on the 1:1 line while the simulated points of SWAT-MG were closer to 1:1 line than SWAT-O.

Figure 8 illustrated the comparison between the daily observed and simulated streamflow. As shown in Figure 8, in the calibration period, SWAT-O and SWAT-MG could simulate the variation tendency of streamflow, nonetheless, the real flow fluctuation is larger than the simulated ones. At the beginning of validation period, the simulated streamflow of SWAT-O presented a steeply decreased trend and obviously underestimated the streamflow, while SWAT-MG could better reproduce the variation of streamflow. Figure 9 illustrated the comparison between the monthly observed and simulated streamflow. Similar to the results of daily scale, SWAT-MG outperformed SWAT-O in both the calibration and validation period. The values of NSE, RSR, PBIAS of SWAT-O and SWAT-MG are listed in

Table 3. Performance for calibration and validation of two different SWAT versions on daily and monthly scales.

	Criteria	NSE		RSR		PBIAS
SWAT Version
		Daily	Monthly	Daily	Monthly	Daily	Monthly
Calibration period
SWAT-O		0.47	0.60	0.73	0.64	0.67	0.67
SWAT-MG		0.58	0.81	0.65	0.43	1.21	1.21
Validation period
SWAT-O		-0.52	-1.88	1.23	1.69	26.6	26.6
SWAT-MG		0.57	0.71	0.65	0.53	-0.93	-0.93

for both the calibration and validation periods at daily and monthly scales. As shown in Table 3, all criteria values of SWAT-MG were higher than SWAT-O, which indicated that SWAT-MG performed better than SWAT-O. The variation of NSE, which emphasizes high flow errors [45], was more apparent than the others. This suggested that the modified groundwater module contribute to improve performance of peak flows simulation, which is in accordance with Figure 7.

Table 4. Different components of simulated streamflow by SWAT-MG (m³/s).

Year	Streamflow	Base Flow			SURFACE FLOW	Base Flow /Streamflow (%)
		Deep Flow	Lower Flow	Upper Flow
1974	2.57	0.13	2.33	0.04	0.07	97.38
1975	2.74	0.17	2.45	0.04	0.08	97.12
1976	2.68	0.18	2.36	0.04	0.08	96.9
1977	2.75	0.19	2.39	0.07	0.11	96.13
1978	2.93	0.19	2.53	0.08	0.13	95.72
1979	2.82	0.2	2.51	0.04	0.07	97.34
1980	2.58	0.19	2.32	0.02	0.05	98.17
1981	2.61	0.19	2.27	0.07	0.08	96.92
1982	2.63	0.19	2.28	0.08	0.08	96.87
1983	2.61	0.19	2.33	0.02	0.06	97.72
1984	2.8	0.2	2.36	0.1	0.13	95.18
1985	2.99	0.21	2.62	0.05	0.11	96.18
Average	2.72	0.19	2.4	0.05	0.09	96.78

shows the annual average values of different streamflow components, simulated by SWAT-MG, from 1970 to 1985. As shown in the Table 4, base flow accounted for over 90% of the streamflow in general. On average, the proportion of base flow to streamflow was up to 96.78%, which reveals that groundwater was the primary source of streamflow, and that conforms to the actual situation of Hailiutu River catchment. Groundwater flow from upper aquifer and surface flow (the summation of surface flow, tile flow and lateral flow) have the 5.2% of the contribution to the streamflow, which were 0.09 m³/s and 0.05 m³/s respectively. Groundwater flow from deep aquifer accounts for 6.83% of the streamflow. Low flow from lower aquifer provides the basic discharge of river which accounts for 87.96% of streamflow and 90.88% of base flow.

The modified SWAT model adopted water depth difference of shallow aquifer and deep aquifer to calculate the amount of water interaction between these two aquifers. Similar to the surface flow, the groundwater from regional aquifer was allowed from upstream flow to downstream of study area. Figure 10 presents the spatial distribution of annual seep water from shallow aquifer simulated by SWAT-MG. The positive value of Figure 10 means the shallow aquifer recharges to deep aquifer while the negative value means the deep aquifer recharges to shallow aquifer. As shown in the Figure 10, the deep aquifer will receive recharge from shallow aquifer in the upstream of study area and recharge to shallow aquifer at the sub-basin which belong to intersection of stream. At the downstream of Hailiutu River basin, the groundwater of regional aquifer was allowed to drain to stream. Therefore, the water depth of deep aquifer of this sub-basin was smaller than the shallow aquifer, and the deep aquifer will receive recharge from shallow aquifer.

As mentioned in Section 2.1, the Hailiutu River is mainly recharged by groundwater, the discharge of Hailiutu River would keep a relatively stable value instead of decreasing to a low value in the dry season (October to May). The contribution of groundwater during rainfall events accounts for approximately 74.8% of the total discharge [32]. The aquifers of the original SWAT model were unable to quickly receive substantial recharge the streamflow during rainfall events while continually recharging to streamflow. The SWAT-MG used an upper aquifer to reflect the fact that streamflow was mainly recharged by groundwater during rainfall events, lower aquifer and regional aquifer to reflect the phenomenon that streamflow was continuously recharged by groundwater, even in the dry period. Moreover, SWAT-MG adopted the difference of water depth between shallow aquifer and deep aquifer to calculate the interactive water amount. According to the simulation results, SWAT-MG outperformed SWAT-O on both daily and monthly scales by all criteria during both calibration and validation periods. The proportion of base flow to total flow was about 96.78% by SWAT-MG, and this value is very close to the value observed by Yang et al. (96.4%) with isotopic technology. The study indicated that modification of the SWAT groundwater module is physically reasonable and its simulation is better than the original SWAT model and calculated high percentage of baseflow (96.78%) in the streamflow of the Hailiutu River catchment can be confirmed by isotopic technology.

Studies have shown that irrigation and river diversion are two factors which affect the variation of streamflow [33,39,41]. In this study, the built SWAT model could simulate the irrigation water of agricultural land through an auto-irrigation module. Yang et al., based on the data before about four decades, calibrated the parameters of ModFlow, then used it to simulate the streamflow of the current period [42]. Yang et al. found that the difference in simulated streamflow and observed streamflow was attributed to the river diversion water [42]. There are many river diversions that were built during the past three decades affects the simulation accuracy of streamflow [39,42]. Therefore, we selected the period of 1970–1985 as a natural period to evaluate the performance of the modified SWAT model. However, if the amount of diverted water is considered in the modified SWAT model, then it can still be applied in the analysis of streamflow variation in the current period since land use change has less impact on streamflow [46].

4. Conclusions

In this paper, a modified groundwater module was proposed by dividing shallow aquifer into the upper aquifer and lower aquifer and integrating all the deep aquifers of a sub-basin into a regional aquifer. The interaction of shallow aquifer and deep aquifer was simulated by water depth difference. The modified and original SWAT model were applied to simulate the streamflow in the semi desert and primarily groundwater recharged Hailiutu River catchment.

SWAT-MG had a better performance than SWAT-O on both daily and monthly scales by all criteria. The groundwater flow from upper aquifer can quickly respond to rainfall events and the lower aquifer could continuously recharge the streamflow. The regional aquifer could reflect the groundwater flow of the confined aquifer in the Hailiutu River basin. Therefore, SWAT-MG performed better in reproducing the groundwater flow processes and achieved a better streamflow simulation accuracy, especially in the peak flow period. With the SWAT-MG, the base flow proportion in the total streamflow was calculated as 96.78%, this value confirmed that the river is mainly recharged by groundwater, and that conformed to the actual situation of Hailiutu River catchment. This value closely corresponds to the value determined by a prior study carried out using isotope technology. In general, the modified SWAT model resulted in an improved representation of the streamflow processes of Hailiutu River catchment.