Using an ETWatch (RS)-UZF-MODFLOW Coupled Model to Optimize Joint Use of Transferred Water and Local Water Sources in a Saline Water Area of the North China Plain

In the saline water area of our research, deep groundwater was over-pumped for agricultural irrigation which resulted in a decline of the deep groundwater level and an increase in the shallow groundwater table. Soil salination was also aggravated due to the strong evapotranspiration (ET) in the shallow groundwater areas, where ET removes water vapor from the unsaturated zone (ETu), and the groundwater (ETg). Joint utilities of multiple water sources of transferred water and local shallow and deep groundwater are essential for reasonable management of irrigation water. However, it is still difficult to distinguish ETu and ETg in coupled management of unsaturated zone and groundwater, which account for the water balance in utilities of multiple water sources in a regional scale. In this paper, we used an RS-based ETWatch model as a source of evapotranspiration data coupled with UZF-MODFLOW, an integrated hydrological model of the unsaturated–saturated zone, to estimate the ETg and ETu on a regional scale. It was shown that the coupled model (ETWatch-UZF-MODFLOW) avoids the influence of ETu on the groundwater balance calculation and improves the accuracy of the groundwater model. The model was used in the simulation and prediction of groundwater level. The eastern North China Plain (NCP) was selected as the study area where shallow groundwater was saline water and deep groundwater cone existed. We compared four different scenarios of irrigation methods, including current irrigation scenario, use of saline water, limited deep groundwater pumping, use of multiple water sources of transferred water and local groundwater. Results indicate that the total ETg for the four scenarios in the study area from 2013 to 2030 is 119 × 108 m3, 81.9 × 108 m3, 85.0 × 108 m3, and 92.3 × 108 m3, respectively, and the proportion of ETg to total ET was 6.85%, 4.79%, 4.97%, 5.37%. However, in regions where the groundwater depth is less than 3 m, ETg accounts for 12% of the total ET, indicating that groundwater was one of the main sources of evapotranspiration in shallow groundwater depth area.

land use/cover, spatio-temporal rainfall and evapotranspiration evaluations [38]. These remote sensing data are often used for coupling of surface water models. However, the RS contribution to groundwater hydrology and groundwater resource evaluation is less distinct, and is therefore less well known [39]. This is because the coupling of RS data to groundwater models requires consideration of the unsaturated zone data for connection. The data in the unsaturated zone are usually not directly monitored by RS technology. Therefore, it is unreasonable to simply couple RS data with groundwater models while ignoring the spatial heterogeneity of the unsaturated zone and its internal moisture variation.
The lowland area of the North China Plain is located in the east-central part of the North China Plain (NCP), which has large areas of low and medium-yielding farmland, and is a key area for regional agricultural development, as well as a major area for the Bohai Sea Granary Science and Technology Demonstration Project to increase grain production. The goal of the Bohai Sea Granary Science and Technology Demonstration Project is to increase grain production in low and medium-yielding fields, and is expected to increase grain production in the Bohai region by five billion kilograms by 2020 [40], which is of great significance for ensuring national food security. However, agricultural irrigation is the key to guaranteeing high and stable grain production, and the eastern area of the NCP is also an area with a severe shortage of freshwater resources, with per capita water resources per unit acreage in the region being only 1/12 of the national level [41], making it one of the most water-scarce regions in China.
Therefore, in this paper, ETWatch, a remote sensing product, is coupled with an integrated hydrological model including unsaturated zone water (UZF package) and groundwater (MODFLOW). Wu et al. [42,43] ever reported a brief description of ETWatch and its application in monitoring regional evapotranspiration. UZF packages use kinematic wave approximation of the one-dimensional Richards' equation through a homogeneous unsaturated zone and can provide a method that is well suited to large-scale models (field scale, regional scale) [44]. PTF functions can provide soil water parameters with fine spatial discretization, which can help to calculate the spatial distribution and variation of water content in unsaturated zone, and thus significantly improve the accuracy of this regional scale coupling model. The coupled model will be applied in the eastern area of the NCP, where deep groundwater was over-exploited for agricultural irrigation since the saline area distributed widely in shallow groundwater.

