Numerical simulation of soil water movement by gravity subsurface hole irrigation Open Access

Dryland agriculture is an essential type of agricultural production in China, and the Loess Plateau is located within the heartland of dryland agriculture in China (Li et al. 2016). Drought, water shortages, and soil erosion are the two bottlenecks limiting the sustainable development of the Loess Plateau and the fundamental reasons for the ecological fragility of the region (Gao et al. 2016). Regulating rainfall–runoff and developing water-saving irrigation technologies are the primary methods to solve the abovementioned issues (Zhao et al. 2009). Rainwater harvesting and storage is an essential measure of rainfall–runoff and hydraulic control engineering. China's Loess Plateau region has a long history of rainwater harvesting and utilization, has accumulated operable rainwater harvesting and utilization techniques, and has built multiple rainwater storage facilities, such as rainwater harvesting cellars and rainwater harvesting ponds (Zhao et al. 2009). Due to the constraints of topography, rainfall resources, and economic conditions, it is no longer practical to use traditional surface irrigation methods for adequate irrigation. Water-saving irrigation methods, such as preserving critical water or life-saving water, can maximize the benefits of irrigation water (Bhatnagar & Srivastava 2003). Currently, water-saving irrigation technology in arid and semiarid regions mainly focuses on sprinkler irrigation (De Wrachien & Lorenzini 2006), surface drip irrigation (Kwon et al. 2020), and subsurface drip irrigation (Kanda et al. 2020b). However, sprinkler, surface drip, and subsurface drip irrigation require relatively large amounts of hydraulic power for irrigation water supply. It is difficult to provide such large amounts of hydraulic power equipment in hilly orchards with poor infrastructure (Bautista-Capetillo et al. 2012). Using manual or water diversion facilities to distribute rainwater collected in rainwater harvesting cellars or ponds and use for sprinkler or drip irrigation will increase operating costs. Therefore, it is necessary to develop a rainwater harvesting and irrigation technology with low kinetic energy, a water supply, and low operating costs – gravity subsurface hole irrigation (Figure 1).

In the rainwater collection of the gravity subsurface hole irrigation system, the rainwater flows into the infiltrating pipes after being filtered by the water holes through the ground. It is introduced into the water delivery pipeline through the infiltrating pipes' vertical pipes, and it is discharged into the buffer tank through the delivery pipeline. After the buffer tank is collected and precipitated, it flows into the cistern by sluice gate 1 (note: the ventilation hole is a device for connecting outside air to ensure air pressure balance). In the process of irrigation, the water is lifted from the cistern into the fertilizer tank through the manual pressure well or an electric pump, flows to the buffer tank through sluice gate 2, then flows into the infiltrating pipe through the water delivery pipeline and vertical pipe, and then penetrates the crop root area through the infiltration surface (infiltrating hole) at the bottom of the infiltrating pipe. In the above two processes, the water in the buffer tank flows into the cistern during rainwater collection. The water in the fertilizer tank flows into the infiltrating pipe during irrigation, which does not need to provide kinetic energy manually. The flow of water depends entirely on the kinetic energy converted from the gravitational potential energy of the water. Therefore, the rainwater collection and irrigation process of gravity subsurface hole irrigation almost all rely on the gravitational potential energy of the system's water flow. This belongs to the semi-automatic ‘active drought resistance’ mode, which aligns with the development concept of rainwater collection and irrigation in dry farming areas. Gravity subsurface hole irrigation is a new technology for rainwater harvesting and irrigation. It combines subsurface drip irrigation and film hole irrigation methods and has similarities with point source irrigation. Both are part infiltration irrigation, but their soil water migration law is different from subsurface drip irrigation (Kanda et al. 2020b) and film hole irrigation (Zhong et al. 2020). Therefore, the study of gravity subsurface hole irrigation can draw on the research ideas of subsurface drip irrigation, film hole irrigation, and other irrigation technologies. Analysis of the characteristics of soil water infiltration patterns determines the dominant factors affecting gravity subsurface hole irrigation. It provides a scientific basis for the design and operation of this irrigation project.

