## Abstract

Vertical line source irrigation (VLSI) is an underground irrigation method suitable for deep-rooted plants. Understanding the characteristics of the soil wetting body of the VLSI was the key to designing this irrigation system. On the basis of experimental verification of the reliability of the HYDRUS simulation results of VLSI under the conditions of soil texture (ST), initial water content (*θ*_{i}), line source buried depth (*B*), line source diameter (*D*) and line source length (*L*), numerical studies of the migration law of the wetting front of VLSI and the distribution characteristics of soil moisture were performed. The wetting front migration (WFM) was mainly influenced by ST, *θ*_{i}, *D* and *L* (*P* < 0.05), while *B* had little effect on WFM (*P* > 0.05). The shape of the soil wetting body changed little under different influencing factors. The water content contour was approximately ‘ellipsoidal’ around the line source. The soil moisture near the line source was close to the saturated moisture content. The moisture content around the line source gradually decreased outward, and the contour lines gradually became dense. According to the simulation results, a prediction model of multiple factors influencing the migration process of the VLSI wetting front was established. The predicted value was in good agreement with the measured value. The results of this research could provide a theoretical basis for further optimizing the combination of VLSI and irrigation elements.

## HIGHLIGHTS

Based on HYDRUS-2D, a multi-factor prediction model was established for the wet front migration process in vertical source irrigation.

The validity of the model was verified by experimental data, and the prediction effect was good.

Numerical studies of the migration law of the wetting front of the vertical line source irrigation and the distribution characteristics of soil moisture were performed.

The wetting front migration was mainly influenced by soil texture, initial water content, line source diameter and line source length (

*P*< 0.05), while line source buried depth had little effect on it (*P*> 0.05).The shapes of the contours of soil moisture content were slightly different, and were all approximate 'ellipsoids' with the axis of the line source.

## ABBREVIATIONS AND SYMBOLS

*α*an empirical parameter inversely related to the air entry value (cm

^{−1})*r*radial coordinate

*z*vertical coordinate

*θ*soil water content (cm

^{3}/cm^{3})*θ*_{s}saturated water content

*θ*_{r}residual water content (cm

^{3}/cm^{3})*θ*_{i}initial water content

*θ*_{f}field capacity (cm

^{3}/cm^{3})*θ*_{w}wilting coefficient (cm

^{3}/cm^{3})*K*(*h*)unsaturated hydraulic conductivity (cm/min)

*K*_{s}saturated hydraulic conductivity (cm/min)

*Se*the effective degree of saturation

*M*_{i}the

*i*th measured value*S*_{i}the

*i*th 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

## INTRODUCTION

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.

## MATERIALS AND METHODS

### Experimental design of VLSI

#### Laboratory experiment

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.

#### Soil sample preparation

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.

ST . | Soil particle size distribution (%) . | γ (g·cm^{−3})
. | θ_{s} (cm^{3}·cm^{−3})
. | θ_{f} (cm^{3}·cm^{−3})
. | K_{s} (cm·min^{−1})
. | ||
---|---|---|---|---|---|---|---|

Sand (0.02–2.00 mm) . | Silt (0.002–0.02 mm) . | Clay (0–0.002 mm) . | |||||

Silt loam | 27.39 | 49.47 | 23.14 | 1.30 | 0.496 | 0.319 | 0.014 |

Loam | 49.72 | 11.66 | 38.62 | 1.35 | 0.479 | 0.271 | 0.018 |

Sandy loam | 78.38 | 12.64 | 8.97 | 1.45 | 0.447 | 0.158 | 0.074 |

Sand | 93.76 | 5.92 | 1.32 | 1.55 | 0.428 | 0.089 | 0.421 |

ST . | Soil particle size distribution (%) . | γ (g·cm^{−3})
. | θ_{s} (cm^{3}·cm^{−3})
. | θ_{f} (cm^{3}·cm^{−3})
. | K_{s} (cm·min^{−1})
. | ||
---|---|---|---|---|---|---|---|

Sand (0.02–2.00 mm) . | Silt (0.002–0.02 mm) . | Clay (0–0.002 mm) . | |||||

Silt loam | 27.39 | 49.47 | 23.14 | 1.30 | 0.496 | 0.319 | 0.014 |

Loam | 49.72 | 11.66 | 38.62 | 1.35 | 0.479 | 0.271 | 0.018 |

Sandy loam | 78.38 | 12.64 | 8.97 | 1.45 | 0.447 | 0.158 | 0.074 |

Sand | 93.76 | 5.92 | 1.32 | 1.55 | 0.428 | 0.089 | 0.421 |

#### Experimental scheme

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.

