Numerical investigation of wetting front migration and soil water distribution under vertical line source irrigation with different influencing factors Open Access

α

an empirical parameter inversely related to the air entry value (cm⁻¹)

r

radial coordinate

z

vertical coordinate

θ

soil water content (cm³/cm³)

θ_s

saturated water content

θ_r

residual water content (cm³/cm³)

θ_i

initial water content

θ_f

field capacity (cm³/cm³)

θ_w

wilting coefficient (cm³/cm³)

K(h)

unsaturated hydraulic conductivity (cm/min)

K_s

saturated hydraulic conductivity (cm/min)

Se

the effective degree of saturation

M_i

the ith measured value

S_i

the ith simulated value

N

the total number

n and m

empirical parameters affecting the shape of the soil water retention curve (m = 1 − 1/n)

X

the horizontal wetting width from point C (cm)

Y

the vertical upward wetting height from point C (cm)

Z

the vertical downward wetting depth from point C (cm)

ST

soil texture

B

line source buried depth

D

line source diameter

L

line source length

WFM

wetting front migration

VLSI

vertical line source irrigation

MAE

mean absolute error

RMSE

root mean square error

PBIAS

percent deviation

NSE

Nash efficiency coefficient

Northwest China has conditions (rich land resources, sufficient sunshine, and a large temperature difference between day and night) that are beneficial for the development of fruit production (Mamitimin et al. 2015; Du et al. 2017). However, rainfall in the area is limited, surface evaporation is strong, and fruit production depends to a large extent on irrigation (Shen et al. 2013; Guo & Shen 2016). Surface irrigation is a common irrigation method in orchards that is simple to operate but easily wastes water (Kang et al. 2003; Deng et al. 2006). In recent years, micro-irrigation techniques such as micro-sprinkler irrigation, surface drip irrigation and subsurface drip irrigation have achieved good savings in water when applied to agricultural production (Pereira et al. 2002; Namara et al. 2007; Seo et al. 2008; Cremades et al. 2015). However, when these techniques are applied to orchard irrigation in Northwest China, they have the following limitations: (1) the climate is dry and it is frequently windy in Northwest China. Water evaporation loss is large, and water use efficiency is low during micro-sprinkler irrigation (Yazar 1984; Stambouli et al. 2012). (2) The quality of water used for irrigation is poor and filtration systems are imperfect in Northwest China. Subsurface drip irrigators are easy to block, and the monitoring and repair of clogged irrigators is difficult (Zhou et al. 2017; Zhang et al. 2020). (3) The root distribution of fruit trees is deep, and the water supply path is long. When using surface drip irrigation, the soil evaporation is large and nutrients easily accumulate on the soil surface (Meshkat et al. 1998, 2000; Yanni et al. 2004). In view of this, it is necessary to develop appropriate irrigation methods to promote the development of the fruit industry and improve the efficient use of water resources.

Vertical line source irrigation (VLSI) is an irrigation method suitable for deep-rooted plants. Its characteristics include a large flow area, high irrigation efficiency, strong adaptability and integration of water and fertilizer (Wang et al. 2012). Its emitter is a plastic, perforated tube with a bottom seal. In this method, the emitter is inserted vertically into the soil and water is supplied directly to the plant roots through the perforations. Zeng et al. (2010) and Cheng et al. (2010) investigated the effects of θ_i and L on the water front migration of VLSI by laboratory experiments. Experimental research is a basic method for scientifically exploring soil water movement under different irrigation methods, but it is time-consuming and laborious. Numerical simulation can be used to calculate soil water movement for different soil properties or design parameters. It is an inexpensive method to determine suitable irrigation technical parameters and optimize the design of an irrigation system.

