Estimating Soil Penetration Resistance of Paddy Soils in the Plastic State Using Physical Properties

Abstract

Soil penetration resistance (SPR) is an important indicator for soil strength which not only affects the growth of crop roots and crop yield but also is crucial in the design and selection of agricultural machinery. The determination of SPR in the laboratory is complex and time-consuming, while measuring SPR on-site shows high uncertainty at different times and locations due to soil heterogeneity. Therefore, this paper investigated the impact of soil parameters on SPR for paddy soils in the plastic state and then established a simple regression model to predict SPR using easy-to-obtain soil physical properties, including clay content, water content and density. Using the combined approaches of central composition rotatable design (CCRD) with response surface methodology (RSM), SPR of 20 soil samples from five paddy fields were measured in the laboratory. The results from the experiments showed that the contribution rate of each single factor to SPR from high to low was soil density, clay content and water content. Statistical analysis for the established equation suggested that the p-value for goodness of fit was significant (p < 0.001) and the p-value for lack of fit was insignificant (p > 0.05); meanwhile, the coefficient of determination (R²) was 0.95, indicating that the model was effective in predicting the SPR. Subsequently, the performance of the regression model was validated by comparing the estimated SPR with in situ field measurements, which showed high accuracy, with percent errors within 10%. Our study successfully proposed a method to estimate SPR using easy-to-measure soil properties that could be obtained from sensors in the soil or field investigations, including soil clay content, water content and wet bulk density.

1. Introduction

Soil compaction resulting from the application of heavy agricultural machinery in intensive agricultural production systems, e.g., field cultivation, fertilization and harvests, has been recognized as being associated with soil degradation and changes in soil mechanical properties [1]. The determination of soil compaction is not only crucial for farmers to assess the necessity of tillage practices but also inevitable for the design and selection of agricultural machinery in different fields for agronomic management practices [2,3].

Soil penetration resistance (SPR) has been regarded as an important indicator to evaluate the degree of soil compaction [4,5], which was found to restrict water uptake by crops [6], affect root proliferation and reduce the growth of crop roots, thereby leading to a reduction in crop yield [7,8]. Souza et al. (2021) stated that high SPR along with large bulk density after the use of heavy agricultural machinery on soil restricted the growth of crop root systems and water availability in deeper soil layers [9]. SPR varied markedly with time and was closely related to changes in soil physical properties [10], such as bulk density [11,12,13], soil water content [13,14,15], soil texture [16,17], soil organic carbon [18], matric potential and degree of saturation [19], etc. Agricultural management practices affect SPR by the disturbance of soil and alteration of soil physical properties. For example, the SPR of paddy soils was found to decrease significantly after puddling due to the reduction of bulk density [11,12]. An increase in tillage intensity contributes to a reduction in SPR by changing the soil bulk density [20], and the installation of subsurface drainage was found to affect soil properties as well as SPR [21]. He et al. (2016) found that SPR was affected by the tine geometry, thickness and penetration depth for sandy loam soil [22].

Although it is common practice to determine the SPR in the laboratory, it is usually expensive and time-consuming, and many researchers have questioned the accuracy of lab measurement for the determination of field parameters due to the differences in loadings [23,24,25]. The measurement of SPR in the field is possible; however, it leads to much uncertainty because the reading of SPR is greatly affected by soil physical properties, which may vary significantly spatially and temporally at different points of one field [26]. Moreover, these approaches are costly due to the need for field samplings and laboratory analyses. Therefore, efforts have been made to use indirect approaches to predict SPR using readily available and costless soil properties for various soil types [9]. Numerous researchers have established different mathematical models to determine SPR using soil physical properties such as moisture content, bulk density, soil texture and organic carbon [27,28,29]. Dexter et al. (2007) presented an exponential equation for the prediction of SPR by an index of soil physical quality and contribution of the pore water to soil strength [17], which showed better performance in predicting SPR than Canarache (1990) and the To and Kay (2005) models [28,29]. However, the model of Dexter et al. (2007) was very complicated, involving the measurement of many soil properties (e.g., bulk density, clay content and several water retention characteristics) [17]. Recently, Zhang et al. (2020) used a mathematical model to estimate SPR from water content, bulk density and shear wave velocity [30]. Whereas the model showed good accuracy (R² = 0.96), the SPR was measured using a small cone penetrometer in a universal test machine laboratory and the model was not validated using data from the in situ measurements, which are necessary to verify the model performance. Gao et al. (2016) proposed a simple model to predict SPR as a function of density, drying and depth in the field on a clay loam soil [31]. However, their model was established using measured data from one soil type and its application on other soil types was not known.

