Influence of Land Use Change on Hydrological Cycle: Application of SWAT to Su-Mi-Huai Area in Beijing, China

The human activities and urbanization process have changed the underlying surface of urban areas, which would affect the recharge of groundwater through rainfall infiltration and may further influence the groundwater environment. Accordingly, it is imperative to investigate the variation of hydrological cycle under the condition of underlying surface change. Based on the high-precision remote sensing data of 2000, 2005, 2010 and 2015, and Soil and Water Assessment Tool (SWAT) model, this work firstly studied the land use change and the corresponding changes in runoff generation mechanism and rainfall infiltration coefficient in Su-Mi-Huai area, Beijing, China. Meanwhile, SWAT-MODFLOW semi-loose coupling model was applied to analyze the water balance in the study area in typical hydrological years. The results showed that the area of the construction land (urban and rural residential land) increased by 1.04 times from 2000 to 2015, which is mainly attributed to the conversion of cultivated land to construction land in the plain area. This change caused the runoff in the area to increase by 7 × 106 m3, the runoff coefficient increased by 17.9%, and the precipitation infiltration coefficient was less than the empirical value determined by lithology. Compared with 2000, the average annual precipitation infiltration coefficient in 2018 decreased by 6.5%. Under the influence of urbanization process, the maximum reduction rate of precipitation infiltration recharge is up to 38%. The study investigated the response of surface runoff and precipitation infiltration recharge to land use change, which can provide helps for water resources managers to coordinate the relationship between land use change and rational water resources planning.


Introduction
Groundwater is an extremely precious natural resource. It plays a key role in the hydrologic cycle, ecological and geological environment as well as economic and social development [1,2]. As an important part of water resources, groundwater is preferred to be used for water supply than surface water due to its characteristics of wide distribution, strong storage capacity, good water quality and easy development [3,4]. Precipitation infiltration is the main natural recharge source of groundwater. In recent years, the development of urbanization in Beijing, China, has greatly changed the underlying surface conditions. The land hardening and vadose thickening brought about by this change make the recharge amount of groundwater obtained from precipitation infiltration decrease continuously.
However, in order to meet the economic development as well as the production and living needs, the exploitation of groundwater is still increasing year by year. To ensure the sustainable development of the city, it is urgent to evaluate the land use change on precipitation infiltration coefficient, runoff and groundwater in this case.
In the application of groundwater modeling, the precipitation infiltration recharge is determined according to the total precipitation and the precipitation infiltration recharge coefficient (α) in a certain area [5,6]. The coefficient mainly depends on the lithology, structure and thickness of the vadose zone as well as the land surface slope [7][8][9]. The natural amplitude of phreatic water level, lysimeter observation, isotope analysis, equilibrium method, chloride tracer and numerical simulation are usually used to calculate α [10][11][12][13][14]. With the rapid economic and social development in recent years, a variety of human activities (such as the acceleration of urbanization process) have contributed to significant changes in land use and land cover, which directly affect the underlying surface of the basin [15], and then changes the rainfall infiltration process.
Recently, researches about the influences of the land use change are abundant, including climate change [16], hydrology and solute budgets [17], streamflow [18], and soil chemical attributes [19]. By analyzing the characteristics of climate change (including spatial distribution and temporal trend of precipitation and temperature) and land use/cover change in Yunnan Province, China during 1961-2011, Shi and Chen [16] concluded that the close relationships of the built-up land with the annual temperature and precipitation changes might be attributed to the developments of urbanization. Based on the lithology of unsaturated zone and the depth to groundwater table in different periods, Meng et al. [6] discussed the temporal and spatial evolution characteristics of precipitation infiltration recharge in the North China Plain in recent 50 years. The study showed that with the increase of water table depth, the precipitation infiltration recharge coefficient increased and reached the peak value at a particular depth (the optimal depth of different lithology is 3 to 6 m). Cherubin et al. [19] conducted the investigation on the effects of land use change on Brazilian tropical soil chemical attributes, indicating that long-term land use change from natural ecosystems to extensive pasture decreases available P, S, Ca, Mg and B contents, and land use change from pasture to sugarcane will lead to increase of nutrients levels and reduction of soil acidity. However, in the area where the groundwater recharge is dominated by the rainfall infiltration, the investigations on the impact of land use change on rainfall infiltration recharge remain limited. Therefore, it is important to study the variation of rainfall infiltration coefficient under the condition of underlying surface change.
The Soil and Water Assessment Tool (SWAT) model can comprehensively consider soil, land use, vegetation, land surface slope, meteorology, reservoir and runoff information to simulate different hydrological cycle processes. Jin et al. [20] improved SWAT model to solve the problem that SWAT could not reflect the land use change during the simulation period. The improved SWAT model can use land use data year by year, and the simulation effect is better than SWAT Model in the Heihe River basin in China where land use and cover change greatly. SWAT model uses land use data in the modeling process, which will affect hydrological processes such as evapotranspiration and groundwater recharge [20]. SWAT model has certain advantages in the calculation of groundwater recharge resources under the condition of underlying surface change. Based on the SWAT hydrological model, Mohamed et al. [21] simulated the responses of streamflow and baseflow to climate and land use change in the Upper West Branch DuPage River basin in Illinois and the Walzem River Basin in Texas, respectively. The modeling results show that changes in streamflow and baseflow in the latter watershed seem to be more affected by urbanization than by climate change. Zhang et al. [22] applied the improved SWAT model to simulate the impact of different land use changes on water balance in the North Johnston River catchment in eastern Australia, showing that urbanization will increase surface runoff. The change of underlying surface will lead to the increase of surface runoff and the decrease of precipitation infiltration recharge, causing the decrease of groundwater recharge amount. Correspondingly, the groundwater level will continue to drop when the groundwater exploitation remains unchanged or continues to increase. In this case, it may lead to surface subsidence and other geological hazards, threaten ecological security, water resources security [23][24][25]. Thus, it is necessary to study the effect of land use change on the rule of runoff, rainfall infiltration coefficient and rainfall infiltration recharge, so as to reasonably plan the relationship between land use and water resources in the context of the social and economic development.
In order to evaluate the influence of underlying surface change on the hydrological cycle, we firstly analyzed the land use change in different periods of the study area, the Su-Mi-Huai District in Beijing, China, and used SWAT model to discuss the influence of underlying surface change on runoff mechanism and rainfall infiltration coefficient. Secondly, the semi-loose coupling model was established by MODFLOW and SWAT, and the results of this model were used to compare the influence of the water balance in the study area before and after land use changes in typical years. This study could provide guidance for the rational planning of land use, the correct prediction of urban water situation and the rational planning of water resources in areas with rapid urbanization.

