Effects of Subsoiling with Different Wing Mounting Heights on Soil Water Infiltration Using HYDRUS-2D Simulations

Abstract

Subsoiling is an essential practice in conservation tillage technology. The amount of disturbed soil at various depths resulting from subsoilers with different parameters has an important effect on soil properties (e.g., bulk density and water infiltration). The information regarding the effects of subsoiling on the characteristics of soil water infiltration is essential for the design of subsoiling tools. In this study, the effects of the wing mounting height (h) (75–155 mm) of the subsoiler on soil disturbance and soil water infiltration were modelled using HYDRUS-2D and validated using field experiments. Results showed that reducing h values resulted in larger soil disturbance area ratios, soil water infiltration rates (f(t)), distances of vertical wetting front movement (DVWs), accumulative infiltrations (AINs), and soil moisture contents at depths of 10–30 cm. The relationships among characteristics of soil water infiltration, h and time (t), were developed. The stable infiltration rates (f_s) varied quadratically with h and the corresponding coefficient of determination (R²) was 0.9869. The Horton model is more suitable for describing the relationship between f(t) and t under the tested soil conditions, as compared with the Kostiakov and Philip models. According to the results of soil water content at different depths from the HYDRUS simulations and field experiments, the developed soil water infiltration model had a good accuracy, as indicated by RMSEs of <0.05, R² values of >0.95, and mean relative errors of <12%. The Above results indicated that increasing the hardpan disturbance by optimizing wing parameters of the subsoiler could improve soil water infiltration characteristics.

1. Introduction

Mechanical subsoiling is one of the key practices in conservation tillage technology. It is mainly employed to break hardpans, eliminate the issue of soil compaction, improve the characteristics of soil water infiltration, and restore the productivity of farmland [1,2,3,4,5]. The wing mounting height of a winged subsoiling tool can greatly affect subsequent soil water infiltration after subsoiling by changing the disturbance behaviour of soils at various layers. The existing research about winged subsoilers has mainly investigated the effects of geometrical or working parameters of subsoiler’s wings on resistances and soil’s dynamic attributes [6,7,8,9]. Spoor and Godwin [6] investigated the effect of wing length on tillage resistance and soil disturbance area. Raper [10] compared the tillage performance of several column-type subsoilers from Kelley Manufacturing Co. (Tifton, GA, USA) and found that adding different wings to the shank more or less increased the specific draft force of subsoilers. A study from Li et al. [7] showed that both draught force and soil disturbance area increased with the addition of two wings to a column subsoiler. Wang et al. [8] investigated soil-winged subsoiler interactions using discrete element method simulations and found that the draft force of the subsoiler in the hardpan was the main source of total draught force in the winged subsoiler. Xia [9] found that wing mounting angle and mounting height significantly affected the shape of the furrow profile disturbed by the winged subsoiler. Reducing top layer soil disturbance and increasing hardpan disturbance during subsoiling operations are often identified as methods to conserve soil and preserve soil water [10,11,12,13,14]. However, differences in the density of the top layer at different locations can be reduced by completely disturbing the top layer soil, which helps to prepare a good seedbed. Moreover, increasing disturbance in the hardpan layer can cause additional energy consumption because of re-compaction from future field operations of agricultural equipment [5]. For a given soil condition, the information regarding the appropriate soil disturbance at various layers that leads to the best water infiltration, which is the premise of designing high-performance subsoilers, is lacking in the literature.