Usually, experiments or numerical simulations can determine the soil water transport pattern around the irrigator. The experimental method is the primary method for studying soil water transport patterns under different influencing factors, but it is also costly and time-consuming (Kanda et al. 2020a). Numerical simulations allow the analysis of soil water movement under different soil characteristics and design parameters, providing a practical and convenient means of determining the appropriate technical parameters for irrigation and optimizing the operation of irrigation systems.

HYDRUS results can better reflect the fundamental laws of soil water movement in terms of point and line source infiltration (Šimůnek et al. 2016). Skaggs et al. (2004) used HYDRUS-2D numerical values to study the soil wetting pattern size and water distribution under drip irrigation and proved that HYDRUS-2D could be used to investigate and design drip irrigation management practices. Based on the reliability of the HYDRUS-2D model, Provenzano (2007) simulated and analyzed the change process of wetted soil bodies under subsurface drip irrigation. The results show that the size of the moist body is a function of time and initial moisture content. Fan et al. (2019) used HYDRUS-2D/3D simulations to investigate the effects of soil texture, initial moisture content, film hole diameter, and depth of water accumulation above the hole on the cumulative infiltration volume (I) and wetting front migration distance (W) during the free infiltration phase of film hole irrigation. I and W estimation models were established, and the reliability of the estimation model was verified using soil box tests. In addition, Naglič et al. (2014) used the HYDRUS-2D/3D model to investigate the factors influencing wetting pattern size for surface drip irrigation, mainly soil texture, drip head flow rate, and initial moisture content, to develop a prediction model for soil wetting pattern size and to assess the accuracy of the model.

Based on the above analysis, this paper will study the irrigation effect of gravity subsurface hole irrigation. Using the combination of HYDRUS-2D/3D numerical simulation and indoor soil box test, the numerical model of soil water movement is established to analyze different soil saturated hydraulic conductivity (K_s), infiltrating hole diameter (D), infiltrating pipe depth (B), and matric potential (Ψ_m). The variation law of cumulative infiltration and wetting front under the combination provides a scientific basis for the rational design and popularization of gravity subsurface hole irrigation.

Figure 2 shows a schematic diagram of the boundary condition calculation area used to simulate the different scenarios in this study.

In all simulation scenarios, the changes in Ψ_m were not considered and were set according to the Ψ_m. The flux boundary was considered on O–O’ because there is no obstacle at the bottom of the infiltrating pipe to block the flow. The bottom boundary adopts free drainage, and the other adopts a zero-flux boundary. The water content of the soil at the bottom of the weep holes was infinitely close to the same value, so data on the wetting front migration distance was recorded from point O.

The simulation used the Galerkin finite element method to spatially discretize the soil profile, simplify the calculation area into rectangular elements, and determine the finite element calculation domain. The transport domain was 200 cm deep (depth) and 80 cm wide (radius), i.e., large enough so that an overlap in water content profiles from neighboring emitters (on other laterals) do not have to be considered. In the process of simulation, the spatial step is 2 cm and the time step is 1 min.

To study the influence degree of factors that may affect gravity subsurface hole irrigation, this paper mainly adopts the numerical method. Using the same irrigation quota (I_q) and referring to the research results of Zhang et al. (2010), the irrigation water volume of each infiltrating pipe is taken as 40 L. Thirty-five scenarios are laid out using single-factor analysis, and the setting of soil matric potential is 30%, 50% and 70% of the wetting soil coefficient (Naglič et al. 2014) (Table 1) to simulate the effects of K_s, D, B, and Ψ_m on I and W for gravity subsurface hole irrigation. Parameters of the soil texture VG-M model were taken from the data of Carsel & Parrish (1988), as shown in Table 2.