Test number . | ST . | θ_{i}/(cm^{3}·cm^{−3})
. | D/cm
. | L/cm
. | B/cm
. | I_{a}/L
. |
---|---|---|---|---|---|---|

T1 | Sand | 50%θ_{f} | 4 | 15 | 30 | 20 |

T2 | Silt loam | 70%θ_{f} | 6 | 20 | 40 | 20 |

T3 | Loam | 50%θ_{f} | 2 | 10 | 30 | 40 |

T4 | Loam | 50%θ_{f} | 6 | 30 | 50 | 40 |

T5 | Loam | 70%θ_{f} | 6 | 10 | 30 | 40 |

T6 | Loam | 70%θ_{f} | 2 | 30 | 50 | 40 |

T7 | Sandy loam | 50%θ_{f} | 2 | 10 | 30 | 40 |

T8 | Sandy loam | 50%θ_{f} | 6 | 30 | 50 | 40 |

T9 | Sandy loam | 70%θ_{f} | 6 | 10 | 30 | 40 |

T10 | Sandy loam | 70%θ_{f} | 2 | 30 | 50 | 40 |

Test number . | ST . | θ_{i}/(cm^{3}·cm^{−3})
. | D/cm
. | L/cm
. | B/cm
. | I_{a}/L
. |
---|---|---|---|---|---|---|

T1 | Sand | 50%θ_{f} | 4 | 15 | 30 | 20 |

T2 | Silt loam | 70%θ_{f} | 6 | 20 | 40 | 20 |

T3 | Loam | 50%θ_{f} | 2 | 10 | 30 | 40 |

T4 | Loam | 50%θ_{f} | 6 | 30 | 50 | 40 |

T5 | Loam | 70%θ_{f} | 6 | 10 | 30 | 40 |

T6 | Loam | 70%θ_{f} | 2 | 30 | 50 | 40 |

T7 | Sandy loam | 50%θ_{f} | 2 | 10 | 30 | 40 |

T8 | Sandy loam | 50%θ_{f} | 6 | 30 | 50 | 40 |

T9 | Sandy loam | 70%θ_{f} | 6 | 10 | 30 | 40 |

T10 | Sandy loam | 70%θ_{f} | 2 | 30 | 50 | 40 |

### Numerical simulation

#### Governing equation

*r*is the radial coordinate;

*z*is the vertical coordinate, specifying that the positive

*z*direction is toward the bottom;

*θ*is the soil water content, cm

^{3}/cm

^{3};

*h*is the pressure head, cm;

*t*is the infiltration time, min; and

*K*(

*h*) is the unsaturated hydraulic conductivity, cm/min.

*S*

_{e}is the effective degree of saturation;

*θ*

_{r}is the residual water content, cm

^{3}/cm

^{3};

*α*is an empirical parameter that is inversely proportional to the intake air value (cm

^{−1}); and

*n*and

*m*are empirical constants affecting the shape of the soil moisture characteristic curve, with

*n*> 1 and

*m*= 1−1/

*n*.

*K*(

*h*) was described using the Mualem model (Mualem 1976):where

*K*

_{s}is the saturated hydraulic conductivity, cm/min.

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.

Soil texture . | Residual water content θ_{r}/(cm^{3}·cm^{−3})
. | Saturated water content θ_{s}/(cm^{3}·cm^{−3})
. | Reciprocal of air entry value α/(cm^{−1})
. | Experienced parameter n
. | Saturated hydraulic conductivity K_{s}/(cm·min^{−1})
. |
---|---|---|---|---|---|

Sand | 0.02 | 0.43 | 0.162 | 2.24 | 0.421 |

Silt loam | 0.02 | 0.50 | 0.014 | 1.51 | 0.014 |

Soil texture . | Residual water content θ_{r}/(cm^{3}·cm^{−3})
. | Saturated water content θ_{s}/(cm^{3}·cm^{−3})
. | Reciprocal of air entry value α/(cm^{−1})
. | Experienced parameter n
. | Saturated hydraulic conductivity K_{s}/(cm·min^{−1})
. |
---|---|---|---|---|---|

Sand | 0.02 | 0.43 | 0.162 | 2.24 | 0.421 |

Silt loam | 0.02 | 0.50 | 0.014 | 1.51 | 0.014 |

#### Initial and boundary conditions

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

#### Model-solving method

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.

### Empirical model description

*et al.*2010; Li & Wang 2011). Therefore, the eigenvalues of WFM in the three directions (vertical upward, horizontal downward and vertical downward from point C) were selected to outline the wetting patterns. Point C was the centre of the line source (Figure 2). An empirical model was proposed to estimate the size of the WFM eigenvalues during VLSI, which considered the influence of

*K*

_{s},

*θ*

