Inclined Block Impacts on Granular Strata: Coupled DEM-FDM Numerical Investigation and Rheological Modelling

In mountain regions, rockfall is a very widespread natural hazard, difficult to be numerically simulated, from the initial detachment of the block to its impact on natural/artificial strata. When a rock block impacts on a granular stratum, the interaction is dominated by inertial effects and is severely affected by block mass, shape, velocity magnitude/direction, and soil geomechanical properties. The kinematic description of the response of the block during the impact requires complex and time-consuming numerical modeling approaches. Moreover, these numerical approaches cannot be easily coupled with classical mass-lumped rockfall analysis tools employed to simulate blocks trajectory. To overcome this problem, a dynamic upscaled visco-plastic rheological model (BIMPAM), suitable for simulating inclined impacts on generally inclined deformable strata, has been conceived in the past by one of the authors. Up to now, due to the absence of suitable numerical and experimental data, BIMPAM has been only validated on experimental large-scale vertical impacts on horizontal strata. In this paper, the impact problem is analyzed, by performing numerical analyses by means of a coupled Discrete Element and Finite Difference numerical model (DEM-FDM), calibrated on available experimental data, to provide (i) a novel numerical dataset highlighting the crucial role played by the impacting velocity inclination on the dynamic interaction between soil stratum and penetrating rigid object and (ii) a detailed micromechanical description identifying the processes responsible for the complex macroscopic behavior. The obtained numerical results have allowed to both validate BIMPAM model, even in case of inclined velocities, and provide a micromechanical interpretation of its constitutive assumptions. The obtained numerical DEM-FDM results have allowed the authors to discuss how the initial velocity inclination affects impact kinematics. An upscaled constitutive relationship (BIMPAM model), previously conceived, has been validated for any impact inclination using the novel dataset obtained in the current work. The problem of DEM-FDM data reliability has been critically tackled proposing an original procedure for the generation/calibration of the coupled DEM-FDM model. The role of grain damaging, when the impact energy is sufficiently large, has been highlighted, from the direct comparison with large scale experimental test results. The obtained numerical DEM-FDM results have allowed the authors to discuss how the initial velocity inclination affects impact kinematics. An upscaled constitutive relationship (BIMPAM model), previously conceived, has been validated for any impact inclination using the novel dataset obtained in the current work. The problem of DEM-FDM data reliability has been critically tackled proposing an original procedure for the generation/calibration of the coupled DEM-FDM model. The role of grain damaging, when the impact energy is sufficiently large, has been highlighted, from the direct comparison with large scale experimental test results.


Introduction
In alpine regions, rockfall is one of the most common natural hazard and causes every year damages to structures/ infrastructures. When rock blocks impact on soil strata, the energy mostly dissipates in the soil. Due to the phenomenon dynamicity, when impact takes place, shockwaves propagate in the medium. The impact phenomenon is governed by (i) impact velocity magnitude/inclination (ii) impacting block characteristics (shape, dimension, mass), (iii) stratum inclination and (iv) in situ granular material relative density (Labiouse and Heidenreich 2009).
In practical applications, for the prediction of falling block trajectories, lumped mass methods are commonly employed (Ritchie 1963). Numerous are the commercial codes based on these methods, in which the impact is described by employing restitution coefficients, defined as the ratio of the final to initial block velocity. Due to the large number of factors influencing the impact process, restitution coefficients are usually empirically calibrated on experimental results (Asteriou et al. 2012).
The study of impacts is also fundamental in the design of sheltering structures, such as embankments or sheltering galleries (Peila et al. 2007;Ronco et al. 2009), commonly tackled by employing pseudo-static approaches and/ or empirical data. In case of sheltering galleries, di Prisco 2007, andCalvetti 2011, have proposed a displacement-based design approach, employing a rheological constitutive relationship, capable of describing the evolution with time of the impacting force (the "BIMPAM" model, di Prisco and Vecchiotti 2006) and on the simulation of the propagation of the shockwave within the granular stratum.
The main objective of this paper consists in describing the role of the block velocity direction in affecting the nature of the impact. This goal is achieved by using a coupled Discrete Element-Finite Difference Method (DEM-FDM) numerical code. The DEM-FDM numerical results are also used to confirm the capability of the BIMPAM model of simulating the system response, even in case of inclined impacts.
From a numerical point of view, an original procedure for the generation of the coupled DEM-FDM model is proposed by the authors. Moreover, the authors have described in detail the calibration procedure of both DEM and FDM model parameters on experimental data.
The paper is structured as follows: in Sect. 2, a stateof-the-art review is proposed and BIMPAM model is briefly summarized (Sect. 2.3); in Sect. 3, the strategies proposed by the authors to calibrate both DEM-FDM and BIMPAM models are described. In this calibration phase, an impact reference case (vertical velocity and horizontal stratum) is employed. A short discussion about the role of granular material uniformity coefficient and grain crushing is reported in Appendices. In Sect. 4, with reference to vertical impacts on horizontal strata, DEM-FDM, BIMPAM and experimental results are compared to validate and discuss the limitations of the two models in predicting impact kinematics (10 experiments characterized by different initial kinetic energies are simulated). In Sect. 5, DEM results are critically discussed from a micromechanical point of view, in particular the temporal evolution of the failure mechanism/block penetration is considered (Sects. 5.1 and 5.2). Finally, in Sect. 6, the role of velocity direction is discussed, by illustrating the results of an additional DEM-FDM numerical campaign. The obtained numerical results are employed to discuss BIMPAM predictions' reliability.

