Grain-scale analysis of proppant crushing and embedment using calibrated discrete element models

Proppant crushing and embedment in hydraulically-induced fractures is a major drawback to the recovery of unconventional oil/gas and geothermal energy production. This study provides a grain-scale analysis of the fracture evolution mechanisms of proppant crushing, rock fracture damage during proppant embedment, the influence of realistic reservoir/fracture fluid on proppant embedment, and the behaviour of proppant packs subjected to in-situ stresses using a discrete element modelling (DEM) approach. The results of this study reveal that the selection of an appropriate proppant type based on the nature of the reservoir formation plays a vital part in quantifying the degree of proppant crushing and embedment within fractures. The utilisation of frac-sand proppants instead of ceramic proppants in shallow soft sedimentary-based siltstone formations reduces proppant embedment up to 88%. However, whatever the depth of the fracture, the injection of ceramic proppants into granite-based geothermal formations is preferred to that of frac-sand proppants due to their lower proppant embedment and greater crush resistance. DEM analysis detected rock-spalling during the proppant embedment process, which ultimately led to the initiation of tensile-dominant secondary fractures in rocks. Fracture initiation, propagation, and coalescence during proppant crushing are analysed using calibrated DEM proppant-rock assemblies. Importantly, this study reveals that the saturation of formation rocks with fracturing/reservoir fluids may cause a significant increase in proppant embedment. Furthermore, proppant crushing, embedment, and re-arrangement mechanisms in proppant packs with different proppant distributions are analysed in this comprehensive numerical study.


Introduction
Proppants used during hydraulic fracturing play a pivotal role in determining the ultimate fracture permeability within fractures [60]. Generally, proppants are used to keep fracture walls open by resisting high closure stresses. In 1947, silica granules were first utilised as a proppant agent, and to date, numerous types of proppants with concentrations up to 20lbm/gal have been used for fracturing treatment (Montgomery and Smith, [42]). With the advancement of fracking in unconventional oil/gas reservoirs, proppants also evolved from traditional silica granules to novel proppant types such as ceramic, resin-coated proppants, frac-sand and many more with various innovative mechanical improvements [7,22,32,35,50,57]. Due to their cost and availability for unconventional oil/gas recovery, ceramic, resin-coated, and frac-sand proppants have become popular types of proppants (Gallagher [19]). However, unforeseen oil/gas recovery reductions experienced from highly-impermeable unconventional deposits have regularly led to questions about proppant performance in fractures [15]. It is now understood that proppant degradation mechanisms may occur in fracture deposits under extreme reservoir conditions, and of these mechanisms, crushing and embedment of proppants require close attention [2,21,44]. The simultaneous processes of mechanical instability of proppants which occur upon resisting high fracture closure stress and proppant grain indentation into the fracture surface can lead to a drastic increase in fracture tortuosity and fracture aperture reductions, lowering the fracture permeability [14]. Although many macro-scale experimental and numerical studies have been conducted to understand the ultimate impact of crushing and embedment, the realistic proppant and rock response upon undergoing crushing and embedment is yet to be fully studied. Therefore, the numerical approach conducted using grain-scale discrete element modelling (DEM) provides comprehensive insights into proppantrock interactions in realistic fracture reservoir environments.
The phenomena of crushing and embedment of proppants in fractures have been investigated using experimental approaches with the assistance of high-stress uniaxial and triaxial test conditions. The studies by Vincent [51] revealed that the crushing of proppants might reduce the fracture conductivity by up to 99%. Laboratory-scale experimental analysis conducted by Bandara et al. [5] on the response of proppant subjected to high stresses revealed that the generation of fines as a result of crushing of ceramic proppant (6% fines), resin-coated sand (35% fines), and frac-sand (51% fines) can reduce proppant pack porosities by 75%, 88% and 98%, respectively. In addition, Alramahi and Sundberg [2] showed that the embedment of high-strength ceramic proppants into shale surfaces may reduce the conductivity by nearly six orders of magnitude at higher stresses. Moreover, experimental analyses by Kumar et al. [31], Akrad et al. [1], Corapcioglu [11], and Rybacki et al. [48] showed that parameters such as rock mineralogical content, type of formation rock, saturation conditions, saturated exposure time, and rock creeping effects may further contribute to proppant embedment. These studies provided an understanding of the mechanical, microstructural, mineralogical, morphological and chemical responses of individual proppants, proppant packs and formation rocks exposed to high-stress reservoir environments. Although these studies have provided insights into the macro-scale analysis of proppant crushing and embedment mechanisms, they were unable to provide non-destructive fracture analysis at a granular level [10]. Therefore, the need to study the proppant fracture mechanisms and embedment into the formation rock has motivated the present study.
Brittle, granular materials like proppants and rocks tend to crush, fracture, and fragment instead of plastically deforming when subjected to high stresses [47]. Therefore, understanding proppant-rock fracture initiation, fracture propagation, and fracture coalescence mechanisms during proppant crushing, and capturing proppant-rock deformation upon experiencing proppant embedment are of crucial importance [39,61]. To date, many finite element models and continuum models have been developed to study grain crushing and embedment mechanisms [26,36,40,58]. However, the inability to accurately observe the evolution processes of grain crushing and deformation during proppant crushing and embedment is a serious drawback [65]. Therefore, discrete element modelling (DEM) provides a better alternative for the quantitative study of this complex behaviour in brittle materials at a granular scale [17,18,37,59]. In the current study, the DEM approach was used to study the complex behaviour of proppant-rock interactions during crushing and embedment. In addition, the crushing and embedment behaviour of different types of proppants and the response of formation rocks are comprehensively analysed. Moreover, the effect of proppant packing and re-arrangement of particles caused by high closure stresses is analysed at a granular-scale.

