Historical evolution and new trends for soil-intruder interaction modeling

Soil is a crucial resource for life on Earth. Every activity, whether natural or man-made, that interacts with the sub or deep soil can affect the land at large scales (e.g. geological risks). Understanding such interactions can help identify more sustainable and less invasive soil penetration, exploration, and monitoring solutions. Over the years, multiple approaches have been used in modeling soil mechanics to reveal soil behavior. This paper reviews the different modeling techniques used to simulate the interaction between a penetrating tool and the soil, following their use over time. Opening with analytical methods, we discuss the limitations that have partially been overcome by the finite element method (FEM). FEM models are capable of simulating more complex conditions and geometries. However, they require the continuum mechanics assumption. Hence, FEM analysis cannot simulate the discrete processes occurring during soil deformation (i.e. the separation and mixing of soil layers, the appearance of cracks, or the flow of soil particles). The discrete element method (DEM) has thus been adopted as a more promising modeling technique. Alongside models, experimental approaches have also been used to describe soil-intruder interactions, complementing or validating simulation results. Recently, bioinspired approaches have been considered promising to improve sustainability and reduce the invasiveness of classical penetration strategies. This review highlights how DEM-based models can help in studying the interaction mechanisms between bioinspired root-like artificial penetrometers and the soil. Bioinspired designs and the merging of multiple analysis approaches can offer new perspectives. These may be pivotal in the design of highly optimized soil robotic explorers capable of adapting their morphology and penetration strategies based on their surrounding conditions.


Introduction
The soil is key to better exploiting natural resources, building structures, and studying geotechnical aspects. Mechanical penetration systems have been developed so that sensorized probes can reach a certain depth in the soil [1]. To dig into soil, drilling devices, vibration systems, pile driving, tunneling systems, and so on are typically used. However, classical excavators create heat during cutting/shearing and use lubricants to remove soil particles. The above techniques are bulky and highly invasive [2]. Soil pollution during excavation is a considerable drawback during monitoring water, minerals, nutrients, and contaminants in the soil [3,4]. Less invasive systems, such as penetrometers, are adopted to monitor and extract the characteristic parameters of the subsoil. These are sensorized metallic probes that are manually, or by mechanical actuators, pushed into the soil from the surface [5]. They typically move following a straight pathway while reaching a certain depth. However, flexible and adaptable probes might be needed when moving underground to avoid collisions and damage in the case of underground heritage and other infrastructures or foundations. To reduce invasiveness and increase dexterity, biologically-inspired solutions have started to be investigated [6][7][8]. Living organisms adopt efficient strategies for digging and exploration, and represent a powerful source of ideas for the design of novel artificial digging systems [9][10][11][12][13]. Many natural organisms, including both animals and plants, can penetrate hard and stiff materials, anchoring in soft soils and moving in loose granular mediums [14,15]. Since artificial systems undergo the same physical laws as natural systems when moving at the same scale in analogous environments, bio-inspiration can provide novel solutions to address geotechnical challenges [13].
Along with many living organisms, plant roots are efficient explorers in soil [16]. They adopt multiple strategies to move in this challenging environment, including mucus and dead cell release (adopted as a lubricant) [17,18], thickening of the body [19], lateral roots and hairs (to improve anchoring) [20], growth from the tip (possibly reducing friction) [4], circumnutation movements (to possibly decrease the soil pressure underneath) [21,22], and exploitation of soil cracks (to move in low-resistance pathways) [23].
Typically, experimental approaches have been adopted to study the penetration strategies used by plant roots with a twofold objective. On the one hand, to verify and validate the efficient employment of the investigated strategy by the biological model and, on the other, to develop novel penetrating devices. However, the right selection and combination of multiple penetrating strategies, the scaling of devices to reach the biological counterpart, and realistic environmental conditions are just a few of the challenges encountered in the implementation and experimentation. In this context, simulating the behavior of a penetrating device in the soil could be a powerful tool to guide the design of efficient artificial systems for underground exploration.
Describing this process appropriately is challenging due to the complex physical processes involved, such as (a) the spatial variation of soil media, (b) the nonlinear behavior of soil material, and (c) the dynamic effects and flow induced in the soil by this interaction [24]. To approach these issues, several modeling techniques have been explored over the years to improve capturing soil behavior during intruder interaction. Abstracting from a specific field of application, we are interested in underlining the limits and benefits of different modeling approaches and discussing their evolution induced by the need to capture a more faithful representation of soil behavior when subjected to a body intrusion or penetration.
Here, we review how the modeling approaches to this issue have evolved over time, the reasons for the switch between one technique and another, and their possible use to complement and support the evolution in geotechnical approaches with bioinspired mechanisms focusing on root-like penetrating models. Figure 1 shows the main contributions we review in this paper, starting with analytical models (detailed in section 2), followed by finite element method (FEM)-based models (described in section 3) and discrete element method (DEM) analyses (discussed in section 4). An overview of root-like models is also reported (section 5), followed by the discussion and conclusions (section 6).

