A unique vadose zone model for shallow aquifers: the Hetao irrigation district, China

Rapid population growth is increasing pressure on the world water resources. Agriculture will require crops to be grown with less water. This is especially the case for the closed Yellow River basin, necessitating a better understanding of the fate of irrigation water in the soil. In this paper, we report on a field experiment and develop a physically based model for the shallow groundwater in the Hetao irrigation district in Inner Mongolia, in the arid middle reaches of the Yellow River. Unlike other approaches, this model recognizes that field capacity is reached when the matric potential is equal to the height above the groundwater table and not by a limiting soil conductivity. The field experiment was carried out in 2016 and 2017. Daily moisture contents at five depths in the top 90 cm and groundwater table depths were measured in two fields with a corn crop. The data collected were used for model calibration and validation. The calibration and validation results show that the model-simulated soil moisture and groundwater depth fitted well. The model can be used in areas with shallow groundwater to optimize irrigation water use and minimize tailwater losses.

1 Introduction

With global climate change and increasing human population, much of the world is facing substantial water shortage (Alcamo et al., 2007). The water crisis has caused widespread concern among public governmental officials and scientists (Guo and Shen, 2016; Oki and Kanae, 2006). Years of rapid population growth have squeezed the world water resources. The available fresh water per capita decreased from 13 400 m³ in 1962 to 5900 m³ in 2014 (World Bank Group, 2019).

Water supply in China is especially stressed. When averaged over the whole country, available water per capita is at the water stress threshold of 1700 m³ yr⁻¹ (Falkenmark, 1989; Brown and Matlock, 2011). It is even less in the arid to semi-arid Yellow River basin that produces 33 % of the total agricultural production in China. To overcome water shortages in the Yellow River basin, crops are irrigated from surface water and groundwater. This irrigation has directly changed the hydrology of the basin. While 50 years ago, the semi-arid North China Plain had springs, shallow groundwater and rivers feeding the Yellow River, at present rivers and springs have dried up where groundwater is used for irrigation (Yang et al., 2015a). At the same time, in arid Inner Mongolia, along the Yellow River, the once-deep groundwater is now within 3 m of the soil surface in the large irrigation projects such as the Hetao irrigation district because of downward percolation of the excess irrigation water that has been applied.

In the Yellow River basin, crop irrigation accounts for 96 % of the total water use (Li et al., 2004). Due to the increased demand for irrigation, the river has stopped flowing downstream for an average of 70 d yr⁻¹ (Hinrichsen, 2002). Saving water upstream in Inner Mongolia by improved management practices mean that more water will be available downstream (Gao et al., 2015). In addition, the Hetao district is suffering from salinization, which leads to the land degradation (Guo et al., 2018; Huang et al., 2018). Salinization is caused by upward migration of water (and salt) from the shallow groundwater table that leads to salt accumulation at the surface (Ren et al., 2016; Yeh and Famiglietti, 2009). Designing improved management practices to save water and decrease salinization can be achieved by field trials or with the aid of a computer simulation mode measuring the fluxes. Field trials are time-consuming and expensive, and only a limited set of water management practices can be investigated. Models can test many management practices; however, the modeling results are often questionable because they have not been validated under local field condition and have not been validated for the future conditions. A combination of field experiments together with models has the benefits of both approaches with few negative effects.

Central to modeling irrigation management practices under shallow groundwater conditions (such as in the Yellow River basin) is simulating the soil moisture content accurately (Batalha et al., 2018, Gleeson et al., 2016; Jasechko and Taylor, 2015; Venkatesh et al., 2011a) because the moisture content plays a critical role in the growth of crops (Rodriguez-Iturbe, 2000), groundwater recharge (Hodnett and Bell, 1986) and upward movement of water to the root zone in areas (Gleeson et al., 2016; Jasechko and Taylor, 2015; Venkatesh et al., 2011a; Batalha et al., 2018). The last effect is unique to shallow groundwater areas where the moisture content and thus the unsaturated conductivity are high and where the drying of the surface soil sets up the hydraulic gradient that causes the upward capillary movement from the shallow groundwater (Kahlown et al., 2005; Liu et al., 2016; Luo and Sophocleous, 2010; Yeh and Famiglietti, 2009). The upward-moving water contains salt that is deposit in the root zone and at the surface.

There is tendency with the ever-increasing computer power to include all processes and the highly heterogeneous field conditions in hydrological models (Asher et al., 2015). In the case of simulating moisture contents, these models become complex and often fully distributed in three dimensions (Cui et al., 2017). Examples of these fully developed models are HYDRUS (Šimůnek et al., 1998), SWAP (Dam et al., 1997) and MODFLOW (Mcdonald and Harbaugh, 2003; Langevin et al., 2017). These models have long run times when applied to scenario simulations for real-world problems. In addition, calibration effort increases exponentially with the number of model parameters (Rosa et al., 2012; Flint et al., 2002). This makes the use of the complex models for real-time management and decision support cumbersome where many model runs are needed (Cui et al., 2017).

To overcome the disadvantages of the full and more complete models, computationally efficient surrogate models have been developed to speed up the modeling process without sacrificing accuracy or detail. Surrogate models are known under several names, such as metamodels, reduced models, model emulators, proxy models and response surfaces (e.g., Razavi et al., 2012a; Asher et al., 2015). We call the complex models “full” or comprehensive models.

Computational efficiency is the main reason for applying surrogate models in place of full models. Other advantages of surrogate models are shortening the time needed for calibration and identifying insensitive and irrelevant parameters in the full models (Young and Ratto, 2011). Most importantly, surrogate models allow investigating structural model uncertainty (Matott and Rabideau, 2008). Finally, surrogate models might be able to deal better with the self-organization of complex systems prevalent in hydrology than the full models (Hoang et al., 2017). For example, full models based on small-scale physics (Kirchner, 2006) cannot necessarily model the repetitive wetting patterns observed in humid watersheds, and for that reason, simple surrogate models often outperform their complex counterparts in predicting runoff when a perched water table is present in sloping terrains (Moges et al., 2017; Hoang et al., 2017).

