Review of Discrete Element Method Simulations of Soil Tillage and Furrow Opening

: In agricultural machinery design and optimization, the discrete element method (DEM) has played a major role due to its ability to speed up the design and manufacturing process by reducing multiple prototyping, testing, and evaluation under experimental conditions. In the ﬁeld of soil dynamics, DEM has been mainly applied in the design and optimization of soil-engaging tools, especially tillage tools and furrow openers. This numerical method is able to capture the dynamic and bulk behaviour of soils and soil–tool interactions. This review focused on the various aspects of the application of DEM in the simulation of tillage and furrow opening for tool design optimization. Different contact models, particle sizes and shapes, and calibration techniques for determining input parameters for tillage and furrow opening research have been reviewed. Discrete element method predictions of furrow proﬁles, disturbed soil surface proﬁles, soil failure, loosening, disturbance parameters, reaction forces, and the various types of soils modelled with DEM have also been highlighted. This pool of information consolidates existing working approaches used in prior studies and helps to identify knowledge gaps which, if addressed, will advance the current soil dynamics modelling capability.


Introduction
In mechanized agriculture, the energy use for soil tillage operations can be as high as 50% of the total energy used in crop production [1,2]. The energy-use efficiency associated with tillage can be increased by improving the design of tillage implements and through their correct operation and settings [3][4][5]. The design optimization of tillage tools and furrow openers conventionally relies on repeated prototyping and evaluation through soil bins and field experimentation. This task is laborious, time-consuming, and expensive [6,7]. In order to reduce the resource intensity involved, various analytical and numerical models for predicting soil-tool interaction and soil forces have been developed. Analytical models are typically based on the universal earthmoving equation (Equation (1)) [8,9]. P = (γd 2 N γ + cdN c + c a dN a + qdN q )w (1) Figure 1. Discrete element method simulation of a narrow tine furrow opener operating in a moist sandy-loam soil. Modified from Barr et al. [21].
The objective of the work reported in this article was to review the various aspects of the application of DEM in the field of soil dynamics by focusing on soil tillage and furrow opening research for tool performance optimization. Contact models for soil particles; DEM particle size and shape considerations; calibration techniques for determining accurate input parameters; predictions of soil failure, particle movement and reaction forces; and the types of soils modelled with DEM have been reviewed. Knowledge gaps that need attention in future research, and that will go some way to advance soil dynamics modelling capability, have also been identified.

Modelling Agricultural Soils with DEM
The discrete element method is used to model soil as a collection of a finite number of individual spherical particles that interact with neighbouring particles and machine parts when subjected to external forces and forced displacement at the soil-tool interface, such as from a soil tillage tool [7]. This process induces the relative motion of particles within the bulk. Contact forces between these particles and their resultant motion are calculated using Newton's Second Law [22]. These calculations involve repeating the same algorithm at each time step of the simulation process, using results from previous calculations as input.

DEM Contact Models
DEM contact models are developed to describe the mechanical and physical interactions of granular particles with neighbouring particles or external objects. The interactions are modelled using equations of motion and contact models expressed as linear, adhesive, and elastoplastic normal contact models, as well as viscosity, tangential force, and torque models. The tangential force and torque models account for friction, rolling, and torsion [23,24]. Physical interactions between particles are expressed via combining functional elements of springs, dampers, and tangential friction. Total contact forces are expressed as the sum of spring (F s ) and damping (F d ) forces. Some commonly used contact models in DEM simulations are listed in Table 1 and reviewed below. These contact models have been implemented in commercially available software such as Bulk Flow Analyst™, Chute Analyst™, Chute Maven ® , DEMpack™, Altair ® EDEM™, ELFEN, GROMOS-96, ITASCA PFC (2D & 3D), LiGGGHTS ® , MIMES, PASSAGE/DEM, Rocky, SimPARTIX ® , StarCCM+, UDEC, 3DEC, and YADE [25,26]. Table 1 also shows the various DEM soft-

Modelling Agricultural Soils with DEM
The discrete element method is used to model soil as a collection of a finite number of individual spherical particles that interact with neighbouring particles and machine parts when subjected to external forces and forced displacement at the soil-tool interface, such as from a soil tillage tool [7]. This process induces the relative motion of particles within the bulk. Contact forces between these particles and their resultant motion are calculated using Newton's Second Law [22]. These calculations involve repeating the same algorithm at each time step of the simulation process, using results from previous calculations as input.

DEM Contact Models
DEM contact models are developed to describe the mechanical and physical interactions of granular particles with neighbouring particles or external objects. The interactions are modelled using equations of motion and contact models expressed as linear, adhesive, and elastoplastic normal contact models, as well as viscosity, tangential force, and torque models. The tangential force and torque models account for friction, rolling, and torsion [23,24]. Physical interactions between particles are expressed via combining functional elements of springs, dampers, and tangential friction. Total contact forces are expressed as the sum of spring (F s ) and damping (F d ) forces. Some commonly used contact models in DEM simulations are listed in Table 1 and reviewed below. These contact models have been implemented in commercially available software such as Bulk Flow Analyst™, Chute Analyst™, Chute Maven ® , DEMpack™, Altair ® EDEM™, ELFEN, GROMOS-96, ITASCA PFC (2D & 3D), LiGGGHTS ® , MIMES, PASSAGE/DEM, Rocky, SimPARTIX ® , StarCCM+, UDEC, 3DEC, and YADE [25,26]. Table 1 also shows the various DEM software that have been used in tillage and furrow opening research and the types of soil the contact models have been used to model. The movement of particles due to the contact forces are governed by Newton's equation of motion for linear and angular motion as expressed by Equations (2) and (3). By solving these equations, the motion of the particles can be determined. For two spherical particles (i = 1,2) of masses m i and radii r i , located at x i in contact, taking F as contact force, g as acceleration due to gravity, I i as the moment of inertia of a particle, ω i as the angular velocity of a particle, and T i as torque due to the tangential component of the contact force:

Linear Spring Contact Model
This contact model is linearly elastic and is the simplest contact model often used to simulate soil-tool interactions (Table 1) [6,27]. A contact force is created between the two spherical particles in contact as described above. The contact force can be decomposed into normal and tangential force components. The overlap at the contact point generates a repulsive contact force and energy dissipation. When an overlap δ n > 0 is formed between the two particles at a relative velocity . δ n in a direction normal to the contact surface, a normal contact force F n is created based on the spring and dashpot models ( Figure 2) such that: Agriculture 2023, 13, x FOR PEER REVIEW 6 of 30 Figure 2. A schematic diagram of the normal (left) and tangential (right) components of the linear spring contact model.