HYDRUS software was developed by vadose zone hydrology scientists van Genuchten and Šimůnek. It has been widely used in agricultural water-saving irrigation research in terms of point source or line source infiltration. The results of its calculations can better reflect the basic law of soil water movement (Šimůnek et al. 2008, 2016). Kandelous & Šimůnek (2010) evaluated the reliability of the HYDRUS-2D simulation of soil water distribution (SWD) in subsurface drip irrigation through laboratory and field experiments. On the basis of experimental verification, Li & Wang (2011) used HYDRUS-2D to simulate the SWD of VLSI under different ST, θ_i, B, D and L at a given irrigation time of 240 min. Li et al. (2013), based on the HYDRUS-3D model, simulated the effects of emitter discharge rate (Q) and L on the SWD of bubbled-root irrigation under 16 L of irrigation amount (I_a). Naglič et al. (2014) analysed the influence of Q and θ_i on the wetting front migration (WFM) for various STs with HYDRUS-2D/3D simulation results. Fan et al. (2018b) adopted HYDRUS-2D software to simulate the effects of ST, θ_i, B, D, and L on the cumulative infiltration of HYDRUS-2D simulation results. Based on the Philip infiltration formula (Philip 1957), a calculation model for the cumulative infiltration of VLSI including the line source seepage area was established. At present, the soil wetting characteristics (including WFM and SWD) have not been studied by either experimental methods or numerical simulations under the different influencing factors of VLSI to find their universal law.

The shape and size of the soil wetting body determine the plant growth status and efficiency of field water use. The size of the soil wetting body can be characterized by the migration distance of the wetting front. Subbaiah (2013) performed a comprehensive review of the prediction model of WFM during surface drip irrigation. Al-Ogaidi et al. (2016) proposed an enhanced empirical model to estimate WFM in surface drip irrigation that performed well in replicating experimental data. Fan et al. (2019) developed an empirical model for predicting the wetted radius and depth of wetting patterns during film hole irrigation using simulation data. The empirical model had a good prediction effect. Douh et al. (2018) presented a semi-empirical model for the simulation of WFM under subsurface drip irrigation and verified its applicability and dependability by comparison with experimental values. The contour line of the wetting front under VLSI has a shape similar to an ‘ellipsoid’. The developed empirical models of SWF for surface drip irrigation, subsurface drip irrigation, and film hole irrigation are not suitable for VLSI. Therefore, it is necessary to develop an empirical model that can predict the SWF during VLSI to provide a reference for the design and management of VLSI systems.

In this study, based on HYDRUS-2D software, a mathematical model of soil water movement in VLSI was established and solved. The reliability of the calculation results was verified by comparison with experimental values. Subsequently, the soil water infiltration process of the VLSI under different combinations of influencing factors (ST, θ_i, B, D, and L) was simulated, and the WFM distance at different times and the SWD after irrigation cut-off were obtained. The characteristics of the soil wetting body in the VLSI were discussed. The simulation data were then used to screen the dominant factors affecting the WFM and an empirical model for predicting the wetting pattern dimensions was developed. Finally, the applicability of the empirical model was evaluated by soil box experiments.

The experimental equipment consisted of five parts: a soil container, a Mariotte bottle, an irrigation emitter, an adjustable-height stand, and a hydraulic hose, as shown in Figure 1.