Most of the reported quantitative models have only been established for upland soils with a few certain soil types, and no studies have been conducted on the determination of SPR for the muddy and wet paddy soils, especially those which are in the plastic state. The measurement of some simple soil physical properties, including soil gravimetric water content, clay content and soil wet bulk density, is less complicated, as these data could be easily obtained from soil buried sensors, literature or field investigations. Therefore, the objective of the current study was to propose a simple model to quickly estimate SPR using the soil physical properties that are easy to obtain for paddy soils with high moisture content that are in the plastic state.

2. Materials and Methods

2.1. Sample Preparation

2.1.1. Sample Collection

The experimental area was located in the Yangtze Delta Plain in eastern coastal China, which has been a traditional rice-growing region for over 7000 years, with a large area of paddy fields, covering 11.5% of the total rice planting area in China [32]. The soil samples were taken from five counties in this study region and the soil physical properties for the five soil samples are shown in

Table 1. Measured soil physical properties.

No.	Sand (%)	Silt (%)	Clay (%)	Soil Texture	Plastic Limit (%)	Liquid Limit (%)	Plasticity Index (%)
1	27.9	49.8	22.3	Loam	24.0	51.8	27.8
2	11.3	62.6	26.1	Silty clay loam	23.5	49.0	25.5
3	5.5	63.7	30.8	Silty clay loam	23.5	49.0	25.5
4	4.0	60.0	36.0	Silty loam	31.5	52.5	21.0
5	2.5	56.9	40.6	Silty clay	22.0	49.5	27.5

. The soil texture of these samples was classified using the standard of the American Society for Testing and Materials, D2487-11 [33], including silty clay, loam, silty clay loam and silty loam, covering most of the typical soil texture types in this region. The soil liquid limit, solid limit, as well as the plasticity index were determined according to the standard test methods of the American Society for Testing and Materials [34]. The initial soil water contents for the five soil samples were all higher than the plastic limit and lower than the liquid limit, as shown in Table 1, indicating that the collected soils were in the plastic state.

2.1.2. Central Composition Rotatable Design

To study the impact of clay content, soil water content and density on SPR, the central composite rotatable design (CCRD) was proposed to schedule the experiments [35]. As shown in

Table 2. Coded levels for clay content (X₁), gravimetric water content (X₂) and wet bulk density (X₃).

Code	X1 (%)	X2 (%)	X3 (g·cm−3)
−1.682	22.3	30	1.8
−1	26.1	34	2.0
0	30.8	40	2.3
1	36.0	46	2.6
1.682	40.6	50	2.8

, the soil samples were classified into five levels based on the CCRD, contributing a total number of 20 remolded soil subsamples with different amounts of clay content, gravimetric water content and density (

Table 3. Experiment scheme and test results (X₁: clay content, %; X₂: gravimetric water content, %; X₃: soil wet bulk density, g cm⁻³; Y: soil penetration resistance).

Run No.	${\mathbf{X}}_1$	${\mathbf{X}}_2$	${\mathbf{X}}_3$	Y/kPa
1	−1	−1	−1	35.4
2	−1	−1	1	72.0
3	−1	1	−1	21.3
4	−1	1	1	51.9
5	1	−1	−1	45.3
6	1	−1	1	78.1
7	1	1	−1	44.6
8	1	1	1	60.4
9	−1.682	0	0	44.3
10	1.682	0	0	58.6
11	0	−1.682	0	40.1
12	0	1.682	0	34.2
13	0	0	−1.682	26.3
14	0	0	1.682	74.1
15	0	0	0	48.3
16	0	0	0	50.9
17	0	0	0	49.2
18	0	0	0	51.3
19	0	0	0	49.2
20	0	0	0	48.3

). The soil subsamples encompassed a wide range of soil physical parameters, which were representative of most soils found in the paddy fields of the Yangtze Delta Plain. The gravimetric water content for all was set between 30% and 50%, suggesting that the remolded soils were still in the plastic state.

