Automating the Calibration of Flat-Jointed Bonded Particle Model Microproperties for the Rewan Sandstone Case Study

The flat joint contact model (FJM) provides significant improvements over its predecessors, the parallel bond and contact bond models, for bonded particle modelling of rocks due to its unique microstructure that allows for the reproduction of the macroscopic compressive–tensile strength ratio, σc/|σt|; internal friction angle, ϕ; and the Hoek–Brown constant, mi. However, the microproperty calibration process is tedious and time-consuming to perform manually due to the various microproperty interdependencies that exist in the FJM. Previous attempts at automating the bonded particle model microproperty calibration process have typically utilised advanced statistical methods, such as artificial neural networks, but they have not yet been widely applied to the FJM over a representative range of confining stresses for calculation of ϕ and mi. In this study, a new method is proposed for automating the FJM microproperty calibration process based on a numerical root-finding algorithm and specific calibration sequencing. The new method is applied to a Rewan Sandstone case study with similar natural porosity to a 2D bonded particle model packed to a low initial mean stress. The resulting FJM microproperties are shown to reproduce both the target macroscopic laboratory properties and a realistic damage evolution, including a normalised crack initiation stress of 0.46 and a normalised crack damage stress of 0.83 coinciding with a reversal of the axial stress–volumetric strain curve in an unconfined compression test simulation. It is also demonstrated that the absolute change in the instantaneous lateral–axial strain ratio (Poisson’s ratio in the linear-elastic phase) provides a reasonable proxy to the acoustic emissions which may be measured in the laboratory. Successful convergence of a novel automated microproperty calibration method for a Rewan Sandstone flat-jointed bonded particle model case study. Reproduction of realistic damage evolution, including crack initiation and crack damage stress thresholds. Microstructural properties including number of flat joint radial elements and installation gap control the macroscopic compressive-tensile strength ratio and dilatancy of the bonded particle model. Successful convergence of a novel automated microproperty calibration method for a Rewan Sandstone flat-jointed bonded particle model case study. Reproduction of realistic damage evolution, including crack initiation and crack damage stress thresholds. Microstructural properties including number of flat joint radial elements and installation gap control the macroscopic compressive-tensile strength ratio and dilatancy of the bonded particle model.