The infiltration process has long been a focus of hydrology and agricultural engineering since it provides the water that is available for plants and groundwater recharge and defines water runoff at soil surface [15,16,17,18]. Many mathematical models have been developed to evaluate the computation of infiltration and they can be classified into physically based models (PBMs), semi-empirical, and empirical models [19,20]. Compared with semi-empirical and empirical models, the physically based models can substantially describe the detailed infiltration process. The Richards equation was derived by the mass conservation law and Darcy’s law and is one of the most commonly used PBMs [19,20]. However, the Richards equation is strongly non-linear and cannot be solved analytically, especially under complex initial and boundary conditions [19,21]. With the development of computer technology, numerical simulations have become effective tools for solving the Richards equation. Based on the finite element method (FEM), the HYDRUS code was developed to solve the Richards equation and has been widely used to simulate water movement in variably saturated media [22,23]. Ebrahimian et al. [22] predicated soil water content, nitrate concentrations, and deep percolation caused by the geometry of the infiltration domain in furrow irrigation using the HYDRUS-2D model. Karandish and Simunek [24] investigated the influence of various water-saving irrigation strategies on maize water footprints using the HYDRUS model. Argyrokastritis et al. [21] investigated the ponded infiltration processes that occurs in agricultural lands irrigated by flooding of their soil surface or under insufficient drainage conditions using the HUDRUS-1D code. Ma et al. [19] modelled water infiltration in a large layered soil column with a modified Green-Ampt model and HYDRUS-1D. The above relevant works mainly focused on the effects of the geometry of the infiltration domain and initial water conditions on soil water infiltration characteristics. The effects of the amount of soil disturbance in different layers due to the variation of soil-engaging tools on water infiltration are absent in the literature.

In this study, soil disturbance area ratio (SDAR) was used to quantify the relative soil disturbance between the hardpan and the top layer. It was defined as the ratio of the disturbance area of the hardpan to that of the top layer. Therefore, the objectives of the study were to: (1) develop a soil water infiltration model using HYDRUS-2D based on the soil properties at various depths, (2) investigate the effects of variation of wing mounting heights of the subsoiler in subsoiling process on soil disturbance amount at various layers and soil water infiltration characteristics, and (3) validate the model using data from double-ring water infiltration tests in the field.

2. Materials and Methods

2.1. Study Site Description

The soil used is a Lou soil (1% clay, 9% silt, 74% sand, and 16% gravel,), which has a granular structure with secondary loess as the parent material and clay as the loamy material [12,13,14]. In this region, summer maize and winter wheat are two common crops. After harvesting the crops, conventional tillage was generally performed before this field test.

Soil cone index (CI) in the tested field was determined using a SC-900 digital display soil compactness meter (Figure 1A). The soil compaction degree and formation of hardpan soil (i.e., Plough sole) can be evaluated using the soil cone index [4,25]. Previous researchers [4,26] have shown that increasing the soil cone index can decrease the elongation and growth of crop roots. The threshold level at which soil strength hinders root elongation varies with plant species, but is usually 1800 kPa or larger [26,27]. The soil cone index in the depth of 16–30 cm from the surface was larger than 1800 kPa (Figure 1B). Therefore, the hardpan layer at the depth of 16–30 cm could limit the development of crop roots, especially for drier soil, as the strength of soil rises as soil dries [27,28,29].

2.2. Experimental Design and Data Collection

2.2.1. Field Experiment

Subsoiling process was performed in the field using subsoilers with wing mounting heights ranging from 75 to 155 mm and a non-winged subsoiler (Figure 2A,B). The working speed and working depth were determined in accordance with the working parameters during subsoiling in local region which were 3 km h⁻¹ and 30 cm, respectively [14,30]. After subsoiling operations, the profiles of disturbed soil under different subsoiling treatments were determined by a profile meter, which is made up of 150 wooden bars with a width of 1 cm. The bars in the profile meter can freely adjust their vertical positions under their own gravity force. The specific measuring method and operations of soil disturbance profiles can refer to the studies from Hang et al. [14] and Chen et al. [29]. Initially, furrow profiles of the soil were traced on engineering graphic paper (grid spacing: 1 mm); the soil disturbance areas of top layer (0–16 cm from surface) and hardpan (16–30 cm from surface) (Figure 2C) were determined by the grid number in the furrow profile and the area of each grid (1 mm²), respectively. Soil disturbance area ratio (SDAR) was used to evaluate the relative disturbance area between different soil layers and it was calculated as follows.