Surrogate models can be classified in two categories (Todini, 2007; Asher et al., 2015): data-driven and physically derived models. Data-driven surrogates analyze relationships between the data available and physically derived surrogates simplifying the underlying physics or reduce numerical resolution. In recent years, most emphasis in the research literature has been on data-driven surrogate approaches (Razavi et al., 2012a). Relatively little research has been published on physically derived approaches. Despite its popularity, data-driven surrogates can be an inefficient and unreliable approach for optimizing complex field situations, especially when data are scarce, such as in groundwater systems (Razavi et al., 2012b). The physically derived surrogates overcome many of the limitations of data-driven approaches and are therefore superior over data-driven methods (Asher et al., 2015).

In the Yellow River basin various water-accounting models have been developed to simulate the soil water content and water fluxes (Xu et al., 2012; Chen et al., 2014; Xue and Ren, 2017; Yang et al., 2017; Ren et al., 2019). Numerical implementations are the finite-element model HYDRUS-1D by Ren et al. (2016) and Luo and Sophocleous (2010) and a finite-difference model by Moiwo et al. (2010). Surrogate models for the North China Plain, where the groundwater is more than 20 m deep, were published by Wang et al. (2001), Kendy et al. (2003), Chen et al. (2010), Ma et al. (2013), Yang et al. (2015, 2017a, b) and Li et al. (2017). In these models, the matric potential is ignored, and the hydraulic potential is equal to the gravity potential; thus the gradient of the hydraulic potential is unity (at least when it is expressed in head units). Under these conditions the water flux becomes negligible when the soil reaches field capacity at −33 kPa (equivalent to −3.3 m in head units), at which point the hydraulic conductivity becomes limiting. These models are not valid for irrigation projects along the Yellow River with shallow groundwater because the matric potential cannot be ignored over the short distance between the water table and the surface of the soil. Since the gravity and matric potential are of the same order, the water moves either down to the groundwater or up from the groundwater to the root zone, depending on the matric potential at the soil (Gardner, 1958; Gardener et al., 1970a, b). In summary, for shallow groundwater at less than 3.3 m from the surface, equilibrium is reached (i.e. fluxes are negligible) when the hydraulic gradient is zero (i.e., matric potential and gravity potential add up to constant value) and thus not when the conductivity becomes limited at a matric potential of −33 kPa.

For the irrigation perimeters with shallow groundwater in the Yellow River basin, we could find only two surrogate models developed by Xue et al. (2018) and Gao et al. (2017a, b). These two models do not consider the dynamics of groundwater depth and matric potential. By including these dynamics, more realistic predictions of moisture contents and upward flow can be obtained and would give better results when extended outside the area they are developed for (Wang and Smith, 2004). The reason is that for areas with shallow groundwater, evapotranspiration sets up a hydraulic gradient that causes the upward capillary water movement to sustain the evapotranspiration demands and crop water use (Kahlown et al., 2005; Liu et al., 2016; Luo and Sophocleous, 2010; Yeh and Famiglietti, 2009).

Advantages of physically driven surrogates are particularly relevant to groundwater studies where water tables are simulated over entire large areas, as shown by Brooks et al. (2007). Despite this, Asher et al. (2015) poses that physically driven methods have not been applied widely to groundwater problems, and even fewer have been applied with the interaction of moisture contents in the vadose zone, which is key in salinization and plant growth of the many cropped irrigated field in arid and semi-arid regions. In these water short areas it is extremely important to develop models that give directions on how to save water. The main objective of this study is, therefore, to develop a novel surrogate model and to validate this approach using experimental data collected in a field with shallow groundwater, where the ultimate goal is to save water in irrigation districts. In addition, sensitive and insensitive model parameters were identified for simulating moisture content in the shallow groundwater area to optimize future data collection efforts. The experimental fields are located in the Hetao irrigation district, Inner Mongolia, China, where in two maize fields, the moisture content and the groundwater table depth were measured over a 2-year period.

The surrogate model developed is a one-dimensional model simulating the moisture content in the root zone using the groundwater depth and information of the soil moisture characteristic curve. It can be easily adapted to the field scale by including the lateral movement of the regional groundwater. However, over short times, lateral movement can be neglected in nearly level areas outside a strip of 5–100 m from the river (Saleh et al., 1989), such as in deltas and lakes (Dam et al., 1997; Kendy et al., 2003).

2 Materials and methods

2.1 Study area

The Hetao irrigation district (HID) is the third-largest irrigation district of China. It covers an area of 1.12×10⁶ ha, of which half is irrigated (Xu et al., 2015). About 5 billion m³ of water is diverted from the Yellow River each year (Xu et al., 2010). The primary irrigation method used is surface flood irrigation (Sun et al., 2013). The groundwater table is very shallow, ranging between 0.5 and 3 m. The overall hydraulic gradient is 0.1 ‰–0.25 ‰ (Ren et al., 2018). Soil salinization is serious, and the chemical composition of groundwater salinity mainly consists of NaCl, K_cl and CaSO₄. The Hetao district has a typical arid continental climate, with high evaporation and low rainfall. The average annual precipitation is 180 mm, and the annual potential evapotranspiration is 2225 mm (Luan et al., 2018). The soil is mainly alluvial deposits with a silty loam texture. It is frozen 5 to 6 months per year, from late November to the middle of May. Maize and wheat are the main food crops, and sunflower is the main cash crop.

Figure 1 Location of the field experiment in Hetao irrigation district. The blue line is the Yellow River.

Table 1 Crop growth stage in 2016 and 2017 for corn growth on the Fenzidi experimental fields in the Hetao district.

2.2 Field experiment and data collection