The purpose of the indoor test is to verify the accuracy of the numerical simulation and the reliability of the empirical model proposed later. The experimental device consists of three parts: a soil tank, a marionette bottle, and the irrigator, as shown in Figure 3. The soil tank was made of 10-mm-thick transparent Plexiglas. To ensure that no intersection occurred during the test, the internal dimensions of the soil tank were designed as 60 × 60 × 100 cm (length × width × height). Several ventilation holes (2 mm diameter) were set in the bottom of the soil tank to prevent the occurrence of air resistance. The diameter of the marionette bottle was 10 cm, and the height was 100 cm. Before the test, water was added to the test soil sample at the set initial water content, mixed evenly, sealed with plastic film, and left to stand for one day. After the soil moisture was uniformly distributed, it was loaded into the soil tank in layers (5 cm) according to the design bulk density. The unit weight of the test soil was as far as possible the same to obtain a uniform soil profile. To facilitate the observation of the soil wetting front transport process, 1/4 of the infiltrating pipe was placed at the corner of the soil tank to ensure that the infiltrating pipe wall was in close contact with the soil. Then the infiltration test was carried out the following day. During the test, the mariotte bottle provided a constant water head. The data were observed and recorded according to the time interval of first dense and then sparse intervals. The readings of the mariotte bottles at different times were the cumulative infiltration of soil. At the same time, the migration distance map of the wetting front was drawn on the transparent soil box with a marker to record the gradual change process of the wetting front. The water supply was stopped when the infiltration rate reached a constant level, and the experiment was completed. Each test was repeated three times to eliminate as many errors as possible.

Soil for the experiment was taken from orchards around Lanzhou at a depth of 0 to 40 cm. Soils used for the experiment were air-dried, rolled, uniformly mixed, and sieved through a 2 mm sieve. The soil particle size distribution in the test soil box was measured using a laser particle size analyzer (clay accounted for 9.85%; silt accounted for 68.45%; sand accounted for 21.70%). The results indicated that the soil was a sandy clay loam. K_s measured by the constant head method was 0.0092 cm min⁻¹, soil bulk density measured by the ring sampler method was 1.35 g cm⁻³, and initial water content measured by the drying method was 0.01 cm³ cm⁻³. The VG model converted the soil moisture content and matric potential to obtain the value of soil matric potential. Three groups were set in the test scheme: treatment 1: D = 12.8 cm, B = 20 cm, Ψ_m = −4,009 cm; treatment 2: D = 11.0 cm, B = 25 cm, Ψ_m = −7,527 cm; treatment 3: D = 14.4 cm, B = 20 cm, Ψ_m = −7,527 cm.

Numerical analysis of the influence of the four influencing factors K_s, D, B, and Ψ_m on I was performed, and the curve of I versus T was drawn, as shown in Figure 4.

