Horizontal penetration of a finite-length intruder in granular materials

In recent years, bio-inspired burrowing robots and other intruder problems in granular media have attracted significant attention. Many of these, especially related to traditional penetration problems in geotechnical engineering, cover vertical penetration. Modelling these types of problems numerically using the discrete element method (DEM) is typically done ignoring gravity by controlling the stresses in the selected representative volume. Additionally, most problems involve infinitely long intruders from a modelling point of view. However, in horizontal penetration there is enough evidence to show that intruders are affected by an uplift force that affects the penetration and needs to be considered. In this paper we use the DEM to demonstrate the impact of considering vertical uplift and gravity for a finite-length intruder penetrating in a dense granular media through particle level and macro-behaviour. Additionally, through the numerical study, important mechanisms emerge during horizontal penetration, including four different distinct stages on the surrounding soil, or the extent of disruption, that are fundamentally distorted when gravity is ignored.


List of symbols
Number of particles with only one or no contacts, respectively n p Ratio of intruder to particle in the chamber center n 1 , n 2 , n 3 , and n 4

Introduction
Many construction activities involve soil penetration with an object that dramatically displaces the soils, e.g., site investigation, pile installation or tunneling, to name but a few. All these activities cause significant disruption to the ground and can have a negative impact on the environment [7]. Such activities also typically consume excessive energy and largely rely on manual operation that could benefit from increased autonomy. In recent years, rapidly growing efforts have concentrated in overcoming the above shortcomings of traditional geotechnical engineering equipment. A major step towards this direction has been using bio-inspired burrowing strategies [15] to develop a new generation of burrowing robots aiming for engineering applications [33]. Inspired by the valve expansion/contraction digging strategy of razor clam, Winter et al. [50] developed the 'RoboClam' robot which uses a significantly reduced amount of energy compared to what would be required for pushing with brute force. It achieves this by using localized fluidization and can dig vertically downwards to a depth of 0.3 m in performed tests. An extension of the robot to burrowing in dry soil was later made [22]. Razor clams do not only burrow downwards into the soil but also upwards out of the soil [6]. Tao et al. [46,47] studied the upward burrowing strategy of razor clams and developed a soft robot that can autonomously burrow out of the soil from a certain depth. The robot achieved upward movement by cyclically inflating and deflating its body.
Robot developments have not been restricted to vertical locomotion though, and unique challenges are presented during horizontal penetration. Maladen et al. [30][31][32] developed a horizontal sand-swimming robot and the Resistive Force Theory (RFT) characterized the interaction between the robot and the ground. Whilst moving, the robot rose slowly due to the drag-induced lift force that is generated during horizontal movement in granular media. The presence of this lift force was also observed in the locomotion of a marine worm-inspired soft robot [39] and of a subterranean burrowing soft robot [35]. To overcome the lift force, strategies need to be employed such as controlling head shapes [30][31][32]. For example, Naclerio et al. [35] counteracted this uplift force by using an asymmetrical tip shape fully coordinated with a downwards air-flow that reduced soil resistance under the tip. Locomotion was achieved by a combination of forwards air-flow and body eversion using a pneumatically inflatable and growing body as initially presented by Hawkes et al. [19].
The underlying mechanism of drag-induced lift force in granular materials has been extensively explored before [14,18,41,43,56]. Gravitational stress gradients around the moving object exert a resulting upwards force that causes lift [18]. However, critically, in those experimental studies, lift force was measured on a vertically restricted intruder where the actual uplift movement was not allowed.
To gain a better understanding of this issue on finitelength burrowing robots, numerical simulations using the discrete element method (DEM) can be helpful. DEM is advantageous to deal with dynamic problems of soil-intruder interaction in granular materials, as it can give simultaneously very precise information about the macroscale behaviour whilst having access to underlying microscale mechanisms [26,37]. DEM-based models have been successfully used to simulate other analogous and more traditional penetration problems in geotechnical engineering such as CPT [1, 8,25], SPT [53][54][55], pile installation [16,28] and bio-inspired self-burrowing explorations [7,21] in granular soils. However, these studies concentrated on vertical penetration of assisted infinite-length penetrometers with no consideration of gravity.
Based on the above, the objective of this study is to investigate the behavior of a finite-length intruder during horizontal penetration into a granular material using DEM. We pay a particular emphasis on gravity which we first hypothesize is critical for accurate horizontal penetration modelling. In the following sections we describe how to build a threedimensional DEM model filled with a calibrated discrete analogue of a representative quartz sand for the penetration. The methodology employed to consider uplift displacement is introduced in detail. The effects of uplift and gravity on horizontal penetration responses are presented from both macroscopic (e.g. tip resistance, uplift magnitude) and microscopic (e.g. particle displacement, contact force, mechanical coordination number) perspectives. The results reveal some interesting features of the unique behavior of finite-length intruders.