$$\mathrm{SDAR}=\frac{A_h}{A_t}$$

(1)

where A_h and A_t stand for disturbance areas of hardpan and top layer, respectively.

Soil water infiltration test was carried out using a double-ring infiltration meter (DIM) (Figure 2D). The specific operations followed these steps: (1) the DIM was placed at the middle of subsoiler path and vertically driven into the soil using a rubber hammer (~10 cm depth); (2) the gap between DIM and the neighbouring soil, created during the hitting process, was filled by soil and the disturbed soil out of the DIM was compacted; (3) the 5 cm high locations from the surface were marked on both inner and outer rings of DIM [23,31,32]; (4) water was quickly injected into the inner and outer rings at the same time, until the water level reached their marked positions; (5) water was continuously injected using a Marriotte bottle for more than 80 min, until the water infiltration rate is basically stable in unit time (i.e., 5 min).

The cutting-ring method was used to collect soil samples before and after subsoiling as well as after water infiltration from seven depths (i.e., 0–10, 10–20, 20–30, 30–40, 40–60, 60–80, 80–100 cm). Soil samples from three locations were weighed wet and oven-dried at 105 °C to a steady weight. The soil bulk density and moisture content (Equation (2)) were determined according to the wet and dry soil weight. The mean of three repetitions was used for both the soil moisture content and bulk density. The soil moisture content at various depths after water infiltration were used for the validation of soil water infiltration model.

$$\mathrm{Soil} \mathrm{moisture} \mathrm{content}=\frac{\mathrm{wet} \mathrm{soil} \mathrm{weight} - \mathrm{dried} \mathrm{soil} \mathrm{weight}}{\mathrm{dried} \mathrm{soil} \mathrm{weight}}$$

(2)

2.2.2. FEM Simulations and Validation

Model Development

HYDRUS-2D software (Version: 2.04) was used for developing a soil water infiltration model. The Galerkin method was used to solve the soil disturbance profiles created by various subsoiling treatments as suggested by Fu et al. [33]. The model development consists of initial and boundary conditions, mesh generation, and soil hydraulic characteristic parameters.

• Initial conditions: The model initial conditions consist of the initial soil moisture content, bulk density, and particle size distribution, which were set as the measured data from above field tests (Table 1).

Table 1. Soil moisture content, density and texture before subsoiling and soil density after subsoiling.

Soil Conditions before Subsoiling						Soil Density after Subsoiling (g cm−³)
Depth (cm)	Moisture Content (%)	Density (g cm−³)	Sand (%)	Silt (%)	Clay (%)	Treatment	Depth (cm)
							0–10	10–20	20–30
0–10	24.2	1.356	87.10	12.10	0.80	h0	1.326	1.355	1.402
10–20	21.5	1.454	89.23	8.05	2.72	h75	1.307	1.264	1.214
20–30	21.6	1.482	87.93	10.03	2.04	h95	1.310	1.284	1.254
30–40	20.4	1.490	90.40	8.15	1.45	h115	1.300	1.304	1.346
40–60	19.8	1.511	92.06	7.06	0.88	h135	1.317	1.291	1.399
60–80	19.5	1.491	91.32	7.50	1.18	h155	1.318	1.287	1.391
80–100	20.7	1.534	87.01	10.45	2.54

Table 1. Soil moisture content, density and texture before subsoiling and soil density after subsoiling.

Soil Conditions before Subsoiling						Soil Density after Subsoiling (g cm−³)
Depth (cm)	Moisture Content (%)	Density (g cm−³)	Sand (%)	Silt (%)	Clay (%)	Treatment	Depth (cm)
							0–10	10–20	20–30