The experiment was carried out in Fenzidi, Bayannur (41^∘9^′ N, 107^∘39^′ E), in the Hetao district in 2016 and 2017 (Fig. 1). In 2016, the experiment was carried out separately in site A (about 3100 m²) and site B (about 7000 m²; Fig. 1). In 2017, Field B was split into Field B1 and B2, and experiments were carried out in these two fields. Field B1 was about 3400 m², and B2 was about 3600 m². Experimental fields were planted both years with maize. The sowing dates were 24 April 2016 and 13 May 2017. The harvest date was 1 October in both 2016 and 2017. The plant growth stages are given in Table 1. The fields were flood irrigated three or four times during the heading and filling stages starting in late June or early July (Table 2).

Figure 2 Daily reference evapotranspiration, ET₀, and precipitation during the growing season in (a) 2016 and (b) 2017.

Table 2 Irrigation scheduling carried out at Fenzidi experimental fields in 2016 and 2017.

Precipitation, air temperature, relative humidity, sunshine duration and wind speed were collected from the weather station at the experimental station. The reference evapotranspiration (ET₀) was calculated based on the FAO Penman–Monteith equation with the daily meteorological data (Allen et al., 1998). Precipitation and ET₀ during the growing season are shown in Fig. 2. The soil moisture was monitored daily in the top 90 cm using HydraProbe soil sensors (Stevens Water Monitoring System Inc., Portland, OR, USA) installed in both experimental fields. Soil moisture was measured at five depths: 0–10, 10–30, 30–50, 50–70 and 70–90 cm. The sensors were connected to data loggers and downloaded via wireless transmission. Calibration was conducted by oven drying soil samples (Wang et al., 2018; Gao et al., 2017a). The groundwater depth was measured by piezometers (HOBO Water Level Logger U20, Onset, Cape Cod, MA, USA) recorded at 30 min intervals.

Soil samples were collected in rings from the same five layers where moisture contents were measured and used for determining soil physical properties including soil moisture at field capacity (θ_fc), soil moisture at saturation (θ_s), dry bulk density (ρ) and saturated hydraulic conductivity (K_s; Table 3). For Field A, B, B1 and B2, the saturated hydraulic conductivity was determined by the constant head method. Field capacity was determined at −33 kPa, and bulk density was determined by oven drying and dividing by the volume of the ring. Soil texture of Field A and B was analyzed with the laser particle size analyzer (Mastersizer 2000, Malvern Instruments Ltd., United Kingdom) in the laboratory and is shown in Table 4. The American soil texture classification was used in this study. The soils vary from silty loam to silty clay loam.

Table 3 Soil physical properties of the Fenzidi experimental fields.

Note: θ_fc is the soil water content at −33 kPa, θ_s is the saturated soil water content, K_s is the saturated hydraulic conductivity and ρ is the bulk density.

2.3 The shallow-aquifer–vadose zone model surrogate model

In developing the shallow-aquifer–vadose zone surrogate model for modeling moisture contents in the vadose zone, we followed the standards of good modeling practice by Jakeman et al. (2006). We made the model as simple as possible, provided justification for our surrogate technique, tested the surrogate model performance and finally provided detail on the method to encourage discussion on the technique that was followed.

2.3.1 Theoretical background

For shallow groundwater (less than 3.3 m deep), the matric potential is a function of depth under equilibrium conditions. Since the soil moisture characteristic curve for each soil is the relationship of moisture content and matric potential, the moisture content is also a function of the depth of the water table under equilibrium conditions.

Soil moisture characteristic curve

There are several formulations describing the soil moisture characteristic curve (Bauters et al., 2000; Brooks and Corey, 1964; Gupta and Larson, 1979; Haverkamp and Parlange, 1986; van Genuchten, 1980); the van Genuchten and Brooks–Corey models are widely used in hydrological and soil sciences. Here, we selected the Brooks–Corey model for its simplicity.

The Brooks–Corey model can be expressed as (Gardner et al., 1970a, b; Mccuen et al., 1981; Williams et al., 1983)

$$\begin{array}{}\text{(1a)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{e}}={\left({\displaystyle \frac{{\mathit{\phi }}_{\mathrm{m}}}{{\mathit{\phi }}_{\mathrm{b}}}}\right)}^{-\mathit{\lambda }}\text{for}\left|{\mathit{\phi }}_{\mathrm{m}}\right|>\left|{\mathit{\phi }}_{\mathrm{b}}\right|,\text{(1b)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{e}}=\mathrm{1}\text{for}\left|{\mathit{\phi }}_{\mathrm{m}}\right|\le \left|{\mathit{\phi }}_{\mathrm{b}}\right|,\end{array}$$

in which S_e is the effective saturation, φ_b is the bubbling pressure (cm), φ_m is matric potential (cm) and λ is the pore size distribution index. The effective saturation is defined as

$$\begin{array}{}\text{(2)}& S_{\mathrm{e}}={\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{d}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{d}}}},\end{array}$$

in which θ is the volumetric moisture content, θ_s is the volumetric saturated moisture content and θ_d is the residual air dry moisture content (all in cm³ cm⁻³). Equation (2) can be simplified to the form by setting θ_d=0:

$$\begin{array}{}\text{(3)}& S_{\mathrm{e}}={\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{s}}}}.\end{array}$$

For cases when the groundwater is close to the surface, under equilibrium conditions when the water flow is negligible (i.e., hydraulic potential is constant with depth), the matric potential can be expressed as height above the water table. For our field experiment the bubbling pressure, φ_b, and the pore size distribution index, λ, in the Brooks–Corey model can be obtained through a trial-and-error procedure by using the measured moisture content and matric potential derived from the groundwater depth after an irrigation event when the equilibrium state was reached and the sum of the gravity potential and matric potential was constant with depth.

2.3.2 Parameters based on soil moisture characteristic curve

The soil of the crop root zone is divided into several soil layers, and each soil layer has its specific soil moisture characteristic curve. After a sufficiently large irrigation and rainfall event, the moisture content is at equilibrium after the drainage stops. After such an event, the soil moisture of the vadose zone stays at the equilibrium moisture content as long as the evapotranspiration is less than upward flux from the groundwater.