Fontainebleau sand analogue
For the granular material, a discrete analogue mimicking the mechanical properties of Fontainebleau sand, was created using unbreakable spherical particles. Particle rotation was fully restricted by fixing all the rotational degrees of freedom of the particles. This simple approach, as a manner to mimic the effect of non-spherical particle shapes, can be traced back to [49] and was successfully employed in previous work showing penetration into granular materials [1, 5,8,54]. More refined particle-scale characteristics such as particle crushing [9], particle shape effects [42] and particle surface roughness [36,40] were not considered in this study for simplicity.
Contacts between particles are assumed to be elasto-plastic. The slip behaviour at contacts is defined by a friction coefficient µ. The simplified Hertz-Mindlin contact model is used to represent non-linear stiffness of each contact, which is controlled by the elastic properties of the material particles, e.g. shear modulus G, and Poisson's ratio v. The contact model properties (G, µ, v) (Table 1) were taken from a calibration given by Ciantia et al. [10,12] for the aforementioned sand. All the numerical tests presented in this study were performed using the DEM code PFC3D [24].
The calibrated parameters were determined via comparisons with two triaxial compression tests at low confining pressure (100 kPa) as reported by Seif El Dine et al. [44]. This low pressure is representative of shallow horizontal penetrations. The numerical tests were performed using a cubical cell of 4 mm in size containing 11,000 particles. Figure 1 shows the PSD used with diameters ranging from 0.1 to 0.4 mm, which were selected to follow the actual PSD of the material. Figure 2 shows the comparison of the DEM results versus the triaxial tests. It shows that the DEM simulation can capture the behavior of the sand in loose and dense state, including dilation effects which is known to be important for penetration problems [2].

Preliminary details of chamber construction
A three-dimensional 0.8 m (x-axis), 0.6 m (y-axis) and 0.6 m (z-axis) cuboid testing chamber was constructed using wall elements. All the chamber walls are rigid and were set to be frictionless. Geometrical model details can be found in Fig. 3 and Table 2. The choice of a cuboid chamber for horizontal penetration was inspired by similar tunneling-related explorations [29], whilst cylindrical chambers are typically used for vertical penetration tests [17,23,38]. In any case, the distance from the intruder to the walls is the critical parameter, rather than the chamber shape.
Discrete elements filling up the chamber have the same contact properties and shape as those used for the particle assembly calibration. However, to obtain a model with a  [54] computationally efficient amount of particles, their size was upscaled applying distinct scaling factors following the particle refinement method (PRM) [11,21,34,45]. If we consider half of the domain, each scaled area has a width of 1/8 T, where T is the full domain's width that has an initial value of 0.6 m (see Fig. 3). By applying this method, the particles at the center of the sample were upscaled with the smallest scaling factor of 39, leading to an intruder/particle ratio of 3.7. This ratio is a key factor in DEM penetration simulations and is larger than the ratios adopted in most of the relevant 3-dimensional DEM penetration studies as summarized in a recent study reported by Chen et al. [7]. Particles further away from the center were upscaled with gradually increasing scaling factors. The ratios between two adjacent zones (n 2 /n 1 = 1.46, n 3 /n 2 = 1.32, n 4 /n 3 = 1.30) are all smaller than 1.5, which is a proper limit ratio between adjacent zones where no obvious particle migration happened as validated by McDowell et al. [34]. The upscaling of particle sizes does not limit the precision of overall response as particle mechanical properties remain unchanged [34].
Mass scaling is not used as this would not effectively reduce computational time as demonstrated by Yousefi and Ng [52], and would also greatly affect the results since we are considering gravity. The radius expansion method (REM) was used to fill the sample. A dense sample with a target porosity of 0.38 was attained after manipulating 134,194 particles. The process of attaining the target porosity started with reducing the inter-particle friction while maintaining an isotropic compression of 5 kPa on all sample walls with the help of a servo-mechanism. After reaching equilibrium at the target porosity, inter-particle friction was reset to the calibrated value. The stress was then ramped up to a target level and the boundary conditions were adjusted for the different cases as specified in the simulation program (Table 3). In all simulations, a local damping of 0.05 [13] was employed and no viscous damping was considered.
The simulated steel intruder consists of a 108 mm long shaft and a cone tip of 60° was modelled using frictional rigid walls. A diameter of 30 mm has been selected to coincide with that of cone penetration test (CPT).