The soil container was made from 10-mm-thick transparent acrylic material to facilitate the visualization of the WFM at different times. The container bottom was perforated (2 mm diameter, 20 mm spacing) to prevent air resistance during irrigation. Some round holes (1 cm in diameter and 5 cm in spacing) were reserved on one side of the container to measure the soil water content after the irrigation cut-off. Considering the axial symmetry of soil water movement in VLSI, we chose a 1/4 cylindrical plastic tube as the irrigation emitter. The bottom of the emitter was sealed, and the L-height cylindrical surface above the bottom had irrigation holes. To observe the shape and migration of the wetting front in the three directions, the gauze-wrapped emitter was placed at a corner of the sidewall while taking the soil hole, ensuring that the emitter was in close contact with both the soil profile and the container walls. An infiltration test was carried out the next day after the completion of soil filling. The water in the irrigation emitter was provided by a Mariotte bottle placed on an adjustable-height stand, and the irrigation emitter and the Mariotte bottle were connected by a hydraulic hose. The size of the Mariotte bottle was 20 cm in diameter and 100 cm in height. According to the literature (Fan et al. 2018a), the soil wetting pattern of vertical moistube irrigation is approximately ‘pear’ shape. The maximum sizes of the soil wetting pattern in the vertical upward, horizontal and vertical downward directions are located at the highest point, the middle point and the lowest point of the moistube, respectively. Therefore, to facilitate the mapping of the wet body outline, the characteristic values of the wet body in the centre (C) of the three lines of vertical downward, horizontal downward and vertical downward are selected. At different times in the experiment, the water level of the Mariotte bottle was recorded, the migration curve of the wetting front was drawn on a plexiglass plate, and then the migration distance of the wetting front was measured with a ruler. After infiltration reached set Ia, the water supply was stopped, the soil was quickly removed from the soil hole, and the soil water content was measured by the drying method (105 °C, 24 h). To minimize the experimental error as much as possible, each experiment was repeated three times, and the results were averaged.

Four soil texture types were used in the experimental study: loam and silt loam in the Lanzhou area of the Loess Plateau and sand and sandy loam in the Wuwei area of the Hexi Corridor. The soil was excavated in the field at a depth of 0–40 cm and was air-dried, rolled and mixed well, and finally passed through a 2 mm circular sieve to produce the experimental soil samples. The particle size distribution, field capacity (θ_f), bulk density (γ), saturated water content (θ_s), and saturated hydraulic conductivity (K_s) of the soil samples were measured by laser particle-size analyser (MS2000), ring knife method, drying method, and constant head permeameter, respectively. The basic physical properties of the soil samples are shown in Table 1. Before the experiment, water was added to the soil sample according to θ_i, mixed thoroughly, sealed with plastic film and allowed to stand for one day to evenly distribute the soil moisture.

Referring to the case of VLSI in the field, 4–8 irrigation emitters could be buried around a fruit tree, A single irrigation emitter can irrigate approximately 20–40 L of water each time. Zhang et al. (2010) and Wang et al. (2012) comprehensively considered the water use efficiency, yield, net benefit and other factors of jujube trees, and pointed out that under the condition of the bubbled-root irrigation of jujube trees in northern Shaanxi, a combination of two emitters per plant and 40 L emitters per plant was the most appropriate arrangement. To fully show the migration distance of the wetting front and water distribution of the wetting pattern in vertical source irrigation, the irrigation water quota was uniformly set at 40 L. The size of the water seepage interface was determined by D and L, the structural parameters of the irrigation emitter. From a practical point of view, D and L were usually set to 2–6 cm and 10–30 cm, respectively (Cheng et al. 2010; Zeng et al. 2010). B is an important design parameter for subsurface irrigation systems, and the value criterion is mainly determined by the location of the plant root distribution. Combined with the variation in L, B is generally 30–50 cm. A suitable θ_i for starting irrigation in the field was approximately 50%θ_f–70%θ_f (Du et al. 2017; Wang et al. 2017). Considering the variation in the above irrigation parameters (ST, θ_i, B, D, and L), we designed laboratory experiments (Table 2). Furthermore, sand and silt loam were used to verify the HYDRUS simulation results, while loam and sandy loam were used to evaluate the WFM model.

The soil–water characteristic curve of the experimental soils was measured by a high-speed centrifuge (SCR-20) and was fitted using the RETC (van Genuchten et al. 1992). K_s was determined by a constant head permeameter (Provenzano 2007). The van Genuchten–Mualem model parameters thus obtained are shown in Table 3.

Figure 2 shows the initial and boundary conditions used in this study to simulate different modelling scenarios.

