Numerical investigation of the influence of banded sand ditches on water infiltration in fine-textured soil

K_ss

Sand ditch soil saturated hydraulic conductivity

K_so

Original homogenous fine-textured soil saturated hydraulic conductivity

W

Sand ditch width

S

Sand ditch spacing

D

Sand ditch depth

SDCS

Sand ditch coarse-textured soil

OHFS

Original homogeneous fine-textured soil

γ

Bulk density

i

Soil infiltration rate (cm/min)

I

Cumulative infiltration (mL)

A

Soil area (cm²)

θ

Soil water content (cm³/cm³)

t

Time (min)

K(h)

Unsaturated hydraulic conductivity (cm /min)

h

Soil water pressure head (cm)

x

Horizontal coordinate (cm)

z

Vertical coordinate (cm)

K_s

Saturated hydraulic conductivity (cm/min)

S_e

Effective degree of saturation

a

Empirical parameter

n and m

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

t₀

Turning point of soil infiltration rate curve under banded sand ditch pattern (min)

RMSE

Root mean square error

NSE

Nash Sutcliffe efficiency

Mi and Si

Observed and predicted values, respectively

M_m

Mean of observed data

N

Total number of observations

CLM

Clay loam*

SLM

Silt loam*

LM1

Loam*

LM2

Loam**

LM3

Loam***

LSD

Loamy sand*

SD1

Sand****

SD2

Sand***

δ

Penetration enhancement rate

Global climate change and intense human activities have made global water security issues increasingly prominent. As a result, society is facing increasingly serious water problems such as water shortages, water pollution, floods, droughts, and other related crises (Abbaspour et al. 2015; Pan et al. 2019). With increased urbanization, which is an important symbol of human social development, water safety has aroused widespread concern from all social classes (Wang & Zhang 2019; Keya & Aysha 2021). With the rapid expansion of urban areas, the original landscape of infiltration underlying the surface (woodland and grassland) is replaced by hardened ground. This phenomenon causes a series of ecological and environmental problems such as impeding rainfall infiltration, increasing surface runoff, urban waterlogging, soil drying, and groundwater-level decline. To reduce the impact of urban construction on the environment, it is necessary to fully utilize urban rainwater resources.

Many countries and regions have introduced the construction of a ‘sponge city’ with natural accumulation, penetration, and purification, which encourages the development of low-impact technologies, such as sponge roads, rain gardens, recessed green spaces, constructed wetlands, and green roofs so that rainwater can be absorbed and utilized on the spot. For local rainwater absorption, the role of the infiltration effect of vegetation, soil, and other natural underlying surfaces is important to maximize the conversion of rainwater into soil water. This process mainly depends on rainfall intensity and the soil infiltration rate. In a single rainfall event, it is impossible to artificially control rainfall intensity and improvement of the rainwater local absorption rate can only be achieved by increasing the soil infiltration rate.

Infiltration characteristics are a fixed property of soil that determines the efficiency of rainwater conversion into soil water. The main factor affecting the soil infiltration rate is soil texture. Furthermore, soil bulk density and soil configuration have a significant impact on soil infiltration characteristics. It was previously determined that setting sand ditches in weak permeability soil enhances water infiltration. Abu-Zreig et al. (2000) and Abu-Zreig & Tamimi (2011) experimented with sand ditch rainfall collection and found that sand ditches can effectively collect rainwater, intercept runoff, reduce sediment transport, prolong the peak reach time, and improve crop production for dry landscapes. Widomski et al. (2010) evaluated the infiltration enhancement effect of a banded sand ditch in an orchard in Olszanka and observed that the infiltration rate of terraces equipped with sand-filled drains increased by an average of 12.6%. Abu-Zreig et al. (2020) confirmed that sand ditches increased water infiltration volume and soil moisture storage, especially in high-density soils, and the presence of sand ditches increased the rainwater infiltration volume by 156%. Barber et al. (2003) studied the effects of ecological ditches (composed of compost, sand, gravel, and perforated drainage pipes) on the peak delay time of urban rainstorms and simulated peak reduction. Ecological ditches can reduce the peak flow of rainstorms by 10–50%. The size of the ditch and media type are key factors for the peak delay time.