Horizontal penetration
Horizontal penetration was initiated from the central horizontal axis and pushed through the sample at a constant rate  of 50 cm/s. Butlanska et al. [4] showed that rates between 2 and 50 cm/s did not change the static penetration resistance observed in DEM modelling of CPT tests. The driving rate led to an inertial number of 0.0054 (< 0.01) indicating that quasi-static conditions could be maintained during the constant penetration [10,12,27].

Consideration of drag-induced uplift
Drag-induced effects have been routinely explored focusing on the measurement of the generated lift forces, where experiments restricted upward displacement [30,31]. This procedure is though unable to reveal the actual uplift movement that would occur during unrestricted advance of bioinspired intruders of finite length. Drag-induced uplift can make pre-defined destination points of the intruders unachievable. It is thus necessary to estimate it. In the DEM model, the intruder is represented with a macro-element using bundled rigid walls. It is anticipated that the intruder will rotate, inclining towards the surface, during penetration. However, as reported by Huang and Tao [20] in a rotational horizontal penetration study, the quantity of inclination during horizontal penetration appears to be negligible when an overburden is applied. Although the effect of rotation on the inclination is not clear, horizontal translation is believed to be the main influencing factor to inclination of the intruder as demonstrated as well in the cases with no overburden in Huang and Tao [20]. Therefore, it is assumed in our study that during horizontal penetration, the intruder does not incline but only translate vertically in the x-z plane. Its vertical displacement is enabled by a user-defined algorithm where the force components acting vertically on the intruder are computed to update its velocity and position at each timestep. The assigned vertical velocity has the following equation: where ̇z t+Δt and ̇z t are the vertical velocities at time (t + ∆t) and t, respectively, ∆t is the time step and m r is an assigned intruder mass. F t tot,z is the resultant vertical force acting on the intruder: where F t tip,z , F t shaft,z , and F t bottom,z are the vertical reaction forces at the intruder tip, along the shaft and at the bottom, respectively, i, j, and k represents a given particle that is in contact with the tip, the shaft and the bottom, respectively, and g is the gravitational acceleration (Fig. 4).
The intruder mass m r is determined from the values of length l p and material density ρ r (Table 2) using steel as reference material. The gravitational acceleration g is set as − 9.8 m/s 2 . Table 3 summarizes the three simulated cases for the horizontal penetration study. In Test T1, drag-induced uplift is restricted following typical testing configurations of horizontal penetration [18,56]. In T1, a K 0 boundary condition is applied with an overburden pressure P 0 = 20 kPa acting on the top boundary with the displacement of other boundaries being fixed. In addition, gravity field is activated after applying this K 0 boundary condition. T2 is the same as T1 but allowing uplift. Comparisons between these two tests are used to examine the effect of uplift on horizontal penetration responses. The effect of gravity is also checked in T3 by setting gravity to 0 from test T2. During the intruder penetration, the specified stress levels in Table 3 were maintained constant using a servo-mechanism.