In all simulation scenarios, the soil moisture content is set according to the initial water content; the upper boundary DE was affected by atmospheric conditions, considering that the surface of the irrigation process was a dry soil layer, and the evaporation amount was small. For the simplified calculation, the zero flux surface was set. The lower boundary FG was not affected by irrigation; it was free drainage, set according to the zero flux surface. The left boundary GH was the central axis of the emitter. AD was the plastic pipe wall, the water was exchanged, and the zero flux surface was set at the right boundary EF, where irrigation did not arrive; it was set according to the zero flux surface. The bottom of the line source was BH sealed, and the boundary was the zero flux boundary; the boundary of the water seepage surface was the full water supply mode, and the water supply was saturated quickly after the start of the water supply and can be treated according to the fixed water head boundary (Li & Wang 2011; Li et al. 2013; Fan et al. 2018b).

Numerical simulation was performed using HYDRUS-2D (Šimůnek et al. 2006). In the process of solving, the time difference was used in the implicit difference format, and the soil profile was spatially dispersed by the Galerkin finite element method. Considering the actual field and calculation accuracy requirements, the finite element calculation domain depth was 100 cm, the width was 50 cm, the space step was 1 cm, and the time step was 0.1 min. The simulation duration was determined by the irrigation quota (40 L). The ST was represented by the van Genuchten–Mualem model parameters, θ_i was set by the initial conditions, while L, D and B were achieved by the boundary conditions.

With the reliability of the numerical model and its solution method verified, 11 simulation schemes were selected to analyse the influence of ST, θ_i, D, L and B on the soil wetting patterns of VLSI as shown in Table 4.

Figure 3 shows the change process of the measured and simulated WFM distance with time in the three directions of vertical upward, horizontal and vertical downward from point C.

Figure 3 shows that the measured values were basically consistent with the calculated values of the model. The measured and simulated values of the wetting front varied nonlinearly with time. The speed of the wetting front showed a trend from fast to slow and gradually tended toward a stable value.

Figure 4 compares the simulated and measured SWD data at the end of irrigation. The curves are simulated values, the dots are measured values, and the soil water content is the volume of water content.

From Figure 4, it can be seen that both the simulated values and the measured values gradually decreased from the line source to the periphery, and the values near the line source basically reached a saturated state.

The statistical characteristics of the simulated and measured values were analysed and are listed in Table 6.

MAE and RMSE are close to 0; PBIAS < ±10; NSE is very close to 1 (NSE ≥ 0.928). The simulation results were in good agreement with the measured data, which proved that the model and the solution method were reliable. The HYDRUS-2D simulation could be used to study the influence of parameters such as the ST, θ_i, L, D and B on the effect of the soil wetting body of the VLSI.

According to the results of the 11 sets of simulation schemes in Table 4, the relationship between the moving distance and time of the wetting front in the three directions was plotted, as shown in Figure 5.

In the early stage of infiltration, the matrix potential gradient between the emitter and the wetting front was large, and the wetting front in all three directions moved faster but slowed down as the wetting front moved away from the emitter and the matrix potential gradient decreased. Because of the gravitational potential, the increase in the Z direction was greater than that in the Y direction, and the increase in the X direction was intermediate between them. Both the wetted Y and the wetted Z included the initial distance L/2 (half of the emitter length), so the values of the wetted Y and Z in the initial stage were greater than that of the wetted X, but with the passage of time, the value of the wetted X would gradually exceed that of the wetted Y. After the wetting front reached the soil surface, the value of the wetted Y remained constant, B − L/2.

Comparing and analysing Figure 5(a)–5(c), the soil texture had a significant influence on the migration of the wetting front (P < 0.01). Coarse-textured soil produced faster migration in all three directions than fine-textured soil. The main reason was that the coarser the soil texture is, the larger the intergranular pores are, the stronger the permeability is, and the faster the wetting front moves. Furthermore, before the wetting front reached the surface, the wetting front of silt loam had only slightly different migration distances in the three directions, while the sandy loam had very different distances. For example, when the irrigation time was 15 h, the difference between the wetted Z and the wetted Y of the silt loam was only 2.5 cm, while that of the sandy loam was 25.7 cm. It should be that the higher the clay content is, the stronger is the capillarity and the more uniform the migration of the wetting front in all directions. In contrast, the rougher the soil texture is, the stronger the action of gravity is, resulting in the vertical downward wetting depth exceeding the vertical upward wetting height.

