Understanding the Role of Shallow Groundwater in Improving Field Water Productivity in Arid Areas

Soil water and salt transport in soil profiles and capillary rise from shallow groundwater are significant seasonal responses that help determine irrigation schedules and agricultural development in arid areas. In this study the Agricultural Water Productivity Model for Shallow Groundwater (AWPM-SG) was modified by adding a soil salinity simulation to precisely describe the soil water and salt cycle, calculating capillary fluxes from shallow groundwater using readily available data, and simulating the effect of soil salinity on crop growth. The model combines an analytical solution of upward flux from groundwater using the Environmental Policy Integrated Climate (EPIC) crop growth model. The modified AWPM-SG was calibrated and validated with a maize field experiment run in 2016 in which predicted soil moisture, soil salinity, groundwater depth, and leaf area index were in agreement with the observations. To investigate the response of the model, various scenarios with varying groundwater depth and groundwater salinity were run. The inhibition of groundwater salinity on crop yield was slightly less than that on crop water use, while the water consumption of maize with a groundwater depth of 1 m is 3% less than that of 2 m, and the yield of maize with groundwater depth of 1 m is only 1% less than that of 2 m, under the groundwater salinity of 2.0 g/L. At the same groundwater depth, the higher the salinity, the greater the corn water productivity, and the smaller the corn irrigation water productivity. Consequently, using modified AWPM-SG in irrigation scheduling will be beneficial to save more water in areas with shallow groundwater.


Introduction
In many semi-arid and arid regions of the world, irrigated agriculture consumes most of the available water. For example, in China, agricultural water use accounts for 60% of the total water use and 90% of the agricultural water is used for irrigation [1,2]. Therefore, developing water-saving agriculture and enhancing water productivity are of great significance for ensuring water security, food security, and ecological security in China [3].
In particular, in the Hetao irrigation district, which is a typical arid area with shallow groundwater and affected by soil salinization, various water-saving measures including lining of main, sub-main and distributor canals have been implemented since 2000 [4,5]. The application of water saving measures has led to a decline in the groundwater table that controls the waterlogging and salinity. Soil salinization is a serious environmental problem [6][7][8]. Healthy soils and healthy land are the basic conditions for the successful implementation and realization of the sustainable development goals related to food, healthy water and climate [9][10][11][12][13]. Therefore, understanding water and salt balance under shallow groundwater is significant for sustainable agriculture.
In districts with shallow groundwater, such as the Hetao Irrigation district, capillary rise from groundwater can be used to supplement surface water irrigation and in closed basins, it can possibly save water for irrigating additional areas. At the same time, the soil acts a filter for groundwater [14]. Then capillary rise from groundwater depends on several factors such as groundwater table depth, soil hydraulic properties, crop growth stage, weather and irrigation [15,16]. According to the entropy model, land-use, soil order and rainfall factors had the highest impact on groundwater potential [17]. Several studies have quantified the amount of water derived from groundwater and found that 20%-40% of the evapotranspiration can be met by capillary rise from water tables at depths of 0.70-1.50 m [18][19][20][21][22][23][24].
With lysmeter experiments, Kahlown and Ashraf [25] found that under the water table at 0.5 m depth, wheat met its entire water requirement from the groundwater and the sunflower's required water absorbed from groundwater is greater than 80%.
Soil salinity is a limiting factor for crop growth, especially within shallow groundwater. Pan et al. used a Soil Water Atmosphere Plant (SWAP) model to simulate soil water-salt balance, summer maize yield and water use efficiency under different irrigation schedules [26]. A modified HYDRUS model was used in simulating polyfluoroalkyl substances (PFAS) transport in the vadose zone [27]. By analysing satellite-based remote sensing images, Yu et al. [8] indicated that groundwater depth was the major controlling factor for regional soil salinity. Metternicht et al. [28] adopted remote sensing to analyze the spatial distribution and temporal changes in soil salinity.
Yield, water productivity (WP), and irrigation water productivity (IWP) vary with both groundwater depth and soil salinity. Mejia et al. [29] conducted a two-year study to investigate the effect of the groundwater table on the corn yield and found a 5%-10% greater yield for corn, and 23% for soybeans at 0.5-0.75 m groundwater depth compared with the same treatment with groundwater depth more than 1 m. Wang et al. [30] investigated the water and salt transport features for saline soil through drip irrigation under film, and found that the increase in irrigation water would expand the standard district for the normal growth of plants. Huo et al. [31] reported that irrigation water productivity (IWP) increased with shallower groundwater tables.
Although experimental results have indicated that shallow groundwater can contribute to crop water use, it is seldom used in irrigation scheduling because measuring the capillary upward flux is complex and cannot be performed routinely. Therefore, models are the only way to obtain fluxes. However, most of the irrigation management models assume that the groundwater is sufficiently deep that the water percolates downwards out of the plant roots [32][33][34]. Several methods that can estimate the upward water such as Hydrus and SWAP are cumbersome because they require spatially varying input data such as the soil water retention curve and unsaturated hydraulic conductivity that are not easily available [35]. Prathapar and Qureshi [36] used the SWAP model to evaluate consequences of deficit irrigation in semi-arid areas with a shallow water table in Pakistan. Other models, such as by Healy and Cook [37] and Ayars et al. [38] are based on measured moisture content or water-table fluctuations and cannot be used in routine applications. Additionally, common numerical methods require a lot of verification and depend on soil parameters.
Because of the intensive data requirements of the above-mentioned models for use in irrigation management with shallow groundwater, the AWPM-SG was established in previous studies to simulate the water cycle in shallow groundwater areas. The objective of this study is, therefore, to (1) develop a modified model that couples soil water and salt and crop growth to analyze the field water and salt cycle and groundwater contribution to crop evapotranspiration; and (2) investigate the effect of groundwater depth and groundwater salinity on yield, crop evapotranspiration, water productivity and irrigation water productivity of maize. Without affecting the maize yield, the irrigation water can be reduced within a reasonable range.

