Determination of Discrete Element Modelling Parameters of a Paddy Soil with a High Moisture Content (>40%)

: Discrete element modelling (DEM) parameters are of great importance for the accurate prediction of soil properties and disturbance. This study aimed to provide an efﬁcient method for accurately determining the DEM parameters of a paddy soil with a high moisture content (45.66%). The DEM parameters of the paddy soil modelled using the Hertz–Mindlin with JKR contact model were determined by using the Plackett–Burman (PB), steepest ascent, and central composite tests. The accuracies of the developed DEM models were evaluated using actual slumping tests. Based on the PB test, the surface energy of soil, coefﬁcients of soil–soil rolling friction, and coefﬁcients of soil–steel static friction exerted larger inﬂuences on the overall relative error ( δ ZH ). The optimization results showed that to achieve a minimum δ ZH (5.96%), the surface energy of soil, the coefﬁcients of soil–soil rolling friction, and the coefﬁcients of soil–steel static friction should be 0.808 J m − 2 , 0.11, and 0.6, respectively. The optimized DEM model had an overall relative error of 7.27% with a coefﬁcient of variation of 1.32%, indicating that the DEM parameters of the calibrated paddy soil had good accuracy.


Introduction
Rice is one of the main grain crops in China, and paddy field cultivation is an important process in rice production [1][2][3]. Precision tillage in paddy fields is conducive to improving water and fertilizer utilization, irrigation uniformity, and rice yield; reducing weeds, pests, diseases and environmental pollution; and promoting sustainable agricultural development [4][5][6]. Understanding the interaction between flat shovels and paddy soil is the key to optimizing flat shovels and to improving soil disturbance in paddy fields. Previous studies on soil-tool interactions can be grouped into three types: experimental, analytical, and numerical methods [7][8][9][10][11][12]. Experimental methods (e.g., field tests and lab soil bin tests) are time-consuming and costly, while analytical methods are limited by simple geometries and soil failure assumptions [10,11]. Numerical methods are able to examine complex tool geometries and more soil dynamic attributes [12]. Discrete element modelling (DEM) is a numerical method used for simulating the mechanical behavior of discontinuous particles (e.g., soil) without limitations on the magnitudes of particle displacement and tool shapes [10,[13][14][15]. In recent years, DEM has been widely used as a modelling tool in agricultural research. Huang et al. [16] developed a subsoiling tool-soil interaction model using DEM, which could predict the soil's microscopic behaviors at different positions with relative errors less than 20%. Hang et al. [17] investigated the effects of tine spacing of a subsoiler on soil disturbance behaviors (e.g., soil particle force) using DEM and found that tine spacing of 400 mm outperformed other spacings. Ding et al. [18] studied the performance of a deep tillage tool as affected by working depth under near-quasistatic conditions (working speed: 0.1 m/s) using DEM simulations. Tanaka et al. [19] simulated the soil loosening process under the action of a vibrating subsoiling tool based on DEM. Tamas and Bernon [20] simulated sweep-soil interactions using DEM simulations and validated them using field tests. Previous DEM studies mainly focused on simulations of disturbance behaviors of soil with a lower moisture content (<25%). By contrast, the soil in paddy fields before rice transplanting has a much higher moisture content (>40%), cohesion, and fluidity, which result in quite different disturbance behaviors. However, the DEM parameters, which are a prerequisite for simulating paddy soil-tool interactions, are absent in previous literature.
The DEM parameters of soil particles have been calibrated using angle of repose (AOR) tests in many previous studies. Wang et al. [21] determined the DEM parameters of soils with various moisture contents (from 0.27% to 22%) using both funneling angle of repose tests and slumping angle of repose tests. Qi et al. [22] obtained the AOR of a sandy loam soil using the lifting cylinder method and calibrated the soil DEM parameters by matching simulated and measured AORs. Barr et al. [23] calibrated the soil rolling friction coefficient and cohesive energy density when the AOR in DEM reached a good match with that measured. Ucgul et al. [24] determined the restitution coefficient between soil particles and yield strength based on the smallest relative error between forces from DEM simulations and measurements. Yang et al. [25] calibrated soil properties by listing a series of values of various soil properties and determined that the relative error between simulated and measured AORs was minimal. The angle of repose test is generally used to calibrate the DEM parameters of granular dryland soil with low moisture content (≤22%) [21][22][23]. However, the moisture content of paddy field soil is generally larger than 40% according to actual drying experiments. The moving behavior of paddy soil is very similar to that of concrete, and soil cohesion plays a key role in soil movement or flow. To characterize the moving behaviors of highly mobile paddy soil (moisture content: 45.66%) and determine the DEM parameters, a slumping test was proposed using a slumping cylinder with performance indices of the slumping value and slumping expansion. In the slumping test, both the horizontal and vertical soil disturbance behaviors could be quantified.
In this study, the DEM parameters' interval and slumping values and the slumping expansion of paddy soil were measured with physical measurements. On this basis, the Plackett-Burman test, steepest ascent test, and Box-Behnken test were used to calibrate the specific values of the DEM parameters. The accuracy of the DEM parameters after calibration was validated by comparing the simulated and measured slumping values and slumping expansion values. It is expected that this result will provide a basis for the selection of the DEM parameters of paddy soil.

Laboratory Test
To characterize the moving behaviors of highly mobile paddy soil (moisture content: 45.66%) and determine the DEM parameters, a slumping test was carried out using a slumping cylinder in the laboratory, as shown in Figure 1.
The measuring method was based on the Chinese standard "GBT50080-2016: Standard test method for performance of ordinary concrete mixtures". The dimensions of the funnel are shown in Figure 1b. Initially, the paddy field soil was poured into a funnel, and the funnel was then raised vertically at a speed of 0.06 m/s (manually controlled) ( Figure 1a). Eventually, soil flowed out of the funnel under the action of gravity. The slumping value (H) and slumping expansion (W) were measured when there was no obvious movement of soil particles. The W was measured according to the maximum diameter of the soil bottom, while the H was calculated by the difference between the funnel height (H f ) and final soil height (h). The measuring method was based on the Chinese standard "GBT50080-2016: Standard test method for performance of ordinary concrete mixtures". The dimensions of the funnel are shown in Figure 1b. Initially, the paddy field soil was poured into a funnel, and the funnel was then raised vertically at a speed of 0.06 m/s (manually controlled) ( Figure 1a). Eventually, soil flowed out of the funnel under the action of gravity. The slumping value (H) and slumping expansion (W) were measured when there was no obvious movement of soil particles. The W was measured according to the maximum diameter of the soil bottom, while the H was calculated by the difference between the funnel height (Hf) and final soil height (h).

DEM Contact Model
The DEM simulations were performed using EDEM2.7 software (DEM Solutions Ltd. Edinburgh, Scotland, UK). The soil tested was a sandy clay soil that was highly cohesive. The Hertz-Mindlin with JKR (Johnson-Kendall-Roberts) contact model is a cohesive contact model that incorporates the influence of van der Waals forces in the contact area into the calculation. The JKR model uses the same calculation method of the Hertz-Mindlin model to calculate tangential elastic force, normal dissipative force, and tangential dissi-

DEM Contact Model
The DEM simulations were performed using EDEM2.7 software (DEM Solutions Ltd. Edinburgh, Scotland, UK). The soil tested was a sandy clay soil that was highly cohesive. The Hertz-Mindlin with JKR (Johnson-Kendall-Roberts) contact model is a cohesive contact model that incorporates the influence of van der Waals forces in the contact area into the calculation. The JKR model uses the same calculation method of the Hertz-Mindlin model to calculate tangential elastic force, normal dissipative force, and tangential dissipative force. The realization of normal elastic contact force is based on the Johnson-Kendall-Roberts theory. Therefore, the tested soil was modelled using the Hertz-Mindlin with JKR model.

Model Parameter Calibration
Soil parameters in DEM mainly consist of two types: material and interaction properties. Material properties are the soil's intrinsic characteristics, e.g., soil bulk density, shear modulus, shape and size of particles, and Poisson's ratio. Interaction properties are characteristics exhibited by the particle in contact with boundaries, surfaces, and other (or similar) particles, e.g., surface energy, coefficients of restitution, and static and rolling fric-tion [21]. The nominal soil particle radius was 2 mm, and the size range was 0.8-1.2 times the nominal particle size, as determined by a combination of preliminary tests and the literature [7,10].
Soil bulk density was determined by measurement. The shear modulus, Poisson's ratio, and density of steel were determined according to previous studies [7,16]. The ranges of shear modulus, surface energy, and Poisson's ratio of soil, as well as the coefficients of restitution, static friction, and rolling friction of soil-soil and soil-steel, were selected from the literature (Table 1) [26][27][28][29][30][31]. The parameters that significantly affect slumping value and slumping expansion of the paddy soil were screened out using the Plackett-Burman (PB) test. The steepest ascent (SA) test was used to further narrow the value ranges of significant parameters [21]. A second-order regression model between indices (i.e., slumping value, slumping expansion, and overall relative error) and significant parameters was then established using a central composite test (CCT). The overall relative error (δ ZH ) was used to evaluate the overall error between simulated and experimental slumping values and slumping expansion. It was calculated as follows: where H s and H e represent the simulated and experimental slumping values in units of mm, and W s and W e represent simulated and experimental slumping expansion in units of mm. Finally, the optimal combination of significance parameters concerning the indices' optimization was calculated according to the second-order regression model.

Laboratory Test Results
The slumping test in the laboratory was repeated five times. As shown in Table 2, the slumping values ranged from 230.17 mm to 238.67 mm with a mean value of 233.5 mm. By contrast, slumping expansions ranged from 344.17 mm to 360.00 mm with a mean value of 350.8 mm. The coefficients of variation of slumping values and slumping expansion from different repetitions were 1.24% and 1.77%, respectively, indicating that the experimental results were reliable.

Plackett-Burman (PB) Test
The results of the Plackett-Burman (PB) test are shown in Table 3. According to the ANOVA outputs (Table 4), the contributions of X 2 , X 6 , and X 8 to overall relative error (δ ZH ) were relatively higher than those of the other parameters, so these three parameters were calibrated using the following steepest ascent (SA) and central composite tests. Table 3. Plackett-Burman test scheme and results. No.

Steepest Ascent (SA) Test
To further investigate the feasible range of soil parameters (i.e., X 2 , X 6 , and X 8 ), the SA test was designed and carried out in accordance with their effects on overall relative error (δ ZH ) (see Table 4). As shown in Table 5, the lowest relative error (9.91%) was at the fourth simulation, implying that more reasonable ranges were 0.60-1.10 for X 2 , 0.05-0.09 for X 6 , and 0.60-1.10 for X 8 . The obtained feasible ranges of soil parameters were then utilized in the following central composite test.

Central Composite Test (CCT)
The specific effects of the surface energy of soil (X 2 ), coefficients of soil-soil rolling friction (X 6 ), and coefficients of soil-steel static friction (X 8 ) on the slumping value (H) and slumping expansion (W), as well as the relative errors between the simulated and measured values of both H and W (i.e., δ 1 and δ 2 ), were investigated using a central composite test. The test design and results are shown in Table 6. The regression equations between the indicators (i.e., H, W, δ 1 , and δ 2 ) and X 2 , X 6 , and X 8 were obtained as follows:  The coefficients of correlation were 0.9778, 0.9904, 0.9778, and 0.9873 for Equations (4)-(7), respectively. Moreover, the residual standard deviations were 1.4993, 2.3696, 0.6421, and 0.5992 for Equations (4)-(7), respectively. The ANOVA outputs (Table 7) showed that Equations (4)- (7) can be used to describe the relationships between indicators (H, W, δ 1 and δ 2 ) and X 2 , X 6 , and X 8 .

Parameter Optimization and Validation
To obtain the optimal combination of soil DEM parameters, the minimum overall relative error (δ ZH ) was used as the target; the constraints were the feasible ranges of soil parameters determined from the SA test. The final optimization results showed that to achieve a minimum δ ZH (5.96%), X 2 , X 6 , and X 8 should be 0.808 J m −2 , 0.11, and 0.6, respectively. To validate the optimization results, five slumping tests were carried out in DEM simulations with the optimized combination of X 2 , X 6 , and X 8 (Figure 2). The results showed that the overall relative error (δ ZH ) was 7.27% with a coefficient of variation of 1.32% ( showed that the overall relative error (δZH) was 7.27% with a coefficient of variation of 1.32% (Table 8), indicating that the calibrated soil DEM parameters had good accuracy.

Conclusions
In this paper, the DEM parameters of a paddy field soil modelled using JKR were determined by a combination of the Plackett-Burman (PB), steepest ascent (SA), and central composite tests. The accuracies of DEM models were evaluated using slumping tests. The following conclusions were drawn: (1) For a given paddy field soil, the surface energy of soil (X2), coefficients of soil-soil rolling friction between (X6), and coefficients of soil-steel static friction (X8) have a greater influence on the overall relative error.

Conclusions
In this paper, the DEM parameters of a paddy field soil modelled using JKR were determined by a combination of the Plackett-Burman (PB), steepest ascent (SA), and central composite tests. The accuracies of DEM models were evaluated using slumping tests. The following conclusions were drawn: (1) For a given paddy field soil, the surface energy of soil (X 2 ), coefficients of soil-soil rolling friction between (X 6 ), and coefficients of soil-steel static friction (X 8 ) have a greater influence on the overall relative error.
(2) An optimal combination of X 2 of 0.808 J m −2 , X 6 of 0.11, and X 8 of 0.6 was obtained given overall relative error minimization (5.96%).
(3) A low relative error (7.27%) was found for simulated overall relative error using a developed DEM model of paddy field soil, indicating that the calibrated DEM parameters were reliable. Paddy soil-flat shovel interactions using DEM simulations will be investigated in our future research.