Simulation configurations
The following test results demonstrate the effects of the examined factors on penetration responses from the perspectives of both macroscale (e.g. tip resistance, uplift displacement) and particle-scale (e.g. contact force, particle displacement).

Horizontal tip resistance
Raw penetrograms show horizontal tip resistance against penetration distance, where different levels of fluctuations can be observed as in Fig. 5. The fluctuations are mainly affected by the ratio of cone to mean particle size n p on the simulations as for low n p the plot becomes very jagged [1]. The fluctuations can be filtered out though to obtain a smooth curve by fitting the raw penetration curves with the following equation: where q t,x (MPa) is the static horizontal tip resistance, h x is the horizontal penetration distance, A and B are fitting parameters. Parameter A (MPa) gives the value of tip resistance at steady state. The raw and filtered resistance values measured in T2 are illustrated against horizontal advance h x in Fig. 5. The horizontal advance was terminated at 70 cm. The parameters of the adjusted penetration curve for all the tests, A, B and goodness of fit, R 2 , are summarized in Table 4.
The filtered penetration curves for all the tests are presented in Fig. 6. The tip resistance values also confirm that particle breakage does not occur for the silica sand particles [51]. The average tip resistances with and without vertical displacement are 6.09 MPa (T2) and 6.23 MPa (T1), respectively, indicating a 2.3% tip resistance reduction if vertical displacement is allowed. This reduction is attributed to loosening effects ahead of the tip induced by vertical displacement. The penetration curve for T3 reveals the effect of gravity compared to T2. Ignoring gravity leads to smaller stress levels in the chamber as shown in Table 3 and thus results in an obvious reduction (22%) in tip resistance (T3).    [18], while the small tip area leads to higher fluctuations. This is clearly shown in Fig. 8 by their probability distributions using all the values from the beginning of the penetration until a penetration of 70 cm is achieved. The forces (F shaft,z , F tip,z and F tot,z ) for the cases with uplift (T2 & T3) and their cumulative distribution functions (CDFs) are plotted in Fig. 9. These cases reveal several common features: (i) the tip force and shaft force present nearly identical magnitudes but opposite directions through the whole horizontal advance, resulting in the total vertical force constantly oscillating around 0. These force trends exhibit a dynamic process of equilibrating vertical forces acting on the intruder by adjusting the intruder's vertical position; (ii) the component forces present both positive and negative values, dominated by unbalanced strong force chains around the intruder. As shown in the CDFs, the shaft force and tip force distributions are in parallel, with a slight gap in between for the case with no gravity (T3). The mean values of the recorded forces for T3 are all approximately 0 (Table 5). However, in T2 the gap becomes more visible and thus reveals a clear effect of the gravity-induced stress gradient. In this case, the mean value of tip force is positive (0.16 kN), while the mean shaft force is negative (− 0.16 kN). Combining with feature (i), it is inferred that generally under gravity, tip force tries to push the intruder upwards and but the movement is quickly suppressed by a reaction shaft force. However, this rebalancing is not instant and the net vertical forces are dominantly positive which result in the characteristic uplift (to be presented in Sect. 3.2 'Displacement-related observations').

