Modeling contribution of shallow groundwater to evapotranspiration and yield of maize in an arid area

Capillary rise from shallow groundwater can decrease the need for irrigation water. However, simple techniques do not exist to quantify the contribution of capillary flux to crop water use. In this study we develop the Agricultural Water Productivity Model for Shallow Groundwater (AWPM-SG) for calculating capillary fluxes from shallow groundwater using readily available data. The model combines an analytical solution of upward flux from groundwater with the EPIC crop growth model. AWPM-SG was calibrated and validated with 2-year lysimetric experiment with maize. Predicted soil moisture, groundwater depth and leaf area index agreed with the observations. To investigate the response of model, various scenarios were run in which the irrigation amount and groundwater depth were varied. Simulations shows that at groundwater depth of 1 m capillary upward supplied 41% of the evapotranspiration. This reduced to 6% at groundwater depth of 2 m. The yield per unit water consumed (water productivity) was nearly constant for 2.3 kg/m3. The yield per unit water applied (irrigation water productivity) increased with decreasing irrigation water because capillary rise made up in part for the lack of irrigation water. Consequently, using AWPM-SG in irrigation scheduling will be beneficial to save more water in areas with shallow groundwater.

with EPIC crop growth model (Williams et al., 2006) to simulate the hydrological processes completely including the soil water flow, water flux at water table and crop growing. The WIPE-EPIC model 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), management (irrigation, growing time) and initial conditions (volumetric soil water content, groundwater depth). The AWPM-SG model needs less demanding data input soil parameters to simulate the groundwater capillary when compared with the SWAP-EPIC model by Xu et al. (2015). The schematic of the AWPM-SG is shown in Figure 2.

Description of crop module
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, Tmx, Tmn are the values of heat units, maximum temperature and minimum temperature in ºC for any day K, and Tb 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
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.
is the value of the maximum of crop leaf area index, 1 , 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 0 is the value of HUI when LAI starts declining.

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 ∆ is the variation of root depth in the i th day (cm); is the root depth in the i th day (cm); is the maximum depth reached by roots (cm).

Crop yield and water use efficiency
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. Ba 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 represent WUE, water use efficiency (kg/m 3 ), and IWP represent 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 ( Figure 1 and Supplementary

Description of soil module
The underground model in this study is a modification of the watershed irrigation potential estimation (WIPE) model designed by Saleh et al. (1989) 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, 1985) of the aquifer and is primarily applicable to shallow aquifers. Precipitation, irrigation, soil properties such as the moisture 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 S1 in the supplementary materials. 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.

Soil water flow
When RD ≤ , water balance is calculated in the zone 1 When mr * RD ≥ mf * RD where mf is field capacity ET = (S15) CAP = 0 (S16) = P + I + CAP + CR + mr * RD * 10 − mf * RD * 10 − ET (S17) (S20) Where 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), is actual evapotranspiration (mm); is the water depth supplied to the root zone from deeper zone due to the root growth (mm); is the water depths that leave the current root zone (mm); CAP is the capillary rise from zone 2 (mm) (Ritchie, 1996)， D 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.

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., 1989) 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 ), 2 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 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 ℎ , and φ is the matric potential calculated as The air entry value is calculated using (Saxton et al 1986) When 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 meter 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 is the extraction rate and dp is the drainable porosity.
Where is water content in the zone 2 (mm); mg is the soil moisture in the zone 2 (cm 3 /cm 3 ).
When RD > , there will be no zone 2, the calculation of zone 1 is similar to zone 2 when RD ≤ .
The input parameters of soil module (Figure 1 and Supplementary Table S3 in the Section 3 of the supplementary materials) 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.

Description of actual evapotranspiration module, ETa
The actual evapotranspiration is calculated and subtracted from soil-water storage.
is a fraction of potential evapotransiration, , which consists of potential evaporation from soil, , and potential transpiration from plants, Where mr is the moisture content of root zone and bt=3 for transpiration and be=0.8 for evaporation. Figure S1 Schematic of the soil profile used in the AWPM-SG Figure S2 Parameters sensitivity analysis for ET, groundwater depth, LAI and soil water content of top 90 cm. Figure  Note: A is the more sensitive parameters for ET (LAImx, ad, be, Tb, mf, ms, mwp), B is the less sensitive parameters for ET (kb, bt, T0, PHU, k2s, ks, dp, C, D0, b); C is the more sensitive parameters for groundwater depth (dp, mf, ms, Tb, PHU), D is the less sensitive parameters for groundwater depth (ks,k2s,D0,b,md,mwp,LAImx,ad,C,T0,kb,be,bt); E is the more sensitive parameters for LAI (LAImx, ad, Tb, T0, PHU), F is the less sensitive parameters for LAI (mf, md, ms, mwp, k2s, ks, dp, C, D0, b, kb, be, bt); G is the more sensitive parameters for soil water content (mf, ms, mwp, LAImx, Tb, T0, be, PHU), H is the less sensitive parameters for soil water content (md, k2s, ks, C, D0, b, dp, ad, kb, bt).    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; 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.