As shown in Figure 4, I increases with increasing T for different combinations of influencing factors (K_s, D, B, and Ψ_m). In a double logarithmic coordinate system, I is linear concerning T, indicating a good power function relationship between the two. K_s has an apparent effect on I, with higher values of K_s resulting in faster soil water infiltration and shorter time to reach the same amount of irrigation under the same D, B, and Ψ_m conditions (Figure 4(a)). For example, it takes 49.33 hours for a clay loam with K_s of 0.0043 cm min⁻¹ (I_q = 40 L) to complete irrigation, while it takes only 1.07 hours for loamy sand (K_s = 0.2432 cm min⁻¹). The times spent on silty loam (K_s = 0.0075 cm min⁻¹), loam (K_s = 0.0173 cm min⁻¹) and sandy loam (K_s = 0.0737 cm min⁻¹) were 27.00, 12.20 and 2.98 hr, respectively. Because K_s is an important physical index reflecting soil infiltration performance, the greater the K_s value, the stronger the soil hydraulic conductivity and the stronger the soil infiltration performance. D significantly influences I, which increases with increasing D at the same moment under the same K_s, B, and Ψ_m conditions (Figure 4(b)). For example, at T = 8 hr, for a small infiltrating hole (D = 8 cm), I was only 18.00 L. When D was increased to 12 and 16 cm, I increased by 49.44% and 97.22% to 26.90 L and 35.50 L, respectively. This increase is because, in the process of gravity subsurface hole irrigation, the bottom surface of the infiltrating hole serves as the infiltration interface for water to enter the soil. The larger D is, the larger the area of the infiltrating hole, the more channels for water to infiltrate into the soil, the faster the infiltration rate, and the more water that enters the soil per unit of time. B substantially influences I, which increases with B at the same moment under the same K_s, D, and Ψ_m conditions (Figure 4(c)). For example, at T = 8 hr, I was only 21.40 L for a low burial depth (B = 30 cm) and increased by 25.70% and 47.20% to 26.90 and 31.50 when B was increased to 40 and 50 cm, respectively. This was due to the irrigation water in the infiltrating pipe, which was self-pressurized from the reservoir. The pressure head at the infiltrating hole was approximately equal to B. Therefore, as B increases, the pressure head increases, resulting in an increase in the pressure potential at the bottom of the infiltrating pipe and thus an increase in I. The effect of Ψ_m on I was small, with a slight increase in I with decreasing Ψ_m at the same time for the same K_s, D, and B conditions (Figure 4(d)). For example, at T = 12 hr, I corresponding to the three substrate potentials (Ψ_m = −3,750 cm, −5,000 cm, and −7,500 cm) was 39.30, 39.40 and 39.50 L, respectively, and they increased by 0.1 L in turn, a small increase. An explanation for this small increase may be that, for the same soil texture, the soil water potential gradient increases with decreasing Ψ_m, which has a facilitating effect on soil water infiltration. However, a decrease in Ψ_m causes a decrease in the unsaturated hydraulic conductivity of the soil and then hinders soil water infiltration. The combined effect of unsaturated soil hydraulic conductivity and the soil water potential gradient weakens the influence of matric potential on soil accumulated infiltration.

A qualitative analysis of the factors influencing I shows that K_s, D, B, Ψ_m, and T all affect I. I increases with increasing K_s, D, and B but decreases slightly with increasing Ψ_m. To accurately describe the quantitative relationships between K_s, D, B, Ψ_m, and T and I, the empirical constants a and b and the values of the coefficient of determination R² were obtained by fitting, based on the results of numerical simulations, using Equation (4), as shown in Table 1.

Using Equation (9), simulation results were fitted to obtain empirical constant a′ coefficient of determination R² values (Table 1). The value of a′ changes considerably, and the method of calculating the average value to determine a′ was not indicative of the actual situation. Further analysis revealed that a′ was influenced by K_s, D, B, and Ψ_m, increasing with K_s, D, and B and decreasing with Ψ_m, and there was a good power function relationship between K_s, D, B, and Ψ_m and a′ (Figure 5).

The Equation (11) coefficient of determination is 0.998, which indicates that the power function continuous multiplication model holds and fits well. As the five influencing factors on the right-hand side of Equation (11) are not uniform in units and are not of the same order of magnitude, it was impossible to directly determine the degree of influence of each on I through the regression coefficients. Standardized regression analysis using SPSS statistical software yielded standardized regression coefficients of 0.31, 0.14, 0.09, −0.04, and 0.78 for K_s, D, B, Ψ_m, and T, respectively, indicating that they all influenced I and the degrees of influence were T > K_s > D > B > Ψ_m in descending order.

Equation (11) was fitted from simulation data and needed to be further verified for accuracy in practical application scenarios. The reliability of Equation (11) was evaluated by three sets of indoor soil tank experiments, and the measured values were plotted against the predicted values, as shown in Figure 6.