Discrete element modelling (DEM)
DEM is a numerical modelling technique with the ability to portray the mechanical behaviour of granular assemblies in terms of particle-like discrete elements such as discs and spheres where the interaction of the particles occurs at particle-to-particle contacts and particle motion is governed by Newton's second Law [12]. In addition, an assembly can be simulated as a bonded particle model (BPM), where bonds between particles can be simulated in order to mimic the brittle behaviour of materials such as rock [45]. Generally, a contact between two particles possesses finite normal and shear stiffnesses, and the mechanical behaviour of the assembly is interpreted by the movement of each particle and the force and moment acting at each contact [46]. Hence, the utilisation of BPM combined with DEM to simulate proppant crushing and embedment behaviour has gained the attention of many researchers [13], Zheng and Tannant [66,67]. However, the accuracy of a DEM simulation relies on selecting an appropriate contact model to represent the mechanical behaviour of the assembly and assigning individual microparameters. In the present study, an existing linear-parallel contact model was selected in Particle Flow Code 3D (PFC3D) to accurately model the mechanical behaviour of granular proppant and rock assemblies. In the linear parallel-bonded contact model, two bond types are available: (1) contact bond (presented in terms of an elastic spring of contact normal k n and shear stiffness k s acting at the point of the contact which permits only force to be transferred) and (2) parallel bond (presented in terms of a set of elastic springs dispersed over a finite section of the contact plane and centred at the contact point which can resist the moment generated as a result of particle rotation) (Lisjak and Grasselli, [38]. Cho et al. [8] show that the contact normal stiffnesses (k n and k s ) remain active upon the breakage of the bond as long as particles stay in contact. However, parallel bond stiffness (k n -parallel bond normal stiffness and k s -parallel bond shear stiffness) is permanently lost when the parallel bond breaks (see Fig. 1). Despite the proppants not being perfect spheres, both proppant types were modelled assuming they were spherical with a diameter of 1 mm (see Fig. 2). In addition, both proppants were modelled using a cluster of homogenous particles. For both proppant types, a uniform particle size of 0.03 mm was selected by iteratively changing the size until no significant variation in the corresponding mechanical strengths was obtained between the laboratory and the calibrated model. Single-grain diametrical compression tests were simulated at an axial loading rate of 0.08 ms -1 to ensure a quasi-static loading environment. All the other required micro-parameters for the calibrated proppants were assigned by iteratively varying them using a trial-and-error method until the stress-strain curve from the numerical model followed a similar pattern to that of the laboratory test specimens. The micro-parameters of the calibrated ceramic and frac-sand proppants are listed in Table 1. A similar proppant calibration technique was adopted in our previous study; Bandara et al. [3]. Figure 3 shows the stress-strain response obtained from the experimentally-tested single grain crushing tests and Fig. 1 a Illustration of a combined parallel and contact bond in PFC3D, and b force-displacement law for a parallel bonded model (redrawn after Lisjak and Grasselli [38], Cho et al. [8] and [46]) the numerically-calibrated ceramic and frac-sand proppant models. Table 2 provides a summary of the strength parameters from the experiments and the calibrated models. According to Fig. 3, for both proppant types, the stress-strain curves from experiments and the DEM model follow a similar trend. Table 2 shows that the differences in the maximum stress (crush strength) and Young's modulus for ceramic (3.48% and 3.73%, respectively), and frac-sand (0.34% and 1.81%, respectively) are minor. These results indicate that calibrated proppant models accurately mimic real proppant behaviour when subjected to increasing stress under diametrical compression conditions. According to Fig. 3a, the stress-strain behaviour of experimentally-tested single proppants inclines downwards during the initial loading stage, which cannot be monitored during the stress-strain response of the numerically-calibrated ceramic proppant. This is a common phenomenon often experienced during experimental testing due to the initial contact between the specimen and the loading plate (De Silva et al. [49]). In Fig. 3a, the accuracy of the calibrated model is represented by shifting and redrawing the projection of the stress-strain response of the numerically calibrated proppant on the elastic region of the experimental results. In Fig. 3b, the stress-strain behaviour of experimentally-tested frac-sand proppant inclines downwards from the stress-strain response of the numerically calibrated proppant during the latter part of the loading stage. This deviation in experimental result can be explained as a result of the coupled effect of grain rotation and asperity damage undergoing in frac sand proppant at higher loading levels.

