1 Introduction

Furrow irrigation is currently the most common surface irrigation practice in the arid region of northwest China. Water flow under furrows laterally infiltrates through the soil to meet the needs of the plants grown on the ridges or raised beds. In traditional furrow irrigation, more water and solute fluxes are transferred in the furrows than planting beds (Bargar et al., 1999; Chen et al., 2011; Zhang et al., 2015), and deep water percolation below the root zone more easily occurred in the furrows than the ridges (Bargar et al., 1999; Abbasi et al., 2003). Additionally, salt concentration within the surface soil below ridges was lower than below the furrows under furrow irrigation with saline water (Chen et al., 2015). Improper furrow irrigation design in the field has higher water deep percolation and lower irrigation uniformity (Moravejalahkami et al., 2009). Hence increased lateral spreading and decreased moisture penetration depth are desirable. Moreover, irrigation of mulched furrows may produce pronounced changes in soil water and solute distribution between furrows and ridges (Bezborodov et al., 2010).

Soil water dynamics under irrigation of mulched furrows are complex, but critical for management of irrigation design, soil storage, and water use efficiency (Tabuada et al., 1995; Abbasi et al., 2003). Both soil physical properties and irrigation design parameters affect water and solute fluxes within the soil media (Zhang et al., 2005; Nie et al., 2009; Zhang et al., 2012). Therefore, these influential variables must be taken into consideration for the optimized design of mulched furrow irrigation systems. Furthermore, effective management of irrigation systems requires accurate estimation of wetted soil volume. The wetted soil depth should be consistent with the anticipated root depth, while wetted soil width should correspond with the spacing of ridges (Zur, 1996). Proper irrigation design can minimize water leaching and nutrient loss (Skaggs et al., 2010). The traditional method of describing furrow infiltration and wetting patterns was investigated using simulation experiments (Nie et al., 2009; Chen et al., 2011). Soil water distribution at selected nodes was measured to draw iso-water content lines (Bargar et al., 1999; Abbasi et al., 2003; Li et al., 2003). Additionally, cumulative infiltration, wetting patterns, and soil water content (SWC) were monitored to investigate soil water dynamics (Strelkoff and Sousa, 1984; Izadi and Wallender, 1985; Trout, 1992; Abbasi et al., 2003). However, knowledge of the effects of soil physical properties and irrigation parameters on the wetted volume was limited due to the restriction of wetting pattern observations in the field. Therefore, specialized and efficient methods which could elaborate soil water dynamics under mulched furrow irrigation systems are necessary.

Mathematical simulation models provide a potential approach for investigating lateral and vertical soil water infiltration (Tabuada et al., 1995; Ebrahimian et al., 2013; Zhang et al., 2013). Many mathematical models have been established based on Richards' equation in a variably saturated, rigid, isotropic porous soil medium (Richards, 1931), which can simulate more flexible representation of the flow domain, boundary conditions, and soil properties (Warrick, 2003). Levenberg-Marquardt optimization algorithm is used to inversely estimate soil hydraulics and/or solute transport parameters. Numerous models have been developed to simulate water flow and solute transport (e.g., bromine and nitrogen) and assess the effects of irrigation schedule and management practices on nitrogen leaching under furrow irrigation (Abbasi et al., 2004; Mailhol et al., 2007; Crevoisier et al., 2008; Ebrahimian et al., 2013). Simulation methods are also less costly and time consuming than conventional field experiments. However, the effect of irrigation of plastic mulched furrows on soil water distribution remains elusive. Mulched treatment retained 40%~60% more water, and saturated hydraulic conductivity of mulched treatment was 123 times greater in unmulched treatment (Humberto and Lal, 2007). The modified soil properties caused by plastic mulching must be considered when estimating the wetted soil volume in the root zone (Chen et al., 2015). Therefore, it is essential to construct a model that represents two-dimensional soil water dynamics under mulched furrow irrigation systems.