• Successful convergence of a novel automated microproperty calibration method for a Rewan Sandstone flatjointed bonded particle model case study. • Reproduction of realistic damage evolution, including crack initiation and crack damage stress thresholds.
• Microstructural properties including number of flat joint radial elements and installation gap control the macroscopic compressive-tensile strength ratio and dilatancy of the bonded particle model. Macroscopic discontinuity friction angle measured from direct shear tests c shr /c all Shear crack proportion at peak stress B CN Bond coordination number σ 3max Maximum confining stress for triaxial tests σ 3n Normalised confining stress for triaxial test stage σ 3 Minor principal stress σ 1 Major principal stress n Number of triaxial stages m i Hoek-Brown constant s Hoek-Brown constant a Hoek Flat joint contact radius λ Flat joint contact model radius multiplier W Bonded particle model assembly width H Bonded particle model assembly height R Model resolution; average number of particles across minimum particle assembly dimension D min Minimum particle diameter D max Maximum particle diameter D avg Average particle diameter D ratio Particle diameter ratio, D max /D min σ 0 Target mean stress for grain-scaling procedure ė QS Quasi-static strain rate ė step Step-strain rate t crit Critical timestep of a bonded particle model assembly v̇Q S Quasi-static loading velocity h 0 Initial bonded particle model specimen dimension in major principal stress axis [a i ,b  Microproperty-macroproperty ratio in iteration i

Introduction
The bonded particle model (BPM) is a subset of the distinct element method (DEM) that represents an intact rock matrix as an assembly of round particles that may independently translate and rotate and which interact at their contact points via a microscopic constitutive law (Potyondy and Cundall 2004). Prior to the application of boundary forces, contacts are initially bonded with finite tensile and shear strengths to represent the cohesion provided by cement and granular interlock in a real rock specimen. In Itasca Consulting Group's Particle Flow Code (PFC), the flat joint contact model (FJM) is the preferred microscopic constitutive law for bonded particle modelling as its unique representation of a contact as a series of elements that may fail independently allows matching of three key macroscopic parameters that its predecessors, the contact bond contact model (CBM) and parallel bond contact model (PBM), could not: (1) the compressive-tensile strength ratio, σ c /|σ t |; (2) the internal friction angle, ϕ; and (3) the Hoek-Brown constant, m i (Potyondy 2012(Potyondy , 2018Wu and Xu 2016). As noted by Cundall (2001), the DEM provides something that traditional continuum methods cannot: a means of modelling the complexity of heterogeneous and anisotropic cohesive-frictional materials without the need to invent macroscopic constitutive laws that often involve many obscure parameters and assumptions. In addition to its use as a tool for studying the failure mechanisms of intact rock, the BPM serves as the base material for the synthetic rock mass (SRM) method, which involves superimposing a discrete fracture network onto the BPM to explicitly simulate a fractured rock mass (Mas Ivars et al. 2011;Pierce et al. 2007). This is particularly useful as a supplementary exercise in the rock mass classification process in which traditional empirically derived classification systems such as the Geological Strength Index (GSI) operate under the assumption of isotropic and homogeneous conditions (Hoek 1994;Hoek et al. 2013Hoek et al. , 1995Hoek et al. , 1992. To date, the BPM and SRM methods have collectively been applied to various rock mechanics applications, including but not limited to: seismicity (Hazzard et al. 2002(Hazzard et al. , 2000Hazzard and Damjanac 2013;Hazzard and Young 2004;Khazaei et al. 2016;Zhang et al. 2017); hydraulic fracturing (Duan et al. 2018;Yoon et al. 2014Yoon et al. , 2015Zhao and Young 2011); scale effects and anisotropy (Bahaaddini et al. 2016;Cundall et al. 2008;Deisman et al. 2010;Esmaieli et al. 2010;Farahmand et al. 2018;Yang et al. 2016;Zhang et al. 2011); true triaxial strength (Duan et al. 2017;Mehranpour and Kulatilake 2016;Schöpfer et al. 2013;Zhang et al. 2019); slope stability (Huang et al. 2015;Kim and Yang 2005;Wang et al. 2003); and underground excavations (Gao et al. 2014;Hadjigeorgiou et al. 2009;Zhang and Stead 2014). Despite this wide range of demonstrated end uses, all such models share a time-consuming preliminary step: calibration of the contact model mechanical properties (microproperties).
As it is generally not possible to measure microproperties directly in the laboratory, calibration is treated as an inverse modelling problem in which the microproperties are iteratively updated until the macroproperties measured from simulated laboratory tests match the corresponding laboratory properties of physical specimens (Chen 2017;Potyondy and Cundall 2004;Shu et al. 2019;Vallejos et al. 2016;Wu and Xu 2016). This is a tedious and time-consuming task to perform manually as it requires potentially hundreds of iterations, each of which requires several minutes to hours depending on the model resolution and contact stiffnesses. While the BPM is a versatile, multi-purpose tool that has the potential to address many long-standing rock mechanics problems, it is computationally expensive and would therefore benefit from the development of an objective, efficient, and repeatable means of automating the microproperty calibration process. Previous attempts at automating the microproperty calibration process have involved the use of advanced statistical methods, such as artificial neural networks (Sun et al. 2013;Tawadrous et al. 2009;Yoon 2007). However, such methods are challenging to automate and have not yet been widely applied to the FJM over a representative range of confining stresses for calculation of ϕ and m i . In this study, we propose a simple trial-and-error method of calibrating the FJM microproperties based on the bisection search method (BSM) numerical root-finding algorithm and appropriate sequencing of the calibration process.

Kinematics
In the general DEM, particles are represented as rigid discs in 2D or spheres in 3D that may independently translate and rotate, interact with other particles at their contact points, and overlap between model cycles for force-displacement updates. As described by Cundall and Strack (1979), particle motion in the DEM is based on Newton's laws of motion and the model state is advanced by a time-centred stepping algorithm. At each model step, a force-displacement law is Overlap of displaced particles resulting in non-zero interparticle force at contact point B at time t 2 = t 0 + 2Δt. Modified from Cundall and Strack 1979 utilised to calculate forces and update velocities resulting from interparticle overlap whereby the force is proportional to the overlap. The force is then converted to stress based on the local contact length in 2D or area in 3D and the resultant stress is then compared with the bond tensile and shear strengths to determine the bond failure state. The implementation of the time-centred stepping algorithm with respect to motion, contact overlap, and updating of displacements and forces is illustrated in Fig. 1. A full description of the model sequence and kinematics is available in the PFC user manual (Itasca Consulting Group 2018).

Micromechanical Properties
The constitutive behaviour of the general DEM is governed by a contact model that locally defines the mechanical properties of each interbody contact. In the BPM, the contact model has an initial cohesive bonded state and a residual frictional unbonded state that is activated when the bond is broken. The FJM is the preferred BPM contact model in PFC and comprises the following microproperties shown in Fig. 2: • Normal-shear stiffness ratio, k n /k s ; • Effective modulus, E c ; • Mean and standard deviation tensile strength, σ b,m and σ b,sd ; • Mean and standard deviation cohesion, c b,m and c b,sd ; and • Bond and grain friction angles, ϕ b and ϕ g respectively, the latter of which is represented in PFC as a friction coefficient, μ, where μ = tan(ϕ g ).
In the bonded state, both strength (σ b,m , σ b,sd , c b,m , c b,sd and ϕ b ) and deformability (k n /k s , E c ) microproperties are active, while only k n /k s , E c , and ϕ g are active in the unbonded state. Note that, to reduce the number of free variables in the system, the same deformability properties are generally specified in the bonded and unbonded states.

Microstructural Properties
In addition to the micromechanical properties, the FJM includes several microstructural properties that enable its unique ability to match σ c /|σ t |, ϕ, and m i . These include the following: • The number of contact radial elements, N r (2D and 3D), and number of contact circumferential elements, N c (3D only); • The bond installation gap ratio, g i /D avg , which defines the particle proximity for contact detection during bond installation; • The gapped and slit proportions, F g and F s , which simulate the initial microcrack conditions; and • The initial surface gap, g 0 , of gapped contacts.
As illustrated in Fig. 3, the bond coordination number, B CN , which is the average number of contacts per particle in the assembly, can be increased from a maximum of 3-4 in a contact-bonded or parallel-bonded BPM to more realistic values of 6-10 for rocks (Oda 1977) by segmenting each contact into multiple elements that may fail independently of each other via N r and N c . Therefore, the influence of the loss of a single element on the residual rolling resistance of a particle in a flat-jointed BPM is less than in a binary contact model such as the CBM or PBM in which the bond is either on or off with no capacity for partial damage. A flat-jointed BPM can therefore sustain greater stress beyond the tensile-controlled crack initiation stress 1 , σ ci , and, subsequently, σ c /|σ t |, ϕ, and m i can be matched. To prevent flat-jointed contacts from overlapping and thus ensure microstructural validity, the radius of each contact is multiplied by a factor, λ, where 0 ≤ λ ≤ 0.58 (Patel 2018). Patel and Martin (2020) showed that setting F g and g 0 > 0 allows the FJM to reproduce the non-linearity of the stress-strain curve prior to the macroscopic crack closure stress, σ cc , which marks the beginning of the linear-elastic phase after the closure of pre-existing microcracks and pore spaces. However, attempting to match σ cc introduces additional complexity and variability to the model. For the purpose of general microproperty calibration where the spatial positions and nature of the initial microcracks affecting σ cc are unknown, it is practical to set F s and F g to 0, with target macroproperties instead scaled to 80% of their laboratory values after Martin et al. (2011) to implicitly simulate the macroscopic influence of the microcracks. Further discussion on the role of pre-existing microcracks represented by the gapped and slit contacts is presented in Sect. 7.2.2 1 .

Fig. 3
Comparison of grain connectivity and bond coordination numbers for: a A real rock with tightly interlocked, angular grains. b A parallel-bonded particle model with a bond installation gap of 0. c A parallel-bonded particle model with a non-zero bond installation gap. d A flat-jointed bonded particle model with a non-zero bond installation gap and two radial elements. Modified from Potyondy (2018) 1 3

Model Generation
The first step in the simulation of laboratory tests with a BPM is the generation and packing of a frictionless, unbonded particle assembly to a target mean stress, σ 0 , which is analogous to the locked-in stresses that may be present in a real rock specimen. If the locked-in stresses are unknown, σ 0 is typically set to 100 kPa to approximate Earth's atmospheric pressure. The key geometric parameters of the particle assembly that must be specified for the packing procedure are the following: • Width, W, and height, H, of the specimen/material vessel; • Model resolution, R, which is the average number of particles across the minimum specimen dimension; • Minimum, maximum, and average particle diameters, D min , D max , and D avg ; and • Maximum-minimum particle diameter ratio, D ratio .
D avg can be calculated as per Eq. (1) where R ≥ 20 after Ding et al. (2014). For uniformly distributed particle diameters, D min and D max can then be calculated as per Eqs. (2, 3) and do not need to be specified directly. In this case, only W, H, R, and D ratio need to be specified As BPMs are typically simulated with particles at the microscopic scale instead of the true mesoscopic scale of rock grains and cement, they are not a direct 1:1 analogue of real rock specimens and matching the macroproperties is therefore a higher priority than explicit reproduction of the grain matrix. For this reason, so long as a microstructurally valid particle assembly is used, it is recommended to select a constant set of geometric and microstructural properties. Calibration should then focus on mathematical convergence of the micromechanical properties and reproduction of a realistic damage evolution identified by σ ci and the crack damage stress, σ cd , coinciding with a pre-peak volumetric strain reversal. The damage evolution is discussed further in Sect. 5.2.
Once the geometric and microstructural parameters of the particle assembly have been selected, packing is performed according to the grain-scaling procedure of Potyondy (2017) which is illustrated in Fig. 4 and involves the following: 1) Randomly placing frictionless particle seeds within a closed material vessel and then upscaling their radii to achieve large initial overlaps. 2) Allowing the assembly to re-arrange within the material vessel to an initially high mean stress and low porosity. 3) Iteratively reducing the radii of all particles by a small amount and cycling the model until the mean stress of the assembly is within a specified tolerance of σ 0 .
At the end of the packing procedure, g i is increased until the specified target B CN is matched. Finally, the contact model is applied and the specimen is saved in two separate states: within the material vessel for compression test simulations and without the material vessel for direct tension test simulations.