Experiments
For model calibration and validation, field experiments by Wang in 2016 in Fenzidi located in the Hetao irrigation district were used. The experimental station was located in an open field and had a representative climate, soil, and groundwater conditions. The Hetao irrigation district, is located in the western part of the Inner Mongolia Autonomous Region. The study site has a typically arid and semi-arid continental climate. The average annual rainfall is 142 mm and falls mainly from June to August. The average annual pan evaporation is 2300 mm. The average number of sunlight hours per month is 266 h. The mean annual temperature is 7°C, with monthly averages of −10.1 in January to 23.8°C in July ( Figure A1 in the Appendix B). Soils begin to freeze in the second half of November to a maximum depth of approximately 1.0-1.3 m, and is completely thawed by the middle of May. Groundwater depth varies between 1.2 and 3.8 m.
The experiment growing maize was carried out by Wang in 2016 [39]. The maize (Kehe-8) commonly grown in Hetao was planted in late April and harvested in late September. Cultivation practices were similar to the practices of farmers and recommended by extension agents. Three times during the growing season 216 and 218 mm of water were applied for the field plots named F1 and F2 ( Table 1). The irrigation amount refers to the local farmers' field. Nitrogen fertilizers were applied at a rate of 207 kg/ha for the first irrigation, and 103 kg/ha for the third irrigation in 2016. During the growing period, the meteorological data, including air temperature, sunshine hours, relative humidity and wind speed were taken from micrometeorological stations in farmland and were used to calculate the FAO (Food and Agricultural Organization of the United Nations) Penman-Monteith reference evapotranspiration (ET 0 ). Daily rainfall was measured using a rain gauge at the experimental site. The rainfall was 55.1 mm, with more than 70% of precipitation occurring in July and August. The main weather data are shown in Figure A1 in the Appendix B.
The particle size distribution was obtained using a laser particle analyzer. The main physical properties of the soil are shown in Table 2. In the root zone, the soil is silty loam, which is the representative soil texture in this region. The dry bulk density was obtained by oven drying undisturbed soil samples of 100 cm 3 at 105 • C for 48 h. Saturated hydraulic conductivity was estimated in eight 100 cm 3 undisturbed soil samples using a constant-head permeameter [40]. Soil water retention curves were determined for each horizon in 100 cm 3 undisturbed soil samples using a pressure membrane apparatus (SEC-1000, Soilmoisture Equipment Corporation, CA, US). The fitted soil hydraulic parameters of the Muchlem-van Genuchten model are presented in Table 3. Note: ms, md are the saturated moisture, residual moisture of 0-90 cm soil, respectively; ks is the saturated hydraulic conductivity; C, b are the content of soil in zone 1 and 2; D0 is the diffusion rate at wilting point; dp is the drainable porosity.

Data
The soil moisture content was measured every day using an Automatic Soil Moisture Monitor (Hydra Probe, Stevens Water Monitoring Systems, Inc., Portland, Oregon, USA). The soil moisture was measured at 20 cm intervals from the surface to a depth of 100 cm. At the same time, it was calibrated with drying method through earth-fetching at regular times.
Soil salinity was measured similarly to the soil water content. It was also measured every day using an automatic soil moisture monitor (Hydra Probe) at 20 cm intervals. Because the soil electrical conductivity (EC) monitored by Hydra Probe was affected by soil water and soil salinity, the soil electrical conductivity is directly related with soil water content. Two references that used Hydra Probe showed the results. In the study of Valdes et al., soil bulk EC (ECb) monitored by Hydra Probe is a function of both water and pore water EC (ECp). Then firstly we transfer the ECb to ECp based on soil water content [41]. Leao et al. reported that soil apparent bulk electrical conductivity is a function of soil water content, pore solution electrical conductivity and solid phase electrical conductivity [42]. We assumed the soil salinity is a function of soil water and soil electrical conductivity monitored by Hydra Probe. Because the monitored soil electrical conductivity needed to be modified according to the soil salinity obtained by earth-fetching, we take the regression as follows: where S is the soil total salinity (g/kg), EC is the soil electrical conductivity monitored by Hydra Probe (ms/cm), and SWC is soil water content (cm 3 /cm 3 ). Figure A2 in the Appendix B shows the contrast between the regression value and measured value. The groundwater depth was monitored daily during the crop growth period by monitoring the piezo metric head of the probe (Hobe 20).
The crop leaf area index (LAI) was measured every 6-12 days using a leaf area metre (YMJ-B). Dry maize yield was determined after harvesting.

Modified Agricultural Water Productivity Model for Shallow Groundwater
In this study, considering the high soil salinity in the Hetao irrigation area, the developed Agricultural Water Productivity Model for Shallow Groundwater (AWPM-SG) model was modified to simulate the actual soil water, soil salinity and groundwater depth. The developed AWPM-SG is the coupling of a crop growth Environmental Policy Integrated Climate (EPIC) model and a soil moisture model of the Watershed Irrigation Potential Estimation (WIPE) root and vadose zone that simulates both upward movement from groundwater and percolation to the groundwater. This was originally developed by Saleh et al. for investigating the impact of irrigation management schemes on groundwater levels in Bangladesh [43,44]. The structure of the modified AWPM-SG is shown in Figure 1. The model simulates water and salt flux on daily time steps and require few parameters (soil hydraulic parameters and crop growing parameters), which are listed in Tables 3 and 4 and. The three parts of the modified AWPM-SG consist of a crop module, an actual evapotranspiration module and a soil module, which are coupled for the first time in the AWPM-SG [15] (Figure A3 drawn by Xiaoyu Gao in the Appendix B). This modified model mainly considers soil salt transport and its effect on crop growth, especially in the shallow groundwater district (Method in Appendix A). The modified model is described in detail in the Appendix A. An overview of modified AWPM-SG is given below.  The crop module of AWPM-SG is mainly based on the Environmental Policy Integrated Climate (EPIC) model. It has been tested and applied widely around the world. Currently the EPIC model has evolved into a comprehensive model capable of simulating photosynthesis, evapotranspiration and other major plant and soil process [45]. The crop module includes the phonological development, crop growth indexes (LAI, biomass, root growing, crop yield) and water productivity (WP, IWP).