The initial moisture content of the soil had a great influence on the migration of the wetting front (P < 0.05), and the migration rate of the wetting front increased with the increase in θ_i in all three directions (Figure 5(b), 5(d) and 5(e)). This was because with the same soil texture, K(h) would increase with the increase in θ_i, while the increased K(h) would promote the flow of water through the pores of the soil. In addition, the larger θ_i is, the higher the water content is in the soil pores before irrigation begins, and less water is required to fill in the soil pores, thus accelerating the migration of the wetting front.

Figure 5(b), 5(f) and 5(g) show that the length of the line source had a significant influence on the migration of the wetting front (P < 0.05). This phenomenon is mainly reflected in the following aspects: (1) The initial values (L/2) of the wetted Z and the wetted Y increased with the increase in L; (2) the constant value (B − L/2) of the wetted Y decreased with the increase in L; (3) the migration rate in the X and Z directions increased with the increase in L. The diameter of the line source had a great influence on the migration of the wetting front in vertical line source irrigation (P < 0.05), and the migration rate of the wetting front increased with the increase in D in all three directions (Figure 5(b), 5(h) and 5(i)). D and L are the indexes of the infiltration channel of the emitter, which determine the infiltration interface of the emitter in size (πD²L/4). The larger D or L is, the larger is the area of the seepage surface, the more channels there are for water to enter the soil and the greater the seepage volume, which accelerates the migration of the wetting front.

By comparing and analysing Figure 5(b), 5(j) and 5(k), the buried depth of the line source had no significant influence on the migration of the wetting front (P > 0.05). When the value of B was small (B = 30 cm), the wetting front could easily reach the soil surface, thus increasing the evaporation of surface water. When B was large (B = 50 cm), there were risks of deep water leakage and surface water deficit. Therefore, the buried depth of the line source should be designed according to the soil texture, root distribution and farming conditions to improve the utilization efficiency of irrigation water.

Using MATLAB software, the contours were plotted of the soil water content of the results of the 11 sets of simulation schemes in Table 4 and are shown in Figure 6.

At the end of irrigation, there was little difference in the soil moisture content contour, which was an approximate ‘ellipsoid’ rotated around the line source. Due to a sufficient water supply, the soil moisture content near the line source was close to saturation. The moisture content of the wetting pattern gradually decreased outward, and the contour lines gradually became dense. The soil texture had a significant effect on the volume of the wetting pattern (P < 0.05); the more clay content in the soil, the smaller the soil wetting pattern, and the higher the average moisture content. The ‘focus’ of the soil moisture content contour moved downward with increasing of soil sand content.

The value of θ_i set in the study was between 50% and 70% of the field water capacity and the change in soil moisture content in this area was small, leading to no significant difference in the volume of the wetting pattern for the three types of soil (P > 0.05) (Figure 6(b), 6(d) and 6(e)). Generally, irrigation was needed when the soil moisture content reached 50%–70% of the field water capacity. The soil moisture content in this area had little influence on the shape and volume of the wetting pattern. In addition, there were variations in the moisture content at each point of the field soil layer, which made statistical analysis difficult. Therefore, it is suggested that the influence of the initial soil moisture content should not be considered in the design of the VLSI system.

The line source length had a significant effect on the volume of the wetting pattern (P < 0.05). With increasing line source length, the volume of the wetting pattern decreased gradually (Figure 6(b), 6(f) and 6(g)). However, the line source diameter had no significant effect on the volume of the wetting pattern (P > 0.05). With increasing line source diameter, the volume of the wetting pattern decreased slightly. This seemed to contradict the observation that the higher L and D were, the faster was the infiltration rate and the faster the migration of the wetting front. In fact, the larger L and D were, the faster the infiltration would be, resulting in a shorter irrigation duration under the same irrigation quota conditions, which produced a decrease rather than an increase in the wetting pattern. From a practical point of view, D and L should be as large as possible to accelerate the rate of water seepage and as small as possible to weaken the impact on the growth of the crop root system and reduce the possibility of surface soil moisture.