Boundary Conditions
The frictionless material vessel walls used for specimen generation and packing act as load platens in compression tests with axial strain induced by converging the walls in the axis of loading at the quasi-static strain rate, ė QS , which is calculated according to Eq. (4) where: • ė step is the step-strain rate (strain per model step) at the quasi-static strain rate, typically set to 1.1e −8 after Zhang and Wong (2014); and • t crit is the global critical timestep of the particle assembly In PFC, t crit is automatically determined during model cycling and is therefore not a free variable. All wall-particle contacts are assigned a simple zero-strength linear contact model with E c set slightly higher (~ 10 to 20%) than the material and k n /k s is equal to 1. The quasi-static loading velocity, v̇Q S , applied to each wall is directly proportional to the quasi-static strain rate and specimen dimension in the axis of loading, h 0 , as per Eq. (5) In confined compression tests, lateral confinement is controlled by a servomechanism that continuously monitors and adjusts the sidewall velocities to maintain the target confining stress. While cumbersome to perform in the laboratory, direct tension tests are more easily simulated numerically and are preferable to Brazilian tensile tests as they are less sensitive to model resolution effects. For the purpose of acquiring a target uniaxial tensile strength, the Brazilian indirect tensile strength can be converted to an equivalent uniaxial tensile strength via an empirical conversion factor of 0.7 after Perras and Diederichs (2014). Axial strain is induced in direct tension tests by assigning two rows of grip particles at opposing ends of the specimen and diverging them at the same loading velocity as the compression tests for direct comparability of results. The general laboratory test simulation configurations are illustrated in Fig. 5. Also shown are the measurement regions which are used to monitor stress-strain quantities in direct tension and unconfined compression tests.