Hertz-Mindlin Contact Model
The Hertz-Mindlin contact model (HMCM), especially when it is used with the parallel bond model (PBM, see Section 2.1.4.2.), is the most popular contact model used by researchers [7,30,[33][34][35]48] to simulate soil-tool interaction in tillage research. However, as an (non-linear) elastic contact model, it fails to predict vertical soil forces accurately (Table 1) [49]. In this model, the contact force consists of a non-linear Hertz component described by the hysteretic spring force-displacement relationship shown in Figure 3 and a damping component (second part of Equation (7)). It is also resolved into normal and tangential components. The HMCM and its parameters are described in Equation (7) to (16) [20,29].
Normal contact force, Fn: where: Equivalent radius, The parameter k n is the normal stiffness, while d n is the damping coefficient. Considering an imaginary rod of radius r = (r 1 + r 2 )/2 between the centres of the two particles and Young's Modulus E: When the tangential component of the contact force, F t > µ t F n , sliding friction occurs. The local friction coefficient µ t = tan∅, where ∅ is the internal friction angle between the particles.
The tangential component of the contact force is also given by: where k t , δ t , d t , and . δ t are tangential components of stiffness, overlap, damping coefficient, and relative velocity.

Hertz-Mindlin Contact Model
The Hertz-Mindlin contact model (HMCM), especially when it is used with the parallel bond model (PBM, see Section 2.1.4), is the most popular contact model used by researchers [7,30,[33][34][35]48] to simulate soil-tool interaction in tillage research. However, as an (non-linear) elastic contact model, it fails to predict vertical soil forces accurately (Table 1) [49]. In this model, the contact force consists of a non-linear Hertz component described by the hysteretic spring force-displacement relationship shown in Figure 3 and a damping component (second part of Equation (7)). It is also resolved into normal and tangential components. The HMCM and its parameters are described in Equation (7) to (16) [20,29].   [29]. The dashed line contrasts the basic linear elastic model relationship.

Hysteretic Spring Contract Model
The hysteretic spring contact model (HSCM) is an elastic-plastic contact model that accounts for the plastic deformation that occurs during the loading and unloading of soil. It makes the particles behave as though they undergo plastic deformation after the load reaches a yield point, as shown in Figure 4 [50]. The main disadvantage of this contact model is that it requires a large number of input parameters, making its setup and calibration of DEM material properties complex (Table 1) [49]. The HSCM comprises two parts: the spring characteristic illustrated in Figure 2 and damping. A comparative study by Ucgul et al. [29] revealed that the HSCM could model soil-tool interaction more accurately than the HMCM. The HSCM has been used to predict soil reaction forces as well as furrow profiles successfully, especially with the linear cohesion model [39][40][41]43]. The governing equations of the HSCM are described in Equation (17) to (20). Normal contact force, F n : k n = 2E eq r eq δ n (8) d n = ln e ln 2 e + π 2 k n m eq (9) where: Equivalent radius, Equivalent Young's modulus, Agriculture 2023, 13, 541 7 of 29 Equivalent particle mass, Tangential contact force, F t : k t = 8G eq r eq δ n (14) d t = ln e ln 2 e + π 2 k t m eq (15) Equivalent shear modulus,

Hysteretic Spring Contract Model
The hysteretic spring contact model (HSCM) is an elastic-plastic contact model that accounts for the plastic deformation that occurs during the loading and unloading of soil. It makes the particles behave as though they undergo plastic deformation after the load reaches a yield point, as shown in Figure 4 [50]. The main disadvantage of this contact model is that it requires a large number of input parameters, making its setup and calibration of DEM material properties complex (Table 1) [49]. The HSCM comprises two parts: the spring characteristic illustrated in Figure 2 and damping. A comparative study by Ucgul et al. [29] revealed that the HSCM could model soil-tool interaction more accurately than the HMCM. The HSCM has been used to predict soil reaction forces as well as furrow profiles successfully, especially with the linear cohesion model [39][40][41]43]. The governing equations of the HSCM are described in Equation (17) to (20).  Normal contact force, Fn: During loading, During unloading and reloading, Normal contact force, F n : During loading, During unloading and reloading, During unloading again, where k 1 and k 2 are the loading and the unloading stiffnesses, respectively, and e is the coefficient of restitution of the particles, and they are related as e = √ k 1 /k 2 . Tangential contact force, F t : where n k is the stiffness factor equal to the tangential stiffness ratio to normal loading stiffness.

Accounting for Cohesion with DEM Contact Models
In reality, agricultural soils exhibit varying levels of cohesion between particles and adhesion to walls and tools they come in contact with. Attractive pressure (that is, cohesive and adhesive forces) are induced due to the capillary effect and water bridge that exists between particles in unsaturated soils [25,51,52]. Thus, a more realistic contact model for agricultural soils should account for cohesion and adhesion. The linear cohesion and parallel bond models have been used in the DEM modelling of agricultural soils (Table 1).

Linear Cohesion Model
When this model is used, a cohesive or adhesive force is added to the normal force component of the contact model used for cohesionless soils. Even though the linear cohesion model itself does not include a tangential component, its addition increases the normal force, which consequently increases frictional force for greater resistance to slippage [41,50,52,53]. The linear cohesion model can be added to any of the three contact models discussed in Section 2.1.1 to Section 2.1.3 above [50]. If F ca (Equation (21)) is the cohesive or adhesive force, then the normal contact force is modified, as shown in Equation (22).
The parameter r c is the contact radius between particles and can be determined using Equation (23). Equation (21) is called the constant cohesion model because the cohesive stressĉ is a constant. The constant cohesion model makes the model particles too sticky [52]. A modification has therefore been proposed, depending on the degree of compression between two adjacent particles. If compression between two adjacent particles is given by Equation (24), then the cohesive stress increases with time t according to Equation (25).

