Study on Hydrologic Effects of Land Use Change Using a Distributed Hydrologic Model in the Dynamic Land Use Mode

Abstract

It is reasonable to simulate the hydrologic cycle in regions with drastic land use change using a distributed hydrologic model in the dynamic land use mode (dynamic mode). A new dynamic mode is introduced into an object-oriented modularized model for basin-scale water cycle simulation (MODCYCLE), a distributed hydrologic model based on sub-watersheds, and the hydrological response unit (HRU). The new mode can linearly interpolate data for the years without land use data and consistently transfer HRU water storage between two adjacent years after a land use data update. The hydrologic cycle simulation of the Sanjiang Plain in China was carried out from 2000 to 2014 in the dynamic mode using land use maps of 2000, 2005, 2010, and 2014. Through calibration and validation, the performance of the model reached a satisfactory level. Replacing the land use data of the calibrated model using that of the year 2000, a comparison model in the static land use mode (static mode) was built (i.e., land use unchanged since 2000). The hydrologic effects of land use change were analyzed using the two models. If the land use pattern remained unchanged from 2000, despite the average annual runoff increasing by 4% and the average annual evapotranspiration decreasing by 4% in this region only, the groundwater storage of the plain areas in 2014 would increase by 4.6 bil. m³ compared to that in 2000, rather than the actual decrease of 4.7 bil. m³. The results show that the fluxes associated with groundwater are obviously more disturbed by land use change in the Sanjiang Plain. This study suggests that the dynamic mode should be used to simulate the hydrologic cycle in regions with drastic land use change, and the consistent transfer of HRU water storage may be considered in the dynamic mode.

1. Introduction

Human activities have caused a large number of natural landscapes (woodland, grassland, wetland, etc.) to be replaced by artificial landscapes (farmland, cities, roads, etc.), resulting in changes in land use/land cover (land use for short) patterns [1,2]. There have been many studies on the impacts of land use/land cover change (land use change for short) on various aspects of the hydrologic cycle. These include analyses of the effects of urbanization; opening up of wasteland for farming; and changes to agricultural land in terms of runoff, evaporation, groundwater recharge, and basic flow under various regional conditions (e.g., global scale, arid, agricultural, urban areas, and watersheds) [3,4,5,6,7,8,9,10,11]. The research methods applied may include statistical analysis [5,12,13,14,15], model simulation, and so on, in which model simulations are widely used in the study of hydrologic effects of land use change [16]. In particular, the emergence of distributed hydrologic models that require land use data provides powerful tools for such studies [17,18,19,20].

At present, the delta method is commonly used to study the hydrologic impacts of land use change using hydrologic models [21,22,23,24,25,26,27,28]. First, a hydrologic model is built using land use data before any changes (initial data). Second, the initial data in the model are replaced by the land use data after any changes (final data), and another model is built. The outputs from the two (or more) models, such as runoff, evapotranspiration, groundwater recharge, discharge, etc., are then compared to analyze the impact of land use change on the hydrologic cycle. However, the delta method is self-contradictory, as it acknowledges that the land use in a watershed is changing, but uses static land use data (or one land use map) in the long-term simulation of the historical hydrologic cycle, which is inconsistent with the reality of land use change in the watershed. In watersheds with little land use change, static land use data can be used, but there is no need to study the hydrologic effects of land use change in such scenarios. There are also contradictions in the model calibration and validation, e.g., using the measured data (such as runoff and groundwater level) from the watershed with continuous land use changes to calibrate the model constructed using one land use map from this watershed. There may be two reasons for these contradictions: one is the lack of land use data, and the other is that the model can only take one land use map as input [17].

As the number of methods of acquiring land use images increases, such as satellites, aircraft, and unmanned aerial vehicles [29,30,31], the constraint of unavailable land use data has gradually been overcome. However, distributed hydrologic models lag behind in expressing dynamic changes in land use, and most of them can only read one land use map in the whole simulation period [27,32,33,34]. In the hydrologic simulation of regions or watersheds where land use changes continuously and severely, such models should be used with caution. Some model developers and scholars have noticed the above shortcomings of hydrologic models and attempted to improve existing hydrologic models by developing new modules and auxiliary tools to enable the land use data to be updated while the model is running [17,35,36,37,38,39,40,41].

Land use data can be used as a kind of time series data, like meteorological data and water use data, and input into the distributed hydrologic model as one of the driving factors, as long as there are sufficient periods (years) of land use maps. Even if there are some years in the simulation period without land use data, methods such as linear interpolation [41] can be used to supplement the data, facilitating a consecutive input data series at the annual scale. However, the key to realizing the dynamic update of land use data when the model is running is the update mechanism (algorithm) of the model. As input data, weather and water use data have fixed formats, which do not increase the burden of the model algorithm, even if they are updated by day or even by hour. On the contrary, the spatial distribution of land use changes with time, which means that the environment of the hydrologic cycle changes, and then causes a change in the hydrologic cycle [42], which is more complex than that caused directly by weather and water use. Therefore, the hydrologic model should not only adapt to these changes in the algorithm, but also adapt to the needs of the algorithm in the data pre-processing step.

At present, the most widely used dynamic update mechanism of land use data in the hydrologic model is based on the hydrologic response unit (HRU), which represents homogeneous areas with a unique land use, soil, and slope class combination [17,36,37,43,44,45]. This mechanism takes the sub-watershed as a unit and changes the fractional area of each HRU in the sub-watershed (the sum of the fractional areas in a sub-watershed is always equal to 1) at the user’s predetermined time node to update the types and areas of the land use. However, this is just the first step for the model in completing the update of land use data, and the next step, which is more important, is water storage transfer; that is, the water storage in HRUs (including snow cover, canopy interception, surface ponding, and soil water content) at the end of the last time interval before land use data update is transferred to the start of the first time interval after the update.