Analytical method
The interaction between penetrating tools and the soil is being widely investigated in different application fields, such as geotechnical engineering problems, agriculture, and the locomotion of robotics systems.
This section reviews how this interaction process has been addressed independently on the considered tool and the application field.
Initially, the tool-soil interactions were investigated using analytical methods based on the passive soil pressure theory, first proposed by Terzaghi [25], and assuming that the soil-tool interaction effect could be modeled as a soil failure pattern.
Along with many others, Hettaratchi et al [26] developed a model to compute the soil load acting on the interface of a penetrating tool to analyze the soil-tool interaction mechanics. The soil load was expressed as a function of the gravitational, cohesive, adhesive properties, and surcharge components of soil reaction per unit width of interface. The proposed model is valid for all typologies of soil failure where deformations occur slowly so that inertia factors may be neglected [26].
In 1975, Durgunoglu and Mitchell [27] analyzed the failure mechanisms associated with the static penetration resistance in cohesionless and cohesive soils. In particular, the authors predicted the penetration resistance encountered by penetrometers differing in tip geometries and surface roughness, and they treated the penetration process as a bearing capacity problem. In particular, the penetration resistance was assumed to depend on cone apex angle and roughness, friction angle of the soil, and penetration depth to base width ratio [27].
McKyes and Ali [28] proposed a threedimensional (3D) model consisting of straight-line failure patterns in the soil to predict the draft forces and the volume of the soil disturbed when subjected to a narrow blade. The authors adopted force factors which, in contrast to those considered by Hettaratchi et al [26], varied with (a) the width-to-depth ratio of the instrument, (b) the rake angle of the instrument, and (c) the friction angle of the soil. In addition, the geometry of the 3D failure pattern in the soil was approximated for different blade shapes and soil strengths, designing the tools based on their draft force requirements and soil cutting efficiency [28].
Perumpral et al [29] suggested a mathematical method to estimate slow-moving narrow-tillage tools in cohesive-frictional soils and, in particular, the draft and the vertical forces the tool is subjected to.
The analytical models reported above have some limitations. As stated by Chi and Kushwaha [30], these models approximate the intruder with one center wedge and two side crescents and, to estimate the soil forces acting on the tool, adopt simple equations, which are developed using limit equilibrium principles. For instance, the models reported in [28,29] consider a plane failure surface at the bottom of the tool to solve the limit equilibrium equation. However, processes occurring in real situations appear very different. When the initial cutting starts, the shear failure induces a failure in front of the penetrating tool. Then, during the cutting, there is a continuous change in (a) the soil's configuration ahead of the tool, (b) the soil-tool interface and (c) the relationship between relative displacement and shear stress [30]. All these properties are not considered in the analytical models described above.
Modeling the effect of the soil-tool interaction as a soil failure pattern is a significant limitation. In fact, the soil failure profile depends on the tool shape, the operating speed, and the soil's physical properties. Analytical models are therefore suitable for predicting the information related to tillage performance (e.g. the draft of an existing tool) but are not sufficiently accurate for designing new tools to improve performance [31]. The tool response and the final soil conditions are affected by e.g. the initial soil conditions, the tool configuration and the operating conditions. Since analytical models consider penetrating tools with a simplified geometric shape, the model's output is significantly affected by this approximation. The tool response is thus evaluated in terms of tool draft, energy requirement, and tool wear in the soil-tool interaction modeling. These characteristics depend on the operating speed and time, as tillage operations are dynamic processes [31]. Furthermore, analytical models (a) do not provide enough information for optimizing the design of tillage tools, because they are not able to predict the fields of stress, displacement, speed, acceleration of soil-tool interactions, and (b) in the assumption of a preliminary soil wedge failure, the soil inertial effect can be analyzed well but there is no well-established theory to consider the soil strength effect on the draft [31]. In this context, cone penetration problems were often investigated. A popular approach was the cavity expansion theory, which assumes that the pressure needed to produce a deep hole in a medium is proportional to the pressure required to expand a cavity of the same volume under the same conditions [32]. Two different typologies of cavity expansion problems can be discussed [32]. One approach assumes the cavity to exist initially and to be characterized by a pressure that is in equilibrium with the stress of the surrounding soil and that, for the cavity expansion, needs to increase until reaching a limit or a steady-state value. By contrast, the other approach does not consider an initial cavity. Therefore, the cavity expansion consists of expanding an existing cavity with a null radius until a final radius approaches infinite. When a cavity expands from an initial radius R 0 to a final radius R f in a certain region, the stresses in the surrounding soil lead to the creation of three different regions [32] as shown in figure 2: (a) a linear elastic zone, where the soil behaves as a linear elastic material because the stresses are small; (b) a non-linear elastic zone, where the soil yields but the stress are not large enough to induce soil failure; (c) a plastic zone, where the stresses are significant and induce a failure condition.
In 1997, Salgado et al [32] developed a theory based on the cavity expansion analysis to estimate the resistance experienced by a penetrometer tip during penetration in a cylindrical sample of silica sand. Good consistency between estimated and measured resistance values was found (the differences were lower than ±30) [32]. Since most studies using the cavity expansion theory were based on an elasticperfectly-plastic model, new models (e.g. camclay and modified camclay) were proposed to take into account the plastic hardening and softening properties of the soil. Moreover, the cone penetration was investigated in different ways (i.e. as a bearing capacity problem, as a steady state problem and using large calibration chambers) [33]. In 2001, Cao et al [34] analyzed the undrained expansion of a cavity in modified camclay by considering a combination of large-strain theory and small-strain theory in the plastic and elastic zones, respectively. The authors accurately predicted both radial and circumferential stresses around the cavity and the pore pressure in the proximity of the cavity wall. When the cone penetration was modeled through the cavity expansion theory, both spherical (e.g. [35]) and cylindrical (e.g. [32]) cavity expansion approaches were used: cylindrical cavity expansion uses cylindrical coordinates, whereas the spherical cavity expansion is formulated using spherical symmetry conditions. In this context, a combined cylindrical-spherical cavity expansion approach was proposed by Yu [36] to compute the cone tip resistance. Later, Salgado and Prezzi [37] evolved a previous developed analysis [32] to compute the cone penetration resistance in the sand, considering the classical cavity expansion analyses based on linear elasticity and perfect plasticity. In 2012, Salgado [38] reviewed the previously adopted approaches (both theoretical and experimental) to investigate the cone penetration boundary value problem (BVP). Theoretically, the cavity expansion analysis (detailed above), the strain path method [39], and numerical analysis (e.g. the Arbitrary Lagrangian-Eulerian (ALE) formulation of the FEM) were used to investigate the cone penetration BVP. Experimentally, it was analyzed in field, chamber and centrifuge [38]. Experimental analyses better clarified some aspects (e.g. dependence of the cone resistance on the penetration rate) with respect to theoretical investigations [38].
In 2013, Li et al [40] investigated the locomotion of a small legged robot on granular soils approximating the net forces acting on the leg in the vertical plane by the linear superposition of resistive forces on infinitesimal leg elements. The acting stresses depended on the intruder depth, orientation, and movement direction. In this model, the researchers took inspiration from the resistive force theory (RFT) method, which was widely adopted to investigate intruders' locomotion. This approach assumes to portion the intruding body into infinitesimal segments, each of which generates thrust and experiences drag, and to compute the net force on a segment by considering the normal and tangential components, which depend only on local properties (i.e. element length, velocity, and orientation) [41]. Figure 3 shows the main idea considered by the RFT: a body in motion propagates a wave in the negative x direction and is propelled forward due to reaction forces it experiences. In particular, each infinitesimal segment of the body moving with velocity v experiences forces dF ⊥,// in both tangent t and normal n directions [41]. If the motion occurs in fluids, these forces are measured through the Stokes' law; if the motion occurs in granular media, these forces are measured experimentally.
The RFT method was popular due to its simplicity. However, it was restricted to analyses considering dry granular media, and it considered the steady state drag forces to characterize the forces experienced by each segment [41]: when locomotion in granular media is analyzed, this approximation is not accurate because the material varies state even when low constant intrusion speed are considered.
In 2018, Kang et al [42] investigated, both experimentally and theoretically, the penetration process done into granular soils by intruders with different geometry. In particular, the authors observed that the penetration of any convex intruders could be understood following a modified Archimedes' law. These findings provide important insights for different application fields, such as cone penetration tests, locomotion of biological and robotic systems in the soil, and digging on Earth and other planets [42].