Parallel Bond Model
An adaptation of the Hertz-Mindlin contact model (HMCM) for cohesive soils is the parallel bond model (PBM) developed by Potyondy and Cundall [54]. The PBM, based on beam theory, uses rectangular (2D) or cylindrical (3D) cement entities as parallel bonds at the point of contact between the two cohesive particles ( Figure 5). This bond is modelled as an elastic beam whose length approaches zero and could be represented by a set of springs uniformly distributed over the contact plane and centred at the contact point [54]. After bond formation, normal and tangential bond forces and moment are calculated in addition to contact forces [50,55,56]. Thus, the bond can withstand or transmit both forces and moments between particles. The bond breaks when its predefined maximum normal or shear strengths are exceeded [57,58]. When no bond exists between particles, the PBM reverts to the HMCM [26]. The PBM is able to model clod formation and the brittle nature of agricultural soils in a more realistic manner [30,59]. It can be used only for particleparticle bonding, not particle-wall (tool) bonding [30,50]. The PBM is the most used model in cohesive soil tillage research [7,30,31,51,[55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70]. Because the base contact model of the PBM is the HMCM, it also fails to predict vertical soil forces accurately as revealed in Table 2 [55,69]. Details of the PBM can be found in Potyondy and Cundall [54]. [52]. A modification has therefore been proposed, depending on the degree of compression between two adjacent particles. If compression between two adjacent particles is given by Equation (24), then the cohesive stress increases with time t according to Equation (25).

Parallel Bond Model
An adaptation of the Hertz-Mindlin contact model (HMCM) for cohesive soils is the parallel bond model (PBM) developed by Potyondy and Cundall [54]. The PBM, based on beam theory, uses rectangular (2D) or cylindrical (3D) cement entities as parallel bonds at the point of contact between the two cohesive particles ( Figure 5). This bond is modelled as an elastic beam whose length approaches zero and could be represented by a set of springs uniformly distributed over the contact plane and centred at the contact point [54]. After bond formation, normal and tangential bond forces and moment are calculated in addition to contact forces [50,55,56]. Thus, the bond can withstand or transmit both forces and moments between particles. The bond breaks when its predefined maximum normal or shear strengths are exceeded [57,58]. When no bond exists between particles, the PBM reverts to the HMCM [26]. The PBM is able to model clod formation and the brittle nature of agricultural soils in a more realistic manner [30,59]. It can be used only for particleparticle bonding, not particle-wall (tool) bonding [30,50]. The PBM is the most used model in cohesive soil tillage research [7,30,31,51,[55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70]. Because the base contact model of the PBM is the HMCM, it also fails to predict vertical soil forces accurately as revealed in Table  2 [55,69]. Details of the PBM can be found in Potyondy and Cundall [54].

Johnson-Kendall-Roberts Cohesion
Another cohesion contact model combined with the HMCM is the Johnson-Kendall-Roberts (JKR) cohesion model [71]. It has been used by researchers such as Cheng et al. [33], Du et al. [36], and Zhai et al. [37] to model cohesive soils in soil tillage simulations with DEM (Table 1). Using this model, the tangential contact force and the normal and tangential damping forces of the HMCM are maintained, while the normal contact force is modified to include cohesion [36]. This modification enables the modelling of strongly adhesive bonds such as exist in dry powders or wet materials (e.g., wet soil). It captures the influence of Van der Waals forces due to contact between two surfaces [50]. A cohesion or adhesion parameter called surface energy is introduced. When this surface energy is zero, the model reverts to the HMCM. The normal contact force in HMCM-JKR is given by Equations (26) and (27), as follows: where a is JKR contact radius and γ is surface energy (J/m 2 ). Sadek et al. [58] n/a n/a n/a n/a n/a

Edinburgh Elasto-Plastic Adhesion Model
Consolidation stress history is one of the main sources of cohesion in cohesive granular materials, and it must be accounted for to accurately model such materials in DEM [78]. The Edinburgh elasto-plastic adhesion model (EEPA) contact model uses a non-linear hysteretic spring model to account for the elastic-plastic contact deformation and an adhesive or cohesive force (pull-off strength) component-acting between dissimilar or similar materials, respectively based on the assumption that this force increases with increasing plastic contact area [50,79]. This model is versatile because, depending on its input parameters, it can be used as either a linear spring model (Figure 6a) or a nonlinear Hertzian spring model (Figure 6b) [79]. Figure 6 shows "a schematic diagram of particle contact and normal force-overlap (f n -δ) curve" for the EEPA contact model. A full description of the EEPA contact model can be found in Morrissey et al. [78]. The EEPA contact model has been used recently for modelling the interaction between tillage tools and agricultural soils [44][45][46].

Edinburgh Elasto-Plastic Adhesion Model
Consolidation stress history is one of the main sources of cohesion in cohesive granular materials, and it must be accounted for to accurately model such materials in DEM [78]. The Edinburgh elasto-plastic adhesion model (EEPA) contact model uses a non-linear hysteretic spring model to account for the elastic-plastic contact deformation and an adhesive or cohesive force (pull-off strength) component-acting between dissimilar or similar materials, respectively based on the assumption that this force increases with increasing plastic contact area [50,79]. This model is versatile because, depending on its input parameters, it can be used as either a linear spring model (Figure 6a) or a non-linear Hertzian spring model (Figure 6b) [79]. Figure 6 shows "a schematic diagram of particle contact and normal force-overlap (fn-δ) curve" for the EEPA contact model. A full description of the EEPA contact model can be found in Morrissey et al. [78]. The EEPA contact model has been used recently for modelling the interaction between tillage tools and agricultural soils [44][45][46].

Particle Size and Shape
Particle shape and size used in DEM significantly affect the necessary simulation time and the accuracy of simulation results [10,28,38]. They are input parameters that should be carefully chosen during calibration for DEM particle interactions to be as close to reality as possible [20,49]. In DEM simulations, it is ideal to use particles of similar sizes to the actual granular materials being modelled. For instance, the actual sizes of agricultural soil particles are relatively small, ranging from several nanometres to about 2 mm for very coarse sand [80]. To model actual particle sizes in DEM requires unrealistically long computation time and impractically high computer processing power [78]. The most timeconsuming part of soil particle DEM simulations is contact detection, and is proportional to the number of particles [81]. For this reason, larger particle sizes than real soil particle sizes are generally employed in DEM [20,38,76]. The larger particles are sometimes implied to represent soil aggregates instead of individual soil particles and somewhat capture the bulk behaviour of a structured soil profile [78].
In reality, soil particles come in various irregular shapes. Thus, particles used in DEM

Particle Size and Shape
Particle shape and size used in DEM significantly affect the necessary simulation time and the accuracy of simulation results [10,28,38]. They are input parameters that should be carefully chosen during calibration for DEM particle interactions to be as close to reality as possible [20,49]. In DEM simulations, it is ideal to use particles of similar sizes to the actual granular materials being modelled. For instance, the actual sizes of agricultural soil particles are relatively small, ranging from several nanometres to about 2 mm for very coarse sand [80]. To model actual particle sizes in DEM requires unrealistically long computation time and impractically high computer processing power [78]. The most timeconsuming part of soil particle DEM simulations is contact detection, and is proportional to the number of particles [81]. For this reason, larger particle sizes than real soil particle sizes are generally employed in DEM [20,38,76]. The larger particles are sometimes implied to represent soil aggregates instead of individual soil particles and somewhat capture the bulk behaviour of a structured soil profile [78].
In reality, soil particles come in various irregular shapes. Thus, particles used in DEM should be not only of a similar size range but also of a similar shape range to actual soil particles to ensure simulations are more representative of realistic bulk behaviour. The basic shape of a DEM particle is a sphere (or circle in 2D modelling) under most DEM codes [20]. Spherical particles approximate and simplify simulations by improving contact detection efficiency and reducing computation time [82]. Usually, irregular (non-spherical) particles are created by clumping a number of spherical particles together, as shown in Figure 7 [28,32,46,83]. This enables the use of spherical particle contact detection algorithms that are simpler and require shorter computation time than those of irregular shapes [24,84].
algorithms that are simpler and require shorter computation time than those shapes [24,84].
Nonetheless, clumped particles also require relatively higher comput than when purely spherical particles are used. Therefore, most studies ado particles to represent soil particles or soil aggregates in DEM simulations. Sp ticle assemblies simulating the soil profile are characterized by considerably nal friction and shear strength than actual soil particles due to lower impa friction. However, this is usually overcome by introducing an arbitrary high tion coefficient to simulate the interlocking tendencies that exist between t shape soil particles [83,85]. . Different associations of spherical DEM particles used to more realistically a effects of soil particle shape on bulk behaviour.

Calibration Techniques for Determining DEM Input Parameters
Running a DEM model involves providing it with such input paramet simulating soil behaviour as close to real soils as possible. Accurate results obtained with accurate input parameters [57,86]. Several approaches exist fo DEM input parameters that accurately represent both soil to soil particle pr soil particle to tool or machine interface properties. The most common calib ods for the former include the angle of repose and hopper discharge, dire triaxial tests, and corresponding in situ soil measurements. The most commo methods for the latter include the inclined plane test, the modified shear tes sponding in situ measurements. All these approaches are focused on bulk re natural stable state and force reactions) of soil under an applied load. After e runs in the laboratory or field, these experiments are then replicated numerica as possible, optimizing parameters iteratively until bulk numerical response field or laboratory measurements [20]. Trial and error methods have tradit relied upon in the past while, more recently, the application of response surf ology (RSM) is demonstrating benefits of significantly reducing the number simulations required for accurate calibration [34,45,[87][88][89][90][91].

Angle of Repose Test
The angle of repose test is used to assess flowability and inter-particle fric soil [92][93][94]. This test is also essential when there is a need for a qualitative a soil surface and furrow profiles in tillage simulations [29]. Various [21,42,84,88,95] used the angle of repose test to calibrate coefficients of static friction between soil particles.
In this test, the soil is allowed to flow by gravity onto a flat surface to Figure 7. Different associations of spherical DEM particles used to more realistically account for the effects of soil particle shape on bulk behaviour.
Nonetheless, clumped particles also require relatively higher computational time than when purely spherical particles are used. Therefore, most studies adopt spherical particles to represent soil particles or soil aggregates in DEM simulations. Spherical particle assemblies simulating the soil profile are characterized by considerably lower internal friction and shear strength than actual soil particles due to lower impact of rolling friction. However, this is usually overcome by introducing an arbitrary high rolling friction coefficient to simulate the interlocking tendencies that exist between the irregular shape soil particles [83,85].

Calibration Techniques for Determining DEM Input Parameters
Running a DEM model involves providing it with such input parameters aimed at simulating soil behaviour as close to real soils as possible. Accurate results can only be obtained with accurate input parameters [57,86]. Several approaches exist for calibrating DEM input parameters that accurately represent both soil to soil particle properties and soil particle to tool or machine interface properties. The most common calibration methods for the former include the angle of repose and hopper discharge, direct shear and triaxial tests, and corresponding in situ soil measurements. The most common calibration methods for the latter include the inclined plane test, the modified shear test, and corresponding in situ measurements. All these approaches are focused on bulk responses (i.e., natural stable state and force reactions) of soil under an applied load. After experimental runs in the laboratory or field, these experiments are then replicated numerically as closely as possible, optimizing parameters iteratively until bulk numerical responses agree with field or laboratory measurements [20]. Trial and error methods have traditionally been relied upon in the past while, more recently, the application of response surface methodology (RSM) is demonstrating benefits of significantly reducing the number of numerical simulations required for accurate calibration [34,45,[87][88][89][90][91].

Angle of Repose Test
The angle of repose test is used to assess flowability and inter-particle friction of loose soil [92][93][94]. This test is also essential when there is a need for a qualitative assessment of soil surface and furrow profiles in tillage simulations [29]. Various researchers [21,42,84,88,95] used the angle of repose test to calibrate coefficients of static and rolling friction between soil particles.
In this test, the soil is allowed to flow by gravity onto a flat surface to form a cone pile. The angle of repose is measured as shown in Figure 8a [29]. Another approach for the angle of repose test is to confine the particles being modelled within the walls of a box, ensuring the top of the particles is levelled. By removing one of the sidewalls, the particles flow to form the angle of repose as shown in Figure 8b [83]. This test is usually used for cohesionless particles and particles with low cohesion with good flowability. However, the general principle is that repeatable observations can be made during key stages of the "angle of repose" experiment. For instance, Roessler and Katterfeld [96] reported successfully calibrating DEM parameters for cohesive soil using the angle of repose test. A cylinder was filled with the cohesive soil, and the cylinder was gradually lifted as shown in Figure 8c. Reproducible phases of soil flow were observed, namely "the build-up of a stable bulk material column, the convex bending of the column, and the beginning of collapse of the column." Aikins et al. [41] observed a reproducible dome-like pile of cohesive soil (Figure 9a,b) and used the results to calibrate soil-soil coefficients of static and rolling friction values.

2023, 13, x FOR PEER REVIEW
cohesionless particles and particles with low cohesion with good flowability. H the general principle is that repeatable observations can be made during key stag "angle of repose" experiment. For instance, Roessler and Katterfeld [96] reported fully calibrating DEM parameters for cohesive soil using the angle of repose test. der was filled with the cohesive soil, and the cylinder was gradually lifted as s Figure 8c. Reproducible phases of soil flow were observed, namely "the build-up ble bulk material column, the convex bending of the column, and the beginning of of the column." Aikins et al. [41] observed a reproducible dome-like pile of cohe (Figure 9a,b) and used the results to calibrate soil-soil coefficients of static and friction values.

Inclined Plane Test
The inclined plane test has been used to determine soil-tool or soil-machin culture 2023, 13, x FOR PEER REVIEW cohesionless particles and particles with low cohesion with good the general principle is that repeatable observations can be made d "angle of repose" experiment. For instance, Roessler and Katterfeld fully calibrating DEM parameters for cohesive soil using the angle der was filled with the cohesive soil, and the cylinder was gradu

Inclined Plane Test
The inclined plane test has been used to determine soil-tool or soil-machine coefficients of static and rolling friction [29]. A schematic diagram of the setup for the inclined plane test is shown in Figure 10. A flat bed of the soil to be modelled is packed into a tray and held on a table with adjustable horizontal inclination. A block of tool material and ball bearings are separately placed on the flat bed, and the table is tilted to an angle Ψ at which the block just starts sliding or the ball just starts rolling down the inclination. The block is used for the determination of the soil-tool coefficient of static friction (µ s ), while the ball bearing is used in the determination of the soil-tool coefficient of rolling friction (µ r ). If the mass of the block is m s , the mass of the ball is m r , and the angles at which sliding and rolling occur are Ψ s and Ψ r , respectively, then the coefficients are calculated according to Equations (28) and (29) [29,97].

Direct Shear Test
The direct shear test is used to determine internal soil parameters namely, and internal friction angle (for soil-to-soil particle interactions). The modified she used to determine the adhesion and external friction angle (for soil to tool or interface properties). These are typically used as direct DEM input parameters to the calibration of other parameters. It has been used for cohesive and adhesive well as cohesionless soils [41,44,53].
This approach uses the normal and shear stresses acting on a column of materials' cross-section. The experimental setup consists of two shear boxes, on on the other and filled with the granular material being modelled. One half is fix the other is made movable horizontally in one direction (Figure 11a). A specified force (Fa) is applied while an increasing horizontal (shearing) force (Fb) is appli movable half till a certain amount of displacement occurs [98]. At that point, the h force would have reached a maximum value and remain constant or slightly in decrease afterward [30]. The experiment is repeated several times with differen force values. Soil-tool coefficient of static friction, µ s = m s g sin Ψ s m s g cos Ψ s = tan Ψ s Soil-tool coefficient of rolling friction,

Direct Shear Test
The direct shear test is used to determine internal soil parameters namely, cohesion and internal friction angle (for soil-to-soil particle interactions). The modified shear test is used to determine the adhesion and external friction angle (for soil to tool or machine interface properties). These are typically used as direct DEM input parameters to support the calibration of other parameters. It has been used for cohesive and adhesive soils, as well as cohesionless soils [41,44,53].
This approach uses the normal and shear stresses acting on a column of granular materials' cross-section. The experimental setup consists of two shear boxes, one placed on the other and filled with the granular material being modelled. One half is fixed while the other is made movable horizontally in one direction (Figure 11a). A specified normal force (F a ) is applied while an increasing horizontal (shearing) force (F b ) is applied to the movable half till a certain amount of displacement occurs [98]. At that point, the horizontal force would have reached a maximum value and remain constant or slightly increase or decrease afterward [30]. The experiment is repeated several times with different normal force values. the calibration of other parameters. It has been used for cohesive and adhesive soils, as well as cohesionless soils [41,44,53].
This approach uses the normal and shear stresses acting on a column of granular materials' cross-section. The experimental setup consists of two shear boxes, one placed on the other and filled with the granular material being modelled. One half is fixed while the other is made movable horizontally in one direction (Figure 11a). A specified normal force (Fa) is applied while an increasing horizontal (shearing) force (Fb) is applied to the movable half till a certain amount of displacement occurs [98]. At that point, the horizontal force would have reached a maximum value and remain constant or slightly increase or decrease afterward [30]. The experiment is repeated several times with different normal force values. In the modified shear box test, the bottom half of the shear box is replaced with a material matching that of the tool or machine part. Normal force and shearing forces are applied the same way as described above. Corresponding normal and shear stresses are plotted as shown in Figure 11b. Given that the cohesive strength of the soil is c, the cross- In the modified shear box test, the bottom half of the shear box is replaced with a material matching that of the tool or machine part. Normal force and shearing forces are applied the same way as described above. Corresponding normal and shear stresses are plotted as shown in Figure 11b. Given that the cohesive strength of the soil is c, the cross-sectional area of the shear box is A, and the internal friction angle is φ, the relationship between the normal force and the horizontal force is displayed in Equation (30).

Triaxial Compression Test
This test can also determine the soil cohesive strength and internal friction angle and used to calibrate DEM input parameters for cohesive-frictional soils [52,72]. A triaxial compression test consists of loading an undisturbed cylindrical soil specimen insulated with an impermeable membrane and subjected to an adjustable confining pressure within a water chamber [18,99,100]. The specimen is then subjected to combined axial (σ 1 ) and radial (σ 3 ) stresses as indicated in Figure 12 until soil failure is achieved [98]. The radial stress (σ 3 ) is first applied around the specimen to a set level via the confining water pressure. An axial strain is then mechanically applied at a controlled rate which generates a corresponding additional deviator stress (q) logged over time and combining with σ 3 to form a resultant axial stress σ 1 . The above steps are repeated several times under increasing radial stress. The plots of the deviator stress (q = σ 1 -σ 3 ) against axial strain identify each deviator stress value at failure and a simple process-for instance using the Mohr circle method-is then used to quantify soil cohesion and internal friction angle [101]. stress (σ3) is first applied around the specimen to a set level via the confining water pressure. An axial strain is then mechanically applied at a controlled rate which generates a corresponding additional deviator stress (q) logged over time and combining with σ3 to form a resultant axial stress σ1. The above steps are repeated several times under increasing radial stress. The plots of the deviator stress (q = σ1 -σ3) against axial strain identify each deviator stress value at failure and a simple process-for instance using the Mohr circle method-is then used to quantify soil cohesion and internal friction angle [101].

In Situ Approaches
Measurements of soil mechanical properties are most accurately done in the laboratory [103]. However, while laboratory methods can precisely measure soil properties, samples may not always be fully representative of field soil conditions due to sampling and handling limitations and time-related changes between sampling and testing. Hence, some researchers have used in situ approaches to measure soil properties for DEM calibration purposes. Kim et al. [44] used an on-site measurement system to determine soil mechanical properties such as shear modulus, Young's modulus, and soil-tool static and rolling friction (Figure 13). On the other hand, Aikins et al. [41] used an on-site cone penetration test (Figure 14) to calibrate Young's modulus of a Black Vertosol.

In Situ Approaches
Measurements of soil mechanical properties are most accurately done in the laboratory [103]. However, while laboratory methods can precisely measure soil properties, samples may not always be fully representative of field soil conditions due to sampling and handling limitations and time-related changes between sampling and testing. Hence, some researchers have used in situ approaches to measure soil properties for DEM calibration purposes. Kim et al. [44] used an on-site measurement system to determine soil mechanical properties such as shear modulus, Young's modulus, and soil-tool static and rolling friction ( Figure 13). On the other hand, Aikins et al. [41] used an on-site cone penetration test ( Figure 14) to calibrate Young's modulus of a Black Vertosol.    Asaf et al. [10] proposed grouser shear and sinkage or penetration tests using wed of different wedge angles and a plate for calibrating DEM contact parameters. Jang e [82] also used a rectangular plate, while Ucgul et al. [38] and Ucgul et al. [29] used circu disc and cone penetration tests to calibrate model parameters. Cheng et al. [33] used a adhesion mass test to determine DEM input parameters of wet clay soil by employing Plackett-Burman test and response surface methodology (RSM) to optimise input para eters. Asaf et al. [10] proposed grouser shear and sinkage or penetration tests using wedges of different wedge angles and a plate for calibrating DEM contact parameters. Jang et al. [82] also used a rectangular plate, while Ucgul et al. [38] and Ucgul et al. [29] used circular disc and cone penetration tests to calibrate model parameters. Cheng et al. [33] used a soil adhesion mass test to determine DEM input parameters of wet clay soil by employing the Plackett-Burman test and response surface methodology (RSM) to optimise input parameters.

Soil Failure and Loosening
Tamas et al. [30] and Barr et al. [21] have revealed the ability of DEM to predict soil rupture and crack propagation, which is an advantage over FEM. Some researchers have used velocity profiles [7,32] or displacement profiles [41,69,104] as soil loosening indicators. Others [21,30] used porosity (in PFC3D Particle Flow Code) or voidage (in EDEM 2.7 TM ), respectively, to measure the degree of particles loosening in DEM. In the work of Tamas et al. [30], for instance, it was found that the DEM modelled soil porosity and soil-break-up resulting from loosening by sweeps increased with both increasing speed and rake angle, which agrees with experimental results.
Identifying particle movement or loosening is mainly used in defining the boundary between disturbed and undisturbed particles to simulate soil failure boundary or furrow profile. Barr et al. [21] argued that using velocity and displacement profiles is based on the assumption that particle movement results in only soil loosening, ignoring the fact that particle movement also occurs during a soil compaction process. The validity of this assumption is therefore limited to tools operating above their critical depth. Additionally, these approaches are open to subjective decisions since a threshold has to be arbitrarily defined to differentiate between the "so-called" loosened and unloosened particles. For example, Murray [69] had to describe loosened particles as having a displacement mag-nitude above 5 mm. Barr et al. [21] instead proposed and used a voidage grid (Figure 15) to define failure boundaries. A voidage grid was applied to the DEM particles after the tillage process was completed and the particles had settled, which reflects experimental practice. Voidage is similar to soil porosity as it measures the proportion of volume not occupied by particles. An empty space will have a voidage of 100%, while a completely filled space will have a voidage of 0%. With V g being grid volume and V p the total volume of particles whose centroids are located within the grid, voidage can be calculated according to Equation (31).
assumption is therefore limited to tools operating above their critical depth. Additiona these approaches are open to subjective decisions since a threshold has to be arbitra defined to differentiate between the "so-called" loosened and unloosened particles. example, Murray [69] had to describe loosened particles as having a displacement mag tude above 5 mm. Barr et al. [21] instead proposed and used a voidage grid (Figure 15 define failure boundaries. A voidage grid was applied to the DEM particles after the age process was completed and the particles had settled, which reflects experimental p tice. Voidage is similar to soil porosity as it measures the proportion of volume not oc pied by particles. An empty space will have a voidage of 100%, while a completely fi space will have a voidage of 0%. With Vg being grid volume and Vp the total volum particles whose centroids are located within the grid, voidage can be calculated accord to Equation (31). Aikins et al. [41] used particle displacement (PD) analysis to determine the loose furrow boundary in DEM. The displacement threshold was not set arbitrarily as was d by Murray [69]. Aikins et al. [41] defined the loosened furrow boundary based on criteria:

= × 100%
1. Minimum particle displacement caused directly by an opener occurs with parti just adjacent to the bottom part of the opener (for wide tines) or particles aligning walls of the slot below critical depth (for narrow tines). 2. To establish a sharp contrast between displaced and undisturbed particles, part locations immediately after particle loosening (i.e., before the particle settle) has to used.
Aikins et al. [41] traced the minimum particle displacement up the profile to prod a loosened furrow boundary as shown in Figure 16a. Figure 16a is a contour plot of Aikins et al. [41] used particle displacement (PD) analysis to determine the loosened furrow boundary in DEM. The displacement threshold was not set arbitrarily as was done by Murray [69]. Aikins et al. [41] defined the loosened furrow boundary based on two criteria:

1.
Minimum particle displacement caused directly by an opener occurs with particles just adjacent to the bottom part of the opener (for wide tines) or particles aligning the walls of the slot below critical depth (for narrow tines).

2.
To establish a sharp contrast between displaced and undisturbed particles, particle locations immediately after particle loosening (i.e., before the particle settle) has to be used.
Aikins et al. [41] traced the minimum particle displacement up the profile to produce a loosened furrow boundary as shown in Figure 16a. Figure 16a is a contour plot of the width and depth of the virtual soil bin profile against displacements (resultant) for each particle within the profile. In some studies [21,32], disturbed soil surface profile after tillage was determined using voidage grid binning and velocity profile. However, Aikins et al. [41] used the profile of the top surface of displaced DEM particles as the disturbed surface profile after the particles had settled because that gives more realistic results and is similar to what actually happens in field experiments. Wang et al. [74] employed another approach using the "clipping" module in EDEM simulation software to define the disturbed soil boundary (Figure 16b). The furrow profile was obtained by connecting the boundaries of the different layers of disturbed soil.
al. [41] used the profile of the top surface of displaced DEM particles as the disturbed surface profile after the particles had settled because that gives more realistic results and is similar to what actually happens in field experiments. Wang et al. [74] employed another approach using the "clipping" module in EDEM simulation software to define the disturbed soil boundary (Figure 16b). The furrow profile was obtained by connecting the boundaries of the different layers of disturbed soil.

Soil Movement and Disturbance Parameters
From the loosened criteria described in the previous section, furrow profile, soil movement, and various soil disturbance parameters have been predicted or determined in DEM simulations with varying levels of relative error (RE) compared to experimental results. Such soil disturbance parameters include the lateral, forward, and upward movement of particles; furrow width at soil surface; loosened furrow cross-sectional area; furrow % backfill and dip area. Barr et al. [21] found an RE of 9% in loosened furrow crosssectional area, 26% in furrow width, 14% in dip area, 0.8% in furrow backfill, 16% in ridge height, and 9% in lateral soil throw in a DEM prediction of soil disturbance parameters with narrow point openers operating in a sandy-loam soil. Barr and Fielke [105] closely predicted lateral soil throw and soil layer mixing using narrow tine openers with 35° and 90° rake angles and a bentleg opener ( Figure 17). Using furrow openers with different rake angles and cutting edge cross-sections, Aikins et al. [41] closely predicted furrow profiles and similar patterns for surface profiles. The majority of DEM predictions of furrow crosssectional area, furrow width, critical depth, and lateral soil throw had an RE of 1% to 20%. However, poor predictions were made for ridge height due to the use of large DEM particles (radius of 5 mm) [41].

Soil Movement and Disturbance Parameters
From the loosened criteria described in the previous section, furrow profile, soil movement, and various soil disturbance parameters have been predicted or determined in DEM simulations with varying levels of relative error (RE) compared to experimental results. Such soil disturbance parameters include the lateral, forward, and upward movement of particles; furrow width at soil surface; loosened furrow cross-sectional area; furrow % backfill and dip area. Barr et al. [21] found an RE of 9% in loosened furrow cross-sectional area, 26% in furrow width, 14% in dip area, 0.8% in furrow backfill, 16% in ridge height, and 9% in lateral soil throw in a DEM prediction of soil disturbance parameters with narrow point openers operating in a sandy-loam soil. Barr and Fielke [105] closely predicted lateral soil throw and soil layer mixing using narrow tine openers with 35 • and 90 • rake angles and a bentleg opener ( Figure 17). Using furrow openers with different rake angles and cutting edge cross-sections, Aikins et al. [41] closely predicted furrow profiles and similar patterns for surface profiles. The majority of DEM predictions of furrow cross-sectional area, furrow width, critical depth, and lateral soil throw had an RE of 1% to 20%. However, poor predictions were made for ridge height due to the use of large DEM particles (radius of 5 mm) [41]. Wang et al. [74] determined soil looseness, furrow width, soil disturbance coe and soil disturbance area ratio with an RE from 3.24% to 41.64% for a winged su Wang et al. [74] determined soil looseness, furrow width, soil disturbance coefficient, and soil disturbance area ratio with an RE from 3.24% to 41.64% for a winged subsoiler and 0.24% to 28.74% for a non-winged subsoiler. There was also a satisfactory agreement in "shape and magnitude" of lateral, forward, and upward displacement of different soil layers resulting from a sweep cultivator between experimental and DEM simulation results [62]. Using a disc with tilt angles from 0 • to 20 • , Murray [69] estimated an average absolute RE of about 10.53% for lateral soil throw. The DEM simulation also revealed, in agreement with the experimental result, that lateral soil throw increased with increasing tilt angle. For a hoe furrow opener, a RE of 14.8% was recorded.
Reduction in the forward movement of soil particles at greater depth has been predicted through DEM simulations [21,32,41]. Other researchers have previously described this phenomenon [9,[106][107][108] through experimental work and analytical models in the concept of critical depth. Below the critical depth, soil movement changes from forward, sideways, and upward directions (generating a loosened crescent-shaped soil failure) to mainly forward and sideways, generating a compaction type failure in a horizontal plane. Barr et al. [21] and Aikins et al. [41] closely predicted the narrowing of furrow down the profile and critical depth with a 90 • rake angle (Figures 15 and 16a). Hang et al. [32] observed a reduction of the forward movement of particles in the inter-row zone between tines with increasing tine spacing in both DEM and experimental results. Greater inter-row zone soil movement with narrower tine spacing was attributed to more intensive interaction between tines. DEM can also simulate greater soil movement with wing attachment to tine cutting tools [31].
Greater soil upheaval observed with a low rake angle opener (35 • ) and increased lateral throw of soil (due to splashing effect) typically found with steeper (90 • ) rake angle openers have been successfully replicated with DEM [21]. However, Ucgul et al. [38] observed that the dynamic height of soil flow under sweep openers was under-predicted by 23% to 35% at speeds of 5 to 12.5 km/h and did not follow the observed shape using spherical particles of 10 mm radii. The prediction was improved by using smaller particles of radii 1.5 mm, which still underpredicted soil flow height by 15%. Using particles of smaller radii, however, considerably increased computation time. Lateral soil throw was also under predicted by about 32% and 9% with radii 10 and 1.5 mm, respectively.
Chen et al. [7] estimated a maximum of 3% RE in furrow cross-sectional area, up to 4% RE in disturbed width, from 14% to 26% more soil (by volume) heaped above the soil surface and 5% to 15% more emptied cross-section below the soil surface. Overall, a close agreement was observed between experimental and simulation results. Hang et al. [32] estimated less than 20% RE between DEM predicted and experimentally determined soil disturbance coefficient and soil looseness. Saunders et al. [76] reported significant underprediction and correlations between measured and predicted furrow area (r = 0.82) and maximum soil throw (r = 0.88) when optimizing the performance of a mouldboard skimmer in a sandy-loam soil. Several furrow parameters from bentleg openers operating in sandy-loam soil, including loosened cross-sectional area (RE = 14.9%), furrow dip area (RE = 14.4%), backfill (RE = 1.8%), ridge height (RE = 16.8%), and lateral soil throw (RE = 14.9%), were accurately predicted using the voidage grid bin approach [39]. Furthermore, Barr et al. [39] observed the same findings as those measured in soil bin investigations, with DEM simulations of bentleg openers also achieving 100% backfill and cancelling furrow spill over.
The DEM has also been used to simulate rotary tiller operations for design optimization. Zhang et al. [109] and Hirasawa et al. [110] closely predicted the height and pattern of soil surface undulations after rotary tillage in DEM. Soil movement pattern during rotary tillage and improving soil layer mixing with tillage depth and travel speed for a rotary tiller have also been closely predicted [36,109]. Cheng et al. [33] recorded an RE of 1.84% in mass of soil that adhered to rotary tiller blades.

Prediction of Tillage Forces
Most attention in the DEM simulation of soil cutting tools has been on predicting soil forces. Table 2 shows relative error (RE) values in draught and vertical force predictions with DEM, travel speed and operating depths for the stated tillage tools, soil types, and DEM contact models. Draught force cycles between peaks at incipient soil failure and troughs at the start of the reloading phase have been well captured in DEM simulations of sweeps by Tamas et al. [30]. In the same DEM study, draught was found to increase with greater speed and rake angle. Draught was predicted with 4% to 12% RE as speed was increased from 0.5 to 2.4 m s −1 with a sweep tine. Bo et al. [31] observed a similar trend between draught force measured for four subsoilers in the soil bin and that obtained through DEM simulation. A winged subsoiler among the four had the highest draught force in both soil bin tests (up to 50% more) and simulations (up to 55% more). With all four subsoilers, DEM predicted draught force with relative errors below 4%.
Additionally, Wang et al. [74] showed that a winged subsoiler operating at a speed and depth of 3 km h −1 and 300 mm, respectively, had an RE of 9.71%, whereas the nonwinged subsoiler obtained an RE of 15.08% when the draught force was compared with the experimental result. Chen et al. [7] observed about 4-31% RE between the draught of experimental and DEM results. A good correlation was obtained between the measured and predicted draught forces (r = 0.95), whereas a more limited correlation was observed for vertical force (r = 0.71). With blunt (R90B) and chamfered (C2S) narrow openers with a 90 • rake angle, a blunt opener with a 45 • rake angle, and a bentleg opener, Aikins et al. [41] predicted draught force with REs of 20%, 22%, 31%, and 5%, respectively. Vertical force was also predicted with 8% and 20% relative error for R90B and C2S, respectively, but poorly for the two other openers.
DEM simulations with a mouldboard plough were also able to simulate the gradual entry of the mouldboard into the soil with a gradual draught increase. For the two soil conditions used in this study, 9% and 2.4% errors in cultivator tool draught were observed for a soft-wet soil and a hard-dry soil, respectively. DEM also closely predicted an increase in draught with increasing depth, with RE ranging from about 3% to 15% [56]. Kim et al. [44] observed that draught force increased with increasing depth with an overall average RE of 7.45%. Tong et al. [73] showed that, across four tillage depths, simulated draught and vertical forces were up to 10% smaller than those measured in the field.
In a DEM prediction of horizontal and vertical soil forces with DEM particles clumped to form different shapes (similar to that shown in Figure 7), Ono et al. [28] obtained the most accurate predictions with the three linearly overlapping spheres. The worst prediction was obtained with simple spherical particles. Ucgul et al. [38] observed a linear increase in draught force against sweep tine width measured experimentally and predicted using DEM with a maximum RE of 8%. Likewise, a non-linear increase in vertical force against width with a maximum RE of 13.7% was recorded. High correlations were recorded between measured and predicted draught forces (r = 0.978) and vertical forces (r = 0.971) with tool speeds from 5 to 12.5 km h −1 and a depth of 70 mm. Prediction of the effect of rake angle on soil forces followed a similar trend (r values of 0.98 and 0.97) and had an RE of 11.6% and 15.2% for draught and vertical forces, respectively. Ucgul et al. [111] again obtained an accurate prediction of draught and vertical forces of a sweep tillage tool at varying speeds and geometry with r values ranging from 0.84 to 0.92. Murray [69] estimated an average RE of 1.86% for draught and 50.7% for vertical force with a flat single disc opener. For rotary tillers, Zhang et al. [109] reported a 12% RE in power consumption, while Du et al. [36] predicted increasing torque with tillage depth (150 to 180 mm) and travel speed (about 2 to 3 km h −1 ).

Soils Modelled in DEM Simulations
Tables 1 and 2 list soil types used in various tillage and furrow opener DEM simulations and their bulk densities, soil water contents, and cohesive strengths. It can be seen that most of the soils modelled with DEM are of sandy to sandy-loam textures. DEM modelling of highly cohesive soils is still relatively scarce in the literature. Bravo et al. [18] used DEM to model highly cohesive clay soil (Vertosol) with cohesive strength of up to about 125 kPa when the soil was highly compacted (bulk density of about 1400 kg m −3 ) and relatively dry (soil water content of about 18%). Aikins et al. [41] also modelled a Vertosol with cohesive strength of 46.4 kPa.
The properties and flow characteristics of sandy soils differ from those of clay soils, which show cohesive and adhesive properties in the presence of sufficient moisture. The successful modelling of a Vertosol and its interaction with tools in DEM by Bravo et al. [18] and Aikins et al. [41] revealed the ability of DEM to model cohesive soils. Although some authors [33,41,43,44,104,112,113] have recently modelled cohesive soils using DEM, more attention is needed in future research to cover the wide spectrum of agricultural soils.

Conclusions
Based on this review, the following conclusions can be drawn:

1.
Even though the Hertz-Mindlin contact model (HMCM) has been used in most DEM studies of tillage and furrow opening, it consistently fails to predict vertical soil force accurately. The Hysteretic Spring contact model (HSCM) can more accurately predict soil forces and particle movement.

2.
Angle of repose, inclined plane, direct shear, triaxial compression, and some in situ tests (grouser shear, plate sinkage, and cone penetration tests) have been used to measure and calibrate DEM input parameters. The angle of repose test has been used mainly for cohesionless soils due to the poor flowability of cohesive soils. However, using results from reproducible phases of the angle of repose experiment, successful calibrations for cohesive soils have been achieved.

3.
Unlike other numerical models, DEM is able to closely predict not only soil forces, but it is also capable of modelling soil failure mechanisms, soil loosening, and soil particle movement. Soil rupture and crack propagation, critical depth, three-dimensional particle movement within the soil profile and lateral particle movement on top of the soil have all been predicted in DEM.

4.
Using voidage or porosity grids to determine loosened furrow cross-sectional profiles has been found to be superior to using particle velocity and displacement profiles. However, some researchers have successfully used a particle displacement approach to determine accurate furrow profiles with a more objective criteria for defining loosened furrow boundary.

5.
Close predictions of draught and vertical forces (≤20%) have been obtained with DEM. These predictions can be improved by using smaller particles of a near-real shape. However, this must be balanced with computation time requirements.
Based on the review conducted, the following recommendations are made for future research: 1.
The Edinburgh elasto-plastic adhesion model (EEPA) has been successfully used to model consolidated or cohesive powders. This contact model is recommended to be studied more extensively for cohesive soils, although some researchers have used it.

2.
Due to pore water pressure within wet and soft soils, coupling DEM and CFD is likely to produce more accurate simulations. This idea can be explored in future research.

3.
A comprehensive analysis of soil disturbance parameters has been successfully done using voidage grids in EDEM ® DEM software. Replication of this approach in other DEM software is recommended. 4.
The criteria introduced by Aikins et al. [41] for defining particle displacement threshold for DEM furrow profile identification need further investigation with particles of smaller radii than the 5 mm used in the study. This approach can provide greater details on the three-dimensional soil translocation process.  Acceleration due to gravity G eq Equivalent shear modulus I i Moment of inertia of a particle k 1 Loading stiffnesses k 2 Unloading stiffnesses k n Normal stiffness k t Tangential component of stiffness m eq Equivalent particle mass m i Mass of spherical particle m r Mass of the ball used in inclined plane test m s Mass of block used in inclined plane test N N factor. Suffixes: γ = gravitational, c = cohesive, a = adhesive, q = surcharge P Soil cutting force (N) q Surcharge stress (Pa) q Deviator stress in triaxial compression test r c Contact radius r eq Equivalent particle radius r i , Radius of spherical particle T i Torque due to the tangential component of the contact force V g Voidage grid volume V p Total volume of particles with centroids within voidage grid w Tool width (m) x i Location of spherical particle γ Specific weight of soil (N m −3 ) ε a Axial strain in triaxial compression test σ 1 Axial stress in triaxial compression test σ 3 Radial stress in triaxial compression test Ψ Inclined plane tilt angle. Subscripts: s = sliding, r = rolling ω i Angular velocity of a particle