Calibration of micro-parameters of rocks
Five different formation rock types were modelled in PFC3D to understand proppant interactions with different formation rocks. The rock types were dry-intact siltstone, 30-day water-saturated siltstone, 30-day NaCl (25%)-saturated siltstone, 30-day oil-saturated siltstone, and dry-intact granite. The densities and porosities for each rock type were measured, and unconfined compression tests were conducted using cylindrical cores of 4 mm in diameter and 8 mm high (similar dimensions were used to calibrate the rocks in PFC3D). All the laboratory compression tests were performed using a 50 kN InstronDX in the Deep Earth Energy Research Laboratory (DEERL) of the Civil Engineering Department at Monash University, Australia.
Similar to the calibration process for proppants, the micro-parameters for the rock DEM models were iteratively varied until the simulations reproduced the laboratory unconfined compression test results. The rocks were modelled as a homogenous agglomerate of particles in a cylindrical shape with dimensions equivalent to the laboratory-tested rock samples (Fig. 4). A uniform particle size of 0.08 mm was selected, using the same technique as was used to identify proppant particle size. A unconfined compression test on the rock models was conducted at a slow axial loading rate of 0.08 ms -1 to ensure a quasistatic loading environment. The appropriate micro-parameters for the individual calibrated rock models were found in a trial-and-error process until the stress-strain curve from the numerical model matched the laboratory-tested UCS results. The micro-parameters of the rock specimens are illustrated in Table 3. Figure 5 illustrates the laboratory and DEM stress-strain responses for each rock type, and Table 4 provides a summary of the results in Fig. 5ae. Table 4 portray excellent agreement between the     by considering the impact of the type of proppant, the type of formation rock, and effect of realistic reservoir environments. A series of DEM simulations using calibrated proppant-rock models was performed to analyse the response of a single proppant placed on the rock surface and subjected to increased stress environments. The set-up of the numerial model in the DEM simulations is shown in Fig. 6. To reduce the computational cost, the thickness of the simulated rock was limited to 2 mm.

Single proppant interaction with sedimentary-based siltstone formation
3.1.1 Crushing and embedment response of proppants with increasing stress levels In this section, the effect of proppant type on interactions with dry-intact siltstone sedimentary-based formations is analysed. Figures 7a and b represents the variation of compressive stress with numerical time-steps for ceramic and frac-sand proppants, respectively. In this analysis, numerical time-steps were used as the variable for the x-axis, instead of axial displacement/strain, to provide better comparative visualisation of the variation of compressive stress. In addition, to accurately study the mechanical response of individual ceramic and frac-sand proppants, the embedment induced as a result of increasing stresses was calculated using Eq. 1, and the results are presented in Fig. 7c. In Eq. 1, the proppant deformation and rock deformation refer to distinct deformation undergoing within proppant and rock, respectively, as a result of the increasing stresses. Thus, to calculate these two parameters, numerical simulation tests were performed on the individual proppant and rock, separately and their respective deformations were recorded with increasing stresses. Then, proppant embedment was calculated by deducting respective proppant deformation and rock deformation from the total deformation at each loading level (total deformation is the deformation we record when conducting DEM simulation of a single proppant-rock test). Similar embedment analysis was adopted by [67] to accurately detect the embedment with load increments.
According to Fig. 7a and b, the ceramic proppant shows a slow stress growth with numerical time steps. However, the frac-sand proppant reaches the failure stress rather rapidly. This behaviour of ceramic proppants is mainly due to the significant plastic deformation which occurs as a result of the dominant proppant embedment mechanism. However, in frac-sand proppant, the stress response is more linearly variable with numerical time-steps due to the dominant crushing process of frac-sand proppant when interacting with siltstone formations at increasing stress levels. According to the results depicted in Fig. 7c, both proppants tend to show an increase in proppant embedment with stress increment. In addition, the results reveal that the ceramic proppant shows a significantly high proppant embedment at each stress level in comparison with fracsand proppant. According to Fig. 7c, at stress levels of 10, 20, 30, 40, and 50 MPa, ceramic proppant shows a  higher embedment-stress gradient can be observed at lower stress levels; however, the rate of proppant embedment declines at higher stresses. Although frac-sand shows significantly lower proppant embedment than ceramic proppants, frac-sand proppants have no potential to withstand any compressive stress beyond 35 MPa, as the frac-sand proppants undergo a significant amount of crushing. Nevertheless, no crushing of ceramic proppant is experienced at increasing stress levels upon interaction with the dryintact siltstone rock formation. However, the following section provides a detailed analysis of proppant-rock fracture initiation, fracture propagation, and fracture coalescence mechanisms when undergoing sand proppant crushing and ceramic proppant embedment while capturing proppant-rock deformation. These findings are consistent with the numerical findings presented by Zheng and Tannant [67] and Zheng et al. [67].

