Probabilistic calibration of discrete element simulations using the sequential quasi-Monte Carlo filter

The calibration of discrete element method (DEM) simulations is typically accomplished in a trial-and-error manner. It generally lacks objectivity and is filled with uncertainties. To deal with these issues, the sequential quasi-Monte Carlo (SQMC) filter is employed as a novel approach to calibrating the DEM models of granular materials. Within the sequential Bayesian framework, the posterior probability density functions (PDFs) of micromechanical parameters, conditioned to the experimentally obtained stress–strain behavior of granular soils, are approximated by independent model trajectories. In this work, two different contact laws are employed in DEM simulations and a granular soil specimen is modeled as polydisperse packing using various numbers of spherical grains. Knowing the evolution of physical states of the material, the proposed probabilistic calibration method can recursively update the posterior PDFs in a five-dimensional parameter space based on the Bayes’ rule. Both the identified parameters and posterior PDFs are analyzed to understand the effect of grain configuration and loading conditions. Numerical predictions using parameter sets with the highest posterior probabilities agree well with the experimental results. The advantage of the SQMC filter lies in the estimation of posterior PDFs, from which the robustness of the selected contact laws, the uncertainties of the micromechanical parameters and their interactions are all analyzed. The micro–macro correlations, which are byproducts of the probabilistic calibration, are extracted to provide insights into the multiscale mechanics of dense granular materials.