The essence of water storage transfer is to initialize the water storage of HRUs after the land use data update. Such HRUs can be called current HRUs, but they may be the same as or different from the HRUs before the update (so-called previous HRUs). However, this moment is not been the beginning of the model simulation period, and the initial water storage of current HRUs should not be assigned at will (the model should not be preheated after each land use data update) but should be inherited from the simulated value of previous HRUs. At present, the general method is to directly transfer the water storage of previous HRUs to the current ones with the same or similar attribute (land use, soil, and slope class combination) after the update (the so-called direct transfer of water storage in this paper), but there may be two problems in this method: one is that new HRUs may appear after the update, and corresponding ones cannot be found in previous HRUs from the same sub-watershed [37]; the other is if there are corresponding ones in previous HRUs, the total water storage (dimension: L3) of the sub-watershed may be changed after the transfer, because the areas of the HRUs have changed despite their water storage (dimension: L) remaining at the previous value, which leads to unreasonable fluctuations of the simulation value and affects the whole simulation effect of the model. To solve the above problems, we introduce a new dynamic update mechanism of land use data to help our hydrologic model achieve the interpolation of years with no land use data and the dynamic update of land use data when the model is running. The new mechanism can reasonably assign initial water storage (or transfer water storage from previous HRUs) to the current HRUs and maintain the total water storage of the sub-watershed before and after land use update, which is called consistent transfer of HRU water storage in this study. The consistent transfer can avoid the sudden change in simulated values caused by the change in the fraction area of the HRUs and ensure the consistency of the hydrologic cycle simulation.

The new mechanism of the dynamic land use update has been realized in the object-oriented modularized model for basin-scale water cycle simulation (MODCYCLE), a hydrologic model based on sub-watersheds and HRUs [46] and closely coupled with the groundwater numerical model [47]. A significant advantage of the model is that it can simulate the groundwater level and flow field in the plain area, overcoming the limitations of groundwater lumped mode in the semidistributed hydrologic model, but it does not reduce the efficiency in simulating the hydrologic cycle of a large region. Previously, the model had been successfully applied in many different watersheds and regions [47,48,49,50,51,52,53], but it had not been verified in the scene of drastic land use change in the dynamic land use mode. In the current study, the improved model was applied to the Sanjiang Plain of China, where land use has been changing significantly [54,55]. It is obviously not advisable to use static land use data to simulate the hydrologic cycle in this region. The objectives of this paper were to (1) establish the hydrologic model of the Sanjiang Plain in the dynamic land use mode and test its performance, (2) establish the hydrologic model of the study area in the static land use mode and compare the hydrologic differences between the two models to clarify the effects of land use change, and (3) compare the differences between the direct and consistent transfer of HRU water storage in the dynamic land use mode to illustrate the necessity of the consistent transfer.

2. Materials

2.1. Study Area

The Sanjiang Plain is located in the northeast of China, where the Heilong River meets two of its tributaries, the Songhua River and the Wusuli River, and it is composed of several watersheds that mainly belong to the Songhua River basin and the Wusuli River basin. Its area is 105,700 km², of which the plain area accounts for 61.2%. It is part of the Northeast Plain, which is the largest in China. In terms of the hydrologic cycle, this region is relatively closed to other regions of the Northeast Plain. As the Sanjiang Plain has similar climatic characteristics, geological conditions, natural landscape, and agricultural management systems, it is often considered as a separate region for management and research purposes (Figure 1).

The Sanjiang Plain is part of the black soil zone of Northeast China, one of the four major black soil zones in the world [56], and it is suitable for a variety of field crops. Large-scale farm exploitation began in the 1950s, since when the cultivated land area has increased from 8200 km² to about 48,077 km² at present (data from the Heilongjiang Statistical Yearbook). Crops have gradually changed from various types to rice and maize. Apart from rice and a small amount of cash crops, almost all crops in the Sanjiang Plain are rainfed. According to the Heilongjiang Statistical Yearbook, paddy fields have been the land use type with the largest increase in area since 2000, accounting for nearly half of the farmland area by 2014 in the Sanjiang Plain.

2.2. Data

As a comprehensive hydrologic model, MODCYCLE requires a digital elevation model (DEM) and information about land use, soil, the river network, meteorology, agricultural management, water conservancy projects, water supply and utilization, and hydrogeological parameters. This section briefly introduces the basic spatial data and land use data.

2.2.1. Basic Spatial Data

The sub-watersheds and main channels are the basic spatial data required for hydrologic simulations in MODCYCLE. These are delineated and generated using a DEM (90 m × 90 m) and a digital real river network (1:250,000). In addition, to incorporate the groundwater numerical simulation, the sub-watersheds need to be further divided into mountain sub-watersheds and plain sub-watersheds according to the topography [47,52]. In our study, the Sanjiang Plain was divided into 1705 sub-watersheds, each of which had a main channel (Figure S1).

The numerical method was used to simulate groundwater in the Sanjiang Low Plain and the Xingkai Lake Plain. Other valley plains are small and scattered, so the lumped method based on sub-watersheds was used. Vertically, the aquifer in the plain area was divided into a shallow aquifer and a deep aquifer with variable thickness; horizontally, each aquifer was cut into 200 rows and 216 columns of uniformly spaced cells (2 × 2 km²) according to the finite difference method. That is, each aquifer has 43,200 cells, with 12,935 effective cells in the plain area (Figure S2).

2.2.2. Land Use Data

The dynamic land use mode in MODCYCLE requires multiple land use maps, so we collected four vector maps interpreted from remote sensing images of land use in 2000, 2005, 2010, and 2014 (1:100,000; data from the Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Figure S3). The main land use types in the Sanjiang Plain include farmland, woodland, wetland, grassland, and residential areas. The proportion of each type is shown in Figure 2. From 2000 to 2014, the area of farmland increased, in which the area proportion of paddy field in the whole region increased rapidly from 9.0% to 22.8%. The areas of woodland and wetland showed different downward trends, while the proportion of grassland and residential areas remained small with some fluctuations. Woodland is mainly distributed in hilly areas, which are difficult to reclaim. Wetland and grassland are mainly distributed in plain areas and tend to be surrounded by farmland.

3. Methods

3.1. Structure and Framework of MODCYCLE

MODCYCLE is an object-oriented modularized model for basin-scale hydrologic cycle simulations. The model is coded in C++ and modularized using object-oriented programming. It can simulate the hydrologic cycle on a long-term and large basin/region scale. MODCYCLE offers good practicality, distributed computing, a strong conceptual and physical mechanism, and a clear and complete cycle path (Figure 3) and accurately reflects the effect of human activities on the hydrologic cycle. There are many distinctive characteristics of the model, including the hierarchical water-balance-checking mechanism, the database platform, support for parallel computing, and tight coupling with a groundwater numerical model [46,47,50,52].