Formation rock damage and fracture evolution during ceramic proppant embedment
The DEM simulation of a single proppant-rock response to increasing stress levels provides a novel approach to understanding the granular level fracture evolution mechanism during proppant embedment. The evolution of the proppant embedment of ceramic proppants with increasing stress conditions is presented in Fig. 8. In the DEM simulation, a crack or fracture visualised as a disc is generated when the bonding between two individual particles breaks.
Here, tensile cracks occurring in the proppant are shown by a red disc, while tensile and shear cracks occurring within rock are shown by blue and pink discs, respectively. At minor stress levels (0.6 MPa), ceramic proppants show the initiation of a small number of tensile cracks at the contact location with the top loading platen. In addition, the generation of tensile dominant fractures can be observed in the siltstone formation, which initiates from its contact interface with the proppant. With the increment of the stress level to 10 MPa, the embedment of ceramic proppant in the siltstone rock causes rock particle spalling, which is typical when the rock surface is indented with high-strength proppants. In addition to fracture aperture reduction due to proppant embedment, the phenomenon of fines generation as a result of rock particle spalling might further reduce fracture conductivity by blocking the path of oil/gas flow [25]. With the increase in ceramic proppant embedment with the increment of stress-levels, the tensile dominant fractures initiated at the proppant-rock contact interface tend to propagate gradually. As illustrated at the fracture evolution figure at 50 MPa, a number of fractures have propagated from the proppant-indented location. Importantly, beyond 50 MPa stress levels, several major fractures coalesce, causing significant damage to the proppant-embedded rock formation. Although ceramic proppant shows an embedment of 323 lm at only 50 MPa, beyond this stress level, the rock surface cannot take any further increase in loading due to the significant damage caused to the rock formation as a result of secondary fracture initiation, propagation and coalescence. Upon reaching this damage stress level (beyond 50 MPa in this instance), ceramic proppant undergoes further embedment at a higher rate until it fully indents into the rock surface. This can be clearly observed in Fig. 7c, as upon reaching the 50 MPa stress level, the embedment-stress gradient of ceramic proppant suddenly increases. Importantly, the initiation and propagation of ceramic proppant embedment-induced secondary fractures in siltstone formations have been revealed in our previous studies [5], 5]. Some of the experimental findings obtained from Bandara et al. [5] are presented in Fig. 9. Figure 9a shows a computed-tomography scanned reconstructed image of secondary fracture initiated and propagated at a proppant-embedded location. Furthermore, Fig. 9b illustrates a 3D-profilometer scanned image at a ceramic proppant-embedded location, where a fracture has initiated upon partial embedment of ceramic proppant in a siltstone formation. The process of secondary fracture initiation and propagation as a result of proppant embedment is mainly due to the layered crystal structural orientation of the kaolinite minerals of siltstone. It has been identified that during the process of proppant embedment, existing hydrogen and van der Waals bonds of layered kaolinite minerals of siltstone specimens collapse and generate fractures [56,41]. These findings further validate the results obtained in the current DEM simulation of the single ceramic proppant-siltstone rock response to increasing stress levels. Moreover, the findings in Fig. 8 reveal that during the ceramic proppant-siltstone rock response to increasing stresses, ceramic proppant does not show a tendency to undergo any crushing. However, small numbers of tensile fractures can be observed near the contact interface between the top-loading platen and the proppant.

Fracture evolution during frac-sand proppant crushing and embedment
In addition, Fig. 8   Granite is the most abundant rock type in geothermal formations [25], which can be enhanced with hydraulic fracturing with proppants. Therefore, in this section, the behaviour of ceramic and sand proppants interacting with granite formations under the influence of increasing stress conditions is analysed in detail. Figure 10a and b provides the single proppant compressive stress response with increasing numerical time-steps for ceramic and frac-sand proppants, respectively. Similar to the findings observed when interacting with siltstone formations, slower stress growth in ceramic proppants is observed compared with frac-sand proppants when proppants are subjected to compression in granite formations. Importantly, according to Fig. 10a, with the increase in loading, ceramic proppant shows a number of additional peak stress points upon reaching its preliminary failure at 123 MPa. According to the results, the ceramic proppant shows initial peak stress at 123 MPa and suddenly drops to 3 MPa, followed by another stress increment curve. In addition, the second load increment curve of the ceramic proppant follows a sawtooth behaviour in the proppant stress-numerical steps response. This is mainly due to the breaking of asperities or the disintegration of the ceramic proppant during the loading process. Interestingly, even after preliminary failure at 123 MPa, ceramic proppant is subjected to further strength increment up to 143 MPa, later followed by full crushing of the proppant. However, for the frac-sand proppant, the loading curve follows a rapid and linear stress response followed by a rapid and catastrophic crushing process at a very low compressive stress of 40 MPa (Fig. 10b). Although both ceramic and frac-sand proppants show a dominant proppant crushing mechanism upon interaction with granite formations, proppant embedment also occurs simultaneously with increasing stress levels.