_{s}−

*θ*

_{i},

*L*and

*D*. The general form of the proposed model was assumed to be as shown in Equations (6)–(8):where

*X*is the horizontal wetting width from point C, cm;

*Y*is the vertical upward wetting height from point C, cm;

*Z*is the vertical downward wetting depth from point C, cm; and

*a*to

*a*

_{5},

*b*to

*b*

_{5}and

*c*to

*c*

_{5}are empirical coefficients.

### Statistical analysis

*M*is the

_{i}*i*th measured value;

*S*is the

_{i}*i*th simulated value; and

*N*is the total number. The closer MAE and RMSE are to 0, −10 ≤ PBIAS ≤ +10, and the closer NSE is to 1, the smaller the difference between the calculated and measured values is, and the better the agreement.

### Simulation schemes for different influencing factors

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.

Scenario . | Soil texture . | Initial water content θ_{i}/(cm^{3}·cm^{−3})
. | Line source length L/cm
. | Line source diameter D/cm
. | Line source buried depth B/cm
. | Irrigation amount I_{a}/L
. |
---|---|---|---|---|---|---|

Case a | Silt loam | 60% θ_{f} | 20 | 4 | 40 | 40 |

Case b | Loam | 60% θ_{f} | 20 | 4 | 40 | 40 |

Case c | Sandy loam | 60% θ_{f} | 20 | 4 | 40 | 40 |

Case d | Loam | 50% θ_{f} | 20 | 4 | 40 | 40 |

Case e | Loam | 70% θ_{f} | 20 | 4 | 40 | 40 |

Case f | Loam | 60% θ_{f} | 10 | 4 | 40 | 40 |

Case g | Loam | 60% θ_{f} | 30 | 4 | 40 | 40 |

Case h | Loam | 60% θ_{f} | 20 | 2 | 40 | 40 |

Case i | Loam | 60% θ_{f} | 20 | 6 | 40 | 40 |

Case j | Loam | 60% θ_{f} | 20 | 4 | 30 | 40 |

Case k | Loam | 60% θ_{f} | 20 | 4 | 50 | 40 |

Scenario . | Soil texture . | Initial water content θ_{i}/(cm^{3}·cm^{−3})
. | Line source length L/cm
. | Line source diameter D/cm
. | Line source buried depth B/cm
. | Irrigation amount I_{a}/L
. |
---|---|---|---|---|---|---|

Case a | Silt loam | 60% θ_{f} | 20 | 4 | 40 | 40 |

Case b | Loam | 60% θ_{f} | 20 | 4 | 40 | 40 |

Case c | Sandy loam | 60% θ_{f} | 20 | 4 | 40 | 40 |

Case d | Loam | 50% θ_{f} | 20 | 4 | 40 | 40 |

Case e | Loam | 70% θ_{f} | 20 | 4 | 40 | 40 |

Case f | Loam | 60% θ_{f} | 10 | 4 | 40 | 40 |

Case g | Loam | 60% θ_{f} | 30 | 4 | 40 | 40 |

Case h | Loam | 60% θ_{f} | 20 | 2 | 40 | 40 |

Case i | Loam | 60% θ_{f} | 20 | 6 | 40 | 40 |

Case j | Loam | 60% θ_{f} | 20 | 4 | 30 | 40 |

Case k | Loam | 60% θ_{f} | 20 | 4 | 50 | 40 |

*θ*

_{f}was obtained using a prediction model established by Rab

*et al.*(2011). The specific expression is:where

*θ*

_{f}is the field capacity, cm

^{3}/cm

^{3}; and

*θ*

_{w}is the wilting coefficient, cm

^{3}/cm

^{3}.

ST . | Residual water content θ_{r}/(cm^{3}·cm^{−3})
. | Saturated water content θ_{s}/(cm^{3}·cm^{−3})
. | Reciprocal of air entry value α/(cm^{−1})
. | Experienced parameter n/(−)
. | Saturated hydraulic conductivity K/(cm·min_{s}^{−1})
. |
---|---|---|---|---|---|

Silt loam | 0.067 | 0.45 | 0.020 | 1.41 | 0.0075 |

Loam | 0.078 | 0.43 | 0.036 | 1.56 | 0.0173 |

Sandy loam | 0.065 | 0.41 | 0.075 | 1.89 | 0.0737 |

ST . | Residual water content θ_{r}/(cm^{3}·cm^{−3})
. | Saturated water content θ_{s}/(cm^{3}·cm^{−3})
. | Reciprocal of air entry value α/(cm^{−1})
. | Experienced parameter n/(−)
. | Saturated hydraulic conductivity K/(cm·min_{s}^{−1})
. |
---|---|---|---|---|---|