Study Area
The study area, located in the northeast of Beijing, China (as shown in Figure 1), is referred to as Su-Mi-Huai area in this work. The total area of the study area is approximately 1078 km 2 . The mountainous area is distributed in the north while the plain area is located in the middle and south, covering an area of about 680 km 2 . The northern mountain boundary of the study area is divided according to the location of the surface watershed. The climate of the area is semi humid and semi-arid continental monsoon climate in typical warm temperate zone, with a corresponding annual precipitation of 589 mm. The precipitation is unevenly distributed within the year and mainly concentrated from June to September, accounting for more than 82% of the annual precipitation. The interannual variation of precipitation is large, and successive dry years often occur. There are four reservoirs, an artificial canal and a river with three branches in this area. These reservoirs are under the state of impoundment all the year around (store water in the reservoir without releasing water to the downstream river). Miyun Reservoir in the northern part of the study area is the largest and located at the mountain pass of the river. The river generally flows from north to south. Chaohe River and Baihe River from Miyun Reservoir converge into Chaobai River in Hecao village, Miyun District. Huaihe River flows into Chaobai River at Niulanshan town in Shunyi district. The vegetation in Miyun Reservoir Basin is mostly composed of weeds, shrubs and forests. When the altitude is higher than 1800 m, the soil type is meadow soil; when the altitude ranges from 900 to 1800 m, the mountain brown soil is dominant; in lower elevations, the soil is mainly leached brown soil.
Water 2020, 12, x FOR PEER REVIEW 3 of 18 rainfall infiltration coefficient and rainfall infiltration recharge, so as to reasonably plan the relationship between land use and water resources in the context of the social and economic development.
In order to evaluate the influence of underlying surface change on the hydrological cycle, we firstly analyzed the land use change in different periods of the study area, the Su-Mi-Huai District in Beijing, China, and used SWAT model to discuss the influence of underlying surface change on runoff mechanism and rainfall infiltration coefficient. Secondly, the semi-loose coupling model was established by MODFLOW and SWAT, and the results of this model were used to compare the influence of the water balance in the study area before and after land use changes in typical years. This study could provide guidance for the rational planning of land use, the correct prediction of urban water situation and the rational planning of water resources in areas with rapid urbanization.