Equilibrium moisture content

The equilibrium soil moisture content, θ_equ, in a layer can be determined by first replacing the matric potential in Eq. (1a) by the matric potential of the layer ${\mathit{\phi }}_{\mathrm{m}}^{z,h}$ that is dependent on the depth of the groundwater and depth of the soil layer, z, e.g.,

$$\begin{array}{}\text{(4)}& {\mathit{\phi }}_{\mathrm{m}}^{z,h}=h-z,\end{array}$$

where ${\mathit{\phi }}_{\mathrm{m}}^{z,h}$ is the matric potential under equilibrium moisture content at a depth z below the surface and h is the depth of the groundwater below the surface:

$$\begin{array}{}\text{(5a)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{eq}}^{z,h}={\mathit{\theta }}_{\mathrm{s}}^z{\left({\displaystyle \frac{h-z}{{\mathit{\phi }}_{\mathrm{b}}^z}}\right)}^{-\mathit{\lambda }}\text{for}⌈h-z⌉>\left|{\mathit{\phi }}_{\mathrm{b}}^z\right|,\text{(5b)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{eq}}^{z,h}={\mathit{\theta }}_{\mathrm{s}}^z\text{for}⌈h-z⌉\le ⌈{\mathit{\phi }}_{\mathrm{b}}^z⌉,\end{array}$$

where ${\mathit{\theta }}_{\mathrm{eq}}^{z,h}$ is the equilibrium soil moisture at the depth z below the surface, while the groundwater depth is h. Note that the superscripts z and h indicate the dependence on the distance from the soil surface, z, and the depth, h, of the groundwater table.

Drainable porosity

The drainable porosity, or specific yield, is defined as the amount of water drained from the soil for a unit decrease in the groundwater table when the soil moisture is at equilibrium. It is a crucial parameter in modeling the moisture content in our case or the amount of runoff for a shallow perched water table when there is rain (Brooks et al., 2007).

Figure 3 Illustration of drainable porosity for a soil moisture characteristic curve with a bubbling pressure of 40 cm. The yellow and the blue lines are the equilibrium moisture contents for the groundwater depth at 130 and 150 cm, respectively. The area between the two lines represents the amount of water for the decrease in groundwater table drained from the profile when the groundwater decreases from 130 to 150 cm.

By subtracting the total moisture content at equilibrium in the profile at the initial water table depth and at the new position one unit lower, we obtain the drainable porosity. For example, the area between the yellow and blue curve is the amount of water drained for a decrease in the water table from 130 to 150 cm (Fig. 3).

The total water content amount of the soil over a prescribed depth with a water table at depth h can be expressed as

$$\begin{array}{}\text{(6)}& W_{\mathrm{eq}}^h=\sum _{j=\mathrm{1}}^n{L}_j\stackrel{\mathrm{‾}}{{\left({\mathit{\theta }}_{\mathrm{eq}}^{z,h}\right)}_j},\end{array}$$

where $\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\mathrm{eq}}^{z,h}}$ is the average equilibrium moisture content of layer j for h taken at the midpoint of the layer, n is the number of layers in the profile and L_j is the height of soil layer j. The drainable porosity, μ^h, with the groundwater at depth h, can simply be found as

$$\begin{array}{}\text{(7)}& {\mathit{\mu }}^h={\displaystyle \frac{W_{\mathrm{eq}}^{h-\mathrm{\Delta }h}-W_{\mathrm{eq}}^{h+\mathrm{\Delta }h}}{\mathrm{2}\mathrm{\Delta }h}},\end{array}$$

where Δh=0.5L_j.

2.3.3 Calculating fluxes in the soil

The model accounts for the downward flux due to the irrigation and rainfall, evapotranspiration by plants and soil, and upward flux from the groundwater to satisfy some or all the evapotranspiration demand by the crop and soil. There are sets of rules implemented in an Excel spreadsheet to calculate the fluxes.

Evapotranspiration

The plant evapotranspiration was calculated in two steps. First the daily reference evapotranspiration (ET₀) was calculated by the Penman–Monteith equation (Allen et al., 1998). We assumed that the moisture content was limiting; therefore the plant evapotranspiration rate was obtained by multiplying the reference evapotranspiration by a crop coefficient (Allen et al., 1998; Sau et al., 2004; DeJonge et al., 2012). Values for the crop coefficients were calibrated according to the water balance in the soil and found to agree with published values for stage of crop development and soil salinity.

On days without rain or irrigation, the evapotranspiration lowers the water table, and the moisture content in the soil decreases due to upward movement of water to the plant roots and soil surface. On days with rain or irrigation, the potential evapotranspiration is subtracted from the irrigation and/or rainfall, and water moves downward.

Upward flux from groundwater

The upward flux from the groundwater, $U_{\mathrm{g}}^h$, is either limited by the potential evapotranspiration or the maximum flux of groundwater. The maximum flux, $U_{\text{g,max}}^h$, depends on the depth of the groundwater, the type of soil moisture characteristic curve and the condition at the surface (Gardner, 1958). These equations have an exponential form (Gardner, 1958; Yang et al., 2011; Zammouri, 2001),

$$\begin{array}{}\text{(8)}& U_{\text{g,max}}^h={\displaystyle \frac{a}{e^{bh}-\mathrm{1}}}\text{for}{U}_{\mathrm{g}}^h\le {\text{ET}}_{\mathrm{p}},\end{array}$$

where a and b are constants and ET_p is the potential evapotranspiration. The upward flux from the groundwater can be written as

$$\begin{array}{}\text{(9)}& U_{\mathrm{g}}^h=\min \left({\text{ET}}_{\mathrm{p}},U_{\text{g,max}}^h\right).\end{array}$$

