Discrete Element Modelling for biocemented sand: effect of calcite distribution at the microscopic scale

. The mechanical e ﬃ ciency of the biocementation process is directly related to the microstructural properties of the biocemented soil, such as, the volume fraction of calcite, its distribution within the pore space (whether localized at the contact between grains or over the grain surfaces) and the contact properties: coordination number, contact surface area, contact orientation, type of contacts (frictional even after treatment, purely cohesive via a calcite bridge or combining friction between particles and cohesion of the localized calcite). Dadda et al, (2018) have used microscopic properties computed from 3D images obtained by X-ray tomography of biocemented sand samples with di ﬀ erent levels of biocementation as an input in current analytical models to estimate the elastic properties (Young’s and shear modulus) and the strength properties (Coulomb cohesion). They pointed out the important role of some microstructural parameters, notably those related to the contact, on such e ﬀ ective parameters. However, the precise evaluation of the e ﬀ ect of microstructural parameters such as the contact surface distribution on the global mechanical behaviour of the soil requires the use of more advanced modelling methods. The paper presents the results of Discrete Element Modelling of triaxial tests with the open source code Yade in which the real microstructural properties of biocemented soil computed on 3D X-ray micro- tomography images are used as input parameters. A particular attention has been paid to take into account the actual distribution of contact surface in the model and not only the average value. It appears that the model is then able to reproduce the evolution of the macroscopic properties (in particular that of the cohesion) with the calcite content.


