2.1. Experimental site
The experiment was conducted at the rangelands of the Isfahan University of Technology, central Iran. The study site is characterized by dry and hot summers and mild and humid winters. The cultivated field with an area of 1000 m2 is located in latitude 32° 43' N and longitude 51° 33' E, 1600 m above sea level. The mean annual temperature is 17.0°C, the mean annual rainfall is 116.9 mm and the mean annual relative humidity of air is 38%. The hottest and coldest months of the year are reported in July and January, respectively. The maximum and minimum absolute temperatures in the area are +48°C and –30°C, respectively (Soltani, 2004). The climate of this region according to Emberger’s climate classification is arid.
2.2. Soil analysis
Soil samples were collected from the 0–50 cm layer below the surface to measure chemical and physical properties. The soil EC and pH were determined in the saturated extract (Slavich and Petterson, 1993). Soil texture was measured by the hydrometer method and soil texture class was determined using the USDA (United States Department of Agriculture) soil texture triangle (Bouyoucos, 1962). The bulk density was determined using the core method (Grossman and Reinsch, 2002), and the porosity of soil samples (f) was calculated using bulk density and particle density (i.e., 2.65 Mg/m3). Field capacity (FC) and permanent wilt point (PWP) were determined using a pressure plate at matric potentials of –33 and –1500 kPa, respectively (Klute, 1986). Physical and chemical properties of the studied soil are presented in Table 1.
Table 1
Physical and chemical properties of the soil
Depth
|
pH
|
EC
|
Sand
|
Clay
|
Silt
|
FC
|
PWP
|
Bulk density
(Mg/m3)
|
Porosity
(m3/100 m3)
|
(cm)
|
(-)
|
(dS/m)
|
(kg/100kg)
|
0–50
|
8.81
|
2.00
|
48.70
|
21.80
|
29.50
|
13.50
|
8.40
|
1.36
|
51.00
|
Note: EC, electrical conductivity; FC, water content at field capacity; PWP, water content at permanent wilting point. |
2.3. Experimental design and treatments
The reservoirs were made with a height of 50 cm and an outer diameter of 11 cm (i.e., a volume of 4749.25 cm3) to irrigate the saplings. We used plaster molds to prepare the wet reservoirs, and then the reservoirs were cooked in a pottery kiln at 900°C. The planting pits, 50 cm×50 cm (width and depth) with a distance of 2 m, were drilled and the saplings together with reservoirs were placed in the pits. A lid was placed on the mouth of each reservoir to avoid evaporation loss. The dimensions of the reservoirs are presented in Table 2.
All the treatments were under the same environmental conditions. In this research, three permeability treatments for the reservoirs were considered. Various permeabilities in the reservoirs were attained by adding different amounts of wheat straw and manure to the reservoirs.
Table 2
Geometric properties of the reservoirs used in this study.
Parameters
|
Reservoir dimensions
|
Thickness (cm)
|
1
|
Height (cm)
|
50
|
Inner diameter (cm)
|
10
|
Outer diameter (cm)
|
11
|
Volume (cm3)
|
4749
|
Area (cm2)
|
1917
|
Tree of Heaven (Ailanthus altissima), Chinaberry (Melia azedarach), White mulberry (Morus alba) and Black locust (Robinia pseudoacacia) were planted in March 2016.
2.4. Determination of water depletion from the reservoirs
The percentage of water loss through the reservoir walls was determined by filling the reservoir installed in the soil with water up to its neck level. After a specific time, the reservoir was filled with the same water up to the same level to replenish the water loss. The percentage of water released per unit of time was calculated by Eq. 1 (Naik et al., 2008):
\(P=\frac{100v}{Vt}\) (1)
where P is the percentage of water loss per unit of time through the reservoir, V is the neck level capacity of the reservoir (cm3), v is the volume of water required to fill up the reservoir again to its neck level between two consecutive fillings (cm3), and t (h) is the time elapsed between two consecutive fillings which was 24 h for all the reservoir.
2.5. Measuring saturated hydraulic conductivity (Ks) of the reservoirs
A modified falling-head method was adapted to measure the saturated hydraulic conductivity (Ks) of the whole reservoirs using tap water. The reservoirs were first submerged, with the reservoir full of water, in a tap water bath for three days for saturation. After that, the reservoir, full of water, was submerged to its neck in a water bucket for Ks measurement. The water level in the bucket was kept constant using overflow. A manometer (length = 100 cm and diameter = 1cm) tube was inserted into a rubber stopper and fitted tightly to the mouth of the reservoir. A schematic diagram of the experimental setup is shown in Fig. 1. Changes in the water head, which is the height of the water level in the manometer tube above the free water surface of the bucket, were monitored with time. The falling-head equation for the calculation of Ks is given by Eq. 2 (Stein, 1990):
\(\text{ln}\left(\frac{h}{{h}_{0}}\right)=\frac{A{K}_{\text{s}}}{aL}t\) (2)
where h0 is initial height of water level in the manometer tube above the free water surface (cm), h is height of water level in the manometer tube at time t (cm), A is the surface area of reservoir (cm2), L is average wall thickness of the reservoir (m), Ks is saturated hydraulic conductivity (m/s), and t is cumulative time (s). A plot of ln(h/h0) versus time, t, gives a straight line. a is cross-sectional area of manometer tube (m2). The Ks of reservoirs can be calculated from the slope of the line if the other terms (i.e., L, A, and a) are known. The average thickness of the reservoir was 1 cm.
2.6. Determination and modeling of soil hydraulic properties
Soil hydraulic properties including soil water characteristic curve and saturated hydraulic conductivity were measured on the intact samples. Intact soil samples were collected from three locations in the field by pressing core samplers of diameter of 5.3 cm and height of 4.5 cm into the soil with minimum soil disturbance. Then, the soil samples were immersed and saturated in water for 48 h. The saturated soil samples were consecutively equilibrated at the matric suctions (h) of 0, 2, 5, 10, 20, 50, 100, 330, 500, 1000, 2000, 5000, 10000, and 15000 cm using a pressure plate (Gee and Ward, 2000). After applying the last pressure, the samples were dried in the oven and weighted. The gravimetric water content values were converted to volumetric water content (θ) by multiplying by the bulk density (ρb). The saturated hydraulic conductivity (Ks) of the soil was measured on the same samples using the constant-head method. The Ks was calculated using the Darcy (1856) equation:
\({K}_{\text{s}}=\frac{V\times L}{A\times t\times \varDelta H}\) (3)
where, Ks is the saturated hydraulic conductivity [LT−1], V is the volume of water outflow [L3], L is the soil column height [L], A is the cross-sectional area of soil [L2], ΔH is the hydraulic potential difference between the two ends of the sample [L] and t is the time difference from t1 to t2 [T].
Soil hydraulic properties were fitted to the van Genuchten-Mualem model (Eqs. 4-6): (Schaap et al, 2006)
(4)
\(K\left({S}_{e}\right)={K}_{s}{S}_{e}^{l}{\left[1+{(1-{S}_{1}^{1/m})}^{m}\right]}^{2}\) (5)
where
\({S}_{e}=\frac{\theta -{\theta }_{r}}{{\theta }_{s}-{\theta }_{r}}. m=1-\frac{1}{n}\) (6)
where θs is saturated water content [L3L−3], θr is residual water content [L3L−3], Se is effective saturation; K(Se) is unsaturated hydraulic conductivity [LT−1], and n and α are shape parameters, and l is the soil pore tortuosity parameter, which is estimated to be about 0.5 on average for many soils.
The Solver tool in MS Excel was used to fit the van Genuchten-Mualem model to the measured data of soil water characteristics curve and optimize its parameters. The optimized hydraulic parameters for the studied soil are presented in Table 3.
Table 3
Optimized van Genuchten-Mualem hydraulic parameters of the studied soil, and reservoir
|
Permeeability
|
Soil hydraulic parameters
|
|
Soil
|
-
|
\({\theta }_{\text{r}}\) (cm3/ cm3)
|
\({\theta }_{\text{s}}\) (cm3/ cm3)
|
n
|
\(\alpha\) (1/cm)\(\)
|
\({K}_{\text{s}}\) (cm/day)
|
l
|
0-10 cm
|
0.068
|
0.381
|
1.09
|
0.008
|
62
|
0.5
|
Reservoir
|
Low
|
0.068
|
0.381
|
1.09
|
0.008
|
0.06
|
0.5
|
Medium
|
0.068
|
0.381
|
1.09
|
0.008
|
0.08
|
0.5
|
High
|
0.068
|
0.381
|
1.09
|
0.008
|
0.2
|
0.5
|
θs: saturated water content, θr: residual water content, Se is effective saturation, n and α are shape parameters, Ks: saturated hydraulic conductivity, and l: soil pore tortuosity parameter.
2.7. Field measurements and numerical modeling of water flow
The soil water content was measured one, three, and six days after irrigation at horizontal distances of 5, 25, and 50 cm from the reservoir and at depths of 0-10, 10-25, 25-40, and 40-75 cm from the soil surface 2 years after establishment sapling with TDR apparatus (TRIME-FM model with 0.1% accuracy) (Fig. 2.).
The water content distribution in the soil was modeled using HYDRUS-2D/3D (Simunek et al., 1999). Assuming a homogeneous and isotropic soil, the governing equation for water flow is the 2D Richards Equation as follows:
\(\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial x}\left[K \left(h\right)\frac{\partial h}{\partial x}\right]+\frac{\partial }{\partial z}\left[K \left(h\right)\frac{\partial h}{\partial z}+K \left(h\right)\right]-S\) (7)
where, θ is volumetric water content [L3L−3], h is pressure head [L], t is time [T], x and z are horizontal and vertical space coordinates, respectively, and Ks is saturated hydraulic conductivity [LT−1]. The S [T−1] as a sink term indicates the rate of water uptake by the roots from the soil.HYDRUS-2D/3D uses the Galerkin finite-element method to solve the transport equations as explained in detail by Simunek et al. (1999).
The domain was defined by an axisymmetrical 2D geometry. The soil and reservoirs were defined as two soil materials. The reservoir were defined as variable head. In addition, boundary conditions on the soil surface were defined as atmospheric boundary condition to consider evaporation from the soil surface. The right and left sides of the modeled environment are no flux, because the range of effect of the reservoirs is up to this boundary. We simulated only the right side of the presumed symmetric profile. Thus, the boundaries of the finite-element mesh are rectangular except on the right edge near the upper right-hand corner where the reservoir is located (Fig. 3). The simulated environment in the model is a range of 100 cm in width and 150 cm in height. The dimensions of the overall mesh were 3 cm and around the reservoir was 1.6 cm.
2.8. Root distribution and water uptake parameters
Plant-root distribution influences soil water and salinity distributions in the root domain under micro-irrigation. Root distribution is quantified by \(\beta\)(r,z) (dimensionless) in 2D geometry as follows:
\(\beta \left(r.z\right)=\left[1-\frac{z}{{z}_{m}}\right]\left[(1-\frac{r}{{r}_{m}})\right]{e}^{-(\frac{{p}_{z}}{{z}_{m}}\left|{z}^{*}-z\right|+\frac{{p}_{r}}{{r}_{m}}\left|{r}^{*}-r\right|)}\) (8)
where z and r are the distances [L] from the plant origin in the directions of z (depth), and r (radius), respectively, and zm and rm are the maximum rooting lengths [L] from the plant origin in the z and r directions, respectively. In Eq. 8, \({r}^{\text{*}}\) and\({ z}^{\text{*}}\), and \({p}_{\text{z}}\) (-) and \({p}_{\text{r}}\) (-) are all empirical parameters controlling spatial root distribution. These empirical parameters were considered to provide zero root water uptake (RWU) beyond zm to account for asymmetrical RWU in the vertical and radial directions and to allow maximum RWU within 0 to zm (Vrugt et al., 2001b). Root distribution parameters by investigating the distribution of plant roots in the field, reviewing resources and with the help of HYDRUS 2D/3D software are shown in Table 4.
Table 4. Root distribution parameters and parameters of the Feddes model (i.e., threshold and critical soil water pressures, h) for the studied plants.
In the absence of osmotic stress, the actual rate of root water uptake, S(h), in any point of the simulation domain is computed according to the multilinear model proposed by Feddes et al. (1978). In the model, it is assumes that S(h) is proportional to the maximum (potential) root uptake rate occurring when water is not limiting plant transpiration:
S(h) = α(h).Sp (9)
where α(h) is a dimensionless prescribed function of soil water pressure head and Sp is the potential (maximal) water uptake rate by roots (cm3 cm−3 d−1).
Parameters of the Feddes model for water uptake by root are as follows: h1 is the pressure head below which plant roots begin to extract water from the soil, h2 is pressure head below which plant roots extract water at the maximum rate, h3.low is limiting pressure head below which plant roots are not able to extract water at the maximum rate at potential transpiration rate of 0.1 cm day−1, h3.high is limiting pressure head at potential transpiration rate of 0.5 cm day−1, h4 is pressure head below which plant roots stop to take up water, which is usually taken at wilting point, r2H is potential transpiration rate at which the limiting pressure head allows extracting water at the maximum rate and r2L is potential transpiration rate at which the limiting pressure head allows extracting water at the minimum rate. Parameters of root water uptake by reviewing the sources, considering the water requirement of the plants and with the help of HYDRUS 2D/3D software are presented in Table 4.
2.9. Evaluation of simulation model predictions
The agreement of HYDRUS-2D/3D software predictions with the measured data was quantified in terms of three statistical measures: the mean bias error (ME), the root square error (RMSE), and the coefficient of determination (R2) defined as follows (Willmott, 1982):
\(\text{M}\text{E}=\sum _{i=1}^{N}({P}_{i}-{O}_{i})/N\) (10)
\(\text{R}\text{M}\text{S}\text{E}=\sqrt{\sum _{i=1}^{N}{({P}_{i}-{O}_{i})}^{2}/N}\) (11)
\({R}^{2}=\frac{\sum _{i=1}^{N}{({P}_{i}-{O}_{i})}^{2}}{\sum _{i=1}^{N}{({O}_{i}-\overline{O})}^{2}}\) (12)
where N is the number of data points, \({P}_{\text{i}}\) is the ith predicted data point, \({O}_{\text{i}}\) is the ith observed data, and \(\overline{O}\) is the mean of observed data. ME can identify potential bias (i.e., underestimation and overestimation) in the predicted values, whereas RMSE and \({R}^{2}\) give overall measures of the goodness-of-fit.