Study Area
The study area, located in the northeast of Beijing, China (as shown in Figure 1), is referred to as Su-Mi-Huai area in this work. The total area of the study area is approximately 1078 km 2 . The mountainous area is distributed in the north while the plain area is located in the middle and south, covering an area of about 680 km 2 . The northern mountain boundary of the study area is divided according to the location of the surface watershed. The climate of the area is semi humid and semiarid continental monsoon climate in typical warm temperate zone, with a corresponding annual precipitation of 589 mm. The precipitation is unevenly distributed within the year and mainly concentrated from June to September, accounting for more than 82% of the annual precipitation. The interannual variation of precipitation is large, and successive dry years often occur. There are four reservoirs, an artificial canal and a river with three branches in this area. These reservoirs are under the state of impoundment all the year around (store water in the reservoir without releasing water to the downstream river). Miyun Reservoir in the northern part of the study area is the largest and located at the mountain pass of the river. The river generally flows from north to south. Chaohe River and Baihe River from Miyun Reservoir converge into Chaobai River in Hecao village, Miyun District. Huaihe River flows into Chaobai River at Niulanshan town in Shunyi district. The vegetation in Miyun Reservoir Basin is mostly composed of weeds, shrubs and forests. When the altitude is higher than 1800 m, the soil type is meadow soil; when the altitude ranges from 900 to 1800 m, the mountain brown soil is dominant; in lower elevations, the soil is mainly leached brown soil.   This area is one of the main water sources in Beijing. The groundwater in the plain area is mainly pore water in loose rock, while in mountain area it is dominated by fissure water in bedrock. From piedmont to plain, the aquifer changes from single layer to multi-layer. The groundwater in the area is mainly supplied by precipitation infiltration. Artificial exploitation is almost the only way of groundwater discharge. In recent years, the overexploitation of groundwater resources led to the increase of the groundwater depth and the thickness of the vadose zone and the change of underlying surface conditions. At the same time, the original overflow zone disappears. The groundwater depth is greater than the evaporation limit (4 m), which means that there is no phreatic water evaporation when the water level is lower than 4 m [26].

