Simulation of Crop Water Demand and Consumption Considering Irrigation Effects Based on Coupled Hydrology‐Crop Growth Model

Hydrological models are widely used in regional drought assessments. Usually, the spatiotemporal development of droughts is evaluated only via drought indexes based on simulated soil moisture. However, the physiological status of crops is key to accurately assessing the effect of drought, and hydrological models either do not consider crop physiological status or oversimplify it. In addition, the effects of reservoirs, irrigation, and crop rotation are generally not considered. To overcome the oversimplification of crop and irrigation modules in hydrological models, the variable infiltration capacity (VIC) model was coupled with the Environmental Policy Integrated Climate model. Irrigation and reservoir modules were introduced into the coupled hydrology‐crop growth (CHC) model. The CHC model was built for the watershed upstream of Heilin hydrological station on the Qingkou River on the border of Jiangsu and Shandong Provinces, China. The model performance was compared with simulation results of the Soil and Water Assessment Tool and VIC, which showed that the CHC model was more accurate. The temporal variation of water consumption and demand of three crop rotations, namely spring peanuts, summer peanuts‐winter wheat, and summer corn‐winter wheat, in the study area was analyzed. Drought area data obtained from Junan County yearbooks were consistent with the simulated water deficit and drought index considering the effect of irrigation. Crop irrigation increased soil water content and evapotranspiration and decreased crop water deficit. Therefore, crop irrigation was considered in hydrological simulations. The CHC model reliably simulated crop water consumption, soil moisture, and drought index, and the model simulated data could be used as basic data for water resources management and drought mitigation measures.