State of the Art
As was previously mentioned, impacts have been studied in the past by performing both experimental tests (Sect. 2.1) and numerical analyses (Sect. 2.2). The dynamic interaction between the block and the granular stratum has been also described by using upscaled lumped models (Sect. 2.3). Although numerous are the papers in the literature on this subject, all the authors, with the exception of Heidenreich and Labiouse 2004, disregard the role of velocity inclination (Sect. 2.4).

Experimental Tests
In the literature are available experimental data concerning both real and small-scale impact tests. Due to the difficulties in imposing an inclined impact velocity, the large-scale experimental test results only concern vertical impacts on horizontal strata (Labiouse et al. 1994;Pichler et al. 2005;Calvetti et al. 2005;Calvetti and di Prisco 2012;Labiouse and Heidenreich 2009). The previously mentioned experimental data obtained by Heideinreich and Labiouse 2004, concerning inclined impacts on granular inclined strata, were obtained in the laboratory by using a small-scale apparatus. Small scale tests are very useful to qualitatively describe the mechanical process but the data are affected by scaling effects. In addition, experimental data concern block kinematic variables, without providing information about micromechanics.

Numerical Analyses
The numerical simulation of impact phenomena requires numerical codes capable of dealing with contact problems involving large displacements and high strain rates. In recent years, large displacements continuum mechanicsbased codes, as the Particle Finite Element Method (Idelsohn et al. 2004) andthe Material Point Method (Suslky et al. 1995), have become popular in the geotechnical scientific community. Unfortunately, these methods have not been applied to study impact problems due to the lack of the implementation of constitutive models suitable for taking the high dynamicity of the phenomenon into account.
In contrast, the discontinuum-mechanics based Discrete Element Method (DEM) (Cundall and Strack 1979), seems to be more suitable to reproduce the impact phenomenon, since the approach is conceived under both large displacements and dynamic conditions. Calvetti et al. 2005 performed 3D DEM tests to reproduce large-scale experimental data concerning vertical impacts on horizontal soil cushions. In particular, they showed that, if the micromechanical parameters are properly calibrated, DEM codes are capable of predicting the evolution of the acceleration of the impacting block, as well as the stress distribution acting on the bottom plate, but not the penetration of the block. Bourrier et al. 2008 andZhang et al. 2017 by performing 2D and 3D DEM analyses, respectively, described the evolution of force chains in the soil stratum induced by the impact of a spherical boulder. In particular, they described the process of wave propagation in the soil cushion induced by the impact, causing the breakage of the initial force chain network. As far as block penetration and peak acceleration are concerned, in Zhang et al. 2017 is evident that DEM numerical simulations tend to overestimate Pichler et al. 2005 measurements. More recently, the influence of the impacting block shape has been also investigated by Shen et al. 2019and Shen et al. 2020whereas Shen et al. 2021 studied the influence of the particle size of the soil.
In the literature, the Discrete Element Method has been employed as a tool for studying impacts, but, in general, what is missing is a discussion on (i) how to properly calibrate the model in order to obtain reliable numerical results and (ii) how to avoid wave reflections associated with the presence of boundaries. These two items are tackled by the authors in Sects. 3 and 4.