On days without rain or irrigation, the soil moisture content is calculated by taking the difference of the equilibrium moisture content associated with the change in depth of groundwater. If in addition the upward flux is less than evapotranspiration, the difference between the upward flux and the evapotranspiration is extracted out of the root zone according to a predetermined distribution, r_j, e.g.,

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\stackrel{\mathrm{‾}}{\left({\mathit{\theta }}^{z,h,t}\right)}}_j={\left(\stackrel{\mathrm{‾}}{{\mathit{\theta }}^{z,h,t-\mathrm{\Delta }t}}\right)}_j+{\stackrel{\mathrm{‾}}{\left({\mathit{\theta }}_{\mathrm{eq}}^{z,h,t}\right)}}_j\\ & -{\stackrel{\mathrm{‾}}{\left({\mathit{\theta }}_{\mathrm{eq}}^{z,h,t-\mathrm{\Delta }t}\right)}}_j-{\displaystyle \frac{r_j\left(K_{\mathrm{c}}{\text{ET}}_{\mathrm{p}}-U_{\mathrm{g}}^h\right)}{L_j}},\end{array}\end{array}$$

where ${\stackrel{\mathrm{‾}}{\left({\mathit{\theta }}^{z,h,t}\right)}}_j$ is the average soil moisture content at time t of layer j, ${\stackrel{\mathrm{‾}}{\left({\mathit{\theta }}_{\mathrm{eq}}^{z,h,t}\right)}}_j$ is the average equilibrium soil moisture content of layer j when the groundwater depth is h at time t, K_c is a reduction factor of the potential evapotranspiration for saline soil water and the canopy, and r_j is the root function that determines that the portion of the evapotranspiration is taken up by the roots in layer j. The value z is taken at the midpoint of layer j. The time t is expressed in days and time, and t−Δt, is the previous day.

The downward flux

The rules for downward flux on days with the effective rain and/or irrigation are relatively simple. If the net flux at the surface (irrigation plus rainfall minus actual evapotranspiration) is greater than that needed to bring the soil up to equilibrium moisture content, the groundwater will be recharged, the distance to soil surface decreases and the moisture content will be equal to the equilibrium moisture content at the new depth. When the groundwater is not recharged, the following water balance is calculated: the rainfall and the irrigation are added to first layer. This layer will be brought up to the equilibrium moisture content, the remaining water fills up the next layer to the equilibrium moisture content and so on. The calculations can be expressed as follows:

$$\begin{array}{}\text{(11)}& {\stackrel{\mathrm{‾}}{\left({\mathit{\theta }}^{z,h,t}\right)}}_j=min\left[{\stackrel{\mathrm{‾}}{\left({\mathit{\theta }}_{\mathrm{eq}}^{z,h,t}\right)}}_j,{\left(\stackrel{\mathrm{‾}}{{\mathit{\theta }}^{z,h,t-\mathrm{\Delta }t}}\right)}_j+{\displaystyle \frac{R_{j-\mathrm{1}}\mathrm{\Delta }t}{L_j}}\right],\end{array}$$

where for j≥2, R_j−1 is the flux from the layer above and equals

$$\begin{array}{}\text{(12)}& R_{j+\mathrm{1}}=R_j-{\displaystyle \frac{\left({\stackrel{\mathrm{‾}}{\left({\mathit{\theta }}_{\mathrm{eq}}^{z,h,t}\right)}}_j-\left(\stackrel{\mathrm{‾}}{{\mathit{\theta }}^{z,h,t-\mathrm{\Delta }t}}\right)\right)L_j}{\mathrm{\Delta }t}}.\end{array}$$

For j=1, R₁ is equal to the rainfall plus the irrigation amounts minus potential evaporation.

Table 4 Soil texture of Field A and B.

Groundwater table depth

The groundwater in Hetao irrigation district has a small hydraulic gradient of 0.10 ‰–0.25 ‰ (Ren et al., 2016). In addition, the soil varies from silt loam to clay loam (Table 4) that has saturated hydraulic conductivity of less than 2 m d⁻¹. This means that the lateral fluxes are small compared the vertical fluxes and can therefore be neglected for the calculation of the groundwater depth. Based on this assumption, the net change in groundwater depth, Δh, can be calculated on days without rainfall or irrigation as

$$\begin{array}{}\text{(13a)}& \mathrm{\Delta }h={\displaystyle \frac{U_{\mathrm{g}}^h}{{\mathit{\mu }}^h}},\end{array}$$

and days with rain or irrigation as

$$\begin{array}{}\text{(13b)}& \mathrm{\Delta }h=-{\displaystyle \frac{R_{\mathrm{5}}}{{\mathit{\mu }}^h}},\end{array}$$

where the upward flux, $U_{\mathrm{g}}^h$, is calculated with Eq. (9), the percolation of the bottom layer, R₅, with Eq. (12) and the drainable porosity, μ^h, with Eq. (7). When the groundwater is close to the surface, the drainable porosity is zero. This would make the change in groundwater infinite. Thus, we limited the maximum decrease in groundwater after the irrigation event to 10–20 cm based on field observations.

2.3.4 Model calibration and validation

The soil moisture contents were measured from 30 May to 25 September in 2016 and 2017. Groundwater depth was observed from 13 June to 26 September in 2016 and 2017. For the convenience of simulation, the period of 13 June to 25 September was set as the simulation period. The model parameters were calibrated with the 2016 data and the validation with data collected in 2017 growing seasons. Soil moisture content of the top 90 cm (0–10, 10–30, 30–50, 50–70 and 70–90 cm) and the groundwater depth were simulated for model calibration and validation.

Relatively few parameters can be calibrated in the shallow-aquifer–vadose zone model. These are the crop coefficient, the K_c value, the two groundwater parameters and the root function. The other input data needed for model were the parameters in the Brooks–Corey equation (e.g., θ_s, θ_d, φ_b and λ) and were obtained by fitting the equation to the soil moisture characteristic curve of each layer of the soil. The saturated moisture content was measured independently as well and agreed with values obtained from the fit. Reference evapotranspiration was calculated directly from observed meteorological data.

