Wetting Patterns in a Subsurface Irrigation System Using Reservoirs of Different Permeabilities: Experimental and HYDRUS-2D/3D Modeling

Precipitation in arid lands is low and often highly variable; therefore, efficient irrigation is the best approach for managing limited water supplies and irregular precipitation events for establishing seedlings. The HYDRUS-2D/3D software is a useful tool for simulating water content distribution in the design of subsurface irrigation systems, determining the proper locations of reservoir and plant, irrigation scheduling and maximizing water use efficiency. In this study, wetting patterns around the reservoirs with different permeabilities (low, medium, and high) were assessed in the field and simulated using the HYDRUS-2D/3D software. Different amounts of animal manure and wheat straw were mixed with clay fractions to produce the reservoirs with different permeabilities. The results showed that the highest soil water content was observed near the reservoir, and it decreased with distancing from the reservoir. In the high permeability treatment (with saturated hydraulic conductivity, Ks = 0.2 cm/day), a large volume of water was released for almost three days, so that more water was deep-percolated and lowered the soil water available to plant roots. In the low and medium permeability treatments (Ks = 0.08, and 0.06 cm/day, respectively), the water slowly leaked, and the maximum soil water content was observed in the 20–40 cm layer over time. The HYDRUS-2D/3D software was able to simulate water flow and soil water content during the growth period of plants, with a high correspondence between measured and simulated soil water contents (R2 = 0.90–0.95), which was higher in the lower discharges.


Introduction
Water scarcity in arid and semi-arid regions is a major concern for establishing plants. High-performance irrigation systems, such as surface and subsurface drip irrigation systems, are often recommended to overcome this problem and to intensely increase the water use efficiency over that of traditional ones (Kandelous and Šimůnek 2010). Subsurface irrigation with a ring-shaped emitter is one of the irrigation techniques in arid lands. Emitters are usually made from rubber, and ring-shaped emitters that have been successfully used in the irrigation design, and are purely empirical. Besides, the current design of the ringshaped emitter does not allow one to easily detect malfunctions because the emitter is fully covered with a permeable textile (Saefuddin et al. 2019).
Several empirical, analytical, and numerical models have been developed for simulating water flow, soil water content pattern, and wetting bulb dimensions in surface/subsurface irrigation systems. Due to the advances in the processing power of computers, and the availability of numerical models simulating water flow and solute transport in the soils, many researchers have become interested in using such models for evaluating water flow in soils under subsurface irrigation systems. Among different models/software, HYDRUS-2D/3D (Šimůnek et al. 2016) is one of the most widely used dynamic and physical-based software for simulating soil water dynamics. One of the advantages of this software is that its input parameters are closely related to soil physical properties, which could be measured either in-situ or in the lab (Karandish and Šimunek 2016). The applicability of HYDRUS-2D/3D to numerous plant/soil scenarios has been demonstrated experimentally, including surface drip (Mailhol et al. 2011;Li et al. 2015), and subsurface drip irrigation systems (Kandelous and Šimůnek 2010;Mailhol et al. 2011). Kandelous et al. (2011) calibrated HYDRUS-2D/3D to simulate water movement from a subsurface drip irrigation system by comparing simulated results with measured soil water contents in several field experiments. Kandelous and Šimůnek (2010) used the HYDRUS-2D/3D to evaluate laboratory and field data involving water movement in a clay loam soil from point-water sources buried at different depths and with different discharge rates. Kandelous and Šimůnek (2010) investigated the capability of HYDRUS-2D/3D to estimate the dimensions of the wetted zone under laboratory and field conditions. The RMSE at different locations varied in the ranges 0.011-0.045 cm 3 /cm 3 and 0.98-4.36 cm for the volumetric water content and wetting dimension, respectively. It was reported that the correspondence between simulations and observations was very good. Dastorani et al. (2010) evaluated the efficiency of surface and subsurface irrigations in dryland environments for planting pistachio trees. They found a considerable difference in efficiency between the two irrigation systems and a relatively higher preference for a subsurface system over the traditional surface method.
Precipitation in arid lands is low and often highly variable; therefore, efficient irrigation is the best approach for managing limited water supplies and irregular precipitation events for establishing seedlings. Knowledge of how water is distributed in the soil is essential for the proper design and management of subsurface irrigation systems and plays a major role in determining the distance between reservoirs and system pressure to supply the water required by the plant. Non-uniform patterns of water distribution around reservoirs have led to the appropriate location of reservoirs as a key factor in determining the effectiveness of irrigation planning based on soil water content. Uncertainty in locating the reservoir can be due to the high sensitivity of point measurements to changes in wetting and soil water distribution patterns and determining the appropriate location. The quantification of soil water condition is influenced by various factors such as climate, soil hydraulic characteristics, and the plant type and distribution of the root system (Soulis and Elmaloglou 2017). The efficiency of reservoir-based irrigation planning depends on the location of the reservoirs and the threshold for irrigation. The question is always in which part of the wetting range the reservoir should be installed. Many recommendations for locating water content measurements are experimental and qualitative that cannot be generalized to all conditions. Furthermore, few studies have been done on plant species that are cultivated in the natural areas, especially in the arid areas of Iran. The main objectives of this study were: 1) To compare the water distribution patterns around reservoirs made with different permeabilities and, select the best reservoir with appropriate permeability, and 2) To assess the efficiency of the HYDRUS-2D/3D software in simulating water content distribution in the subsurface irrigation. The results of the present study will be applicable to biological restoration in arid and semi-arid regions.