Description of Soil Module Considering Soil Salinity
Because EPIC cannot simulate the crop growth above the ground without the soil data, EPIC was coupled with the modified watershed irrigation potential estimation model (WIPE) to simulate the whole water recycle. The WIPE model was designed by Saleh et al. [44] to study the impact of irrigation management schemes on groundwater levels in Bangladesh. The model divides the soil profile into four zones namely the actual root zone, potential root zone, transmission zone, and the saturated zone as shown in Figure A3 in the Appendix B. The soil texture in zone 1 and 2 are same. The zone 1 is calculated using the water and salt balance method. The WIPE model simulates both downward recharge when the soil moisture in zone 2 is greater than field capacity and upward movement when the soil moisture is below field capacity using an analytical solution developed by Gardner [46]. This is a one-dimensional model employing the Thornthwaite-Mather procedure to calculate the recharge below the root zone to the aquifer and is primarily applicable to shallow aquifers. Full details of the model are provided in Appendix A (modified Agricultural Water Productivity Model for Shallow Groundwater).

Actual Evapotranspiration Module Considering Soil Salinity
Actual evapotranspiration is an input to the EPIC model. We used the model of Kendy et al. that previously was able to calculate the actual evapotranspiration from soil-water and salt storage in North China Plain with good accuracy [33]. In this model, the ratio of E p to T p depends upon the development stage of the leaf canopy, moisture content, soil salinity and root development. The input data are the daily leaf area index, LAI simulated by the EPIC model and simulated soil moisture and soil salinity of actual root zone (zone 1) by WIPE model and water balance.

Model Calibration and Validation
The simulation period for maize lasted from late April to late September in 2016 using the observed initial soil water content and groundwater depth subject to the imposed irrigation schedule and fertilizer applications. The S1 data were used for calibration and S2 to validate the model. Soil moisture and soil salinity in the top 90 cm (zones 1 and 2), groundwater depth, and LAI were simulated. The soil hydraulic parameters (md, ms, C, α and ks) and the crop parameters were calibrated (Tables 3 and 4). During the calibration, we set the parameters of the model according to the measured data and recommended values. We then analyzed both their sensitivity and uncertainty of parameters and found the sensitive parameters for soil water, groundwater depth and crop LAI, such as maximum leaf area index (LAI mx ), mf and so on. Then we adjusted the parameters as their sensitivity to make the simulation result of the model closer to the measured data. Finally, the calibrated parameters were used to validate the model using the S2 data. The default values of the EPIC model for maize were used as initial values for simulating crop growth [47]. The default value of the maximum rooting depth was 90 cm, as measured in the experiment by Wang [39]. Initial soil water content and groundwater depth were specified according to the measurements. A sensitivity and uncertainty analysis of parameters was performed once more for soil water content, soil salinity, groundwater depth and crop LAI.
The upper boundary condition was determined by the actual evaporation and transpiration rates, and the irrigation and precipitation were fluxed. A no-flux boundary condition was specified at the bottom of the column.
The mean relative error (MRE), root mean square error (RMSE), Nash and Sutcliffe model efficiency (NSE), coefficient of determination (R 2 ) and coefficient of regression (b) were used to quantify the model-fitting performance for both calibration and validation processes [48]. A MRE close to 0 indicates good model predictions. A RMSE value close to 0 indicates good model predictions. NSE = 1.0 represents a perfect fit, NSEclose to 0 represents the predicted values near the averaged measurement, and negative NSEvalues indicate that the mean observed value is a better predictor than the simulated value. An R 2 value close to 1 indicates good model predictions. A b value close to 1 indicates good model predictions [49].

Sensitivity Analysis of Parameters
By increasing or decreasing the percentage of the parameters, the variation of ET, groundwater depth, LAI, soil water and soil salinity were analyzed at 90 cm. Then the sensitivity of each parameter on ET, groundwater depth, LAI and soil water and salinity of 90 cm can be obtained ( Figure A4).

Uncertainty Analysis of Parameters
The d-factor was used to analyze the uncertainty of the parameters [50,51]. The d-factor is indicative of the average distance between the upper and lower confidence intervals (in this study the 95% prediction interval). The d-factor was calculated as follows: where d x is the average distance between the lower X Li and upper limits X Ui of the confidence interval and σ x is the standard deviation of the observed data. n is the number of data points. Larger d-factors lead to larger uncertainty.

Scenarios for Evaluation of Water Use and Yield
After calibration and validation of the modified AWPM-SG, the model was used to determine the response of groundwater depth (GWD) and groundwater salinity on water use, maize yield, WP and IWP. To do so the modified AWPM-SG was run for ten scenarios with 7 GWD levels (100, 150, 200, 250, 300, 350 and 400 cm) and 3 groundwater salinity levels (S, S1, S2-groundwater salinity of 1 g/L, 2 g/Land 4 g/L) were simulated. The rainfall and irrigation schedule for the simulation treatment was similar to that in the F1 plot used for calibration of the model. The initial soil water, groundwater depth soil texture and fertilizer application are set as the data of F1 plot. The meteorology data was based on the data in 2016.
The response of groundwater depth (GWD) and groundwater salinity on water use, maize yield, WP and IWP were investigated through running model with MATLAB. And then regression between WP, IWP and groundwater depth under different groundwater salinity was carried out using Origin 2017.