Due to the high infiltration rate of sand ditches, rainfall, and runoff can quickly enter the soil profile through the sand ditch to limit water evaporation and quickly store precipitation resources in deep soil, which help solve the water demand of plants in arid periods. Sand ditch rainfall harvesting technology provides an important method for alleviating the negative hydrological impact of urbanization. When applied to road green belts, urban road runoff, rain gardens, and green roofs in cities, these ditches increase the soil infiltration rate and cities obtain a higher rainwater adsorption rate. This has a major impact on solving the water shortage problem in the hydropenia season and flood problems in the rainy season. Therefore, to thoroughly implement sand ditch technology in urban rainwater management, reasonably set sand ditch parameters, and avoid unnecessary consumption of human and financial resources, it is necessary to quantitatively analyze the change of the soil infiltration rate after setting sand ditches. Numerical simulation is an effective method for the analysis of a large amount of data. HYDRUS can be used to simulate the infiltration process of water in different layers of soil under different water accumulation conditions (Šimůnek et al. 2018). Wang et al. (2014) simulated and analyzed the infiltration characteristics of fine soil with a coarse sand interlayer under the conditions of pounding rainfall with HYDRUS software, and the grainy sand interlayer with higher matrix suction promoted rather than inhibited soil infiltration. Sasidharan et al. (2020) used HYDRUS software to conduct a numerical simulation and compared the groundwater recharge effects of drywall (DWS) and infiltration bases (IBS) on the penetration of tight clay layers. The authors found that drywall facilitates the rapid infiltration and recharge better than infiltration bases. Hilten et al. (2008) and Xie et al. (2019) evaluated the effect of green roofing on rainwater runoff with the help of the HYDRUS simulation. Saito et al. (2012) utilized experimental results to confirm that HYDRUS can reliably simulate soil water movement in sand ditches and evaluated the water storage capacity of sandy ditches.

Based on the HYDRUS model, this paper choose the soil water movement equation and defined the solution conditions under a banded sand ditch pattern. The effects of soil parameters (sand ditch coarse-textured soil (SDCS) and original homogeneous fine-textured soil (OHFS)) and sand ditch structural parameters (sand ditch spacing (S), sand ditch width (W), sand ditch depth (D)) on the soil infiltration rate were simulated. This simulation data was used to analyze the influence degree of five factors on the soil infiltration rate. Then, the calculation model of the soil infiltration rate was established. The reliability of the calculation model was verified through laboratory experiments.

The laboratory test was used to evaluate the reliability of the soil infiltration rate prediction model under a banded sand ditch pattern. The tested loam and fine sand were from Jingtai County, Baiyin City, China, and the silty loam and coarse sand were from Qilihe District, Lanzhou City. The soil was collected from 0 to 40 cm underground. Soils used for the experiment were air-dried, rolled, uniformly mixed, and sieved through a 2 mm sieve. The bulk density (γ) of the soil used in the experiment was determined by a ring knife method, and the soil saturated hydraulic conductivity was measured using the constant head method.

Figure 1 illustrates the equipment for the laboratory experiments, and Figure 2 shows the device diagram of a laboratory experiment. The experimental device consists of four parts as follows: an adjustable bracket, a soil tank, a water supply system, and a drainage system. The water supply system is composed of a water bucket, inlet pipe, peristaltic pump, and outlet pipe; the drainage system consists of an outflow pipe and a measuring cylinder. Both soil tanks are composed of 10 mm thick plexiglass bonds, and the internal dimensions of the two soil tanks are 30 × 10 × 100 cm and 40 × 10 × 100 cm (length × width × height). Vent holes occur at the bottom of the soil tank at 2 mm diameter to prevent air resistance during infiltration. The diameter of the water bucket was 10 cm and the height was 100 cm. The outflow pipe inlet was set 2 cm above the soil surface to maintain a constant 2 cm head on the soil surface during infiltration.

