Numerical investigation of wetting front migration and soil water distribution under vertical line source irrigation with different in ﬂ uencing factors

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 veri ﬁ cation 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 about 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 in ﬂ uenced 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 in ﬂ uencing 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 decreases outward, and the contour lines were gradually dense. According to the simulation results, a prediction model of multiple factors in ﬂ uencing 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. the migration process of the vertical line source irrigation wetting front was established. The model was validated by experimental data, and the prediction effect was good.


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. ; Du et al. ).
However, rainfall in the area is limited, surface evaporation is strong, and fruit production depends to a large extent on irrigation (Shen et al. ; Guo & Shen ). Surface irrigation is a common irrigation method in orchards that is simple to operate but easily wastes water (Kang et  The shape and size of the soil wetting body determine the plant growth status and efficiency of field water use.

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 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 observate 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 irri- 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 a 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 1 day to evenly distribute the soil moisture.

Experimental scheme
Referring to the case of VLSI in the field, 4-8 irrigation emit-   (Table 2). Furthermore, sand and silt loam were used to verify the HYDRUS simulation results, while loam and sandy loam were used to evaluate of the WFM model.

Governing equation
It is assumed that the soil is an isotropic homogeneous porous medium and that there is no air resistance that affects water movement. Ignoring the influence of the temperature field, the diffusion of water from the seepage hole into the soil and then to the surroundings is a spatial three-dimensional infiltration process. The basic equation used in this study is the Richards equation, as follows: where 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.
The soil water retention was described using the van Genuchten equation (Van Genuchten ).
where S e is the effective degree of saturation; θ r is the residual water contents, 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; Silt ( 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  and was fitted using the RETC (Van Genuchten et al. ). K s was determined by a constant head permeameter (Provenzano ). The van Genuchten-Mualem model parameters thus obtained are shown in Table 3. In summary, the initial conditions can be expressed as:

Initial and boundary conditions
The boundary conditions can be expressed as:

Model solving method
Numerical simulation was performed using HYDRUS-2D (Šimunek 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

Empirical model description
The shape of the soil wetting body under VLSI was approximately an ellipsoid surrounding the line source (Zeng et al.

;
Li & Wang ). Therefore, the eigenvalues of WFM in  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
The mean absolute error (MAE), root mean square error (RMSE), percent deviation (PBIAS) and Nash efficiency coefficient (NSE) were used to statistically analyse the agreement between the calculated and measured values of the wetting pattern dimensions. These statistical parameters are defined as follows: where M i is the ith measured value; S i is the ith 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, eleven 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.
The van Genuchten-Mualem model parameters for different STs were taken from Carsel and Parrish (1988) data, as shown in Table 5.
where θ f is the field capacity, cm 3 /cm 3 ; and θ w is the wilting coefficient, cm 3 /cm 3 . Nagličet al. () showed that θ w can be expressed by θ r , i.e., Combining (13) and (14), available   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

Wetting front position during irrigation
According to the results of 11 sets of simulation schemes in Table 4, the relationship between the moving distance and time of the wetting front in 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  Table 4, respectively.
Notes: X ¼ Horizontal wetted width from point C; Y ¼ vertical wetted height from point C; Z ¼ vertical wetted depth from point C.
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 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 gravity action, 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 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. The larger D or L is, the larger the area of the seepage surface, the more channels 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 Soil water distribution at the end of irrigation Using MATLAB software, the contours of the soil water content of the results of the 11 sets of simulation schemes are plotted 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 The fitting coefficient and the determination coefficient of fitting regression line are listed in Table 7.  Table 8