Experimental Site
The experiment was conducted at the campus of the Isfahan University of Technology, central Iran. The study site is natural rangeland characterized by dry and hot summers and mild and humid winters. The cultivated field with an area of 1000 m 2 is located at latitude 32° 43' N and longitude 51° 33' E, 1600 m above sea level. The mean annual temperature in the study area 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. The climate of this region according to Emberger's climate classifications is arid.

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/m 3 ). 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). The physical and chemical properties of the studied soil are presented in Table 1.

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 cm 3 ) 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.

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): where P is the percentage of water loss per unit of time through the reservoir, V is the neck level capacity of the reservoir (cm 3 ), v is the volume of water required to fill up the reservoir again to its neck level between two consecutive fillings (cm 3 ), and t (h) is the time elapsed between two consecutive fillings. t (h) was 24 h for all the reservoirs.

Measuring Saturated Hydraulic Conductivity (K s ) of the Reservoirs
A modified falling-head method was adapted to measure the saturated hydraulic conductivity (K s ) 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 K s measurement. The water level in the bucket was kept constant using overflow. A manometer (length = 100 cm and diameter = 1 cm) 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, A schematic diagram of the falling-head permeameter was used to measure the saturated hydraulic conductivity of the reservoirs 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 K s is given by Eq.
(2) (Stein 1998): where h 0 is the initial height of water level in the manometer tube above the free water surface (cm), h is the height of water level in the manometer tube at time t (cm), A is the surface area of the reservoir (cm 2 ), L is the average wall thickness of the reservoir (cm), K s is saturated hydraulic conductivity (cm/day), and t is cumulative time (day). A plot of ln (h/h 0 ) versus time, t, gives a straight line. A is the cross-sectional area of a manometer tube (m 2 ). The K s 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.

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 a 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, 10,000, and 15,000 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 (K s ) of the soil was measured on the same samples using the constant-head method. The K s was calculated using the Darcy (1856) where θ s is saturated water content (cm 3 /cm −3 ), θ r is residual water content (cm 3 /cm −3 ), S e is effective saturation; K(S e ) is unsaturated hydraulic conductivity (cm/day), and n and α (1/cm) 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.

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 from the soil layers of 0-10, 10-25, 25-40, and 40-75 cm for 2 years after plant establishment using TDR apparatus (TRIME-FM model with 0.1% accuracy) (Fig. 2).
The water content distribution in the soil was modeled using HYDRUS-2D/3D (Šimůnek et al. 1999). Assuming a homogeneous and isotropic soil, the governing equation for water flow is the 2D Richards' equation as follows: where, θ is volumetric water content (cm 3 /cm −3 ), h is pressure head (cm), t is time (day), x and z are horizontal and vertical space coordinates, respectively, and K s is saturated hydraulic conductivity (cm/day). The S (1/day) as a sink term indicates the rate of water uptake by the roots from the soil. The HYDRUS-2D/3D uses the Galerkin finite-element method to solve the transport equations, as explained in detail by Šimůnek et al. (1999). The van Genuchten parameters (θ s , θ r , n, α and K s ) were optimized by the inverse modeling in HYDRUS-2D/3D. The domain was defined by an axisymmetric 2D geometry. The soil and reservoirs were defined as two soil materials. The inlet boundary condition of reservoirs was defined as variable head. In addition, boundary conditions on the soil surface were defined as atmospheric boundary condition to consider the evaporation from the soil surface. The right and left sides of the modeled environment were set no flux, because the water flow was symmetrical around the vertical line passing through the center of the reservoirs and wetting front would not reach the left lateral 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 finite elements (FE) size of the overall mesh was 3 cm and around the reservoir was 1.6 cm.

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 (r,z) (dimensionless) in 2D geometry as follows: Where z and r are the distances (cm) from the plant origin in the directions of z (depth), and r (radius), respectively, and z m and r m are the maximum rooting lengths (cm) from the plant origin in the z and r directions, respectively. In Eq. (8), r * and z * , and p z (-) and p r (-) are all empirical parameters controlling spatial root distribution. These empirical parameters were considered to provide zero root water uptake (RWU) beyond z m to account for asymmetrical RWU in the vertical and radial directions and to allow maximum RWU within 0 to z m (Vrugt et al. 2001). 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.
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: where α(h) is a dimensionless function of soil water pressure head and S p is the potential (maximal) water uptake rate by roots (cm 3 cm −3 d −1 ).
Parameters of the Feddes model for water uptake by root are as follows: h 1 is the pressure head below which plant roots begin to extract water from the soil, h 2 is pressure head below which plant roots extract water at the maximum rate, h 3.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 , h 3.high is limiting pressure head at potential transpiration rate of 0.5 cm day −1 , h 4 is pressure head below which plant roots stop to take up water, which is usually taken at wilting point, r 2H is potential transpiration rate at which the limiting pressure head allows extracting water at the maximum rate and r 2L 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.