Introduction
Bio-cementation focuses on the bio-transformation of sand to sandstone to improve its mechanical properties, using microorganisms such as urease bacteria. [1] and [2] are the first to have used the potential biocatalyst to cement soils. Since then, research in this field has developed in several ways. As a result, bio-mediated methods of strengthening granular soil through mineral precipitation, have emerged as a promising means of soil improvement [3,4]. Recently, the behaviour of the bio-cemented sand has been experimentally explored in several studies [5][6][7]. These experimental studies are being accompanied by some numerical studies.
To design an efficient calcite bio-mineralization scheme, it is important to be able to predict the improvement in soil properties resulting from calcite precipitation. It is quite clear that the mechanical efficiency of the bio-cementation process is directly related to the microstructural properties of the bio-cemented sand, such as, the volume fraction of calcite, its distribution within the pore space (localized at the contact between grains, over the grain surfaces) * e-mail: moha.abbas1996@gmail.com and the contacts properties: coordination number, contact surface area, contacts orientation, types of contact. [8] have used microscopic properties computed from 3D images obtained by X-ray tomography of bio-cemented sand as an input in current analytical models to estimate the elastic properties (Young's and shear modulus) and the strength properties (Coulomb cohesion). However, the effect of these microstructural parameters such as the contact surface distribution on the global mechanical behaviour of the soil remains hypothesized or vague and still needs some clarification using advanced modeling methods.
The objective of the present work is to contribute in explaining the bio-cemented sand behaviour by using the real microstructural properties of bio-cemented sand (such as contact surface area distribution...) computed on 3D X-ray micro-tomography images as an input in a numerical DEM triaxial test model to explore the influence of these parameters on the microscopic behavior of the bio-cemented sand.

Methodology for modeling bio-cemented sand
The proposed method is implemented on Yade [9], an extensible open-source framework for discrete numerical models, based on Discrete Element Method. A simulation run in Yade operates according to a classic DEM algorithm. The basic steps include the detection of collision between discrete elements (DE also referred to as particles), the computation of interaction forces between these DE based on constitutive laws, the computation of particles' accelerations based on Newton's second law, and the time integration of these accelerations for the determination of the particles' next positions. These steps are repeated until the simulation is over.

Constitutive Law
For the bio-cemented sand modeling, the constitutive law used is a cohesive law within Yade presented by [10]. The law is for linear traction, compression, bending and twisting, using both cohesive and frictional interactions, and defining failure criterion by Mohr-Coulomb plasticity surface, thus adding adhesion and moment to the traditional Cundall's linear elastic-plastic law that has three parameters: tangential stiffness, normal stiffness, and friction coefficient [9]. This law considers spring-type stiffnesses: considering two spheres each of an elastic modulus E and diameter D = 2R, the normal stiffness K n and the tangential stiffness K s are the defined by: where ν is the stiffness ratio which is a parameter of the model. Adhesions related to the normal and the shear components are calculated from the normal cohesion C n and the shear cohesion C s . These two values should be given as parameters. For two particles of radii R 1 and R 2 , the adhesion will be: a n = C n min(R 1 , R 2 ) 2 , a s = C s min(R 1 , R 2 ) 2

Simulation's parameters for the untreated sand
The untreated sand is cohesion-less, so the cohesion parameters C n and C s defined in 2.1 are considered null. According to the chosen constitutive law, the tangential stiffness K s and the normal stiffness K n are mainly dependent on elastic modulus E. As a result, the material parameters which control the behaviour are the elastic modulus E and the contact friction angle ϕ. By running several simulations of drained triaxial tests, we observed that E changes significantly the macroscopic elastic modulus of the sample. On the other hand, the contact friction angle ϕ affects only the peak stress and the residual stress. More precisely, it is proportional to both stresses, but still has much greater effect on the peak stress than the residual stress that reaches a certain maximum stress at a certain ϕ max .

Material parameters
The bio-cemented sand modeling is based on resembling the cementation effects by adding cohesion and applying certain constraints to the existing untreated sand model. The bio-cemented sand is made up of a group of particles connected to each other in a way that they can neither slide nor rotate as long as they are cemented to each other; which means that the friction has no effect anymore as long as the cementation bonds exist. The constitutive law chosen for this modeling in 2.1 is then adapted to describe such behaviour between particles. Mainly, the friction will be ignored as long as the adhesion is active on each contact, and when the maximum tensile or maximum shear force is reached, the cohesive link is broken and the behaviour will be purely friction afterward. The bond between particles created by this contact law has no weight (i.e. all the mass is concentrated in the particles and is therefore conserved when the bond breaks). In order to estimate the normal cohesion parameter, we are going to relate the adhesion to the cohesive contact surface areas that were computed on 3D X-ray microtomography in [7] for each sample. [11] have developed an expression of the tensile strength of cohesive material using the volume averaging technique inside Voronoi cells taking into account the capillary forces (suction). In order to describe our material, we propose to replace the capillary force by the bonding force induced by the contact cementation. This inter-particle force, i.e. normal adhesion, can be defined by: where S c is the cohesive contact surface area of a cemented bond and σ ten is its tensile strength, as introduced in [12]. Generally, in granular materials, the tensile strength depends on the intrinsic strength of the cement material and its adhesion at the grain surface. According to [8], the tensile strength σ ten was calibrated in order to reflect the experimental results obtained for the biocemented samples in which the cohesion was measured by drained triaxial tests on the same bio-cemented sand. This calibration yielded to σ ten = 2.75 MPa . Considering that we are studying the same samples, the value of the intrinsic tensile strength is then taken as σ ten = 2.75 MPa for all the simulations. From the inter-particle force and the normal adhesion formulas and by placing f n = a n , we can then determine the normal cohesion C n based on the formula: The shear cohesion parameter was determined as a function of the normal cohesion C n . To do so, the normal stiffness parameter S n and the shear stiffness parameter S τ were computed using the formulas established by [13]. They depend on the shear modulus and Poisson's ratio of both the grains and the cement layer in addition to the radius of the sand grains and the radius of the cement layer as defined and used in [14,15] . The results are shown in Table 1. The shear cohesion C s was then calculated as follow: Based on our calculations of this ratio for each subsample (Table 1), we can consider that the ratio S τ /S n is almost constant and equal to 2.4 for all the treated samples considering. Thus, the shear cohesion was then expressed as: Since a s /a n = C s /C n , the shear adhesion a s can be also expressed in terms of the normal adhesion a n using the same ratio S τ /S n as follow: a s = 2.4 a n

Cohesion contact area distribution and percentage of cohesive contacts
The cohesive surface area distributions were extracted from the 3D X-ray micro-tomography images by [7], in which the distribution of the contact area before and after cementation for each treated case were already computed based on these images [7]. Using these distributions, we then build up the distribution of cohesive contact area S c by subtracting the contact surface area before treatment from the contact surface area after treatment for each contact ( Figure 2). From the resulting distribution, the mean cohesive contact areaS c and the corresponding standard deviation σ std were calculated for each case. The best fitting is the log-normal distribution ( Figure 3) for all the cases which is then built up using the meanS c and the standard deviation σ std for each treated case. The values from the log-normal distribution, corresponding to the cohesive contact area S c , are then randomly spread on all the contacts in the DEM triaxial model for the fully cemented cases. In this way, the actual cohesive contact area distribution is maintained for these subsamples. However, in order to maintain the real cohesive  contact area distribution for partially cemented subsample, the percentage of cohesion should be taken into account, which means that the values extracted from the corresponding log-normal fitting can not be assigned directly to all the contacts. To do so, first a certain percentage of the contacts is assigned a zero cohesion, which reflects the frictional contacts remaining after the bio-cementation process, and then for the rest of the contacts which are cohesive contacts, values from the corresponding log-normal distribution representing the cohesive contact area S c are assigned in a random manner in a similar way to the fully cemented cases. Using this procedure, we will keep both the actual frictional to cohesive contacts ratio for the partially cemented subsamples and the actual cohesive contact area distribution for the cohesive contacts within these subsamples.

Setting up the simulations
The triaxial cylindrical mold used in the experiments in [7] was replaced by a 3D rectangular parallelepiped by applying preservation of volume and neglecting corner effects ( Figure 1). It was then scaled in a way to fit the simulations by using a scaling factor of 30. In addition, rigid walls boundary conditions were used in which the walls act as dynamic bodies. The sample contains exactly 1602 spherical particles having a mean diameter 210 µm ± 20%. The particles were generated using a random dense pack [9], which spreads them randomly inside the specified volume. The sample generated was then used for all the different simulations, that is to say that the sample used for each of the treated cases was exactly the same as the sample used   for the untreated simulation, having the same number of particles and the same radius and initial position for each particle.

Untreated sand simulation
The Young's modulus parameter E of the particles and the contact friction angle ϕ should be assigned taking into consideration their effect on the macroscopic behaviour as discussed in 2.2. The aim is to reproduce the experimental results for the untreated case by changing the two aforementioned parameters. The values assigned for these parameters are then used as base values for the bio-cemented sand simulations. In other words, these two parameters were used as tuning parameters. The experimental results for the untreated sample under a confining pressure of 50 KPa were extracted from [7], and after several tries, the best fit was reached using a Young's modulus parameter E = 80 MPa and a contact friction angle ϕ = 15 • . The simulations were then repeated using 25, 75 and 100 KPa confining pressures and the corresponding simulations results versus the experimental results along with the 50 KPa confining pressure case are plotted in Figure 4. From Figure 4, we can see that the general behaviour of the untreated sand was reproduced. The residual stress was almost the same for the different confining pressures, but the peak stress varied more from the experimental peak values as the confining pressure increased. In addition, the macroscopic elastic modulus obtained is almost the average value of the macroscopic elastic modulus experimentally obtained for different confining pressures which for instance varied for different confining pressures. As a conclusion, the base values of the Young's modulus parameter E and the contact friction angle ϕ which will be taken into account for the bio-cemented sand simulations are 80 MPa and 15 • respectively.

Fully cemented sand -High calcite content
For the fully cemented sample simulation, the data related to the subsample 13MB in [7] was considered. The two simulation's parameters, the Young's modulus parameter Axial strain(%) Experiment CP= 25KPa Simulation CP= 25KPa Experiment CP= 50KPa Simulation CP= 50KPa Experiment CP=100KPa Simulation CP=100KPa Figure 5. High calcite content subsample 13MB simulation: Deviatoric stress versus axial strain E and the contact friction angle ϕ, of the untreated sand sample discussed in 2.2 are part of the parameters to be assigned in the treated sand sample simulations. The treated sample can be seen as an untreated sand sample with certain calcite layers between its particles connecting them to each other, that is cementing the sample. Since the cementation process was done on the same untreated samples modeled in 3.2, the values of the parameters related to the untreated sample simulations can be considered as base values for our simulations in this case, i.e. values that we can not go below. For Young's modulus parameter E, a certain modification is certainly needed because obvioulsy the normal stiffness increases with the increase of percentage of calcite in the sample, which was then considered 1 GPa. For the contact friction angle ϕ, the base value was considered for all the treated samples that is ϕ = 15 • . This is because the particles of the treated sample are the same as the particles of the untreated sample. Although the precipitation of calcite on the particles' surfaces after treatment will increase the roughness of these surfaces, this increase in roughness will be neglected in our study. In Yade, cohesion is a material property not a contact property, but adhesion is a contact property. As a result, the normal adhesion a n was assigned based on the values of the cohesive contact area for each contact as stated in 2.3.2. Initially, the normal cohesion was assigned the mean cohesive value calculated according to the formula mentioned in 2.3.1, then the normal adhesion was assigned for each contact as stated in 2.3.2, a way used to assign different cohesion for the same material, that is different adhesion. By observing the results in Figure 5, we can see that the general behaviour of the high calcite content sand was reproduced. The residual stress was almost the same for different confining pressures, and the peak stress values were close to the experimental values. For instance, the numerical peak stress when the confining pressure was 100 KPa was less than experimental peak for the same case, while both the numerical and the experimental peak were almost the same when the confining pressure was 50 KPa, but the numerical peak was greater than the experimental peak when the confining pressure was 25 KPa. In other words, the peak stresses in the simulations are closer to each other Axial strain(%) Experiment CP= 25KPa Simulation CP= 25KPa Experiment CP= 50KPa Simulation CP= 50KPa Experiment CP=100KPa Simulation CP=100KPa Figure 6. Medium calcite content subsample 13TT simulation: Deviatoric stress versus axial strain than the experimental peak stresses. In addition, we can observe from the curve shape after the peak that the general behaviour is more ductile numerically, while in reality, it is more brittle. This difference can be partly attributed to the DEM model (for instance the considered boundary conditions) but can also be explained by the procedure used in the experiments. Still, we can say that the general behavior was reproduced for the high calcite content sample using the real cohesive contact area distribution. The same was observed after repeating the simulations using a medium calcite content sample ( Figure 6).

Partially cemented sand -Low calcite content
For the partially cemented sample simulation, the data related to the subsample 2T in [7] was considered. The Young's modulus parameter was taken 300 MPa while for the contact friction angle ϕ, the base value was taken that is 15 • . The normal adhesion a n was assigned based on the values of the cohesive contact area for each contact as stated in 2.3.2 for partially cemented sand. The percentage of cohesive contacts was varied. Several cases are considered: 25%, 50%, 75% & 100% of the contacts are cohesive. This is because there exist an important percentage of purely frictional contacts for the low calcite content samples according to [7]. We then studied how the behaviour of the low calcite content sample changed with the change of the cohesive contact fraction. Figure 7 shows more clearly the effect of the change of the percentage of cohesion. We can observe how the behaviour is changing from an almost purely frictional behaviour when the percentage of cohesion is 25% to an almost purely cohesive behaviour when the percentage of cohesion is 75%, to an obvious cohesive behaviour when the all the contacts are cohesive, i.e. percentage of cohesion is 100%. We can say that any sample having a certain percent of cohesion will be stiffer than a purely frictional sample. It becomes more and more stiff as the percentage of cohesion increases until that percent becomes 100% turning it into a purely cohesive sample.. It is important to mention that the contact friction angle ϕ which affects the peak stress and the residual stress according to 2.2 can play an important role in reproducing the exact experimental behaviour for the partially cemented sand.

Mechanical Properties
Some of the mechanical properties of the bio-cemented sand were then studied; the peak friction angle, the residual friction angle and Coulomb cohesion which were alll determined using the Mohr's circle plots. The simulation results are listed in Table 2 while the experimental results extracted from [7] are listed in Table 3. We can observe that the peak friction angle and the residual friction angle increase as the percentage of calcite by mass increases which is in accordance with the experimental results in [7]. Nevertheless, it appears that the difference between peak and residual angles is smaller in the numerical simulation than in the experiments while the experimental residual angle is satisfyingly reproduced by the DEM model. The increase in the peak friction angle is explained by the increasing strength of each bond with the cementa- tion level. The evolution of the residual friction angle can be explained by the failure of distributed contact bonds that induce the creation of clusters of particles (with an irregular and therefore more frictional geometry). Groups of bonds connecting the particles to each other are broken through out the triaxial samples, still another group of bonds is not broken and after certain strain, the sample reaches a stable state when the unbroken bonds that are still connecting some particles are no more breaking, as if the particles that are still connected to each other form one new particle of greater size (Figure 1). The cohesion also increases importantly as the percentage of calcite by mass increases which is in accordance with the experimental results too. However, unlike both the peak friction angle and the residual angle results, the numerical results are greater than the experimental results. This is because the peak stresses are closer to each other numerically than they are experimentally, which leads to the decrease of the slope of the plotted tangents and thus the increase of the vertical axis intersection value. Still in general, the results can be considered similar to the experimental results.

Summary and Conclusion
The objective of this paper is the numerical modeling and investigation of the bio-cemented sand behaviour by using the real microstructural properties of bio-cemented sand which were computed on 3D X-ray micro-tomography images (distribution of the contact surface area, type of contacts). These properties were used as input parameters in Yade, an open source D.E.M. platform, in order to explore the influence of these parameters on the microscopic behavior of the bio-cemented sand. The main focus was to take into account the real distribution of contact surface in the model instead of the average value, and to study the effect of mixed cohesive and frictional contacts in the macrostructural behaviour. A cohesive frictional constitutive law for contacts was chosen for the simulations then the cohesive model was built up by relating the cohesion parameters in the model to the cohesive contact areas, then using a log-normal fit for the distribution so to assign randomly the values to all contacts. Untreated simulations were then performed in order to obtain base values for the Young's modulus parameter and the contact friction angle. These values were then adapted and used in the fully cemented samples simulations and the results were similar to the experimental results, still more ductile numerically. It appears that the model was then able to reproduce the evolution of the macroscopic properties with the calcite content. The simulations were then performed on partially cemented sample using different percentages of cohesive contacts. The results showed that the higher is the percentage of cohesion contacts, the stiffer is the sample and the more cohesivelike is the post-peak behaviour. These simulations made us understand the evolution of behaviour of the partially cemented sand. The peak and residual friction angles were then calculated, and a noticeable increase in both values was observed as the percentage of calcite increases in the sample in accordance with the experimental result. This can be explained by the formation of groups of particles that act as new particles of bigger sizes within the sample. Coulomb's cohesion also showed an important increase as the percentage of calcite increases which was also in accordance with the experimental results.