For better understanding the model fitting performance, statistical indicators were used to evaluate the hydrological model goodness of fit (Ritter and Muñoz-Carpena, 2013). The statistical indicators including the mean relative error (MRE; Dawson et al., 2006), the root-mean-square error (RMSE; Abrahart and See, 2000; Bowden et al., 2002), the Nash–Sutcliffe efficiency coefficient (NSE; Nash and Suscliff, 1970), the regression coefficient (b; Xu et al., 2015), the determination coefficient (R²) and the regression slope (Krause et al., 2005) were used to qualify the model fitting performance during the model calibration and validation in this study. These statistical indicators can be expressed as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }\text{MRE}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\displaystyle \frac{(P_i-O_i)}{O_i}}\times \mathrm{100}\mathit{\%},\text{(15)}& {\displaystyle }& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N(P_i-O_i{)}^{\mathrm{2}}},\text{(16)}& {\displaystyle }& {\displaystyle }\text{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(P_i-O_i{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^N(O_i-\stackrel{\mathrm{‾}}{O}{)}^{\mathrm{2}}}},\text{(17)}& {\displaystyle }& {\displaystyle }b={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{O}_i\times P_i}{{\sum }_{i=\mathrm{1}}^N{O}_i^{\mathrm{2}}}},\text{(18)}& {\displaystyle }& {\displaystyle }{R}^{\mathrm{2}}={\left[{\displaystyle \frac{{\sum }_i^N(O_i-\stackrel{\mathrm{‾}}{O})(P_i-\stackrel{\mathrm{‾}}{P})}{{\left[{\sum }_i^N(O_i-\stackrel{\mathrm{‾}}{O})\right]}^{\mathrm{0.5}}{\left[{\sum }_{i=\mathrm{1}}^N(P_i-\stackrel{\mathrm{‾}}{P})\right]}^{\mathrm{0.5}}}}\right]}^{\mathrm{2}},\end{array}$$

where N is the total number of observations, O_i and P_i are the ith observed and predicted values ($i=\mathrm{1},\mathrm{2},\mathrm{\dots },N$), and $\stackrel{\mathrm{‾}}{O}$ and $\stackrel{\mathrm{‾}}{P}$ are the mean observed values and mean predicted values, respectively. For the MRE and RMSE, the values closest to 0 indicate good model predictions. NSE = 1.0 means a perfect fit, and the negative NSE values indicate that the mean observed value is a better predictor than the simulated value (Moriasi et al., 2007). For b and R², the values closest to 1 indicate good model predictions.

3 Results

In this section, we present first the 2016 and 2017 experimental observations of the Fenzidi experimental fields in the Hetao irrigation district (Fig. 1). This is followed by the calibration and validation of the shallow-aquifer–vadose zone model of moisture content in each of the five layers and the groundwater table depth.

Figure 4 Simulated and observed groundwater depth during the growing period for the Fenzidi experimental fields in the Hetao irrigation district: (a, b) calibration in 2016 and (c, d) validation in 2017. (Note: additional irrigation means the irrigation recharge from the adjacent field which leads to the water table rise and is not planned.)

3.1 Results of the field experiment

The total precipitation at the experimental during growing season was 62 mm in 2016 and 67 mm in 2017. The maximum daily rainfall was 23 mm in July 2017 (Fig. 2). The reference evapotranspiration varied between 1 to 5.5 mm d⁻¹, and the total ET₀ was 517 and 442 mm in the growing seasons during 2016 and 2017, respectively (Fig. 2). Daily observation consisted of groundwater depth (blue points; Fig. 4) and soil moisture content at five soil depths up to 90 cm (blue points; Fig. 5) for Field A and B in 2016 and Field B1 and B2 in 2017.

Figure 5 Simulated and observed soil moisture content for five soil depths during the growing period for the Fenzidi experimental fields in the Hetao irrigation district: (a, b) calibration in 2016 and (c, d) validation in 2017.

3.1.1 Groundwater observations

In 2016, the groundwater depth was generally more than 100 cm, except during the last two irrigation events in Field B, when it reached a depth of 72 cm for 1 or 2 d (Fig. 4). In 2017, groundwater tables were slightly closer to the surface than in 2016, especially in Field B2. The minimum groundwater depth was 61 cm on 21 June 2017 in Field B2 after an irrigation event.

In general, groundwater rose during an irrigation event and then decreased slowly due to upward movement of water to the plant roots to meet the transpiration demand. However, in the beginning of the growing season, we can see that the water table increased without an irrigation event. This occurred in Field A on 24 June 2016 and in Field B1 and B2 on 20 June 2017 (Fig. 4). This is curious and could be due to water originating from irrigation in a nearby field.

The water table at the end of the period of observation on 24 September 2016 is approximately 2 m deep, whereas on 15 June 2017, the depth decreased to around 125 cm. This is due to an irrigation application after the crops were harvested to leach the salt from the surface to deeper in the profile, bringing the water table up to near the surface. Evapotranspiration during the winter is small but sufficient for bringing the water table down. There was also a rainfall event on 5 June 2017 of 13 mm (Fig. 2) before the water table was measured, increasing the water level.

3.1.2 Soil moisture

Moisture contents are shown for the five layers and the two fields for 2016 and 2017 in Fig. 5. The moisture contents were near saturation when irrigation water was added and subsequently decreased (Fig. 5). For example, the soil moisture content changed in the 0–10 cm layer from 0.26 cm³ cm⁻³ to 0.42 cm³ cm⁻³ after the irrigation on 13 July 2016 in Field A and then gradually decreased to 0.34 cm³ cm⁻³. The moisture content decreased faster in the 10–30 cm depth than at any other depth for Field A, B and B1 but not for Field B2. The moisture content in Field A also showed a decrease at the 50–70 cm depth. For all plots, the moisture content at the 70–90 cm depth stayed nearly constant and only decreased during the growing season when the water table decreased below the 150 cm depth (Fig. 5). In Field A, the initial moisture content when the observation started was less than saturation; then after the first irrigation, it remained close to the saturated moisture content.