Evaluation of the Modified AWPM-SG
Data from the experimental study conducted in 2016 at S1 and S2 sites of the Fenzidi experimental station in the Hetao irrigation district, were used to calibrate and validate the modified AWPM-SG. The maize (Kehe-8) commonly grown in the Hetao irrigation district was planted in late April and harvested in late September. Experiments were carried out in two replicates.
The experimental data of F1 were used for calibration and F2 for validation. The soil data and crop data input are shown in Tables 3 and 4. The accuracy of the fit test is shown in Table 5. In the case where the observed data at F1 and F2 sites was remained nearly the same as the simulated values, while the data of validation treatment at F2 sites gave unrealistic results, which was due to the abnormal growth of crops at the F2 site. The soil water content of the top 90 cm (zones 1 and 2) was used to calibrate and validate the model. Simulated mean soil water content at the F1 site in the 90 cm soil zone was generally satisfactorily simulated with a mean relative error (MRE) and Nash and Sutcliffe model efficiency (NSE) of −1.07% and 0.49, respectively ( Figure 2A). The accuracy of fit was reasonable for most statistics presented in Table 5. At the S2 site the soil moisture was well predicted with the calibrated input data (Figure 2A, Table 5) including an NSE of 0.28 and R 2 of 1 ( Table 5). For calibration and validation, the averaged RMSE value was 2.11 cm.

Soil Salinity: Calibration and Validation
The soil salinity of the top 90 cm (zones 1 and 2) was used to calibrate and validate the model. Simulated mean soil salinity at the F1 site in the 90 cm soil zone was generally satisfactorily simulated with the MRE and NSE of −2.57% and 0.85 ( Figure 2B). The goodness of fit was reasonable for most statistics presented in Table 5. At the S2 site, the soil salinity was well predicted with the calibrated input data (Figure 2A, Table 5), including an NSE of 0.17 and b of 0.92 (Table 5). For calibration and validation, the averaged RMSE value was 0.3 g/L. Table: Calibration and Validation Figure 2C shows that the modified AWPM-SG simulations of the groundwater table followed the same trend as that observed during the calibration period, as indicated by the goodness of fit parameters with an R 2 of approximately 0.6, and RMSE value of 0.2 m (Table 5). For the validation period, the simulated data in S2 ( Figure 2C(a)) had a similar accuracy with the MRE of −2.63% (Table 5). Figure 2D show that the simulated LAI (calculated with Equations (A3) and (A6) in the Appendix A) followed the same trend as that observed during the calibration period. The average NSE and R 2 for LAI were 0.91 and 0.97, respectively, indicating a good fit ( Table 5). The average RMSE value for calibration and validation was 0.62 cm 2 /cm 2 . The measured LAI in F2 was less than that in F1, which can have many causes that were not included in the model such as the low temperature and snow in the seeding stage. LAI was therefore not as well predicted by the modified AWPM-SG for the validation period as for the calibration period with an MRE of 14.03% ( Figure 2D).

Parameter Sensitivity and Uncertainty Analysis
To investigate the effect of calibrated parameters on the crop evapotranspiration, groundwater depth, LAI, soil water content, crop yield and soil salinity of 90 cm soil, the parameters sensitivity and uncertainty were analyzed together ( Figure A4). When the parameters varied, the six indexes were linearly related to the change in the parameters, except for the Tb, PHU and C. For the crop evapotranspiration, the LAI mx affects the ET most obviously, and then the be and T b significantly affected the ET ( Figure A4A). The other parameters affected the ET slightly when ∆ET/ET was less than 5% as the parameters increased by 25%. Most of the parameters had little impact on the ET. The initial LAI mx is from the recommended value of the EPIC. Through the calibration and sensitivity of the parameters the final LAI was determined. Detailed information is shown in Figure A4 in the Appendix B.
For the groundwater depth, mf, dp, T b and PHU are the main impact factors. The others had almost no influence on the groundwater depth ( Figure A4B). mf was measured in the field experiment using undisturbed samples. The initial dp was obtained by the recommended value of the local government, and Tb, PHU was obtained by the recommended value of EPIC. Detailed information is shown in Figure A4 in the Appendix B.
Similar to the ET, the main impact factors on LAI are LAI mx , Tb, PHU, T0 and ad. The other parameters had little effect on the LAI. The initial Tb, PHU, T0 and ad are recommended by EPIC. Detailed information is shown in Figure A4 in Appendix B.
For the soil water content of 90 cm, parameters such as mf, ms, mwp, LAI mx , PHU, Tb affect the soil water content more obviously than the other parameters. mf, ms, and mwp were measured in the experiment. The initial LAI mx , PHU, and T b are recommended by EPIC. Detailed information is shown in Figure A4 in the Appendix B.
For the crop yield, the main impact factors are BE, T0, Tb, DLAI and LAI mx . The other parameters had little effect on the yield. The initial BE, T0, Tb, DLAI and LAI mx are recommended by EPIC. Detailed information is shown in Figure A4 in the Appendix B.
For a soil salinity of 90 cm, parameters such as mf, B, and ms affect the soil salinity more obviously than the other parameters. mf, ms were measured in the experiment. B is recommended by FAO-56. Detailed information is shown in Figure A4 in the Appendix B.
Combining the observed data, the uncertainty analysis of the parameters of groundwater depth, LAI, soil water content and soil salinity were analyzed ( Table 6). Note: In the uncertainty analysis, larger d-factor represents larger uncertainty of parameters. mf, ms, md, mwp are the field capacity, saturated moisture, residual moisture, wilting moisture of 0-90 cm soil, respectively; k2s is the saturated hydraulic conductivity in zone 1; ks is the saturated hydraulic conductivity in zone 2; C, b are the content of soil in zone 1 and 2; D0 is the diffusion rate at wilting point; dp is the drainable porosity; LAImx is the maximum of crop leaf area index; ad is a parameter that governs LAI decline rate for crop; DLAI is ratio that period from leaf area beginning to decline to whole growing period; Tb is the minimun temperature for plant growth, T0 is the optimal temperature for plant growth; kb dimensionless canopy extinction coefficient; be, bt are the parameters for evaporation and transpiration; PHU is total potential heat units required for crop maturation.
Through the uncertainty analysis of the parameters, the d-factors calculated by Equation (2) ranged from 0 to 1.62 with an average value of 0.17. A larger d-factor represents larger uncertainty. From the data, it can be inferred that for the groundwater depth, LAI, soil water content and soil salinity, the most uncertain parameters are the same with the most sensitive parameter. The most uncertain parameters are mf, LAI mx , mf, and mf for groundwater depth, LAI, soil water content and soil salinity, respectively. The trend of uncertainty of the parameters is basically similar to the sensitivity. Detailed information is shown in Table 6.