Introduction
The discrete element method (DEM) captures the collective behavior of a granular material by tracking the kinematics of the constituent grains [9]. From just a few micromechanical parameters, DEM can provide comprehensive cross-scale insights [2,6,15,42] that are difficult to obtain in either stateof-the-art experiments or sophisticated continuum models. However, automated and systematic calibration of these parameters against macroscopic experimental measurements is still lacking, and often takes significant effort from DEM analysts.
Assuming homogeneous macroscopic deformation, the effective elastic properties of an ideal granular packing can be derived analytically from contact mechanics theories, micromechanical parameters and the microstructure of 4 Graduate School for International Development and Cooperation, Hiroshima University, 1-5-1, Kagamiyama, Higashi-Hiroshima 739-8529, Japan the packing [26,31,35,46], which to some extent facilitates the calibration of DEM models with analytical macromicro relations. Experimental studies were also conducted to improve the theories based on measured micromechanical [13,37] and microstructural behaviors [19,36] of grains and granular assemblies. Nevertheless, calibrating contact laws using microscopic measurements to reproduce the micromechanics at contacts does not necessarily recover the mechanical behavior at the macro-scale. Furthermore, in industrial applications, physically based contact laws (e.g., the classical linear (CL) [9] and Hertz-Mindlin (HM) noslip contact laws [23]) are often coupled with fictitious non-physical parameters, e.g., rolling and twisting stiffnesses [21,22,52], which makes the calibration a daunting task. To tackle the issues mentioned above, an "inversemodeling calibration" approach [1,27] is needed to infer the micromechanical parameters from macroscopic experimental measurements.
The conventional calibration procedure for DEM simulations of granular materials employs a "one at a time" analysis of the micromechanical parameters. Many parametric studies were conducted to derive micro-macro interpolation charts for various contact laws and materials using the "one at a time" approach. For DEM simulations of rocks [7,27,50] in which intergranular forces depend linearly on relative displacements between adjoining grains, a linear correlation between the bulk Young's modulus and the normal contact stiffness was identified, whereas a nonlinear correlation was found for sandy soils [38]. Given a constant value for the normal contact stiffness, both Young's modulus and Poisson's ratio were found to be linearly related to the tangential contact stiffness [4,38], despite the fact that normal and tangential stiffnesses can jointly affect the deformability of a granular material. The micromechanical parameters that characterize deformability (e.g., contact stiffnesses) and yield (e.g., intergranular friction angle), however, are generally believed to be uncoupled, and thus can be calibrated separately. This has led to many parametric studies into the macroscopic internal friction angle which was proved to have a nonlinear dependency on the intergranular friction angle [7,45,49], irrespective of the contact stiffnesses.
Because several micromechanical parameters can collectively determine the bulk behavior of a granular material, a parameter set identified for a DEM model with a known grain configuration can be regarded as one of the numerous solutions to the parameter identification problem. Among a handful of systematic approaches to calibrating DEM models, the design of experiments (DOE) methods are efficient in searching for possible solutions in the multi-dimensional parameter space, with a manageable number of DEM simulations [18,24,39,54]. Hanley et al. [18] applied DOE to calibrate DEM models of crushable agglomerates. The interaction between the key parameters was considered by the orthogonal arrays designed using the Taguchi methods. Yoon [54] developed a two-step optimization process in which a DOE method (Plackett-Burman design) was first applied to select the parameters that have the largest impacts on macroscopic material properties. The statistical micro-macro correlations were then determined by running additional DEM simulations. Despite a good agreement on the instantaneous macroscopic material properties (e.g., compressive strength), the transient response in the simulation (e.g., stress-strain curves) did not agree well with the experimental data. The efficiency of DOE methods largely depends on the prior knowledge of the interaction between parameters, which is needed to prepare DOE samples, but is also generally limited due to the diversity of granular materials.
The aforementioned calibration approaches aim to calibrate micromechanical parameters against the macroscopic material properties (e.g., Young's modulus, peak-and criticalstate friction angles), rather than the transient behavior of the bulk material. This is very likely to hinder the predictive capacity of DEM models. Local rather than global solutions to the parameter identification problem might be discovered and adopted in the models. This leads to deviations from the transient experimental response in the simulations. In addition, the dependency of granular materials on stress and fabric history cannot be accounted for in these approaches. To capture the elastoplasticity of granular materials, the experimental data measured throughout the entire loading history needs to be fully considered within the parameter identification procedure [41].
The sequential data assimilation techniques [12,34] can be applied to solve the "inverse-modeling calibration" problems and to overcome the above-mentioned difficulties. For the current inverse problem, the sequential quasi-Monte Carlo (SQMC) filter 1 and sequential importance sampling (SIS) are implemented, which can jointly account for the effects of loading history on the elastoplastic behavior of granular materials [43]. The SQMC filter applies the recursive formula of sequential Bayesian estimation. Experimental data obtained step by step during a loading process are assimilated into DEM models to approximate the evolution of posterior PDFs in the corresponding parameter space. This "inversemodeling calibration" approach is expedient, because the synergy of the SQMC and SIS algorithms is well-justified for nonlinear and non-Gaussian geomechanical problems, as demonstrated in the applications of these methods to inverse finite element analyses [33,43]. A few innovative calibration methods for DEM or molecular dynamics simulations were recently proposed using Bayesian approaches [3,40,53]. Hadjidoukas et al. [16,17] employed the Transitional Markov Chain Monte Carlo algorithm to find optimal parameter values in DEM simulations of silo discharge. By assuming uniform prior distributions for all micromechanical parameters, the parameter uncertainties were quantified and a Bayesian model selection framework was provided for two-dimensional (2D) monodisperse granular systems. The SQMC filter implemented in this study requires no assumptions about the prior distributions of the model states. The Bayesian calibration is conducted for three-dimensional (3D) DEM simulations of polydisperse granular materials under triaxial compression, which in turn reveals the posterior PDFs that help in assessing the robustness of the contact laws. To the authors' knowledge, this work is the first attempt to develop a sequential data assimilation procedure for calibrating DEM models over the transient behavior of bulk granular materials. For the current implementation, both SQMC and SIS algorithms are implemented in the open-source DEM package YADE [44].
The remainder of this paper is organized as follows. Section 2 outlines the contact laws and grain configurations which are the key ingredients of the DEM models. The fundamentals of SQMC and SIS, and their implementation, are detailed in Sect. 3, followed by the descriptions of the posterior PDFs resulting from the transient experimental data in Sects. 4 and 5. Section 6 discusses the robustness of the proposed method. Conclusions are drawn in Sect. 7.

DEM simulations of granular materials
The open-source DEM package YADE [44] is applied to perform 3D DEM simulations of dense granular materials under drained triaxial compression. The stochastic responses obtained from the simulations using quasi-randomly sampled parameter values, i.e., samples, are processed through the SQMC and SIS algorithms to approximate the posterior PDFs, conditioned to the experimental data. YADE models granular materials as packings of solid grains with simplified geometries and vanishingly small intergranular overlaps. It tracks the trajectory of each grain within the explicit time integration scheme, based upon the net force resulting from the intergranular forces. In the present work, the robustness of the SQMC filter for calibrating DEM models is examined by modeling a representative volume of cohesionless granular soil using dense packings of polydisperse spherical grains. The simple HM and CL contact laws are selected to describe the contact force-displacement relationships in the normal and tangential directions. To account for the effects of grain shape and roughness, rolling and twisting moment transfer is allowed at the contact points of adjoining spheres. Both inter-granular tangential forces and rolling/twisting moments are assumed to be bounded by Coulomb type yield criteria. Each DEM simulation of the representative volume of the granular soil under drained triaxial compression is performed in a strain-controlled quasistatic manner. For each calculation cycle of an incremental loading, a global damping ratio of 0.2 is adopted to mimic the presence of a background fluid and ensure numerical stability. In the subsequent cycles, which are run to dissipate the total kinetic energy before the next incremental loading, the global damping ratio is raised to 0.9 until the quasistatic macroscopic variables (see Sect. 3.1) can be extracted.

Micromechanical contact laws and parameters
Two different contact laws are employed in the present work, i.e., the classical linear force-displacement model and the nonlinear model which combines the formulations for Hertzian normal and Mindlin's simplified tangential stiffnesses. For two adjoining spheres with a negligible normal overlap u n and a tangential relative displacement du s at the contact point, the intergranular normal force F n and tangential force increment dF s can be related to the normal stiffness k n and tangential stiffness k s respectively. The general forms of these force-displacement models are given by where k n and k s are defined differently according to the contact laws. In YADE, the linear normal and tangential spring stiffnesses are computed from the harmonic averages of the contact-level Young's moduli E c , Poisson's ratios ν c and the radii associated with the two spherical grains. Adopting identical values for E c and ν c of both grains, the expressions for k n and k s are simplified as: where the equivalent radius r * is defined as 1/(1/r 1 + 1/r 2 ), and r 1 and r 2 denote the radius of the two grains in contact. By adopting identical E c and ν c for the contacting spherical grains the formulations of k n and k s given by the HM contact law read: where the values of k n and k s depend nonlinearly on u n , which is updated in every calculation cycle of DEM sim- ulations. The respective expressions for the normal and tangential stiffnesses according to the two contact laws implemented in YADE are summarized in Table 1.
Because the grain shapes are simplified as rigid spheres, rolling resistance between adjoining grains is needed to mimic the effects of grain shape and surface roughness. One additional phenomenological parameter, namely the rolling stiffness k m , is thus introduced, which linearly relates the rolling/twisting moment in the contact area with the relative rotation/torsion between two contacting grains. Note that rolling and twisting are unified in one moment-rotation equation for simplicity. In YADE [32,44], the rolling stiffness defined in the implementations of the HM and CL forcedisplacement models varies slightly (see Table 1). A constant rolling stiffness k m = β m is implemented in the former, whereas the latter uses β m as a dimensionless coefficient to calculate k m from the radii and contact tangential stiffness of the two contacting grains.
Friction criteria are imposed on intergranular tangential forces and rolling/twisting moments to account for sliding between contacting grains and, on the macro-sale, plastic deformation in the granular soil. Both failure criteria adopt the Mohr-Coulomb type formulation. The maximum tangential force F s is calculated as the intergranular friction coefficient tan(φ μ ) multiplied by the magnitude of the normal force F n , whereas the rolling/twisting moment is constrained by the rolling friction coefficient η m , F n and the characteristic length related to the radii (see Table 1). Despite being phenomenological, rolling stiffness and friction are generally needed to obtain internal friction angles close to experimentally measured values for dense granular soils, without making an extra computational effort to capture the kinematics related to irregular grain shapes and surface roughness. Note that although the contact laws employed in the current work are implemented differently, the micromechanical parameters to be identified in the Bayesian calibration are the same, i.e., the Young's modulus E c and Poisson's ratio ν c of contacting grains, the rolling stiffness and friction coefficient β m and η m , and the intergranular friction angle φ μ .

Scale and resolution of DEM granular packings
In the present work, the DEM granular packings are prepared to model a representative volume of Toyoura sand, such that the packings can be repetitively stacked up to construct a large-scale simulation domain. The micromechanical force and contact networks calculated with DEM are typically averaged over the repeated representative volumes to extract macroscopic variables such as stress and fabric tensors [5]. Alternatively, these packings can be embedded at the Gauss integration points of a FEM mesh to derive the local material responses in a hierarchical FEM×DEM multiscale model [15,30]. It is well known that the microstructure of a granular material (e.g., coordination number and anisotropy [28,29]) plays a key role in the macroscopic constitutive behavior. Therefore, the number of constituent spherical grains N g , namely, the resolution of the DEM packings needs to be sufficiently large in order to ensure a realistic representation of the microstructure. On the other hand, in those cases in which the DEM packings are designated to return the material responses at a local scale, N g should not be excessively large so as to suppress localized cracks and failures in the microstructure.
In the light of the above considerations, grain configurations consisting of various numbers of spherical grains are created and investigated. For each grain configuration, a cloud of spherical grains with the same density ρ g = 2650 kg/m 3 is first generated in a cuboidal periodic cell. The diameters of the spherical grains are sampled from a scaled grain size distribution of Toyoura sand. Initial isotropic compression is then applied to consolidate each cloud into a dense cuboidal packing (50 mm × 50 mm × 100 mm) with an initial void ratio of e 0 = 0.68, the same as those in the experiments [47]. During the consolidation stage, a very high Young's modulus is adopted to create grain configurations with negligible overlaps between grains. While maintaining a low confining pressure ( p = 100 Pa) by using a servo-controlled periodic boundary condition, the initial intergranular friction is tuned down slightly whenever quasistatic states are reached. After a couple of cycles of reducing the friction and minimizing unbalanced forces, various "stress-free" dense packings, which are "identical" in the sense of the initial void ratio, are generated. Among all the randomly generated dense packings three (N g = 1000, 2000 and 5000) are selected as illustrated in Fig. 1, because of their relatively smooth and spherical contact orientation diagrams after the initial isotropic compaction. Note that N g = 1000 is the minimum number of spherical grains needed to create "stress-free" dense packings, whose contact orientation distributions are uniform and statistically stable.

State space model and state estimation
The SQMC filter is a sequential data assimilation method that is known to be particularly preferable for merging sparse experimental data into prognostic models of nonlinear and non-Gaussian systems. It takes advantage of the recursive formula of the sequential Bayesian framework to approximate the posterior probability distributions using an ensemble of sampling points in the multi-dimensional parameter space R n . When applied to estimate the expectations of the state variables together with the importance sampling, the recursive Bayesian filter is capable of tracking the evolution of the probability on each sample, conditioned to the sequentially measured experimental data. In fact, the process of updating the posterior PDFs based on the existing and newly obtained knowledge from experimental investigations is highly suited for calibrating the DEM models of granular materials which are sensible to stress history. Given the Toyoura sand specimen (e 0 = 0.68) and its equivalent DEM granular packings, the physical measurements on the specimen and the numerical predictions resulting from the DEM models can be described in a nonlinear and non-Gaussian state space model [25]: where the state vector x t consists of three independent variables which characterize the triaxial responses of the DEM packing, namely, stress ratio σ a /σ r , radial strain ε r and volumetric strain ε v at a discrete data assimilation step t, whereas the observation vector y t is directly measured in the drained triaxial compression experiments [47]; ν t and ω t are the system error and the observation error respectively, whose PDFs follow normal distributions with zero means. F denotes the operator that represents the nonlinear dynamic model (i.e., the DEM model). The current state of F depends on all preceding states of the dynamic model, which evolves with physical constraints like externally applied loads. The nonlinear observation operator H is reduced to an identity matrix of size three in the current problem.

SQMC filter
The SQMC filter approximates posterior PDFs via a set of samples (referred to as an ensemble) and the associated weights [25,43]. From the ensemble {x (9) where N p is the number of samples, δ is the Dirac delta function, and the superscript (i) indicates the ith sample x (i) t−1|t−1 and its associated weight w (i) t−1 . It should be noted that all the weights are no less than zero and the summation of them must be one, i.e., 0 ≤ w (i) t−1 ≤ 1 and With the help of the ensemble sampled using the quasi-Monte Carlo method, the recursive formulas for approximating the onestep-ahead forecast distribution p(x t |y 1:t−1 ) and the filtered distribution p(x t |y 1:t ) are simplified as follows: p(x t |y 1:t ) = p(y t |x t ) p(x t |y 1:t−1 ) p(y t |y 1:t−1 ) Given that the observation system is assumed to be linear, the likelihood p(y t |x t|t−1 ) of the ith sample reads: where m is the dimension of the state vector, M t is a predetermined covariance matrix of the observation error, and H t is the observation operator H in matrix form. From the previous weight w (i) t−1 and the current regulated weightw (i) t on each sample (see Eq. 12), the new weight w t (i) at the data assimilation step t can be updated as:

Sampling method and SQMC filtering procedure
The SIS algorithm keeps track of the evolution of the initially generated samples, rather than repeatedly performing resampling from the updated posterior distribution. It allows for calibration against the transient behavior of the modeled system, e.g., the stress-strain response of granular soils, through the filtering and weight updating steps. The initial sampling and the SQMC filtering procedure for solving the current parameter identification problem are summarized in the four steps below.

Initialization:
Generate an ensemble of realizations {x 0 (1) , x 0 (2) , . . . , x 0 (N p ) } by initializing dynamic models (DEM simulations) using parameter sets sampled from a lowdiscrepancy sequence in R n . n is the number of the parameters to be identified and is equal to five in the present work.

Prediction:
Update the state of each realization x t (i) from the corresponding DEM simulation run at the data assimilation step t.

Filtering:
Given the experimental data y t , calculate the weights on } that represents the "fitness" of the dynamic models to the physical system using Eqs. (12) and (13). 4. Weight updating: Approximate filtered distribution p(x t |y 1:t ) with the updated ensemble and the associated weights. Return to Step 2 with p(x t |y 1: } for prediction and filtering at the next data assimilation step. The SQMC filter makes no assumptions on how the prior probabilities of the dynamic model distribute over the prior samples. Therefore, a sufficiently large ensemble size is required for the posterior PDFs to stabilize within the available data assimilation steps. For the current numerical models, which have five micromechanical parameters, a total number of 2000 samples (N p = 2000) can efficiently estimate the posterior PDFs at acceptable computational cost. The fact that the DEM simulations are run in parallel on a multi-core system with the open-source DEM package YADE also helps in reducing the total computational time needed in the prediction step.

Numerical and experimental triaxial compression tests: the dynamic model and the physical system
The calibrations of the DEM simulations of Toyoura sand under drained triaxial compression are adopted as an example, to demonstrate the capability of the SQMC and SIS algorithms in evaluating the parameter uncertainties of the DEM models. Based upon the wide ranges assumed for the micromechanical parameters in Table 2, a five-dimensional Halton sequence [10] is generated and employed as the ensemble which contains 2000 samples (i.e., combinations of parameters E c , ν c , β m , η m , and φ μ ). The values held within each sample are correspondingly assigned to the micromechanical parameters of the CL and HM contact laws. The HM and CL contact laws are applied respectively to the three "stress-free" initial packings (N g = 1000, 2000 and 5000), which are created as the computational representative volumes of the Toyoura sand specimen (e 0 = 0.68). Drained triaxial compression on the Toyoura sand specimen under various confining pressures are simulated in a two-phase loading program: the packings are first isotropically com- pressed to prescribed confining pressures σ c (see Table 3), and then sheared along triaxial loading paths in a quasistatic manner. The DEM simulations are strain-controlled, so that the state variables can be extracted at the exact macroscopic strain levels where the experimental measurements were made. The probabilistic calibrations of six DEM models (combinations of two contact laws and three grain configurations) are conducted under various confining pressures as listed in Table 3. Each probabilistic calibration makes use of 2000 deterministic DEM simulations run with the same set of samples (N p = 2000). As the simulations undergo each incremental load, the extracted state variables and the corresponding experimental measurements are processed through the one-step-ahead forecasting and filtering steps (Eqs. [10][11]. As a result, the evolution of the posterior PDFs can be captured over the loading history.

Posterior probabilities of parameters in two DEM models
To understand the robustness of the proposed calibration approach, the effect of grain configurations and loading conditions on the posterior probability distributions are investigated (see Table 3). The evolutions of the posterior probabilities of the micromechanical parameters are first presented for both contact laws shown in Fig. 3. The goal is to present an global picture of numerous solutions to the current inverse problem and their evolutions over the loading history. The dependencies of the posterior PDFs on the selected grain configurations and loading conditions are then analyzed in Sects. 5.1 and 5.2, by fixing either grain configurations or loading conditions in the DEM simulations. Kernel density estimators are applied to postprocess the discrete samples and their associated importance weights into smooth posterior PDFs. Note that the goal is not to reveal the probability density at a specific point in the parameter space, but to inter-pret the robustness of the contact laws from the posterior distributions.

Identified micromechanical parameters
To keep track on how the expectations of the state variables converge to their true values, the quasi-randomly generated parameters are weighted by the updated posterior probabilities at each data assimilation step (axial strain increment dε a = 0.1%). The evolutions of the weighted averages, referred as identified parameters, are shown in Fig. 2. The posterior probabilities obtained from the three DEM packings (N g = 1000, 2000 and 5000) and the two contact laws are compared in Figs. 2 and 4, to investigate the effect of the DEM models on the calibration process. The DEM packings undergo the exact triaxial loading history as the experimental specimen at a confining pressure of σ c = 0.2 MPa. After recursively updating the weights through sufficient steps of data assimilation, the SQMC filter eventually leads to less pronounced variation of the posterior probability distributions. Figure 2 shows that the identified parameters reach plateaus after approximately 60 data assimilation steps. The combinations of parameters that are estimated to give the highest posterior probabilities are considered to be the calibration results. Although the evolutions of the identified parameters become stationary eventually for both contact laws and all three DEM packings, the convergence rate appears to be lower with larger fluctuations for the DEM models governed by the HM contact law (dashed curves). Note that the true values, which the identified parameters converge to, are very close between the DEM models that only differ in their grain configurations. The agreement between the identified true values obtained with various grain configuration is even better for the DEM models which adopt the CL contact law. The identified φ μ and η m for the simulations conducted with the HM and CL laws differ slightly, whereas the differences between the identified values for E c , ν c and β m of the two contact laws are more pronounced. This is because the same Coulomb type friction is applied on both tangential forces and rolling/twisting moments. It appears that the true values identified for the strength parameters φ μ and η m are independent of the normal and tangential stiffnesses, although they are formulated differently in the two contact laws.

Posterior PDFs of micromechanical parameters
A key advantage of the SQMC filter for parameter identification is the ability to capture the posterior PDFs in the multi-dimensional parameter space. Because of the impor- Therefore, the calibration must be able to take the loading history into account, which is only possible within a sequential Bayesian framework using the SIS. The subplots CL(a)-(e) and HM(a)-(e) in Fig. 3 show a typical evolution of the weights, calculated using Eqs. 12-14, associated with the five micromechanical parameters conditioned to the experimental triaxial behavior of Toyoura sand. The subplots CL(a), (e) and HM(a), (e) of Fig. 3 clearly show that the posterior probability distributions projected on the E c and φ μ axes gradually shift from uniform to Gaussian-like for both contact laws, as more experimental data becomes available. It is surprising to see that the true values for E c and φ μ already become identifiable within the first 20 data assimilation steps (ε a = 2%). The weights associated with ν c in Fig. 3 CL(b) and HM(b) grow into bimodal and multimodal distributions at approximately the 40th data assimilation step (ε a = 4%), which corresponds to the strain where σ a /σ r reaches the peak in the experiment. After 60 data assimilation steps, there are no significant changes in the posterior distributions over E c , ν c and φ μ . Further data assimilation only causes a few samples to gain more weights. This follows the concept of Bayesian updating, but leads to the degeneracy problem, i.e., very few samples having relatively large weights. Refining the range of ν c or using a larger ensemble size would help better capture the bimodal/multimodal distribution. Although suitable values for β m and η m are identified, the posterior probability distributions appear to be more scattered and do not evolve with consistent shapes. The fact that these distributions can hardly be described by a Gaussian or mixtures of Gaussians might be due to the nonphysical nature of the rolling/twisting parameters. Note that the weights in Fig. 3 HM(a)-(e) distribute around the true values with bigger standard deviations compared with those shown in Fig. 3 CL(a)-(e), although the uncertainty assumed for the observation error is identical in all cases. This suggests that, given the current ensemble of prior samples, the HM contact law has a higher model robustness in predicting the triaxial behavior of granular soils than the CL contact law.
To better illustrate the dependency of the posterior probability distributions on grain configurations and loading conditions, the discrete weights as in Fig. 3 are further processed through kernel density estimation to obtain the smooth density functions shown in Figs. 4 and 6. The posterior PDFs obtained using different DEM packings (N g = 1000, 2000 and 5000) at the 60th data assimilation step (ε a = 6%) are plotted together in Fig. 4. Note that after the 60th step, there are no further changes in the evolutions of the identified parameters, except for β m and η m . As can be expected from Fig. 2, the posterior PDFs for the different DEM mod-els governed by the same contact law generally collapse, as shown in Fig. 4. The density functions are almost identical for the DEM models using the CL contact law, regardless of the numbers of grains N g in the packings, whereas those obtained using the HM contact law differ slightly depending on N g . One of the reasons might be that the rolling stiffness implemented for the HM contact law is not scaled by grain radii which increases as the number of grains in the packing decreases. The implementation of the rolling/twisting  governed by CL and HM contact laws, at the 60th data assimilation step scheme for the CL contact law, on the other hand, calculates the rolling stiffness from the radii of the contacting grains, the tangential contact stiffness and dimensionless factor β m . Nevertheless, the agreement between the posterior PDFs using the same contact law is generally good, and in line with most large-scale DEM simulations that use scaled grain size distributions to represent granular soils [6,14,20]. The agreement of the posterior distributions in Fig. 4 verifies a scaling rule [48]: the macroscopic behavior of a granular material can be recovered from a "prototype" DEM packing with a minimal number of grains, as long as the same bulk properties, such as initial void ratios and stress states, are ensured for the DEM packing. Based on the scaling rule, fast probabilistic calibration of micromechanical parameters can be first conducted with the "prototype" DEM packing. Duplicates of the "prototype" packing can then be either employed as representative volume elements for FEM × DEM multiscale modeling, or stacked up to assemble a large-scale DEM simulation domain.

Identified micromechanical parameters
It is known that the elastoplastic behavior of granular soils strongly depends on the loading history experienced. With DEM, the plastic deformation of granular soils can be easily captured via irreversible change of the microstructure subjected to external loads. The questions are how the loading history can be taken into account during the calibration process such that the correct elastoplastic behavior is reproduced at the macro-scale, and how the identified micromechanical parameters of both contact laws would be affected by a range of loading conditions. Figure 5 shows

Posterior PDFs of micromechanical parameters
In a similar fashion to the posterior probability distributions described in Sect. 5.1, the discrete approximation of posterior probabilities are smoothed using Gaussian kernels. A closer look at the posterior PDFs in Fig. 6 CL(a) and HM(a) shows that the PDF projected over E c of the CL contact law shifts slightly to the right as σ c increases, whereas those of the HM contact law collapse between 10 9.5 and 10 10 Fig. 6 HM(e). In fact, there are still considerably large overlaps between the distributions obtained at various confining pressures, where the probability densities are relatively high. The above findings suggest that for a given grain configuration, the nonlinear HM contact law is more robust at predicting the triaxial behavior of granular soils under various confining pressures. This means that applying a certain combination of parameters for various confining pressures, in the case of the HM law, will not lead to a large discrepancy between the numerical predictions and experimental data, as long as the parameters are selected around the peaks of the joint distributions. These findings are reasonable, because the HM contact law is analytically derived from the elasticity of two contacting spheres which provides the nonlinear dependency of contact stiffnesses on pressure. Using a well calibrated combination of parameters, the HM contact law should be able to more accurately model a granular material subjected to a wide range of external loads, until the effect of grain crushing is no longer negligible.

Numerical predictions versus experimental data
The DEM simulation results reproduced with the most probable parameters identified by the SQMC and SIS algorithms are plotted with the experimental data in Figs. 7 and 8. Good agreement between the numerical predictions and experimental data is achieved as shown in Figs. 7 and 8, regardless of the contact laws and the numbers of spherical grains in the DEM simulations. The combinations of parameter values identified to be the most probable in reproducing the experimental data at their respective confining pressures are listed in Tables 4 and 5. As can be seen in the tables, the intergranular friction angle which gives the best agreement in Fig. 8 decreases with an increase of the confining pressure for both contact laws. This dependency can be attributed to the crushing behavior of Toyoura sand in the experiments, which neither of the two contact laws can mimic in DEM simulations using smooth spherical grains. It should be noted that although the SQMC filter can give the most optimal combinations of parameters as listed in Tables 4 and 5, it is fundamentally different from any optimization techniques, because the probabilistic calibration is conducted by evaluating the parameter and model uncertainties through the posterior PDFs within the Bayesian framework. Each optimal solution presented here is simply one of many possible solutions to the parameter identification problem. The following section aims to interpret the interactions between micromechanical parameters from the posterior PDFs in the five-dimensional parameter space. The statistical micro-  macro correlations are also established following the same concept mentioned above.

Correlations between micro-and macro-mechanical parameters
The posterior PDFs obtained from the DEM simulations (N g = 1000 and σ c = 0.2 MPa) using both contact laws at three selected data assimilation steps are illustrated in Fig. 9a-c. At these characteristic steps, the granular material shows distinctive macroscopic physical states, namely, purely elastic, maximum volumetric-contraction and peak stress ratio states, as shown in Fig. 7. The objective is to analyze the evolution of the posterior PDFs in the respective parameter spaces, so that better prior knowledge of the interactions between the micromechanical parameters can be obtained. The color bars are omitted in Fig. 9 because only the relative values matter in solving the parameter identification problem. The plots in the diagonal of Fig. 9 present the marginals of the posterior PDFs of the five micromechanical parameters, i.e., E c , ν c , β m , η m , and φ μ , respectively. The projections of the continuous posterior PDFs obtained using the HM and CL contact laws are shown in the below-and above-diagonal panels, and colored by their respective densities. The blue and red color schemes are adopted for the distributions obtained using the CL and HM contact laws. Figure 9a shows that at a very early stage the posterior distributions and the peaks associated with E c and φ μ of Table 5 Most probable sets of micromechanical parameters identified for the HM contact law both DEM models are already identifiable. While these two parameters seem to be correlated, the other parameters ν c , β m and η m show less pronounced or insignificant correlations with E c and φ μ , irrespective of the contact laws. The 2D projections of the posterior PDFs associated with E c are almost aligned into straight bands perpendicular to the E c axis. The physical system, which evolves from the initial elastic to maximum volumetric-contraction states, causes the peaks over φ μ to shift towards larger values for both contact laws, as shown in Fig. 9b. However, the true values of E c are almost the same as those in the initial calibration stage, as shown in Fig. 9a. In Fig. 9a, the distributions in both 2D E c -φ μ projections become increasingly parallel to the E c axis, forming into long tails as E c increases. These long tails gradually disappear, as the physical system evolves from the purely elastic to maximum volumetric-contraction states, as can be seen for both contact laws in Fig. 9b. This indicates that if, by trial and error, the identified E c happens to be located at these local optima, which could happen because E c is commonly calibrated against the bulk Young's modulus, the resulting numerical simulations will not yield good agreement with the experimental results, even though φ μ could be well determined by the bulk shear strength. While the physical system is approaching the peak stress ratio state, the distributions keep shrinking towards the true values, with no evident shift in the parameter spaces. It should be noted that the calibrations at least for E c and φ μ of both contact laws have already finished at the maximum volumetric-contraction state. The Fig. 9 Posterior PDFs approximated by the SQMC and SIS algorithms at the data assimilation steps that correspond to three characteristic states in the stress-strain response. Diagonal panels indicate marginals of the joint posterior PDFs, whereas the below-and above-diagonal panels present 2D projections of the posterior PDFs obtained using the HM (red color scheme) and CL (blue color scheme) contact laws respectively. The Therefore, it can be concluded that the micromechanical parameters E c , ν c and φ μ have the predominant effect on the calibration of the DEM models. Once these three parameters are calibrated with high accuracy (posterior probability), the true values for the rolling parameters β m and η m can be identified by tuning the evolution of post-peak stress ratio over axial or deviatoric strain. It is worth noting that this technique is often adopted in the literature, because rolling stiffness and rolling friction can significantly affect the bulk shear strength, but have a negligible influence on the initial elastic behavior and dilatancy of a dense granular material [38].
A very important byproduct of the proposed calibration procedure is the large number of simulations that are needed to capture the complete picture of the posterior PDFs over the explored parameter space. From the huge amount of simulation data, although many of them largely deviate from the experimental results, it is possible to extract meaningful universal trends between the micro-and macro-mechanical parameters. From the numerically predicted stress-strain curves, the bulk Young's modulus and Poisson's ratio, E bulk and ν bulk are determined as the secant values at the devi-atoric strain ε d = 1%. The internal friction angle ϕ bulk is simply calculated from the peak stress ratio. Figure 10a, b respectively show the samples in the parameter subspaces constructed with the micromechanical parameters of the CL and HM contact laws and the corresponding bulk properties calculated directly from the numerical stress-strain curves. The samples in Fig. 10 are colored by the probability density functions estimated with Gaussian kernels.
From the statistics in Fig. 10a, b clear micro-macro correlations between the micromechanical parameters and the bulk properties of granular materials can be identified. The correlations between E bulk and E c are almost linear and piecewise linear for the CL and HM contact laws respectively in loglog scales. Although the correlation between E bulk and φ μ is not as significant as that between E bulk and E c , a clear dependency of the bulk Young's modulus on the intergranular friction can be seen from the samples with high probability densities. As shown in Fig. 10a, the bulk Poisson's ratio ν bulk predicted by the CL contact law might have unrealistic values above the physical limit. The statistics show that ν bulk is more likely to be unrealistic with the increase of E c and φ μ . All the bulk Poisson's ratios predicted by the HM contact law, on the other hand, fall within the physical range between 0 and 0.5. From the subplots that characterize the dependency of ν bulk on the micromechanical parameters in Fig. 10a, b, it can be observed that ν bulk is nonlinearly correlated to both E c and φ μ in the case of the HM contact law, whereas this trend is barely noticeable for the CL contact law. The ϕ bulk -φ μ correlations appear to be nonlinear and are very similar for both the CL and HM contact laws, except that the samples are more scattered in the former than in the latter. In both ϕ bulkφ μ subplots, there exist upper bounds on the possible values of ϕ bulk depending on the selection of φ μ . In the case of the HM contact law, a lower bound is also present which further narrows down the choices for the initial guess of φ μ . Nevertheless, the samples become more scattered and the two bounds diverge increasingly as φ μ increases, meaning that the calibration against large ϕ bulk is generally more difficult than that against small ϕ bulk . Because of the non-physical nature of the rolling parameters, no significant correlations can be found from the evenly distributed samples in the relevant 2D micro-macro parameter spaces. This means that it is impractical to establish meta-models relevant to the rolling parameters, which may undermine the efficiency of an optimization process. The SQMC filter can resolve this difficulty without initially knowing the interactions between the parameters being identified, which usually comes at the price of increased computational cost. Note that the SQMC filter or other sequential Monte Carlo methods can be applied iteratively using the knowledge of the interactions obtained in the filtering step. In such a way, the computational cost could be significantly reduced. It can also be observed in Fig. 10a, b that the data gathered from the DEM simulations using the HM contact law are less scattered that those using the CL contact law, which suggests that the HM contact law is more robust than the CL contact law within the parameter space currently explored.

Conclusions
A novel probabilistic calibration approach is proposed for the DEM simulations of granular soils. The SQMC and SIS algorithms are implemented within the open-source DEM package YADE. The micromechanical parameters of the contact laws are successfully calibrated against the stress-strain behavior of Toyoura sand in drained triaxial compression conditions at various confining pressures. Compared with general optimization methods, the synergy of the SQMC and SIS algorithms can estimate the evolution of poste-rior probability distribution over a loading history in a high dimensional parameter space. From the distribution of the estimated probability density functions, one can easily evaluate the robustness and uncertainties of DEM models and interactions between the micromechanical parameters.
The probabilistic calibrations of six DEM models (combinations of two contact laws and three grain configurations) are conducted considering confining pressures ranging from 0.2 to 2.0 MPa. Despite some fluctuations in the evolution of the identified parameters during the calibration, stabilized posterior probability distributions are obtained for both contact laws. The distributions of E c and φ μ evolve from being uniform to Gaussian-like, whereas the distributions of the other parameters can hardly be approximated by mixtures of Gaussians. The effect of grain configuration on the posterior PDFs obtained using either the CL or HM contact law under the same loading condition appears to be negligible. It is proved that the macroscopic behavior of a granular material can be recovered from a "prototype" DEM packing with a minimal number of grains, as long as the same bulk properties, such as initial void ratios and stress states, are ensured for the DEM packing. For a predetermined grain configuration (e.g., N g = 1000), the posterior PDFs obtained using the HM contact law under various loading conditions generally collapse near the consistent peaks in the parameter space. These findings indicate that although the scaling rule holds for both contact laws (perfectly for the CL contact law), the HM contact law is more robust in predicting the triaxial behavior of a dense granular material at a wide range of confining pressures. The scaling rule is of great importance, because in most laboratory-and industrial-scale applications the number of physical grains can easily exceed several millions. One reasonable approach to DEM modeling at these scales is to apply the scaling rule to reduce the computational effort. However, the up-scaling has to be carried out in such a way that the model resolution is still sufficient with negligible effects on the predicted bulk behaviors.
The SQMC and SIS algorithms are implemented as a separate calibration toolbox independent of the DEM package. It is straightforward to apply the proposed calibration approach to other DEM codes. Current work involves probabilistic calibration of dense granular materials under quasistatic loading conditions. In future research, it will be worth investigating the applicability of the SQMC and SIS to DEM calibrations against more complex behavior of granular materials, such as rockfalls [11], avalanches [51] and silo discharges [8]. Ongoing work on probabilistic calibration of DEM simulations is directed towards the development of an iterative SQMC filter, which is capable of resampling micromechanical parameters from prior probability distributions updated by the previously conducted probabilistic calibrations. The iterative version of the SQMC filter would allow the probabilistic calibration to learn from its preceding experience and explore a parameter space with different resolutions. By putting more samples in highly probable parameter subspaces, the computational cost can be greatly reduced.