2.1.3. Soil Sample Preparation

As the soil clay content was fixed for all soil samples (22.3%, 26.1%, 30.8%, 36.0% and 40.6%; see Table 2), the 20 subsamples with different levels of clay content and density were remolded in the laboratory. The calculated amount of water was added to the soil samples to obtain subsamples with different amounts of gravimetric water content (30%, 34%, 40%, 46% and 50%). After complete mixture with the water, the soil subsamples were finally remolded according to following the equation:

$$m=S\times h\times \rho$$

(1)

where m is the mass of soil added to the soil container (g); h is the height of soil in the container for the soil samples (m); $\rho$ is the target density (g cm⁻³); S is the bottom area of the soil container (m²). All soil samples were sealed in plastic bags for 24 h before the SPR tests.

2.2. Soil Penetration Resistance Tests in the Laboratory

The TE-3 soil penetrometer (Nanjing Soil Instrument Factory Co., Ltd, Nanjing, China.) with a measuring distance of 0.2 m was used to measure the SPR. As shown in Figure 1, the soil penetrometer consisted of a measuring head, a measuring rod, a pedal, a handle, a screw rod, a recording pen, a recording tape and a bolt. The height and diameter of the soil container were 0.25 and 0.2 m, respectively. The container with remolded soil samples inside was placed on the self-made aluminum alloy base (Figure 1). The experiments were conducted following these five steps:

Step 1.

The measuring rod was raised until the zero mark of the screw rod was horizontal with the top of the recording disk;

Step 2.

The soil penetrometer was placed vertically in the container, and the handle was fixed with one hand while the manual handle was slowly rotated with the other hand at a constant speed to start the measurement;

Step 3.

The measuring rod was pressed down and the cone probe was wedged into the soil driven by the bevel gear; then, the internal spring was pressed and deformed with data recordings on the paper tape. The ratio between the real soil penetration depth and the value on the recording tape was 2:1;

Step 4.

When the screw on the screw rod collided with the recording disk, the recording pen (separated from the recording paper) was lowered down and the soil penetrometer was pulled out, after which one measurement was completed. Every soil sample was evenly distributed into three subsamples; thus, the measurement was repeated three times for every soil sample.

Step 5.

Repeat steps 1 to 4 for all the 20 subsamples. After the measurement, the soil samples were immediately taken in soil cores and stored in sealed bags for initial water content and density measurements.

2.3. Data Processing and Analysis

2.3.1. Model Establishment

The SPR was expressed by a quadratic polynomial regression using soil clay content (X₁), gravimetric water content (X₂) and density (X₃) as follows:

$$\begin{array}{l}Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+{\beta }_{12}{X}_1{X}_2+{\beta }_{13}{X}_1{X}_3+{\beta }_{23}{X}_2{X}_3+{\beta }_{11}{X}_1^2+{\beta }_{22}{X}_2^2\\ +{\beta }_{33}{X}_3^2\end{array}$$

(2)