FEM
Limitations shown by the analytical models are partly overcome by the FEM, which represents a suitable numerical technique for analyzing complex engineering problems characterized by geometric and material nonlinearity, as in the case study of the soiltool interaction [31]. In the following, section 3.1 summarizes the constitutive relations that have been proposed for FEM model formulation, followed by section 3.2, where FEM models more strictly related to soil-intruder interaction are reviewed.

Soil constitutive relation
The soil constitutive relation is at the base of FEM models formulation. Many characteristics of soil response cannot be captured through elastic models. Indeed, most soils experience important volume changes even when subjected only to shear stress variation. When the shear stress becomes sufficiently large, they reach a state of continuing shearing with no additional changes in stresses at large strains [43]. This behavior can be well captured by an elasticperfectly-plastic description of the soil response, which divides the strain increments occurring with stress variations into elastic and plastic components [43]. In particular, these models consider a region of stress space that can be reached elastically. When the boundary of this region is reached, the material yields at constant stress [43]: the boundary of the elastic area is known as the yield surface, and a yield function mathematically describes it.
Among the elastoplastic constitutive models, camclay and modified camclay models were developed for predicting the response of soil samples subjected to triaxial tests [44]. These models were applied mainly in the context of the loading of clay samples and geotechnical constructions on clay [44]. In 1990, Wood [44] provided a stress-dilatancy relationship and an expression of plastic potential describing the yield loci for the soil. Considering this yield loci and identical plastic potentials in the framework of volumetric hardening elastic-plastic soil behavior models lead to recreating the original camclay model, which was developed by Roscoe and Schofield [45]. The modified camclay model represents a later and more adopted version of the earlier camclay model and is developed using four main ingredients: elastic properties, yield surface, plastic potential, and hardening rule [44]. Moreover, the modified camclay model is based on two main assumptions [44]: (a) any change in mean effective stress (p) occurs together with recoverable changes in volume, and (b) recoverable shear strains accompany any change in deviator stress (q). Circular or elliptic shapes are considered in the yield locus in the p-q plane, and the yield locus is assumed to expand at constant speed. This model provides the methodology to verify if plastic strains are occurring and expressions of elastic and plastic stress-strain responses, which are generical and valid for all effective stress paths in the p-q plane. In this context, it is important to highlight the concept of critical state: it refers to a condition in which the plastic shearing could continue indefinitely without implying variations in volume or effective stresses (i.e. condition of perfect plasticity) [44]. An example of elastoplastic constitutive models is given by the model developed by Yao et al [46]. They predicted the mechanical behavior of granular soils (i.e. stress-strain relation, positive/negative dilatancy, and strain hardening/ softening), and validated the modeling results by comparison with experimental findings obtained from tests on Cambria sand and Coarse-grained material. The model performed well in a large range of stresses [46]. Among the non-linear hypoelastic models, the Duncan-Chang model (1970) [47] was widely adopted to describe the soil-tool interaction because simple and easy to implement.
In addition, several models based on the theory of hypoplasticity were developed. Hypoplasticity represented an alternative to the classical theory of elasto-plasticity (e.g. a different concept of yielding) [48]. Among hypoplastic constitutive models, Wu et al [49] proposed a model without the elastoplastic theory (i.e. yield surface, plastic potential, flow and hardening rules, and decomposition of deformation in elastic and plastic parts). In particular, the authors included the effect of void ratio and stress level by integrating the critical state in the model and validated it by performing several numerical experiments finding a good consistence for reproducing the main behavior of granular materials [49]. With respect to previous constitutive models, hypoplastic models are simpler to formulate and apply [50], despite the need to include a broader range of application fields (e.g. cohesive materials), and result in a more accurate numerical implementation [50].
In 2004, Dafalias and Manzari [51] proposed some improvements to a previously developed sand plasticity model: the dilatancy expression as a function of a fabric-dilatancy tensor was introduced to overcome the low accuracy observed in the simulation of reverse loading when low effective stresses were considered despite the bounding surface formulation; then, the plastic potential was expressed as a function of a modified Lode angle in the multiaxial generalization to simulate the deviatoric plastic strain under radial monotonic loading properly. The authors considered the simple triaxial and the general multiaxial formulation to be connected systematically and investigated several density and pressure values in both drained and undrained loading conditions [51]. Dafalias and Manzari [51] highlighted the importance of including a fabric-dilatancy tensor to have a realistic description of soil response in reverse loading and provided a concept of fabric-dilatancy tensor that can be considered generical and, therefore, appliable in other constitutive models. This model belongs to a category of constitutive models known as Simple ANIsotropic SAND, which were developed for the critical state soil mechanics and the bounding surface plasticity [52]. In this context, Taiebat and Dafalias [52] proposed a new constitutive model to simulate both drained and undrained behaviors of sand in a wide range of confining pressures. This model successfully simulated the behaviors in both conditions, and respect to previous works, it induced plastic deformations under constant stress-ratio loading, and it did not require initial values of stress ratio [52].