0–10	24.2	1.356	87.10	12.10	0.80	h0	1.326	1.355	1.402
10–20	21.5	1.454	89.23	8.05	2.72	h75	1.307	1.264	1.214
20–30	21.6	1.482	87.93	10.03	2.04	h95	1.310	1.284	1.254
30–40	20.4	1.490	90.40	8.15	1.45	h115	1.300	1.304	1.346
40–60	19.8	1.511	92.06	7.06	0.88	h135	1.317	1.291	1.399
60–80	19.5	1.491	91.32	7.50	1.18	h155	1.318	1.287	1.391
80–100	20.7	1.534	87.01	10.45	2.54

2.

Boundary conditions: A 60 cm wide (i.e., outer diameter of DIM) disturbance area was set as the constant pressure head (5 cm) at the centre of the soil ridge in accordance with the water infiltration test in the field. Atmospheric boundary condition was used for other locations at the soil surface (Figure 3). Both lateral boundaries were considered to be zero flux faces. The bottom boundary was set as a free drainage boundary as the groundwater depth was relatively large (30–50 m) [34] and the groundwater which moved into the test zone was neglected.

3.

Mesh generation: The FE-Mesh module of HYDRUS was used to generate mesh with size of 50 mm (Figure 3C). To improve the accuracy of the model, the triangular mesh size of the disturbed region was appropriately encrypted (i.e., 5 mm). The initial, maximum, and minimum time steps were set as 0.001 min, 10 min, and 0.0001 min, respectively.

4.

Soil hydraulic characteristic parameters: These parameters were determined using the “Rosetta” module of HYDRUS based on the measured soil bulk density and particle sizes at various depths, including residual moisture content, saturated moisture content, reciprocal of air inlet, shape parameter, and saturated hydraulic conductivity.

Figure 3. Simulations of soil water infiltration test: (A) diagram showing the disturbance area; (B) formation of disturbance profile; (C) mesh generation; (D) distance of vertical wetting front movement (DVW).

Figure 3. Simulations of soil water infiltration test: (A) diagram showing the disturbance area; (B) formation of disturbance profile; (C) mesh generation; (D) distance of vertical wetting front movement (DVW).

Model Application

The developed model was used to investigate the effect of wing mounting height of subsoiler on soil water infiltration characteristics (i.e., water infiltration rate, accumulative infiltration, distance of vertical wetting front movement) and moisture content of soil at various depths after water infiltration.

Developing appropriate relationships between soil water characteristics and time are conducive to better understanding hydrologic process under various conditions. Current models which describe soil water infiltration rate mainly consist of empirical models, semiempirical models, and the models based on physical significance; moreover, Kostiakov model (Equation (3)), Philip model (Equation (4)) and Horton model (Equation (5)) are the representative models of above three types of models, respectively.

f(t) = ϛt^−ϵ

(3)

f(t) = 0.5υt^−0.5 + f_s

(4)

f(t) = f_s + Δfe^−Γt

(5)

where f(t) is water infiltration rate in real time, mm min⁻¹; f_s is stable infiltration rate, mm min⁻¹; Δf is the difference between initial and stable infiltration rate, mm min⁻¹; t was infiltration time, min; ϛ and ϵ are empirical constants; υ is water absorption rate, mm min⁻¹; Γ is model parameter.

Based on the soil water infiltration results of HYDRUS, three models were used to compare and evaluate the water infiltration process, and then a water infiltration rate model, suitable for the tested soil conditions, was obtained.

Model Validation

To validate the accuracy of the HYDRUS model, soil moisture contents at various depths from both the model and the field experiment were compared and analysed using root mean square error (RMSE) and coefficient of determination (R²). The RMSE and R² are generally calculated as follows [33,35,36].

$$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(M{\left(o\right)}_i-M{\left(s\right)}_i\right)}^2}$$

(6)

where M(o)_i and M(s)_i (I = 1…n) stand for simulated and measured soil moisture content, respectively (%).

$${\mathrm{R}}^2=1-\frac{\sum {\left({M\left(o\right)}_i-{M\left(s\right)}_i\right)}^2}{\sum {\left({M\left(o\right)}_i-M\left(o\right)\right)}^2}$$

(7)