Particle-scale contact forces
The contact force network (Fig. 10, captured at h x = 50 cm) provides a microscale perspective on intruder-soil interaction. Strong force chains carry a large proportion of the stress in the material, while weak force chains are able to show particle rearrangement resulted from the disturbance of objects. Figure 10 illustrates 3D contact force vectors in the examined tests represented in their planar projection along a 2 cm thick vertical section located centrally in the chamber. Forces exceeding the whole ensemble average (λ) but smaller than λ + 5σ are plotted in dark grey while they are in black if CF > λ + 5σ where σ is the standard deviation. The forces smaller than the average force are plotted in orange. This color provides better visualization of the weak contact force network behind the tail. Forces below 10% of the average value are not represented. The statistics (λ and σ) obtained from test T1 are used as references for other plots.
The magnitudes of contact forces at the tip decrease with the decrease of measured tip resistance, e.g. T3 with the smallest resistance values (Fig. 6) presents thinner and looser contact forces. The spatial distribution of contact forces is also of great interest. It reveals several significant common features: (i) the strong force network clearly concentrates on the intruder tip, (ii) the region of the shaft shows relatively small forces, (iii) the weak force network concentrates behind the intruder, and (iv) a large empty area appears right behind the tail. The first two features are consistent with observations seen during vertical penetration [3], while the last two are unique observations for short intruders driven horizontally.
To examine the effect of gravity on contact force distribution, we compare the two tests with (T2, Fig. 10b) and without (T3, Fig. 10c) gravity. The extent of the weak force chains behind the tail reveals significant differences. T3 in Fig. 10c shows a symmetric contact force network with the empty area behind the tail prolonged to the entry point of the chamber (left boundary). In contrast, the test with gravity in T2 (Fig. 10b) shows clear asymmetric weak contact force Fig. 9 Vertical tip, shaft and total force evolution (a) and cumulative distribution function (CDF) of the forces (b) for T2: VD, g ≠ 0; Vertical tip, shaft and total force evolution (c) and cumulative distribution function (CDF) of the forces (d) for T3: VD, g = 0 chains with respect to the x-axis. The forces towards the left boundary have gradually recovered as particles behind the tail fall rapidly to the empty space under gravity, resulting in a faster soil densification process than those without gravity. Figure 11 illustrates the drag-induced uplift displacement of the intruder h z recorded in T2 and T3 versus horizontal The tests experience nearly identical rapid uplift at the initial stage of penetration, because initially there is no shaft forces to counteract the uplift forces from the tip. This is slowed down as the shaft forces become larger due to larger penetration after approximately h x /l p ~ 0.1 (Fig. 11b). At a penetration of h x /l p ~ 0.4, both lines start to separate due to the different shaft forces acting on both cases as shown in Fig. 9. Then, in T2 the intruder moves vertically at an approximately constant velocity v z of 0.6 cm/s. The velocity was calculated using the vertical movement at final penetration (70 cm, h x /l p = 5.22) h z2 , at full penetration point h z1 and the elapsed time. The values of h z2 and h z1 are given in Table 4. At the final penetration point in the figure, the intruder has been lifted up by 1.33 cm in T2, corresponding to roughly half of the intruder's diameter. In the case with no gravity (T3) the uplift behavior is much less notable. This confirms that gravity is an essential factor to correctly predict the uplift, which coincides with Guillard et al. [18].

Particle displacements
The advance of the intruder unavoidably triggers rearrangement of particles in its vicinity. Different from vertically driven long intruders, the horizontally driven short intruder in this study presents several unique particle rearrangement characteristics associated with gravity-related mechanisms, e.g. lift effects at the tip and shaft and soil collapse behind the tail after the intruder passes through. These phenomena are visualized by illustrating 3D contours of particle displacement magnitudes accumulated from the initiation of penetration in a vertical plane (Fig. 12). Particles with displacement magnitudes less than 0.003 m are not displayed for a better visualization of large displacement zones. To better reflect the status of particle movement, scaled particle arrow velocity field (> 0.01 m/s) is also displayed in conjunction with displacement field.
In the displacement contours of all the cases, there appears to be a small empty area right behind the tail after the intruder passes through, which coincides with the observations in Fig. 10. When gravity is deactivated, the empty area is relatively large and the particle displacements are symmetrical around the horizontal center line (Fig. 12c). Particle densification driven by the confining pressure acting on chamber boundaries takes place slowly to fill the empty space in the case of no gravity. In contrast, if gravity is present, the empty area is rapidly filled with particles that are collapsing downwards as confirmed by the considerable particle arrow velocities behind the tail (Fig. 12a, b). Right behind the tail, the empty areas in both tests become much smaller but still visible, which may be attributed to the rapid horizontal advance and arching effects (Fig. 13) that delay the immediate collapse of particles. The particles right below the penetration path are compressed by the intruder.
The particles with gravity ( Fig. 12a, b) displace asymmetrically around the center line with more displacements in the upper side, particularly when vertical displacement is allowed (T2). This asymmetry is caused by distinct mechanisms associated with relative position of particles to the intruder. For instance, for particles behind the intruder tail, asymmetry is mainly attributed to the soil collapse behavior due to gravity. Soil collapse leads to the largest particle displacement in the specific area starting from the left boundary and approaching the intruder tail. However, around the shaft and tip, asymmetry occurs under the influence of gravitational stress gradient and is enhanced by the intruder uplift behavior. The distance of affected particles around the shaft and tip is up to 4 times the intruder diameter if gravity is considered. Conversely, this range reaches only ~ 3 times the intruder diameter in the case of no gravity. No localized surface deformation was supposed to be observed from the displacement contours, since most of the particles close to the boundaries were not displaced by the horizontal penetration behaviour. This confirms the sufficiency of the distance from top surface to intruder to eliminate boundary effects and the appropriateness of using rigid wall as top boundary.