In this study, irrigation of mulched furrows in a cross-sectional ridge-furrow configuration was investigated in a soil chamber, which avoided the added complexity caused by vegetation and site-specific factors. The objectives of this study were to: (1) analyze the physical infiltration process and effects of experimental variables on wetting patterns; (2) develop wetting pattern propagation equations to predict the wetted lateral and downward distances; (3) establish a mathematical model to simulate soil water flow under mulched furrow irrigation systems.

2 Materials and methods 2.1 Soil properties

Two-dimensional infiltration characteristics of two soil types under mulched furrow irrigation were tested in a cuboid soil chamber. Silt loam soil and sandy loam soil (US Department of Agriculture) were collected at 10~30 cm depth from fallow fields. Soil textural fractions were analyzed using the Laser Mastersizer 2000 (Malvern Instruments, Malvern, England). Bulk density was determined by ring soil samples. Saturated soil water content (SWC) of the soil samples was measured using the oven-dried method. Saturated hydraulic conductivity was estimated using the ROSETTA code (Schaap et al., 2001). Some physical properties of the field-collected soil samples are summarized in Table 1.

2.2 Experimental setup

The experimental setup consisted of a cuboid soil chamber and water supply system (a Mariotte flask) (Figure 1). The cuboid soil chamber, made from 10-mm-thick plexiglass material, was 70 cm long, 5 cm wide, and 70 cm high. The bottom of the soil chamber included numerous 2-mm parallel air vents for ventilation. The soil wetting pattern was inscribed on two selected vertical surfaces of the soil chamber to validate the experimental accuracy. The Mariotte flask with a flexible hose was connected to a rubber hole on the soil chamber to supply water. This system maintained a constant water level under mulched furrows. Previous studies have found that the experimental chamber shape has negligible impact on the water flow pattern (Zhang et al., 2005; Nie et al., 2009). The middle of the furrow cross-section and ridge cross-section are zero flux boundaries in homogeneous soil conditions. Therefore, only half of the soil water distribution with mulched furrow irrigation (i.e., from the middle of the furrow to the middle of the ridge) was simulated in the soil chamber.

2.3 Experimental process

The soil was air-dried at ambient temperature immediately after collection from the field. The air-dried soil was passed through a 2-mm sieve and prepared for packing into the soil chamber (Figure 1). Before packing, the soil was uniformly watered using a watering pot until the intended initial SWC (Table 2) was reached. The prepared soil was loaded and compacted into the soil chamber in 5-cm increments to obtain the required bulk density (Table 2) and homogeneous soil profiles. The packed soil in the entire chamber was allowed to equilibrate for 24 hours to ensure uniform distribution of the intended initial SWC throughout the chamber (Table 2). In order to avoid soil evaporation, the soil was covered with a polyethylene sheet before the experiment.

To investigate the effect of initial SWC, irrigation water level, and soil texture on two-dimensional infiltration under mulched furrow irrigation, six experimental treatments with two replications were conducted in the soil chamber (Table 2). The manipulated variables in the experimental treatments were determined using actual field values. The packed soil was excavated into an isosceles-trapezoid-shaped furrow and shaped into a ridge-furrow configuration (Figure 2). The furrow was mulched with plexiglass material and sealed using silicone rubber. The furrow was 15 cm below the ridge, with slope coefficient of 1. The furrow water level, the standing water height measured from the furrow bottom, was measured using a steel ruler. The water level was maintained by adjusting the Mariotte flask at the beginning of experimental treatments. Each experiment lasted for 120 min to 240 min to ensure free drainage at the bottom boundary. During each experiment, positions of the moving water front in horizontal and vertical directions were recorded with a marker on the two selected vertical surfaces of the soil chamber every 20 min. The cumulative infiltration was measured with readings of the graduated Mariotte flask water volume every 5 min. At the end of the experiment, the soil chamber was placed horizontally and opened from one vertical surface of the soil chamber for immediate soil sampling. The sampling layout was at horizontal intervals of 5 cm, starting 2.5 cm from the center of the furrow and moving outward to the edge of the wetted area. The soil was also vertically sampled in 5 cm intervals, starting 2.5 cm from the soil surface and moving downward to the wetted front (Figure 2). Soil samples were collected from outside the wetted soil to confirm the initial SWC. The sampled soil was oven dried at 105 ℃ for 8 hours to determine the gravimetric soil water content. The two-dimensional wetting patterns were characterized by the lateral and downward distances (Figure 2). The lateral distance (R_l) was measured from the foot of the ridge shoulder to the edge of the wetting front, and the downward distance (R_d) was measured from the soil surface to the wetting front.