Macroproperty Measurement
The key macroproperties for microproperty calibration are as follows: • Poisson's ratio, ν; • Young's modulus, E; • Uniaxial tensile strength, |σ t |; The ϕ fitted to a minimum of five confined compression test data points using the equations of Hoek and Brown (1997) and Hoek et al. (2002) is used for direct compatibility v QS Load platen Lateral wall servomechanism   Hoek and Brown (1997) and Hoek et al. (2002): For n ≥ 5; s = 1; a = 0.5; and σ 3max equal to 25% of σ c with the bond friction angle, ϕ b . To enable automation of the iterative calibration process, it is necessary to calculate macroproperties from measured stress-strain quantities during model cycling and not with external post-processing tools. In compression tests, stress-strain quantities are measured from wall displacements and force responses. However, as there are no lateral walls in unconfined compression tests and neither lateral nor axial walls in direct tension tests, it is necessary to also measure quantities from measurement regions installed within the specimen. Measurement regions average the stress-strain quantities measured at the contact points of all particles with centroids inside the region and are therefore analogous to internal stress-strain probes as previously shown in Fig. 5. While |σ t | and σ c are simply identified as the maximum stress quantities in their respective laboratory test simulations, ν, E, and ϕ are calculated according to the methods summarised in Table 1. For laboratory specimens, secant values of ν and E are typically calculated over the range from 0 to 50% of the peak stress in an unconfined compression test (Brady and Brown 2006), and this range is also used for BPM microproperty calibration. Further discussion on the measurement of E and ν with respect to damage thresholds is provided in Sect. 6.1.

Application of the Bisection Search Method to Microproperty Calibration
The BSM, also known as the binary search method (Burden and Faires 2011), is a numerical root-finding algorithm that iteratively converges a search range, [a i ,b i ], on a target, T, by: (1) taking the output, c i , as the average of [a i ,b i ]; (2) identifying whether c i is less than or greater than T; and (3) discarding the unused portion of the search range and repeating steps 1 and 2 until c i is within tolerance, ε, of T. This general convergence process is illustrated in Fig. 6. While the BSM is the slowest of the numerical root-finding algorithms, it is also the most reliable, simplest to implement numerically, and is guaranteed to converge on a solution provided one exists within the initial search bounds.
With respect to calibration of BPM microproperties, the BSM is performed on the ratio of the microproperty, m i , to its corresponding macroproperty, M i , rather than absolute values of m i , to ensure that the algorithm is generic Fig. 6 Convergence of a search range on a target using the bisection search method for a range of rock types of variable strength and stiffness. Moreover, it provides a means of directly comparing the proportionality of calibrated microproperties for different rock types. Therefore, [a i ,b i ] and c i are, respectively, the search range and bisected value of the microproperty-macroproperty ratio, R i , while T is the value of R i that produces the target value of M i in the relevant laboratory test simulation. The generic application of the BSM to the calibration of a particular microproperty is illustrated in Fig. 7.
A limitation of the BSM is that the initial search range, [a 0 ,b 0 ], must be specified. This requires initial parameters to be assumed, which is not ideal for the development of an objective and generic microproperty calibration algorithm. To overcome this problem, the first two iterations of the calibration sequence are used to establish initial error bounds and estimate [a 0 ,b 0 ] for initiation of the BSM in iteration 3. In iteration 1, R i is set equal to 1, such that the microproperty is directly proportional to the corresponding target macroproperty ( Table 2). The exception is k n /k s , which is directly set equal to 1. In iteration 2, a linear microproperty-macroproperty relationship is temporarily assumed and R 2 is set equal to R 1 multiplied by the resulting gradient as per Eq. (6) For E c and σ b,m , this often results in calibration being achieved in the second iteration due to the approximately linear microproperty-macroproperty relationships. However, the non-linear relationships between c b,m /σ b,m, k n /k s , c b,m , and ϕ b and their corresponding macroproperties generally require more than two iterations to converge. For iteration 3, which is the first iteration in which the BSM is invoked, where e 2 is the error resulting from iteration 2. The BSM convergence process is then applied from iteration 4 until calibration is achieved. Fig. 7 Application of the bisection search method to microproperty calibration

