Textural and Mineralogical Controls on Rock Strength Elucidated Using a Discrete Element Method Numerical Laboratory

: Numerical modelling techniques such as the discrete element method are now well established and extensively used in many applications including solid earth geoscience, materials science, geotechnical engineering and rock mechanics. The potential for this technique in understanding comminution mechanisms has been identiﬁed as highly promising. This work utilizes the discrete element method as a numerical laboratory to conduct investigations relevant to comminution that would otherwise be costly or time-consuming to perform in the ﬁeld or laboratory. A benchmark numerical model for impact breakage of rock specimens is ﬁrst established and validated against results of controlled laboratory experiments. Thereafter, the model is utilized to systematically investigate the potential dependency of ore breakage properties upon the prevalence of pre-existing fractures, as well as the mineralogical composition of the ore. These numerical experiments serve to highlight the potential for quantitatively relating the mechanical response of ore to its textural and mineralogical characteristics. Tandem utilization of numerical and laboratory experimentation to formulate and test hypotheses is a promising avenue to illuminate such relationships.


Introduction
The Discrete Element method (DEM) is a well-established numerical technique for simulating and studying the flow of granular materials as well as the fracture of rocks and rock-like materials [1][2][3]. A notable advantage is the ability to simulate at a desired scale for different applications [4,5]. Its application to serve as a "numerical laboratory" for predicting rock response or fracture, whilst varying micro-structural parameters and loading conditions, has been demonstrated and comprehensively tested in rock engineering and fracture mechanics [6][7][8][9]. The concept of the DEM as a numerical laboratory entails prescribing the boundary conditions and specifying material properties in the model, and then systemically varying these to investigate the variability of macroscopic measurable quantities [10]. This approach yields quantitative hypotheses that are amenable to subsequent experimental verification [11]. Areas of application where this has been demonstrated include materials science [12], rock mechanics [11,[13][14][15], seismology [9] and mining sciences [16].
In minerals processing, the DEM has been demonstrated as a viable tool for simulating breakage of ore and equipment operation at known conditions [17][18][19][20]. The DEM has proven useful for gaining understanding of the underlying causes for operational challenges and to assist in the design of comminution equipment [21]. However, to date, the DEM has not been extensively applied as a numerical laboratory in comminution research, DEM has not been extensively applied as a numerical laboratory in comminution research, specifically to investigate the phenomenology of rock breakage and its dependence upon the textural characteristics of ore.
Research towards this end typically relies upon laboratory impact breakage experiments studies involving individual rock specimens [22]. The Short Impact Load Cell (SILC) has arisen as the preferred device for such studies [22]. The SILC is a hybrid of the drop weight tester and split Hopkinson bar [23] and is a shortened version of the ultrafast load cell [23]. Figure 1a shows a schematic of the SILC apparatus, while Figure 1b displays typical force-time profiles obtained from laboratory experiments. It consists of a pneumatic gripper that clamps and releases a steel ball at a prescribed height above a vertically oriented cylindrical steel rod on which the rock specimen rests. The steel rod is instrumented with strain gauges to measure the changing forces acting on rock specimens during breakage events. Analysis of the force time-histories permits estimation of the rock stiffness, fracture force or tensile strength, and fracture energy [24][25][26].
SILC investigations have shown that factors such as rock shape [29] and size [30], pre-weakening and mineralogical composition [22] cause variations in the measured mechanical response during impact breakage experiments [22,27,30]. Rock shape markedly contributes to the variability of mechanical properties measured in SILC experiments, potentially masking the influence of intrinsic textural characteristics of the ore. This has prompted research employing specimens of controlled shape and size [31,32], as well as synthetic rock specimens of homogeneous mineral composition [30,33]. Such studies have confirmed that the mechanical properties measured by the SILC apparatus depend systematically upon the textural characteristics of the rock specimens. Furthermore, these experiments show that the peak fracture force as measured by the SILC is proportional to the tensile strength of the specimens [30][31][32], offering the possibility to harmonize geotechnical and comminution mechanical parameters.
Having demonstrated that mechanical properties and textural characteristics are inter-related [31], the next logical step is to illuminate the quantitative nature of this relationship. A purely empirical approach would involve extensive laboratory investigations in which the textural characteristics of the rock specimens are carefully measured, or synthetic specimens are prepared with strictly controlled internal structure. Such an experimental campaign is likely to be laborious and expensive, with the ever-present concern that any resulting relationships may apply only to the specimens examined.
The primary objective of the presented research is to demonstrate that the bonded particle DEM may be employed as a numerical laboratory to investigate textural controls on comminution that are difficult to access via purely laboratory-based investigations.  [26]. (b) Typical force-time profiles obtained in SILC experiments [24,27,28].
SILC investigations have shown that factors such as rock shape [29] and size [30], preweakening and mineralogical composition [22] cause variations in the measured mechanical response during impact breakage experiments [22,27,30]. Rock shape markedly contributes to the variability of mechanical properties measured in SILC experiments, potentially masking the influence of intrinsic textural characteristics of the ore. This has prompted research employing specimens of controlled shape and size [31,32], as well as synthetic rock specimens of homogeneous mineral composition [30,33]. Such studies have confirmed that the mechanical properties measured by the SILC apparatus depend systematically upon the textural characteristics of the rock specimens. Furthermore, these experiments show that the peak fracture force as measured by the SILC is proportional to the tensile strength of the specimens [30][31][32], offering the possibility to harmonize geotechnical and comminution mechanical parameters.
Having demonstrated that mechanical properties and textural characteristics are interrelated [31], the next logical step is to illuminate the quantitative nature of this relationship. A purely empirical approach would involve extensive laboratory investigations in which the textural characteristics of the rock specimens are carefully measured, or synthetic specimens are prepared with strictly controlled internal structure. Such an experimental campaign is likely to be laborious and expensive, with the ever-present concern that any resulting relationships may apply only to the specimens examined.
The primary objective of the presented research is to demonstrate that the bonded particle DEM may be employed as a numerical laboratory to investigate textural controls on comminution that are difficult to access via purely laboratory-based investigations. Numerical simulations are conducted to identify potential relationships between mechanical properties and textural characteristics. Such numerically obtained relationships are treated as hypotheses, to be tested via subsequent controlled laboratory experiments. In the fol-lowing, two such numerical hypotheses are presented, governing the relationship between peak fracture force and the density of microcracks or mineralogical composition within rock specimens. For a binary phase mineralogical texture, it is shown that fracture force decreases proportionately with an increase in the concentration of the weaker component. Fracture force similarly decreases linearly with an increase in the surface area of initial fractures within pre-weakened rock specimens.
The following section describes the numerical methodology. This entails the description of the numerical SILC apparatus, model validation as well as the construction of numerical rock specimens with prescribed fracture networks or mineralogical composition. The subsequent section presents results of numerical experiments involving these damaged or heterogeneous numerical specimens. The paper concludes with a discussion of the findings, guidelines for future research and experimental validation and the primary conclusions.

Materials and Methods
A bonded particle DEM model [34] is employed for the numerical experiments presented herein. Bonded particle DEM models represent rock specimens as an assembly of DEM spheres. Adjacent spheres are bound together via cylindrical linear-elastic beams. Deformation of this bonded assembly induces forces and moments between the bound DEM spheres, which are used to incrementally update the positions and orientations of the spheres. When the shear and normal stresses within an individual beam exceed a prescribed failure criterion, the beam interaction is replaced with a frictional elastic-repulsive interaction. This method mimics the spontaneous formation of frictional fracture surfaces during crack propagation within brittle materials.
Past research has demonstrated that such a formulation reproduces key phenomenology in fracture mechanics such as Mode I tensile crack propagation [35] and wing-crack formation [2,36]. This research employs ESyS-Particle [37], an open-source software package that implements a bonded particle DEM. A detailed mathematical description of the bonded particle formulation of ESyS-Particle can be found elsewhere [38,39]. The software has been successfully utilized to study the behaviour of rock and rock-like materials in the past [15,[40][41][42][43][44].

Numerical Reproduction of the SILC Apparatus
A numerical reproduction of the SILC apparatus has been constructed, as illustrated in Figure 2. The steel ball is represented as a single, large DEM sphere, initially located a small distance above the DEM rock specimen and prescribed with a constant velocity commensurate with having fallen under gravitational acceleration from a user-specified height above the specimen.
The specimen itself is an assembly of non-overlapping DEM spheres filling a cylindrical volume, constructed using a volume-filling algorithm [45], as implemented in the GenGeo software that accompanies ESyS-Particle. The method involves firstly inserting non-overlapping spheres at random locations within the cylindrical volume, with radii in a specified range. Thereafter, four existing objects (spheres and/or volume boundaries) are selected and a new sphere is inserted to touch each of these, so long as the radius of the candidate sphere resides within the specified range. This sphere insertion procedure is repeated until the algorithm fails repeatedly to find a hole of sufficient size to insert new spheres. The resultant assembly consists of non-overlapping spheres that fill the prescribed volume with a final porosity that decreases as the range of sphere radii increases. For the simulations considered herein, the smallest spheres are one quarter the size of the largest spheres, resulting in a porosity of approximately 35%. The porosity of the DEM sphere packings is a major factor influencing the relationship between microscopic parameters governing interactions and the macroscopic mechanical response of bonded assemblies of DEM spheres. Consequently, the porosity of all numerical specimens is held fixed in this investigation.
sphere packings is a major factor influencing the relationship between microscopic parameters governing interactions and the macroscopic mechanical response of bonded assemblies of DEM spheres. Consequently, the porosity of all numerical specimens is held fixed in this investigation. The lower-most spheres comprising the specimen undergo elastic interactions with a planar substrate that remains stationary throughout simulations. Simple linear, elastic repulsive forces act between such spheres and the substrate. The net force acting on this substrate is recorded throughout simulations, producing force time-histories akin to that measured during laboratory SILC experiments.
At commencement of each simulation, all pairs of adjacent DEM spheres are identified and bound together via linear-elastic beam interactions. Four model parameters govern the mechanical properties of these beams and, by extension, the macroscopic mechanical properties of the specimen. These four parameters are: the elastic beam modulus (Yb), beam Poisson's ratio (νb), cohesive strength (Cb) and tangent of the internal angle of friction of the beams (tan θ).
In fracture mechanics investigations using the bonded particle DEM, it is typical to conduct simulations of uniaxial compression tests and modify these four parameters by trial-and-error, to match the measured geotechnical properties (Young's modulus, Yc, and unconfined compressive strength, σ ) of a particular rock type. In this study, an alternative, parametric approach has been employed to calibrate the mechanical properties of numerical rock specimens. Each of the four model parameters were systematically varied to construct empirical calibration relationships between these parameters and measured macroscopic geotechnical properties. The details of this calibration procedure are provided in Appendix A and the so-obtained calibration relationships are as follows: These calibration relationships are valid only for DEM sphere assemblies constructed using the packing algorithm described above, with a ratio of sphere radii satisfying Rmin:Rmax = 1:4 and the diameter of the largest DEM sphere less than one tenth the length The lower-most spheres comprising the specimen undergo elastic interactions with a planar substrate that remains stationary throughout simulations. Simple linear, elastic repulsive forces act between such spheres and the substrate. The net force acting on this substrate is recorded throughout simulations, producing force time-histories akin to that measured during laboratory SILC experiments.
At commencement of each simulation, all pairs of adjacent DEM spheres are identified and bound together via linear-elastic beam interactions. Four model parameters govern the mechanical properties of these beams and, by extension, the macroscopic mechanical properties of the specimen. These four parameters are: the elastic beam modulus (Y b ), beam Poisson's ratio (ν b ), cohesive strength (C b ) and tangent of the internal angle of friction of the beams (tan θ).
In fracture mechanics investigations using the bonded particle DEM, it is typical to conduct simulations of uniaxial compression tests and modify these four parameters by trial-and-error, to match the measured geotechnical properties (Young's modulus, Y c , and unconfined compressive strength, σ c ) of a particular rock type. In this study, an alternative, parametric approach has been employed to calibrate the mechanical properties of numerical rock specimens. Each of the four model parameters were systematically varied to construct empirical calibration relationships between these parameters and measured macroscopic geotechnical properties. The details of this calibration procedure are provided in Appendix A and the so-obtained calibration relationships are as follows: These calibration relationships are valid only for DEM sphere assemblies constructed using the packing algorithm described above, with a ratio of sphere radii satisfying Rmin:Rmax = 1:4 and the diameter of the largest DEM sphere less than one tenth the length of the smallest linear dimension of the numerical rock specimen. The length of the cylindrical specimen is twice its diameter, ensuring that the well-known aspect-ratio dependency of mechanical properties is avoided. Adherence to these constraints ensures that all numerical specimens constructed and calibrated according to these relationships have statistically identical macroscopic geotechnical properties, regardless of the precise locations and sizes of individual spheres comprising the specimens.

Model Validation
The purpose of a numerical laboratory is to quantitatively predict measurable responses of a modelled system, under changing initial and boundary conditions. To validate such a numerical laboratory, it is insufficient to simply demonstrate that the model reproduces the results of an individual laboratory test case; the typical method employed to validate DEM models in the literature. The numerical laboratory must demonstrably predict the results of laboratory experiments across a broad range of experimental conditions.
The herein proposed SILC numerical laboratory is validated with reference to the experimental results of Barbosa et al. [30], who conducted laboratory SILC experiments using 3D-printed synthetic rock specimens. These authors manufactured synthetic rocks using silica-sand and jet binder and processed in a three-dimensional (3D) printer programmed with the desired geometry (cylinder) and size. Specimens were produced with 20% binder saturation and dried at 80 • C for about 24 h [33]. The length to diameter ratio (L:D) of all specimens varied between 1:1 to 2:1. Figure 3 shows a picture of different sizes of the 3D printed cylinders that were manufactured for the laboratory studies. A subset of the manufactured cylinders, with diameters of 10.2 mm, were subjected to uniaxial compression tests, yielding measurements of the Young's modulus and uniaxial compressive strength of the synthetic specimens. Thereafter, the authors [30] conducted SILC experiments and reported the variation in peak fracture force, as a function of the diameter of the cylindrical specimens. These results comprise the target for the validation exercise. and sizes of individual spheres comprising the specimens.

Model Validation
The purpose of a numerical laboratory is to quantitatively predict measurable sponses of a modelled system, under changing initial and boundary conditions. To v date such a numerical laboratory, it is insufficient to simply demonstrate that the mo reproduces the results of an individual laboratory test case; the typical method employ to validate DEM models in the literature. The numerical laboratory must demonstra predict the results of laboratory experiments across a broad range of experimental con tions.
The herein proposed SILC numerical laboratory is validated with reference to experimental results of Barbosa et al. [30], who conducted laboratory SILC experime using 3D-printed synthetic rock specimens. These authors manufactured synthetic ro using silica-sand and jet binder and processed in a three-dimensional (3D) printer p grammed with the desired geometry (cylinder) and size. Specimens were produced w 20% binder saturation and dried at 80 °C for about 24 h [33]. The length to diameter ra (L:D) of all specimens varied between 1:1 to 2:1. Figure 3 shows a picture of different si of the 3D printed cylinders that were manufactured for the laboratory studies. A sub of the manufactured cylinders, with diameters of 10.2mm, were subjected to uniaxial co pression tests, yielding measurements of the Young's modulus and uniaxial compress strength of the synthetic specimens. Thereafter, the authors [30] conducted SILC expe ments and reported the variation in peak fracture force, as a function of the diameter the cylindrical specimens. These results comprise the target for the validation exercise Validation commences with calibration of the numerical model parameters, to ma the measured geotechnical properties of the 3D-printed. The above calibration relatio ships (Equations (1) and (2)) are employed for this purpose and the obtained model rameters are summarized in Table 1. The geotechnical properties measured in the expe ments, as well as the obtained values for the calibrated numerical specimens, are summ rized in Table 2. According to this table, a comparison of the Young's modulus and uniaxial compressive strength from both experiment and numerical simulation shows t there is no statistically significant difference based on the calculated standard deviati As such, the numerically calibrated parameters are considered a close match to the exp iment.
The standard deviation of mechanical properties for the numerical simulations is s nificantly lower than that of the experiments, demonstrating that the numerical mode an "idealized" representation of rocks that have heterogeneous structure yet homoge ous mechanical properties. This is a distinct advantage when employing the bonded p ticle DEM as a numerical laboratory. Specimens may be constructed with statistica Validation commences with calibration of the numerical model parameters, to match the measured geotechnical properties of the 3D-printed. The above calibration relationships (Equations (1) and (2)) are employed for this purpose and the obtained model parameters are summarized in Table 1. The geotechnical properties measured in the experiments, as well as the obtained values for the calibrated numerical specimens, are summarized in Table 2. According to this table, a comparison of the Young's modulus and the uniaxial compressive strength from both experiment and numerical simulation shows that there is no statistically significant difference based on the calculated standard deviation. As such, the numerically calibrated parameters are considered a close match to the experiment.
The standard deviation of mechanical properties for the numerical simulations is significantly lower than that of the experiments, demonstrating that the numerical model is an "idealized" representation of rocks that have heterogeneous structure yet homogeneous mechanical properties. This is a distinct advantage when employing the bonded particle DEM as a numerical laboratory. Specimens may be constructed with statistically identical mechanical properties, avoiding the inherent uncertainties arising from the greater variability of rock specimens employed in laboratory studies.
Having calibrated the numerical model parameters, a suite of numerical SILC simulations was conducted, to measure the variation in peak fracture force for cylindrical specimens of differing diameters, in accordance with the experimental conditions employed by Barbosa et al. [30]. For each specimen size, five different numerical samples were fractured using the numerical SILC apparatus. Each numerical sample contained approximately 40,000 DEM spheres. The measured peak fracture force for each of these was then averaged and the standard deviation calculated. The results are illustrated in Figure 4, in comparison with the experimental results of Barbosa et al. [30]. A statistical T-test was conducted to test whether the experimental and numerical results reflect samples from statistically different distributions. For the 2-tailed T-test, a p-value of 0.035 confirms that the two sets of results are statistically identical. The numerical laboratory matches not only the measured fracture force at the specimen size used for geotechnical characterization (10.2 mm) but also correctly predicts the observed change in peak fracture force as the specimen size is varied. Furthermore, visual inspection of the fragmentation obtained in the numerical simulations also closely matches that reported for the laboratory studies, as shown in Figure 5a,b. For both the laboratory experiments and the numerical simulations, the specimens are observed to cleave into two large fragments, with some small debris between the primary facture surfaces.
It is emphasized that no additional calibration of the numerical model was conducted in order to obtain agreement between the laboratory and numerical SILC results; calibration was only based on the geotechnical properties as reported by Barbosa et al. [30] for specimens with a diameter of 10.2 mm. This is considered a robust validation of the capacity of the numerical SILC model to quantitatively predict the mechanical response of brittle-elastic materials under a broad range of geometrical, initial and boundary conditions, a requisite feature for the model to be considered a numerical laboratory.

Construction of Pre-Weakened Numerical Rock Specimens
Pre-weakened numerical rock specimens are constructed using a geometrical met wherein pre-existing cracks are represented as a set of triangulated planes that are domly distributed within the specimen volume. For this study, cracks are constraine be square surfaces whose side length ranges between 0.8 mm and 2.0 mm. The orie tions of these crack surfaces are randomized in all three coordinate directions, to repre pre-weakening due a variety of sequential, independent comminution and materials h dling stages, as is typical of ore samples within a minerals processing circuit.
These crack surfaces are imposed as additional internal boundaries during const tion of the DEM specimens. The sphere packing algorithm of GenGeo ensures that serted spheres do not cross these surfaces, forming geometrically flat contacts betw spheres either side of the crack surfaces. These DEM spheres are initially unbonded, dergoing only linear elastic repulsion and frictional interactions should they come

Construction of Pre-Weakened Numerical Rock Specimens
Pre-weakened numerical rock specimens are constructed using a geometrical method wherein pre-existing cracks are represented as a set of triangulated planes that are randomly distributed within the specimen volume. For this study, cracks are constrained to be square surfaces whose side length ranges between 0.8 mm and 2.0 mm. The orientations of these crack surfaces are randomized in all three coordinate directions, to represent pre-weakening due a variety of sequential, independent comminution and materials handling stages, as is typical of ore samples within a minerals processing circuit.
These crack surfaces are imposed as additional internal boundaries during construction of the DEM specimens. The sphere packing algorithm of GenGeo ensures that inserted spheres do not cross these surfaces, forming geometrically flat contacts between spheres either side of the crack surfaces. These DEM spheres are initially unbonded, undergoing only linear elastic repulsion and frictional interactions should they come into

Construction of Pre-Weakened Numerical Rock Specimens
Pre-weakened numerical rock specimens are constructed using a geometrical method wherein pre-existing cracks are represented as a set of triangulated planes that are randomly distributed within the specimen volume. For this study, cracks are constrained to be square surfaces whose side length ranges between 0.8 mm and 2.0 mm. The orientations of these crack surfaces are randomized in all three coordinate directions, to represent pre-weakening due a variety of sequential, independent comminution and materials handling stages, as is typical of ore samples within a minerals processing circuit.
These crack surfaces are imposed as additional internal boundaries during construction of the DEM specimens. The sphere packing algorithm of GenGeo ensures that inserted spheres do not cross these surfaces, forming geometrically flat contacts between spheres either side of the crack surfaces. These DEM spheres are initially unbonded, undergoing only linear elastic repulsion and frictional interactions should they come into contact. This approach captures the essential features of pre-existing cracks: frictional surfaces that constrain deformation in their vicinity to be parallel and perpendicular to the fracture surface.
Illustrations of homogeneous and pre-weakened DEM rock specimens are shown in Figure 6. Figure 6a,b show the homogeneous specimen without micro-cracks and its bond networks, respectively, while Figure 6c shows a thin section from the bond network around the middle of the specimen in the Z-direction. Figure 6d and e show the specimen with micro-cracks and its bond networks, respectively, while Figure 6f shows a thin section from the bond network. The presence of some of the cracks denoted by holes within the bond network are highlighted with red ellipses. contact. This approach captures the essential features of pre-existing cracks: frictional surfaces that constrain deformation in their vicinity to be parallel and perpendicular to the fracture surface.
Illustrations of homogeneous and pre-weakened DEM rock specimens are shown in Figure 6. Figure 6a,b show the homogeneous specimen without micro-cracks and its bond networks, respectively, while Figure 6c shows a thin section from the bond network around the middle of the specimen in the Z-direction. Figure 6d and e show the specimen with micro-cracks and its bond networks, respectively, while Figure 6f shows a thin section from the bond network. The presence of some of the cracks denoted by holes within the bond network are highlighted with red ellipses. For the numerical SILC experiments, pre-weakened specimens with differing levels of damage were constructed. This is achieved by prescribing the total number of pre-existing cracks, Nc = 0, 30, 90, 180 and 300. For each of these damage levels, 30 DEM rock specimens with different random crack distributions were constructed.

Construction of Numerical Rock Specimens with Mineralogical Composition
The bonded particle DEM approach is further adapted to investigate the effect of mineralogical composition of the numerical rock specimens. In this study, binary specimens comprised of two distinct mineral phases are considered; one phase representing a valuable mineral and the other gangue. This simplifies the construction and reduces the complexity of analyzing results obtained for more complex, multi-mineral textures.
The numerical rock specimens are constructed such that the mechanical properties of one of the two mineral phases may be weakened with respect to the other mineral phase. For the purposes of discussion, it is assumed that the softer component is the "valuable mineral" whilst the other phase is the "gangue". The steps involved to build the numerical rock specimen with the desired mineralogy are as follows: • First generate a homogeneous rock specimen. For the numerical SILC experiments, pre-weakened specimens with differing levels of damage were constructed. This is achieved by prescribing the total number of pre-existing cracks, Nc = 0, 30, 90, 180 and 300. For each of these damage levels, 30 DEM rock specimens with different random crack distributions were constructed.

Construction of Numerical Rock Specimens with Mineralogical Composition
The bonded particle DEM approach is further adapted to investigate the effect of mineralogical composition of the numerical rock specimens. In this study, binary specimens comprised of two distinct mineral phases are considered; one phase representing a valuable mineral and the other gangue. This simplifies the construction and reduces the complexity of analyzing results obtained for more complex, multi-mineral textures.
The numerical rock specimens are constructed such that the mechanical properties of one of the two mineral phases may be weakened with respect to the other mineral phase. For the purposes of discussion, it is assumed that the softer component is the "valuable mineral" whilst the other phase is the "gangue". The steps involved to build the numerical rock specimen with the desired mineralogy are as follows:

•
First generate a homogeneous rock specimen. • Specify random "seed points" within the homogeneous specimen which represent the centroids of individual mineral grains.
• Group together spheres of the homogeneous rock specimen closest to each of these seed points; a Voronoi tessellation procedure. • Randomly assign each group of spheres to be either "valuable" or "gangue" in accordance with a prescribed probability, denoting the mineral grade of the specimen.
These four steps are illustrated in Figure 7. Figure 7a shows the initial homogeneous cylindrical specimen, while Figure 7b shows the insertion of seed points (coloured spheres). Grouping of DEM spheres into individual mineral grains, based upon their proximity to the nearest seed point, is shown in Figure 7c. Figure 7d shows the final two-phase specimen in which a percentage of the mineral grains are selected as representing the valuable mineral phase (red spheres) and gangue (light grey spheres). A 20% grade is illustrated in this case. • Specify random "seed points" within the homogeneous specimen which represent the centroids of individual mineral grains.

•
Group together spheres of the homogeneous rock specimen closest to each of these seed points; a Voronoi tessellation procedure.

•
Randomly assign each group of spheres to be either "valuable" or "gangue" in accordance with a prescribed probability, denoting the mineral grade of the specimen.
These four steps are illustrated in Figure 7. Figure 7a shows the initial homogeneous cylindrical specimen, while Figure 7b shows the insertion of seed points (coloured spheres). Grouping of DEM spheres into individual mineral grains, based upon their proximity to the nearest seed point, is shown in Figure 7c. Figure 7d shows the final twophase specimen in which a percentage of the mineral grains are selected as representing the valuable mineral phase (red spheres) and gangue (light grey spheres). A 20% grade is illustrated in this case. Figure 7. The four stages of creating specimens with a binary mineralogical texture. (a) homogeneous cylindrical specimen (b) homogenous specimen with seed points (coloured spheres) (c) DEMspheres close to each seed points are grouped together (d) grouped spheres assigned to be either "valuable"(red spheres) or "gangue" (grey spheres). For this study, six different mineral grades are examined: 0% (homogeneous gangue), 10%, 25%, 50%, 75% and 100% (homogenous valuable). Figure 8a-f depicts the typical distribution of valuable minerals within the numerical rock specimens at different mineral grades. These were randomly distributed for each specimen.
With respect to mechanical properties of these DEM mineralogical rock specimens, two key assumptions are made:  For this study, six different mineral grades are examined: 0% (homogeneous gangue), 10%, 25%, 50%, 75% and 100% (homogenous valuable). Figure 8a-f depicts the typical distribution of valuable minerals within the numerical rock specimens at different mineral grades. These were randomly distributed for each specimen.
With respect to mechanical properties of these DEM mineralogical rock specimens, two key assumptions are made:

Results
This section describes the results of two numerical experiments, one which systematically explores the influence of pre-weakening in the DEM rock specimen and a second which investigates the mechanical response of the DEM rock specimen as the mineralogical composition changes.

Breakage Behaviour under the Influence of Pre-Existing Cracks
Breakage investigations were conducted on three random specimens at each pre-existing crack level to predict the fracture behaviours. Predictions of the force-time histories and fracture patterns are presented from Figures 9-14.
An examination of the force-time histories in these Figures indicate that homogeneous rock specimens with 0 cracks (Figure 9a) and specimens with the least amount preexisting cracks, i.e., 30 cracks (Figure 10a) exhibit the highest consistency of force-time histories. Histories are fairly consistent in the presence of 90 and 180 cracks as shown in Figures 11a and 12a, respectively, while in the presence of 300 cracks (Figure 13a), the force-time history shows more variability.
The fracture patterns for three numerical specimens at each damage level are shown in Figure 9b-d to Figure 13b-d. The key observation here is also the consistency of the fracture pattern obtained at 0 crack level with specimens showing a tendency of splitting into almost two equal fragments (Figure 9b-d). A similar split of rock specimens into two fragments was observed in the presence of 30 cracks into the rock specimen, where some inconsistencies in the symmetry of the cleavage begin to appear. This can be seen with the fracture line of specimen 1 (Figure 10b), moving slightly more diagonally to the left-hand side whereas the fracture lines of specimens 2 and 3 tend to pass through the center of the specimen (Figure 10c,d). This behaviour can be attributed to the presence of micro-cracks driving the failure of crack propagation in variable directions within the rock specimen.
As the crack density intensifies, the irregularity of the fracture patterns within rock specimens become more apparent. These are shown with the fracture patterns at 90 cracks (Figure 11b-d), 180 cracks (Figure 12b-d) and 300 cracks (Figure 13b-d). Among these, fracture patterns at 180 cracks and 300 cracks show the highest irregularities. The irregularity of fracturing is further emphasized by examining the relative sizes of fragments produced. Figure 14 illustrates fragment masses for differing levels of damage. Relative

Results
This section describes the results of two numerical experiments, one which systematically explores the influence of pre-weakening in the DEM rock specimen and a second which investigates the mechanical response of the DEM rock specimen as the mineralogical composition changes.

Breakage Behaviour under the Influence of Pre-Existing Cracks
Breakage investigations were conducted on three random specimens at each preexisting crack level to predict the fracture behaviours. Predictions of the force-time histories and fracture patterns are presented from  Figures 11a and 12a, respectively, while in the presence of 300 cracks (Figure 13a), the force-time history shows more variability.
The fracture patterns for three numerical specimens at each damage level are shown in Figure 9b-d to Figure 13b-d. The key observation here is also the consistency of the fracture pattern obtained at 0 crack level with specimens showing a tendency of splitting into almost two equal fragments (Figure 9b-d). A similar split of rock specimens into two fragments was observed in the presence of 30 cracks into the rock specimen, where some inconsistencies in the symmetry of the cleavage begin to appear. This can be seen with the fracture line of specimen 1 (Figure 10b), moving slightly more diagonally to the left-hand side whereas the fracture lines of specimens 2 and 3 tend to pass through the center of the specimen (Figure 10c,d). This behaviour can be attributed to the presence of micro-cracks driving the failure of crack propagation in variable directions within the rock specimen.
As the crack density intensifies, the irregularity of the fracture patterns within rock specimens become more apparent. These are shown with the fracture patterns at 90 cracks (Figure 11b-d), 180 cracks (Figure 12b-d) and 300 cracks (Figure 13b-d). Among these, fracture patterns at 180 cracks and 300 cracks show the highest irregularities. The irregularity of fracturing is further emphasized by examining the relative sizes of fragments produced. Figure 14 illustrates fragment masses for differing levels of damage. Relative fragment masses are calculated by identifying all DEM spheres forming a connected cluster and accumulating the masses of these spheres. It is the fragment mass relative to the original specimen. As the damage level increases, the number and range of masses of fragments increases markedly, as compared with pristine specimens (0% cracks) or mildly damaged specimens. a b c d fragment masses are calculated by identifying all DEM spheres forming a connected cluster and accumulating the masses of these spheres. It is the fragment mass relative to the original specimen. As the damage level increases, the number and range of masses of fragments increases markedly, as compared with pristine specimens (0% cracks) or mildly damaged specimens.    fragment masses are calculated by identifying all DEM spheres forming a connected cluster and accumulating the masses of these spheres. It is the fragment mass relative to the original specimen. As the damage level increases, the number and range of masses of fragments increases markedly, as compared with pristine specimens (0% cracks) or mildly damaged specimens.    original specimen. As the damage level increases, the number and range of masses of frag ments increases markedly, as compared with pristine specimens (0% cracks) or mildly damaged specimens.          These results may also be examined with reference to the initial fracture surface area of pre-existing cracks, as illustrated in Figure 15. The mathematical expression for the calculation of the initial fracture surface area is given by Equation (3).
where, = initial fracture surface area (mm 2 ) = minimum length of crack (mm) = maximum length of crack (mm) = number of cracks. It is evident that the peak fracture force for the DEM SILC simulations decreases linearly as the initial fracture surface area increases. This represents a numerically derived hypothesis that is amenable to experimental verification. These results may also be examined with reference to the initial fracture surface area of pre-existing cracks, as illustrated in Figure 15. The mathematical expression for the calculation of the initial fracture surface area is given by Equation (3).
where, FSA = initial fracture surface area (mm 2 ) L min = minimum length of crack (mm) L max = maximum length of crack (mm) Nc = number of cracks. It is evident that the peak fracture force for the DEM SILC simulations decreases linearly as the initial fracture surface area increases. This represents a numerically derived hypothesis that is amenable to experimental verification.

Mechanical Response under the Influence of Varying Mineralogical Composition
In the results that follow, the mineral grade is systematically varied and the resultant change in fracture force measured. Six different mineral grades are examined: 0% (homogeneous gangue), 10%, 25%, 50%, 75% and 100% (homogenous valuable). Thirty SILC simulations are conducted for randomly differing specimens for each level of mineral grade, to provide representative statistics on fracture force versus mineral grade. Figure 16 shows a typical force-time history for the different grades. While the peak force decreases as the grade increases, it is further observed that there tends to be increasing oscillation of the force time-history in the presence of the valuable mineral. These oscillations are noticeable at even the lowest grade considered (10%). These oscillations are likely a consequence of increased incidence of elastic wave reflections at internal boundaries between the two mineral phases. The disappearance of such oscillations for specimens comprising 100% grade supports this conjecture.

Mechanical Response under the Influence of Varying Mineralogical Composition
In the results that follow, the mineral grade is systematically varied and the resultant change in fracture force measured. Six different mineral grades are examined: 0% (homogeneous gangue), 10%, 25%, 50%, 75% and 100% (homogenous valuable). Thirty SILC simulations are conducted for randomly differing specimens for each level of mineral grade, to provide representative statistics on fracture force versus mineral grade. Figure 16 shows a typical force-time history for the different grades. While the peak force decreases as the grade increases, it is further observed that there tends to be increasing oscillation of the force time-history in the presence of the valuable mineral. These oscillations are noticeable at even the lowest grade considered (10%). These oscillations are likely a consequence of increased incidence of elastic wave reflections at internal boundaries between the two mineral phases. The disappearance of such oscillations for specimens comprising 100% grade supports this conjecture.

Mechanical Response under the Influence of Varying Mineralogical Composition
In the results that follow, the mineral grade is systematically varied and the resultan change in fracture force measured. Six different mineral grades are examined: 0% (homo geneous gangue), 10%, 25%, 50%, 75% and 100% (homogenous valuable). Thirty SILC simulations are conducted for randomly differing specimens for each level of minera grade, to provide representative statistics on fracture force versus mineral grade. Figure 16 shows a typical force-time history for the different grades. While the pea force decreases as the grade increases, it is further observed that there tends to be increas ing oscillation of the force time-history in the presence of the valuable mineral. These os cillations are noticeable at even the lowest grade considered (10%). These oscillations ar likely a consequence of increased incidence of elastic wave reflections at internal bound aries between the two mineral phases. The disappearance of such oscillations for speci mens comprising 100% grade supports this conjecture.    Homogeneous specimens with zero grade, i.e., without any inherent weak-phase mineral, have the highest fracture force of approximately 1.2 kN. The force decreases as the composition of the grade increases, reducing to 0.1 kN for specimens which are entirely composed of the weaker phase (i.e., 100% grade).  Figure 17 shows the variation of the predicted fracture force as the grade increases. Homogeneous specimens with zero grade, i.e., without any inherent weak-phase mineral, have the highest fracture force of approximately 1.2 kN. The force decreases as the composition of the grade increases, reducing to 0.1 kN for specimens which are entirely composed of the weaker phase (i.e., 100% grade). Preferential exposure (liberation) of mineral phases is also examined. Figure 18a-r summarize the locations and distributions of the valuable mineral, the fracture paths and locations within the rock specimen as well as the fragmentation obtained for different mineral compositions. At 0% grade (Figure 18a), the fracture path travels almost vertically downward (Figure 18b) with the specimen spitting into almost two halves (Figure 18c). For 10% grade (Figure 18d), fracture paths tend to flow around the valuable minerals (Figure 18e) which results in an irregular fragmentation (Figure 18f). At 25% and 50% grade composition (Figure 18g,j), the fracture paths and directions are controlled and directed along the densely populated regions of valuable minerals which give rise to tributaries within the fracture patterns (Figure 18h,k) and associated irregular fragmentation of the specimens (Figure 18i,l). However, at 75% and 100% grade (Figure 18m,p), when the rock specimen is predominantly composed of valuable minerals, there is no evidence for preferential fracture propagation. Instead, the fracture locations are randomly distributed (Figure 18n,q). This anomalous behaviour is interpreted as a deviation of the mechanical properties of the specimens from distinctly brittle-elastic to a more ductile nature. A lack of distinct through-going fractures, as well as localized indentation and minor fragmentation at the contact points (Figure 18o,r) are further evidence of the ductility of specimens comprised of predominantly weaker minerals. Preferential exposure (liberation) of mineral phases is also examined. Figure 18a-r summarize the locations and distributions of the valuable mineral, the fracture paths and locations within the rock specimen as well as the fragmentation obtained for different mineral compositions. At 0% grade (Figure 18a), the fracture path travels almost vertically downward (Figure 18b) with the specimen spitting into almost two halves (Figure 18c). For 10% grade (Figure 18d), fracture paths tend to flow around the valuable minerals ( Figure 18e) which results in an irregular fragmentation (Figure 18f). At 25% and 50% grade composition (Figure 18g,j), the fracture paths and directions are controlled and directed along the densely populated regions of valuable minerals which give rise to tributaries within the fracture patterns (Figure 18h,k) and associated irregular fragmentation of the specimens (Figure 18i,l). However, at 75% and 100% grade (Figure 18m,p), when the rock specimen is predominantly composed of valuable minerals, there is no evidence for preferential fracture propagation. Instead, the fracture locations are randomly distributed (Figure 18n,q). This anomalous behaviour is interpreted as a deviation of the mechanical properties of the specimens from distinctly brittle-elastic to a more ductile nature. A lack of distinct through-going fractures, as well as localized indentation and minor fragmentation at the contact points (Figure 18o,r) are further evidence of the ductility of specimens comprised of predominantly weaker minerals. Figure 18. Locations and distributions of valuables within rock specimens, fracture locations and paths within rock specimens, fragmentation of specimens. (a) rock specimen with 0% mineral grade, (b) Fracture locations and patterns at 0% mineral grade, (c) Fragment mass at 0% mineral grade, (d) rock specimen with 10% mineral grade, (e) Fracture locations and patterns at 10% mineral grade, (f) Fragment mass at 10% mineral grade, (g) rock specimen with 25% mineral grade, (h) Fracture locations and patterns at 25% mineral grade, (i) Fragment mass at 25% mineral grade, (j) rock specimen with 50% mineral grade, (k) Fracture locations and patterns at 50% mineral grade, (l) Fragment mass at 50% mineral grade, (m) rock specimen with 75% mineral grade, (n) Fracture locations and patterns at 75% mineral grade, (o) Fragment mass at 75% mineral grade, (p) rock specimen with 100% mineral grade, (q) Fracture locations and patterns at 100% mineral grade and (r) Fragment mass at 100% mineral grade

Discussion
The primary objective of the research described herein is to demonstrate that the bonded particle DEM may be employed as a numerical laboratory to investigate textural controls on comminution. The approach involves first validating that the DEM model for a given experimental apparatus, the SILC in this instance, can quantitatively predict the results of laboratory control experiments. The validation exercise conducted herein demonstrates that the proposed DEM SILC model predicts the variation in peak fracture force as both the impact energy (drop height) and physical dimensions of specimens are varied. Qualitative analysis further demonstrates that the DEM SILC model predicts fracture patterns and fragmentation, in reasonable agreement with previously reported laboratory observations.
Having established the validity of the DEM SILC model, it is employed to conduct numerical control experiments that are difficult or tedious to conduct in the laboratory. This work focusses upon two such numerical control experiments, examining the variation in peak fracture force of cylindrical specimens with prescribed pre-existing cracks or differing binary mineral composition. Whilst simplified, these numerical control experiments yield quantitative predictions on the relationship between peak fracture force and textural characteristics of rock.

Breakage Behaviour under the Influence of Pre-Existing Cracks
The effect of pre-damage was investigated by specifying five levels of pre-damage within the numerical simulations. These pre-weakened rock specimens were simplified representations of rock specimens with inherent flaws or pre-weakening, with the added Figure 18. Locations and distributions of valuables within rock specimens, fracture locations and paths within rock specimens, fragmentation of specimens. (a) rock specimen with 0% mineral grade, (b) Fracture locations and patterns at 0% mineral grade, (c) Fragment mass at 0% mineral grade, (d) rock specimen with 10% mineral grade, (e) Fracture locations and patterns at 10% mineral grade, (f) Fragment mass at 10% mineral grade, (g) rock specimen with 25% mineral grade, (h) Fracture locations and patterns at 25% mineral grade, (i) Fragment mass at 25% mineral grade, (j) rock specimen with 50% mineral grade, (k) Fracture locations and patterns at 50% mineral grade, (l) Fragment mass at 50% mineral grade, (m) rock specimen with 75% mineral grade, (n) Fracture locations and patterns at 75% mineral grade, (o) Fragment mass at 75% mineral grade, (p) rock specimen with 100% mineral grade, (q) Fracture locations and patterns at 100% mineral grade and (r) Fragment mass at 100% mineral grade.

Discussion
The primary objective of the research described herein is to demonstrate that the bonded particle DEM may be employed as a numerical laboratory to investigate textural controls on comminution. The approach involves first validating that the DEM model for a given experimental apparatus, the SILC in this instance, can quantitatively predict the results of laboratory control experiments. The validation exercise conducted herein demonstrates that the proposed DEM SILC model predicts the variation in peak fracture force as both the impact energy (drop height) and physical dimensions of specimens are varied. Qualitative analysis further demonstrates that the DEM SILC model predicts fracture patterns and fragmentation, in reasonable agreement with previously reported laboratory observations.
Having established the validity of the DEM SILC model, it is employed to conduct numerical control experiments that are difficult or tedious to conduct in the laboratory. This work focusses upon two such numerical control experiments, examining the variation in peak fracture force of cylindrical specimens with prescribed pre-existing cracks or differing binary mineral composition. Whilst simplified, these numerical control experiments yield quantitative predictions on the relationship between peak fracture force and textural characteristics of rock.

Breakage Behaviour under the Influence of Pre-Existing Cracks
The effect of pre-damage was investigated by specifying five levels of pre-damage within the numerical simulations. These pre-weakened rock specimens were simplified representations of rock specimens with inherent flaws or pre-weakening, with the added benefit of being consistent replications for control and statistical reliability through nu-merical experiments. Results indicated that the fracture patterns became increasingly irregular as the level of pre-damage increased, as more opportunities for the propagation of cracks occurred within specimens. From the perspective of fracture mechanics, this result is not unexpected. Stress concentrates at the tips of pre-existing cracks, forming preferential pathways for crack propagation that would be otherwise absent from a pristine rock specimen. These stress concentrations promote reactivation of pre-existing cracks rather than absorbing surface energy to form new fracture surfaces. Consequently, the required fracture force decreased as the initial fracture surface area of the specimens increased by virtue of the pre-damage. Results of the simulations highlighted this to be a linear relationship (Figure 15), which was not unexpected given that the elasticity of all specimens remained unchanged [34].
The reported linear relationship between peak fracture force and initial fracture surface area is a numerically derived hypothesis, amenable to laboratory verification. Such verification might be achieved through a variety of means. One approach might be to prepare 3D-printed specimens with prescribed internal fractures matching the fracture networks employed in the numerical experiments. Breakage of these synthetic, damaged specimens in the SILC will provide measurements of peak fracture force that may be directly compared with the numerical predictions. An alternative procedure for verifying this numerical hypothesis might involve controlled pre-weakening of rock specimens via compression or low-energy impact, non-destructive imaging of damaged specimens to estimate internal fracture surface area, followed by breakage in a SILC device. Such an experimental campaign should be sufficient to verify whether the reported linear relationship for peak fracture force versus fracture surface area is a reasonable hypothesis.

Breakage Behaviour under the Influence of Varying Mineralogical Composition
The second examination of the potential application of the proposed methodology was the consideration of mineralogy. For simplification, the construction of specimens specified two mineral phases, a representation of typical ore bodies known to exhibit such distinct behaviour between the mineral and gangue. The grade was systematically varied using mechanical properties to create specimens over five varieties spanning homogenous gangue to homogenous mineral.
The variability in force time histories and nature of fracture were found to be most consistent when closer to homogeneity of either form. Greater variability was observed in the fracture force as the grade increased from completely gangue to become binary, decreasing subsequently as the grade increased to 100% mineral. The fracture force decreased with an increase in ore grade, with fracture occurring predominately within regions in which the weaker mineral phase was present. This continued to occur until the numerical experiments with 100% grade, at which the change in mechanical properties from brittleelastic to more ductile behaviour meant that specimens yielded readily under low force with no discernible fracture pattern.
The fracture of the rock specimen along the path of the mineral suggests that crack propagation is enhanced under the presence of a weaker phase. This leads to a more exposed surface area (liberation) of the weak phase. These findings demonstrate notable application for investigations of mineral liberation, in which the goal of exposing the highest surface area of the valuable material is the desired outcome. Prior work has examined the manner of preferential breakage through laboratory testing on common variabilities of ore [46]. It is worth highlighting that the main difficulty of these studies lies in the impracticability of being able to facilitate closely controlled experiments. This is mainly due to the inherent variability of specimens in internal characteristics such as texture and shape, as well as controlling the external environment they have been subject to and establishing a repeatable mode of testing. Such issues have been evident in the inconsistencies found in the breakage of ores even considered to be homogenous [47]. The methodology exhibited in this study can be utilized as a means of examining a wide range of effects associated with simple to complex mineralogical structures, with an added advantage of being able to consistently replicate both the specimen characteristics and the manner of impact. This would subsequently enable for the identification of specific scenarios to be closely studied in greater detail via experiments in which consistent specimens can be generated, e.g., by 3-d printing [33] or selective preparation of rock specimens whose mineralogical composition is measured via non-destructive imaging techniques.

Conclusions
In this work, DEM was utilized as a numerical laboratory to examine the effect of specimen pre-weakening and varying mineralogical composition on fracture behaviour. The methodology to calibrate BPM parameters to meet macroscopic mechanical properties was demonstrated to be capable of generating results consistent with standard geotechnical tests. It was further shown that empirical calibration relationships could be used to tune these across a range of properties characteristic of common rock types.
For the first application, the effect of pre-damage on the fracture force, nature and extent of fracture was investigated. An inverse linear relationship was exhibited between the fracture force and fracture surface area, highlighting the dominant role of pre-weakening in decreasing the amount of force required to cause fracture while increasing the resulting extent. This observation offers a hypothesis which warrants further scrutiny through further simulations coupled with closely controlled experiments which can investigate the optimal form of pre-damage which maximizes the efficiency of breakage in terms of energy expended to produce a given fracture surface area.
Secondly, this paper also assessed the potential of the numerical approach for studying the effect of ore texture on fracture behaviour. Ore texture is increasingly being considered as the next frontier in the design of flexible circuits of the future, the goal being to advantageously utilize information regarding mineralogy as a proxy to identify the most efficient manner of processing. Obtaining statistically representative information on the patterns of different textures is a crucial step toward this goal. Numerical experiments with a binary mineralogical composition indicated that fracture of specimens occurred mainly along the path of the weaker phase. This highlighted a promising opportunity to facilitate further studies on mineralogical characteristics that can be exploited to enhance liberation. While the texture considered was a simplified scenario, future simulations can consider more complex architectures and a range of parameters consistent with measurements from prior work. This can be supplemented by closely controlled laboratory tests to derive meaningful relationships that can be put to use toward quantitatively characterizing fracture using texture.
Author Contributions: T.O.: conceptualization, data curation, formal analysis, investigation, methodology, writing original draft; L.B.: Conceptualization; funding acquisition, project administration, resources, supervision, writing original draft, reviewing and editing manuscript; D.W.: conceptualization, data curation, methodology, supervision, reviewing and editing manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Acknowledgments:
Simulations were performed using High performance computers provided by the University of Cape Town's ICTS Computing team (http://hpc.uct.ac.za). Details of the computational resources that were allocated are as follows: cluster type Dell C6420, 40 cores, 294GB RAM. For the research presented herein, each simulation employed 5 processor cores, with typical execution times of between four and five hours. Multiple simulations were conducted simultaneously for parametric studies, thus reducing the real time required to conduct the numerical experiments.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Calibration of BPM Parameters
In the rock mechanics literature, it is typical to calibrate the macro-mechanical properties of BPM models utilising numerical reproductions of standard geotechnical tests such as the uniaxial compression test (UCT) and either direct or indirect tension tests [48,49]. These model parameters are then applied to the DEM environment specific to the experimental apparatus under investigation. Following this, numerical results are compared against similar laboratory tests to verify that the selected model parameters produce results reasonably consistent with experimentation without the need for further tuning and refinement [49,50].
A range of calibration methodologies have been identified in literature, including trial and error, design of experiments and machine learning or neural networks approaches [51][52][53]. For this study, an approach called the Systematic Parametric Study (SPS) was used to calibrate BPM parameters against standard measures of macro-mechanical properties. The SPS approach, which takes a form similar to design of experiments, entails systematically varying model parameters against macroscopic mechanical properties. These are then regressed against a fitting function to describe an empirical relationship between the model parameters and macroscopic mechanical properties.
Four BPM parameters were calibrated: elastic beam modulus (Y b ), beam Poisson's ratio (ν b ), cohesive strength (C b ) and tangent of the internal angle of friction (tan θ). The first two parameters govern the elastic response of the DEM specimen, whilst the latter two govern the onset and nature of the brittle failure. UCT simulations were conducted to establish the relationship between these four parameters and DEM specimen properties such as the macroscopic Young's modulus and unconfined compressive strength. The macroscopic elastic properties of a given specimen could subsequently be tuned to a desired value using the established calibration relationships.

Appendix A.1. UCS Simulations
A set of Uniaxial Compressive Strength (UCS) simulations were conducted on a numerical cylindrical-drill core specimen as shown in Figure A1, performed according to International Society of Rock Mechanics (ISRM) and American Society for Testing and Materials (ASTM) geotechnical standards [54][55][56][57]. A cylindrical geometry with an L:D ratio of 2:1 (20 mm:10 mm) was packed with bonded DEM-spheres with a radius ranging between 0.08 mm to 0.32 mm using the volume filling algorithm [45]. Planar walls were placed on at the top and bottom of the cylinder as shown in Figure 3 and prescribed with a downward speed from the top and an upward speed from the bottom of 0.125 m/s to impart compressive stress. The total instantaneous force on the DEM-spheres was recorded and the axial stress was calculated by dividing the instantaneous total force by the initial cross-sectional area of the specimen, in accordance with the standard laboratory procedure. The corresponding strain was determined by monitoring the mean vertical position at 25 and 75% of the height of the specimen as indicated in the Figure. This was calculated by dividing the relative vertical distance by the initial distance between these layers of DEM-spheres.  Figure A2 illustrates a typical stress-strain curve obtained from the numerical UCT, which is observed to follow conventional experimental results. The Young's modulus of Figure A1. A typical DEM specimen employed for uniaxial compression calibration simulations. Figure A2 illustrates a typical stress-strain curve obtained from the numerical UCT, which is observed to follow conventional experimental results. The Young's modulus of the rock was determined by calculating the slope of the linear portion of the stress-strain curve within zone 1. The peak point between zones 2 and 3 of the graph determines the unconfined compressive strength of the specimen. Figure A1. A typical DEM specimen employed for uniaxial compression calibration simulations. Figure A2 illustrates a typical stress-strain curve obtained from the numerical UCT, which is observed to follow conventional experimental results. The Young's modulus of the rock was determined by calculating the slope of the linear portion of the stress-strain curve within zone 1. The peak point between zones 2 and 3 of the graph determines the unconfined compressive strength of the specimen. Figure A2. Axial stress vs. axial strain and curves for a typical uniaxial compression simulation with an illustration of the calculated Young's modulus and unconfined compressive strength.  Figure A3 shows a graph of the macroscopic Young's modulus against elastic beam modulus at different Poisson's ratios. It is observed that there is a linear increase in Y c as the Y b increases and an incremental vertical spread of the data points, i.e., different slopes. This suggests that the Y c is influenced by both Y b and v b . Similar trends can also be observed in the graph of UCS versus the C b at different values of θ as depicted in Figure A4. The elastic beam parameters (Y b and v b ) were found to have no influence on the UCS. Likewise, the beam strength properties (C b and θ) did not influence the macroscopic Young's modulus of the DEM specimens. angle (tanθ) Figure A3 shows a graph of the macroscopic Young's modulus against elastic beam modulus at different Poisson's ratios. It is observed that there is a linear increase in Yc as the Yb increases and an incremental vertical spread of the data points, i.e., different slopes. This suggests that the Yc is influenced by both Yb and vb. Similar trends can also be observed in the graph of UCS versus the Cb at different values of θ as depicted in Figure A4. The elastic beam parameters (Yb and vb) were found to have no influence on the UCS. Likewise, the beam strength properties (Cb and θ) did not influence the macroscopic Young's modulus of the DEM specimens.  In order to isolate and study the effect of vb independently, Yc was divided by Yb and plotted against vb in Figure A3. This collapsed the data points into an almost linear (monotonic) trend. A similar technique was followed by plotting UCS/Cb against tan θ which resulted in a similar trend.
From inspecting the correlation between the BPM and macroscopic parameters, fitting functions in Equations (A1) and (A2) were proposed. In order to isolate and study the effect of v b independently, Y c was divided by Y b and plotted against v b in Figure A3. This collapsed the data points into an almost linear (monotonic) trend. A similar technique was followed by plotting UCS/C b against tan θ which resulted in a similar trend.
From inspecting the correlation between the BPM and macroscopic parameters, fitting functions in Equations (A1) and (A2) were proposed.
K i (i = 1, 2, 3 and 4) are the constants determined from the empirical fit. These reflect different elements of the spatial relationship of connecting adjacent DEM-spheres within the bonded rock specimens [35].
Non-linear regressions of the simulation data points were conducted using the SciPy optimize curve fit function [58] which returned the calculated parameters (K-values) and the coefficient of determination (R 2 ) as summarized in Table A2. Observance of the fitted model and high R-squared suggested the correlations were reasonable to adopt. These constants were substituted into Equations (1) and (2) to establish the relationships between the BPM and macroscopic mechanical properties and selected as a basis for BPM calibration of rock specimens. Figures A5 and A6 show the fitted empirical relationships against collapsed data points for the two scenarios, while Figures A7 and A8 summarize the fitted relationships across the range of tests on a surface dot-plot.