Fracture evolution during ceramic proppant crushing and embedment
In this section, the mechanisms of ceramic proppant embedment and crushing, fracture generation, propagation, and coalescence during proppant embedment and crushing under increasing stress conditions are discussed. Figure 11 shows the DEM-simulated fracture evolution results for ceramic proppant interaction with the granite formation. According to Fig. 11, during the initial loading stage, tensile fractures appear near the proppant-top loading platen contact interface. In addition, the gradual and slow indentation of ceramic proppant into the granite formation during the initial loading phase induces shear fractures near the proppant-rock contact area. The generation of shear fractures on the granite formation rock upon the interaction of ceramic proppant may be due to ceramic proppant rotation and sliding on the stiff granite surface. This ceramic proppant behaviour of sliding on stiff granite formations can be clearly seen in the first row of Fig. 11.  to significant asperity changes within the proppant as well as within the rock surface. To accurately visualise the asperity damage in the proppant assembly, the DEM results of the fracture propagation pattern of the proppant are presented in the third and fourth rows of Fig. 11. According to the results, upon reaching 100 MPa, the cluster of particles starts to fragment. Further increment of loading up to 120 MPa causes additional tensile fractures to initiate from the top and bottom contact regions of the ceramic proppant. Importantly, upon reaching 123 MPa, rock particle spalling from the granite surface can be detected. Similar to the findings described in Sect. 3.1.2, with the initiation of particle spalling from the rock surface, ceramic proppants tend to embed into the formation at a higher rate, eventually causing significant damage to the granite near the contact with the ceramic proppant. The sudden unloading phase identified after reaching 123 MPa in Fig. 10a is due to the initiation of particle spalling from the granite surface. However, with the further increase in loading, the phenomenon of particle spalling increases substantially, leading to a rapid increment in the rate of proppant embedment, as detected in the second loading curve. In addition, simultaneously to the proppant embedment process, a significant increase in proppant fragmentation can be observed during the second loading curve, as illustrated in the fourth row of Fig. 11. Upon reaching 100 MPa during the second loading curve, a substantial number of tensile-dominant fractures can be seen near the top and bottom contact points of the ceramic proppant. However, further increase in loading up to 140 MPa permits the sudden propagation of distinct vertical split fractures converging from the top and bottom of the ceramic proppant contact points, making the ceramic proppants mechanically unstable to withstand any additional load.

Fracture evolution during frac-sand proppant crushing
As described in Sect. 3.2.1, frac-sand proppants undergo a dominant crushing mechanism in granite formations when exposed to elevated stress conditions, which is shown in Fig. 12. According to the results, upon initial compression of sand proppant, stress concentrates at the proppant contact points with top loading platen and rock that act as flaw initiation zones [11]. As a result, a significant number of tensile-dominant fractures initiate within the sand proppant assembly, even at low stress of 20 MPa. Further increase in stress to 30 MPa results in major diametrical cracks propagating and converging from the top to the bottom of the sand proppant assembly. Moreover, as the sand approaches its failure strength, sand granules start to form fractures with the rapid propagation of multiple cracks. On the verge of sand proppant failure, two major fractures can be observed almost perpendicular to each other, with additional micro-cracks coalesced from the two tensile fracture planes. A similar fracturing mechanism was identified by Zhao et al. [64] by analysing the single sand fracture mechanism using X-ray micro-tomography (see Fig. 13). Note that the current DEM simulations were performed based on the assumption that sand proppants are perfectly spherical. The realistic fracturing mechanism of non-spherical sand particles may slightly differ in terms of the number of fragments created.