Calibration Sequence
The global calibration is performed as a sequence of local calibrations in series and parameter interdependencies are therefore not directly addressed within a BSM sequence for a particular microproperty-macroproperty pair. This is in contrast to global optimization methods that test a wide array of microproperty combinations and return the combination that produces the minimum error as the calibrated set of microproperties. While global optimization methods directly address parameter interdependencies, they require a significantly greater number of iterations than the proposed BSM-based procedure as they inherently include more fringe cases. The calibration sequence presented in Fig. 8 was informed by the sensitivity analyses of Wu and Xu (2016), Vallejos et al. (2016), and Chen (2017) and was designed to minimise the influence of the various microproperty interdependencies. Specifically, the microproperty controlling the proportion of shear and tensile strengths and hence the macroscopic failure mechanism, c b,m /σ b,m , is calibrated first for an arbitrarily low value of σ b,m as it affects the calibrated values of all other microproperties. Next, the microproperties controlling the deformability, k n /k s and E c , are calibrated. As k n /k s significantly affects both ν and E, while E c only affects E, k n /k s is calibrated before E c . Finally, the microproperties controlling the strength, σ b,m and c b,m , are calibrated separately. As σ b,m controls the initial extension strain-dominated portion of the stress-strain curve and hence σ ci , it is calibrated before c b,m which primarily affects the shear-controlled pre-peak plastic portion of the stress-strain curve from σ ci through to σ cd . Although c b,m /σ b,m is already calibrated at the beginning of the sequence, changes to other microproperties can slightly affect σ c /|σ t | and it is therefore necessary to adjust c b,m again at the end of the sequence. As ϕ b affects all macroproperties, it is performed as an outer loop within which all other microproperties are recalibrated in each iteration. It is noted that Vallejos et al. (2016) recommend matching ϕ to ϕ b , while Chen (2017) reports that ϕ is Fig. 8 Recommended sequence for microproperty calibration more sensitive to ϕ g . In the proposed sequence, ϕ g is treated as a fixed parameter equal to the discontinuity friction angle, ϕ j , measured from laboratory direct shear tests, and ϕ is a variable parameter that is matched to ϕ b .

Calibration Results
The proposed automated calibration algorithm was applied to the calibration of FJM microproperties in PFC for a Rewan Sandstone material for which the target macroproperties were derived from the median values of a large rock testing database consisting of 997 laboratory test results (Tsang et al. 2018). Material generation was carried out according to the grain-scaling procedure and all microproperties were calibrated to a tolerance of 3%. The packing parameters for the material generation procedure including the microstructural properties are summarised in Table 3. Although a D ratio of 1.66 is commonly used in the BPM literature based on the Lac du Bonnet Granite BPM of Potyondy and Cundall (2004), the D ratio of 2.0 reported by Sharrock et al. (2009) for a BPM representing the Hawkesbury Sandstone was adopted for the Rewan Sandstone case study. The target and calibrated macroproperties for the Rewan Sandstone are summarised in Table 4, the calibrated microproperties and microproperty-macroproperty ratios are summarised in Table 5, and the calibration record is presented in Fig. 9.

Damage Evolution
A unique feature of the BPM is the ability to visualise microscopic compressive and tensile forces at particle contacts. When these forces overcome the bond strength in either tension or shear, the bond fails and represents a new microcrack (Fig. 10). The progressive failure of bonds in a BPM is analogous to the damage evolution of laboratory specimens and allows for the identification of two key pre-peak damage thresholds: the crack initiation stress, σ ci , and the crack damage stress, σ cd , which Diederichs (2003) suggests can be taken, respectively, as the lower and upper bound operational strengths. Typically occurring at around 42-47% of the peak stress in laboratory UCS tests, σ ci marks the end of the linear-elastic phase and the onset of irrecoverable plastic strain (Nicksiar and Martin 2013). σ cd occurs later at around 78-90% of peak stress (Pepe et al. 2018;Xue et al. 2014), marking the end of the stable cracking process and onset of unstable cracking. That is, if the applied load is held constant at σ cd , cracking will continue and the specimen will eventually fail. Any gain in applied load from σ cd through to the peak stress is therefore time-dependent as a function of the boundary conditions and σ cd can accordingly be taken as the long-term yield stress (Hudson et al. 1972;Martin 1997;Martin and Chandler 1994).