Evaluation of 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 (R 2 ) defined as follows (Willmott 1982): where N is the number of data points, P i is the i th predicted data point, O i is the i th observed data, and 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.

Water Depletion from the Reservoirs
The percentage of water depletion of the reservoirs with different permeabilities is presented in Table 5. We defined three treatments, namely, reservoirs with low (22, and 24%), medium (30%), and high permeability (48, and 50%) in this study according to the percentage of water depletion. At high cooking temperature, organic materials were completely burned, and different pore spaces were produced in their places. Therefore, this resulted in different permeabilities of the reservoirs concerning the quantity of the organic matter. Therefore, it was expected that the organic materials had minimal nutritional and direct influences on the survival and growth of seedlings. Plants with a high transpiration demand can be suitably irrigated with R 4 and R 5 reservoirs because of their high depletion rates. Plants with a low transpiration demand can be suitably irrigated with the conventional reservoirs made of clay, and low percentages of manure and wheat straw (i.e., R 1 and R 2 ) due to their low depletion rates. It seems that the use of high percentages of manure and wheat straw might induce the formation of large pores and microcracks in the walls of reservoirs, and therefore, result in a higher water depletion rate. The saturated hydraulic conductivity of reservoirs is also presented in Table 5. The K s varied from 0.06 (low permeability) to 0.2 (high permeability) cm/day. In other words, as the porosity and macropores of the pottery wall increases, the hydraulic conductivity increases. The pattern of soil water distribution is affected by the physical and hydraulic properties of clay pots as well as the physical properties of soil (Ghorbani Vaghei et al. 2016). Moreover, the water flow rate through the clay pots is affected by some factors including the hydrostatic pressure, the K s of the clay pot material, wall thickness, surface area, soil type, plant species, and the rate of evapotranspiration (Vasudevan et al. 2011). Figure 3 shows the laboratory-measured and predicted data of soil water characteristic curve. The soil water characteristic curve is an important soil hydraulic function which indicates the relationship between soil water content and matric suction. It affects water flow and storage, and consequently wetting pattern in the soil. This curve in the wet range (i.e., low matric suctions) is mainly affected by soil structure and macropores and in the dry range (i.e., high matric suctions) is mostly affected by soil texture and micropores. The van Genuchten-Mualem model has an excellent fit to the measured data.