Study Area
The study area is located in the Cangzhou city which is in the eastern North China Plain (NCP) (37 • 30 -39 • 0 N, 115 • 42 -117 • 58 E) ( Figure 1). The slope of the area varies from 1/8000 to 1/15,000. The ground elevation is 2-15 m. The climate is a warm temperate sub-humid continental monsoon climate with a perennial average temperature of 12.2 • . It is cold and dry in winter and hot and rainy in summer. The mean annual precipitation in the area varies from 400 to 500 mm, more than 80% of which is received during July-September period. The mean annual evaporation is 1919 mm, with 34.5%, 36.9%, 20.7%, and 7.9% occurring in spring, summer, autumn, and winter, respectively [45].
The predominant sedimentary strata in the eastern NCP area are quaternary loose sedimentary strata with a thickness of about 350-550 m, with each stratum being a single layer (Figure 1). The aquifer is mainly silt and fine sand. The main stratigraphy is composed of unconsolidated sediments of Q4 origin. Four (I, II, III, and IV) aquifer layers consist of silty clay and clay [46]. Holocene Series aquifer group (depth of the aquifer floor is about 40 m, salty water type), Upper Pleistocene Series aquifer group (depth of the aquifer floor is about 120-250 m, salty water type), Middle Pleistocene Series III aquifer group (depth of the aquifer floor is 350-430 m, freshwater type), Lower Pleistocene Series aquifer group (depth of the aquifer floor is 350-550 m, freshwater type) ( Figure 2). Shallow aquifer includes Layers I and II, which is mainly composed of salty water and the degree of exploitation is Water 2020, 12, 3361 4 of 25 relatively low. Deep aquifer includes Layers III and IV with freshwater and it is the main pumping layer in the study area [46][47][48]. The predominant sedimentary strata in the eastern NCP area are quaternary loose sedimentary strata with a thickness of about 350-550 m, with each stratum being a single layer ( Figure 1). The aquifer is mainly silt and fine sand. The main stratigraphy is composed of unconsolidated sediments of Q4 origin. Four (I, II, III, and IV) aquifer layers consist of silty clay and clay [46]. Holocene Series aquifer group (depth of the aquifer floor is about 40 m, salty water type), Upper Pleistocene Series aquifer group (depth of the aquifer floor is about 120-250 m, salty water type), Middle Pleistocene Series III aquifer group (depth of the aquifer floor is 350-430 m, freshwater type), Lower Pleistocene Series aquifer group (depth of the aquifer floor is 350-550 m, freshwater type) ( Figure 2). Shallow aquifer includes Layers I and II, which is mainly composed of salty water and the degree of exploitation is relatively low. Deep aquifer includes Layers III and IV with freshwater and it is the main pumping layer in the study area [46][47][48].   The volumes of groundwater exploitation in 2011 and 2012 were 10.75 × 10 8 m 3 and 10.70 × 10 8 m 3 , respectively, 73% and 72% of which were pumped for agricultural exploitation. The amounts of groundwater exploited from shallow aquifers in these two years are 3.5 × 10 8 m 3 and 3.2 × 10 8 m 3 , respectively, accounting for about 32.5% and 29.9% of total groundwater exploitation. Due to the widely Water 2020, 12, 3361 5 of 25 distributed salty water, the extent of exploitation of shallow groundwater was low. Groundwater exploitation in the study area is mainly concentrated in the deep aquifer group.
The main crops in the study area are wheat-maize and cotton, and their distribution is shown in Figure 3. The volumes of groundwater exploitation in 2011 and 2012 were 10.75 × 10 8 m 3 and 10.70 × 10 8 m 3 , respectively, 73% and 72% of which were pumped for agricultural exploitation. The amounts of groundwater exploited from shallow aquifers in these two years are 3.5 × 10 8 m 3 and 3.2 × 10 8 m 3 , respectively, accounting for about 32.5% and 29.9% of total groundwater exploitation. Due to the widely distributed salty water, the extent of exploitation of shallow groundwater was low. Groundwater exploitation in the study area is mainly concentrated in the deep aquifer group.
The main crops in the study area are wheat-maize and cotton, and their distribution is shown in Figure 3.

Model Description
In this paper, remote sensing data and pedo transfer function parameters proposed by Zhang and Marcal [50] at the regional scale were coupled with the groundwater model to calculate the ETg ( Figure 4). ETWatch, an ET monitoring system based on remote sensing technology, is a data source for ET. In addition, there is the modflow-unsaturated zone flow (UZF) package to determine the proportion of ETg and ETu in ET, respectively ( Figure 5).This coupled model was used in the simulation and prediction of groundwater irrigation in the eastern NCP to test performance. ETWatch and Pedo transfer function parameters are both the work of others and are available to us as mature tif-formatted data products; we use ArcGIS as a tool to slice these spatial data according to a grid of groundwater models and calculate each grid for each stress period. Pedo transfer function parameters are global data and can be used for any regional study; while ETWatch is a regional product, other remote evapotranspiration data products such as MODIS, SEBAL, etc., can be used instead of ETWatch for model coupling.

Model Description
In this paper, remote sensing data and pedo transfer function parameters proposed by Zhang and Marcal [50] at the regional scale were coupled with the groundwater model to calculate the ET g ( Figure 4). ETWatch, an ET monitoring system based on remote sensing technology, is a data source for ET. In addition, there is the modflow-unsaturated zone flow (UZF) package to determine the proportion of ET g and ET u in ET, respectively ( Figure 5). This coupled model was used in the simulation and prediction of groundwater irrigation in the eastern NCP to test performance. ETWatch and Pedo transfer function parameters are both the work of others and are available to us as mature tif-formatted data products; we use ArcGIS as a tool to slice these spatial data according to a grid of groundwater models and calculate each grid for each stress period. Pedo transfer function parameters are global data and can be used for any regional study; while ETWatch is a regional product, other remote evapotranspiration data products such as MODIS, SEBAL, etc., can be used instead of ETWatch for model coupling.

ETWatch
A brief description of ETWatch and an application for monitoring regional evapotranspiration were reported by Wu et al. [42,43]. ETWatch integrates SEBAL and surface energy balance system (SEBS) to estimate surface fluxes under clear-sky conditions together with the Penman-Monteith approach to calculate daily ET, based on a surface resistance model, meteorological data and surface parameters from remote sensing. A brief description of ETWatch is given in the following text, and the main procedure is demonstrated in Figure 6.