Impact of rock-saturation condition on mechanism of proppant embedment
According to the results reported in Sects. 3.1.1 and 3.2.1, the occurrence of proppant embedment is noteworthy even in dry siltstone and granite formations. In realistic hydraulic fracturing reservoir conditions, formation rocks are often exposed to various types of reservoir and fracture fluids over a long period. Therefore, in such circumstances, the impact of proppant embedment may be even critical and have a noteworthy impact on oil/gas recovery from sedimentary-based unconventional reservoir formations. Therefore, this section quantifies the effect of ceramic proppant embedment in siltstone formations exposed to different saturation conditions (30 days water-saturated siltstone, 30 days NaCl (25%) saturated siltstone, 30 days oil-saturated siltstone). DEM simulations of compressing a single proppant on a rock face were conducted to understand the influence of reservoir/fracture fluid saturation on proppant embedment. The simulation results under different saturation conditions are compared with the simulation results of dry-intact siltstone specimens with ceramic proppant. In this section, the type of proppants is limited to ceramic, as the effect of proppant embedment becomes critical when using ceramic proppant compared with lessstiff sand proppant. The embedment into formation rock by a single proppant was calculated using Eq. 1. Figure 14a represents the embedment obtained from DEM simulations for different saturated siltstone formations with increasing stress levels. It is evident that regardless of the saturation condition, fracture/reservoir fluid saturation of siltstone formation rock has a significant impact on the increment of proppant embedment compared with dry-intact siltstone formation. A maximum ceramic proppant embedment of 613 lm was detected for watersaturated siltstone formation rock at 50 MPa stress level, while, under similar stress conditions, dry-intact proppant showed an embedment of only 323 lm. In addition, for both dry-intact and water-saturated siltstone formations, maximum stress-level to withstand embedment was 50 MPa. For 25% NaCl saturated and oil-saturated siltstone specimens, the maximum stress level that could withstand embedment further reduced to 25 and 20 MPa, respectively. Beyond the maximum stress levels, significant damage was experienced by the formation rock due to ceramic proppant embedment, and the formation rock is no longer able to withstand any further increase in stress. Proppant embedment in dry-intact siltstone formation shows a reduction in its increment rate upon reaching 10 MPa stress. In contrast, for the water-saturated, 25% NaCl-saturated and oil-saturated siltstone specimens, a higher incremental rate in proppant embedment can be detected with increasing stress levels. For instance, even at a very low confining stress of 20 MPa, high proppant embedment amounts of 323.9 (33% increment compared to embedment in dry-intact siltstone at 20 MPa), 426.1 (74% increment), and 625.7 lm (155.6% increment) are shown by water-saturated, 25% NaCl saturated, and oil-saturated siltstone specimens, respectively.
During hydraulic fracturing, various forms of fracturing fluids can be injected to rock formations. These fluids   [20,54]. The interaction of formation rocks with injected fracture fluids can result in the softening of fracture surfaces and weakening of formation rock due to the reduction in the rock's Young's modulus [33]. Importantly, the effect of rock's Young's modulus on proppant embedment has been subjected to extensive study by a number of authors [2,11,67]. According to these authors, the decrease in Young's modulus of a formation leads to a substantial increase in proppant embedment. Therefore, in order to accurately quantify proppant embedment with fracture fluid interaction in the current DEM analysis, Young's modulus variations in calibrated siltstone assemblies are presented in Fig. 14b. In addition, the mechanical strength variation (unconfined compressive strength) of each calibrated saturated-siltstone formation is illustrated in the same figure. According to the results, the reduction in the Young's modulus and the UCS of the calibrated saturated-siltstone show an inverse correlation to ceramic proppant embedment, consistent with the findings of Alramahi and Sundberg [2] and Zheng et al. [67]. In reality, siltstone saturation with different fracture/reservoir fluids can lead to several microstructural and mineralogical alterations in rock specimens. The experimentally-tested siltstone specimens used for the DEM calibration process in the present study have mainly kaolinite and quartz-rich mineralogical compositions. In addition, the experimentally-tested siltstone specimens have higher volume pore structures with porosities ranging up to 20%. Therefore, the saturation of siltstone specimens causes significant strength degradation to its microstructure, as clay minerals and cement materials interact with surrounding water bodies by means of absorption, ion exchange, and swelling mechanisms [63]. In addition, the existence of a large volume pore structure within siltstone specimens further accelerates fracture surface softening. Furthermore, the interaction of NaCl with the kaolinite and quartz minerals of siltstone leads to mineral corrosion, causing formation softening and strength weakening [16,53]. Similarly, with the interaction of oil and siltstone, the organic compounds present in oil reduce siltstone's mechanical stability due to chemical weathering [6].

Proppant embedment-induced formation rock damage
In this section, we evaluate proppant embedment-induced formation rock damage in saturated siltstone specimens. Figure 15 shows the DEM results of fracture damage in calibrated siltstone assemblies at their respective maximum failure stress. In addition, the tensile and shear fractures generated during proppant embedment-induced formation damage are separately recorded. Figure 16 represents the variation of tensile and shear fracture counts produced by DEM-calibrated rock assemblies under different saturation conditions. Similar to the findings reported in the previous sections, significant formation damage occurs in the process of proppant embedment, and the extent of formation damage is greater in saturated siltstone specimens than dry siltstone specimens. According to Fig. 15, in the DEM simulation with dry siltstone, proppant embedment at the  directions at the onset of the respective failure stress. In addition, the fractures propagated in multiple directions coalesce from the proppant-indented location, forming a fracture network within all the saturated siltstone specimens. Therefore, it is evident that asperity damage in saturated siltstone specimens is substantially greater, which can negatively affect oil/gas recovery. Moreover, according to Fig. 16, the total induced fracture count increases to 7280, 6738, 9129 in water-saturated, 25% NaCl saturated, and oil-saturated specimens, respectively, compared with dry siltstone specimens (4503 fracture counts). The higher fracture counts portrayed by oil-saturated siltstone specimens correlate to the higher embedment shown by the ceramic proppant even at a very low-stress level. In addition, the fracture mechanism illustrated by all the simulated siltstone specimens shows a tensile-dominant fracturing mode. The small number of shear fractures shown by all the specimens may result from frictional slip between fracture surfaces during fracture propagation and coalescence.