2.4 Mathematical simulation 2.4.1 Numerical model We simulated two-dimensional soil water flow using HYDRUS-2D (Šimůnek et al., 2006). The measurements of water infiltration under mulched furrow irrigation in the soil chamber were used to verify the accuracy of HYDRUS-2D simulated predictions of soil water distribution. Water infiltration and redistribution in the cross-sectional ridge-furrow configuration was assumed to be a two-dimensional isothermal Darcian water flow in a variably saturated rigid porous medium (Abbasi et al., 2004; Mailhol et al., 2007; Crevoisier et al., 2008). Subsurface water flow was described using Richards' equation (1931) in its two-dimensional form:

$\frac{{\partial \theta }}{{\partial t}} = \frac{\partial }{{\partial x}}\left({K(h)\frac{{\partial h}}{{\partial x}}} \right) + \frac{\partial }{{\partial z}}\left({K(h)\frac{{\partial h}}{{\partial z}}} \right) + \frac{{\partial K(h)}}{{\partial z}}$

(1)

where θ is volumetric water content (cm³/cm³); h is soil water pressure head (cm); t is time (min); K(h) is unsaturated hydraulic conductivity (cm/min); x is the horizontal space coordinate (cm); z is the vertical space coordinate. The unsaturated hydraulic conductivity function, K(h), was described using the capillary model of Mualem (1976). The soil water retention curve, θ(h), was described using the closed-form equation of van Genuchten (1980) as follows:

$K(h) = {K_s}{S_e}^l{[1 - {(1 - {S_e}^{\frac{1}{m}})^m}]^2}$

(2)