Considering SDCS, OHFS, S, W and D factors, four groups of laboratory experiments were designed and test schemes are shown in Table 1. First, the air-dried fine-textured soil was loaded into soil tanks in layers according to the soil bulk density. To prevent the sidewall soil from collapsing, the original soil was supported with a thin metal plate and sand was loaded into the sand ditch in layers according to the bulk density. The thickness of each soil layer was 5 cm and a roughening treatment was conducted between layers. Finally, the metal plate was slowly pulled up. The soil box filled with soil samples was deposited for 24 h and then the testing began. Water was added to the predetermined 2 cm constant head at the beginning of the test. Moreover, the peristaltic pump was switched to supply water according to the set flow rates (q) and the excess water flowed out into the measuring cylinder. During experiments, cylinder data were recorded at the predetermined time until infiltration was stable. Table 1 shows design schemes of the laboratory experiments.

where θ_r is the residual water content (cm³/cm³), θ_s is the saturated water content (cm³/cm³), K_s is the saturated hydraulic conductivity (cm/min), S_e is the effective degree of saturation, α is an empirical parameter (1/cm), and m and n are empirical constants (–).

Figure 3 shows the structure of the soil profile in the banded sand ditch pattern. Sand ditches of the same width and depth are arranged in parallel with the same interval. Therefore, this is an axisymmetric two-dimensional plane problem. Considering the symmetry of the sand ditch, half of the control area of the sand ditch was selected as the simulation area (Figure 4). The initial matrix of test soils are used as the initial conditions. In most cases, rainwater cannot infiltrate and there will be a thin layer of water on the upper boundary of the soil where the constant head is set (h = 2 cm). Free drainage BC was used along the bottom boundary without the influence of groundwater and rainfall. The contact boundary between sand and the original soil is not affected by the outside and so it is set as a free drainage BC and a no-flux boundary for all remaining boundaries. When solving, the Galerkin finite element method was used to disperse the soil profile in space, and the implicit difference scheme was used to disperse the profile in time; the simulated soil depth was 100 cm, the time step was 0.1 min, the space step was 1 cm, and the total simulation time was 120 min.

To establish a soil infiltration rate prediction model of banded sand ditches that is more applicable and truly reflects the infiltration process, 10 soil parameters from the published literature were used for numerical simulation and five of these soils were selected from waterlogged cities (Carsel & Parrish 1988; Young et al. 2002; Ma et al. 2011; Fan et al. 2016; Li et al. 2018). The parameters of the 10 soils (VG-M) are listed in Table 2.

To set reasonable simulation schemes, the appropriate size range of the sand ditch is selected according to the existing data (Abu-Zreig & Tamimi 2011; Abu-Zreig et al. 2020; Saito et al. 2012). Combinations of 10 types of soil (Table 3), 5 types of spacing (S = 40, 50, 60, 70, and 80 cm), 5 widths (W = 6, 8, 10, 12, and 14 cm), and 5 depths (D = 10, 15, 20, 25, and 30 cm) were used to form 51 groups of simulation schemes. These designs included 5 groups of original soil infiltration contrast schemes without a sand ditch, 21 groups of single factor simulation schemes, which can be used to analyze single influencing factors, and 25 representative orthogonal simulation schemes. A total of 46 simulation scheme groups were used to explain the change of the soil infiltration rate under a banded sand ditch mode and these are listed in Table 3.

To illustrate the influence of single factors changing the original soil infiltration rate under the banded sand ditch pattern, the numerical simulation results of 21 groups of single factor schemes were used for analysis.

Under conditions of S = 50 cm, W = 10 cm, and D = 20 cm, the dynamic change process of the soil infiltration rate under different combinations of SDCS and OHFS is simulated and shown in Figure 5. To compare and analyze the infiltration characteristics of the banded sand ditch and the original homogeneous fine-textured soil, the original control soils without sand ditch infiltration rate curves were added to Figure 5(a) and 5(b).