Upscaled Modeling
An alternative approach to simulate block-granular soil stratum interaction, assuming the block mass to be concentrated in the center of gravity and both kinematic and static variables to be described with respect to its position, has been proposed by di Prisco and Vecchiotti (2006). This consists in an upscaled constitutive law putting in relation both generalized kinematic and loading variables, conceived for reproducing not only vertical impacts of spherical blocks on horizontal homogeneous strata, but also inclined impacts along inclined slopes. This model is inspired to the macroelement theory (Nova and Montrasio 1991), developed to calculate settlements of shallow foundations under complex loading. Recently, to extend its applicability, the model has been modified by introducing both the boulder rotational degree of freedom and the toppling mechanism, dominant when the shape of the block is not ideally spherical (Dattola et al. 2021). BIMPAM is a simple tool to be used in the engineering practice to calculate both impact forces and block penetration.
In this paper, the modified version of Dattola et al. 2021 is employed. For the sake of brevity, the governing equations will be not hereafter reported, but a synthetic description of the model characteristics and constitutive parameters is provided. BIMPAM uses as input data: initial impact velocity ( V 0 ), initial velocity inclination with respect to the stratum normal vector ( 0 ), slope inclination, boulder radius, aspect ratio and mass. It provides as output: block trajectory, translational and rotational block velocities and accelerations.
In BIMPAM, constitutive equations are formulated in terms of generalized displacements/rotations and generalized forces/moments. In the one-dimensional case, the model may be schematized as in Fig. 1, where the viscoplastic slider simulates the development of irreversible strains in the near field, whereas viscous damper and elastic spring the mechanical response in the far-field, where the dissipation is due to elastic wave traveling in an infinitely extended spatial domain.
The constitutive relationship of the elastic spring is assumed to be non-linear, taking into consideration both geometrical evolution of the block footprint, increasing during penetration, and increase in the soil confining pressure.
The visco-plastic slider constitutive relationship is based on the plastic model conceived by Nova and Montrasio 1991 with reference to rigid footings subject to inclined and   eccentric loads, assuming the development of a generalized failure mechanism under the penetrating block. In case of inclined impacts, when superficial sliding may occur, the penetration failure mechanism is inhibited and the viscoplastic element is substituted by a plastic slider obeying to a Mohr-Coulomb perfectly plastic law. Unfortunately, at the moment, due to the lack of a suitable experimental and numerical dataset, BIMPAM has been validated only in case of vertical impacts of spherical blocks on horizontal strata (Table 1).

Discussion
In Table 1, the studies available in the literature are classified according to (i) the approach employed, (ii) input and (iii) output data. As far as input data are concerned, in Table 1 the studies, focusing on stratum conditions (in terms of inclination and material mechanical properties) and block initial conditions (in terms of impact energy and block velocity inclination), are identified. As far as output data are concerned, in Table 1, the studies are classified according to the block kinematic description provided and to the availability of micromechanical information.
From the analysis of the data reported in Table 1, it is quite evident that, up to now, in the literature, much effort has been dedicated to study the influence of impact energy on the dynamic block-stratum interaction, whereas data concerning the role of impact velocity inclination, and granular stratum properties are missing.

Strategies of Models Generation/ Calibration
As was previously mentioned, to obtain reliable results, fundamental is the use of robust generation and calibration strategies. The generation procedure of the DEM-FDM model is discussed in Sect.