Fig. 9
Calibration record for the Rewan Sandstone case study bonded particle model with 3% acceptance tolerance 1 3 In the laboratory, pre-peak damage thresholds can be identified from volumetric strain and acoustic emission (AE) monitoring (Brace et al. 1966;Eberhardt et al. 1998;Lockner 1993;Nicksiar and Martin 2012). In a BPM, tensile and shear crack events can be tracked directly without the need for complex AE waveform analysis. Additionally, as demonstrated by Diederichs et al. (2004), the proportion of the lateral and axial strains, ε ratio , which is Poisson's ratio if taken in the linear-elastic phase, can also be used as a damage indicator. This is because the rate of lateral dilation increases relative to the rate of axial contraction with increasing plastic strain and the associated reductions in axial and volumetric stiffness. As the volumetric strain reversal represents the point at which the rate of extensile lateral dilation overtakes the rate of compressive axial contraction and hence the onset of unstable cracking, there is a corresponding increase in the magnitude and volatility of the rate of change of ε ratio for which the absolute change, |Δε ratio |, can be used as a proxy to AE events.
The damage evolution of the calibrated Rewan Sandstone case study BPM is presented in Fig. 11 with respect to the stress-strain and crack-strain curves. Brittle behaviour is observed in which predominantly axially aligned tensile cracks progressively increase in intensity from σ ci through to σ cd beyond which inclined shear cracks initiate to link the existing tensile cracks and form a macroscopic axial splitting failure mechanism. The damage thresholds, σ ci and σ cd , occur at 46% and 83% of the peak stress, respectively, where σ ci was taken at the first major perturbation in the plot of |Δε ratio | rather than the first crack event, i.e., the onset of systematic cracking; and σ cd was taken at the volumetric strain reversal. An onset of micro-shear cracking coinciding with σ cd is observed, progressing to a final shear crack proportion at peak stress, c shr /c all , of 7.6%. These results demonstrate that the automated calibration procedure produces not just matching of target macroproperties but also a realistic damage evolution for the Rewan Sandstone case study. However, it is noted that the typical porosity of Rewan Sandstone is 10-15%, which is similar to the porosity of 16% for the 2D BPM packed by the grain-scaling procedure with a low initial target mean stress of 100 kPa. Further research is required to determine whether the BPM is suitable for materials with varying porosities or whether alternative DEM models such as the clumped model (Cho et al. 2007), grainbased model (Bahrani and Kaiser 2016;Lan et al. 2010;Saadat and Taheri 2020;Wong and Zhang 2018), or bonded block model (Purvance and Garza-Cruz 2020; Sinha and Walton 2020; Turichshev and Hadjigeorgiou 2017) produce more realistic damage evolutions for these materials.

Measurement of Deformability Properties
By definition, the true Young's modulus and Poisson's ratio are the tangent properties measured in the linearelastic phase of an unconfined compression test. For rocks, measuring the secant properties over the range from 0 to 50% of peak stress is a practical compromise when the crack initiation stress, and hence, the end of the linearelastic phase is unknown. For typical rocks in which the crack initiation stress is around 42-47% of peak stress, this means that 3-8% of the plastic phase in which the axial stress-strain rate is decelerating due to stiffness loss is included in the measurement of deformability properties. This is compounded by the stiffness gain in the preelastic phase relating to crack closure and the net effect of the two phenomena is that the secant properties can be significantly less than the tangent properties (Fig. 12). In a BPM with F g = 0, the crack closure phenomenon is not reproduced. The initial acceleration in the axial stress-strain rate relating to pore closure and stiffness gain does not exist and the resulting secant line is steeper than Fig. 11 Stress-strain and damage evolution for the Rewan Sandstone case study bonded particle model in a simulated UCS test it otherwise would be. Therefore, where possible, it may be prudent to use the tangent deformability properties as the target macroproperties for calibration of deformability microproperties if the crack closure phenomenon is not reproduced. The tangent deformability properties are also immune to interdependencies between c b,m /σ b,m , ϕ b , and ϕ g that can prolong the unstable cracking phase and peak stress, resulting in a reduced normalised crack damage stress, σ cd /σ c , and hence more of the plastic phase included in the secant deformability property calculation range.

Influence of Microstructural Properties
In the 2D flat-jointed BPM, the key microstructural properties affecting the mechanical behaviour are g i /D avg and N r , which control B CN , and F s and F g , which control the initial microcrack intensity. In the proposed automated BSMbased calibration scheme, the calibration is undertaken with respect to the micromechanical properties (ϕ b , k n /k s , E c , σ b,m , c b,m , c b,m /σ b,m ) for a constant set of microstructural properties. However, the selected set of microstructural properties (g i /D avg , N r , F s , F g ) must be capable of reproducing realistic rock behaviour as a precondition for the calibration process. Discussion and guidance on the selection of microstructural properties are thus provided herein.