Figure 5 shows that there is an obvious turning point in the soil infiltration rate curve in the banded sand ditch pattern that we defined as ‘t₀.’ The values of ‘t₀’ change with changes in SDCS and OHFS. During the infiltration process, the initial value of ‘i’ was relatively large. With increased time, ‘i’ rapidly decreases, and when the infiltration progresses to ‘t₀,’ ‘i’ will rapidly decrease again and then gradually stabilize.

Figure 5(a) shows that when only changing the texture of OHFS, the more delicate the OHFS, the higher the penetration enhancement rate (δ) (defined as the percentage increase in ‘i’ in the banded sand ditch mode compared with ‘i’ in the same period without the sand ditch mode). For example, at the end of the simulation time (2 h) and when the original homogeneous soil CLM sets the sand ditch soil conditions of SD2 as S = 50 cm, W = 10 cm, and D = 20 cm, ‘δ’ is 102%; in contrast, LM3 has relatively lower clay content under the same conditions and ‘δ’ is only 19.8%. This is because when the water infiltration starts, water simultaneously moves into the SDCS and OHFS, but because the saturated hydraulic conductivity (K_ss) of SDCS is large, the SDCS will soon saturate. While the saturated hydraulic conductivity (K_so) of OHFS is small, the wet front moves slow and in a short time it is unable to reach the bottom of the sand ditch. Thus, a lateral seepage surface appears at the interface between the sand ditch and the original soil (the size of the seepage surface is the sum of the width of the sand ditch and the depth of the sand ditch minus the OHFS wet front migration distance) and causes the SDCS water to move through the lateral seepage surface in the OHFS. The finer the texture of OHFS, the smaller the K_so and the slower the migration speed of the wet front, the larger the lateral seepage surface, the greater the lateral seepage volume, and the more obvious the effect of improving permeability. In addition to the above phenomena, the effects of OHFS on ‘t₀’ are all manifested when there is a greater clay content in OHFS. In addition, the smaller ‘t₀’ is because before the water reaches the bottom of the sand ditch and while the water infiltrates rapidly in the vertical direction in the SDCS, it also infiltrates laterally to the OHFS. The higher the soil clay content, the weaker the water conductivity and the less water can penetrate the OHFS from the SDCS side, and the SDCS becomes easy to saturate.

Figure 5(b) shows that when OHFS, S, W, and D are fixed, the higher the sand content of SDCS and the smaller ‘t₀.’ For example, in LM1 with the SD2 ditch and LM1 with SD4, ‘t₀’ is 15 and 3 min, respectively, because the higher the sand content in the SDCS, then the larger K_ss, the faster the water migration in the SDCS, and the shorter the time required for water to reach the bottom of the sandy ditch. The influence of SDCS on ‘i’ mainly occurs before ‘t₀.’ The higher the sand content of SDCS, the greater the initial value of ‘i.’ For example, in the first minute of infiltration (for LM1 with the LSD ditch), the value of ‘i’ is only 0.294 cm/min, while for the case of LM1 with the SD4 ditch, ‘i’ increases to 0.764 cm/min. After infiltration reaches ‘t₀,’ then ‘i’ rapidly decreases. The higher the sand content of SDCS, the faster the decay of ‘i.’ When time reaches 70 minutes, the influence of the SDCS texture type on ‘i’ is weakened and its value stabilizes because after the water reaches the bottom of the sand ditch and it becomes saturated, then it cannot increase infiltration by storing excess water and can only increase the infiltration of the OHFS through the lateral seepage surface, which results in a rapid decrease in ‘δ.’ As time passes, the OHFS wetting front reaches the bottom of the sand ditch and the lateral seepage surface disappears. Infiltration in the banded sand ditch pattern is transformed into a vertical one-dimensional infiltration in OHFS, and the texture type of SDCS has no effect on ‘i.’

Figure 6 shows the variation in the soil infiltration rate under the banded sand ditch pattern, which is under the constant soil conditions of LM1 with the SD2 ditch by changing the sand ditch structure (the size of parameters S, W, and D).