The burial depth of the line source directly changed the moisture distribution position of the wetting pattern, which had a significant effect on the volume of the wetting pattern (P < 0.05). With increasing B, the volume of the wetting pattern gradually increased (Figure 6(b), 6(j) and 6(k)). The difference in the wetting pattern with different burial depths of line sources was mainly reflected in the upper part of the line sources. When B was small, the wetting pattern was located in shallow soil. When the wetting pattern moved downward, it gradually moved away from the soil surface (B = 50 cm).

The calculated and measured values of the moving distance of the wetting front in the three directions from point C under VLSI are compared as shown in Figure 7.

The MAE, RMSE, PBIAS and NSE values for the measured and calculated values are presented in Table 8.

MAE, RMSE, and PBIAS values ranged from 0.770 to 2.310 L, 1.093 to 2.910 L, and 2.050% to 3.131%, respectively. Meanwhile, NSE values were very close to 1.0. The distance of the WFM in each direction of point C of the VLSI calculated by the formula is only slightly different from the measured value. This shows that the established prediction model can effectively calculate the moving distance of the wetting front.

In previous research (Fan et al. 2018b), the authors adopted HYDRUS-2D software and simulated the effects of ST, θ_i, B, D, and L on the cumulative infiltration of VLSI. Based on the Philip infiltration formula (Philip 1957), a calculation model for the cumulative infiltration of VLSI including the line source seepage area was established. However, the soil wetting characteristics (including WFM and SWD) have not been studied by either experimental methods or numerical simulations under the different influencing factors of VLSI to find their universal law. Based on previous research on the effects of ST, θ_i, B, D, and L on the cumulative infiltration of VLSI, this study qualitatively analysed the characteristics of soil wetting in VLSI under different influencing factors.

Based on the HYDRUS-2D model, a mathematical model of soil water movement in VLSI was established, and the reliability of the numerical simulation method was verified. On this basis, the effects of ST, θ_i, L, D and B on the characteristics of the wetting pattern were simulated and analysed. The following conclusions were drawn: the dynamic change law and trend of wetting bodies under different influencing factors were basically the same. During the early stage of irrigation, the WFM rate was in the order of horizontal > vertical downward > vertical upward; as time passed, the WFM rate showed a trend of vertical downward > horizontal direction > vertical upward. The ST had a significant influence on the WFM distance. The coarser the grain content of the soil, the faster the wetting front moved. The values of L and D had little influence on the WFM distance. The larger L and D were, the faster the wetting front moved; θ_i and B had a weak influence on the WFM distance. Based on a comprehensive analysis, ST was the dominant factor affecting the distance of the wetting front of the VLSI. Under different influencing factors, the shapes of the contour of soil moisture content were slightly different, but all were approximate ‘ellipsoids’. The soil moisture content near the line source was close to saturation, the soil moisture content gradually decreased from the inside to the outside, and the contour lines gradually became dense. When the same I_a was used, the wetting body was mainly affected by ST, L and B; θ_i and D had little effect on it.

In this study, the WFM and SWD of the VLSI under different influencing factors were determined, and a prediction model of the WFM was established. Comprehensive optimization of the irrigation parameters to reduce surface evaporation and deep leakage of soil water, and effective matching of the soil water with the plant roots to improve the benefits of the VLSI system, still require further research.

This research was supported by the National Natural Science Foundation of China (No. 51969013), the Natural Science Foundation of Gansu Province, China (No. 18JR3RA144), and Hongliu Supporting Discipline of Lanzhou University of Technology.

All relevant data are included in the paper or its Supplementary Information.