Generation and Test Description
The numerical analyses have been performed by combining Discrete Element Method (PFC3D, Itasca 2018) and Finite Difference Method (FLAC3D, Itasca 2019) numerical codes. The coupling strategy, proposed by Itasca (2018) with reference to quasi-static penetrations, has been used in this work to allow the dissipation of the compressive wave propagating from the impact point due to the dynamicity of the phenomenon (Calvetti and di Prisco 2007).
The impacting block is modeled as a rigid sphere of radius R and mass m . The free fall phase is not simulated. The impacting block is generated on the top of the granular stratum ( Fig. 2a) with an initial, inclined (Fig. 2b), translational velocity V 0 ( V 0 = √ 2gH , being H the falling height and g the gravity acceleration). The measurement circles, introduced for micromechanical analyses (Sect. 5), are reported in Fig. 2c.
The soil stratum close to the impact point (near field), where irreversibilities develop and large deformations are expected, is a cubic sub-domain, modeled as an assembly of polydisperse spherical particles (Shen et al. 2019;Bourrier et al. 2008;Calvetti et al. 2005). The employed grain size distribution is given in Fig. 3. A linear force-displacement contact law, set in series with a slider obeying Coulomb's friction law along tangential direction, has been used. Alternatively, to what suggested by other authors introducing a rolling resistance (Marchelli et al. 2020;Zhang et al. 2017;Wensrich and Katterfeld. 2012), in this paper, to roughly mimic the mechanical response of a real granular material with non-spherical grains, particle rotations have been inhibited. This is a simplified approach, successfully employed in previous works (Calvetti 2008;Ciantia et al. 2019Ciantia et al. , 2016Redaelli et al 2021;Calvetti et al. 2019), to simulate the mechanical behavior of granular materials. Particle crushing effects and numerical (local or viscous) damping are not accounted for (grain crushing is discussed in Appendix B), and friction is considered as the unique source of energy dissipation in the DEM discretized subdomain.
Energy radiation is allowed by the presence of the FDM external sub-domain (Fig. 2). Since significantly small displacements are expected in the far field, an isotropic elastic constitutive law is there assumed. To allow elastic wave absorption, independent dashpots along normal and tangential directions (Lysmer and Kuhlemeyer 1969) are introduced on the external FDM subdomain boundaries (Fig. 2b). The FDM sub-domain is discretized by means of a structured finite difference grid.
The generation procedure requires the following five stages: i) The FDM region is generated.
ii) The DEM subdomain is filled with particles by employing pluviation method. During gravity deposition, a very low initial interparticle friction coefficient 0 is assigned to all particles. Its value is chosen to achieve a sufficiently dense stratum. During this phase, the compliance of the homogenous material in the FDM sub-domain is assumed to be null. At the end of this stage, the authors evaluate in the DEM subdomain the value of both earth trust coefficient K 0 and stratum density . iii) The final elastic properties are assigned to the material in the FDM sub-domain. The geostatic state of stress is there initialized by using a K 0 procedure, while at the boundaries of the FDM region displacements are imposed to be null. This phase stops when inertia forces become negligible with respect to gravitational forces. At the end of this phase, equilibrium is reached even at the interface between the two spatial sub-domains. iv) Dashpots are introduced on outer lateral and bottom surfaces of the FDM subdomain and the interparticle friction coefficient is increased up to the prescribed value. The displacements at the edges of the boundaries are still maintained constrained to avoid the rigid translation of the whole domain. v) The spherical boulder is generated in contact with the granular cushion and its initial velocity is assigned.
To avoid any dependence of numerical results on model dimensions and mesh size, the authors have performed a sensitivity analysis. The size of the cubic DEM sub-domain has been chosen equal to b = 4m to obtain the plastified zones under the penetrating block to be fully contained. The FDM subdomain dimensions are fixed to optimize the dampers position (Kellezi 1998). The authors have chosen L = 6m (Fig. 2b). As far as the mesh size is concerned, a good compromise between solution convergence and computational costs was individuated in l = 0.5m.