Soil-intruder FEM models
Unlike analytical models, FEM models are able to predict the draft without a preliminary assumption on the soil failure pattern and can compute the fields of displacement, stress, speed, and acceleration of soiltool interactions [31].
When a tool moves in the soil, the soil-tool interface shifts continuously, introducing a discontinuity in the displacement field. For this reason, an idealized joint element (i.e. interface element) is considered between the tool and the soil when the tool response is simulated [31].
The joint element was first proposed by Goodman et al [53] and then adopted in FEM models that simulated interacting systems. Goodman et al proposed a two-dimensional (2D) joint-element with zero width, four nodes, and eight displacement degrees of freedom enabling normal and tangential deformations [31,53].
Desai et al [54] proposed substituting the joint at zero width with a thin-layer element. This solution (a) enables both 2D and 3D problems to be represented, (b) provides an improved definition of normal and shear behavior, (c) presents various deformation modes if appropriate stress redistribution schemes are included, and (d) it is a more realistic representation than the zero thickness element since a thin layer of soil is often involved in the interaction behavior [54].
In 1995, a thin interface element was also recommended by other researchers [55] to overcome the limitations observed when zero thickness interface elements are considered [53] such as: tangential force oscillations, due to inappropriate quadrature schemes, which affect the element's performance; the use of extremely large stiffness parameters which lead to a lack of accuracy; an insufficient mesh fineness which induces inaccurate interface stress predictions [55].
In 1977, the FEM was employed to model the soiltool interaction by Yong and Hanna [56]. The authors developed a 2D FEM model to stimulate soil cutting under a wide blade. The model predicts the stress and the deformation fields in the soil and the tangential and normal pressures that are developed at the soiltool interface, with a good agreement between modeling and experimental results [56].
Later, several 3D finite element analyses were developed to describe the soil-tool interaction. Among these, Chi and Kushwaha [30] analyzed the 3D soil failure when subjected to the action of a narrow blade. Figure 4 shows the finite element mesh for the blade. The authors considered interface elements between the soil and the tool interface so that both drift and lift forces were not majorly affected by the adhesion and friction reacting between the soil and the blade's surface. A hyperbolic model was also used to approximate the shear stress-relative displacement between the soil and the blade due to the reduced computing time [30]. In this work [30], the authors found that the soil stress increased with the blade displacement, and the stress level was also found outside the failure surface during loading. To overcome the latter issue, the stresses were corrected in the normal direction with respect to the failure surface following the procedure proposed in [57,58]. Considering a small displacement in relation to the nodes on the interface at each increment, reaction forces are computed [30]: a force-displacement relationship is obtained and the draft is observed to increase with the displacement. The estimated draft force behavior was found to be in agreement with the experimental finding. In particular, good consistency was observed concerning the failure zone and displacement field [30]. The authors found that the initial soil failure can occur either around the tool tip alone or simultaneously around the tool tip, the tool edge and the soil surface. The failure of the soil structure thus occurs when a complete failure surface is formed from the tool tip to the soil surface [30]. Compared to the previously discussed analytical models, the numerical model developed by Chi and Kushwaha [30] presents more realistic results.
Upadhyaya et al [59] highlighted the ability of FEM to analyze complex engineering problems involving geometric and material nonlinearities (e.g. soil-machine or soil-plant interaction problems). However, many complex processes (e.g. the separation and mixing of soil layers, the appearance of cracks, the flow of soil particles, etc) occur during the soil deformation and none of these discrete processes can be adequately modeled with the finite element analysis, which appear to be more suitable for continuum studies [59]. Moreover, traditional FEM does not provide accurate results when investigating large deformation problems since it is a gridbased method. Therefore, it suffers from grid distortion leading to an inaccurate solution or failure of computation [60]. Therefore, other methods were preferred to simulate large deformation problems. Among these, the smoothed particle hydrodynamics (SPH) method was proposed [60]. This method is Lagrangian meshfree and considers particles carrying field variables and moving with the material velocity. Bui et al [60] adopted the SPH for simulating large deformation problems of non-cohesive and cohesive soils and agreed well with FEM solutions and experiments. Besides the SPH method, the total Lagrangian and the update Lagrangian approaches represent another option to simulate large deformations. However, since the calculation needs to stop if a few elements within the mesh become importantly distorted, approaches based on the ALE were proposed [61]. In this context, three methods were popular [61]: the remeshing and interpolation technique with small strain, the efficient ALE approach, and the coupled Eulerian-Lagrangian approach. In 2015, Wang et al [61] analyzed the performances shown by the three approaches reported above by considering different geotechnical problems (i.e. pipeline penetration and lateral movements in clay) and found that, despite the three methods showed different limiting aspects, their results were similar for the quasi-static penetration problems.
Recently, the material point method (MPM) [62] has been attracting attention for its ability to model very large deformation problems and for its straightforward implementation respect to other continuum meshless methods [63]. In the MPM, the information on material is considered to be in correspondence of points that are scattered through the problem domain, and it is mapped to the nodes of a finite element mesh; a calculation is performed, and the resulting nodal displacements are mapped back to the material points, and the cycle restarts again [63]. The MPM is advantageous because the background mesh can be reused in its original undistorted form in each cycle, and all the deformation is recorded in the material points so that a distorted mesh does not occur and large deformations can be accommodated [63].