where M(o) is the mean of the simulated soil moisture content M(o)_i (I = 1…n) (%).

3. Results and Discussion

3.1. Effect of Wing Mounting Height (h) on Soil Disturbance Area Ratio

The soil disturbance area ratios (SDARs) were 0.4463, 0.4271, 0.3694, 0.3572, 0.3488, and 0.3215 for subsoilers with an h of 75, 95, 115, 135, 155 mm, and a non-winged subsoiler, respectively. Increasing h gave smaller SDARs which corresponded to lower quality of hardpan disturbance. In addition, the soil disturbance area ratio from the non-winged subsoiler was smaller than those obtained from the winged subsoilers. For the non-winged subsoiler, the tillage depth/tool width ratio was 7.5 due to the much smaller tool width of 40 mm, and only the soil near the soil surface was disturbed in a crescent manner. This led to much less soil disturbance in the hardpan, which can explain the much smaller values of SDAR for the non-winged subsoiler.

3.2. Effect of Wing Mounting Height on Soil Water Infiltration Rate

As shown in Figure 4A, with the increase in time, soil water infiltration rates in real time (f(t)), at various wing mounting heights (h), initially reduced rapidly and then gradually reduced, and eventually stable infiltration rates (f_s) were obtained. The f_s varied quadratically with h and the corresponding coefficient of determination (R²) was 0.9869 (Figure 4B). The high value of R² indicates that the quadratic function can be used to describe the relationships between f_s and h. Moreover, decreasing h from 155 to 75 mm gave larger values of soil disturbance area, soil disturbance area ratio, f(t), and f_s, which implied that reducing the wing mounting height of a subsoiling tool would improve the water infiltration rate. Additionally, both the f(t) and f_s of the non-winged subsoiler (i.e., h0) were much smaller than those of a winged subsoiler with various h values. The above results agreed well with the study from Yao [23] who found that larger soil disturbance areas resulted in higher soil water infiltration rates.

To describe the relationship between f(t) and t (i.e., time), values of f(t) at various subsoiling treatments (h0–h155) were firstly fitted using three typical models (i.e., Kostiakov, Philip, and Horton) based on MATLAB R2014b software. As shown in

Table 2. Results of fitting of soil water infiltration rate with time under different conditions.

Treatment	Kostiakov			Philip			Horton
	ϛ	ϵ	R2	υ/ (mm/min)	fs/ (mm/min)	R2	fs/ (mm/min)	Δf/ (mm/min)	Γ	R2
h0	5.138	0.123	0.834	6.356	2.725	0.962	3.277	3.739	0.345	0.979
h75	6.723	0.116	0.831	7.824	3.75	0.963	4.470	4.331	0.326	0.968
h95	6.793	0.129	0.824	8.542	3.524	0.955	4.328	5.200	0.380	0.972
h115	6.660	0.123	0.827	8.066	3.566	0.959	4.314	4.808	0.366	0.972
h135	6.646	0.135	0.834	8.606	3.311	0.958	4.111	5.267	0.383	0.980
h155	6.035	0.138	0.815	8.130	2.936	0.947	3.708	5.220	0.407	0.980

, the variation in the range of R² for fitted equations were 0.815–0.834, 0.947–0.963, and 0.968–0.980, respectively, for the Kostiakov, Philip, and Horton models. The higher values of variation in the range of 0.968–0.980 indicated that the Horton model is more suitable to describe the relationship between f(t) and t under the tested soil conditions.

The stable water infiltration rate (f_s), the difference between initial and stable water infiltration rate (Δf), and the model parameter (Γ) in the Horton model were fitted to the h, respectively, as shown in

Table 3. Fitting results of Horton model parameters.

Fitted Equations	R2
fs(h) = −0.0002h2 + 0.0264h + 3.2815	0.987
Δf(h) = 0.0103h + 3.7757	0.839
Γ(h) = 0.5947sin(0.009957h + 0.306) + 0.3819sin(0.2085h + 8.976)	1.000