Groundwater Recharge on soil Water in Root Zone
The water flux at the lower boundary (F) of zone 2 is the exchange between groundwater and soil water, as shown in Figure 3. Obviously, groundwater recharge decreases with deeper groundwater, which has been proven in previous studies [15,16]. The results show that groundwater recharge to soil water decreases slightly with increasing groundwater salinity under full irrigation when groundwater depth is shallower than 3 m, while F with groundwater mineralization being 4 g/L is less than that with 1 g/L by 1.67 mm. When groundwater is deeper than 3 m, the effect of groundwater mineralization on F can be neglected, which is related to the relation between groundwater depth and F [31]. In general, regardless of which groundwater depth, the influence of groundwater salinity on groundwater rise is relatively small. Even when the groundwater depth is 1 m, the groundwater rise with a groundwater salinity of 4 g/L is only 2.7% less than that with the groundwater salinity of 1 g/L. In the case of other groundwater depths, the decline rate is smaller. Figure 3. Water fluxes at water table during the growing period with groundwater depth ranging from 1 m to 4 m and groundwater salinity ranging from 1 to 4 g/L. Note: S, S1 and S2 refer to the scenarios with groundwater salinity of 1, 2 and 4 g/L, respectively.

Soil Salt Content in Actual Root Zone
The variation of average salt in the actual root zone with groundwater depth and mineralization changes throughout the growing period is shown in Figure 4. The results show that with increasing groundwater salinity, the average salt content of the actual root zone increases gradually during the growth period under constant groundwater depth. This increase varies with different groundwater depths. For example, when the groundwater depths are 1 m and 2.5 m, the average salt contents of the actual root zone under groundwater salinity of 1, 2, and 4 g/L are 1.07, 1.15, 1.30 kg/m 2 and 1.06, 1.07, 1.10 kg/m 2 , respectively. The average salt content of the actual root zone increases only by 0.001 kg/m 2 under a groundwater depth of 4 m. Therefore, the groundwater depth directly affects the influence of groundwater salinity on the salt content of the actual root zone, which can be indicated in the above text that the groundwater capillary rise and groundwater salinity to the root zone increase with groundwater shallower.
In constrast, the average salt content of the actual root zone decreases with groundwater deeper under the same groundwater salinity, while the average salt content of the actual root zone ranges from 1.30 kg/m 2 to 1.07 kg/m 2 with groundwater depths varying from 1 m to 4 m under groundwater salinity of 4 g/L. Therefore, extensive study on groundwater depth and groundwater salinity is also an important part of a rational irrigation system [22].

Crop Evapotranspiration
The response of groundwater salinity to crop evapotranspiration under the same irrigation amount is shown in Figure 5A. On one hand, crop evapotranspiration decreases with increasing groundwater salinity under the same groundwater depth. The effect of groundwater salinity on crop evapotranspiration can be neglected when groundwater is deeper than 3 m. For example, crop evapotranspiration increases by 11.53 mm with groundwater salinity from 1 g/L to 4 g/L under groundwater depth of 1.5 m, but this increase is only 0.78 mm under groundwater depth of 3.0 m.  On the other hand, the groundwater recharge to soil water and crop evapotranspiration increases when groundwater is shallower, which indicates that the influence of groundwater depth on crop evapotranspiration is greater than that by average salt of the actual root zone. Xu et al. [48] reported that the groundwater contribution to crop growth was significant when the depth of the groundwater table was less than 1.50 m but was irrelevant for depths greater than 2.00 m. As shown in Figure 5, the crop evapotranspiration increases by 9 mm with groundwater depth ranging from 1.5 m to 1 m under groundwater salinity of 1 g/L, while the average salt of the actual root zone increases with groundwater shallower. With increasing groundwater salinity, the influence of groundwater depth on crop evapotranspiration is almost the same as that of the average salt of the actual root zone. For instance, crop evapotranspiration with groundwater depths of 1 m and 1.5 m under groundwater salinity of 4 g/L are similar.

Maize Yield
The effect of groundwater salinity on crop yield under the same irrigation amount is similar to that on crop evapotranspiration, as shown in Figure 5B. The influence of different groundwater salinities on maize yield mainly comes from two aspects. On one hand, crop yield decreases with increasing groundwater salinity under the same groundwater depth. When groundwater depth is greater than 3 m, this influence can be neglected, while maize yield was 13.07 Mg/ha, 13.06 Mg/ha and 13.04 Mg/ha with groundwater salinity of 1.0 g/L, 2.0 g/L and 4.0 g/L under groundwater depth of 3 m [25].
On the other hand, under same groundwater salinity, maize yield with shallower groundwater is greater than that with deeper groundwater. Maize yield is affected by groundwater depth greater than that by soil salt under low groundwater salinity. For example, maize yield with groundwater depth of 1.0 m under groundwater salinity of 1.0 g/L is greater by 0.1 Mg/ha than that with groundwater depth of 1.5 m, but this increase is only 0.02 Mg/ha under groundwater salinity of 4.0 g/L.
According to the effect of groundwater salinity on maize yield and crop evapotranspiration, the effect of groundwater salinity on the yield is slightly less than that on the crop evapotranspiration when the groundwater salinity is less. Based on the above analysis, when the salinity of groundwater is 1.0 g/L, the water consumption of maize under groundwater depth of 1.0 m increased by 3.6% compared to that under groundwater depth of 2.0 m, but the yield of maize under groundwater depth of 1 m was 1.6% larger than that of 2 m. Therefore, the effect of groundwater salinity on crop yield is slightly less than that on crop water consumption [52].