Mechanical response of proppant pack to increasing stress levels
The arrangement of proppants within fractures after proppant injection with fracture fluid is a complex process heavily dependent on many dynamics, including particle type, morphology, geometry, existing inter-particle forces and particle-environmental forces (e.g. particle-fluid interactions, buoyancy, drag and Magnus lift forces) [27]. However, it is also crucial to study the behaviour of proppants within a proppant assembly or a pack and its response to stress conditions as it plays a significant role in ultimate oil/gas recovery. Therefore, this section provides a comprehensive analysis of proppant pack behaviour and its response to increasing stress.
To enable the accurate analysis of the effect of proppant concentration on the mechanisms of proppant pack embedment, crushing, and proppant re-arrangement, DEM models of four different proppant arrangements: partialmonolayer proppant arrangement (the results from Sect. 3 were obtained for purposes of comparison), uniformmonolayer proppant arrangement, uniform-multilayer proppant arrangement with two proppant layers, uniformmultilayer proppant arrangement with three proppant layers were simulated in PFC3D (see Fig. 17). Each proppant layer consists of three by three proppant distributions (a total of nine proppants). Moreover, in order to prepare two and three layers of proppants, each proppant layer is placed one above the other in a cubic particle arrangement. Importantly, to minimise computational time, the particle size in bonded particle agglomerates representing proppant (ceramic and frac-sand) was increased to 0.09 mm.

Effect of proppant concentration on proppant pack embedment
Proppant embedment = Total deformation ÀProppant deformation -Rock deformation The effect of proppant concentration on proppant embedment is subjected to analysis in this section. According to Sect. 3, proppant embedment is the dominant proppant degradation mechanism for three situations, including ceramic proppants embedding to siltstone, ceramic proppants embedding to granite, and sand proppant embedding to siltstone formation. Hence, these three combinations of proppant-rock models were used to further study the proppant embedment mechanisms under different proppant mass concentrations. The proppant embedment variations with compressive stress for ceramic proppantsiltstone, ceramic proppant-granite, and sand proppantsiltstone formations are illustrated in Fig. 18a, b, and c, respectively. Here, the proppant embedment of the DEM simulations was calculated using Eq. 2. A similar method adopted in Eq. 1 was utilised to calculate proppant deformation and rock deformation in Eq. 2 as well. However, an additional parameter, proppant pack deformation (re-arrangement/deformation occurring within proppant packs when subjected to increasing stresses) was calculated for Based on the DEM simulation results, an increase in proppant concentration within a fracture appears to reduce proppant embedment. For instance, at a compressive stress of 30 MPa, changing the proppant concentration from a partial monolayer to uniform monolayer (1 layer), two uniform layers, and three uniform layers of proppants in ceramic-siltstone assemblies cause an embedment change from 256.2 to 22.1 (91.4% reduction), to 16.4 (93.6% reduction), and to 10.8 lm (95.8% reduction), respectively. In addition, a substantial drop in embedment is found for the simulation of ceramic proppant embedding to granite at 40 MPa. Embedment changes from 62.7 to 8.8 (85.9% reduction) and to 4.7 lm (92.4% reduction) when the proppant concentration varies from partial monolayer to uniform monolayer (1 layer) and to two uniform layers, respectively. Similarly, in the simulation with sand proppant and siltstone at 25 MPa, embedment varies from 51.7 to 15.8 (69.4%), and to 2.9 lm (94.3%) when partial monolayer of proppants shifts to uniform monolayer (one layer) and two uniform layers, respectively. For both ceramic-granite and sand-siltstone simulations, no proppant embedment is shown at any stress level when the fracture consists of a uniform multilayer with three layers of proppants. Therefore, the results confirm that a higher proppant concentration within fractures greatly contributes to the minimisation of the effect of proppant embedment. Importantly, the DEM results of the current analysis agree well with the experimental studies conducted by Bandara et al. [6], Wen et al. [55], and Volk et al. [52], who experimentally confirmed that proppant embedment is drastically curtailed with the increase in proppant concentration. This is primarily due to the higher single-point loading acting on the proppant-rock contact interface in partial and uniform monolayer proppant arrangements compared with multilayer proppant packing arrangements. To accurately visualise this, the force chain between the proppant and rock interface of the DEM simulation of the uniform monolayer (1 layer) and uniform multilayer arrangements at 150 N are analysed. Here, the centre proppant in the bottom proppant layer is represented in Fig. 19. According to Fig. 19, it is evident that larger contact forces are visible at the proppant-rock interface of proppants in the uniform-monolayer proppant arrangement