Table 5 Fitted Brooks–Corey parameters for the soil moisture characteristic curve.

It is interesting that while the soil profile was saturated (Fig. 4), the groundwater table was between 75–100 cm (Fig. 5). Before equilibrium moisture content was reached the water table was likely near the surface during the irrigation event. Because the drainable porosity was extremely small, even a minimum amount of evapotranspiration or drainage would cause the water table to decrease to roughly the height of the capillary fringe equal to the bubbling pressure, φ_b, in Eq. (5). The values of bubbling pressure are listed in Table 5.

Figure 6 Soil moisture characteristic curve of the four experiment fields for the Fenzidi experimental fields. The red line is the fit with the Brooks–Corey equation.

3.1.3 Soil moisture characteristic curve

In 2016 and 2017, the observed reduced moisture contents were plotted versus the height above the water table for the five soil layers of the two field sites in Fig. 6. These plots were used to define the soil moisture characteristic curves, which were of critical importance in simulating the moisture contents.

To define the soil moisture characteristic curve, the Brooks–Corey equation (Eq. 1) was fitted through the points closest to saturation at each matric potential representing the equilibrium conditions after an irrigation event. The fitted parameter values are shown in Table 5. Points to the left of the soil moisture characteristic curve are a result of evapotranspiration drying out the soil when the upward movement of water was insufficient for replenishing the moisture content in these layers, and thus matric potential and groundwater depth were not in equilibrium. In addition, the few points to the right indicate that the soil moisture was greater than the equilibrium moisture content. Many of the outlier soil moisture contents occurred in the layer from 0–10 cm, indicating that the soil was still draining after a rainfall event shortly before the measurements. Thus, the soil was not at the equilibrium moisture content.

The saturated moisture contents in Table 5 agree in general with those measured in Table 3 but are not exact. This is not a surprise, as the alluvial soil deposited by the rivers with layers varies over short distances. The variation within the field was also obvious from the soil's physical measurements. Field B1 and B2 are within Field B. The soil's physical properties of the various layers (Table 4) were not the same for the three sites, clearly showing the variability within the field.

Generally, large values of a pore size index coefficient λ are for sandy soils, and lower values are for clay soils (Bahmani and Bayram, 2018). We find this to be true for our site: for example, in Field A, λ=0.23 corresponds to a sandy layer with only 8 % clay in the 30–50 cm layer (Tables 4 and 5). In the 70–90 cm layer of Field B, the λ=0.07 corresponds to the clay layer of 23 % clay. In addition, bubbling pressure, φ_b, is greater for soils with a large clay content (Bahmani and Bayram, 2018). This is demonstrated for Field A in the 10–30 cm layer, where the bubbling pressure of 75 cm corresponded with the clay layer of 20 % clay. However, the correspondence between Tables 4 and 5 is not always perfect. This is especially obvious for the layer of 70–90 cm in Field A, where the values in Table 5 clearly indicate that the soil has a dense clay layer; however, the soil description in Table 4 shows that the soil is 39 % sand. This is due to the alluvial deposition patterns with changes in soil texture over short distances, as mentioned before.

3.2 Modeling results

The four parameters that can be calibrated in the shallow-aquifer–vadose zone model are the crop coefficient K_c value and the root function, both related to removal of water by the atmosphere, and the two groundwater parameters that determine the upward movement of water from the groundwater.

Table 6 Calibrated parameter values of the shallow-aquifer–vadose zone model.

3.2.1 Calibration of the parameters related to moisture content

The first step in the calibration was to fit the K_c value from the water balance. From the moisture contents and the groundwater depth, we can calculate approximately the amount of water lost to evapotranspiration. By comparing these values to the reference evapotranspiration calculated with the Penman–Monteith equation, we found that initially during the early stages the crop coefficient was 0.3 until the filling stage and then increased to 0.7 during the filling stage to the maturing stage (Table 6). These values are in accordance with the findings of Katerji et al. (2003) that salinity reduces the evapotranspiration (Katerji et al., 2003). According to the observed total salt content, the mean total salt content of experiment field in the 0–100 cm soil layer during the crop growth period was 2.29 g kg⁻¹ in Field A, 1.79 g kg⁻¹ in Field B, 2.33 g kg⁻¹ in Field B1 and 2.09 g kg⁻¹ in Field B2.

Table 7 Model statistics for calibration of the shallow aquifer model in 2016. The mean relative error (MRE), root-mean-square error (RMSE), regression slope, coefficient of determination (R²) and regression coefficient (b) are listed. Note: SWC is the soil water content, and GWD is the groundwater depth.

The second step was calibrating the moisture content by adapting the root function indicating the layers from which the water was taken up. Calibration was done manually by trial and error. We found that we could use the same root function for Field A, B, B1 and B2 (Table 6). The calibrated soil moisture contents of the five soil layers for the two fields in general are in agreement with the measured values in 2016 (Fig. 5a, b), with the coefficient of determination R² ranging between 0.48 to 0.94 with slopes of around 1, the MRE ranging between −9.38 % and 6.96 %, and the RMSE varying from 0.01 to 0.04 cm³ cm⁻³ for the five layers (Table 7). Finally, the parameters behaved physically and realistically, as water was extracted from shallow layers when the groundwater was close to the surface and from the deeper layers when the groundwater and the associated capillary fringe went down.

3.2.2 Validation of the parameters related to moisture content

The moisture contents predicted by the shallow-aquifer–vadose zone model were validated with the 2017 data in Field B1 and B2. Although the validation statistics of the five layers were slightly worse than for calibration in Table 7, the overall fit was still good, as shown in Fig. 5c and d. The coefficient of determination varied between 0.39 and 0.90. The MRE varied between −9.34 % and 19.48 %, and the mean RMSE range was from 0.01 to 0.07 cm³ cm⁻³ for the five soil layers (Table 7).

3.2.3 Calibration of the parameters related to groundwater depth