Water Productivity and Irrigation Water Productivity
Water productivity is the ratio of crop yield to evapotranspiration. The average water productivity is approximately 2.88 kg/m 3 , while water productivity was 2.83 kg/m 3 and 2.90 kg/m 3 under groundwater depths of 1 m and 3 m, respectively. The WP declined when the groundwater depth was greater than 3.5 m.
Irrigation water productivity (IWP) is the ratio of maize yield to the amount of irrigation water. Results show that the IWP of maize range from 4.0 kg/m 3 to 4.2 kg/m 3 . The IWP decreases with groundwater deeper, while IWP decreases from 4.16 kg/m 3 to 4.04 kg/m 3 with groundwater depth ranging from 1 m to 4 m under groundwater salinity of 1 g/L.

Contribution of Groundwater Recharge on Crop Water Use
The response of groundwater salinity to groundwater recharge to crop water use is shown in Figure 6. Similar to the findings of Luo and Sophocleous [23], ratios of seasonal groundwater evaporation to seasonal potential ET were plotted against the depth to the water table. When the groundwater was deeper than 3.0 m, the crop ET decreased significantly due to a lack of groundwater upward flux to crop growth under deeper groundwater. When groundwater depth ranges from 1 m to 2 m, the groundwater contribution to crop water use increases slightly with increasing groundwater salinity under the same groundwater depth. As shown in Figure 6, groundwater contributions to crop evapotranspiration with groundwater salinity of 1 g/L is lower than that with 4.0 g/L by 0.15% under groundwater depth of 1 m. This indicates that the influence of groundwater salinity on crop evapotranspiration is slightly greater than that on groundwater recharge. Then, the effect of groundwater depth on groundwater contribution is the important factor. Because of the lower effect of groundwater on soil water and crop evapotranspiration with deeper groundwater, the influence of groundwater salinity on soil water and crop evapotranspiration will decrease [25,53]. Babajimopoulos et al. found that the specific field conditions about 3.6 mm/day of the water in the root zone originated from the shallow water table amounting to about 18% of the water, which was transpired by the maize with the water table observed at a mean depth of 0.58 m below soil surface [54]. Figure 6. Groundwater contribution to water use (F/ET) during the growing period with groundwater depth ranging from 1 m to 4 m and groundwater salinity ranging from 1 g/L to 4 g/L. Note: S, S1 and S2 refer to the scenarios with groundwater salinity of 1 g/L, 2 g/L and 4 g/L, respectively.

Relationship between WP, IWP and Groundwater Salinity, Groundwater Depth
Relationship between WP and groundwater depth and salinity is shown in Figure 7A. The WP declined when the groundwater depth was greater than 3.5 m. However, WP increased mildly with increasing groundwater salinity, especially with shallow groundwater. For example, WP increased by 0.01 kg/m 3 with groundwater salinity ranging from 1.0 g/L to 4.0 g/L, which indicates that the effect of salt on crop evapotranspiration is more significant than that on crop yield. The result is similar to the study of Jiang et al. [52] which reported that water productivity irrigated by brackish water is slightly higher than that irrigated by fresh water.
Relationship between IWP and groundwater depth and salinity under same irrigation amounts is shown in Figure 7B. The effect of groundwater on IWP is different from that on WP. Results show that the IWP of maize decreases from 4.16 kg/m 3 to 4.04 kg/m 3 with groundwater depth ranging from 1 m to 4 m under groundwater salinity of 1 g/L, which is similar to the results of Sun et al. [55]. However, IWP increases gradually with deeper groundwater under groundwater salinity of 4 g/L, which is due to the decline of soil salt from groundwater decline with deeper groundwater. For example, IWP ranges from 4.02 kg/m 3 with groundwater depth varying from 1 m to 4 m. Therefore, IWP increases with deeper groundwater when groundwater salinity decreases to a certain extent.
In contrast, IWP decreases with increasing groundwater salinity at the same groundwater depth, especially for shallow groundwater. As shown in the Figure 7, when the groundwater depth is 1 m and 1.5 m, the IWP of maize is 4.16 kg/m 3 , 4.12 kg/m 3 , 4.02 kg/m 3 and 4.13 kg/m 3 , 4.10 kg/m 3 and 4.03 kg/m 3 under groundwater salinity conditions of 1.0 g/L, 2.0 g/L and 4.0 g/L, respectively. It is further proved that the shallower the groundwater depth, the groundwater salinity has a certain inhibition on the yield and has a greater impact on irrigation water productivity. The groundwater salinity is affected by irrigation water quality, soil prosperity, land use and management and so on [56,57]. Relationship between water productivity, irrigation water productivity and groundwater depth, groundwater salinity. Note: (A) Relationship between water productivity and groundwater depth, groundwater salinity. (B) Relationship between irrigation water productivity and groundwater depth, groundwater salinity. Note: S, S1 and S2 refer to the scenarios with groundwater salinity of 1, 2 and 4 g/L, respectively.
In addition, the change of WP caused by groundwater level fluctuation is higher than that by groundwater salinity. The WP increased by 0.01 kg/m 3 with groundwater salinity ranging from 1.0 g/L to 4.0 g/L under groundwater depth of 1.0 m, while the WP increased by 0.03 kg/m 3 with groundwater depth ranging from 1.0 m to 1.5 m under groundwater salinity of 1.0 g/L. However, the change of IWP caused by groundwater level fluctuation and groundwater salinity has not obvious difference. Above results can be attributed to that the effect of capillary rise on crop evapotranspiration is more significant than that on crop yield under shallow groundwater [31,58].