DEM
Solutions to the issues encountered in continuum mechanics appear to be provided by the DEM [64].
DEM models can simulate (a) the mechanical behavior of a system composed of discrete particles with arbitrary shape and (b) the interaction between this system and a body, which can be rigid or flexible [64]. When the soil-tool interaction is considered, the soil typically represents the component of the system to be analyzed, and is modeled with discrete particles (i.e. granules) that displace independently from one another and interact only at contact points. The DEM appears to be a suitable method for modeling the nonlinear behavior of the soil and the soil-tool interaction due to its ability to consider the soil failure, deformation, and translocation [64].
The majority of DEM models are based on the approach proposed by Cundall and Strack [65], which derives from a study previously developed for studying rock mechanics problems [66]. In a granular medium, the interaction among distinct particles represents a transient problem with states of equilibrium developing whenever the internal forces are balanced [65]. The particles equilibrium forces and displacements are computed by tracing the movement of each particle. Propagations of disturbances in the medium induce particle movements. These disturbances originate at the boundaries and proceed with a speed depending on the physical properties of the medium [65]. The DEM describes this dynamic process by considering small time steps so that, during a single time step, the velocity and acceleration values can be assumed to be constant and the disturbances are transmitted from one disc to the discs in contact with it. For each time step, the resultant force acting on a disc is then computed considering only the interactions with its neighbors [65]. DEM models can properly represent many behaviors observed in granular materials (e.g. the shear strength-dilatancy, pressure, density dependency, jamming, and flow) depending on the individual particle interactions [67]. DEM simulations can be performed in 2D and 3D environments. 2D simulations are computationally cheaper because a smaller number of degrees of freedom is assigned to each particle (three in 2D and six in 3D), and the number of contacts per particle is smaller since the system is restricted to in-plane contacts [67]. Moreover, in 2D simulations, the visualization of results is easier, and the commercial 2D DEM codes are cheaper than 3D codes [67]. However, respect to 3D simulations, 2D simulations are less accurate in capturing the soil behavior since they are mostly capable of providing qualitative information. Indeed, in most cases, considering a lower number of degrees of freedom affects the mechanics making the results not easily verifiable. In DEM simulations, soil particles are represented in different sizes and shapes. Since DEM simulations are computationally expensive, scaling up the particle size is often needed, affecting the model's accuracy [64]. Moreover, simplifications are introduced in the shape used to represent the particles. Disks or ellipses and spheres or ellipsoids are used in 2D and 3D environments, respectively.
Using circular and spherical particles is often preferred due to their computational simplicity [68][69][70]. Using ellipses/ellipsoids implies the need to solve in each time increment at each contact point a nonlinear equation [67]. Often, a compromise is to adopt 'clumps' or 'clusters' of particles by bonding them together or bringing the particles together into a rigid agglomerate [67]. Figure 5 shows an example of granular assembly modeled through 3D DEM simulations using spherical particles [71].
It is important to highlight that representing particles with circular and spherical shapes leads to important rolling respect to real granular materials. This issue can be overcome by using polygons and polyhedrons. However, this implies a more complicated and expensive detection of contacts and computation of forces and torques [72]. Another solution considers aggregates or clumps of particles (disks or spheres in 2D and 3D, respectively) [72]. In this case, a large number of spheres is required, implying an important complexity in evaluating the contribution of non-convex surface of clumps. Therefore, it is often preferred to fix the particle's rotation or consider the transfer of a moment between elements in the local constitutive law. In 3D simulations, the transfer of moment can be considered, for instance, respect to the rolling inter-particle relative rotation. Indeed, at the contact point between two particles, a contact model with a rolling resistance is used to reproduce the mechanical behavior of assemblies of subrounded and subangular particles [72]: this model simulates the inter-particle normal contact response with a constant stiffness spring and the shear and moment responses with constant stiffness springs and sliders. Several studies analyzed the effect of including a rolling resistance at the contact region on the macroscopic behavior and found that the macroscopic internal friction angle increases due to the rolling elastic stiffness increase but without affecting the dilatancy importantly [72].
In DEM models, interactions among particles are governed by physical laws and are computed using contact models that consider the contact among particles to take place in correspondence with punctual regions, i.e. individual points. After detecting contact points, all the forces acting on a particle are known, and the particle position and orientation are computed using Newton's second law of motion [70]. As stated earlier, the DEM assumes that each particle transmits the contact forces only to close neighbors within the chosen time step. These calculations are repeated until the system reaches a defined balance or a target simulation run time [70]. Different contact models have been proposed to describe the mechanical relation between particles [70]: (a) the linear spring contact model is based on the work by Cundall [66] and considers only frictional linear elastic collisions, and the stiffness and damping parameters are determined for each material as a constant; (b) the Hertz-Mindlin contact model assumes the deformation at the contact point to be nonlinear elastic, and stiffness and damping coefficients are calculated using relative displacement-based equations; (c) the hysteretic spring contact model (HSCM) is based on the work by Walton and Braun [73] and considers the plastic deformation behavior in contact mechanics equations. Ucgul et al [70] concluded that the HSCM provides a more realistic force and soil movement prediction in DEM simulation rather than assuming a simpler elastic model. DEM simulations are widely used in many studies concerning soil operations.
Regarding 2D DEM models, good results are observed in the analysis by Tanaka et al [74], who modeled in 2D the interaction between the soil and a vibrating subsoil to find the best shape for the subsoiler shanks and obtained highly consistent modeling and experimental results [64,74]. In 2001, Masson and Martinez [75] used 2D DEM simulations for investigating loose and dense specimens of granular material when subjected to direct shear tests. At a macroscopic scale, the authors observed the typical features of the shear response of granular materials (i.e. a perfectly plasticity state, peak stress, and a dilatant/contraction behavior for the dense/loose sample) [75]. At a microscopic scale, the authors evaluated the particle displacements and the particle rotations, which provided a useful indicator of strain localization. The directions of the microscopic stress correlated to the strain tensors and the microstructure-induced anisotropy [75]. Shmulevitch et al [76] built a 2D model to describe the performance of several soil bulldozer blades in low-cohesive soil, achieving a good agreement with experimental results [76]. In addition, DEM models have been adopted to simulate cone penetration problems to better understand granular material behavior and its influence on cone penetration resistance [77]. 2D DEM models of cone penetration tests were developed by [78][79][80]. In these 2D models, kinematic constraints differ significantly from those in real granular materials, and, therefore, it is only possible to compare modeling results qualitatively with experimental findings [77][78][79][80].
Regarding simulations in 3D, Horner et al [81] developed a 3D DEM model to simulate soil deformation when subjected to plowing [24,81]. In 2002, Hofstetter [82] developed a 3D DEM analysis to study the interaction between a bulldozer blade and the soil by estimating the horizontal force in agreement with the experimental values and underestimating the vertical force with respect to the value obtained experimentally [24,82]. In 2006, Cui and O'Sullivan [83] analyzed the response of granular material in the direct shear apparatus at macro-scale and microscale levels by DEM simulations in 3D and compared the modeling results with experimental findings obtained through physical tests finding a good consistency. The micro-scale interactions observed in the DEM simulations represented well the interactions in the physical tests [83]. Mak et al developed a 3D DEM to simulate the soil cutting process of a simple soil-engaging tool [84]. In 2011, Arroyo et al [85] adopted 3D DEM models to build a virtual calibration chamber representing each particle with a spherical shape and preventing particle rotation for modeling the rotational resistance of nonspherical particles [77,85]. In 2012, McDowel et al [77] used DEM models to simulate in 3D the cone penetration of granular materials. Filling the whole chamber with small sized particles involves highly time-consuming simulations. The authors [77] therefore reduced the computational time by assuming (a) smaller particles next to the cone and (b) larger particles further away from the cone. Grains were organized in layers to prevent small particles from migrating into voids located between larger particles. McDowell et al [77] found that the cone penetration resistance was similar to that obtained considering a chamber filled with small particles, in agreement with the values found experimentally [77]. In 2019, Khosravi et al [86] adopted the DEM to simulate in 3D cone penetration tests and, in particular, to investigate to what extent the tip resistance and the friction sleeve shear stress values are affected by modeling parameters. The achieved modeling results were in agreement with the values found experimentally and showed that (a) parameters affecting the soil friction angle (e.g. interparticle friction coefficient, rolling resistance coefficient) have a strong influence on the tip resistance and a less important effect on the friction sleeve shear stress, and (b) variations in the probe-particle friction coefficient along the friction sleeve produce significant changes in the friction sleeve shear stress values but have a negligible effect on the tip resistance value [86].
Furthermore, DEM simulations are used to simulate how the crushing of grains affects the mechanical behavior.
Lin et al [87] used the DEM to model the effect of grain crushing. They evaluated the model's performance by simulating biaxial tests in both 2D and 3D environments proposing a good compromise between accuracy and computational cost. For developing this DEM model, the authors combined two numerical methods that are frequently used to model grain crushing: the replacement method, which uses a pre-defined failure criterion to check if particles break and that, if met, implies that the grain is deleted and substituted by smaller grains; and the agglomerate method, which considers crushable grains made of connected smaller grains and checks the grain crushing by checking the strength of the bonds [87]. Both the replacement method and agglomerate method have been widely used. However, both show some limiting aspects [87]. For instance, the replacement method cannot accurately model grains crushing due to the predefinition of fragmentation type. Meanwhile, the agglomerate method is more accurate, but it requires many particles, implying highly timeconsuming simulations.