Optimal Values of Soil Hydraulic Parameters from Inverse Modeling
Soil hydraulic parameters were optimized in the inverse modeling by HYDRUS-2D/3D. The HYDRUS-2D/3D software has been widely used for inverse modeling of water flow. The estimated soil hydraulic parameters by inverse solution are presented in Table 6. The optimized values of soil hydraulic parameters can be compared with those estimated by fitting the van Genuchten-Mualem model to the laboratory-measured data. The optimized θ s values (0.50 cm 3 cm −3 ) were higher than the measured value (0.38 cm 3 cm −3 ). The parameter n determines the shape of soil water characteristic curve, and its optimized values (average 1.56) were higher than the fitted value (1.09). The parameter α for the model fitting to the measured data was equal to 0.008 cm −1 which is greater than its optimized value (i.e., 0.0018 cm −1 ). Thus, the unsaturated zone occurs at lower matric suctions for the fitted set of parameters than the optimized set of parameters. The measured K s (62 cm/ day) was greater than the optimized value (average 35 cm/day). The differences between the optimized and fitted sets of parameters can be explained by the fact that the optimized parameters are obtained by using the dynamic soil water content distribution in different parts of the soil profile and taking into account both wetting and drying processes, but in the laboratory, the hydraulic properties of the soil samples were measured at equilibrium condition (e.g., in pressure plate) or at steady-state condition (e.g., constant-head method).