Conclusions
The influence of different groundwater salinity on soil salinity, groundwater recharge, crop water use, yield, water use efficiency and other parameters was investigated using the modified AWPM-SG. When the groundwater depth was more than 3 m, the groundwater recharge to soil water had a negative value. Under the same groundwater salinity, soil salinity in the root zone increased with groundwater shallower and maize growth was affected significantly. Under the same salinity, the smaller the groundwater depth, the smaller the water productivity and irrigation water productivity of maize.
Because groundwater salinity has little effect on groundwater recharge and great restraint on the crop water use, the groundwater salinity is higher under the same groundwater depth. Then the contribution of groundwater recharge to crop water consumption increased slowly. In conclusion, the inhibition of soil salt content on maize yield is slightly less than that on maize water consumption. Therefore, under the same groundwater depth, the higher the salinity, the greater the corn water productivity, and the smaller the corn irrigation water productivity. In summary, the effect of groundwater salinity on crop water productivity and irrigation water productivity is not significant. In the future, the modified AWPM-SG can be promoted for similar areas with shallow groundwater.
Author Contributions: Conceptualization, Z.H. and Z.Q.; methodology, X.G. and P.T.; software, X.G.; formal analysis, X.G. and P.T.; data curation, X.G. and S.Q.; writing-original draft preparation, X.G. and P.T.; writing-review and editing, Z.Q.; visualization, X.G.; project administration, X.G.; funding acquisition, Z.Q. and X.G. All authors have read and agreed to the published version of the manuscript. Acknowledgments: The authors would like to thank the editor and all the reviewers for their insightful comments and constructive suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.

A.1. Modified AWPM-SG Model Overview
AWPM-SG incorporates a new crop growth model on WIPE's original code. It combines the WIPE model (Saleh et al., [44]) with EPIC crop growth model (Williams et al., [47]) to simulate the hydrological processes completely including the soil water flow, water flux at water table and crop growing. In addition, the modified AWPM-SG considers the cycle of soil salinity. The modified AWPM-SG requires various inputs for soil (field capacity, permanent wilting point, residual moisture, saturated water content, saturated hydraulic conductivity), groundwater (drainable porosity), daily weather (minimum and maximum temperature, solar radiation, rainfall), crop (crop-specific base temperature, value of the maximum of crop leaf area index, parameter that governs LAI decline rate for crop, the value of HUI when LAI starts declining, harvest index, reduction in yield per increase in the electrical conductivity of saturation extract (ECe), electrical conductivity of the saturation extract at the threshold of ECe when crop yield first reduces below Ym, a yield response factor), management (irrigation, growing time) and initial conditions (volumetric soil water content, groundwater depth). The modified AWPM-SG needs less demanding data input soil parameters to simulate the groundwater capillary when compared with the SWAP-EPIC model by Xu et al. [48]. The schematic of the modified AWPM-SG is shown in Figure 2.

A.1.1. Description of Crop Module-EPIC Considering Soil Salinity
The EPIC plant growth model was developed to estimate soil productivity as affected by erosion throughout the U.S. The processes simulated include leaf-area index; conversion of biomass; above ground mass; economic yield; root growth and water use.

Phonological Development
In EPIC, phonological development of the crop is based on daily heat unit accumulation. It is computed using the equation: where HU, T mx , T mn are the values of heat units, maximum temperature and minimum temperature in • C for any day K, and T b is the crop-specific base temperature in • C. A heat unit index (HUI) ranging from 0 at planting to 1 at physiological maturity is calculated as follows: where HUI is the heat unit index for day i and PHU is the potential heat units required for maturity. The value of PHU may be provided by the user or calculated by the model from normal planting and harvest dates.

Potential Growth
Interception of solar radiation is estimated with a Beer's law equation (William et al., [43]): where PAR is intercepted photosynthetic active radiation in MJ/m 2 , RA is solar radiation in MJ/m 2 , LAI is the leaf area index, and subscript i indocates the day of the year. The constant, 0.5, is used to convert solar radiation to photosynthetically-active radiation. Experimental studies indicate that the extinction coefficient varies with foliage characteristics, sun angle, row spacing, row direction, and latitude. The value used in EPIC (0.65) is representative of crops with narrow row spacings. A somewhat smaller value (0.4-0.6) might be appropriate for tropical areas in which average sun angle is higher and for wide row spacings. Using Monteith's approach, potential increase in biomass for a day can be estimated with the following equation: where ∆B p,j is the daily potential increase in biomass in t/ha, BE is the crop parameter for converting energy to biomass in (kg/ha)/(MJ/m 2 ).

A.1.2. Leaf Area Index
Leaf area index (LAI) is simulated as function of heat units, crop stress and crop development stages. From emergence to the start of leaf decline, LAI is estimated with the equation: where HUF is the heat unit factor, REG is the value of the minimum crop stress factor. LAI mx is the value of the maximum of crop leaf area index, ab 1 , ab 2 are the crop parameters. From the start of leaf decline to the end of the growing season, LAI is estimated with the equation where ad is a parameter that governs LAI decline rate for crop and HUI 0 is the value of HUI when LAI starts declining.

A.1.3. Effect of Environmental Stress on Biomass Growth
When an environmental stress factors (soil water, temperature, nitrogen, Phosphorus and so on) is less than 1, actual crop biomass growth is calculated as follows: where ∆B i is the actual growth of crop biomass in the i-th day (t/hm 2 ); REG is the value of the minimum crop stress factor.
where WS is water stress factor; TS is temperature stress factor; SS is soil salinity stress factor; T p is the potential crop evapotranspiration (mm/d); T a is the actual crop evapotranspiration (mm/d); TG is the average daily temperature ( • C); T b is the base temperature ( • C); T 0 is the optimal temperature ( • C); EC e is mean electrical conductivity of the saturation extract for the root zone (ms/cm); EC ethreshold is electrical conductivity of the saturation extract at the threshold of EC e when crop yield first reduces below maximum expected crop yield (Ym (ms/cm)); Ky is a yield response factor; B is the reduction in yield per increase in EC e [%/(ms/cm)], when crop yield first reduces below Ym; EC 1:5 is electrical conductivity of soil solution with soil ratio water of 1:5.

A.1.4. Root Growth
Root length is simulated as function of the heat units and the maximum of root length. The root length will be maximum before the mature stage: where ∆RD i is the variation of root depth in the ith day (cm); RD i is the root depth in the ith day (cm); RD mx is the maximum depth reached by roots (cm).