A bioinspired perspective
Bioinspiration is an old approach that has recently become more common in geotechnical engineering [13]. Bioinspired robotics support soil investigations by developing autonomous devices capable of moving in clutter environments (e.g. [8,[88][89][90][91][92][93]). These systems are often adopted as robophysical models to analyze the complexity of body-environment interactions which otherwise would be challenging to predict and embed in simulated environments [94]. Plant roots have been extensively used as inspiration for developing penetration mechanisms with the twofold goal of proving that strategies adopted by the biological model are effective and energy advantageous for both natural and artificial systems [4,22,95].
However, understanding root-soil interactions remains challenging with (a) direct observations and analyses of the biological model, (b) physical implements, or (c) mathematical modeling. The first approach is limited by the intrinsic difficulties of working with, observing, and disambiguating the behavior of natural roots. Physical implements are often at larger scale than natural systems, possibly causing discrepancy in the performances of the natural vs the artificial models due to their different scales. Mathematical modeling is still limited by the unknown laws undergoing the root-soil interaction.
In fact, these approaches should be adopted in synergy to compensate and build on the knowledge from all three.
Here, we briefly summarize the key achievements of modeling borrowing intruders, implementing plant-inspired penetration strategies, and plantrelated investigations, coupled, when available, with experimental results.
Greacen and co-workers developed a mathematical model where the root growth was approximated by expanding spherical and cylindrical cavities [96][97][98]. The authors estimated the stresses developed around the cavities and observed higher pressures in the proximity of the penetrating root, with the highest values closer to spheres rather than cylinders [19,96,97]. However, the model did not account for several soil mechanical properties, which thus affected its validity and performance [97,99].
In 1993, Faure [100] used FEM analyses to investigate the pressures exerted by the root growth in the soil. The author considered a linearly elastic soil model, with the root tip located at a depth of 10 cm. As the root grew, a circular plastic region appeared around the root tip, which was wider in stronger soils [100]. The author also observed high soil compaction in a radial direction around the root tip, exerting a compression force with constant direction, independent of the soil properties (i.e. depending only on the biological material) [100]. However, the analysis did not involve the time domain, and, therefore, penetration and deformation rates could not be defined.
Kirby and Bengough [19] adopted FEM to imitate the growth of pea roots in sandy loam and clay loam at different densities and estimated the stresses the root-inspired penetrating body was exposed to while penetrating. Root penetration was approximated by a hollow soil cylinder sliding upward against a thin inner cylinder (i.e. the root), and minimum and maximum friction cases were analyzed at the interface between the root and soil [19]. The authors found that (a) the highest axial stress occurred in relation to the tip, which was even higher when the root was frictionless, and (b) both the axial and shear stresses were three times higher when penetrating in sandy loam [19]. A penetrating root with a smaller diameter was also subjected to higher axial stresses than those roots with a larger diameter. This indirectly suggests that root thickening may be an advantage to overcome mechanical constraints. In this work, the authors aimed to understand growth limits in roots [19]. However, in the model they did not include a differential elongation rate at the flanks of the root behind the tip, nor did they account for the steady-state Figure 6. Example of approach used to model the penetration process into soil: root penetration considered as a mechanical inclusion problem for a cylinder with finite length (i.e. the root) in a matrix with infinite length (i.e. the soil packing). Reprinted from [101], Copyright (2020), with permission from Elsevier. condition (i.e. the state in which the predicted stresses and strains vary slightly) which is reached after some penetration depth [19].
In 2018, Mishra et al [95] analyzed how the performance of a penetrometer is affected by its tip morphology. Adopting the Zea mays roots as a model, the authors developed a method to extract the root tip morphology and fabricate a bio-inspired penetrating probe to compare its performance with probes with different tip shapes. Penetration tests were performed in sandy loam soil, where the probes were pushed from the top down to a 10 mm depth at a constant speed of 10 mm min −1 , and numerical simulations based on the DEM were also performed [95]. The soil was approximated as a packing composed of homogeneously sized spheres with mechanical properties typical of sand. The probes were modeled as rigid bodies penetrating at a higher speed (10 mm s −1 ) than experimental settings in order to reduce the computational time [95]. The results of the DEM numerical simulations were in agreement with the experimental findings. In fact, probes with conic or root-like tip shapes were found to be more efficient than others in terms of penetration force and effort. However, when probes with conic and root-like shapes were compared, the root-like tip showed a better performance experimentally. However this advantage was not observed numerically due to the small geometric difference and the particle diameter considered. The authors demonstrated that tip geometry might be crucial for designing penetration systems [95].
Calusi et al [101] developed a mathematical model to investigate mechanical stresses acting on plant roots during their penetration into soil. The authors approached the penetration phenomenon as a mechanical inclusion problem for a cylinder with finite length (i.e. the root) in a matrix characterized by infinite length (i.e. the soil packing) (figure 6). The axial growth was analyzed, assuming it to occur only at the tip region of the embedded cylinder [101]. The proposed mathematical model is based on continuum mechanics and assumes both the cylinder and the matrix to be linearly elastic, homogeneous, and isotropic materials. The modeling results were compared to those obtained experimentally, analyzing the growth process of a Zea mays primary root in artificial soils with different concentrations of Phytagel and in real soils that differ in compactness. The results from simulations were consistent with the experimental ones [101]. Key findings included: (a) mechanical responses depended on the mechanical properties of the environment, (b) the scaling parameter, which is related to the input energy, increased in artificial soil with a higher concentration of Phytagel and real soil with higher compactness, and (c) the root elongation decreased in artificial soils with a high concentration of nutrients. However, in the latter case, the radial expansion observed in the model underestimated the experimental finding. The mathematical model did not discriminate between cell division and cell elongation, and did not involve the mechanical properties of the root tissue and soil medium, or the hydraulics of root growth [101].
In 2021, using 3D DEM, Chen et al [102] evaluated the self-penetration of a bioinspired probe, which adopts an 'anchor-tip' burrowing strategy, such as that adopted by earthworms or plant roots when radially expanding their bodies to increase body-soil friction and enhance penetration. The authors evaluated the effect of the probe geometry and soil depth on the ability of the bio-inspired probe to self-penetrate. Each DEM simulation consisted of three stages: first, the probe was pseudo-statically pushed into the granular soil, then an anchor was thickened, and lastly, the probe tip and the anchor were moved in opposite directions. The authors found that (a) the self-penetration ability (i.e. the ability to generate the anchoring forces needed for in-depth penetration) of the probe was influenced by the probe configuration and soil conditions, (b) these forces had more impact when generated by longer anchors and with greater expansion magnitudes and friction coefficients, and (c) the penetration resistance decreased temporarily due to the thickened anchor [102]. The model did not consider the effect of the medium density on the self-penetration ability of the probe. However, Martinez et al found that coarse-grained soils are more challenging [103].

Discussion and conclusions
This paper provides an overview of the different modeling techniques adopted to simulate the interaction between a penetrating tool and the soil to highlight their advantages and limitations.
Initially, analytical models based on Terzaghi's study [25] were widely adopted. However, several aspects limited the performance of these models: they involved a preliminary assumption regarding the soil failure pattern and simplified tool shape, and thus were unable to predict fields of stress, displacement, speed, and acceleration of the soil-tool interactions. The FEM partly overcame these limitations, providing more realistic results than those achieved with analytical models [31]. However, FEM analyses were found to be more suitable for continuum studies and, therefore, not for describing typical processes occurring during soil deformation (i.e. the separation and mixing of soil layers, the appearance of cracks, the flow of soil particles) [59].
By contrast, the DEM represents a more suitable solution to study these aspects [64].
The state of the art of the DEM numerical models used to study soil-tool interactions seems promising in the design of artificial penetrating systems. Despite early investigations dating back to the 1970s [65,66], which proved their effectiveness in soil-intruder analysis, the extensive use of DEMs was only possible with the introduction of high-speed processors (in the late 1990s), which solved complex problems in high discrete domains, including soil penetration scenarios. Today, parameters such as the resistance exerted by the soil on a penetrating body can be predicted accurately, providing an essential contribution in the design of robotic systems capable of efficiently penetrating the soil.
The further improvement of modeling approaches will have an impact on both engineering and biology, by not only driving novel design solutions for artificial systems but also helping in understanding the phenomena behind living beings and their interactions with the surrounding environment.

Data availability statement
No new data were created or analyzed in this study.