Mechanical coordination number
In granular materials, the average coordination number is usually defined as Z = 2C/N, where C is the number of contacts and N is the number of particles. However, in numerical simulations, there are possibly some particles with zero contacts and some particles that only have a contact with one neighboring particle. These particles do not make contributions to the stable state of stress. Hence Thornton [48] defined a mechanical average coordination number: where N 1 and N 0 are the number of particles with only one or no contacts respectively.
The particles contained within the upper and lower zones of two cylindrical measurement areas (Fig. 14) with 5 cm radius and 60% and 30% of the sample length, respectively, are used to measure the evolution of Z m . The values of both upper and lower parts (Z m (up) and Z m (low)) from both measurement areas are shown in Fig. 15 for T2. The evolution follows consistent trends, which can be classified into four phases (I, II, III and IV) depending on the intruder's penetration. The features of each phase are described below: (I) At the initial stage of penetration, Z m (up) and Z m (low) both fall abruptly because of the effect of dilatancy in dense granular samples [48, 53,55]. Fig. 12 Particle displacements (70% transparency) and arrow velocities (scaled by magnitude) in a vertical slice of 2 cm thickness at h x = 50 cm. Particles with displacement magnitude less than 0.003 m and velocity magnitude less than 0.01 m/s are not illustrated Fig. 13 Contact force network for particles within a 2 cm thick slice at x = 30 cm. The model state corresponds to a state where the intruder has exited the chamber boundary. Force legends are the same as Fig. 10 A slight gap is observed between the two values induced by vertical stress gradient and uplift. A greater value of Z m (low) indicates a denser particle arrangement in the lower part. As the intruder advances, the abrupt falls of Z m gradually flatten before full penetration at h x /l p = 1. (II) After full penetration, the separation between upper and lower parts becomes clearly visible. Z m (up) keeps decreasing during the entire phase. Conversely, Z m (low) starts to flatten quickly, even showing a slight increase (Fig. 15a). The continuous decrease of Z m (up) is caused by not only the loosening effect of intruder penetration, but also the fact that soil behind the tail starts to collapse, making Z m (up) smaller and thus the soil looser. In contrast, the compression effect of the intruder (Fig. 12) on soil particles in front of and below the intruder dominates the evolution of Z m (low) in this phase. (III) After reaching a penetration equal to the measurement area, Z m (up) begins to increase due to the initiation of densification process in the entire area under gravity and the applied boundary stress. Z m (low) keep recovering and evolve in parallel with Z m (up). In this phase, gravity is acting to reach a new equilibrium in the measured area after the influence of the intruder. (IV) In phase IV, instead of recovering to the initial particle contact arrangement, both curves quickly evolve into plateaus after about 1.5l p ahead of the measurement areas maintaining a similar gap in between as in phase III. This indicates the loosening effect of finite intruders in granular materials which, for this dense particle arrangement, is larger for the upper part of the soil compared to the lower part. Figure 16 shows the Z m evolution patterns in the upper and lower parts of all the examined tests. T1 presents similar evolution trends to T2, characterizing the general evolution of Z m in the samples with gravity. From the comparisons between T1 and T2, we find that the main effect of uplift lies in phase III and IV of Z m (up), which shows a faster recovery and larger residual in T2. The T3 test with zero gravity shows continuous reduction throughout the whole penetration process within the chamber (phase I, II and III) due to continuous creation of the large cavities behind the intruder and slow densification process taking placing behind the tail. In Phase III it begins recovering but remains much looser than the other tests, showing that the disturbance is much greater. This again shows another important factor supporting the inclusion of gravity for modelling horizontal penetration. Figure 6 has shown the differences of horizontal resistances observed between the cases with and without gravity. It is interesting to further explore the underlying reasons behind the phenomenon from the angle of stress levels at the penetration level. Figure 17 shows the vertical stress values at the top, middle and bottom levels for both cases with and without gravity. Although the stress levels at the top boundary are the same, the stress levels at the middle and bottom with gravity are greater than those without gravity. Particularly, the vertical stress level at the penetration level (red dash line in Fig. 17) seems to be a dominating factor in horizontal resistance. To verify this point, we run a model T4 which is without gravity and confined by 24.5 kPa on the top boundary. This stress level is identical to that at the penetrometer level in the samples T1 and T2 with gravity. Figure 18 shows that T1 and T4 produce very similar horizontal resistance results, validating that stress level at the penetration dominates soil resistance.