$\theta \left(h \right) = \left\{ \begin{array}{l} {\theta _r} + \frac{{{\theta _s} - {\theta _r}}}{{{{\left({1 + {{\left| {\alpha h} \right|}^n}} \right)}^m}}}\;\;\;h < 0\\ {\theta _s}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;h \ge 0 \end{array} \right.$

(3)

where the effective saturation (S_e) is calculated according to the formula:

${S_e} = \frac{{\theta - {\theta _r}}}{{{\theta _s} - {\theta _r}}}, m = 1 - \frac{1}{n}, n > 1$

(4)

θ_r is residual soil water content (cm³/cm³); θ_s is saturated soil water content (cm³/cm³); K_s is saturated hydraulic conductivity (cm/min); m and n are empirical parameters related to the pore-size distribution; α is an empirical constant that is inversely related to the air-entry pressure value (min^-1); and l is the empirical shape parameter. 2.4.2 Initial and boundary conditions We simulated a presumed two-dimensional (2D) symmetric vertical soil profile (Mailhol et al., 2007; Ebrahimian et al., 2013). Initial and boundary conditions for the simulation of cross-sectional ridge-furrow configuration are presented in Figure 3. The simulated flow domain, for which the numerical solution was obtained, was 70 cm in width (Figure 3, BC section) and 70 cm in depth (Figure 3, DC section). The initial condition in the simulation domain was intended initial SWC throughout the soil chamber. During water application, the upper boundary condition along GF was kept at a constant pressure head (h₀). The mulched furrow AG was zero flux. Atmospheric condition was applied to the soil surface ED and EF. Mean daily soil evaporation used for the atmospheric condition was 5.7 mm/d (Zhang et al., 2010). Zero flux boundary conditions were assumed in middle furrow cross-section AB and ridge cross-section CD. Water flux at the bottom boundary (BC section) was simulated as free drainage. 2.4.3 Parameter estimations and simulation processes The simulations of water flux were subjected to various uncertainties resulting from parameterization of soil hydraulic parameters (Loos et al., 2007). To analyze the effect of uncertainty in soil hydraulic parameters on the model output, two parameterization approaches were compared. In simulation method I (NNP), soil hydraulic parameters were derived using the pedotransfer function of the ROSETTA software package (Schaap et al., 2001) via the neural network based on data of soil particle size distribution and bulk density (Table 1). In simulation method II (Inverse), an inverse approach was used to estimate soil hydraulic parameters using a Levenberg-Marquardt optimization algorithm implemented in HYDRUS-2D. The hydraulic functions described in Equations (2) and (3) contain six unknown parameters (θ_r, θ_s, α, n, K_s, and l). Experiments 1 and 4 were used to optimize the soil hydraulic parameters of the van Genuchten-Mualem model in silt loam soil and sandy loam soil. Cumulative infiltration was used in the objective function during the Inverse approach. Parameters n and K_s were simultaneously optimized with other constant parameters (θ_r, θ_s, and α, ). The parameter l was set as 0.5. Experiments 2, 3, 5, and 6 were conducted to validate the accuracy of estimated soil hydraulic parameters in silt loam soil and sandy loam soil. The geometry of experimental ridge-furrow configuration and water level were required for the HYDRUS-2D mathematical simulation. In order to identify the effects of unmulched and mulched furrows on cumulative infiltration and wetting patterns, simulations 1 and 2 used in HYDRUS-2D model are listed in Table 3. 2.5 Statistical analysis

Performance of the HYDRUS-2D model was assessed using the root mean square error (RMSE), which provides a quantitative comparison of the goodness-of-fit between simulated and observed values (Patel and Rajput, 2008; Kandelous and Šimůnek, 2010). The RMSE was calculated using the formula:

$RMSE = {\left[ {\frac{1}{n}\sum\limits_{i = 1}^n {{{(I{{(o)}_i} - I{{(s)}_i})}^2}} } \right]^{\frac{1}{2}}}$

(5)

where n is the total number of data in each experiment; I(o)_i is the ith observed data; I(s)_i is the ith simulated data.

The parameters of empirical wetting patterns propagation equations were estimated by using OriginLab 8.0 software (OriginLab Corporation, Northampton, MA, USA). Statistical analysis was performed using the software program SPSS, ver. 17.0 (SPSS Inc., Chicago, IL, USA). We prepared contour maps of the spatial distribution of soil water content increment (SWCI) in the SURFER package (Golden Software Version 8.0, 2002), using Kriging with a linear variogram.

3 Results and discussion 3.1 Soil water distribution

Soil water distribution under mulched furrow irrigation determines the optimum spacing and number of planting rows. The wetted soil volume achieved through irrigation should match the distribution of the root system. Soil water content increment (SWCI) in a semi ridge-furrow configuration for all the experimental treatments is illustrated in Figure 4. Mulched furrow irrigation clearly increased soil moisture lateral spreading, and resulted in excessively deep water infiltration under ridge shoulders and ridges. Capillary force initially drove lateral subsurface water flow, which wetted the ridges. Downward flow trends over time occurred due to gravitational forces. The volume of wetted soil tended to increase with an increase in initial SWC in both silt loam soil (Figures 4a and 4b) and sandy loam soil (Figures 4d and 4e). Previous studies have reported that soil hydraulic conductivity increased and unsaturated wetted zone can be easily saturated when initial SWC was greater (Enciso-Medina et al., 1998; Nie et al., 2009; Li and Liu, 2011). Similar results have also been reported in field measurements and model simulations under surface irrigation (Zhang et al., 2005; Zhang et al., 2013). Therefore, irrigation amount should take into account the initial SWC because soil moisture affects crop germination and growth. Higher furrow water level would trigger greater wetted soil volume in both silt loam soil (Figures 4a and 4c) and sandy loam soil (Figures 4d and 4f). Higher furrow water level produced a larger surface infiltration area and gravitational potential on furrow, which in turn increased the wetted soil volume. Previous studies show a positive correlation between cumulative infiltration and wetted perimeter without overlapped wetting patterns in stagnant blocked furrows (Strelkoff and Sousa, 1984; Izadi and Wallender, 1985; Abbasi et al., 2003). Irrigation with a higher furrow water level increased the uniformity of water distribution throughout the root zone. The volume of wetted soil was greater in sandy loam soil than silt loam soil, due to the greater saturated hydraulic conductivity of sandy loam soil (Table 1). Soil texture was the primary factor that determined the soil matric potential gradient to drive subsurface water movement.

3.2 Movement of the wetting distances

Ponding infiltration under cross-sectional irrigation of mulched furrows is regarded as a two-dimensional lateral and downward water fluxes (Figure 2). The lateral and downward wetted distances after each infiltration event are listed in Table 4. The wetting distances significantly differed between the two soil types. The lateral distance in silt loam soil was greater than the downward distance. The R_l/R_d ratio was more than 1 in Experiments 1-3. However, in sandy loam soil, the lateral distance never exceeded the downward distance. The average R_l/R_d ratio was 0.88 in Experiments 4-6.

Soil water movement under mulched furrow irrigation was mainly driven by capillary and gravitational forces. The wetting patterns moved outward by the relative effect of gravitational over or below capillary driving forces. Lateral distance, downward distance, and R_l/R_d as a function of infiltration time in silt loam soil and sandy loam soil are presented in Figure 5. The wetting front initially moved laterally and downward rapidly, and then slowed down with the increase in infiltration time. Figure 5 illustrates the exponential relationship between the wetted lateral and downward distances and t^0.5 for all experimental treatments. The variation of lateral distance (R_l, cm) with infiltration time (t, min) is described by the equation:

${R_l} = {L_1}{t^{0.5}} + {L_2}$

(6)

A correlative analysis between downward distance (R_d, cm) and infiltration time is described by:

${R_d} = {L_3}{t^{0.5}}$

(7)

Similarly, the variation of R_l/R_d with infiltration time (t, min) is described by:

${R_l}/{R_d} = {L_4}{t^ - }^{0.5} + L{}_5$

(8)

where L₁ to L₅ are the fitting coefficients.

The R² values of wetted lateral distance, downward distance, and R_l/R_d for all experimental treatments were more than 0.87. The regression parameters in Equations (6)-(8) were all significant at the P < 0.01 level (Figure 5). Thus, Equations (6) and (7) can be used to estimate the wetted lateral and downward distances for mulched furrow irrigation. The greater L₁ and L₃ in Experiments 4-6 illustrated that the wetting front moved faster in sandy loam soil than in silt loam soil. The asymptotic values of L₅ in silt loam soil and sandy loam soil were 0.93 and 0.62, respectively (Figure 5). In silt loam soil, the downward wetting distance was relatively equal to the lateral distance. This result is comparable to findings of Nie et al. (2009), who reported a linear relation between the wetting distances and t_0.5in furrow irrigation. As indicated in Equations (6) and (7), the fitting coefficients must be calibrated before these relationships can effectively be used to design the mulched furrow irrigation system.

3.3 Model validation

Before modeling two-dimensional infiltration under mulched furrow irrigation, model calibrations were performed using experiment 1 and 4 for silt loam soil and sandy loam soil, respectively. In Figure 6, the optimized model parameters based on the inverse optimization (Inverse) were clearly different from the neural network parameter estimates (NNP). The soil water retention curve, θ(h), significantly differed between the two parameterization approaches. With the same pressure head, the Inverse simulation estimated lower soil water content than NNP. The estimated soil hydraulic parameters were compared with the parameter ranges documented in previous publications (Šimůnek et al., 2006; Nie et al., 2009).

The soil hydraulic parameters estimated through the HYDRUS-2D model were validated using data from Experiments 2 and 3 for silt loam soil and Experiments 5 and 6 for sandy loam soil. The accuracy of model parameters was further assessed using RMSE. Cumulative infiltration and downward distance RMSE values showed a better goodness-of-fit for the Inverse methods relative to NNP (Table 5). The observed and predicted cumulative infiltration over time in the two soil types is presented in Figure 7. The HYDRUS-2D model with the inverse parameters effectively simulated cumulative infiltration and ridge, ridge shoulder, and furrow soil water content for both soil types (Figure 8). There were significant discrepancies between measured soil water content and predictions using the NNP method, with lower predicted soil water contents than the observed values. Therefore, the HYDRUS-2D simulations with the neural network-estimated parameters did not provide a good description of soil water movement. Saturated soil hydraulic conductivity was likely influenced by effects of soil repacking on the soil structure in the chambers used for these experiments. The hydraulic conductivity predicted by NNP was lower than the Inverse method predictions because of the effects of entrapped or dissolved air in the soil (Abbasi et al., 2004). Discrepancies between observed and predicted soil water contents are consistent with similar discrepancies reported in other studies (Abbasi et al., 2004; Ebrahimian et al., 2013), because no hysteresis of soil water characteristic curve was set in the simulations. Therefore, the inverse optimized parameters were more accurate than the estimated parameters based on neural network. These results confirm that the HYDRUS-2D model can be used to accurately simulate two-dimensional soil water dynamics under mulched furrow irrigation.

3.4 Effects of mulched furrows relative to unmulched furrows on soil water distribution

Field experiments and model simulations under furrow irrigation found that deep water and solute percolation below the root zone more easily occurred in furrows than ridges (Bargar et al., 1999; Abbasi et al., 2004; Chen et al., 2011). In order to increase the wetted lateral spreading, the effects of mulched and unmulched furrows on soil water distribution were simulated using the HYDRUS-2D model with inverse parameters. The wetting pattern under unmulched furrow irrigation differed from under mulched furrow irrigation in silt loam soil (Figure 9). The water infiltration contact area only reached the ridge shoulder position under mulched furrows, whereas the contact area with the furrow and ridge shoulder was greater under unmulched furrows. The greater wetted perimeter and gravitational potential under unmulched furrows, thus the more cumulative infiltration was produced. Lateral distances predicted by both simulation cases were relatively equal (Table 6). However, downward distance under unmulched treatment increased by 32.14% compared to under mulched treatment. The R_l/R_d ratio for mulched furrows was greater than for the unmulched furrows. These results confirm that deep water and solute percolation more easily occurred under unmulched furrow irrigation. The center of the wetted soil volume moved on ridge shoulder under mulched furrow treatment (Figure 9). The R_l was greater than R_d under mulched furrows (Table 6), illustrates that the effect of capillary forces driving lateral infiltration exceeded the gravitational forces driving downward infiltration. The speed of irrigated water flow under mulched furrows is faster, which promote irrigation uniformity between ridges and furrows and reduce soil evaporation (Zhang et al., 2005). Hence irrigation of film-mulched furrows is an effective practice.

4 Conclusions

Mulched furrow irrigation clearly increased soil moisture lateral spreading on ridge shoulders and ridges, due to enhancement of capillary driving force. Increasing initial SWC and irrigation water level resulted in an increase of the wetted soil volume. Simple regression equations describing the relationships between wetting front and elapsed time can be used to estimate the wetted lateral and downward distances, which conveniently provide guidelines and criteria for irrigation design. The optimized parameters by the Inverse approach were accurate, and the observed and simulated values were very close in the mathematical model. Based on model calibrations, HYDRUS-2D model provided a basis to evaluate the effect of mulched furrows on soil water distribution. There were significant differences in the wetting pattern and soil water distribution between unmulched and mulched furrows. The downward wetted distance under mulched furrows decreased by 32.14% compared to under unmulched furrows. Film-mulched furrows is recommended for furrow irrigation. The information provided by the HYDRUS-2D model is helpful for field mulched furrow irrigation design and management.

Acknowledgments:

This work is supported by National Natural Science Foundation of China (NO. 41401036), China Postdoctoral Science Foundation (NO. 2015T81070, 2014M560818), and West Light Foundation of the Chinese Academy of Sciences. We would like to thank Christine Verhille at the University of British Columbia for her assistance with English language and grammatical editing of the manuscript. We also would like to thank reviewers and the editors for their constructive comments on this manuscript.