Influence of the Bond Coordination Number
Parameter sensitivity analyses by Scholtès and Donzé (2013) and Wu and Xu (2016) revealed that σ c /|σ t | increases linearly with increasing B CN . This implies that, for a particular rock type with a particular characteristic σ c /|σ t |, there is a maximum B CN beyond which the minimum σ c /|σ t | of the BPM assembly is higher than the characteristic σ c /|σ t |, such that the microproperties cannot be matched irrespective of the calibration scheme. This affects the damage thresholds, σ ci and σ cd , and the measurement of the secant deformability properties. This in turn implies that there are maximum values of g i /D avg and N r for a particular rock type. For the Rewan Sandstone case study, the influence of B CN on the microproperties returned by the BSM-based automated calibration scheme was studied by varying g i /D avg and N r in isolation and re-running the calibration procedure for an acceptance tolerance of 3%. The results of the sensitivity analysis are presented in Fig. 13 and demonstrate the following: • With increasing N r for a constant g i /D avg of 0, the calibrated k n /k s and σ b,m increase slightly; c b,m and c b,m /σ b,m increase significantly; E c decreases slightly; and ϕ b decreases moderately. • With increasing g i /D avg for a constant N r of 2, the calibrated k n /k s increases significantly, while all other micro- properties experience only slight increases or decreases. It is noted, however, that for the case of g i /D avg = 0.05, c b,m /σ b,m first increased before decreasing at higher values of g i /D avg . This range coincides with the greatest increase in the B CN with diminishing returns observed thereafter as there are fewer unbonded contacts remaining. • For the combination of N r = 2 and g i /D avg = 0.20, the calibration procedure could not converge due to an excessively high minimum σ c /|σ t |. N r = 2 and g i /D avg = 0.18 corresponding to a target B CN of 10 therefore represent an upper bound for the Rewan Sandstone material.
These results indicate that g i /D avg primarily controls the stiffness ratio and hence the dilatancy of the specimen, whereas N r primarily controls the compressive-tensile strength ratio due to the reducing influence of the loss of each bonded element on the residual moment resistance of a particle with increasing N r . While g i /D avg and N r both increase the B CN , they do so in different ways. g i /D avg increases the B CN by extending the effective interaction range of proximal particles for bond installation (Fig. 14a). With respect to the particle geometry, g i /D avg increases B CN in the radial direction resulting in an overconnected, more interlocked particle assembly. Increasing g i /D avg can therefore be considered as analogous to increasing grain angularity. Although g i /D avg was in use in earlier parallel-bonded BPMs, it could not sufficiently increase the B CN to represent rock. This led to the introduction of the FJM by Potyondy (2012) to partition the overconnected contacts into multiple elements that can fail independently of each other via N r (and N c in the 3D BPM). Increasing N r therefore increases the B CN tangentially with respect to the particle geometry (Fig. 14b). While the partitioning of contact elements via N r achieves numerical matching of the B CN , it is noted that in a mechanical sense, there are diminishing returns with increasing N r as the area, and hence, load-bearing capacity of each element is proportional to (1/N r )λr c . Moreover, the partitioning of pre-existing contacts into multiple elements does not address the problem of the spatial distribution of intergranular contacts around the grain boundary in a real Implicit increase of B CN by partitioning pre-existing contacts into multiple elements rock assembly such as that shown in Fig. 3a. The numerical matching of B CN in the FJM is therefore implicit. Where grain angularity and interlocking are known to contribute significantly to macroscopic strength-typically in more brittle rocks of low porosity and high σ c /|σ t |-the use of non-uniform particle shapes via a clumped, grain-based, or bonded block model may be more suitable than the flatjointed BPM. For moderately porous rocks such as the Rewan Sandstone where macroscopic strength is instead controlled primarily by interstitial cement, the implicit representation of B CN in the FJM is generally adequate.

Influence of Slit and Gapped Contacts
Initial microcracks in rock reduce the macroscopic stiffness, such that less stress is transferred to the specimen per increment of strain and, subsequently, lower peak strengths are achieved (Martin and Stimpson 1994). In the FJM, initial microcracks are represented by slit and gapped contacts with their proportions specified by F s and F g , respectively. Slit contacts are unbonded contacts with g 0 = 0 and are mechanically equivalent to closed microcracks, such that friction is mobilised under compressive loading, but no crack closure effect is observed in the stress-strain curve. Where the initial non-linear stiffness associated with the crack closure effect is an important model feature, gapped contacts may be used. Gapped contacts are unbonded contacts with g 0 > 0 and are mechanically equivalent to open microcracks. Patel and Martin (2020) showed for a BPM representing a low-porosity Lac du Bonnet granite that the use of slit contacts enabled reproduction of the tensile-compressive elastic bimodularity effect which they reasoned was an indicator of the initial microcrack condition due to the differential stiffness response of microcracks under tensile and compressive loading. By increasing F s from 0 to 0.36, they found that |σ t | and σ c reduced by 40.8% and 24.9%, respectively, demonstrating that microcracks in a BPM have a greater effect on tensile strength than compressive strength. They also showed that the normalised crack initiation stress, σ ci /σ c , increased from 0.29 at F s = 0 to a more realistic value of 0.41 at F s = 0.36. It is noted, however, that a realistic σ ci /σ c of 0.46 was observed for the Rewan Sandstone case study with F s and F g both set to 0. The reproduction of a realistic σ ci /σ c is believed to be related to the good match between the porosity of the BPM and the natural porosity of the Rewan Sandstone. It is also noted that the BPM modelling of Patel and Martin (2020) was performed in 3D for which the use of spherical particles results in an even higher porosity than its 2D counterpart.
Irrespective of the rock type, it is recommended to set F s and F g to 0 for the purpose of automating the microproperty calibration as the use of randomly positioned slit and gapped contacts introduces a significant degree of variability that is detrimental to the ability of the calibration process to converge on a solution. Moreover, unless micro-Computed Tomography (μCT) is performed on the rock specimen (Ghamgosar et al. 2015;Knackstedt et al. 2006), the nature and spatial distribution of initial microcracks are generally unknown, and the specification of microcracks in a BPM is therefore not possible to validate directly. In the Rewan Sandstone case study, the target macroproperties were downgraded to 80% of the laboratory-measured values after Martin et al. (2011) to implicitly account for the influence of microcracks on the macroscopic strength and stiffness. Where explicit inclusion of microcracks in the BPM is required, it is recommended to follow a procedure similar to that of Patel and Martin (2020) whereby: (1) the microproperties are first calibrated to the target intact properties without microcracks; (2) an increasing proportion of slit contacts are introduced to downgrade the macroscopic strength and stiffness; and (3) gapped contacts are introduced to simulate the crack closure effect by increasing g 0 of the slit contacts until σ cc is matched.