. The R² values of the fitted equations were >0.83, indicating a high fitting degree. The fitted equations of h in Table 3 were then substituted into the Horton model, and the relationship among soil water infiltration rate (f(t,h)), wing mounting height (h), and time (t) can be obtained, as shown in Equation (8).

f(t,h) = −0.0002h² + 0.0264h + 3.2815 +
(0.0103h + 3.7757)e^{−[0.5947sin(0.009957h + 0.306) + 0.3819sin(0.2085h + 8.976)]t}

(8)

3.3. Effect of Wing Mounting Height on Distance of Vertical Wetting Front Movement

As shown in Figure 5, during the first 20 min of soil water infiltration, the distance of vertical wetting front movement (DVW) showed no obvious variation under different h values. With the further increase in time (≥25 min), the DVW showed a decreasing trend with the increase in h. The DVW corresponding to the subsoiler without wings (h0) was significantly lower than that of other subsoiling treatments.

A study from Yao [23] demonstrated that DVW and time of water infiltration (t) are approximately satisfied as follows.

DVW = b₁ + b₂t^0.5 + b₃t

(9)

where b₁, b₂, b₃ are all constants.

This function (i.e., Equation (9)) was used to fit the DVW under different h values, and the fitting results under each subsoiling treatment (h) were obtained, as shown in

Table 4. Fitting parameters of DVW with time (t).

Treatment	b1	b2	b3	R2
h0	121.37	−4.35	13.42	1.000
h75	−629.50	268.57	−7.53	0.998
h95	−702.22	288.45	−9.22	0.996
h115	−744.81	299.94	−10.06	0.995
h135	−618.04	264.51	−7.89	0.995
h155	−502.89	213.77	−3.70	0.994

. The R² values of the fitting results were all greater than 0.994, indicating that the fitting results are reliable. According to the results in Table 4, h was fitted for b₁, b₂, and b₃, respectively, and the relationships between h and b₁, b₂, and b₃ were established, respectively, as shown in

Table 5. Coefficients of variation in the fitted equations of DVW with h.

Fitted Equations	R2
b1(h) = 0.0781h2 − 16.178h + 124.84	0.995
b2(h) = −0.0284h2 + 5.8292h − 5.2937	0.998
b3(h) = 0.0021h2 − 0.4463h + 13.512	0.999

. Similarly, the R² values were all greater than 0.995 and the fitting was valid.

By substituting b₁(h), b₂(h), and b₃(h) into the DVW expression, the relationship among DVW, wing mounting height (h), and time (t) can be established as follows.

DVW = (0.0021h² + 0.4463h + 13.512)t + (−0.0284h² + 5.8292h − 5.2937)t^0.5 +
0.0781h² − 16.178h + 124.84

(10)

3.4. Effect of Wing Mounting Height on Accumulative Infiltration

As shown in Figure 6, with the increase in time (i.e., t), the values of accumulative infiltration (AIN) gradually increased for all subsoiling treatments. However, the increase rate of AIN was different to some extent, i.e., the increase rate of AIN from h of >135 mm was much smaller than these from smaller h values. The AIN for h of >135 mm and h0 after 60 min was significantly lower than those of other treatments.

The AIN and infiltration time (t) can be approximated using the following function.

AIN = a₁t² + a₂t + a₃

(11)

Table 6. Fitting equations of accumulative infiltration (AIN) with time (t).

Treatment	Fitted Equations	R2
h0	AIN = −0.0067t2 + 3.9004t + 6.0851	0.999
h75	AIN = −0.0091t2 + 5.2844t + 8.3945	0.999
h95	AIN = −0.0093t2 + 5.1688t + 8.7256	0.998
h115	AIN = −0.0089t2 + 5.1218t + 8.366	0.999
h135	AIN = −0.0087t2 + 4.9069t + 8.5697	0.998
h155	AIN = −0.0081t2 + 4.4461t + 8.2953	0.998