A.1.5. Crop Yield and Water Productivity
The crop yield considering the water stress can be expressed as the equation: where WSYF is a parameter expressing the sensitivity of harvest index to drought, WS is the water stress factor. HI is the harvest index. B a is the crop biomass above ground and calculated by LAI and BE (a conversion factor of crop transferring the energy to biomass).
Water productivity and irrigation water productivity was computed using the following equation: where WP represents the water use efficiency (WUE, kg/m 3 ), and IWP represents irrigation water use efficiency (kg/m 3 ). Y denotes the final grain yield (kg/ha), ET is the total ET (mm), I is the irrigation amount from planting to harvest (mm). The input parameters of crop module (Table 4) include daily maximum, minimum and mean value of temperature, T max , T min and T mean ( • C); crop-specific base temperature, Tb ( • C); the potential heat units required for maturity, PHU; the maximum of crop leaf area index, LAI max ; the maximum of crop root, RD mx ; the crop parameters, ab1, ab2; the parameter that governs LAI decline rate for crop, ad; the value of HUI when LAI starts declining, HUI0; the parameter expressing the sensitivity of harvest index to drought, WYSF; the harvest index, HI ( Table 4). The output is daily leaf area index, LAI; crop height, h; root growth, R; crop yield, Y and water productivity, WP, IWP.

A.2. Description of Soil Module-WIPE Considering Soil Salinity
The underground model in this study is a modification of the watershed irrigation potential estimation (WIPE) model designed by Saleh et al. [44] to study the impact of irrigation management schemes on groundwater levels. This is a one-dimensional model employing the Thornthwaite-Mather procedure to calculate the recharge (Steenhuis and van Molen [59]) of the aquifer and is primarily applicable to shallow aquifers. Precipitation, irrigation, soil properties such as the moisture content, soil salinity content and the hydraulic conductivity and the initial groundwater level are required to run this model.
The model begins by dividing the soil profile into four zones namely the current root zone, future root zone, transmission zone, and the saturated zone over an impermeable bed as shown in Figure A3 in the Appendix B. The zone 1 is the zone occupied by roots; the zone 2 is the zone that is not currently occupied by the roots but will be so after their complete development; the zone 3 is the unsaturated transition zone below the root zone with lower boundary at water table and the thickness of this layer varies in time according to extraction/evaporation and recharge; the zone 4 is the saturated zone and is regarded as the water table. The zone 3 is always at the constant moisture content and is equal to the saturated moisture content minus the drainable porosity.

A.2.1. Soil Water Flow
When RD i ≤ RD mx , water balance is calculated in the zone 1: Otherwise : P w = 0 (A26) Wr i+1 = Wr i + P + I + CAP + CR − ET − P w (A29) where Wr is water content in the zone 1 (mm); mr, mg are the soil moisture in the zone 1, 2 (cm 3 /cm 3 ); P is precipitation (mm), I is irrigation (mm), ET is actual evapotranspiration (mm); CR is the water depth supplied to the root zone from deeper zone due to the root growth (mm); P w is the water depths that leave the current root zone (mm); CAP is the capillary rise from zone 2 (mm) (Ritchie, 1996), D r is the averaged diffusivity of zone 1 (cm 2 /day); D0 is the diffusivity at wilt pointing of zone 1 (cm 2 /day); b is the empirical parameter of soil. In the study, the soil texture of zone 1 and zone 2 are same.

A.2.2. Water Balance Calculation of Zone 2:
When mg ≥ mf (field capacity) where mf is the redistribution moisture content of root zone and the flux J is given by (Saleh et al., [44]) which is always directed downwards (mm). Here ms is the saturated moisture content of root zone (cm 3 /cm 3 ), md is the air dry moisture content of root zone (cm 3 /cm 3 ), k 2 s is the saturated hydraulic conductivity of root zone (mm/day), and C is a constant to 13. In this condition, there will be no upward evaporation flux from the aquifer so flux = −J. When mg < mf there will not be any downward flux so that J = 0. However, the upward evaporative flux from the aquifer will be non-zero and is a function of depth to water table from soil surface as given by Gardner [46]: where ks is the saturated hydraulic conductivity of transmission zone (mm/day), h is the depth to the water table (mm), αis the diffusivity coefficient which is the inverse of air entry ϕ h , and ϕ is the matric potential calculated as: The air entry value is calculated using (Saxton et al. [60]): When the water table is closer to soil surface the flux will be maximum and as water table goes down the flux will decrease. The limiting depth at which the flux becomes zero is approximately 4.5 m below ground level.
Irrigation using groundwater is simulated by extracting water from the aquifer and adding it to root zone. The water table depth is updated as: where ex t is the extraction rate and dp is the drainable porosity.
Wg i+1 = Wg i + P w − CAP − CR − J + u (A36) where Wg is water content in the zone 2 (mm); mg is the soil moisture in the zone 2 (cm 3 /cm 3 ). When RD i > RD mx , there will be no zone 2, the calculation of zone 1 is similar to zone 2 when RD i ≤ RD mx . The input parameters of soil module (Table 3) include Saturated moisture, ms; Air-dry moisture, md; Saturated hydraulic conductivity, Ks, constant C, b; the thickness of zone 1 and 2 (maximum root depth), the diffusivity at wilt pointing, D0; initial soil moisture and groundwater depth. The output is daily soil moisture and groundwater depth.

A.3. Description of Actual Evapotranspiration Considering Soil Salinity
The actual evapotranspiration is calculated and subtracted from soil-water storage. ET a is a fraction of potential evapotransiration, ET p , which consists of potential evaporation from soil, E p , and potential transpiration from plants, T p . The ratio of E p to T p depends upon the development stage of the leaf canopy, expressed as τ, the dimensionless fraction of incident beam radiation that penetrates the canopy (Campbell and Norman [61]): kb is the dimensionless canopy extinction coefficient, with a value of about 0.82 (Stockle [62]) and LAI is leaf-area index, daily values of which can be obtained by simulation of EPIC. Accordingly, ET p is allocated to: Actual evapotranspiration, ET (mm), can be limited by the availability of water and salt in the soil. Thus total actual evaporation and transpiration from soil are modeled as Kendy et al. [33] calculate the crop evapotranspiration using this method in the north of China: where mr is the moisture content of root zone and bt = 3 for transpiration and be = 0.8 for evaporation.