Introduction
Hydrological models have been widely used in regional drought assessments . Generally, the soil moisture simulated from these models is used to calculate the drought index for spatial and temporal drought assessments. However, drought cannot be accurately estimated if soil moisture alone is considered, because the soil moisture in agricultural land is greatly affected by the crop type and growth phase. Therefore, crop growth dynamics should be considered to attain accurate agricultural drought assessments (Hou et al., 2007). Nevertheless, commonly used hydrological models, such as the variable infiltration capacity (VIC) model, oversimplify crop growth simulation. For example, in the VIC model, the leaf area index (LAI) value remains constant throughout a month, and that of the same month each year does not vary with changes in the external environment (Siad et al., 2019;Tsarouchi et al., 2014). The root length also does not change during the modeling period. These oversimplifications affect the accuracy of transpiration and soil water content simulations. Consequently, to improve the performance of simulation of crop water demand and consumption, hydrological models can be integrated with crop models (Casanova & Judge, 2008;Y. Li et al., 2014;Maruyama & Kuwagata, 2010;Van den Hoof et al., 2011). In this study, crop water consumption was defined as the amount of water consumed by a crop (Supit et al., 2010). Crop water demand refers to the water required by the crop under suitable growth conditions, and water deficit is the difference between crop water demand and consumption (Kersebaum et al., 2016;Supit et al., 2010).
Crop water consumption is often estimated based on crop evapotranspiration (Shao et al., 2019), of which solar radiation is the main driver, and can be calculated based on the theory of energy balance and water vapor diffusion. Penman's equation (Penman, 1948) is widely used to calculate crop water consumption and has frequently been revised, and the Penman-Monteith equation was proposed (Monteith, 1965). However, because the Penman-Monteith model ignores soil evaporation (Gardiol et al., 2003), it cannot be used for sparse vegetation and crop evapotranspiration. The Food and Agriculture Organization of the United Nations (FAO) published a revised evapotranspiration equation (Allen et al., 1998;Guo, 1997) and the revised equation considers the effect of soil evaporation. Compared with global lysimeter observations, the Penman-Monteith equation is a universal model for accurately and reliably calculating transpiration volume and is widely applied (Anzhi, 2001). The dual-source model (Shuttleworth-Wallace model;Shuttleworth & Wallace, 1985) was also developed based on the Penman-Monteith model. This model regarded canopy and soil as two independent and interactive sources of water vapor, and a sparse vegetation evapotranspiration model was developed. However, this model only calculates the total transpiration instead of the water consumption in different soil layers. Moreover, this method was developed for a single point and cannot be applied in regional simulations (Tong et al., 2004). When regional crop water consumption is calculated, the crop water consumption at various single points needs to be interpolated (Kersebaum et al., 2016), which increases the errors caused by the spatial heterogeneity of meteorological factors and crop water consumption (Tong et al., 2004).
With the development of new measurement technologies, remote sensing is widely used to estimate crop water consumption (Olivera Rodriguez et al., 2021). Many researchers (Dangwal et al., 2015;Hoffmann et al., 2016;Jackson et al., 1981;Y. J. Wang et al., 2014) have calculated crop water consumption based on remote sensing data and estimation equations, including the crop coefficient, energy balance, water balance (Vosselen et al., 2005), and water vapor diffusion (Tong et al., 2004) methods. The crop coefficient method uses the Penman equation to calculate crop water demand and the remote sensing vegetation index to calculate the soil moisture index. Water demand is multiplied by the moisture index to obtain actual crop water consumption (Zhang, 2014). In the energy balance method, the distribution of crop water consumption is calculated based on satellite remote sensing data and energy balance equations (H. Li et al., 2008). This method considers the moisture status of vegetation and involves physical mechanisms. However, these methods ignore the characteristics of different crops at different growth stages (McDaniel, 2015;Yang et al., 2010;Zhang, 2014), and some are only suitable for areas with high vegetation cover. Remote sensing data also have the disadvantage of time discontinuity, and it is difficult to simulate the dynamic process of crop water consumption.
To overcome the shortcomings of temporal discontinuities in remote sensing data, some researchers have used crop growth models to simulate crop water consumption (Camargo & Kemanian, 2016;Z. Wang et al., 2013). In some models, such as the Agricultural Production Systems Simulator crop model (Gaydon 10.1029/2020MS002360 3 of 28 et al., 2017, Cropping Systems Simulation model (Stöckle et al., 2003), Soil Water Atmosphere Plant model (Mokhtari et al., 2018), and World Food Studies model (Cheng et al., 2016), empirical parameters need to be adjusted by incorporating experimental data, and it is difficult to apply them to large-scale areas.
The Soil and Water Assessment Tool (SWAT) (Arnold et al., 2012) includes a simplified Environmental Policy Integrated Climate (EPIC) crop model, in which potential transpiration is calculated using the Penman-Monteith equation (Jensen & Allen, 2016), and evapotranspiration is divided into soil evaporation and crop transpiration based on the aboveground biomass. Evapotranspiration is further divided into different soil layers according to soil depth. When the soil water content is less than the available soil water content multiplied by 0.25, crop water consumption is suppressed. However, water stress occurs at 0.7 times the available soil water content in some crops (Allen et al., 1998). The threshold should be revised based on crop species and transpiration. Therefore, the value of 0.25 in the EPIC model is not reasonable (Luo et al., 2008). McNally et al. (2015) assumed that evapotranspiration was affected by water stress based on the ratio of soil moisture to critical soil water content but did not consider the effects of crop type, potential evapotranspiration, and wilting moisture content. The Simulation Dual Kc (SIMDualKc) model (Rosa et al., 2016) uses the double-crop coefficient method to calculate crop evapotranspiration and divides crop evapotranspiration into crop transpiration and soil evaporation, where the former is the potential crop transpiration multiplied by the water stress coefficient. This coefficient was proposed by Allen et al. (1998) and is more reasonable, as it considers the effects of crop type, transpiration, and wilting moisture content.
The integration of crop and hydrological models to improve crop irrigation estimation has been widely evaluated at field or local scales (Balkovic et al., 2013;Faivre et al., 2004). For example, Y. Li et al. (2014) coupled the World Food Studies and HYDRUS-1D models to simulate wheat production in arid regions, but the application area was small, and these authors did not specify crop irrigation simulation. Nicola et al. (2005) coupled the dynamic vegetation and land surface models to simulate soil moisture and LAI, but they used point-scale data to verify the simulation results, and the rationality of area drought analysis could not easily be explained. Additionally, these models (Y. Li et al., 2014;Nicola et al., 2005) cannot simulate macroscopic hydrological processes and do not consider hydrological effects on crop growth and irrigation. The Community Land Model 5.0, developed by Lawrence et al. (2019), considers hydrological effects on crop growth and irrigation. However, the crop module was simplified due to the limited crop types; the water stress effect did not consider the joint influence of soil water, evapotranspiration, and crop types; the crop areas of each grid cell were divided into irrigated and rainfed fractions according to the data set of areas equipped for irrigation (Portmann et al., 2010), which did not reflect the actual situation (see Section 2.4); CLM5.0 did not consider the effect of reservoir and pond on irrigation. Fatichi et al. (2012) applied an ecohydrological model to a region smaller than 1 km 2 , which did not take into account the spatial heterogeneity of soil and meteorological data (Sawada et al., 2014).
Regional studies have been conducted, for example, Sawada et al. (2014) simulated large-scale agricultural drought based on the coupled SiB2 crop model (Sellers et al., 1996) and the geomorphology-based hydrological model (Dawen et al., 2004). Zhang et al. (2019) coupled a simple hydrological model with a crop model to simulate crop yields for agricultural drought assessment. However, the resolution of crop model parameters was low, for example, the LAI had a 16-km resolution (Sawada et al., 2014). The crop model of Zhang et al. (2019) was uniformly applied to the whole city. In addition, these authors neglected the effects of reservoirs and did not consider crop rotation. Most previous studies did not consider differences in crop rotation type in a single grid. The resolution of coupled crop and hydrological models, which could possibly be used for regional or local studies, is still limited (Sawada et al., 2014;Tsarouchi et al., 2014). Furthermore, existing studies have several limitations and unconsidered factors, for example, sparse irrigation data (Vosselen et al., 2005), and the effects of reservoirs, irrigation, and crop rotation (J. Li et al., 2016;S. Liu et al., 2018;Osborne et al., 2007).
To overcome these limitations, a coupled hydrology-crop growth (CHC) model was developed to simulate crop water demand and consumption in the Qingkou River basin in eastern China. Crop rotation was considered for each grid, and irrigation and reservoir modules were employed to incorporate the impact of irrigation on crop growth. To improve the simulation of crop water consumption, the crop water stress simulation method recommended by the FAO was introduced to improve the transpiration module in the crop model. The objectives of this study were (a) to validate the CHC model performance; and (b) to analyze the validity of temporal variations in crop water consumption; and (c) to assess the impacts of irrigation on crop water consumption, hydrological processes, and soil water content at the watershed level.

Integration of Hydrological and Crop Models
Crop growth is important for the simulation of crop water demand and consumption. Commonly used hydrological models, such as the VIC model, oversimplify the crop growth simulation. However, such simplification affects the accuracy of simulation of transpiration and soil water content. To address this issue, a hydrological model (VIC) was coupled with a crop model (EPIC; Figure 1; Zhang et al., 2021). The VIC model transferred the simulated soil moisture to the EPIC model at the beginning of the modeling period. The EPIC model calculated the LAI, root length, and evapotranspiration in this period and transferred the evapotranspiration to the VIC model. The latter then completed the water budget for the modeling period based on the transferred evapotranspiration to obtain the soil moisture content at the beginning of the next period. This simulation was repeated. Because errors existed in the EPIC model's simulation of evapotranspiration (Camargo & Kemanian, 2016;Z. Wang et al., 2013), the transpiration module was improved using the soil water stress simulation method recommended by the FAO (Allen et al., 1998) to improve the accuracy of simulation of crop water demand and consumption. When the CHC model was run, meteorological data and other parameters were fed into the VIC model ( Figure 1), which was initially run to drive the crop growth module to calculate the crop irrigation water demand. The gridded routing model was simultaneously run based on the VIC output to drive the reservoir module to simulate the reservoir water budget. The irrigation module calculated the actual irrigation water, based on the reservoir water volume and irrigation water demand. Subsequently, the crop growth model continued to simulate crop growth based on the actual irrigation water volume and output the daily crop water consumption and demand in the improved evapotranspiration module. Finally, a long series of hydrology and crop growth simulations were completed. The hydrological and crop models served as the basis of irrigation and reservoir modules (Sections 2.2-2.5).
The main goal of this study was to simulate crop water requirements and consumption and irrigation effects, the key point of this study is the simulation of crop growth, and the simulation of other vegetation growth would not affect the simulation of crop growth. Therefore, the VIC model was used to simulate the hydrologic process of the naturally vegetated area.

VIC Model
The VIC model (Liang et al., 1994) is a large-scale, semi-distributed hydrological model. In this study, the variable infiltration curve in the VIC model was used to obtain the infiltration amount, and the Arno model was used to obtain the base flow. The soil water content was calculated by combining it with the water balance equation. There are reported hydrologic parameters of VIC in China (Wu et al., 2011). To consider the heterogeneity effects of land conditions, the VIC model subdivides each merged grid's land cover into an arbitrary number of "tiles," each corresponding to the fraction of the cell covered by the particular land cover. The hydrologic process was simulated in each land cover type of the merged grid and the evapotranspiration, soil moisture, and runoff of the merged grid are the weighted average based on the area ratio of each land cover type. However, the VIC model does not involve plant growth simulation. The seasonal variation of leaf area is generalized with a series of discrete values to simulate monthly transpiration, and the leaf area in the same month does not change from year to year. The evapotranspiration of agricultural land includes the net evaporation from canopy interception and transpiration of crops but omits soil evaporation. Therefore, when crop LAI is low, the simulated evapotranspiration is lower than the actual evapotranspiration. To increase the accuracy of simulated evapotranspiration, the crop model needs to be coupled with the VIC model. The crop model calculates soil evaporation and transpiration from the crop canopy. Because the crop model simulates crop growth, the LAI can be updated daily, whereas the VIC retains a constant monthly LAI. This can improve the accuracy of simulation of crop land evapotranspiration in the VIC model.

Crop Growth Model
The simplified EPIC model was developed using the SWAT model. The EPIC model uses a unified approach to simulate more than 100 crop types and is widely used. Different crop types have different parameters, such as planting date, root growth, and LAI. Therefore, different crop growth stages are simulated with different results. Rich data and documents are available for this model, which requires little data and can easily be applied to large-scale simulations. The EPIC model uses temperature, soil water, and nutrient data to simulate potential evapotranspiration, maximum crop transpiration, actual crop transpiration, actual soil evaporation, plant height, LAI, biomass, root length, and yield. The simulation method is described by Neitsch et al. (2011). This model also simulates the effect of water stress on crop growth, quantified by the ratio of actual to potential plant transpiration (Neitsch et al., 2011): where wstrs is the water stress for a given day; E t is the maximum plant transpiration on a given day (mm; Equation 2), which is the crop water demand, and contributes to potential evapotranspiration; E t,act is the actual transpiration on a given day (mm), which is the crop water consumption; w actualup is the total plant water uptake for a given day (mm); λ is the latent heat of vaporization (MJ kg −1 ); Δ is the slope of the saturation vapor pressure-temperature curve, de/dT (kPa °C −1 ); H net is the net radiation (MJ m −2 d −1 ); G is the heat flux density to the ground (MJ m −2 d −1 ); γ is the psychrometric constant (kPa °C −1 ); K 1 is a dimension coefficient; ρ air is the air density (kg m −3 ); P is the atmospheric pressure (kPa); o z E e is the saturation vapor pressure of air at height z (kPa); e z is the water vapor pressure of air at height z (kPa); r a is the diffusion resistance of the air layer (s m −1 ); r c is the plant canopy resistance (s m −1 ; Equation 3); g l is the maximum conductance of a single leaf (m s −1 ); LAI is the LAI of the canopy (Equation 4); LAI max is the maximum LAI of the canopy; fr LAImx is the fraction of the plant's maximum LAI corresponding to a given fraction of potential heat units for the plant (Equation 5); fr PHU is the fraction of potential heat units accumulated for the plant on a given day in the initial growth period; and l 1 and l 2 are the shape coefficients. In the calculation of total plant water uptake, as the water content of soil decreases, the soil water becomes increasingly difficult for the plant to extract, and the crop water uptake is modified (Neitsch et al., 2011): where  , up ly E w is the adjusted amount of water absorbed by crops from the ly th soil layer (mm);  , up ly E w is the amount of water absorbed by crops from the ly th soil layer before correction (mm); SW ly is the soil water content of the ly th soil layer (mm); AWC ly is the available water capacity of the ly th soil layer (mm); WP ly is the water content of soil layer ly at the wilting point (mm); w actualup,ly is the actual water uptake for the layer (mm); and n is the total number of soil layers. However, this model cannot fully simulate the effect of water stress in the study area, because 0.25AWC ly is lower than the water content at the wilting point, and the soil water cannot reach 0.25AWC ly . Uniyal et al. (2019) and Marek et al. (2016) indicated that errors in plant stress functions lead to higher evapotranspiration. As a result, the water stress effect cannot be considered in the simulation of crop growth (Luo et al., 2008;Marek et al., 2017).
The transpiration module of the SIMDualKc (Rosa et al., 2016) could improve the evapotranspiration module of the EPIC model (X. Liu et al., 2015). Equation 6 was replaced by Equation 9 to consider the effect of water stress on crop water consumption: , , up ly up ly s ly w w K where K s,ly is the water stress coefficient of soil layer ly calculated using Equation 10: where TAW and RAW are the total and readily available soil water (mm), respectively; D r is the water content of root zone depletion (mm); and p is the depletion fraction at the initiation of stress: where p normal is applied for ET c = 5 mm day −1 ; ET c is the crop evapotranspiration; and p normal is defined by Allen et al. (1998). To consider the effects of crop rotation, crop rotation simulations were introduced into the crop growth model. More than one crop could be planted in some grids.

Irrigation Module
Each merged grid was divided into irrigation and non-irrigation areas (Haddeland et al., 2006) according to the irrigation ratio data from yearbooks (Cheng, 1996(Cheng, -2016. If the soil water content was below some threshold, crops experienced water stress, and crop transpiration decreased accordingly. If there was sufficient water for irrigation, water was applied to the soil via irrigation (Haddeland et al., 2006) until the soil water content reached field capacity. The difference between field capacity and actual soil water content is the amount of I r in the irrigation water required by the crop. In practical production, irrigation water is obtained from rivers, groundwater, and reservoirs. However, the study area is close to the coast, and to prevent seawater from intruding, the local government prohibits residents from extracting groundwater for irrigation. Therefore, the simulated irrigation water comes only from rivers and reservoirs. When this water was less than that required for irrigation, crops were insufficiently irrigated. To calculate the water in the reservoirs, the reservoir module was coupled with the gridded routing module. The actual irrigation water was then simulated based on Equation 12: where I a is the actual irrigation water (mm); I r is the irrigation water demand (mm); W is the water stored in the water body (mm); and e i is the irrigation efficiency, which considers the effect of seepage loss on irrigation. The seepage loss increases the soil moisture of canals or ditches. Since the canal or ditch area ratio in the merged grid is small and does not substantially influence the average soil moisture of a merged grid, the irrigation module does not consider the effect of increased soil moisture in canals or ditches. In this study, the irrigation efficiency e i was 0.49 according to Cao et al. (2016).
When irrigation was simulated, the proportion of irrigated areas of each merged grid needed to be determined first. Several data sets of irrigated areas exist, including the MIRCA 2000 data set (Portmann et al., 2010). This data set specifies irrigable grids at an 8 × 8 km 2 scale, which does not reflect the actual situation. Because this grid is large, ensuring that each point in the irrigable grid belonging to the irrigable area can be difficult (van Dijk et al., 2018). Siebert et al. (2013) published a global irrigation map, but the irrigation area ratio was higher than the crop area ratio in many grids of Jiangsu Province. Thenkabail et al. (2009) and Siddiqui et al. (2016) also published irrigated area maps, but the study area comprised 100% irrigated area, which was incorrect according to the local statistical yearbook (Cheng, 1996(Cheng, -2016. Crops are typically irrigated from the closest irrigation source (Xu, 2014). Therefore, the irrigation module assumed that the closer the crops are to rivers and reservoirs, the more easily they are irrigated. Hence, the locations of rivers and reservoirs needed to be determined first. Based on digital elevation model (DEM) data and ArcGIS software, the cumulative area of the convergence of each point was calculated. When the accumulated area of confluence exceeded some threshold at one point, the point belonged to a river or reservoir. This threshold was determined based on the remote sensing images. Then, the distance l between each crop grid and its nearest river grid was calculated. The distance l of all crop grids was collected and arranged from small to large. Subsequently, the area ratio of crop grids within the distances was calculated based on Equation 13: where m is the number of crop grids whose distance to the nearest water area was less than l; and n is the total number of crop grids ( Figure 2). Figure 2 shows that the closer crops were to the river, the greater the distribution density of the farmland was. This indicated that the assumption of irrigation from the nearest river was rational, because the cropland density near the river was higher than that far from the river. Finally, when P was approximately equal to the effective irrigated area ratio in the yearbook (Cheng, 1996(Cheng, -2016, the corresponding distance l was used to determine the irrigated areas. In these areas in the simulation, irrigation water was obtained from rivers or reservoirs through ditches or canals. According to the effective irrigation area ratio of 61.7% in the study area (Song, 2004), the longest distance between crops and the nearest rivers was 0.77 km (Figure 2). Based on this distance, the general distribution of irrigable farmland in the region was determined.

Reservoir and Routing Modules
Reservoir water volume is an important source of irrigation water in the irrigation module. To determine the amount of irrigation water from a reservoir, a reservoir module was built to simulate the storage of reservoirs (Sun, 2016). Numerous reservoir module algorithms have been designed (Haddeland et al., 2006;Hanasaki et al., 2006Hanasaki et al., , 2018, and the following water balance equation was used to build the reservoir module (Neitsch et al., 2011): where V is the volume of water stored in the reservoir at the end of the modeling period (m 3 ); V stored is the volume of water stored in the reservoir at the beginning of the modeling period (m 3 ); V flowin is the volume of water entering the reservoir during the modeling period (m 3 ); V flowout is the volume of water flowing out of the reservoir during the modeling period (m 3 ); V pcp is the precipitation volume during the modeling period (m 3 ); V evap is the evaporation from the reservoir during the modeling period (m 3 ); and V seep is the water lost from the reservoir by seepage during the modeling period (m 3 ). Because most reservoirs in the study area were small, V was taken as the reservoir capacity when the calculated reservoir water volume was greater than the reservoir capacity.
The reservoir inflow V flowin was determined using the flexible large-scale hydrological routing model (G. . In this model, the division method of the grid and sub-watershed merged grids was introduced in the river network division module, which extracted more realistic river network data. A high-resolution DEM and sub-basin kinematic wave routing were combined to incorporate the sub-basin damping effect, which provided the flexible large-scale hydrological routing model with a more reasonable structure and refined characteristics. The diffusion wave analysis method was used for river routing.

Study Area
The study area (182 km 2 ) is located in the upper watershed of the Qingkou River above the Heilin hydrological station in northern Jiangsu Province, China (118°44´17″-118°55´10″E, 34°59´37″-35°9´37″N; Figure 3). The watershed has a semi-humid continental climate, with abundant sunshine and humid rainy summers. It is characterized by a warm and temperate monsoon season and has cold and dry winters. The annual mean precipitation and temperature are 952 mm and 14.2°C, respectively. The planting system yields one to two harvests per year. The major crops are peanuts, rice, corn, and wheat, with peanuts covering the largest area.

Data
The meteorological data (including daily rainfall and air temperature data from 1976 to 2016) used for the Heilin Basin were obtained from Heilin, Qingshuijian, and Zhubian hydrological stations. The discharge data (including daily flow data in the flood season from 1978 to 2016) were obtained from Heilin hydrological station at the outlet of the study area. The soil water content data of wheat cropland were acquired from Doushan hydrological station in Shandong Province. Soil moisture gravimetric observations are made by applying the oven-drying method every 10 days (i.e., on the 1, 11, and 21 days of each month at 00:00 UTC) from 2006 to 2016. The MODIS evapotranspiration products (Vermote, 2015) comprised the evapotranspiration volume of the global 500 × 500 m 2 grid every 8 days, and the average evapotranspiration volume of the study area in 2001-2016 was calculated. Water surface evaporation observational data were obtained from Qingshuijian hydrological station. No published crop spatial distribution data exist, and the data were obtained using a decision tree model based on Landsat8 data (Figure 4). The crop rotations included wheat-peanut, wheat-corn, and spring peanut. The DEM data (Farr & Kobrick, 2000) were obtained from the Shuttle Radar Topography Mission and had a resolution of 30 m. The water area data were acquired from the 30-m resolution water area data product developed by the European Commission Joint Research Center (Pekel et al., 2016). The grid vegetation type data were derived from the global 1 km land cover type database developed by the University of Maryland. The LAI of vegetation other than crops was derived from the Land Data Assimilation System (Pisek & Chen, 2007;Wiegand et al., 1979), and the LAI of crops was simulated using the EPIC model. The soil parameters of the VIC model were calculated based on the global soil database and soil parameter equations (Reynolds et al., 2000;Saxton et al., 1986). The parameters of different crop types were obtained from the crop parameter library of the SWAT model. The hydrological model parameters were calibrated using the measured discharge data.

Calibration and Validation of CHC Model
The CHC model was run with a 5-km resolution and daily time step. A trial-and-error method and auto-optimization procedures for calibration were used based on Mbonimpa et al. (2015) and Rosenbrock (1960 1976-1977, 1978-2008, and 2009-2016, respectively. The Nash-Sutcliffe efficiency of the calibration and verification periods was 0.78 and 0.76, respectively, and the PBIAS (percent bias) of these periods was −2.1% and −11.1%, respectively. The model performance was good (Daggupati et al., 2016). The correlation coefficient between the simulated and MODIS evapotranspiration was 0.75, the PBIAS was 11.2%, and RMSE was 0.6 mm/d. The reservoir discharge influenced the flow at the outlet of the study area and evapotranspiration. The assessment indices of discharge and evapotranspiration indicated that the simulation had a high accuracy. Therefore, the reservoir simulation was rational. The correlation coefficient of soil moisture was 0.65, which was significantly correlated (Fang, 2011), with a root mean square error of 0.023. Although the objective function did not consider the effects of soil moisture and evapotranspiration, the above results showed that the calibration using discharge was reasonable. where NSE is the Nash-Sutcliffe efficiency; f t is the model simulation value at time t, t = 1,2,…,T; y t is the observed value at time t; and E y is the average value of the observed values at all times.
To illustrate the superiority of the CHC model, the SWAT model was built for the watershed upstream of Heilin hydrological station, and the necessary data are described in Section 2.7. The SWAT performance was compared with the CHC model performance (Table 1). The CHC model's accuracy was better than that of the SWAT model. One reason for this was that only daily precipitation data could be input into the SWAT model, whereas six-hour precipitation data could be fed into the VIC model. Moreover, the SWAT model outputs the total soil moisture of all the layers, whereas the VIC model outputs the soil moisture of each layer individually. Finally, the evapotranspiration module of the CHC model was improved using the water stress method (Allen et al., 1998), which increased the correlation coefficient between the simulated and MODIS evapotranspiration.
To verify the simulation performance of cropland evapotranspiration, two merged grids with high proportion of cropland (more than 80%) were chosen and their correlation coefficients between the evapotranspiration data of CHC simulation and MODIS were compared with that of VIC simulation (shown in Table 2). The CHC error was smaller than that of VIC in the merged grids with high proportion of cropland, which showed that the model's performance improved in the cropland hydrologic simulation after the coupling of the VIC model with the crop model.

Time Series of Crop Water Demand, Consumption, and Deficit
Figure 5 presents the time series of crop water demand, consumption, and deficit, including the three main crop rotations of the study area: spring peanut, summer peanut-winter wheat, and summer corn-winter wheat. The amount of water varied with different crop rotations. For example, the minimum water demand of summer peanuts occurred in 2006, whereas that of spring peanuts occurred in 1979, and the maximum water deficit of summer peanuts occurred in 2001, which differed from that of spring peanuts. This difference was mainly due to the different growth time of spring and summer peanuts. Furthermore, the months during which the leaf area reached its maximum differed, and the larger the leaf area was, the greater the water consumption became. Different water stress and root parameters were relevant to different crops, which led to varying root depth and crop water uptake at the same time. Therefore, crop water differed between different crop types.
The high water demand in some years indicated hot and sunny weather, but this did not necessarily mean large water deficits, because constant short-term precipitation occurred to satisfy crop water consumption. The high water consumption in some years indicated not only hot weather but also high precipitation and irrigation water, which could meet the demands of normal crop growth. Large water deficits corresponded to high water demand and low water consumption in some years, indicating the probability of drought, when the weather was hot but precipitation and irrigation water were too low for normal crop growth. The water consumption of spring peanuts peaked in 2008, between 1995 and 2015, but the area affected by drought was 0 in Junan County (Cheng, 1996(Cheng, -2016, indicating sufficient rainfall and irrigation water in 2008, which enabled normal crop growth. This finding was consistent with the small water deficit of spring peanuts in 2008. In 2001, water deficits in spring peanuts occurred widely, and 27,000 hm 2 of land was affected by drought in Junan County (Cheng, 1996(Cheng, -2016, indicating a severe drought. For the peanut-wheat crop rotation, the maximum water deficit occurred in 2002, and the drought-affected area was 31,401 hm 2 , the greatest drought area was in 1990-2015. Comparison between water deficits and drought areas showed that the water deficit simulation was reasonable. Seasonal differences in crop water variation occurred throughout the year. For example, after the sowing of spring peanuts from May onwards, the water consumption increased until it peaked in July due to the continuous increase in leaf area (Figure 6a). Regardless of the rotation method, the water consumption peaked in July, because the crop leaf area was at its maximum and the weather was hot.
Water demand and consumption also often had similar trends and peaked at a similar time. However, occasionally in summer, the changes in water demand and consumption had opposite trends. When the water demand peaked, the water consumption reached its minimum, and the water Note. ET, evapotranspiration; NSE, Nash-Sutcliffe efficiency.   deficit often peaked correspondingly. During summers with high temperatures and low rainfall, the air humidity gradually decreased, causing water demand to increase. Furthermore, the soil water content gradually decreased, resulting in increased water stress and decreased crop water consumption. This situation was alleviated by rainfall events, which increased water availability and decreased water stress, resulting in increased water consumption. Therefore, point A was the minimum point in Figure 6c. After rainfall, the air humidity often increased and the temperature decreased slightly, reducing the water demand. Therefore, point B was the maximum point in Figure 6c. The water deficit is the difference between water demand and consumption. Thus, the water deficit corresponded to the maximum value (Figure 6c), which was also the maximum value throughout the year.

PBIAS (%) RMSE
In spring, the maximum water deficit corresponded to the maximum water demand and consumption. Because the temperature during spring gradually increased, water demand increased correspondingly, but the soil moisture content was not too low during the rainy season. Furthermore, the water stress effect was small. The effect of the lower water stress coefficient on water consumption was not as high as that of increasing ET 0 . Hence, water consumption continued to increase. However, influenced by water stress, water consumption increased at a much slower rate than water demand did, and this difference was reflected by the increasing water deficit, which eventually peaked.
Furthermore, water consumption and demand differed between different crops (Figure 7). The water demand, consumption, and deficits of corn were high, and those of wheat were low. This was mainly related to the crop LAI max . The LAI max of corn was greater than that of wheat. When LAI max increased, the LAI in Equation 4 increased accordingly. The plant canopy resistance r c obtained using Equation 3 was lower, and the transpiration was higher based on Equation 2. Consequently, the water consumption, demand, and deficits of corn were greater.
Different crops exhibited different water changes. The water demand, consumption, and deficits of peanuts and corn increased sharply after sowing and gradually declined during harvesting. In contrast, wheat water demand and consumption initially increased and then declined before January. After January, the same pattern was observed. This was mainly due to the dual effects of climate conditions and changes in leaf area during crop growth. The leaf area of peanuts and corn gradually increased as these crops grew, and they experienced the transition from spring to summer, which was accompanied by gradually increasing temperatures. Therefore, according to Equations 2, 3, 7-9, and 16, the crop water demand, consumption, and deficits increased. Maturing peanuts and corn exhibited withering leaves and decreased leaf area. The weather also gradually cooled during the abovementioned transition. The cooler weather and withered leaves caused water demand, consumption, and deficits to gradually decrease according to Equations 2, 3, 7-9, and 16.
where T is the daily average temperature (°C). The other variables are described below Equation 2. The change in winter-wheat leaf area was consistent with the above, but because of the long growth period of this crop and complex weather changes, it underwent a transition from autumn to winter and then to spring. The temperature and precipitation initially decreased and then increased, whereas the leaf area first increased and then decreased. Therefore, the increased water demand and consumption (hereinafter water demand and consumption are collectively referred to as water volume) during the seedling stage (October to early November) was related to the leaf area, and the impact of temperature and precipitation was small. Because the leaf area increased during the seedling stage, the relative change was large. Water volume in the tillering period (early November to mid-December) gradually decreased mainly because of gradually decreasing temperature and precipitation. The changes were substantial, resulting in decreased soil water content and potential transpiration and causing the water volume to decrease further. The crop leaf area remained relatively constant. Therefore, the leaf area was no longer a major factor.
Wheat growth slowed during mid-December to late February because of the cooler temperatures and lower water levels in winter. Figure 7c shows low water deficits from October to February related to the low temperatures and small leaf area during this time, and crop water demand and consumption were low, but soil water content was high. Therefore, the degree of water stress was lower. During the renewed growth stage (late February to early April), wheat underwent a transition from winter to spring, and the temperature gradually increased, leading to a gradual increase in potential evapotranspiration and rainfall. Subsequently, during the rainy season, wheat experienced increasing water availability. The mild climate conditions promoted leaf growth (increased leaf area), and the combined effect of climate and leaf growth resulted in gradually increasing water volume. Finally, wheat grew and matured, undergoing jointing, heading, milk ripening, and maturation stages. In the later stages, the leaves gradually withered, and the rainy season was over. Although the temperature increased constantly and the water demand per unit leaf area increased, the total leaf area and total crop water demand and consumption decreased. The precipitation also decreased, leading to increased water stress. Therefore, leaf area became the main influencing factor at this stage because of the large relative changes.
Generally, the water volume of different crops varied greatly, which can be explained by climate and crop growth. The explanations show that the simulation of crop water demand and consumption is reasonable. Crops influenced the external environment, for example, through increased total evapotranspiration from the land to the atmosphere, and vice versa. During the cold and dry winter, crops stopped growing, but during spring, the leaf area increased, which enhanced the total evapotranspiration from the land to the atmosphere. Therefore, interactions between crops and atmosphere existed. These interactions included the two-way effects of crops on the atmosphere and vice versa and were also interdependent (Osborne et al., 2007). Some studies (Y. Siad et al., 2019;Tsarouchi et al., 2014) have emphasized the necessity of coupling crop and hydrological models. The above results showed that increased soil moisture due to increased precipitation affected leaf growth and, subsequently, crop water consumption. In addition to precipitation, irrigation influenced these changes, and the impact of irrigation could not be ignored, and there were also two-way interactions, which would be further analyzed in the Section 4.

Validation of Drought Events Simulated by the CHC Model Considering Irrigation Effects
The CHC model can simulate the process of agricultural drought, but the water deficit cannot be used for the quantification of drought severity, because the water deficit of one day only shows the effect of drought of one day on the crop growth and the reported drought area is a cumulative impact result of drought. Besides, different crops have different levels of crop water consumption, the crop water deficit is variable for different crop types, and the crop water deficit cannot be compared for different crop types. As a result, a drought index called CWAPI (Crop Water Anomaly Percentage Index) was proposed to identify the drought area every year. This drought index can consider the effect of irrigation and can be used to compare the drought state for different crop types. This index was calculated using the following equations: where E t is the maximum plant transpiration on a given day of j (mm); E t,act is the actual transpiration on a given day of j (mm); CWA j is the Crop Water Anomaly index on a given day of j; k i is the weight of CWA j ; CWAI i is the weighted average of CWA of the period i; and m is the number of each year, m = 1,2,…,n.
To validate the irrigation effect on the CWAPI, a drought index called SMAPI (Soil Moisture Anomaly Percentage Index) was compared with the CWAPI, because the SMAPI cannot consider the irrigation effect. The SMAPI calculation and drought identification method can be found in Mao et al. (2016) and Wu et al. (2011).
To validate the effect of irrigation on the drought in the Qingkou River, we calculated the SMAPI and CWA-PI and identified the drought area each year based on the SMAPI and CWAPI. Meanwhile, the reported drought area caused by the drought loss in the study area was collected and the reported drought area was compared with the identified drought area based on the SMAPI and CWAPI. To neglect the effect of the agricultural acreage variation, the ratio of the drought area and agricultural acreage was calculated.
In 1995,2004,2005,2008, and 2013, the reported drought area ratio is zero or small and the identified drought area ratio based on the SMAPI and CWAPI both is zero. In 2006, the reported drought area ratio is 2.4%, but the drought area based on the SMAPI is 88.0%, while the drought area based on the CWAPI is 0%. This result illustrated that CWAPI can reveal the drought state better, because CWAPI considers the effect of irrigation and SMAPI cannot consider the irrigation effect. The similar phenomenons occurred in 1992, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2009, 2011, and 2012. 4. Discussion

Effects of Crop Irrigation on Crop Water Demand, Consumption, and Deficit
To explain the effects of irrigation on crop water consumption, demand, and deficit, the difference of crop water simulated by the CHC model with irrigation module (crops and irrigation, CI) and the CHC model without irrigation module (crops and no irrigation, CNI) in the spring peanut of Qingkou River is shown in Figure 8. Irrigation increased the crop water consumption and demand and the difference of crop water consumption was larger than that of crop water demand. Because irrigation increased the soil moisture (see Section 4.3), the water stress effect decreased, and water consumption increased. Besides, irrigation also increased the LAI of crops due to the larger soil moisture and the water demand and consumption increased. The crop water consumption was mainly affected by the water stress effect and the larger LAI while the crop water demand was mainly affected by the larger LAI. As a result, the difference of crop water consumption between CI and CNI was larger than that of crop water demand and the crop water deficit decreased.

Effect of Crop Irrigation on Hydrological Factors in Typical Years
To analyze the effects of irrigation on the hydrological factors in typical years, daily average soil water content, evapotranspiration, runoff, water demand, water deficits, and flow under high-flow, normal-flow, and low-flow years (1990, 2009, and 2002, respectively) were calculated using the VIC model (no crops and no irrigation, NCNI), the CHC model without irrigation module (crops and no irrigation, CNI), and the CHC model with irrigation module (crops and irrigation, CI). Table 3 shows that soil moisture, transpiration, and evapotranspiration of all soil layers in the CI simulation were higher than those in the CNI simulation in the high-flow year, indicating that irrigation increased soil moisture and evapotranspiration. This finding was consistent with that of previous studies (Chen et al., 2018;Haddeland et al., 2006;Leng et al., 2017). Irrigation had a lower impact on runoff, water demand, and water deficit. The surface soil moisture in the NCNI simulation was higher than that of the crops, and the bottom soil moisture was lower, which was related to the water consumption distribution in each soil layer in the VIC model. The agricultural land in the VIC model was simulated using corn as the default. Corn had a larger leaf area than most other crops did, which enhanced evapotranspiration and reduced runoff. In addition, the VIC model did not consider the small leaf area of seedlings. Tables 4 and 5 show that the relevant conclusions were similar to the above in the normal-flow and low-flow years, but the difference in the hydrological elements between irrigation and non-irrigation scenarios in the low-flow year was greater than that in normal-flow and high-flow years, including differences in evapotranspiration, runoff, and water demand. This was because more irrigation water was used during the low-flow year, and the lower water content of non-irrigated soil resulted in greater water stress.
Regarding the daily average flow of the outlet of the study area, Tables 3-5 indicate that during normal-flow and low-flow years, irrigation increased evaporation and water consumption. Therefore, the irrigation flow decreased. In high-flow years, the amount of irrigation was low, but it nonetheless increased soil water content. Hence, runoff increased, leading to slight increases in flow. The water content of the first two soil layers in the VIC model was higher, but that of the third layer was lower, which led to a lower base flow. Consequently, the overall flow was lower. The difference between the scenarios with no crops and those with crops increased with increasing annual precipitation.

Effects of Crop Irrigation on Changes in Soil Moisture
To analyze the impact of irrigation on the simulation of soil moisture content and the validity of the irrigation module, the soil moisture values in the NCNI, CNI, and CI simulations were compared in the low-flow, normal-flow, and high-flow years (2002, 2009, and 1990, respectively; Figure 9). Figure 9a shows that the soil moisture in the CI simulation was similar to that in the CNI simulation. This was mainly because a high amount of irrigation water was not required in the high-flow year. During winter, crops were irrigated  because of the lower rainfall. Hence, the soil moisture in the CI simulation was higher than that in the CNI simulation in winter. The soil moisture simulated by the VIC model was generally higher than that simulated by the CHC model. The crop water uptake was mainly concentrated in the third soil layer in the VIC model, with lower water uptake from the surface soil layer. In addition, the monthly leaf area in the VIC model remained constant, and cases of small leaf area in the crop seedling stage were ignored. Therefore, the CHC model simulation was more reasonable. Figure 9b shows that the difference between soil moisture simulations in the second soil layer by the CI and CNI was greater than that in the first layer, and the difference in the duration was also longer. This was mainly because the second soil layer was thicker than the first layer, and the saturated water content of the second layer was much higher than that of the first layer. The water supply of the second soil layer was affected by the first soil layer, and the recharge rate of the former was slower. Therefore, the saturated water content in the second soil layer was reached over a longer period. The soil moisture simulated by different methods was similar only after the soil water content reached saturation. As a result, the difference and its duration in the second soil layer were greater than those in the first layer. This also showed that the CHC model simulation was more reasonable.
The soil moisture in the second layer simulated by the CHC model was lower than that simulated by the VIC model during the non-flood season but higher than that simulated by the VIC model during the flood season. This was mainly because the evaporation from the second layer in the non-flood season simulated by the CHC model was higher than that by the VIC model, which was related to the difference in the evaporation proportion between different soil layers in the two models. Figure 9c shows that the third layer was more strongly affected by water accumulation. Hence, the difference in the duration of soil moisture between the CI and CNI simulations was greater. The greater difference between the VIC and CHC models in the third layer was mainly due to different parameter values of the third layer. The soil moisture simulated by the VIC model was lower than that simulated by the CHC model in the third layer because the third layer in the VIC model was less thick than that in the CHC model because of the different calibrations of each model. According to the simulation results, the gap in soil moisture gradually widened between the third soil layer in the VIC and CHC models since the first year. The crops in the VIC model mainly consumed soil water from the third layer, and those in the CHC model mainly consumed water from the first layer. The water consumption varied throughout the year. Baseflow was proportional to the soil water content of the third layer. Therefore, they would eventually enter a different equilibrium state (when the input and output would be nearly equal). Figures 10 and 11 show the soil moisture contents of the first, second, and third layers in the NCNI, CNI, and CI simulations during the normal-flow and low-flow years in the watershed upstream of Heilin hydrological station. The differences between the soil water contents in the three simulations in the normal-flow and low-flow years were similar to those in the high-flow year. Comparison of the irrigation water in the high-flow, normal-flow, and low-flow years revealed that the irrigation water in the high-flow year was lower than that in the normal-flow and low-flow years. Because precipitation was higher during the high-flow year, the soil water content was greater, and less irrigation water was required. The duration of irrigation time overlapped in these three years during January to March, which is the non-flood season, with less rainfall, lower soil water content, and greater irrigation water demand. Irrigation had a greater impact on     soil water content in the low-flow year. However, if precipitation occurred after irrigation, the difference in soil moisture between the different years would be eliminated. The impact of irrigation on the second layer had hysteresis and amplification effects, and the latter was stronger in the low-flow year. Because less rain fell during the low-flow year, the probability of soil water content reaching field capacity was lower. Hence, it was difficult to eliminate the difference in soil water content between the CI and CNI simulations. As long as this difference in soil water contents existed, the differences would continue to occur and accumulate.

Effects of Crop Irrigation on the Correlation Coefficients Between the Simulation and Observation
To discuss if including irrigation would improve the skill of the model, the CHC model was built in a larger area including the Qingkou River. This area mainly includes Jiangsu province. The parameters were determined based on the references (Liu, 2017;Mao et al., 2016;Wu et al., 2007). The soil moisture and ET were calculated in the Jiangsu province with and without the irrigation module respectively and the difference between the correlation coefficient (R) of MODIS ET and evapotranspiration simulated by the CHC model with irrigation module and the R-value of MODIS ET and the evapotranspiration simulated by the CHC model without the irrigation module is shown in Figure 12. The R-value was improved in 56.3% of all the merged grids in Jiangsu province. The maximum of the R difference was 0.25 while the minimum of the R difference was −0.06. The R improvement of soil moisture using the irrigation module is shown in Figures 13 and 14. The R-value improved in 59.6% of all the merged grids with the soil moisture monitoring stations. The maximum of the R difference was 0.28, while the minimum of the R difference was −0.04. The correlation coefficient differences of different values are evenly distributed in the study area. These results showed that the irrigation module improved the performance of evapotranspiration and soil moisture simulation in more than half of the merged grids.

Conclusions
To overcome the oversimplification of crop and irrigation modules in hydrological models, the VIC model was coupled with crop and gridded routing modules to simulate the hydrology and crop growth. The reservoir and irrigation modules were introduced to incorporate the impact of irrigation on crop growth. The transpiration simulation module of the crop model was improved, based on the water stress simulation method recommended by the Food and Agriculture Organization of the United Nations (FAO). A CHC model was built to simulate soil moisture and transpiration in the upstream watershed of the Qingkou River at Heilin hydrological station on the border of Jiangsu and Shandong Provinces of China. The accuracy of the CHC model simulation results was higher than that of the SWAT simulation results, and the CHC model's performance improved in the cropland hydrologic simulation after coupling of the VIC model with the crop model. The crop water consumption and drought events simulated by the CHC model considering irrigation effects were validated and the impact of crop irrigation was analyzed. The following conclusions were drawn: 1. The CHC performance was validated using the SWAT results, MODIS data, and drought event records. The CHC model's accuracy was better than that of the SWAT model in terms of discharge, soil moisture, and evapotranspiration; the CHC model's performance improved in the cropland evapotranspiration after the coupling of the VIC model with the crop model; the simulated ratio of the drought area and agricultural acreage considering the irrigation effect is close to the reported drought area ratio. 2. The crop water temporal variation in the three crop rotations in the watershed upstream of Heilin hydrological station and its validity were analyzed. The simulated water deficit was verified using statistical data from the drought-affected area in Junan County. Seasonal differences in crop water were reasonable throughout the year due to the different leaf area and weather. Furthermore, difference validity between crop water demand, consumption, and deficits was analyzed. In summer, meteorological and soil water conditions cause the crop water demand to be minimum to correspond to maximum water consumption, which generated the maximum water deficit. The water gaps between different crops were mainly related to the two-way interactions between crops and atmosphere. 3. The effect of irrigation on the crop water deficit and cumulative value of hydrological elements in the study basin was analyzed. Irrigation increases the crop water demand and consumption and decreases the crop water deficit. Irrigation increased soil water content and evapotranspiration, consistent with the findings of previous studies. The VIC model's proportion of water uptake in each soil layer increased the surface soil water content, and the default corn (which represented all the crops) resulted in lower soil water content in the third layer. The amount of irrigation water was high in the low-flow year, and when no irrigation occurred, water stress strongly affected crop growth, resulting in greater differences between hydrological factors with and without irrigation. 4. The effect of irrigation on soil water content and evapotranspiration was analyzed. The influence of irrigation on the soil water content of the second layer had hysteresis and amplification effects, and the latter was stronger in the low-flow year. Due to the lower soil water content and higher water demand in the low-flow year, the irrigation volume was higher than that in the high-flow year. The irrigation

Data Availability Statement
The Landsat 8 data were downloaded from http://www.gscloud.cn. The DEM data from the Shuttle Radar Topography Mission were obtained from https://earthexplorer.usgs.gov/. The water area data were downloaded according to Pekel et al. (2016). The grid vegetation type data were obtained according to Hansen et al. (2000). The global soil map was downloaded according to Reynolds et al. (2000). The parameters of the different crop types were obtained from https://swat.tamu.edu/.