Water Distribution and Wetting Pattern in the Soil
The wetting patterns in the soil and the spatial distribution of soil water and matric potential depends on soil hydraulic properties, emitter discharge rates, spacing and placement of emitters, irrigation amount and frequency, water uptake rates, and root distribution patterns (Gardenas et al. 2005). As shown in Figs. 4, 5, and 6, the water slowly moved from the reservoir with low permeability, but it moved faster and in all directions from the reservoir with high permeability. In other words, with increasing the hydraulic conductivity of the reservoir, the dimensions of the wetted area around the reservoir increased, and the water content moved faster. As it can be seen, in the low permeability treatment, the water content values are at the lowest values after 24 h, the water content values slowly increased after three days, and finally, the soil around of the plant was saturated and water is slowly provided to the plant. The trend was the same in the moderate permeability treatment, but in the high permeability treatment, a large volume of water was discharged for almost three days, so the soil water content after one day was greater than the other two treatments, and more water was deep-percolated and lowered the soil water available to plant roots. One of the important goals of pressurized irrigation systems, especially drip irrigation systems, is to prevent deep percolation and water outflow from the plant root area (Khorami et al. 2013). The irrigation system should be designed so that water is evenly distributed in the root zone. In the low permeability and medium permeability treatments, it is observed that the maximum soil water content occurred in the 20-40 cm layer overtime. The soil water content after three days of irrigation in the low permeability treatment ranged from 0.25 to 0.14 cm 3 cm −3 . In the moderate permeability treatment, it changed from 0.25 to 0.13 cm 3 cm −3 , and finally, in the high permeability treatment, the water content varied from 0.26 to 0.12 cm 3 cm −3 . Such differences might be related to different discharge rates and hydraulic conductivities of the used reservoirs. Li et al. (2004)  subsurface irrigation depended on the soil type and layering, the reservoir outlet flow rate, the time of cessation of irrigation, etc. In their research, they examined the distribution of the wetting pattern in the subsurface irrigation for two loamy and sandy soils, which can show the shape of the moisture distribution in the horizontal and vertical sections with exponential functions. They stated that with increasing discharge rate, the horizontal distribution of water content would increase, and with decreasing discharge rate, the water content distribution would increase in the vertical direction. The distribution of soil water content around the reservoir must be known for the proper management of subsurface systems to wet the plant root zone uniformly. It is also aimed to increase the efficiency of the water/fertilizer use, and to maintain a dry soil surface for low water losses due to evaporation. The horizontal and vertical expansion of the water content profiles is a function of the flow rate and the operating time of the reservoir (Mohammadzade et al. 2015). The water content distribution at different distances and depths and for different treatments at one, three, and six days after irrigation in different permeability treatments are presented in Figs. 4, 5, and 6 and using Surfer 10.0 software and the Kriging method. The subsurface irrigation system is based on the use  of continuous water content, which determines the proper and uniform irrigation. In general, the highest soil water content was observed near the reservoir, and the water content decreases with increasing distance from the reservoir.

Comparison of Measured and Simulated Soil Water Content Data
The calibration results of measured and simulated values of soil water content with HYDRUS-2D/3D software in reservoirs with low, medium and high permeabilities are shown in Fig. 7. HYDRUS-2D/3D software was able to simulate the water flow and soil water content with high values of coefficient of determination (R 2 = 0.90-0.95). This shows the high accuracy of simulation of soil water content pattern in the two-dimension, which was higher in the lower discharges (Table 7). Kandelous and Šimůnek (2010) simulated the water movement in a subsurface drip irrigation system and concluded that the correspondence between simulations and observations was very good. Nayebloie et al. (2015) simulated the 2D water content distribution in the subsurface drip irrigation and stated that soil water flow can be reasonably simulated while root water uptake and evaporation processes are both actively going on.

Conclusions
In this study, reservoirs with different permeabilities were made by using manure and wheat straw, so that the saturated hydraulic conductivity (K s ) of reservoirs varied from 0.06 to 0.2 cm/day. Then, the reservoirs were used in the field for reservoir irrigation of saplings and HYDRUS-2D/3D software was used to simulate wetting pattern in the soil profile around the reservoir.
With increasing the hydraulic conductivity of the reservoirs, the dimensions of the wetted area increased, and the water content moved faster in the soil. The water slowly released from the reservoirs with low and medium permeabilities, but in the high permeability treatment, a large volume of water was released for almost three days, so that more water was deep-percolated which lowered the soil water available to plant roots leading to adverse effects on the plant growth. In the low and medium permeability treatments, it was observed that the maximum soil water content occurred in the 20-40 cm layer overtime. Plants with a high transpiration demand can be suitably irrigated with the reservoirs of their high depletion rates. Plants with a low transpiration demand can be suitably irrigated with the conventional reservoirs made of clay, and low percentages of manure and wheat straw due to their low depletion rates. The HYDRUS-2D/3D software was able to simulate soil water content during the growth period of plants with a high correspondence between the measured and predicted data. This finding showed the high accuracy of simulating the water content pattern in the two-dimension by HYDRUS-2D/3D. It is suggested that the results of HYDRUS-2D/3D simulations and different irrigation scenarios be used to determine the water content distribution and suitable places for planting seedlings.