Calibration of Post-peak Behaviour
In the automated calibration sequence presented in Fig. 8, ϕ g is set equal to ϕ j and is held constant throughout the calibration procedure. The use of a constant grain friction angle (or friction coefficient) is common practice in the BPM literature and was found to produce reasonable post-peak behaviour for the Rewan Sandstone case study. If, however, it were necessary to more closely match specific post-peak behaviour or a particular critical strain value, ϕ g could be selected to match the gradients of the post-peak portions of the stress-strain curves as suggested by Wu and Xu (2016). Furthermore, if a variable ϕ g were to be implemented, the degradation of surface asperities and frictional strength of Fig. 15 Performance of the bisection search calibration algorithm for the Rewan Sandstone case study microcracks as a function of local shear strain during the post-peak phase may be better simulated.

Algorithm Performance
The BSM was utilised in this study as it was the simplest and most reliable numerical root-finding algorithm that could be readily implemented in Itasca's FISH programming language. However, it is relatively slow compared to other rootfinding methods such as the Secant or Brent methods which may be accessed via the Python interface that is available in PFC and later versions. These methods were not utilised in this study but could potentially improve the efficiency of future versions of the calibration algorithm. The performance of the BSM-based algorithm with respect to calibration precision is presented in Fig. 15 for acceptance tolerances of 0.1%, 1%, 3%, and 5% and shows an exponentially reducing efficiency for more precise tolerance. The calibrations were performed on a desktop computer with an Intel i7-4770 K 3.50 GHz processor and 32 GB of memory with a typical runtime of 8 min per iteration for the Rewan Sandstone case study material loaded at the quasi-static strain rate with a model resolution of 20. Therefore, the calibration procedure can be expected to converge on a solution within approximately 13 h for an acceptance tolerance of 3% versus 28 h for an acceptance tolerance of 1%. No appreciable improvement in performance is achieved by increasing the acceptance tolerance beyond 3% while reducing the tolerance below 1% results in a large increase in convergence time. A tolerance of 1-3% is therefore recommended depending on the convergence time required by the user.
If a more precise tolerance is required, a global optimization method may be utilised as a secondary calibration step upon completion of the initial calibration. For maximum efficiency, the residual errors of the BSM procedure could be used to inform the initial bounds of the global optimization method and the set of microproperties producing the minimum global error would then be taken as the final calibrated set of properties. Unlike the sequential procedure proposed in this study, the optimization procedure would calibrate all microproperties in parallel and it would therefore be essential to minimise the initial parameter bounds.

Conclusions
The calibration of microproperties for the flat joint contact model is significantly more complex than its predecessors, the contact bond and parallel bond contact models, due to its unique microstructural features that allow for matching of the macroscopic compressive-tensile strength ratio, σ c /|σ t |; internal friction angle, ϕ; and Hoek-Brown constant, m i . In this study, a method for automating the calibration of FJM microproperties based on the bisection search method numerical root-finding algorithm and specific calibration sequencing is proposed. Application of the calibration scheme to a Rewan Sandstone case study demonstrated that both the target macroproperties and a realistic damage evolution were reproduced, including a normalised crack initiation stress of 0.46 and a normalised crack damage stress of 0.83 coinciding with a volumetric strain reversal.
Successful application of the automated calibration procedure is dependent on the selection of a valid set of microstructural properties which are held constant. These include the number of radial elements, N r ; the bond installation gap ratio, g i /D avg ; and the slit and gapped contact proportions, F s and F g . Sensitivity analyses of the automatically calibrated microproperties for varying g i /D avg and N r revealed that g i /D avg primarily controls the dilatancy, while N r primarily controls the compressive-tensile strength ratio, σ c /|σ t |. For the purpose of automating the calibration process, it is recommended to set F s and F g to 0 and instead capture the influence of microcracks via an empirical scale factor. Where the inclusion of explicit microcracks is required, it is recommended to follow a secondary calibration procedure similar to Patel and Martin (2020) whereby the proportion of slit and gapped contacts are increased from 0 until the required macroscopic strength and stiffness reductions and crack closure effects are reproduced.
Finally, it is noted that the particle assembly for the Rewan Sandstone case study was generated using the grainscaling method of Potyondy (2017) with a low target initial mean stress equal to the Earth's atmospheric pressure. This resulted in a porosity of 16% and provided a good match with the natural porosity of around 10-15% for the Rewan Sandstone. This is believed to be an important factor in the ability of the automated calibrated scheme to reproduce a realistic damage evolution. Future research is therefore recommended to focus on calibration of microproperties for rock materials with varying porosities and whether alternative DEMs such as the clumped, grain-based, and blockbased models may produce more realistic damage evolutions for these materials.