ETWatch
A brief description of ETWatch and an application for monitoring regional evapotranspiration were reported by Wu et al. [42,43]. ETWatch integrates SEBAL and surface energy balance system

ETWatch
A brief description of ETWatch and an application for monitoring regional evapotranspiration were reported by Wu et al. [42,43]. ETWatch integrates SEBAL and surface energy balance system (SEBS) to estimate surface fluxes under clear-sky conditions together with the Penman-Monteith approach to calculate daily ET, based on a surface resistance model, meteorological data and surface parameters from remote sensing. A brief description of ETWatch is given in the following text, and the main procedure is demonstrated in Figure 6.

UZF Package and Pedo transfer Functions Parameters
Unsaturated-zone flow processes are simulated using the UZF1 package of MODFLOW [53]. UZF1 uses the method of characteristics for solving a kinematic wave approximation for onedimensional downward vertical flow in the unsaturated zone that is derived by neglecting the diffusive term in Richards' equation. UZF1 also relies on the assumption of uniform hydraulic properties in the unsaturated zone for each vertical column of model cells [54]. The unsaturated-zone flow equation solved by UZF1 that is coupled to the groundwater flow equation is where θ is volumetric water content in the unsaturated zone [L 3 ·L −3 ], K(θ) is the vertical hydraulic conductivity as a function of θ[L·T −1 ], and qET is the rate of ET removal from the unsaturated zone above the water table (expressed per unit depth [L·T −1 ·L −1 ]). The kinematic-wave equation provides an efficient method for simulating unsaturated-zone flow, ET, and recharge over regional scales. If ET potential is not satisfied by unsaturated-zone water, then MODFLOW-NWT also can simulate ET derived from groundwater upfluxon the basis of a linear function of water table depth [55]. Additional detail of the approach used in UZF1 and its advantages are described in Niswonger et al. [53] and Niswonger and Prudic [56].
Eastern NCP was discretized into 3 layers, according to the lithology of sediments and aquifer characteristics. The first layer was the salty water aquifer. The second and third aquifer was the confined aquifer. The grid size of the simulation area is set as 1500 m ×1500 m, with a total subdivision into 111 rows ×126 columns, 7201 active units. The calibration period of model simulation is from January 2011 to December 2011. Each stress period is the corresponding natural month and the time step is one day. The validation period of the model is from January 2012 to December 2012, and each stress period is the corresponding natural month, and the time step is one day.

UZF Package and Pedo Transfer Functions Parameters
Unsaturated-zone flow processes are simulated using the UZF1 package of MODFLOW [53]. UZF1 uses the method of characteristics for solving a kinematic wave approximation for one-dimensional downward vertical flow in the unsaturated zone that is derived by neglecting the diffusive term in Richards' equation. UZF1 also relies on the assumption of uniform hydraulic properties in the unsaturated zone for each vertical column of model cells [54]. The unsaturated-zone flow equation solved by UZF1 that is coupled to the groundwater flow equation is where θ is volumetric water content in the unsaturated zone [L 3 ·L −3 ], K(θ) is the vertical hydraulic conductivity as a function of θ[L·T −1 ], and q ET is the rate of ET removal from the unsaturated zone above the water table (expressed per unit depth [L·T −1 ·L −1 ]). The kinematic-wave equation provides an efficient method for simulating unsaturated-zone flow, ET, and recharge over regional scales. If ET potential is not satisfied by unsaturated-zone water, then MODFLOW-NWT also can simulate ET derived from groundwater upfluxon the basis of a linear function of water table depth [55]. Additional detail of the approach used in UZF1 and its advantages are described in Niswonger et al. [53] and Niswonger and Prudic [56]. Eastern NCP was discretized into 3 layers, according to the lithology of sediments and aquifer characteristics. The first layer was the salty water aquifer. The second and third aquifer was the confined aquifer. The grid size of the simulation area is set as 1500 m × 1500 m, with a total subdivision into 111 rows × 126 columns, 7201 active units. The calibration period of model simulation is from January 2011 to December 2011. Each stress period is the corresponding natural month and the time step is one day. The validation period of the model is from January 2012 to December 2012, and each stress period is the corresponding natural month, and the time step is one day.