where Y is the response value; X₁, X₂ and X₃ are the coded values of each index; and ${\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3,{\beta }_{12},{\beta }_{13},{\beta }_{23},{\beta }_{11},{\beta }_{22}\mathrm{and}{\beta }_{33}$ are the coefficients of the terms in the regression equation.

The results from lab experiments for 20 soil samples were processed and analyzed using the software Design-Expert. The analysis of variance (ANOVA) test was conducted to assess the significance (p ≤ 0.05 level) of all the coefficients in each equation, and the coefficients that were not significant were excluded from the equation. Meanwhile, p-values from both goodness of fit and lack of fit tests were determined, and the determination coefficient (R²) was computed for the regression model for the prediction of SPR.

2.3.2. Response Surface Methodology (RSM)

The response surface methodology (RSM) was used to analyze the results from the quadratic polynomial regression. The RSM is both a statistical and a mathematical technique that has been widely used for process optimization. It allows the explanation of the combined effect of several factors on a response value and it is commonly used with CCRD for second-order response surface models [36]. In the current study, the response surface graph is a 3-D surface graph of SPR responses to soil clay content, gravimetric water content and soil wet bulk density. Then, the 2-D contour plots were projected to provide better analysis of the interactive effect of two factors on the SPR.

2.4. Computation of the SPR

The average SPR (kPa) of each soil sample was determined using the following strategy. As shown in Figure 2, the curve on the recording paper represented the relationship between the deformation of spring and penetration depth at different soil depths. The Y-axis represented the deformation of the spring (x, cm), and the X-axis represented the depth of penetration (h, cm). Subsequently, the average deformation of spring and SPR was computed as follows:

$$\Delta x=\frac{S}{h}$$

(3)

$$P=\frac{\Delta x}{lA}\times k$$

(4)

where $\Delta x$ is the average deformation of the spring (cm); S is the surrounded area of the curves (cm²); h is the depth of the measurement (cm); P is the average SPR (kPa); l is the maximum spring deformation (cm); A is the bottom area of the penetrometer; k is the measuring distance of the spring (kN).

The surrounded area of the curve for each measurement was solved using MATLAB (MathWorks Inc., Natick, MA, USA), before which the closed area outside the curve was blackened, scanned and saved (Figure 3). As shown in Figure 3, a closed square of 1 × 1 cm was drawn on the paper as a reference area, while the closed area 2 was the area to be calculated in Equation (3). The pixels of both areas were computed by MATLAB; therefore, the area of closed area 2 was solved using the proportions between the pixels of the two areas. Finally, the average SPR of each soil sample was obtained through Equation (4).

2.5. Verification for Model Accuracy

As the paddy soils were taken from the fields to the laboratory and the SPR was measured indoors, it was necessary to verify the performance of the established mathematical model in estimating SPR by comparison with in situ measurements on different soils. Therefore, in situ SPR was measured with the soil penetrometer at depth from 0 to 20 cm from the same five paddy fields and an additional field following similar steps to those introduced in Section 2.4 (Figure 1). Subsequently, the in situ measured SPR was compared with the model estimations from soil physical properties using statistical methods of percent bias and coefficient of determination.

3. Results and Discussion

3.1. Experimental Results

The results from experimental tests are shown in Table 3, based on which the statistical analysis was conducted using Design-Expert (

Table 4. Results of statistical analysis for the experiment.

Coefficient	Value	P	F
${\beta }_0$	49.37	<0.0001	22.41
${\beta }_1$	5.26	0.0013	19.67
${\beta }_2$	−4.58	0.0032	14.90
${\beta }_3$	14.37	<0.0001	146.7
${\beta }_{12}$	1.98	0.2313	1.62
${\beta }_{13}$	−2.32	0.1644	2.25
${\beta }_{23}$	−2.87	0.0932	3.44
${\beta }_{11}$	1.74	0.1617	2.28
${\beta }_{22}$	−3.31	0.0167	8.22
${\beta }_{33}$	1.30	0.2855	1.27
P1	<0.0001
P2	0.08
RMSD (kPa)	4.9
R2	0.95

P: p-value for each coefficient; F: F-value for each coefficient; P₁: p-value for the goodness of fit; P₂: p-value for lack of fit; RMSD: root mean square deviation; R²: coefficient of determination.

).

The coefficients of all terms in the quadratic polynomial regression equations are shown in Table 4, in which ${\beta }_0, {\beta }_1, {\beta }_2 and {\beta }_3$ were extremely sensitive (p < 0.005 level), and ${\beta }_{22}$ was sensitive (p < 0.05 level) to the response value of SPR according to the ANOVA test. The remaining coefficients were all insignificant and thus were excluded from the model. Therefore, the quadratic polynomial regression equation for SPR was established as follows:

$$Y=49.37+5.26X_1-4.58X_2+14.37X_3-3.31X_2^2$$

(5)

The p-value from the goodness of fit test was highly significant at p < 0.0001 level and the p-value from the lack of fit test was insignificant at $p\ge 0.05$ level. The root mean square deviation was 20.4 kPa and the coefficient of determination (R²) for the model was 0.95, indicating high accuracy of the established regression equation for estimating the SPR using soil physical properties.

3.2. Impact of Single Factors on SPR

The impact of each factor on SPR was assessed by reducing the dimensionality of the prediction models, during which one factor was assessed while the others were fixed at zero levels in the regression equation. Subsequently, the curves of response values for SPR as affected by each single factor are shown in Figure 4.

As shown in Table 4 and Figure 4, the soil wet bulk density was the primary factor that was positively and linearly associated with SPR. When soil wet bulk density changed from 1.8 to 2.8 g cm⁻³, the SPR increased accordingly from 23.9 to 77.3 kPa. Our results are in line with many published works that suggested positive relationships between soil density and SPR for different soil types. For example, Ayers and Perumpral (1982) showed direct and consistent relationships between soil bulk density and core index, which was a measure of the SPR, for five different soil types with the soil water content below 20% [37]. They indicated that when the soil water content was below 6.7%, SPR increased exponentially with dry density, while it tended to increase linearly with density when the soil water content was between 8.8% and 11.8%. A similar finding has been reported by Vaz et al. (2001), who showed that SPR increased with bulk density (1.2–1.6 g cm⁻³) at the same soil water content (15–30%), and its impact on SPR decreased with higher water content [13].

Similar to soil density, SPR was also positively associated with soil clay content. The impact of clay content on SPR is comparable to previous studies on upland soils. Elbanna and Witney (1987) plotted three curves to describe the relationship between SPR and clay ratio when the moisture content was 20%, 30% and 40%, indicating that SPR increased linearly with clay ratio for soils with different clay contents (from loamy sand to heavy clay soil) [38]. Hummel et al. (2004) used the clay fraction of soil as a major variable for predicting SPR, suggesting that clay content was a significant parameter in predicting SPR for three soil types (Drummer silty clay loam, Plainfield sandy loam and Shoals silt loam) with soil water content below 25% [39].

By contrast, the slope of soil water content and SPR curve was lowest among the three curves, indicating that soil water content has the least impact on SPR as compared with soil density and clay content. As illustrated by Figure 4, SPR was negatively associated with soil water content, which decreased from 78.3 to 28.2 kPa when soil water content increased from 30% to 50% (coded value from −1.682 to 1.682). The reduction in SPR with the increase in soil water content is attributed to the more plastic state of soil through lubrication by water, as well as the rearrangement of soil particles [9]. Many have reported strong correlations between SPR and soil water—for example, the SPR characteristic curve that describes a positive correlation between soil SPR and matric suction [40,41]. Vaz et al. (2001) found that SPR increases exponentially with decreasing water content (15–30%) [13], which was in line with Upadhyaya et al. (1982) [27].

It should be noted that some researchers have pointed out that soil bulk density was a less significant driver of SPR than soil water content and clay content [39]. However, our study indicated that soil density was the primary contributor to SPR, because in the current study, we used soil wet density, which is easier to measure, while the other studies used dry soil bulk density. Therefore, the soil wet bulk density in the current study contributed to the combined effects from both dry bulk density and moisture.

3.3. Analysis of Interactive Factors

In order to gain a better understanding of the results, the predicted models are presented in Figure 5 as 3-D response surface plots and 2-D contour maps.

Keeping one factor among clay content, gravimetric water content and wet bulk density as constant, the interactive impact of other two factors on SPR was analyzed using RSM (Figure 5).

Figure 5a illustrated that when soil density was 2.3 g cm⁻³, SPR was positively associated with clay content, and the maximum value was reached when the coded value of clay content was within 1–1.682 (clay content between 36% and 40%). With the increase in water content, SPR increased first and reached the maximum value when the soil water content was between −1.682 and −1 (30–40%); then, it showed a decreasing trend when soil water content continued to increase. The contour lines suggested that the clay content had a greater impact on SPR as compared to water content, as the slope of SPR was greater in the direction of clay content than water content.

When soil water content was 40%, the SPR increased with clay content and density; meanwhile, its maximum value was achieved when clay content was within 36–40% (coded value 1–1.682) and density was between 2.6 and 2.8 g·cm⁻³ (coded value 1–1.682). The contour map showed (Figure 5b) that the variation rate of SPR was faster in the direction of density than clay content, indicating a greater impact of soil density on SPR than clay content.

When soil clay content was kept at 31%, SPR was positively associated with density, and it increased with soil water content at first but then declined. The maximum SPR fell in the area where water content was within 36–40% (coded value 1–1.682) and density was between 2.6 and 2.8 g·cm⁻³ (coded value 1–1.682). Similarly, soil density contributed more to the SPR than water content.

According to the F and P-values in Table 3 and slopes of contour lines in Figure 5, soil density was the primary factor that affected the SPR, which is followed by clay content and water content. Considering the impact of interactive factors on SPR, their contribution rates from high to low were X₁X₃ > X₂X₃ > X₁X₃. However, p-values of all the coefficients for interactive factors were all insignificant in determining SPR; thus, these coefficients were deleted in the regression models.

3.4. Model Validation

Table 5. Comparison between field measured and model predicted soil penetration resistance (SPR). PBIAS: percent bias.

Location	Clay (%)	Water (%)	Density (g cm−3)	SPR (kPa)
				Measured	Estimated	PBIAS
Hangzhou	22	22.2	1.8	18.6 ± 1.2	16.8	9.7%
Wuhu	26	36.7	2.2	38.6 ± 1.7	34.4	8.5%
Jiaxing	31	50.4	2.0	15.0 ± 0.8	23.1	9.4%
Huzhou	36	40.1	1.9	40.6 ± 1.2	36.7	9.6%
Yixing	40	38.5	2.1	53.1 ± 1.0	51.1	3.8%
Huzhou	35	22.4	2.4	66.3 ± 1.8	62.7	5.4%

presents the measured in situ resistance forces at the depth of 5, 10, 15 and 20 cm from six paddy fields in six different counties of the study region. The SPRs obtained from in situ measurements were then compared with the estimations from the model using soil clay content, water content and density. The percent bias between the measured and predicted SPRs ranges from 3.8% to 9.7%, values which are all within 10%. The comparisons of the measured and model estimated SPR were plotted as the dots in Figure 6, which fitted well in the linear regression line (R² = 0.96). The low bias for each prediction as compared with the measurement as well as the high determination of coefficient suggests that the established model is accurate and efficient in estimating the SPR using readily available soil physical properties.

4. Conclusions

This paper investigated the impact of soil clay content (X₁), water content (X₂) and density (X₃) on SPRs for paddy soils in the plastic state. SPR of the 20 soil samples were tested in laboratory based on the CCRD and analyzed using RSM. Results from the ANOVA test suggested that, though all three soil parameters had a significant impact on SPR, soil density was the primary factor, which was followed by water content and density. The impact of interactive factors on SPR from high to low was ordered as X₁X₃ > X₂X₃ > X₁X₃; however, none of the three interactive factors was significant in determining the SPR. Subsequently, we quantified the SPR using the three soil physical parameters by a regression equation, which showed good accuracy, with significant goodness of fit at the p < 0.001 level and insignificant lack of fit at the p > 0.05 level. Moreover, the coefficient of determination (R²) for the model was 0.95, indicating that the model was effective in predicting the SPR. Furthermore, the performance of the established model for SPR prediction was validated with the field experiment, in which the percent bias between field-measured and model-computed SPR were within 3.8% and 9.7%. Our results proved that the SPR for paddy soils in the plastic state could be quickly estimated through the mathematical method using field-derived soil parameters with high accuracy.