As shown in Figure 6, the model estimates of I were consistent with the measured values of the experiment, and the P-value calculated using the t-test, P > 0.05, indicates that there was no significant difference between the two. The estimated and measured values were analyzed using Equations (7) and (8). The results showed that MAE ranged from 0.009 to 0.019 L and NSE ≥ 0.998. This indicated that the overall error between the model estimates and the experimentally measured values was small. The estimated model can accurately predict the I values under each combination of influencing factors.

During irrigation, soil water movement was accompanied by the passage of the wetting front. From the 35 sets of simulation scenarios, single-factor comparison scenarios with different combinations of I, K_s, D, B, and Ψ_m influencing factors were selected to map the changes in W, as shown in Figure 7.

As seen in Figure 7, there was a slight variation in the shape of the soil wetting pattern. The contour lines are approximately ‘ellipsoidal,’ and the soil wetting pattern size in three directions follows the pattern of vertical downwards > horizontal direction > vertical upwards (with point 0 being the starting point for the transport of the wetting front). There were differences in the effects of I, K_s, D, B, and Ψ_m on the soil wetting pattern size. Under the same K_s, D, B, and Ψ_m conditions, W tended to increase in all directions as I increased, with the most significant increase in the vertically downward direction, followed by the horizontal direction, and the smallest in the vertically upward direction (Figure 7(a)). The more significant I was, the more water entered the soil, the further the water was transported around, and the larger the wetting pattern volume. In addition, the gravitational potential promotes downward infiltration and inhibits the upward absorption of soil water. For the same I, D, B, and Ψ_m conditions, W in the vertically upward and horizontal directions gradually decreased with increasing K_s, while W in the vertically downward direction gradually increased (Figure 7(b)). This may be because the K_s of the soil characterized the soil texture type in the simulation scenario. For coarse-textured soil with a high K_s, large soil pores are high, and the gravitational potential is more pronounced during infiltration than the Ψ_m. For fine-textured soils with low K_s, the number of tiny soil pores is high, and the effect of the Ψ_m is more pronounced than the gravitational potential during infiltration. Therefore, soil water emerges under the combined effect of gravitational potential and Ψ_m, and the characteristic pattern is shown in Figure 7(b). For the same I, K_s, B, and Ψ_m conditions, W in all directions tends to decrease as D increases (Figure 7(c)). This phenomenon seems to contradict that the more considerable the D, the greater the water infiltration rate, and the faster the transport of wetting fronts. Yet the more significant the D, the faster the water infiltration, and the shorter the irrigation time under the same irrigation quota conditions, therefore producing the phenomenon of the wetting pattern size decreasing and not increasing. For the same I, K_s, D, and Ψ_m conditions, W in all directions tends to decrease as B increases (Figure 7(d)). This was due to an increase in the value of B, which represents an increase in the soil pressure potential at the infiltrating pipe. The pressure potential promotes rapid infiltration of soil water, resulting in a reduction in irrigation time for the same irrigation quota and a decrease in wetting pattern size instead of an increase in size. For the same I, K_s, D, and B conditions, W in all directions tends to increase slightly with increasing Ψ_m (Figure 7(e)). This may be due to the slight change in Ψ_m during this interval, resulting in a more negligible effect of the Ψ_m gradient. In addition, as Ψ_m increases, less water is required to fill the pores, and the water being reduced continues to diffuse outwards, resulting in a larger wetting pattern volume.

In the above qualitative analysis of the influencing factors, W, I, K_s, D, B, and Ψ_m all influence W_j, and the degree of influence varies. To accurately describe the quantitative relationship between I, K_s, D, B, Ψ_m, and W_j, the empirical constants c_j, d_j, and the value of the coefficient of determination R² were obtained by fitting based on numerical simulation data using Equation (5). This is shown in Table 3.

Based on Equations (12)–(14), the simulated data were fitted again to obtain the empirical constant c_i′ and the coefficient of determination R² values, as shown in Table 3.