describes the detailed fitting results under different h values with a coefficient of determination (R²) of >0.998. In the same way, h was fitted using the coefficients of t under different subsoiling treatments (R² values were all greater than 0.97), and then substituted into Equation (10) to obtain the mathematical model between AIN, t, and h (Equation (12)).

AIN = (−0.0002h² + 0.0453h + 6.1008)t² + (−0.0002h² + 0.0308h + 3.9044h)t
+ 3 × 10⁻⁷h² − 5×10⁻⁵h − 0.0067

(12)

3.5. Effect of Wing Mounting Height on Soil Moisture Content

After the free infiltration of water into the soil for three hours, the difference in soil moisture content in the top 10 cm below the surface and depths of >30 cm under various h values was small (Figure 7). By contrast, the soil moisture content in the depth of 10–30 cm, when h ≤ 95 mm, was significantly higher than those when h > 95 mm. Both the present research and a study from Yao [23] showed that there was a significant positive correlation between soil water infiltration rate and soil disturbance area. The above results are attributed to the fact that a lower wing mounting height can break more hardpan soil and improve the fragmentation degree of hardpan. As a result, the water can seep out more quickly into the hardpan section. The smaller the wing mounting height of the subsoiler, the faster the soil water infiltration, and the larger the distance of vertical wetting front movement (DVW) and accumulative infiltration (AIN).

3.6. Model Validation

As shown in Figure 8, the simulated (curves) and measured (scattered points) values of soil water contents at various depths (especially the range above 80 cm) after 3 h of free infiltration of soil water were basically consistent. According to the results for soil water contents at different depths from FEM simulations (curve) and field experiments (scatter), RMSEs were lower than 0.05, R² values were higher than 0.95, and mean relative errors were less than 12% (

Table 7. Root mean square error (RMSE), determination coefficient (R²) and mean relative error.

Item	h0	h75	h95	h115	h135	h155
RMSE	0.032	0.026	0.031	0.032	0.036	0.031
R2	0.948	0.982	0.970	0.964	0.956	0.968
Mean relative error (%)	11.03	8.19	10.58	11.87	11.38	9.9

). This indicated that the developed soil water infiltration model had a good accuracy. These errors can be attributed to the following aspects: (1) the soil water infiltration model did not consider the roots and organic matter in the soil; and (2) there are more or less differences in the bulk density and moisture content of soil at different locations at the same depth in the actual field.

4. Conclusions

The effects of the wing mounting height (h) (75–155 mm) of a subsoiler on soil disturbance and soil water infiltration were modelled using HYDRUS-2D and validated using field experiments. The following conclusions were drawn:

• Reducing h values gave larger soil disturbance area ratios, soil water infiltration rates f(t), distances of vertical wetting front movement (DVW), accumulative infiltration (AIN), and soil moisture content at depths of 10–30 cm.
• The relationships among characteristics of soil water infiltration, h and time (t), were developed. The stable infiltration rates (f_s) varied quadratically with h and the corresponding coefficient of determination (R²) was 0.9869.
• The Horton model is more suitable to describe the relationship between f(t) and t under the tested soil conditions, as compared with the Kostiakov and Philip models.
• Overall, reducing the h can improve the accumulative infiltration (AIN) and distance of vertical wetting front movement (DVW) with an increase in time of water infiltration. The relationships among DVW, wing mounting height (h), and time (t), and among AIN, wing mounting height (h) and time (t) were established.
• According to the results for soil water contents at different depths from FEM simulations and field experiments, RMSEs were lower than 0.05 and R² values were higher than 0.95, and mean relative errors were less than 12%. The developed soil water infiltration model had a good accuracy.
• Given the fact that increasing the hardpan disturbance by reducing wing mounting height of the subsoiler could improve soil water infiltration characteristics, it is recommended to appropriately reduce the wing mounting height of the subsoiler before subsoiling. It should be noted that the results obtained in this study are limited to only one soil type (Lou soil) and future studies will be needed to consider more soil types.