The model first divides the watershed/region (one region may involve multiple watersheds) into different sub-watersheds using a digital elevation model (DEM). Each sub-watershed has a main channel, and the hydraulic relationship between sub-watersheds is constructed from the tree topology of the main channel. According to the superimposed distribution of land use, soil, and agricultural management, each sub-watershed is further divided into several HRUs. These are not the same as a piece of land. In fact, they are the aggregate of unique land use, soil type, and agricultural management patterns in the sub-watershed. HRUs do not have spatial location attributes in the sub-watershed of the model; thus, they are relatively independent and have no hydraulic connection with each other. The channel system in each sub-watershed consists of a sub-channel and a main channel. The sub-channels collect the water yield from the HRUs, which is partly transported to the natural waters (depressions, wetlands, and ponds) in the sub-watershed and partly to the main channel. The main channels and reservoirs together constitute the river network system of the model.

The groundwater system is divided into shallow and deep layers. The model only simulates the interaction between shallow and deep groundwater within a sub-watershed. Groundwater in different sub-watersheds is independent, and there is no lateral groundwater exchange between sub-watersheds. This groundwater system may be reasonable in mountainous areas, where surface watersheds are usually identical to the groundwater watersheds (except in karst mountainous areas) [57,58]. However, the water parting of sub-watersheds in plain areas is not obvious, and the groundwater aquifers in different sub-watersheds are interlinked; thus, the lateral movement of groundwater cannot be ignored. To handle this issue, MODCYCLE is closely coupled with a groundwater numerical model, controlled by the partial differential equation of three-dimensional groundwater flow based on the finite difference method, on the basis of the original hydrologic cycle simulation framework, and separates the large and complete plain areas from watersheds/regions to carry out groundwater numerical simulations [47,52].

3.2. Dynamic Land Use Mode

After adding the dynamic update mechanism of land use data, MODCYCLE has two land use data input modes: static and dynamic. These are used to simulate the hydrologic cycle of a watershed with little/no land use change and severe land use change, respectively. The static mode uses one land use map to construct an HRU dataset, whereas the dynamic mode uses multiple land use maps to construct the same number of HRU datasets. In the dynamic mode, the land use data may not be available every year during the simulation period; thus, HRU interpolation is needed for years without land use data. The area of an HRU in year i can be interpolated using the following formula:

$$A_i=\{\begin{array}{l}{A}_a, if no land use data before the yea{r}_i\\ A_b+(A_a-A_b)\cdot \frac{yea{r}_i-yea{r}_b}{yea{r}_a-yea{r}_b}\\ A_b, if no land use data after the yea{r}_i\end{array},$$

(1)

where A_i is the area of an HRU in year i (year_i) without land use data (m²) and A_a and A_b are the areas of the same HRUs in the same sub-watershed in the most recent year with land use data after (year_a) and before (year_b) year i, respectively (m²). If there are no land use data before or after year i, i.e., the years at either end of the data series, the HRU areas should not be subjected to linear extension in order to avoid negative values.

The sub-watersheds remain unchanged once delineated in the simulated watershed, so the model takes the sub-watershed as the basic unit for the transfer of the HRU water storage. It should be emphasized that the model requires the division of the soil profiles of different soil types into the same number of layers and thickness to ensure that the soil water content of different soil types can be transferred via one-to-one correspondence in soil layers. Assuming that the model divides all soil types into 1–L layers from top to bottom, the components that need to be transferred for HRU water storage include the soil water content in layers 1–L, surface ponding, snow cover, canopy interception, and other water storage components. The calculation process for each sub-watershed and water storage component is the same. An example is described below.

The HRUs and their areas in a sub-watershed in year i may change in year i + 1, but one of the following two scenarios will always occur: an HRU will have a reduced area (reduced HRU) or the same/increased area (non-reduced HRU). The former includes HRUs in year i that disappear in year i + 1, and the latter includes HRUs that do not exist in year i but appear in year i + 1 (new HRUs) and those that have no change in area. The model classifies HRUs in the two years according to the above classification method (Figure 4).

It is assumed that there are M reduced HRUs and N non-reduced HRUs, where the sum of the reduced area of the former equals the sum of the increased area of the latter. The water storage component of reduced HRUs in year i + 1 is

$$W_{i+1}(m)=W_i(m),$$

(2)

where W_i+1(m) and W_i(m) are the water storage components of the m-th reduced HRU (m = 1, 2… M) at the beginning of the simulation in year i + 1 and at the end of the simulation in year i, respectively (mm). The latter is the final simulation result in the i-th year of the model. In other words, the water storage components of the reduced HRUs are directly inherited from the previous HRUs.

Because the HRUs lose their spatial attributes in the process of delineation and interpolation, it is impossible to determine the transformation relationship between different land use types. Therefore, the model first accumulates the water storage components of the reduced part of the reduced HRUs and then distributes the accumulated water storage component to the increased part of the non-reduced HRUs according to the proportion in area. The accumulated water storage component of the reduced part of the reduced HRUs is calculated by the following formula:

$$V_i={\displaystyle \sum _{m=1}^M{W}_i(m)\cdot [A_i(m)-A_{i+1}(m)]},$$

(3)

where V_i is the accumulated water storage component of the reduced part of the reduced HRUs (10⁻³ m³); A_i(m) and A_i+1(m) are the areas of the m-th reduced HRU in years i and i + 1, respectively (m²); and W_i(m) is the water storage component of the m-th reduced HRU (m = 1, 2…M) at the end of the simulation in year i (mm). The accumulated water storage component is distributed to the increased part of the non-reduced HRUs as follows:

$$V_i(n)=V_i\cdot \frac{A_{i+1}(n)-A_i(n)}{{\displaystyle \sum _{n=1}^N[A_{i+1}(n)-A_i(n)]}},$$

(4)

where V_i(n) is the water storage component distributed to the increased part of the n-th non-reduced HRU (n = 1, 2…N) (10⁻³ m³), V_i is the accumulated water storage component of the reduced part of the reduced HRUs (10⁻³ m³), and A_i+1(n) and A_i(n) are the areas of the n-th non-reduced HRU in years i + 1 and i, respectively (m²).

The water storage component of a non-reduced HRU in year i + 1 may be divided into two parts: one is from the accumulated water storage component of the reduced part of the reduced HRUs, i.e., V_i(n), and the other is from the water storage component of the corresponding HRU in year i, i.e., W_i(n). Both are distributed to the non-reduced HRU in year i + 1 using an area-weighted average, which can be calculated as follows:

$$W_{i+1}(n)=\frac{V_i(n)+W_i(n)\cdot A_i(n)}{A_{i+1}(n)},$$

(5)

where W_i+1(n) and W_i(n) are the water storage components of the n-th non-reduced HRU (n = 1, 2…N) at the beginning of the simulation in year i + 1 and at the end of the simulation in the year i, respectively (mm); V_i(n) is the water storage component distributed to the increased part of the n-th non-reduced HRU from the accumulated water storage component of the reduced part of the reduced HRUs (10⁻³ m³); and A_i+1(n) and A_i(n) are the areas of the n-th non-reduced HRU in years i + 1 and i, respectively (m²).

In summary, to ensure that the total water storage of each sub-watershed remains consistent before and after land use changes, all kinds of HRUs of a sub-watershed in both years are divided into two categories to initialize the water storage component of each HRU at the beginning of the year i + 1. Combining Equations (2) and (5), we obtain a formula for calculating the water storage component of each HRU at the beginning of year i + 1 for a sub-watershed:

$$W_{i+1}(HRU)=\{\begin{array}{l}{W}_{i+1}(m)\\ W_{i+1}(n)\end{array}$$

(6)

The model applies the above consistency calculation process (also called initialization process or consistent transfer) to other water storage components and other sub-watersheds, thus maintaining the water storage across the whole watershed as the land use is updated between adjacent years. The typical case in Figure 4 is represented according to the above process of water storage transfer, as shown in Figure 5.

3.3. Model Performance Evaluation

The hydrologic model is capable of simulating the actual hydrologic cycle and predicting the future trend of the hydrologic cycle by the processes of parameter calibration and model validation, but statistical performance measures and corresponding performance evaluation criteria are needed in these processes [59]; that is, the effect that the model achieves in the simulation of measured data indicates that the model has the above abilities. At present, many statistical performance measures are applied, including Nash–Sutcliffe efficiency (NSE) [60], th ecoefficient of determination (R²), root-mean-square deviation (RMSD), the index of agreement (d) [61], percentage bias (PBIAS) [62], and so on. These statistical performance measures cooperate with graphic performance measures to support model performance evaluation. In the current study, two statistical performance measures, the NSE and R², were used to evaluate the model performance in terms of streamflow simulation. Other performances of the model, groundwater level, and irrigation amount simulation were evaluated by R². The expressions for computing the NES and R² are shown in Equations (7) and (8). Performance evaluation criteria referenced the corresponding quantitative thresholds recommended by Moriasi [59].

$$NSE=1-\frac{{\displaystyle \sum _{i=1}^n{(Q_{m,i}-Q_{s,i})}^2}}{{\displaystyle \sum _{i=1}^n{(Q_{m,i}-{\overline{Q}}_m)}^2}},$$

(7)

$$R^2=\frac{{\left[{\displaystyle \sum _{i=1}^n(O_i-\overline{O})(P_i-\overline{P})}\right]}^2}{{\displaystyle \sum _{i=1}^n{(O_i-\overline{O})}^2\cdot {\displaystyle \sum _{i=1}^n{(P_i-\overline{P})}^2}}},$$

(8)

where Q_m,i is the i-th observed value of streamflow, Q_m is the mean value of the observed streamflow, and Q_s,i is the i-th simulated value of streamflow. O_i is the i-th observed value; O is the mean value of the observed values; P_i is the i-th simulated value; and P is the mean value of the simulated values for streamflow, groundwater level, and irrigation amount.

3.4. Analysis of Land Use Change and Its Hydrological Effect

By comparing the simulated output of MODCYCLE in the dynamic and static modes of land use, this study shows the applicability and rationality of the dynamic mode in the area with dramatic land use change and the impact of the continuous land use change on the main hydrologic cycle components in the study area. In the hydrologic effect analysis section, the differences between the two land use change patterns (one for the dynamic mode and the other for the static mode) are first analyzed to support the subsequent comparative analysis of hydrologic components.

Taking sub-watersheds as individual units, the comparison between the two land use change patterns is carried out by examining the relative difference in land use (LUD_R) between the average area of a certain type of land use in the dynamic mode (LU_D) and the corresponding type of land use area in the static mode (LU_S):

$$\mathrm{LU}{\mathrm{D}}_{\mathrm{R}}=(\mathrm{L}{\mathrm{U}}_{\mathrm{D}}-\mathrm{L}{\mathrm{U}}_{\mathrm{S}})/A(sub),$$

(9)

where A(sub) is the area of the sub-watershed for the considered LU (km²).

In the hydrologic effect analysis, this study mainly analyzes the differences in runoff, evapotranspiration, and groundwater between the two modes. Among them, comparisons between the two modes in terms of the average annual runoff yield and evapotranspiration at the sub-watershed scale are calculated by the following formulas:

$$\mathrm{RF}{\mathrm{D}}_{\mathrm{R}}=(\mathrm{R}{\mathrm{F}}_{\mathrm{D}}-\mathrm{R}{\mathrm{F}}_{\mathrm{S}})/\mathrm{R}{\mathrm{F}}_{\mathrm{S}},$$

(10)

$$\mathrm{E}\mathrm{T}{\mathrm{D}}_{\mathrm{R}}=(\mathrm{E}{\mathrm{T}}_{\mathrm{D}}-\mathrm{E}{\mathrm{T}}_{\mathrm{S}})/\mathrm{E}{\mathrm{T}}_{\mathrm{S}},$$

(11)

where RFD_R is the relative difference in the average annual runoff yield between the dynamic and the static land use mode (-), and RF_D and RF_S are the average annual runoff yield of the sub-watershed in the dynamic mode and the static mode, respectively (mm). ETD_R is the relative difference in the average annual evapotranspiration between the dynamic and the static land use mode (-), and ET_D and ET_S are the average annual evapotranspiration of the sub-watershed in the dynamic mode and the static mode, respectively (mm).

4. Results

Using the available data, we first calibrated and validated the hydrologic model with the dynamic land use mode. Holding other data and parameters constant, the model was then converted to the static land use mode, and the land use data were replaced by the land use map of 2000. The differences in runoff, evapotranspiration, and groundwater in the Sanjiang Plain in the dynamic mode and the static mode were analyzed to clarify the influence of land use changes on the hydrologic cycle. Finally, we analyzed the difference in hydrologic cycle simulation between the consistent transfer and the direct transfer of water storage in the process of land use update.

4.1. Model Calibration and Validation

Many aspects of comprehensive distributed hydrologic models can be calibrated [63,64,65]. However, the measured data used for calibration and validation are often limited; measured streamflow data are the most prevalent because of their accessibility. Coupling MODCYCLE with the groundwater numerical model enables the simulation of groundwater levels, so the measured groundwater depth can be used for verification. In addition, a large number of rice irrigation observation experiments were carried out at several experimental rice-planting stations in the Sanjiang Plain, and some representative measured rice irrigation data were obtained. These can be used to verify the irrigation simulation for paddy fields in the model.

The water yield of the Sanjiang Plain only accounts for a small part of the flow through the main stream of Songhua River. Therefore, it is not appropriate to use the measured flow data from the main stream to calibrate and validate the model. We used the measured flow data from six hydrologic stations on the representative tributaries of the Songhua River and the Wusuli River and the measured groundwater depth data from 21 representative observation wells in the plain area to calibrate (2000–2008) and validate (2009–2014) the model (measured data from hydrologic bureaus of the Songhua River basin management organization and Heilongjiang).

The model can simulate the historical hydrologic cycle through parameter calibration (model calibration). Model calibration approaches include automatic calibration, manual calibration, and hybrid calibration [66,67]. This study adopted the manual calibration approach based on the experience of modelers.

Uncertainty is common in current hydrologic models [68], and parameter uncertainty presents significant challenges to model calibration [69]. According to the calibration experience of previous studies, several key parameters sensitive to runoff and groundwater response were selected through trial calculations, and the ranges of these parameters were continuously narrowed to determine the optimal parameter combination for optimizing the simulations. The Sanjiang Plain is composed of several watersheds, each of which has different underlying surface characteristics, so the sensitivity ranking and final calibration values of the model parameters are different. Considering the differences in each watershed, the key parameters of the model for the whole Sanjiang Plain were sequenced according to their sensitivity, and the maximum and minimum values were determined by model calibration (

Table 1. Key parameters used for hydrologic cycle and groundwater simulations at the representative hydrologic stations and groundwater observation wells.

Variable	Description (Units)	Recommended Range	Actual Min./Max. Values Used	Hydrologic Process
Hydrologic cycle simulation
MXSURPOND	Maximum depth of surface ponding (mm)	0.0–150.0	1–100	Runoff
ALPHA_BF	Baseflow alpha factor (days)	0.0–1.0	0.05–0.08	Groundwater
ESCO	Soil evaporation compensation factor (-)	0.01–1.0	0.9–0.92	Evaporation
SURLAG	Surface runoff lag coefficient (-)	1.0–24.0	5.0–5.0	Runoff
SOL_AWC	Available water capacity of the soil layer (-)	0.0–1.0	0.01–0.25	Soil water
SOL_K	Saturated hydraulic conductivity (mm/h)	0.0–25.0	0.018–25	Soil water
GWDMN	Threshold water level in shallow aquifer for baseflow (m)	0.0–5.0	2.5–2.5	Groundwater
Groundwater numerical simulation
HY	Hydraulic conductivity (m/d)		0.25–7.5	Groundwater
SC1	Primary storage coefficient 1 (-)		0.008–0.175	Groundwater
SC2	Primary storage coefficient 2 (-)		0.004–0.175	Groundwater

).

The observed streamflow and simulated streamflow at the Baoquanling hydrologic station on the Wutong River, a tributary of the Songhua River, and the Baoqing hydrologic station on the Naoli River, a tributary of the Wusuli River, are shown in Figure 6. This figure also presents the observed and simulated groundwater depths at two representative observation wells located in the Sanjiang Low Plain and the Xingkaihu Plain. Model calibration and validation were evaluated using the coefficient of determination (R²) and the Nash–Sutcliffe efficiency (NSE). As shown in the charts, the streamflow simulations at the two stations were satisfactory according to the criteria recommended by Moriasi [59]. The groundwater depth fitting between the observations and simulation results were only evaluated using R², and the performance was relatively poor compared with that of the streamflow simulations. Fitting the simulated average groundwater depth from a 2 km grid cell to the measured values at the point scale of the observation wells effectively conceals many small-scale factors affecting the groundwater level of a single well, such as nearby groundwater exploitation and the uneven rainfall distribution. These factors inevitably affect the accuracy of the simulations. However, the simulated period and range of groundwater depth fluctuations generally reflect the actual situation.

We also compared the simulated data with the measured data at all hydrologic stations and groundwater observation wells. Figure 7a shows the results for the monthly flow. According to the criteria suggested by Moriasi [59], both the NSE and R² were satisfactory. Figure 7b shows that the overall simulation effect of groundwater depth at all observation wells at the end of the month from 2002 to 2014 was much better than that of a single well, and the R² value reached 0.90. The fluctuation in groundwater levels varied at different observation wells. The scatter point group determined by the simulated values and the observed values from a single well is distributed in different locations in Figure 7b. Although the R² value of a single scatter point group is relatively small, all scatter point groups were distributed on the diagonal 1:1 line, which greatly improved the R² value. In general, the simulation of groundwater level in the whole plain area achieved satisfactory results.

Paddy fields represent the most prominent type of land use change in the Sanjiang Plain. Because surface water and pumped groundwater are used for large-scale irrigation, the increase in the area of paddy fields may have greatly changed the hydrologic cycle in the Sanjiang Plain. Therefore, the demanded irrigation of paddy fields should be accurately simulated, which is completed in MODCYCLE using a method similar to that of Tsuchiya [70]. We collected the annual irrigation amount entering paddy fields from eight experimental irrigation stations in the Sanjiang Low Plain from 2008 to 2012. A comparison between the measured annual irrigation amount and the simulated value at each experimental station shows that the R² value reaches 0.61, a satisfactory simulation effect (Figure 8).

In summary, through the model calibration and validation in terms of streamflow and groundwater level, the model can simulate the historical hydrologic cycle and predict the future hydrologic cycle. The verification of paddy field irrigation further demonstrates the model’s ability to simulate artificial interventions in the hydrologic cycle.

4.2. Hydrologic Effects of Land Use Change

The model was used to simulate the hydrologic cycle of the Sanjiang Plain from 2000 to 2014 in two modes: the continuous yearly change in land use (actual or close to actual dynamic mode) and no change in land use since 2000 (assumed static mode). The effects of land use changes on the hydrologic cycle were analyzed in terms of runoff, evapotranspiration, and groundwater.

4.2.1. Comparison of the Two Land Use Change Patterns

To enhance the comparison of the hydrologic effects, the differences between the two land use change patterns must be clarified. Although actual land use data were only obtained for four out of 15 years, and those for other years were interpolated, the dynamic land use mode provides a better reflection of the real land use trend in the Sanjiang Plain than the static land use mode. According to Equation (9) in the Methods section, if LUD_R is greater than 0, the area of this type of land use increased (or decreased in some years and increased in other years but increased overall) from 2000 to 2014; otherwise, the area of this land use type decreased. Changes in the four main land use types in the Sanjiang Plain are shown in Figure 9.

In the plain areas, farmland was the main land use type, and the conversion from dry land to paddy fields was the most remarkable; in the middle and northeast of the Sanjiang Low Plain, dry land areas were observed to increase. There was not a lot of woodland in the plain areas, and that which existed was mainly shrubland, but there was still a decreasing trend in the northeastern part of the Sanjiang Low Plain. Wetland along rivers tended to increase, whereas that in low-lying areas tended to decrease. Hilly areas were dominated by woodland, with relatively small changes in land use. Some woodland in valley plain areas was converted to dry land.

4.2.2. Runoff Effect

Since 2000, land use in the Sanjiang Plain has changed continuously, especially in terms of dry land and paddy fields. If land use had not changed since then, the average annual runoff yield and the annual runoff process in the Sanjiang Plain would inevitably be different from the actual values. The distribution of RFD_R at the sub-watershed scale is shown in Figure 10a. The runoff in sub-watersheds dominated by farmland is generally reduced; that is, the increase in farmland and the transformation from dry land to paddy fields enhanced the utilization of precipitation in farmland and reduced the occurrence of runoff. Of course, there may also be waterlogging hazards in some areas, with the need for drainage leading to an increase in runoff, such as the Xingkai Lake Plain. Runoff did not change much in hilly areas and sometimes decreased because part of the woodland was reclaimed for farmland in valley plain areas. The annual process of runoff also diverged between the two land use modes, as shown in Figure 10b. Before 2005, the two modes nearly had similar annual runoff, but after 2005, the difference increased. Overall, in the dynamic mode, the farmland area gradually increased, with paddy field areas increasing sharply, which resulted in the corresponding reduction in runoff.

4.2.3. Evapotranspiration Effect

Evapotranspiration (ET) is an important flux in the hydrologic cycle and typically changes with the intensification of human activities. In this section, we compare the temporal and spatial differences in the ET (mainly including canopy interception evaporation, snow sublimation, ponding evaporation, vegetation transpiration, and soil evaporation on HRUs) in the two modes to illustrate the impact of land use change on regional evapotranspiration. The spatial distribution of ETD_R in the Sanjiang Plain at the sub-watershed scale is shown in Figure 11a. The map shows that the ET increased significantly in plain areas, where paddy fields expanded dramatically in the dynamic land use mode. In hilly areas, the ET also differed between the two modes, mainly in valley plains where the land use changed most significantly. Figure 11b compares the annual ET process in the whole region. The annual ET deviation of the two modes gradually enlarged, with the ET clearly higher in the dynamic mode. The main reason for this difference is the change from dry farmland to paddy fields.

4.2.4. Groundwater Effect

Groundwater is an important source of fresh water in the Sanjiang Plain, and it also plays an important role in maintaining ecological stability. The large-scale cultivation of rice means that the overexploitation of groundwater is serious in areas with insufficient surface water and poor groundwater recharge conditions. If the pattern of land use in 2000 was maintained, the balance of groundwater and groundwater levels would not have caused the current problems.

We compared the main groundwater recharge, discharge, cumulative storage change, and groundwater level in the dynamic and static land use modes (Figure 12). Driven by the increase in paddy field area, the total groundwater recharge, including rainfall and irrigation percolation recharge and surface water seepage recharge, and groundwater exploitation (mainly paddy field irrigation) increased significantly in the dynamic mode. However, there may be a critical value of groundwater exploitation below which groundwater evaporation and baseflow do not change significantly, such as the exploitation before 2006. If the groundwater was continuously exploited on a large scale for many years, even an increase in precipitation would not be sufficient to compensate for the deficit in groundwater. In this situation, the groundwater evaporation and baseflow would decrease, and it would be difficult for groundwater storage to recover to the level of 2000. As can be seen from Figure 12f, in the eastern part of the Sanjiang Low Plain, the groundwater level in the dynamic mode was at least 5 m lower than that in the static mode at the end of 2014, forming a large-scale groundwater level depression. This region is the main rice-producing area, in which there is insufficient surface water, and the surface is covered with 3–17 m of silty clay, which is not conducive to groundwater recharge. Once the groundwater falls, it is difficult to be recovered. Rice cultivation in the Xingkai Lake Plain also developed very fast, but this was mainly irrigated by surface water and had better groundwater recharge conditions, resulting in relatively stable groundwater levels.

4.3. Difference between Consistent and Direct Transfer of HRU Water Storage

HRU water storage transfer is one of the important steps of dynamic land use update in MODCYCLE. Its function is to initialize water storage of the new year’s HRUs after land use update and maintain the consistency of the total water storage of each sub-watershed with that of the previous year. This step can prevent unacceptable sudden changes in hydrologic cycle simulation and is also called the consistent transfer of HRU water storage, which was introduced in the Methods section. At present, the direct transfer of HRU water storage is most widely used; that is, after the land use update, the water storage of the current HRU comes directly from the previous HRU with the same or similar attributes. However, there may be two defects in the direct transfer: one is that the current HRUs may not find the corresponding ones in the previous HRUs of the same sub-watershed, and the other is that the total water storage of the sub-watersheds may change greatly after the transfer. This section comparatively analyzes the differences between consistent and direct transfer in water storage simulation.

Due to the first defect of the direct transfer, we modified the code of MODCYCLE for water storage transfer so that the HRUs with corresponding ones can directly inherit the water storage of the previous HRUs (dimension: L), while the HRUs without corresponding ones adopt the initial value at the beginning of the model simulation period. As a feasible direct transfer, it is compared with the consistent transfer.

The time of land use update is the first day of each year in MODCYCLE, during the winter in the Sanjiang Plain and when the water storage of HRUs is mainly snow and soil moisture. The daily variation of snow water storage in the watershed scale is significant, which means that it cannot give prominence to the sudden change in water storage caused by direct transfer. The soil moisture is generally stable in winter because it is in the freezing state, so the soil water content is taken as the object of comparison between the two water storage transfers (Figure 13).

The soil water content in the Sanjiang Plain is relatively stable in winter, but there are differences among different land use types, even though the differences are much smaller than those in other seasons [71]. In the hydrologic cycle simulation based on the dynamic land use mode, if the transfer of HRU water storage cannot be handled properly after land use update, it leads to sudden changes in soil water content simulation, which will be more pronounced in the areas with dramatic land use change. Figure 13 shows that the direct transfer results in varying degrees of sudden changes in soil water content on January 1 of each year in the whole region, the plain areas, and the areas with rapidly expanding paddy fields. In the whole region and the plain areas, the sudden changes before 2006 are less significant than those after 2006 (note: the increase in the paddy field area was relatively small before 2006, as shown in Figure 2), but in the areas with rapidly expanding paddy fields, the relatively large sudden changes appear from the second year, which indicates that the increase in the paddy field area is one of the important factors aggravating the sudden changes. The differences in the soil water content between the two transfer methods are similar in the whole region and the plain areas (green lines in Figure 13a,b); only the total soil water content of the former is much greater than that of the latter, which indicates that the drastic change in land use in the plain areas is the main reason for the sudden changes in soil water content under the direct transfer.

The direct transfer of HRU water storage not only causes sudden changes but also results in durative variation in soil moisture simulation after land use update, which can be seen from the differences between the two transfer methods. In the three spatial scales, the simulated values of soil water content in the first year are the same under the two transfer methods, but after the first land use update, the differences appear, the maximum value of which is up to 1.5% of the total soil water content simulated using the consistent transfer. Although the magnitude of the differences caused by the direct transfer is not enough to affect the hydrologic cycle simulation of the whole large-scale region, the simulation error is magnified in the small-scale watershed with severe land use change. For example, in sub-watershed No. 385, with an area of 273.5 km², the dry land area increased by 66.2 km² and the wetland area decreased by 84.2 km² from 2000 to 2014. The sudden changes caused by the direct transfer and their durative impact make the maximum difference in soil water content reach 8.7% of the total volume, and the average value also reaches 3.9%.

In terms of the water balance in the whole region, the direct transfer mainly affects soil evaporation, vegetation transpiration, surface water outflow, and groundwater storage by changing HRU water storage (

Table 2. Comparison of the average annual water balance from 2002 to 2014 under the consistent transfer and the direct transfer in the Sanjiang Plain (10⁹ m³).

Hydrologic Fluxes and Storage Changes	Consistent	Direct
Recharges
Precipitation on non-surface water	56.31	56.31
Precipitation on surface water	0.12	0.12
Surface water inflow a	46.73	46.73
Groundwater boundary inflow	0.00	0.00
Outside water diversion for irrigation b	0.75	0.75
Outside water diversion for other water supply c	0.75	0.75
Total recharge	104.65	104.65
Discharges
Interception evaporation	2.70	2.70
Snow sublimation	0.18	0.18
Ponding evaporation	7.34	7.34
Soil evaporation	23.55	23.58
Vegetation transpiration	11.44	11.43
Surface water evaporation	1.76	1.76
Surface water outflow	54.81	54.72
Groundwater boundary outflow	0.00	0.00
Evaporation of irrigation system	1.42	1.42
Other water consumption d	1.21	1.21
Total discharge	104.41	104.34
Storage changes
Storage change in soil water e	0.65	0.77
Storage change in surface water	−0.10	−0.10
Storage change in groundwater	−0.31	−0.36
Total storage change	0.24	0.31

Note: ^a Surface water inflow refers to the flow from the upper reaches of the Songhua River into the Sanjiang Plain. ^b Outside water diversion involves artificial transfers from outside of the Sanjiang Plain. ^c Other water supply includes industrial, domestic, and eco-environmental water supply. ^d Other water consumption includes industrial, domestic, and eco-environmental water consumption. ^e Soil water (or HRU water) includes canopy interception, snow cover, ponding, and soil profile moisture.

). The direct transfer has little impact on the water recharge and discharge of the whole region, but it has a significant impact on water storage change, especially on the storage change in soil water, which is 18.5% larger than that under the consistent transfer.

The above comparative analysis shows that if the dynamic land use mode is used to simulate the hydrologic cycle in the area with severe land use change, special attention should be paid to the water storage transfer after the land use update. Otherwise, unexpected errors in hydrologic cycle simulation can occur, and the uncertainty of the simulation may increase. The consistent transfer of HRU water storage is the recommended method in this study, and it can be used in the dynamic land use mode of the hydrologic cycle model with a similar structure.

5. Discussion

At present, the distributed and semidistributed hydrologic models are mainly based on static land use data input [27,33,34,72]. There are doubts regarding the use of runoff and groundwater-level data monitored under the environment of dramatic land use change to calibrate and validate the model under the static land use mode. In this case, the function of the land use dynamic update is useful for a hydrologic model. The hydrologic model of the Sanjiang Plain from 2000 to 2014 was constructed under the dynamic mode of MODCYCLE using land use maps of 2000, 2005, 2010, and 2014. The model’s simulation on the actual hydrologic process reaches a satisfactory level through calibration and validation, which indicates that the calibrated parameters approximately reflect the actual land surface, vadose zone, and hydrogeological characteristics of the Sanjiang Plain. There are many concerns about whether the dynamic mode is better than the static mode for a hydrologic model. Therefore, we used the land use data of 2000 to replace the land use data of the above calibrated model, and we established a static land use hydrologic model under the same parameters and other data. By comparing the runoff, evapotranspiration, and groundwater data from the two models, we found that there are significant differences, especially in groundwater. By the end of 2014, the difference in groundwater storage in the plain areas between the two models reaches 9.3 bil. m³, and the difference in groundwater level in the northeastern part of the Sanjiang Low Plain is more than 5 m. Thus, the hydrologic process simulated under the static mode is far from the actual hydrologic process simulated under the dynamic mode. This indicates that the static mode is not as good as the dynamic mode in simulating the hydrologic cycle of the Sanjiang Plain; at least, it is difficult to do so only using the land use map of 2000.

Because of the lack of a dynamic land use mode, many hydrologic models have to be used to study the hydrologic effects of land use change by the delta method [21,22,23,24,25,26,27,28]. In addition to abovementioned doubt regarding model calibration, the delta method mostly uses the annual mean values to study the hydrologic effects of land use change [38]. Because the models are all hypothetical, and none of them really reflects the hydrologic process under the environment of continuous land use change, the comparison of time process values between different models is practically meaningless. The hydrologic model based on the dynamic land use mode solves the above problems. The model under the dynamic mode of the Sanjiang Plain is representative of the actual hydrologic cycle rather than a hypothetical model. Thus, we assume that, since 2000, the local government has taken active land use control measures to avoid the disorderly growth of paddy fields, strengthen the protection of wetlands and woodlands, and generally maintain the land use pattern at that time. What would be the state of the hydrologic cycle of the Sanjiang Plain over the 15 years? Would it be better than the hydrologic cycle problem under the actual land use change scenario, and how much would it improve? These questions are the reasons for the establishment of the hydrologic model using the land use of 2000. Comparisons of annual mean and time process values between the dynamic mode and the static mode are more practical than those between hypothetical models of the delta method, and the former can support the formulation of local land use management measures. If the future land use patterns can be predicted, we can use these land use data to establish models under the dynamic or the static mode and compare them with the calibrated model to study the hydrologic cycle effects under the future land use policy, which is also more practical than the model comparison assumed by the delta method.

The advantages of MODCYCLE compared to other hydrologic models in the dynamic mode are also concerns of researchers. At present, we have not found a distributed hydrologic model with the ability of dynamic land use update. In semidistributed hydrologic models, the Soil and Water Assessment Tool (SWAT) has this ability [43], but it does not consider the consistent transfer of HRU water storage in the sub-watershed before and after land use update, nor can it interpolate the years without land use data. The comparison of hydrologic cycle simulation between MODCYCLE and SWAT under the dynamic mode will be our follow-up work, but we compared the differences between the consistent transfer and the direct transfer of HRU water storage under the dynamic mode of MODCYCLE. In the simulation of soil water content, the direct transfer leads to sudden changes and subsequent continuous variations. The smaller the area of the watershed and the more severe its land use change, the more significant the sudden changes and subsequent impacts. In summary, the direct transfer leads to unexpected errors and increases uncertainty in the hydrologic simulation in areas with dramatic land use change, so it is suggested to adopt the consistent transfer of HRU water storage in the dynamic mode.

Although the model has made some progress in expressing continuous land use change, it still has some limitations. The response of soil structure and land surface properties to land use change has a time lag, which often needs several years of gradual development [73]. Land use change also affects land–atmosphere interactions and feedback, as well as causing climate variation [74,75], which, in turn, changes the land cover pattern. These have yet not been considered in the dynamic mode and may assist in the further improvement of the model. The model uses the linear interpolation method to supplement the years without data. However, land use changes non-linearly. In the absence of land use data, this method is only an expedient measure. Fortunately, five-year land use updates are sufficient in most cases [17]. The uncertainty of land use data is one of the sources of model uncertainty. Land use dynamic update may further increase the uncertainty of the model, so the evaluation of model uncertainty may be our next research subject.

6. Conclusions

A new dynamic land use mode is introduced into MODCYCLE, a hydrologic model based on sub-watersheds and HRUs. The new mode can linearly interpolate data for the years without land use data and transfer HRU water storage between two adjacent years consistently. The consistent transfer of HRU water storage represents the key function of the dynamic land use mode, which can maintain the stability of the total water storage of the sub-watershed before and after the land use data update and avoid the sudden changes in water storage caused by the direct transfer in the simulation. The hydrologic cycle simulation of the Sanjiang Plain from 2000 to 2014 was carried out under the dynamic land use mode using the land use maps of 2000, 2005, 2010, and 2014. Through calibration and validation, the performance of the model reached a satisfactory level. As a comparison for analyzing the impact of land use change on the hydrologic cycle, another model under the static land use mode was established by replacing the land use maps of the calibrated model using that of 2000.

The differences in simulated runoff, evapotranspiration, and groundwater between the dynamic land use mode and the static land use mode were compared. The results show that the increase in farmland and the transformation of dry land to paddy fields enhanced the utilization of precipitation and reduced the average annual yield of runoff by 4% compared to the static mode. Under the dynamic land use mode, evapotranspiration increased significantly in the low-plain areas with paddy field expansion; some differences also appeared between the two modes in hilly areas, but these were relatively small. On the whole, maintaining the land use pattern of 2000 in the Sanjiang Plain would reduce average annual evapotranspiration by 4% compared with the actual land use change. Driven by the increase in paddy field areas, total groundwater recharge and exploitation increased remarkably, and phreatic evaporation and baseflow decreased gradually in the late simulation period; the groundwater storage in 2014 decreased by 4.7 bil. m³ compared with that in 2000, which would increase by 4.6 bil. m³ in the static mode. Additionally, a large-scale groundwater level depression formed in the northeastern part of the Sanjiang Low Plain, and this was consistent with the actual depression distribution. From another point of view, the comparative analysis results of the two modes show that the dynamic mode is helpful in improving the simulation accuracy of the hydrologic model in the areas with dramatic land use change. In terms of analyzing the hydrologic response of land use changes, hydrologic comparisons between the actual dynamic land use mode and the hypothetical static land use mode are more practical and reasonable than those between the two hypothetical land use modes.