Calibration
Model parameters have been calibrated on one reference experimental test obtained in Listolade (Calvetti and di Prisco 2012), involving a spherical block (radius R = 45cm , mass m = 850kg ), falling from height H = 4.9m , corresponding to an initial velocity V 0 = 9.8m∕s , impacting on a real rockfall sheltering gallery protected by a dense 2m thick granular stratum.
During the test, the evolution with time t of the block vertical acceleration a z (dashed line Fig. 3a) was measured by means of an accelerometer embedded in the sphere. Both block vertical velocity and displacement were derived by integrating a z (dashed line Fig. 3b and c). The granular stratum was placed on a reinforced concrete slab, on which different load cells were placed to measure the local stress increments. From the test results analysis, the mean compression wave velocity in the soil c = 250 − 300m∕s was measured. The particle size distribution of the soil cushion was unknown. The soil stratum thickness was sufficient to neglect the effect of boundary wave reflection on the previously mentioned experimental curves.
The mechanical behavior of the material in the FDM subdomain is a function of ( Table 2): (i) material density and at-rest earth thrust coefficient K 0 (imposed to coincide, for the sake of homogeneity, with the values obtained in the DEM subdomain in the second phase of the model generation) and (ii) elastic parameters ( E, ). The Poisson ratio is assumed equal to 0.25 , whereas the elastic modulus has been calibrated on the experimentally measured compression wave velocity c: In the light of the results of an ad hoc numerical campaign, the authors have derived that the dynamic mechanical response of the DEM stratum is instead mainly governed by (i) micromechanical parameters ( K n , k s , , normal, tangential stiffness and friction coefficient, respectively) and (ii) initial void ratio e 0 (Table 3). Indeed, the numerical results are not affected by the PSD choice, if and only if Coefficient of Uniformity ( C u = D 60 ∕D 10 , Fig. 4) is between 1.25 and f R∕D 50 and D 50 < 0.2R , being D 50 defined in Fig. 4 (Appendix A). The density of the particles has been imposed to be equal to s = 2690kg∕m 3 , whereas the chosen particle size distribution (PSD) to be described by the function of Fig. 4, satisfying the constraint mentioned here above.
To assign e 0 , the authors employed the definition of the relative density ( Dr 0 = e max − e 0 ∕(e max − e min ) ) and, therefore, had to perform a series of numerical tests, simulating the pluviation procedure, by changing falling height  and 0 to assess e max and e min . Since the soil stratum experimentally tested (Calvetti and di Prisco 2012) was very dense, the authors have chosen e 0 to obtain Dr 0 ≅ 90% ( e 0 = 0.538).
As far as micromechanical parameters are concerned, it is worth mentioning that damping is absent and k s ∕K n is imposed, as was suggested by Calvetti and Emerlault 1999 to obtain realistic values for the material Poisson ratio, equal to 0.25.
The two remaining micromechanical parameters ( K n and ) are calibrated (Table 3) from experimental data: A) K n ∕D 50 is related to the propagation of the compression wave velocity. To obtain the value experimentally measured, K n ∕D 50 was fixed equal to 90MPa. B) The value affects the displacement versus time curve (Fig. 5). To satisfactorily fit the experimental results, the authors have chosen = 0.35.
Once calibrated, the numerical model captures satisfactorily the experimental results (Fig. 3) in terms of time evolution of block acceleration a z (Fig. 3a), velocity v z (Fig. 3b) and displacement u z (Fig. 3c).
The standard geotechnical parameters are the soil specific weight , the critical state soil friction angle ', the initial relative density D r0 of the soil stratum.
The elastic parameters are the non-dimensional Janbu's coefficients k E , n E (Janbu 1963) and the Poisson's ratio v whereas the damping parameters are N , T and R (Sieffert and Cevaer 1991).
The visco-plastic parameters are , , , vt , vr , v , c v , Δ 1 and Δ 2 . , and govern the shape of the yield function, vt and vr introduce in the hardening rule the dependence on the generalized tangential displacement and rotation, respectively, while v , c v , Δ 1 , and Δ 2 govern the shape of the viscous nucleus.  Between these 18 parameters, the 3 elastic and the 3 damping parameters have been fixed according to Dattola et al. 2021, whereas the 3 geotechnical standard parameters ( , ',D r0 ) have been fixed equal to those characterizing the DEM subdomain. Among the 9 visco-plastic parameters, 5 ( , , , vt , vr ) are fixed and correspond to those employed in di Vecchiotti 2006 andDattola et al. 2021 whereas 4 ( v , c v ,Δ 1 ,Δ 2 ) have been calibrated in order to reproduce the reference impact test.
The comparison between experimental data and BIM-PAM model predictions for the reference impact test, is reported in Fig. 3 in terms of temporal evolution of a z , v z and u z . The agreement is more than satisfactory.

Validation
Both DEM-FDM and BIMPAM models have been calibrated on the reference test results, characterized by a relatively small value of the initial boulder kinetic energy and a vertical velocity direction. Here below both models are employed to reproduce the experimental data of Calvetti and di Prisco 2012, concerning different falling heights (Fig. 6).
As was expected, the higher the falling height H , the larger both the peak acceleration a z,peak (Fig. 6a) and the maximum penetration u z,max (Fig. 6b).
As is evident, for low energy impacts ( H < 10m ), the agreement among experimental, DEM-FDM and BIM-PAM results is very satisfactory. In case of higher energy impacts ( H>10 m), DEM-FDM model tends to overestimate both a z,peak and u z,max . The overestimation of a z,peak depends, according to the authors, by a limit of DEM-FDM model not accounting grain damage or crushing for. This item is discussed in Appendix B. The overestimation of the final block penetration is, according to the authors, a consequence of the inability of the DEM numerical model of correctly reproduce the dissipation of energy. This item is discussed in Sect. 5. On the other side, BIMPAM predictions remain very accurate also for high energy impacts, despite the occurrence of grains damage phenomena is not explicitly taken into account by the rheological model. This is due to the choice of the viscous nucleus introduced by di Prisco and Vecchiotti 2006 in the formulation of BIMPAM inspired to projectile penetration experimental data (Boguslavskii et al. 1996).

Micromechanical Discussion
In this section, to exploit the additional information obtained by performing the numerical DEM analyses, the reference impact test is described in detail: in particular, the temporal evolution of both failure mechanism (Sect. 5.1) and block penetration (Sect. 5.2) is discussed.

Failure Mechanism
One of the fundamental hypotheses of BIMPAM model concerns the development of a generalized failure mechanism, evolving with time, under the penetrating block. To discuss such a hypothesis, here below the authors analyze DEM numerical results from both kinematic and static points of view. In Fig. 7 particle velocity vectors are plotted at the time instants a-d of Fig. 3a. For the sake of clarity, being the problem axisymmetric, only a vertical cross section of the model is considered in Fig. 7a-d. At the very beginning (instant A, t = 0.005 s), particles next to the impact point experience a very large acceleration (Fig. 7a). The wave front progressively propagates spherically from the impact point (Fig. 7b). Subsequently, while the particles under the penetrating block are pushed downward, an upward movement is laterally observed (Fig. 7c, d). In Fig. 7c and d, a sort of generalized failure mechanism clearly appears, becoming dominant with respect to the propagation wave.
In Fig. 8, the evolution of vertical stress z (Fig. 8a) and stress ratio k = x ∕ z (Fig. 8b) are plotted for different spherical measurement domains of radius equal to 10 D 50 placed at a relative depth z∕R = − 1.33 and different normalized distances x∕R = 0, 1.11, 2.22, 3.33 from the vertical line passing Fig. 7 Field of particle velocities under the impact point at selected time instants (Fig. 3a) for the DEM-FDM reference test ( x − z section) through the impact point (Fig. 2c). The temporal evolution of vertical stresses allows to capture the spherical propagation of the compression wave (Fig. 8a). In Fig. 8b, the evolution of k with time is compared with both active and passive earth thrust coefficients k a and k p respectively, evaluated considering the macroscopic critical friction angle ( ' = 35° obtained from DEM compression triaxial test results performed on the same material, here omitted for the sake of brevity). As is evident, under the impact point, k decreases, and a pseudo-active stress state is attained. On the contrary, beside the impact point, a significant increase in k is observed and a pseudopassive stress condition is reached.  These results confirm the development of a generalized failure mechanism under the impact point, the lack of a clear localization pattern is instead justified by the continuous penetration of the block inside the stratum which is associated with a continuous increase in the contact surface between block and soil.

Block Penetration Assessment
As was previously observed, when the impact energy is sufficiently low, the penetration of the block u z,max (Fig. 6b) is correctly reproduced by both BIMPAM and DEM-FDM models. In contrast, for higher impact energies, the DEM-FDM model seems to overestimate the experimental datum. To better understand the reason of this discrepancy, in Fig. 9a the block penetration versus time response for two different falling heights ( H = 4.9 − 38.6m ) is compared. The solid red line, corresponding to the DEM-FDM prediction, is characterized by an initial rapid penetration, followed by a plateau and a subsequent unexpected increase in u z . In contrast, experimental test results are characterized by a block rebound. The numerical plateau, not observed in case of H = 4.9m , is associated with a strong reduction in normal stresses under the impact point and an almost nullification of contact coordination number Z , already observed by Bourrier et al. 2008 andBourrier et al. 2010 (Fig. 9b).
According to the authors, this numerical response is due to the incapability of the system of correctly reproducing energy dissipation. Some authors have tried to overcome this problem by adding normal and tangential dashpots ("viscous damping") at particles contacts (Malone and Xu 2008). Although this strategy allows to both increase the dissipated energy and better reproduce the rebound of the block even for high energy impacts, the use of damping is not physically based and requires the definition of additional artificial viscous parameters (O' Sullivan 2011).

Inclined Impacts
As was mentioned in Sect. 2.4, inclined impacts have been very rarely studied in the literature; in particular, both numerical and large-scale experimental test results relative to inclined impacts are not available. In this section, DEM-FDM results are employed to test the reliability of BIMPAM model. Impacts characterized by the same initial velocity modulus ( V 0 = 9.80m∕s corresponding to a falling height H = 4.9m ) and different 0 (Fig. 2b) values, ranging between 0 and 60°, have been simulated.
In Fig. 10, the results, obtained for 0 = 20, 40, 60 • with DEM-FDM and BIMPAM models, are compared in terms of temporal evolution of a z , a x , v z , v x , u z and u x . The agreement seems to be satisfactory. As 0 increases, the peak of vertical acceleration (Fig. 10a) decreases, whereas horizontal acceleration increases (Fig. 10b). The peak instants of time of both a z and a x increase with 0 . The final value of the vertical velocity is significantly smaller than the initial one (Fig. 10c), whereas a lower reduction in v x is observed (Fig. 10d). In all the cases, the acceleration peak time simulated by BIMPAM is slightly larger than the one predicted by the DEM-FDM model. DEM-FDM model predicts a vertical penetration larger than the one predicted by BIMPAM (Fig. 10e). This discrepancy should be probably not only a signature of the already observed phenomenon but is also a consequence of the fact that the block enters in contact with a too small number of particles.
The comparison in terms of block trajectory and load path is illustrated in Fig. 11a and b, respectively. The two models predict very similarly the progressive deviation of the trajectory due to the interaction with the soil stratum. Analogously, even the generalized loading paths of the two models are very similar: the irregularity of the generalized loading path obtained using the DEM-FDM model is due to the numerical oscillations also evident in Fig. 10.
Finally, in Fig. 12a and b, the results for the complete set of inclinations are compared in terms of a z,peak ∕a z,0 and a x,peak ∕a z,0 where a z,0 is the peak vertical acceleration in the case of vertical impact. In Fig. 12c and d, the comparison is done in terms of vertical e z and horizontal e x restitution coefficients, as defined in Heideinreich and Labiouse 2004. The agreement between the two models' predictions is more than satisfactory. Both models predict, in agreement with what experimentally observed by Heideinreich and Labiouse 2004, restitution coefficients ( e z and e x ) to be very slightly dependent on the initial velocity direction.
To discuss the role of rotational energy in Fig. 13 for all the tests the ratio E kin,rot,max ∕E kin,in , being E kin,rot,max the maximum rotational kinetic energy of the block and E kin0 the initial kinetic energy ( E kin,in = 1 2 mV 2 0 ) is plotted. As is evident, even for the large inclinations, the rotational contribution remains negligible with respect to the translational one, testifying that, when the shape of the block is perfectly spherical, its rotation does not significantly affect the impact kinematics.

Conclusions
In this paper, inclined impacts of spherical rock blocks on horizontal granular strata have been investigated by using a coupled DEM-FDM numerical model. The obtained numerical DEM-FDM results have allowed to (i) discuss how the initial velocity inclination affects impact kinematics and (ii) provide a dataset to discuss the predictions of a previously conceived upscaled constitutive relationship (BIMPAM model) for any impact inclination. The agreement between DEM-FDM data and BIM-PAM predictions is very good, testifying the potentiality of the upscaled model in simulating dynamic soil-block interactions.
DEM-FDM results have also been discussed to confirm a fundamental hypothesis of BIMPAM, according to which, in case of relative dense strata, a generalized failure mechanism develops under the block footprint.
The problem of DEM-FDM data reliability has been critically tackled and the authors have proposed an original procedure for the generation of the coupled DEM-FDM model and for the calibration of the constitutive/contact parameters, based on large scale experimental test results. The authors have highlighted the role of grain size distribution in affecting numerical results and their dispersion.
In case of vertical impacts on horizontal strata, owing to a direct comparison with large scale experimental test results, the authors have also highlighted the role of grain damaging when the impact energy is sufficiently large.

Appendix A
See Fig. 14.   Fig. 13 Comparison between DEM-FDM numerical predictions (red crosses) and BIMPAM results (blue circles) in terms of the energetic ratio between the maximum rotational kinetic energy and the initial kinetic energy (Color figure online) A sensitivity analysis has been carried out to highlight the influence of grain size distribution on DEM-FDM numerical results. In particular, the authors assessed the dependence of numerical results on (i) the dimensions of grains ( D 50 ) and (ii) the material polydispersity ( C u ). To this aim, the authors have used the constitutive parameters of Tables 2 and 3, imposed in all the models analyzed a value of relative density Dr 0 ≅ 90% and considered only one falling height ( H = 4.9m ) and one block radius ( R = 0.45m).
(i) The authors have assessed that for D 50 ≤ 0.2R , since size effects are inhibited and the dispersion of results is acceptable, the results are not significantly affected by particle dimensions (the corresponding results are not reported for the sake of brevity) providing that a scaling rule as a function of particle size is adopted for contact stiffness (the ratio K n ∕D 50 must remain constant, as stated by Calvetti 2008); (ii) The results of the parametric analysis on material polydispersion have allowed the authors to derive that the system mechanical response depends very slightly on grain size distribution, provided that 1.25 ≤ C u ≤ f R∕D 50 (Fig. 14). The first threshold ensures the grain size distribution to be sufficiently polydisperse to mimic the response of a real granular material, whereas the second limitation prevents the presence in the material of too big grains. For this reason, this second limitation is not fixed but is a function of the ratio R∕D 50 employed, in the case taken into account f R∕D 50 = 5 = 2.

Appendix B
See Fig. 15 and Table 5. In Sect. 4, the authors have shown that DEM results seem to overestimate a z,peak in case of high energy impacts (Fig. 6a). This discrepancy may be justified by observing that DEM-FDM simulations do not account for particle damage phenomena, such as grain crushing, typically occurring in real materials when subject to large stress values. To clarify this aspect, an additional numerical campaign has been performed by accounting for different vertical impact tests with both rigid and crushable grains.
The authors have simulated grain crushing by using the breakage criteria and splitting procedure proposed by Ciantia et al. 2015, where an Hertzian contact rule was used. The reference material is in this case the Fointainbleau Sand (Yang et al. 2010), for which the micromechanical parameters have been calibrated by Ciantia 2021. An upscaling factor of 250 for the size of the particles has been employed. As in Sect. 3.1.1, an external elastic FDM sub-domain, whose parameters are reported in Table 5, is introduced. By following the strategy described in Sect. 3.1.1, the geometric dimensions of the numerical model have been fixed equal to b = 3m and L = 6m.
As is evident in Fig. 15, the numerical results obtained by assuming grains to be crushable (filled circles) are more satisfactorily fitted by the empirical relationship introduced by Labiouse et al. 1994 and verified for large impact energies by Calvetti and di Prisco 2007.

Fig. 15
Comparison between the predicted DEM-FDM peak acceleration for the different falling heights considering rigid (empty circles) and crushable particles (filled circles). The solid line represents the empirical correlation proposed by Labiouse et al. 1994