Figure 6 shows that like the influence of SDCS and OHFS on ‘I,’ ‘t₀’ also exists under the influence of sand ditch structure parameters (S, W, and D). Using ‘t₀’ as the boundary, ‘i’ has a rapid decline in two time periods and gradually stabilizes. Figure 6(a) shows that S does not affect ‘t₀.’ In the case of LM1 with the SD2 ditch (W = 10 cm and D = 20 cm), the ‘t₀’ values corresponding to different ‘S’ values (40, 50, 60, 70, and 80 cm) are all 15 min. Furthermore, as ‘S’ increases both ‘i’ and ‘δ’ decrease. The reason is that under certain W conditions, the increase in S indirectly reduces the proportion of sand ditches in the soil profile (W/S) and thereby weakens the extreme permeability of the sand ditch. In addition, the wider the ‘W,’ the weaker the influence of lateral infiltration of water in the sand ditch on vertical infiltration and the shorter the time for water to reach the bottom of the sand ditch. As ‘D’ increases, ‘t₀’ also gradually increases while ‘i’ and ‘δ’ also significantly increase (Figure 6(c)). The reason is that as ‘D’ becomes deeper, then it takes for the water to reach the bottom of the sandy ditch. In addition, as ‘D’ becomes deeper, the larger the lateral contact surface between the sand ditch and the original soil, the greater the amount of lateral infiltration of water in the sand ditch, and the more noticeable the effect of increasing infiltration.

The saturated hydraulic conductivity of the soil is affected by the soil texture and plays an important role in the engineering design of irrigation and drainage systems. To simplify the construction of the model, K_ss and K_so are utilized to characterize the coarse-textured soil and the original homogeneous fine-textured soil difference. The sensitivity analysis of the five factors that may affect ‘t₀’ and the results show that the values of ‘P’ for K_ss, K_so, and D are less than 0.01, which indicates that these factors have a significant impact on ‘t₀.’ The effect of ‘S’ is small (P = 0.13) and can be ignored. The influence degree of each factor is ‘K_ss > D>K_so > W.’ To quantitatively analyze the functional relationship between ‘t₀’ and ‘K_ss, D, K_so, and W,’ except for the factors to be analyzed, other factors are controlled so that they remain unchanged and a trend chart of ‘t₀’ is created under the influence of this factor (Figure 7).

Using numerical simulation data to evaluate the simplified prediction model that characterizes ‘t₀,’ we compared and analyzed the calculated value of the simplified prediction model with the numerical simulation value. The comparison result is shown in Figure 8.

Figure 8 shows that points of the ‘t₀’ simulated value and the model-calculated value coordinates are distributed and are near the 1:1 line, which indicates that the correlation is good and is used to calculate the error-index parameters. The RMSE is 2.07 min and the NSE value is 0.962, which indicates that the ‘t₀’ prediction model has good accuracy. The significance of the model is to estimate the time when water reaches the bottom of the sand ditch to build a more accurate prediction model of ‘i.’

Using a significance analysis of each influencing factor, it was concluded that during the t ≤ t₀ stage, each factor's influence degree of ‘i’ can be described as ‘K_ss>W+S>K_so>D’. Furthermore, ‘K_ss’ and ‘W’ have a significant impact on ‘i.’ Similarly, during the t > t₀ stage, each factor's influence degree for ‘i’ can be described as ‘K_so > D + S > K_ss > W,’ and ‘K_so’ and ‘D’ have a significant impact on ‘i’. Using Equation (10) and considering ‘t₀’ as the demarcation, 46 sets of simulation plans (Table 3) are applied to the recursion data fitting to calculate the infiltration rate versus time. Then, values of β₁, β₂, ε₁, ε₂, and R² are obtained and results are shown in Table 4.