Silt loam | 0.067 | 0.45 | 0.020 | 1.41 | 0.0075 |

Loam | 0.078 | 0.43 | 0.036 | 1.56 | 0.0173 |

Sandy loam | 0.065 | 0.41 | 0.075 | 1.89 | 0.0737 |

*et al.*(2014) showed that

*θ*

_{w}can be expressed by

*θ*

_{r}, i.e.,Combining (13) and (14), available

## DISCUSSION AND RESULTS

### Comparison of HYDRUS simulations with experimental observations

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.

ST . | Wetting front migration distance . | Soil water content distribution . | ||||||
---|---|---|---|---|---|---|---|---|

MAE (cm) . | RMSE (cm) . | PBIAS (%) . | NSE (-) . | MAE (cm^{3}/cm^{3})
. | RMSE (cm^{3}/cm^{3})
. | PBIAS (%) . | NSE (−) . | |

Sand | 0.820 | 0.963 | 1.481 | 0.981 | 0.019 | 0.025 | −4.033 | 0.928 |

Silt loam | 0.743 | 1.021 | 2.420 | 0.984 | 0.015 | 0.020 | −4.087 | 0.963 |

ST . | Wetting front migration distance . | Soil water content distribution . | ||||||
---|---|---|---|---|---|---|---|---|

MAE (cm) . | RMSE (cm) . | PBIAS (%) . | NSE (-) . | MAE (cm^{3}/cm^{3})
. | RMSE (cm^{3}/cm^{3})
. | PBIAS (%) . | NSE (−) . | |

Sand | 0.820 | 0.963 | 1.481 | 0.981 | 0.019 | 0.025 | −4.033 | 0.928 |

Silt loam | 0.743 | 1.021 | 2.420 | 0.984 | 0.015 | 0.020 | −4.087 | 0.963 |

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.

### Wetting front position during irrigation

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*^{2}*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.

### Soil water distribution at the end of irrigation

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

### An empirical model for predicting wetting front migration

#### Establishment of empirical model

*K*

_{s},

*D*,

*L*, and (

*θ*

_{s}−

*θ*

_{i}). Fitting parameters a–a5, b–b5 and c–c5 in the different directions were obtained by Origin9.0 fitting. The determination coefficients of the fitting regression line were 0.99, 0.98, 0.98, respectively, all close to 1, indicating that the fitting effect was good.

Horizontal of point C . | Upward of point C . | Downward of point C . | |||
---|---|---|---|---|---|

a | 3.8 | b | 3.9 | c | 8.1 |

a1 | 0.25 | b1 | 0.17 | c1 | 0.52 |

a2 | 0.14 | b2 | 0.16 | c2 | 0.35 |

a3 | 0.14 | b3 | −0.04 | c3 | 0.08 |

a4 | 0.08 | b4 | 0.08 | c4 | 0.18 |

a5 | 0.34 | b5 | 0.32 | c5 | 0.43 |

R^{2} | 0.99 | R^{2} | 0.98 | R^{2} | 0.98 |

Horizontal of point C . | Upward of point C . | Downward of point C . | |||
---|---|---|---|---|---|

a | 3.8 | b | 3.9 | c | 8.1 |

a1 | 0.25 | b1 | 0.17 | c1 | 0.52 |

a2 | 0.14 | b2 | 0.16 | c2 | 0.35 |

a3 | 0.14 | b3 | −0.04 | c3 | 0.08 |

a4 | 0.08 | b4 | 0.08 | c4 | 0.18 |

a5 | 0.34 | b5 | 0.32 | c5 | 0.43 |

R^{2} | 0.99 | R^{2} | 0.98 | R^{2} | 0.98 |

#### Evaluation of empirical model

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.

WFC . | MAE (cm) . | RMAE (cm) . | PBIAS (%) . | NSE (−) . |
---|---|---|---|---|

Horizontal of point C | 1.018 | 1.607 | 2.050 | 0.979 |

Upward of point C | 0.770 | 1.093 | 3.131 | 0.978 |

Downward of point C | 2.310 | 2.910 | 2.374 | 0.959 |

WFC . | MAE (cm) . | RMAE (cm) . | PBIAS (%) . | NSE (−) . |
---|---|---|---|---|

Horizontal of point C | 1.018 | 1.607 | 2.050 | 0.979 |

Upward of point C | 0.770 | 1.093 | 3.131 | 0.978 |

Downward of point C | 2.310 | 2.910 | 2.374 | 0.959 |

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.

## CONCLUSIONS

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.

## ACKNOWLEDGEMENTS

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.

## DATA AVAILABILITY STATEMENT

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

## REFERENCES

*Research Report No. EPA/600/2-91/065*