MODFLOW-NWT
MODFLOW-NWT is a Newton-Raphson formulation for MODFLOW-2005 to improve solution of unconfined groundwater-flow problems. MODFLOW-NWT is a standalone program that is intended for solving problems involving drying and rewetting nonlinearities of the unconfined groundwater-flow Water 2020, 12, 3361 8 of 25 equation. The model is based on solving the 3D partial differential equation of groundwater flow using the finite-difference method. The governing equation for groundwater flow is as follows [57], where K xx , K yy , K zz are values of hydraulic conductivity along the x, y, z coordinate axes, which are assumed to be parallel to the major axes of hydraulic conductivity (L·T −1 ); h is the groundwater level (L); W are the source/sink terms, with W < 0 for flow out of the groundwater system, and W > 0 for flow into the system (T −1 ); S s is the specific storage of the porous material (L −1 ); t is time (T

Data
The data required for model construction include elevation data, meteorological data, groundwater table data, land use data, crop planting data, pumping data, and hydrogeological parameters.  [49,56]. The initial values of hydrogeological parameters are mainly set with reference to the results of previous research [46,[57][58][59][60].
In addition, Zhang and Marcal [50] applied PTF (Pedo transfer functions) to a high-resolution (1 km) global database of soil properties and generated the corresponding high-resolution maps, including estimated Kosugi parameters, saturated hydraulic conductivity, derivatives such as field moisture capacity, wilting point, plant available water, etc. In this paper, the soil parameters in the UZF model, including saturated water content θ s and residual water content θ r , were from the research results of Zhang and Marcal. They provided a download of the results at https://doi.org/10.7910/DVN/UI5LCE.

Initial Conditions and Lateral Boundary
Initial groundwater level: The initial groundwater level of the aquifer was obtained using the Kriging interpolation method, based on groundwater level data from January 2011. The lower limit of plant available water (wilting coefficient) for the study area was used as the initial soil water content for the UZF calculation. The lower limit value of plant available water (wilting coefficient) diagram was obtained from the calculation results of Zhang and Marcal [50].
Lateral boundary conditions: The shallow groundwater in the study area generally flows from southwest to northeast, and there is a groundwater table depression cone in the western part of Suning county and the southern part of Nanpi county, where the groundwater near the cone converges to the cone center; the deep groundwater flows mainly from all directions to the complex cone center Qingxianand Cangxian counties (Figure 1), and according to the groundwater flow map of the study area, Darcy's law is used to determine the lateral recharge and discharge of the boundary. The finite-difference method is used in the model, so the boundaries are regular meshes. The large size of the regular meshes will make the mesh boundaries vary greatly from the actual boundaries, so we refine the boundary meshes to improve the accuracy of the boundaries. Since the boundaries are mostly constant head boundaries and general head boundaries, the refinement of the mesh boundaries does not make it more difficult to set the boundaries when the values are given. We set up some wells along the general head boundary and conduct pumping tests to get more accurate hydrogeological parameters to calculate the boundary flow.

Sinks and Sources
Recharge from Rainfall Infiltration is the most important recharge in the study area, which is mainly affected by rainfall, shallow aquifer's depth and the lithology of the vadose zone. Groundwater recharge from field irrigation was calculated as follows, where Q pre is the amount of precipitation recharge (m 3 ), α i is the precipitation infiltration coefficient corresponding to region i, P i is the amount of precipitation in region i (m), F i is the calculated area of area i (m 2 ). The rainfall infiltration coefficient is influenced by many factors, such as precipitation characteristics, soil infiltration performance, ground surface conditions, and submerged depth, etc. The choice of the rainfall infiltration coefficient in this study is mainly based on the previous research results of the hydrogeological zones of the North China Plain [46] (Table 1).  We used the runoff data from hydrology stations to calculate the amount of leakage by the water balance method, and inputted leakage data into the river package.

Groundwater Pumping Data
The pumping data were from the hydrology department of the local government; they included annual pumping data of 20 counties in eastern NCP.

Calibration and Validation
The monthly groundwater level (GWL) data from 59 observation wells (during the period 2011-2012) were used to calibrate and validate the model. The simulated GWLs were thus aggregated into the monthly values for comparison. The observed daily dataset in 2011 and in 2012 were applied, respectively, for model calibration and validation. Parameters referring to the hydraulic conductivity and specific yield were calibrated through an iterative process using a trial-and-error method; they were relatively uncertain in the study area. The amount of water pumped from the wells were not calibrated either; the only pumping data we obtained was the overall water consumption of the farmland in a region, but the distribution of wells we cannot specify, so we completed the calibration process by redistributing the amount of water pumped from each pumping well. The root mean square (RMS), normalized root mean square (NRMS), and the correlation coefficient were used as indicators of goodness of fit.

Scenarios Setting
In the eastern NCP, the most traditional and dominant irrigation method was pumping the deep aquifer's groundwater. Pumping deep aquifer for irrigation has led to deep aquifer groundwater levels continuously declining and the formation of large groundwater cones. The joint utilization of multiple water sources in agriculture was to solve the water shortage of eastern NCP.
The government has limited or forbidden exploitation of deep groundwater in some groundwater cone area during recent years. Additionally, the transfer of the Yellow River water to the lowland area of the NCP is part of South-to-North Water Transfer Project [48]. Much water transferred from the Yellow River into the lowland area will be used for agricultural irrigation in the future. In addition to transferred water, combined use of salty water from shallow aquifers with fresh water from deep aquifers for irrigation is a very good method. It can not only reduce the shallow groundwater level and alleviate soil salinization, but also save deep freshwater resources.
The use of saline water for irrigation purposes is increasingly popular in some countries and regions including India, China, Italy, and the USA [61][62][63][64]. Particularly, drip irrigation with saline water was widely applied to reclaim saline soil for planting cotton in the semiarid and arid saline-alkali region of China [65,66], and for constructing ecological landscapes in coastal regions [67,68].
Therefore, four scenarios of agricultural water use are set up to simulate the influence of various agricultural water use methods on groundwater. These scenarios considered current water consumption, limiting pumping of deep groundwater and combined utilization of salty water and transferred water ( Table 2). For scenario 1, only groundwater is used. According to the amount and location of groundwater pumping in 2012, the change in groundwater level was predicted.
For scenario 2, saline water and deep groundwater are used. According to previous research on saline water irrigation, saline water of shallow aquifer can be used to replace deep aquifer's fresh water during wheat irrigation [69]. It was estimated that the amount of deep aquifer's fresh water replaced by shallow saline water was 1.65 × 10 8 m 3 if the planted area of winter wheat in the study area replaced irrigation of deep fresh groundwater with shallow saline water.
For scenario 3, transferred water and groundwater are used. The amount of various water-saving measures (such as agricultural planting structure adjustment, water and fertilizer integration technology, conservation tillage project, sprinkling irrigation and drip irrigation) and external water diversion (from the Yellow River, Zhanghe river and reservoirs) was used to replace the local groundwater irrigation( Table 2).Water saved by agricultural water-saving measures 0.88 (10 8 m 3 ), deep aquifer's fresh water replaced by shallow saline water 0.46 (10 8 m 3 ).
For scenario 4, saline water, surface water and deep groundwater are used in combination. The water-saving measures in scenario 2 and 3 are adopted, and the amount of industrial and domestic groundwater consumption to be replaced by the South-to-North water diversion project is calculated. Total industrial and domestic water consumption is 267 × 10 8 m 3 . For scenario 2, saline water and deep groundwater are used. According to previous research on saline water irrigation, saline water of shallow aquifer can be used to replace deep aquifer's fresh water during wheat irrigation [69]. It was estimated that the amount of deep aquifer's fresh water replaced by shallow saline water was 1.65 × 10 8 m 3 if the planted area of winter wheat in the study area replaced irrigation of deep fresh groundwater with shallow saline water.

Model Calibration and Validation
For scenario 3, transferred water and groundwater are used. The amount of various watersaving measures (such as agricultural planting structure adjustment, water and fertilizer integration technology, conservation tillage project, sprinkling irrigation and drip irrigation) and external water diversion (from the Yellow River, Zhanghe river and reservoirs) was used to replace the local groundwater irrigation (Table 2).Water saved by agricultural water-saving measures 0.88 (10 8 m 3 ), deep aquifer's fresh water replaced by shallow saline water 0.46 (10 8 m 3 ).
For scenario 4, saline water, surface water and deep groundwater are used in combination. The water-saving measures in scenario 2 and 3 are adopted, and the amount of industrial and domestic groundwater consumption to be replaced by the South-to-North water diversion project is calculated. Total industrial and domestic water consumption is 267 × 10 8 m 3 .

Model Calibration and Validation
The  Comparison between the simulated and observed GWLs for some observation wells during the simulation period (2011-2012) is described in Figure 8. These observation wells are located in different areas of the study area. Results showed that the simulated GWLs matched well with the observed ones both for calibration and validation. This indicates that the groundwater model coupled with the ETWatch has good performance results. Inter annual patterns of groundwater level variation (water level rise after the rainy season) are well fitted with small errors. The discrepancy between Comparison between the simulated and observed GWLs for some observation wells during the simulation period (2011-2012) is described in Figure 8. These observation wells are located in different areas of the study area. Results showed that the simulated GWLs matched well with the observed ones both for calibration and validation. This indicates that the groundwater model coupled with the ETWatch has good performance results. Inter annual patterns of groundwater level variation (water level rise after the rainy season) are well fitted with small errors. The discrepancy between simulated and observed data is primarily due to the simulated data, but rather the inaccuracy of the observed data. The observed groundwater level fluctuates several times during the year, mainly because we do not have dedicated observation wells, some of the groundwater level probes are installed in wells near residential areas, and although these areas are centrally water supplied, the wells are still used (for non-agricultural purposes). Therefore, there are temporary fluctuations in the observed groundwater level data, whereas we mainly simulate the overall trend of the groundwater level over the entire hydrological year, and these temporary fluctuations can be ignored because such short-term fluctuations do not affect the pattern of the groundwater level over the entire hydrological year. In addition to these small fluctuations that cannot be fitted, the overall discrepancy between simulated and observed data may be due to the accuracy of our initial data. wells are still used (for non-agricultural purposes). Therefore, there are temporary fluctuations in the observed groundwater level data, whereas we mainly simulate the overall trend of the groundwater level over the entire hydrological year, and these temporary fluctuations can be ignored because such short-term fluctuations do not affect the pattern of the groundwater level over the entire hydrological year. In addition to these small fluctuations that cannot be fitted, the overall discrepancy between simulated and observed data may be due to the accuracy of our initial data.  Spatial distribution of the main calibrated hydrogeological parameters (i.e., Kx and Sy) is presented in Figure 9. Spatial distribution of the main calibrated hydrogeological parameters (i.e., K x and S y ) is presented in Figure 9. We have many years of experimental experience in this area and have accumulated a large number of measured data to calculate the value of the conductivity coefficient. In addition, we use a pumping test to calculate the hydraulic conductivity of the area without data, and we interpolate the hydraulic conductivity of this area according to the measured hydraulic conductivity and the collected data. The hydraulic conductivity in this region is basically isotropic, so KX = KY. The reason for setting Kx = Ky = 10Kz is that the model has three aquifers, and the value of Kz averages the hydraulic conductivity of the aquifer and the aquitard; we define the value of Kz in this way because we do not consider the vertical flow within the aquifer.

Prediction of Groundwater Dynamics and Simulation with Different Irrigation Scenarios
By using the calibrated ETWatch-UZF-MODFLOW model, the groundwater level of eastern NCP from 2013 to 2020 was run to show the dynamics of shallow groundwater affected by different utilization of agricultural irrigation ways. Figure 10 shows the shallow groundwater table in 2018, 2024, and 2030 for four scenarios, respectively. The flow direction of shallow groundwater is mainly from Southwest to Northeast. In the western part of the area, the groundwater depth is morethan5 m, while the groundwater depth in the eastern coastal area is less than 5 m, and the shallowest part is less than 1 m. Ground water depth and evapotranspiration are closely related.
The simulation results show that the groundwater table varies little. Groundwater table of Scenarios 2 and 4 declined due to the combined irrigation of shallow saline groundwater in these two scenarios. Additionally, groundwater table in the eastern area in 2030 is about 1-2 m lower than that in 2012. For Scenarios 1 and 3, intensive irrigation (deep groundwater and surface water sources) results in a further increase in the shallow groundwater level in the eastern region, which will accelerate the local soil salinization extent (Figures 10 and 11). . Spatial distribution of the calibrated shallow K x (a), deep K x (b) and S y (c) in aquifers, with K x = K y = 10 K z (m d −1 ) (K x , K y , K z are the values of hydraulic conductivity along the x, y and z coordinate axes, and S y is the specific yield. Other calibration parameters include storativity (d) and rainfall infiltration coefficient (e).
We have many years of experimental experience in this area and have accumulated a large number of measured data to calculate the value of the conductivity coefficient. In addition, we use a pumping test to calculate the hydraulic conductivity of the area without data, and we interpolate the hydraulic conductivity of this area according to the measured hydraulic conductivity and the collected data. The hydraulic conductivity in this region is basically isotropic, so K X = K Y . The reason for setting Kx = Ky = 10Kz is that the model has three aquifers, and the value of Kz averages the hydraulic conductivity of the aquifer and the aquitard; we define the value of Kz in this way because we do not consider the vertical flow within the aquifer.

Prediction of Groundwater Dynamics and Simulation with Different Irrigation Scenarios
By using the calibrated ETWatch-UZF-MODFLOW model, the groundwater level of eastern NCP from 2013 to 2020 was run to show the dynamics of shallow groundwater affected by different utilization of agricultural irrigation ways. Figure 10 shows the shallow groundwater table in 2018, 2024, and 2030 for four scenarios, respectively. The flow direction of shallow groundwater is mainly from Southwest to Northeast. In the western part of the area, the groundwater depth is more than 5 m, while the groundwater depth in the eastern coastal area is less than 5 m, and the shallowest part is less than 1 m. Ground water depth and evapotranspiration are closely related.  The simulation results show that the groundwater table varies little. Groundwater table of Scenarios 2 and 4 declined due to the combined irrigation of shallow saline groundwater in these two scenarios. Additionally, groundwater table in the eastern area in 2030 is about 1-2 m lower than that in 2012. For Scenarios 1 and 3, intensive irrigation (deep groundwater and surface water sources) results in a further increase in the shallow groundwater level in the eastern region, which will accelerate the local soil salinization extent (Figures 10 and 11).

Evapotranspiration Lost from Unsaturated Zone and Groundwater
The water balance of unsaturated zone and saturated zone can be calculated by using the combined ETWatch-UZF-MODFLOW model. This was used for estimating the water loss of the study area in detail.
For the water balance in the unsaturated zone, the main input sources of water are rainfall and irrigation. The total infiltration amounts of the unsaturated zone under the four scenarios are almost the same, 1630 × 10 8 m 3 , while the leakage from the unsaturated zone to groundwater is about 0.45 × 10 8 m 3 . It indicates that water in the unsaturated zone is mainly lost by evapotranspiration and the majority of irrigation water is used to meet crop evapotranspiration. Since a large proportion of this

Evapotranspiration Lost from Unsaturated Zone and Groundwater
The water balance of unsaturated zone and saturated zone can be calculated by using the combined ETWatch-UZF-MODFLOW model. This was used for estimating the water loss of the study area in detail.
For the water balance in the unsaturated zone, the main input sources of water are rainfall and irrigation. The total infiltration amounts of the unsaturated zone under the four scenarios are almost the same, 1630 × 10 8 m 3 , while the leakage from the unsaturated zone to groundwater is about 0.45 × 10 8 m 3 . It indicates that water in the unsaturated zone is mainly lost by evapotranspiration and the majority of irrigation water is used to meet crop evapotranspiration. Since a large proportion of this freshwater for soil evaporation and crop transpiration comes from deep groundwater, i.e., evapotranspiration from agricultural cultivation leads to a negative deep groundwater balance, it cannot be fundamentally altered by current water conservation measures. As a result, the current mode of agricultural production will inevitably lead to a continuous decline in local groundwater.
The total evapotranspiration in the study area for the four scenarios during the simulation period is 1740 × 10 8 m 3 , 1710 × 10 8 m 3 , 1710 × 10 8 m 3 , 1720 × 10 8 m 3 , of which evapotranspiration lost from the unsaturated zone is 1620 × 10 8 m 3 , 1630 × 10 8 m 3 , 1620 × 10 8 m 3 , 1630 × 10 8 m 3 , and evapotranspiration from groundwater is 119 × 10 8 m 3 , 81.9 × 10 8 m 3 , 85.0 × 10 8 m 3 , and 92.3 × 10 8 m 3 , respectively. Evapotranspiration from groundwater accounted for 6.85%, 4.79%, 4.97%, and 5.37% of the total evapotranspiration, respectively. This shows that the proportion of evapotranspiration from groundwater is not high. This is mainly due to the relatively large span of groundwater depth in the evaluation area. In the west part of the study area, groundwater table also declined (lower than 5 m) due to the over-pumping of deep groundwater. This depth was larger than the extent depth of evaporation of groundwater (Wang et al., 2008). However, groundwater table in the east part was shallow, which accelerated the evaporation from groundwater. Additionally, soil salination was much deteriorated, where combined utilization of transferred water, salty water and deep fresh water is necessary.
To further investigate groundwater evapotranspiration in the west of the study area, an offshore subarea with a groundwater depth of less than 3 m was delineated in the study area. Under the irrigation conditions in Scenario 3, the total groundwater evapotranspiration (ET a ) in this subarea is 499 × 10 8 m 3 ; the evapotranspiration from the unsaturated zone (ET u ) and from groundwater (ET g ) are 440 × 10 8 m 3 and 595 × 10 8 m 3 , respectively. There was12% of ET a that was lost by groundwater ( Figure 12). The large percentage was caused by the local water table depth, cultivated acreage, planting structure, irrigation and rainfall. At the same time, these data can provide a reference basis for the construction of groundwater model and water balance calculation in the coastal region of the North China Plain. freshwater for soil evaporation and crop transpiration comes from deep groundwater, i.e., evapotranspiration from agricultural cultivation leads to a negative deep groundwater balance, it cannot be fundamentally altered by current water conservation measures. As a result, the current mode of agricultural production will inevitably lead to a continuous decline in local groundwater. The total evapotranspiration in the study area for the four scenarios during the simulation period is 1740 × 10 8 m 3 , 1710 × 10 8 m 3 , 1710 × 10 8 m 3 , 1720 × 10 8 m 3 , of which evapotranspiration lost from the unsaturated zone is 1620 × 10 8 m 3 , 1630 × 10 8 m 3 , 1620 × 10 8 m 3 , 1630 × 10 8 m 3 , and evapotranspiration from groundwater is 119 × 10 8 m 3 , 81.9 × 10 8 m 3 , 85.0 × 10 8 m 3 , and 92.3 × 10 8 m 3 , respectively. Evapotranspiration from groundwater accounted for 6.85%, 4.79%, 4.97%, and 5.37% of the total evapotranspiration, respectively. This shows that the proportion of evapotranspiration from groundwater is not high. This is mainly due to the relatively large span of groundwater depth in the evaluation area. In the west part of the study area, groundwater table also declined (lower than 5 m) due to the over-pumping of deep groundwater. This depth was larger than the extent depth of evaporation of groundwater (Wang et al., 2008). However, groundwater table in the east part was shallow, which accelerated the evaporation from groundwater. Additionally, soil salination was much deteriorated, where combined utilization of transferred water, salty water and deep fresh water is necessary.
To further investigate groundwater evapotranspiration in the west of the study area, an offshore subarea with a groundwater depth of less than 3 m was delineated in the study area. Under the irrigation conditions in Scenario 3, the total groundwater evapotranspiration (ETa) in this subarea is 499 × 10 8 m 3 ; the evapotranspiration from the unsaturated zone (ETu) and from groundwater (ETg) are 440 × 10 8 m 3 and 595 × 10 8 m 3 , respectively. There was12% of ETa that was lost by groundwater ( Figure 12). The large percentage was caused by the local water table depth, cultivated acreage, planting structure, irrigation and rainfall. At the same time, these data can provide a reference basis for the construction of groundwater model and water balance calculation in the coastal region of the North China Plain.

Optimal Agricultural Irrigation Ways Based on Analysis of Groundwater Balance
Results of groundwater balance was calculated and demonstrated in Figure 13. The equilibriums of the four scenarios are −327 × 10 8 m 3 , −287 × 10 8 m 3 , −272 × 10 8 m 3 , and −257 × 10 8 m 3 , respectively. All scenarios show negative water equilibrium during 2013-2030, suggesting the consuming of water resources in this region. This reflects the impact of the various water conservation measures on water balance. Comparing with the current ways of water utilization, Scenario 2, Scenario 3, and Scenario 4 can save 40.4 × 10 8 m 3 , 55 × 10 8 m 3 , and 69.9 × 10 8 m 3 of water over 18 years, respectively. Joint use of multiple water sources of Scenario 4 saves a maximum of 69.9 × 10 8 m 3 ; it is still smaller than the total deficit of −257 × 10 8 m 3 .Therefore, the trend of groundwater level decline continues due to deep groundwater abstraction. It shows that the total abstraction of groundwater for the four scenarios are 225 × 10 8 m 3 , 221 × resources in this region. This reflects the impact of the various water conservation measures on water balance. Comparing with the current ways of water utilization, Scenario 2, Scenario 3, and Scenario 4 can save 40.4 × 10 8 m 3 , 55 × 10 8 m 3 , and 69.9 × 10 8 m 3 of water over 18 years, respectively. Joint use of multiple water sources of Scenario 4 saves a maximum of 69.9 × 10 8 m 3 ; it is still smaller than the total deficit of −257 × 10 8 m 3 .Therefore, the trend of groundwater level decline continues due to deep groundwater abstraction. It shows that the total abstraction of groundwater for the four scenarios are 225 × 10 8 m 3 , 221 × 10 8 m 3 (including shallow groundwater pumping of 1.45 × 10 8 m 3 ), 203 × 10 8 m 3 , 180 × 10 8 m 3 (including shallow groundwater pumping of 1.45 × 10 8 m 3 ) (Figure 13).
Although the sources of irrigation water for the four scenarios vary greatly, the total amount of irrigation water is the same for all four scenarios under the current cultivated acreage and cropping structure (to ensure normal crop growth), and the similarity of the simulation results for the four scenarios indicates that the total amount of irrigation determines the overall trend of the local groundwater flow field (especially the shallow aquifer). From the simulation results, the difference in the flow field between the four scenarios is not very large, and the direction of groundwater flow is basically the same for all four scenarios. Scenario 4can slow down the rate of decline of deep groundwater and reduce the shallow groundwater table in the study area. Therefore, it was referred to as the optimal irrigation method to be used in the eastern NCP, even though the groundwater decline continues. Figure 13. Simulated changes of groundwater balance terms in the eastern NPC basin during the whole period(2013-2030).CHD IN is the recharge from constant head boundary, RIVERIN is the recharge from Rivers, GHB IN is the recharge from general head boundary, UZF RECHARGE is the recharge from unsaturated zone, GW ET is the groundwater evapotranspiration, CHDOUT is the groundwater discharge to constant head boundary, WELLS OUT is the groundwater exploitation, SURFACE LEAKAGE is the discharge to surface, RIVER LEAGAGE is the discharge to rivers, and ΔS is the storage change of groundwater system. PS:GW ET, CHD OUT,WELLS OUT,SURFACE LEAKAGE,RIVER LEAGAGE and ΔS are all negative. Figure 13. Simulated changes of groundwater balance terms in the eastern NPC basin during the whole period (2013-2030).CHD IN is the recharge from constant head boundary, RIVERIN is the recharge from Rivers, GHB IN is the recharge from general head boundary, UZF RECHARGE is the recharge from unsaturated zone, GW ET is the groundwater evapotranspiration, CHDOUT is the groundwater discharge to constant head boundary, WELLS OUT is the groundwater exploitation, SURFACE LEAKAGE is the discharge to surface, RIVER LEAGAGE is the discharge to rivers, and ∆S is the storage change of groundwater system. PS:GW ET, CHD OUT, WELLS OUT, SURFACE LEAKAGE, RIVER LEAGAGE and ∆S are all negative.
Although the sources of irrigation water for the four scenarios vary greatly, the total amount of irrigation water is the same for all four scenarios under the current cultivated acreage and cropping structure (to ensure normal crop growth), and the similarity of the simulation results for the four scenarios indicates that the total amount of irrigation determines the overall trend of the local groundwater flow field (especially the shallow aquifer). From the simulation results, the difference in the flow field between the four scenarios is not very large, and the direction of groundwater flow is basically the same for all four scenarios. Scenario 4 can slow down the rate of decline of deep groundwater and reduce the shallow groundwater table in the study area. Therefore, it was referred to as the optimal irrigation method to be used in the eastern NCP, even though the groundwater decline continues.

The Strengths and Limitations of the Coupled Model
We have developed an unsaturated-saturated zone coupled model; the important value of this model is its coupling of remote sensing data to the integrated hydrological model and remote sensing data, which has an irreplaceable advantage in reflecting the spatial heterogeneity of regional models, so the coupled model improves the accuracy of the hydrological model at aregional scale. The limitation of the model is that the PTF functions can only specify the parameters of the shallow unsaturated zone, which makes this model not effective in improving the accuracy of water level calculation in thick vadose zone area. The accuracy of the ETWatch data and the PTF function parameters we use are both at the kilometer level, so our model grid setting is also at the kilometer level. The data accuracy limited the refinement of the model.

Conclusions
In this paper, we developed an ETWatch-UZF-MODFLOW coupled model to predict the sustainability of groundwater under the current irrigation strategy (Scenario 1) in the eastern North China Plain. It also aimed to evaluate the effectiveness of three different combined multiple water use irrigation strategies (Scenarios 2, 3 and 4). During the construction of the model, remote sensing evapotranspiration data (ETWatch) were used instead of meteorological station scatter data, and the PTF function calculations were used to support the parameters required for the calculation of groundwater budget in the unsaturated zone. These data integration methods provide a new approach to the construction of integrated unsaturated zone-groundwater hydrological models at the regional scale and under irrigation conditions. The calculated results obtained from the ETWatch and PTF function can reflect the evapotranspiration on a regional scale and spatial heterogeneity of soil parameters. Meanwhile, the UZF package can distinguish ET g and ET u , avoiding the influence of ET u on groundwater balance calculations and further improving the groundwater balance calculations for shallow aquifers, and these calculations can significantly improve the fitting accuracy of the groundwater model.