Effect of proppant concentration on proppant pack crushing
In this section, an analysis of the influence of proppant concentration on the proppant pack crushing effect is presented. As Sect. 3 indicated, the crushing phenomenon is the most dominant proppant degradation mechanism in the sand proppant-granite assembly. Hence, a sand proppantgranite DEM assembly is utilised to analyse the influence of proppant concentration on proppant crushing. Figure  and Herskovits et al. [24], and these researchers further confirmed the resistance of proppant crushing in multilayer proppant packs compared to monolayer proppant arrangements.

Effect of proppant concentration on proppant pack re-arrangement
In addition to proppant crushing, with an increase in the number of proppant layers, proppant packs tend to undergo a significant amount of proppant re-arrangement and elastic and plastic proppant pack deformation with increasing stress conditions [62]. The DEM results of sand proppant re-arrangement with increasing stress levels are illustrated in Fig. 24. Generally, with increased stress, proppant packs undergo a higher strain hardening process due to proppant re-arrangement into voids, resulting in a denser pack arrangement, ultimately increasing proppant-proppant contacts and limiting proppant motion [43]. According to As proppants are initially modelled in a cubic arrangement, proppants maintain their stable cubic arrangement at lower stress levels. However, with the increase in stress levels, proppants in multilayer proppant packs transfer to a more stable rhombohedral proppant pack arrangement. However, prior to shifting to the final stable rhombohedral arrangement, proppants are likely to undergo a transitional position analogous to a hexagonal arrangement. Similar findings have been reported for spherical granular packs subjected to stress conditions by Li et al. [34]. Importantly, all these findings are based on the assumption that all the proppants are spherical and equal in size. However, in reality, proppant arrangement may vary compared with the present DEM analysis, as proppants often have many imperfections and irregularities [67].

Discussion on proppant selection for different rock formations/closure stress
In general, shale and siltstone formations are the most abundant sedimentary-based rock formations which contain trapped unconventional oil and gas resources. In addition, granite formations often exist as the formation rock for geothermal resources. Therefore, for the recovery of unconventional oil/gas and geothermal energy, attention must be given to the behaviour of proppants in tight rock formations if production recovery targets are to be met. The series of DEM proppant-rock numerical simulations in the present study gives detailed insights into the mechanisms of proppant crushing, embedment, and re-arrangement in fractured reservoir environments. According to the results of this study, the selection of an appropriate proppant type based on the type of reservoir formation plays an important role in determining the extent of proppant crushing and embedment in the fractures. The results reveal  Importantly, the DEM analysis reported here shows that ceramic proppants have the potential to withstand high in-situ stresses up to 140 MPa without undergoing full crushing. Furthermore, the current study reveals that ceramic proppants often undergo very high embedment in formation rocks on exposure to reservoir and fracture fluid saturation. However, in realistic reservoir conditions, the mechanism of proppant crushing and embedment may be severely intensified due to consistent proppant-rock interaction and saturation with various fracture/reservoir fluids, high temperatures followed by cool conditions, and high closure stresses. Therefore, it is crucial to take these factors into consideration in the synthesis of proppants and the selection of appropriate proppant types prior to proppant injection. Importantly, the DEM analysis indicates the importance of higher proppant concentrations in fractures owing to their ability to minimise proppant crushing and embedment. The achievement of higher proppant concentrations in fractures with better proppant distribution patterns (preferably proppant packs with multilayer packing arrangements) can be managed by the process of graded proppant injection. The injection of finer proppants into fractures followed by larger proppants allows finer particles to initially penetrate into deep fractures, while larger proppants assist in forming a multilayer proppant packing arrangement with better proppant coverage. The mechanism of graded proppant injection has gained the attention of many researchers  and this was comprehensively analysed in our previous study with the aid of experimental approaches [4].

Conclusions
The mechanisms of proppant crushing and embedment in formation rocks have been subjected to comprehensive analysis with the aid of a numerically-simulated DEM approach. The influence of proppant type, type of formation rock, and reservoir/fracture fluid saturation condition on crushing, embedment, and proppant rock damage were extensively studied using a single proppant-rock model in PFC3D. In addition, the proppant crushing, embedment, and re-arrangement undergone in multilayer proppant packs with increasing stress levels were evaluated and discussed. The following conclusions can be drawn: • The presence of a multilayer proppant pack leads to lower proppant crushing effects as in-situ stresses tend to redistribute and reduce the single-point loading on a single proppant. In addition, with the increment of stress, proppant packs tend to undergo a greater strain hardening process due to proppant re-arrangement into voids.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visithttp://creativecommons. org/licenses/by/4.0/.