The ranges of variation in the values of c₁′, c₂′, and c₃′ were all large and could not be determined by averaging. Further analysis revealed that c_j′ was influenced by K_s, D, B, and Ψ_m; c₁′ increased with K_s and Ψ_m and decreased with D and B, c₂′ increased with Ψ_m and decreased with K_s, D, and B, c₃′ tended to decrease with K_s, D, B, and Ψ_m, and all had good power function relationships between K_s, D, B, and Ψ_m and a′ (Figure 8).

Equations (18)–(20) have a coefficient of determination R² ≥ 0.975, indicating that the power function continuous multiplication model fits with high accuracy and can accurately describe the quantitative relationship between W from gravity subsurface hole irrigation and the influencing factors. Through comparing the standardized regression coefficients (Table 4) of the factors in the same direction, it is known that the influence of the factors in the vertically downward direction was, in descending order, I > K_s >D > B > Ψ_m. The influence of the factors in the horizontal direction was, in descending order, I > D >K_s > B > Ψ_m. The influence of the factors in the vertically upward direction was, in descending order, I > D > K_s > B > Ψ_m. The influence of I on W was the most significant in all three directions, and the influence of Ψ_m was the least significant.

To evaluate the accuracy of Equations (18)–(20), the measured values were plotted against the predicted values using three sets of indoor soil tank test data (Figure 9) and the two indicators of MAE and NSE in three directions were obtained by regression analysis.

As seen in Figure 9, the model estimates of W in the three directions for the three groups of indoor soil tank test data showed the same trend as the measured values of the test. The P-values were calculated using a t-test with P > 0.05, indicating that both were not significantly different and were in agreement. The estimated and measured values were analyzed using Equations (7) and (8). The results show that the MAE in the vertically downward direction ranges from 0.096 to 0.320 cm and NSE from 0.979 to 0.998. MAE in the horizontal direction ranges from 0.011 to 0.060 cm and NSE from 0.998 to 0.999. Vertically upward MAE is between 0.164 and 0.351 cm and NSE between 0.963 and 0.989. The prediction effect of the W estimation model in the three directions was acceptable and can accurately predict W under the combination of various influencing factors.

In this paper, under the background of gravity subsurface hole irrigation, a water-saving irrigation mode, 35 simulation scenarios are set up considering different influencing factors. Using numerical simulation as the primary technical method, the dynamic variation rules of cumulative infiltration and wetting body size under different influencing factors are studied. In this research, 35 simulation scenarios were set up to numerically investigate the dynamics of I and wetting pattern size dimensions under different combinations of influencing factors. I was positively correlated with K_s, D, B, and T and negatively correlated with Ψ_m. The five influencing factors on I were, in descending order, T, K_s, D, B, and Ψ_m. Starting from the edge of the weep hole, W in all directions increases with I and Ψ_m and decreases with an increases in D and B. Under the combined effect of gravity and Ψ_m, W in the vertically downward direction increases as the K_s of the soil increases, while W in the vertically upward and horizontal directions decreases. The five influences on W in the vertically downward direction were, in descending order, I, K_s, D, B, and Ψ_m, while K_s and D were swapped in the vertically upward and horizontal directions. Based on the simulation results, a model for estimating ‘I’ in the form of power function continuous multiplication was developed that incorporated T, K_s, D, B, and Ψ_m. A model for estimating W in three directions (vertically downward, vertically upward, and horizontally) was developed that incorporated I, K_s, D, B, and Ψ_m. The model's reliability was evaluated using three sets of measured data, with the MAE close to 0 and the NSE close to 1. This indicates that the model is highly accurate and that it is feasible to calculate the I and W for gravity subsurface hole irrigation soils.

Yanwei Fan designed the experiments and wrote the paper; Chunyan Zhu conducted experiments and analyzed the data; and Guilin Bai and Tianhua Ma and Zhenchang Wang revised the paper.

This research was supported by the National Natural Science Foundation of China (No. 51969013) and the National Natural Science Foundation of Gansu Province (No. 21JR7RA225).

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

The authors declare there is no conflict.