SWAT Model
The input data needed for SWAT model are mainly divided into two categories: spatial data and attribute data are shown in Table 1. The former includes the digital elevation model (DEM), land use and soil type, while the latter includes hydrological data, meteorological data (daily rainfall, daily average wind speed, daily relative humidity, solar radiation and daily maximum and minimum temperature) [27,28]. Among them, DEM data is used to (1) generate river network, (2) determine watershed outlet and (3) divide sub-basin. The hydrological response unit (HRU) is generated by using the superimposition of the land use type map, soil type distribution map and slope data. Soil data is used to calculate and create the soil databases needed for the SWAT model.  1 The data set is provided by Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences. (http://www.gscloud.cn). 2 The data set is provided by Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn). 3 The data set is provided by Food and Agriculture Organization of the United Nations (FAO) (http://www.fao.org). 4 The data set is provided by Beijing Hydrological Center, China.
The groundwater recharge in SWAT model is calculated by the Green and Ampt formula [29]. The Soil Conservation Service (SCS) runoff curve method is used to estimate surface runoff in SWAT model. The SCS calculation formula [30] is as follows: where Q sur f is the cumulative runoff (mm), R day is the rainfall depth of one day (mm), S is retention parameters (mm). Only when R day > 0.2S, surface runoff is generated.
We processed the remote sensing image data (Table 1) of different periods to get the land use transfer matrix applied in land use analysis. The specific measures are as follows: firstly, we open the land use map of two periods in ARCGIS at the same time and fuse the two maps according to the land use field. Then, the land use maps of the two periods are superimposed and cross analyzed to obtain the intersection layer. We assign values to the attribute table of the intersection layer and  calculate its area and then export the attribute table. Finally, we open the attribute table through Excel,  and perform the Pivot Table and Pivot Chart operations to get the land use transfer matrix [31].

MODFLOW Model
MODFLOW is a three-dimensional, physically based, distributed finite-difference groundwater model for variably saturated subsurface systems [32]. In MODFLOW, if the change of water density is not considered, the three-dimensional flow of groundwater in porous media can be expressed by the following partial differential equation: where S s is water storage rate of water bearing medium (1/m), h is water level elevation of aquifer (m), K x , K y , K z is the permeability coefficient of x, y and z directions respectively (m/d), W is the source and sink term of aquifer (1/d). MODFLOW software has a main program and several relatively independent subprograms. Each subroutine has several modules, and each module completes part of the numerical simulation. This modular structure is quite conducive to the increase and decrease of data and the establishment of later coupling model. In this study, the optimized RCH subprogram package and the newly developed RAW (recharge and well) subprogram by Han et al. [33] are selected to process the source and sink data. These two subroutines are the principle of area data processing. Precipitation infiltration recharge, agricultural irrigation supply and partial mining are carried out in the form of area. In the simulation, the study area can be divided into different sub regions to deal with the related recharge and discharge. This method is helpful to use them to process the replenishment and discharge data, and to replace the water supplied to aquifer in SWAT model with the same areal distribution.

The Semi-loose Coupling
This work mainly uses three parts to complete the semi loose coupling of surface water and groundwater: (1) surface water runoff model based on SWAT2005 software, (2) groundwater flow model based on MODFLOW-2005 and (3) semi-loose coupling model of surface water and groundwater based on GIS platform for data transmission. The process of establishing SWAT-MODFLOW semi-loose coupling model is described as follows: (1) the groundwater recharge in the plain area calculated by SWAT model is used to replace the rainfall infiltration amount in MODFLOW model, (2) the optimized RCH sub-package is selected to deal with precipitation infiltration recharge and (3) it is substituted into model as the area recharge item. In addition, the MODFLOW model's lateral recharge in the piedmont area is replaced by the mountainous aquifer recharge calculated using the SWAT model, and the recharge is substituted into the model in the form of well (WELL).
The key problem in the establishment of semi-loose coupling model is that the basic calculation unit of MODFLOW model is finite difference grid (CELL), which has spatial information. The basic calculation unit of SWAT model is hydrological response unit (HRU), which is the basis of model simulation calculation [34,35]. Each HRU represents the combination of specific land use, soil type and slope in the watershed. They have vector information, but their data cannot be directly transferred to the grid (CELL). For this reason, we established an interactive interface of HRU-CELL through ARCGIS software [36]. The specific situation is shown in Figure 2. Through this interactive interface, each hydrological response unit corresponds to one or more grids, and the values on each HRU are transferred to the corresponding CELL, that is, the semi-loose coupling of SWAT-MODFLOW is realized.

Model Calibration and Validation
The SWAT model takes the day as the time step, and uses the sequence uncertainty fitting method (SUFI2) in the SWAT-CUP (Calibration and Uncertainty Programs) 2012 to identify and analyze the sensitivity of the model parameters and determine the parameter range [37]. The SUFI2 algorithm judges the sensitivity of parameters by t-test and p-factor method. Further, the parameter is more sensitive with a larger absolute value of t. The p value represents the significance of t value. The smaller the p factor, the more significant the parameter sensitivity.
The represents the ratio of the observed data variance which is interpreted by the model as follow [27,39]: The appraises a standardized difference between the simulation results and the observed results as follow [27,39]: where and are the observed and computed values, respectively, and are the averages of the observed and computed values, respectively. When the values of and are gradually close to 1, it indicates that the simulation value of the model is closer to the measured value, that is, the less variation in the error, and the better simulation effect [40]. In general, the accuracy of monthly runoff simulation can be satisfied when the values of and are above 0.6 [41].

Precipitation Frequency Analysis and Precipitation Infiltration Coefficient
The amount of precipitation controls the balance and dynamics of groundwater system in the study area and determines the quantity of groundwater recharge resources. Therefore, we should analyze the precipitation characteristics and predict the variation trend of rainfall infiltration.
In this study, the precipitation guarantee rate was firstly analyzed. According to the precipitation data of many years, the rainfall series , , …, , …, can be obtained by arranging the rainfall from large to small. The guarantee rate (cumulative frequency) of precipitation is defined as follows:

Model Calibration and Validation
The SWAT model takes the day as the time step, and uses the sequence uncertainty fitting method (SUFI2) in the SWAT-CUP (Calibration and Uncertainty Programs) 2012 to identify and analyze the sensitivity of the model parameters and determine the parameter range [37]. The SUFI2 algorithm judges the sensitivity of parameters by t-test and p-factor method. Further, the parameter is more sensitive with a larger absolute value of t. The p value represents the significance of t value. The smaller the p factor, the more significant the parameter sensitivity.
The R 2 represents the ratio of the observed data variance which is interpreted by the model as follow [27,39]: The E ns appraises a standardized difference between the simulation results and the observed results as follow [27,39]: where Q o and Q c are the observed and computed values, respectively, Q o and Q c are the averages of the observed and computed values, respectively. When the values of E ns and R 2 are gradually close to 1, it indicates that the simulation value of the model is closer to the measured value, that is, the less variation in the error, and the better simulation effect [40]. In general, the accuracy of monthly runoff simulation can be satisfied when the values of E ns and R 2 are above 0.6 [41].

Precipitation Frequency Analysis and Precipitation Infiltration Coefficient
The amount of precipitation controls the balance and dynamics of groundwater system in the study area and determines the quantity of groundwater recharge resources. Therefore, we should analyze the precipitation characteristics and predict the variation trend of rainfall infiltration.
In this study, the precipitation guarantee rate was firstly analyzed. According to the precipitation data of many years, the rainfall series R 1 , R 2 , . . . , R i , . . . , R n can be obtained by arranging the rainfall from large to small. The guarantee rate P i (cumulative frequency) of precipitation R i is defined as follows: where P i is greater than or equal to the frequency of occurrence of precipitation R i , i is the serial number of precipitations in order of magnitude, n is the number of years of statistics. Precipitation infiltration coefficient is the ratio of precipitation recharge to precipitation. This coefficient is an important hydrological parameter in the study of groundwater resources estimation and mutual transformation between atmospheric water, surface water and groundwater.
The precipitation infiltration recharge coefficient changes greatly due to the comprehensive influence of natural factors and human factors. The annual precipitation infiltration coefficient can be calculated as follows [42]: where α is precipitation infiltration coefficient, PERC is groundwater recharge (mm), PREC is precipitation (mm). The empirical value of the precipitation infiltration coefficient is determined by previous studies based on factors such as the lithology, structure and thickness of the unsaturated zone [7][8][9].

The Results of Model Calibration and Validation
SuZhuang station, a hydrological observation station in the study area, is selected to calibrate the runoff parameters based on the monthly measured runoff. The SuZhuang station is located at the downstream outlet of the basin, which can reflect the runoff from the entire watershed to the outlet. The runoff data of this station is provided by Beijing Hydrological Center, China. The collected runoff data is divided into two parts: the first part (from January 2000 to December 2012) is used for model calibration, and the second part (from January 2013 to December 2016) is used for model validation. The comparison of simulated and observed monthly streamflow during calibration and validation periods is shown in Figure 3. During the calibration period, the values of R 2 and E ns are 0.64 and 0.60, respectively. During the validation period, the values of R 2 and E ns are 0.78 and 0.70, respectively. The values of R 2 and E ns in the calibration period and the validation period are both greater than 0.6, which meets the accuracy requirements, indicating that the calculation results of SWAT model are reliable.
The results of MODFLOW model and SWAT-MODFLOW model have been presented in our previous study, which indicate that the simulation results of models are reliable [43].

Analysis on Land Use Change
In order to fully demonstrate the direction of land type transfer and supplementary sources, we used the land use transfer matrix to conduct the research. Table 2 shows the land use transfer matrix of the study area from 2000 to 2015. It can be seen that there are eight types of land use in the study area, mainly including cultivated land, woodland, rural area and urban area.
From 2000 to 2015, the area of cultivated land in the study region decreased by 22.85%. Rural residential land increased by 44.22%. The area of urban residential land has become 2.81 times the area of the same land type in 2000.
Rural residential land and urban residential land are mainly converted from cultivated land. The area of grassland, woodland and waters decreased to varying degrees, and the percentage of reduction was 44.2%, 11.5% and 19.1%, respectively. The expansion of land used for urban construction in the study area is more significant, which is mainly manifested in the increase of urban and rural residential land and the decrease of arable land.
Water 2020, 12, x FOR PEER REVIEW 8 of 18 Figure 3. Runoff fitting diagram of SuZhuang station, adapted from our previous study [43].

Analysis on Land Use Change
In order to fully demonstrate the direction of land type transfer and supplementary sources, we used the land use transfer matrix to conduct the research. Table 2 shows the land use transfer matrix of the study area from 2000 to 2015. It can be seen that there are eight types of land use in the study area, mainly including cultivated land, woodland, rural area and urban area.
From 2000 to 2015, the area of cultivated land in the study region decreased by 22.85%. Rural residential land increased by 44.22%. The area of urban residential land has become 2.81 times the area of the same land type in 2000.
Rural residential land and urban residential land are mainly converted from cultivated land. The area of grassland, woodland and waters decreased to varying degrees, and the percentage of reduction was 44.2%, 11.5% and 19.1%, respectively. The expansion of land used for urban construction in the study area is more significant, which is mainly manifested in the increase of urban and rural residential land and the decrease of arable land.   Figure 4 shows the land use types of the study area in 2000 and 2015, the zoning of the empirical value of precipitation infiltration coefficient in plain area and the rainfall infiltration coefficient zoning of the study area based on SWAT surface runoff simulation calculation using the land use map of 2015.

Spatial Variation Analysis of Precipitation Infiltration Coefficient
As the process of urbanization accelerates, land use types have changed, and the precipitation infiltration coefficient will also change accordingly. Figure 4b,d shows that the precipitation infiltration coefficient calculated by SWAT model is smaller than the empirical value of precipitation infiltration coefficient in plain area. In the mountainous area, the rainfall infiltration coefficient calculated by SWAT model is greater than 0.33. The maximum range of precipitation infiltration coefficient obtained by the two methods is in the northern plain area, but the maximum range of empirical value is 0.01-0.17 larger than the SWAT calculation result. infiltration coefficient calculated by SWAT model is smaller than the empirical value of precipitation infiltration coefficient in plain area. In the mountainous area, the rainfall infiltration coefficient calculated by SWAT model is greater than 0.33. The maximum range of precipitation infiltration coefficient obtained by the two methods is in the northern plain area, but the maximum range of empirical value is 0.01-0.17 larger than the SWAT calculation result.   Figure 4b, the areas with large empirical value of precipitation infiltration recharge coefficients are concentrated in the upper part of the alluvial fan, where the lithology is mainly gravel and coarse sand and the precipitation infiltration coefficient ranges from 0.56 to 0.65. The precipitation infiltration coefficient in the middle part of the plain area is greater than 0.33, while it is less than 0.33 in the northern mountain and the downstream boundary. In 2000, the land use types in the study area were relatively primitive, and large-scale urbanization had not started. As it can be seen from Figure 4a, the main land use type of the plain area was cultivated land, accounting for 69.24% of the plain area. Rural and urban residential areas accounted for 16.19% and 5.14% of the plain area, respectively. The woodland is mainly distributed in the northern mountainous area, accounting for 21.35% of the total area.

As shown in
According to the Table 2 and Figure 4c,d, the proportion of cultivated land, the main land use type in the region, decreased from 51.29% to 39.56% by 2015 due to the acceleration of urbanization. In the plain area of the upper reaches of the river, the precipitation infiltration coefficient of the area with more cultivated land is relatively large, ranging between 0.48 and 0.55, less than the empirical value of precipitation coefficient.
Compared with the land-use type maps in Figure 4, it can be clearly seen that rural residential  Figure 4b, the areas with large empirical value of precipitation infiltration recharge coefficients are concentrated in the upper part of the alluvial fan, where the lithology is mainly gravel and coarse sand and the precipitation infiltration coefficient ranges from 0.56 to 0.65. The precipitation infiltration coefficient in the middle part of the plain area is greater than 0.33, while it is less than 0.33 in the northern mountain and the downstream boundary. In 2000, the land use types in the study area were relatively primitive, and large-scale urbanization had not started. As it can be seen from Figure 4a, the main land use type of the plain area was cultivated land, accounting for 69.24% of the plain area. Rural and urban residential areas accounted for 16.19% and 5.14% of the plain area, respectively. The woodland is mainly distributed in the northern mountainous area, accounting for 21.35% of the total area. According to the Table 2 and Figure 4c,d, the proportion of cultivated land, the main land use type in the region, decreased from 51.29% to 39.56% by 2015 due to the acceleration of urbanization. In the plain area of the upper reaches of the river, the precipitation infiltration coefficient of the area with more cultivated land is relatively large, ranging between 0.48 and 0.55, less than the empirical value of precipitation coefficient.

As shown in
Compared with the land-use type maps in Figure 4, it can be clearly seen that rural residential area and urban residential area have increased significantly. In the areas where urban residential land increased (such as the western, northeastern and southwestern riverside areas of the study area), the precipitation infiltration coefficient decreased 0.18-0.44. The precipitation infiltration coefficient of the concentrated area of urban residential land in the study area is generally less than 0.32 due to the ground hardening caused by urbanization.  Figure 6 presents the area changes of the main land use types in the five periods as well as the corresponding annual precipitation infiltration coefficient.

Time Variation Analysis of Precipitation Infiltration Coefficient
According to the analysis of Figure 6, from 2000 to 2018, the average annual precipitation infiltration coefficient shows a decreasing trend with the change of land use, which is negatively correlated with the change of construction land area (urban area and rural area), but positively correlated with the change of cultivated land area. From 2000 to 2018, the average annual precipitation infiltration coefficient decreased by 6.5%. During the years from 2000 to 2005, the average annual precipitation coefficient decreased significantly, by 3.4%. During this period, the area of construction land increased, such as urban area, rural area and traffic land, which increased by 41.3%, 13.1% and 25.5%, respectively. The runoff coefficient is the ratio of surface runoff to rainfall. It can comprehensively reflect the influence of natural geographical factors on the relationship between precipitation and runoff. The average annual runoff and runoff coefficient under each land use scenario are calculated. The annual average value of surface runoff in the study area is increasing continuously. The average runoff coefficient in 2015 was 17.9% higher than 2000.
From 2000 to 2005, the change range of land use was the largest, and the annual increment of surface runoff was the largest, which was 5 × 10 6 m 3 . Comparing each year's surface production flow and annual average surface production flow under the four phases of land use, the changes are the same, both showing an increasing trend. The added value is the largest during the period of the fastest land use change (2000)(2001)(2002)(2003)(2004)(2005). From 2000 to 2015, when the land use type changed from cultivated land to construction land, the runoff coefficient increased with the increase of average runoff yield, and the increase rate of runoff coefficient was greater than that of surface runoff. Water 2020, 12, x FOR PEER REVIEW 11 of 18  According to the analysis of Figure 6, from 2000 to 2018, the average annual precipitation infiltration coefficient shows a decreasing trend with the change of land use, which is negatively correlated with the change of construction land area (urban area and rural area), but positively correlated with the change of cultivated land area. From 2000 to 2018, the average annual precipitation infiltration coefficient decreased by 6.5%. During the years from 2000 to 2005, the average annual precipitation coefficient decreased significantly, by 3.4%. During this period, the area of construction land increased, such as urban area, rural area and traffic land, which increased by41.3%, 13.1% and 25.5%, respectively.  According to the analysis of Figure 6, from 2000 to 2018, the average annual precipitation infiltration coefficient shows a decreasing trend with the change of land use, which is negatively correlated with the change of construction land area (urban area and rural area), but positively correlated with the change of cultivated land area. From 2000 to 2018, the average annual precipitation infiltration coefficient decreased by 6.5%. During the years from 2000 to 2005, the average annual precipitation coefficient decreased significantly, by 3.4%. During this period, the area of construction land increased, such as urban area, rural area and traffic land, which increased by41.3%, 13.1% and 25.5%, respectively. yield under the four land use types has the same change trend, and the runoff yield fluctuates with rainfall. In 2012, the surface runoff was the largest, ranging from 82 × 10 6 m 3 to 94 × 10 6 m 3 . The surface runoff in 2003 and 2009 was the smallest, ranging from 22 × 10 6 m 3 to 28 × 10 6 m 3 . According to the precipitation histogram, the year with the largest precipitation is 2012, and the year with the minimum precipitation is 2003 and 2009. The peak and valley values of runoff are the same as those of precipitation. The runoff coefficient is the ratio of surface runoff to rainfall. It can comprehensively reflect the influence of natural geographical factors on the relationship between precipitation and runoff. The average annual runoff and runoff coefficient under each land use scenario are calculated. The annual average value of surface runoff in the study area is increasing continuously. The average runoff coefficient in 2015 was 17.9% higher than 2000.
From 2000 to 2005, the change range of land use was the largest, and the annual increment of surface runoff was the largest, which was 5 × 10 6 m 3 . Comparing each year's surface production flow and annual average surface production flow under the four phases of land use, the changes are the same, both showing an increasing trend. The added value is the largest during the period of the fastest land use change (2000)(2001)(2002)(2003)(2004)(2005). From 2000 to 2015, when the land use type changed from cultivated land to construction land, the runoff coefficient increased with the increase of average runoff yield, and the increase rate of runoff coefficient was greater than that of surface runoff.

Analysis of Water Balance in Typical Years
In this study, we used the rainfall series data from 1956 to 2016 to calculate the guarantee rate of each year (Appendix A) and count the frequency of precipitation during this period (Table 3). It can be seen that the extremely dry year is 2000 (the guarantee rate is 95%) in which, the rainfall is 313.9 mm. The general dry year is 2013 (the guarantee rate is 75%), and the corresponding rainfall is 465.4 mm. The normal year is 1970 (the guarantee rate is 50%), the simulation period representative year is

Analysis of Water Balance in Typical Years
In this study, we used the rainfall series data from 1956 to 2016 to calculate the guarantee rate of each year (Appendix A) and count the frequency of precipitation during this period (Table 3). It can be seen that the extremely dry year is 2000 (the guarantee rate is 95%) in which, the rainfall is 313.9 mm. The general dry year is 2013 (the guarantee rate is 75%), and the corresponding rainfall is 465.4 mm.
The normal year is 1970 (the guarantee rate is 50%), the simulation period representative year is 2010, the rainfall is 555.7 mm. The rainy year is 2012 (the guarantee rate is 15%), the rainfall is 753.6 mm, and the statistical annual average rainfall is 588.8 mm. The water balance of representative years is selected in this section to study the impact of land use change on the overall water balance in the study area. As shown in Table 4, the phreatic aquifer water equilibria of the study area in three characteristic years before and after land use change are selected, which are calculated by MODFLOW model and SWAT-MODFLOW semi-loose coupling model [43] respectively.
In the rainy year, the calculated precipitation infiltration replenishment after considering the land use change is less than the calculation result when the land use change is not considered, and the reduction rate is 16.5%. In the normal year, under the condition of land use change, the amount of precipitation infiltration recharge decreased by 33.6%. In the dry year, under the condition of land use change, the amount of precipitation infiltration recharge decreased by 38%. In the rainy year, the calculated precipitation infiltration replenishment after considering the land use change is less than the calculation result when the land use change is not considered, and the reduction rate is 16.5%. In the normal year, under the condition of land use change, the amount of precipitation infiltration recharge decreased by 33.6%. In the dry year, under the condition of land use change, the amount of precipitation infiltration recharge decreased by 38%.
The results show that the impact of land use change on water balance is mainly reflected in the aspect of precipitation infiltration recharge [44,45]. Land hardening leads to a decrease in the amount of groundwater recharge by precipitation infiltration, which is especially obvious in the year of decreasing atmospheric precipitation, that is, the response in dry year is stronger than that in normal year, and the response in rainy year is the weakest. Although the investigation conducted in this work is based on a small-scale model, its results are similar to those obtained from large-scale studies [44][45][46]. Urbanization increases runoff and reduces groundwater recharge, posting a challenge for catchment management. The study on the change of flood peak in river urban agglomeration [47,48] and the change of regional climate [49,50] are related to the rapid urban growth. Groundwater is very important for irrigation and global food security. When the supply of groundwater is reduced, the groundwater level drop caused by irrigation cannot be replenished in time, which will affect the aquatic ecology of rivers [51]. Accordingly, it is important to reasonably plan the land use of a town even a region, and comprehensively consider its impact on runoff and groundwater infiltration.
The methodology used in this work can quantify the influence of land use change on hydrological cycle in a watershed, including precipitation infiltration coefficient, surface runoff and water balance, which is helpful to provide guidance for water resource management. It is different from the qualitative understanding obtained by remote sensing data or other traditional methods, that is, the logarithmic relationship between urbanization rate and runoff [52]. In areas strongly affected by human activities, it may not be realistic to obtain detailed hydrological parameters to establish a full coupling model of surface-water and groundwater for water resources evaluation and hydrological process research. The combination of the remote sensing data, SWAT model and MODFLOW model make it possible. The application of this comprehensive method in Su-Mi-Huai area demonstrates the simulation results are credible. Specifically, the remote sensing data of different periods were used to calculate the degree of urbanization construction and the transformation of land use types. The SWAT and MODFLOW model are semi loosely coupled to quantify the amount of water resources affected by land use change. In order to ensure the sustainable development of an area, we should not only balance the relationship between urban development and water resources allocation but also consider the effect of urbanization on hydrological cycle. Certainly, this requires the coordinated planning and joint efforts of land planning, water resources management and other departments.

Conclusions
Based on the Landsat 8 high-precision remote sensing data and the SWAT hydrological model, this study analyzes the change law of the falling water infiltration coefficient and surface runoff in the study area under the land use change scenario and analyzes the effect of land use change on the water balance of the study area using the MODFLOW model. The conclusions as follows are reached.
From 2000 to 2015, the land use change was mainly manifested in the transformation from cultivated land to urban and rural residential land in Su-Mi-Huai area in Beijing, China.
By 2015, the precipitation infiltration coefficient value in the area is smaller than its empirical value. In the area of increasing urban residential land, the precipitation infiltration coefficient decreased by 0.18-0.44. From 2000 to 2018, the average annual precipitation infiltration coefficient in the study area has gradually decreased, which is positively correlated with the change of cultivated land area and negatively correlated with the change of construction land.
Under the conditions of land use from 2000 to 2015, the annual average value of surface runoff and the runoff coefficient in the area showed an increasing trend. During the period of the fastest land use change (2000)(2001)(2002)(2003)(2004)(2005), the increase values of surface runoff and annual average surface runoff of each year are the largest.
Changes in land use mainly reduce precipitation infiltration replenishment. In the dry year, precipitation infiltration recharge has the strongest response to land use changes, with a reduction rate of 38%.
This study examined the influence of hydrology cycle due to the land use change primarily caused by urbanization. The current study provides useful implications for the rational management of water resources under the situation of land use change caused by urbanization. The sustainable development of cities must consider the impact of land use change on the hydrological cycle and find a reasonable balance between urban planning and water resources management.