Table 4 shows that the coefficient of determination of the regression line fitted by all simulation schemes is close to 1 (R² ≥ 0.96). These results indicate a good quantitative relationship between the infiltration rate and various influencing factors, and the Kostiakov infiltration formula can be used to quantitatively describe the change in the soil infiltration rate with time under a banded sand ditch pattern. The ε₁ and ε₂ exponent values show a small fluctuation range around their average value. Regarding results of the parameter sensitivity analysis, slight fluctuations in the values of ‘ε₁’ and ‘ε₂’ have little effect on the infiltration rate. To simplify this calculation, the average values of fitting parameters ε₁ and ε₂ can be used (0.33 and 0.65. respectively) to fix them. However, values of β₁ and β₂ greatly fluctuate. In this case, the average value cannot simply be used as a fixed value and further analysis is needed.

Substituting ε₁ = 0.34 and ε₂ = 0.63 into Equation (10), using ‘t₀’ as the demarcation point, and again fitting the numerical results of 46 sets of simulation schemes (Table 4), we can obtain new values of ‘β₁’ and ‘β₂.’ Next, quantitative analyses of the functional relationship between β₁, β_2, and K_so, K_ss, S, W, and D are shown in Figure 9.

Figure 9 shows that throughout the infiltration stage, K_ss, K_so, W, D, and S have a good linear relationship with β₁ and β₂.

To verify the reliability of the empirical model, Figure 10 shows that the calculated value of the estimated model is compared and analyzed with a measured value of the indoor test. The error analysis results are listed in Table 5.

Figure 10 shows that the calculated value agrees with the measured value. With Formula (6) and Formula (7), the estimation model is statistically analyzed and RMSE values are between 0.0085 and 0.0290 cm/min, which are close to 0, and the NSE value is between 0.824 and 0.989, which is close to 1. These results suggest that the model has a good prediction effect and can be used to estimate the soil infiltration rate under the banded sand ditches model. Moreover, the error mainly occurs in the advanced stage of infiltration (approximately at the stable infiltration stage), which is caused by the error of the Kostiakov formula used in this study. In addition, the Kostiakov formula is an empirical model with a simple form, which causes fitting results to be incompatible with measured values at the beginning and end. Specifically, the fitted infiltration rate is often too large at the initial moment of infiltration, and the fitted infiltration rate is lower than the measured value at the stage of stable infiltration.

Based on the HYDRUS model, a mathematical model of soil water movement under a banded sand ditch pattern was established, and 46 scenarios were set to numerically investigate the effects of soil texture parameters and sand ditch structure parameters on the soil infiltration rate under this pattern. When compared with the original homogeneous fine-textured soil infiltration process, the following conclusions were drawn from this information. The banded sand ditch shows a significant increase in the infiltration effect and there is a turning point ‘t₀’ in the curve of the soil infiltration rate with time. In addition, ‘t₀’ is not related to S but it increases with increasing K_so and D and decreases with increasing K_ss and W. Therefore, there is a determined power function relationship (R² = 0.958) between ‘t₀’ and K_ss, K_so, D and W. The Kostiakov formula can be used to describe the infiltration process of banded sand ditch furrow soil. Because the Kostiakov formula has a small index change, the average values are 0.34 (t ≤ t₀ stage) and 0.63 (t > t₀ stage) when determining the formula parameters. Coefficients of the Kostiakov formula and K_ss, K_so, W, S, and D satisfy the linear relationship. A model for estimating the soil infiltration rate under a banded sand ditch pattern was proposed, and the reliability of the model was verified by laboratory experiments (RMSE values are between 0.0085 and 0.0290 cm/min and the NSE value is between 0.827 and 0.989), which can be used to predict the soil infiltration rate under a banded sand ditch pattern. In this research, only five extremely important influencing factors (K_ss, K_so, W, S, and D) were analyzed when quantitatively describing the permeation effect under a sand ditch model and effects of rainfall intensity, soil bulk density, and soil initial moisture content on the soil infiltration rate under a sand ditch mode were not considered. Therefore, in future research, other influencing factors should be examined, which would improve the model established in this study and result in a more accurate description of the change process of the soil infiltration rate under a strip sand trench mode.

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.