Influence of other parameters on horizontal penetration
In this study, we have deeply investigated two important factors on horizontal penetration behaviours: gravity and vertical displacement. Effects of other parameters such as intruder length, gravity level, and overburden stress level are believed to also be important and should be covered in future studies. Here to stimulate that, we anticipate the possible effects of these parameters for future validation: • Intruder length for longer intruders, there will be more forces generated at the shaft to resist the intruder's uplift  Table 3 so uplift will be more limited. Equally, this is likely to reduce any rotational effects of the intruder. • Gravity level with the increase of gravity level, the uplift will become more significant due to the increase of vertical stress gradient along the intruder. A lower gravity should result in a reduction of uplift and penetration resistance. • Stress level with the increase of stress level applied on the top surface, penetration tip resistance will increase, while uplift and intruder rotation will decrease.

Conclusions
In this study, we have presented a series of DEM-based models simulating horizontal penetration of a small-scale intruder and associated uplift behaviors in a dense sand sample. Drag-induced uplift has been considered using a userdefined algorithm that updates intruder velocity accordingly from computed lift forces acting on it. Macroscale results have been presented in appropriate conjunction with microscale observations to show the influence of uplift and gravity during horizontal penetration process. The correct modelling of horizontal penetration in a dense granular material using DEM must include gravity because otherwise: • The well-known uplift observed during horizontal penetration is not captured because the stress gradient is less pronounced around the intruder. • The granular material disruption is not correctly captured. In the absence of gravity, the material disturbance is symmetrical around the intruder which is not realistic. On the one hand, the disrupted area of soil above the intruder is larger when considering gravity but on the other hand, the loosening of both above and below the intruder, observed by the coordination number, is larger in the case of no-gravity modelling. Both aspects of behavior are critical for the future design of actuating mechanisms in robot for navigation in granular materials and its interaction with the surroundings. • Critically, the predicted penetration resistance forces are 22% lower in the case of no-gravity being considered.
Additionally, outside the issue of gravity and considering the case with gravity only, we have identified four distinct zones in respect to the granular material disturbance in relation to intruder penetration. The particle disturbance of an area reached a constant value once the intruder tip was approximately 1.5l p ahead of the measurement area.
To the best of the authors' knowledge, these are the first reported three-dimensional DEM-based simulations of finite-length intruders in granular materials taking uplift behavior and gravity into account. All these considerations potentially support the development and correct DEM modelling of new robotic devices. To provide more robust support, more penetration scenarios allowing movement in more degree of freedom such as tilting are to be carried out in future studies.
Acknowledgements The authors are grateful for the PFC software support from Dr. Marcos Arroyo at Universitat Politècnica de Catalunya.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.