The final step was to calibrate the groundwater table coefficients with the 2016 data for both fields. We found that for fields not in the same location (e.g., A, B), the subsurface was sufficiently different; thus the same set of parameters could not be used (Table 6). The difference between the calibrated parameters for the two fields was small (Table 6). The measured and simulated groundwater depths were in good agreement with the chosen set of parameters (Fig. 4a, b), with the coefficient of determination R² being 0.67 for Field A and 0.85 for Field B (Table 7). Only from 15 July to 24 July did the observed water table on Field B decrease slower than the simulated water table. This is partly related to the fact that the properties of the soil below 90 cm were not measured, and the assumption was made that the soil moisture characteristic curve below 90 cm was the same as that from 70–90 cm. Thus the drainable porosity of the soil, which is a very sensitive parameter, might be different than what was used in the model. Another reason might be that the equation for upward movement might be too simple. Other statistical indicators show a good fit as well (Table 7).

3.2.4 Validation of the parameters related to groundwater depth

Since Field B1 and B2 are in the same location as Field B, we used the same set of groundwater parameters for the three fields (Table 6). The resulting fit between observed and predicted daily groundwater depths for Field B1 and B2 in 2017 was better than for the calibration in 2016 (Fig. 4c, d), with R² values of 0.84 for Field B1 and 0.86 for Field B2 (Table 7). In both cases, the slope of the regression line was close to 1. The other statistics indicated a good fit as well (Table 7), with the MRE being −0.05 for Field B1 and −0.02 for Field B2, the RMSE being 18.02 cm for Field B1 and 16.95 cm for Field B2, and the regression coefficient b being 0.94 and 1 for Field B1 and B2, respectively. The general agreement between the measured and simulated groundwater depth suggests that the two parameters are adequate, and the model can be used as a tool to simulate the change of the groundwater depth.

4 Discussion

In this paper, a novel surrogate model was developed for irrigation systems where the groundwater is close to the surface. The model uses the soil moisture characteristic curve to derive the drainable porosity and to predict the moisture contents in the soil. It is based on a definition of field capacity that is used less often (or equilibrium moisture content, as it is called in this paper) based on the observation that the flow becomes negligible when the hydraulic gradient is zero. In other words, the system is in equilibrium when the sum of the matric potential and the gravity potential is constant. Thus, when we chose the groundwater level as the reference point for the gravity potential, the matric potential is equal to the height above the groundwater. This is different from other applications of Darcy's law, where the groundwater is below 3.3 m. In these cases, groundwater movement stops when the conductivity becomes negligible at −33 kPa or 3.3 m in head units. The hydraulic conductivity value above −33 kPa (3.3 m in head units) does not limit the system, reaching equilibrium for daily time steps. No need therefore exists to measure this parameter in great detail for surrogate models. The opposite is true for the soil moisture characteristic curve for determining the spatial distribution of moisture content with depth above the groundwater.

In general, this surrogate model simulated the soil moisture content in each soil layer well, especially when compared to other models that attempted the soil moisture contents in the Yellow River basin such as the North China Plain (Kendy et al., 2003) and the Hetao irrigation district in Gao et al. (2017b) during the crop growth period. Our simulation results suggest that the reduction factor of the potential evaporation for soil saline K_c and root function parameters, together with the information of the soil moisture characteristic curves, can be used to adequately predict the soil moisture content. To predict the groundwater depth, two additional parameters are needed for the exponential function that defines the upward movement of groundwater.

The simulations, together with the observed data, indicated that information about the soil is very important to obtain the exact moisture content in the soil. However, generalized soil moisture characteristic curves for each soil type can be used in the simulation and will not result in great differences in water use by plants, since percolation to deeper layers was negligible, and thus the only loss of water was by evapotranspiration independent of the soil moisture content.

Finally, in the simulations we did not consider the influence of crop type and the influence of crop growth on soil moisture and groundwater depth. It would be of interest to investigate in future work whether the simulations would be improved by considering the dynamic crop characteristics during the growing season (Singh et al., 2018; Talebizadeh et al., 2018). A mature crop model, such as the EPIC model (Williams et al., 1989), which needs relatively few parameters, will certainly help to predict the crop yield but might not change the water use predictions. Actually, the EPIC model was already applied to the Hetao irrigation district by many researchers to analyze the crop growth during the crop growth period (Jia et al., 2012; Xu et al., 2015).

5 Conclusions

A novel surrogate vadose zone model for an irrigated area with a shallow aquifer was developed to simulate the fluctuation of groundwater depth and soil moisture during the crop growth stage in the shallow groundwater district. To validate and calibrate the surrogate model we carried out a 2-year field experiment in the Hetao irrigation district in upper Mongolia with groundwater close to the surface. Using meteorological data and the soil moisture characteristic curve and upward capillary movement, the surrogate model predicted the soil water content with depth and groundwater height on a daily time step with acceptable accuracy during validation and was an improvement of two previous models applied in the Hetao district that could predict the overall water content in the root zone but not the distribution with depth.

The surrogate modeling results show that after an irrigation event, as long as the upward flux from the groundwater to the root zone was greater than the plant evapotranspiration rate, the moisture contents in the vadose zone could be found directly from the soil moisture characteristic curve by equating the depth to the groundwater with the absolute value of the matric potential. When the plant evapotranspiration rate exceeded the upward movement, moisture contents would be indicated by groundwater depth and were predicted by a root zone function. Another finding was that the daily moisture contents were simulated without using the unsaturated hydraulic conductivity function in the surrogate model. For a daily time step, equilibrium (defined as the hydraulic potential being constant) in moisture contents in the profile was attained; thus precise unsaturated conductivity was not needed. Of course, for shorter time steps, for predicting the transient fluxes and groundwater, the conductivity function is needed. For management purposes a daily time step is acceptable.

Future improvement to this model will focus on coupling the EPIC model and applying it to simulate other crops and other locations with a shallow groundwater table. The surrogate model should also be compared with a “full” model to test the conditions under which the surrogate model will fall short.