Experimental and Numerical Investigation of Characteristics of Highly Heterogeneous Rock Mechanical Responses in Tight Sandy Conglomerate Reservoir Rock Under Tri-axial Compression

Due to the use of horizontal wells and hydraulic fracturing, commercial tight oil production from some tight sandy conglomerate reservoirs has been achieved. Since the widely distributed gravels in the sandy matrix in conglomerate reservoir rocks are harder than the matrix, the rock mechanical response in conglomerates under compression is highly heterogeneous. This increases the complexity of understanding the hydraulic fracturing behaviors in conglomerate reservoirs. Previous tri-axial compression tests provided the stress-strain relationships of conglomerate samples as a whole, and the stress and strain in the gravels and in the sandy matrix were not investigated due to the limitation of the compression test lab. This study presents tri-axial test results for a conglomerate sample cored from a reservoir that has been economically developed. Lab results are then used to calibrate the numerical model for the simulation of the tri-axial compression process. Numerical results indicate that the elastic modulus and size of gravels have significant impacts on the axial stresses and axial strains in the conglomerate. Stress concentrations are observed in gravels due to the heterogeneous mechanical properties in the conglomerate. The reorientation of the maximum horizontal principal stress is quantified to study the mechanisms of the interaction types between hydraulic fractures and gravels embedded in the tight sandy matrix.


INTRODUCTION
Unconventional oil and gas resources in low permeability reservoirs have attracted tremendous attention in the upstream in the petroleum industry (Xie et al., 2020;Zhao et al., 2020). Due to the low permeability in such reservoirs, it usually requires horizontal wells with hydraulic fracturing technologies to obtain commercial production (Zhang et al., 2017;Tang et al., 2018;Cheng et al., 2019;Dong et al., 2021). Recently, the high potential of hydrocarbon production in tight sandy conglomerate reservoirs has been proved by large-scale field production involving horizontal wells with hydraulic fracturing in several locations such as the Junggar Basin and Songliao Basin in China (Feng et al., 2013;Xiang and Zhang 2015). However, due to the presence of gravels in the tight sandy matrix, the rock mechanical characteristics are rather heterogeneous. Since the hydraulic fracturing quality is closely related to the complex mechanical behaviors in conglomerates, it is important to understand how the conglomerates deform under the compression induced by reservoir stimulation processes.
Experimental tests are widely used in rock mechanical analyses. In such tests, the exertion of confining pressure and axial stress can establish the in-situ stress conditions in the subsurface. They help to obtain essential rock mechanical parameters and constitutive relationships using rock samples, which are especially critical for the optimization of hydraulic fracturing parameters. Kluge et al. (2020) proposed a novel shear test strategy to correlate the permeability evolution and microfaults. Their test was established based on an MTS 815 tri-axial compression cell and concluded that the induced fractures increase the permeability by two-three orders of magnitude. Based on a modified tri-axial compression test with changing confining pressure and axial stress, it is found out that stress paths of loading and unloading directly govern the brittleness, elastoplasticity, and tensile and shear failure mechanisms in tight rocks (Guo et al., 2019). The anisotropy of deformation in tri-axial tests is discussed by Togashi et al. (2017), where a novel method with only a single test on one sample is proposed. Thus, it is more economical and convenient to quantify the deformation anisotropy. The advantage of this method is that there is no need to carry out multiple tests on various samples. Similarly, Aghababaei et al. (2019) also proposed a multi-stage strategy to obtain key strengths and failure behaviors using reduced numbers of samples to decrease the time and cost involved in tri-axial compression tests. However, they pointed out that the newly developed experimental method yields lower strength measurements than the single stage method. In order to understand the fracture patterns, Baumgarten and Konietzky (2013) conducted conventional single stage compression tests, multi-stage tests and continuous failure state tests. Both uni-axial and tri-axial tests were employed. They indicated that the combination of lab tests and numerical simulations helped to improve the understanding of post-failure behaviors.
The presence of gravels in a relatively homogeneous matrix leads to complexity in evaluating the rock mechanical behaviors, which have attracted attention from many researchers. Kumara et al. (2013) carried out a study on the stress and strain relationships in sand-gravel mixtures. This study provides insight into the deformation characteristics in sandy conglomerate reservoirs. They found out that the shape of gravels and the percentage of sand in the mixture both affect the stress-strain curves, and 30% of sand leads to the highest stress-strain curve indicating the highest strength. Akram et al. (2019) employed uni-axial, tri-axial, and Brazilian tests to investigate the effects of specimen size, clast size, cement matrix, and clast properties on the strength and deformation   Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 735208 patterns in conglomerate rocks. They concluded that with the increase of particle size, peak strength and Young's modulus decrease. Zhou et al. (2020) further experimentally investigated the effects of sphericity, where sphericity is negatively related to gravel content and the plasticity of the rock is positively correlated with heterogeneity. They emphasized the effect of gravel on the emergence of micro-cracks on the edges of gravels. In a rock mechanical test, Shi et al. (2013) reported elastic parameters based on acoustic measurements in 80 samples from the Triassic and Permian Formations in the conglomerate reservoirs in the Mahu Sag, Junggar Basin in Xinjiang, China. The dynamic modulus range is between 28.4 and 32.21 GPa while the dynamic Poisson's ratios are between 0.2055 and 0.2858. The correlation between P and S waves is good while the correlation between dynamic and static parameters is poor. A major importance of tri-axial tests is that they provide strengths, elastic parameters, and failure patterns for rock samples. It is widely accepted that the quality of hydraulic fracturing in tight and heterogeneous is jointly governed by stress shadows, in-situ stress, brittleness, fracability, clay and organic matter, and fracturing parameters (Guo et al., 2018a;Guo et al., 2018b;Dahi Taleghani et al., 2018;Li et al., 2018;Lecampion et al., 2018;Guo et al., 2019;Mao et al., 2020;Hou et al., 2019;Wang et al., 2019;Chen et al., 2021). To better characterize the hydraulic fracture network, tri-axial tests can provide essential measurements for rock mechanical properties (Guo et al., 2020;Duan et al., 2021). In addition, gravels in conglomerates introduce heterogeneities and increase the complexity of hydraulic fracture propagation. Some efforts were made to investigate these phenomena. Heterogeneities related to topology and geometry were discussed and it was observed that structured fracture networks have less dispersion than disordered networks, and the corresponding production performance can also be affected (Hyman and Jimenez-Martinez 2018;Zhi et al., 2021). The layered formation is deemed as a typical mechanical heterogeneity that leads to complex fracture propagation behaviors (Yue et al., 2019;Tan et al., 2020;Tan et al., 2021). For mechanical heterogeneities characterized by stratified formations, fracture tip locations, modulus values, and each layer's height percentage are the key parameters affecting hydraulic fracture paths and geometries.
Based on the literature review, it is noted that tri-axial tests have been widely studied, while the mechanical responses within the highly heterogeneous sandy conglomerates require further investigation. In this study, a novel workflow including tri-axial testing and the finite element analysis of mechanical responses within the testing sample is proposed. Then, this study provides lab-calibrated distributions of stresses and strains in the sandy matrix, in the gravels, and in the interface between gravels and matrix. The effects of gravel size, gravel spacing, and gravel stiffness on stress concentration are also quantified. The heterogeneity of mechanical responses discussed in this article provides a reference for the hydraulic fracturing interaction with gravels and insights for hydraulic fracture parameter design and optimization.

TRI-AXIAL COMPRESSION TEST
In the experimental study, a tight sandy conglomerate rock sample taken from the Lower Triassic Baikouquan Formation in Mahu Sag, Junggar Basin in northwestern China is used for the tri-axial compression test. 1 × 10 8 tons of reserves were already reported in the low-porosity and low-permeability reservoirs. Continuous and commercialized tight oil production has been achieved in this area after the use of horizontal wells and hydraulic fracturing. To further improve the recovery in this area, optimization of drilling and hydraulic fracturing parameters is required where rock mechanical tests can provide necessary data for the optimization. Due to the fandeltaic sedimentary environment in the sag, large-scale accumulation is favorable. Tight sandy matrix and gravels are the main components of the reservoir rock, and the clastic gravels are generally quartz and feldspar. Therefore, the stiffness of the gravels is higher than the tight matrix. Gravels in this area have a maximum diameter of 24.8 mm and have good sphericity. Correlations between P wave velocity and S wave velocity are good, and one established correlation  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 735208 is V s 0.3692 V p + 851.3 (Shi et al., 2018;Qin and Yang, 2019;Zhou et al., 2020;Qian et al., 2021). Figure 1 is the tight sandy conglomerate sample taken from the formation, where strong heterogeneity can be observed with gravels embedded in the matrix. In the sample, the gravel sizes range between 2 and 12 mm. The tight sandy sample is from a reservoir containing sandstones, mudstones, and conglomerates. Gavels show relatively good roundness and mainly contain feldspars and quartzes. The support mode is a combination of the matrix and gravels . The tri-axial compression test is conducted in the Rapid Triaxial Testing System GCTS RTR-1500 ( Figure 2). This platform can provide a maximum pressure of 140 MPa, a maximum compressive force of 1,500 kN, a maximum tensile force of 820 kN, and a maximum temperature of 150°C. A confining pressure of 40 MPa is used in the test with continuous increases in the axial loading. Figure 3 presents the experimental results of the relationship between deviatoric stress, radial strain, and axial strain. Before reaching the peak strength, a linear correlation is observed. The size of the core is a standard sample with a diameter of 25 mm and a length of 50 mm. The weight is 53.5 g. They are the essential parameters to be used in the following finite element analysis.

NUMERICAL STUDY
The tri-axial compression test is capable of providing key elastic parameters. However, the parameters are obtained based on the compression of the entire conglomerate while the distribution of stress and strain in the rock cannot be directly examined by the compression test. Since the heterogeneity in the conglomerate sample is strong, it is not reasonable to assume homogeneous and uniform mechanical responses to compression in the sample. Therefore, a numerical model based on finite element methods and momentum balance in the stress tensor is introduced to quantify the mechanical behaviors in the sample.

Mathematical Model
Cauchy stress tensors are used to describe the stress components in the 3D domain. The momentum balance between stress components and traction boundaries is expressed as: where σ is the stress tensor; t is the traction boundary.
Since the compression test results exhibit strong linear elasticity, linear elastic materials are used for the sandy matrix and the gravels. Based on Hooke's law: where σ 0 is the initial stress; C is the elasticity tensor; ε el is the elastic strain; ε is the total strain; ε inel is the inelastic strain.
The total strain can also be expressed as: FIGURE 6 | Schematic of the simplified geometry for conglomerate compression simulation.
FIGURE 7 | Matching between lab and modeling results for linear elasticity. Note that only linear regimes are studied.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 735208 The elasticity tensor is related to Young's modulus and Poisson's ratio as: Although heterogeneity in the rock mechanical properties is considered in this study, isotropy is used in the assumption. Therefore, a highly symmetric elasticity matrix is obtained: where E is the Young's modulus; ] is the Poisson's ratio.
In order to simulate the tri-axial compression process in the laboratory, a three-dimensional geometry is established as in Figure 4. Three types of boundary conditions are used in the geometry. The first type is the axial stress exerted at the top of the cylindrical domain. It is used to represent the time-dependent axial load in the GCTS testing system. The second type is the confining stress or confining pressure exerted radially. This is used to represent the confinement in the testing system. The third type is the fixed bottom boundary.
For the first and second types of stress boundaries, their effects are written as: where n is the normal vector and P is the nominal stress.
For the third type of fixed boundary, the effect is written as: u 0 (9) Note that the boundary tractions are time-dependent as functions of time. They represent the loading process in the testing system during the compression test.
To build up the numerical model, three-dimensional tetrahedral cells are used in the mesh. Local grid refinements

Base Case
In the base case, the mathematical model is calibrated with lab data. In the model, the time-dependent axial stress and confining stress boundaries are from the tri-axial compression test records as in Figure 5. Note that it only reports loading after the confining stress is fully increased to the target of 40 MPa and the process for the increase in loading for the confining stress is not plotted. This is used in Eqs 7, 8 in the model. An assumption is made in the simulation to reasonably simplify the process to quantify the heterogeneity. Based on the observation of the tight conglomerate core and statistical review carried out by Liu et al. (2018) in the same formation, uniformly distributed spheres as in Figure 6 are used to represent the gravels in the sandy matrix. Therefore, quantitative analyses can be better focused on the size, spacing, and elastic parameters of the gravels.
In general, four layers of gravels are uniformly distributed vertically. In each layer, five gravels are uniformed placed. In the base case, the diameter of each spherical gravel is 2 mm. In each layer, the spacing between two neighboring gravels is 4.25 mm. The vertical spacing between two neighboring layers is 8 mm.
Before detailed numerical analyses, a calibration of the parameters used in the model is carried out. After the calibration of rock mechanical parameters of the gravels and matrix, axial and radial strains from the model are used for matching purposes. In Figure 7, the matching between the lab results and the modeling results for axial strain, radial strain, and deviatoric stress is achieved. Note that only the results in the linear deformation regimes are used as this study does not consider failure mechanisms. Only results before the peak strength are used. Based on the match, the linear elastic assumption used in the modeling study can be verified.
The mechanical parameters used in the numerical model are calibrated in the matching process. In the calibration, the Young's modulus of the sandy matrix is 17 GPa; the Poisson's ratio of the sandy matrix is 0.41; the Young's modulus of the gravel is 47 GPa; the Poisson's ratio of the gravel is 0.31.
Distributions of the axial strain, the axial stress, and the orientation of the maximum principal stress in the horizontal direction in the base case are plotted. They are plotted at the end of the elastic deformation when the axial strain is at its maximum. In Figure 8, the two-dimensional distributions of axial strains in planes Z 10, 20, 30, and 40 mm are presented. The deformation in gravels is generally smaller than that in the matrix, as gravels have higher stiffnesses. Distinct strains are observed as the effects of the boundaries between gravels and the matrix are dominant. In the sandy matrix, areas between neighboring gravels also experience relatively higher axial strains compared to areas near the boundary of the rock sample. This is the stress concentration caused by non-uniform deformation in gravels and in the matrix during the tri-axial compression process. Therefore, it is found out that the stress concentration is only observed in between gravels. The strain magnitudes are higher at areas closer to the fixed bottom boundary of the rock sample. This is because the fixed bottom leads to higher changes in deformation at nearby areas, while the top and the confining boundaries are traction boundaries that allow for movements.
In Figure 9, one-dimensional distributions of axial strains along X 0, 3.125, and 6.25 mm at four X-Y planes in Figure 8 are plotted. These lines are selected as they penetrate the three gravels, zero gravels, and one gravels as it moves away from the center of the sample. Vertical lines in the one-dimensional plots represent the boundaries between gravels and the sandy matrix. Note that in Figure 9, only half of the lines are presented due to symmetry. At X 0mm, it is clear that the axial strain magnitudes are lower in the gravels than in the matrix. The strains in the gravels and between gravels are positively correlated with the distance to the bottom boundary, while the differences in strains near the boundary of the sample are not distinct except for Z 10 mm. At X 3.125mm, the strain distributions between gravels are discussed. Areas closer to the center of the sample have higher axial strains. This indicates that the stress concentration is more noticeable in areas more surrounded by gravels. At X 6.25mm, the effect of gravels on decreased axial strains is again exhibited. The changes in the strain are relatively sharp at the gravel-matrix boundary due to the changes in elastic parameters. The axial strains outside the gravels are around 0.35%. Figure 10 plots the axial stress distributions at X-Y planes with different Z values. Strong heterogeneities in axial stresses are also observed. In general, stress concentrations are located within the gravels as they are typically harder than the sandy matrix and have greater elastic moduli. The elevated axial stresses in gravels at Z 10 mm are the greatest as they are closer to the fixed bottom boundary which allows for limited space for deformation. Another observation is that the interfaces between gravels and the matrix endure decreases in axial stresses, which is caused by the sharp differences in elastic moduli between two different materials.
To better present the heterogeneity in axial stresses, Figure 11 describes the one-dimensional axial stress distributions at X 0, 3.125, and 6.25 mm at four different X-Y planes. Due to symmetry, only half of the domain is plotted. It is noticed that Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 735208 the axial stress at and near the center of the gravel in the plane tends to be greater than that in the outer gravel. Significant decreases in axial stresses are observed at the boundaries between gravels and the matrix. The magnitudes of the decreases can be as great as 15 MPa. For the one-dimensional distribution at X 3.125 mm with no presence of gravels, oscillations in stresses are observed due to the presence of gravels at nearby areas, while the general trends are flat compared to X 0 and 6.25 mm.
In addition to axial strains and axial stresses in the matrix and gravels, the orientation of the maximum horizontal stresses is quantitatively analyzed. This is of special significance as previous studies have indicated that the interaction between heterogeneous gravels and hydraulic fractures is complicated, and hydraulic fractures in tight sandy conglomerate reservoirs can penetrate the gravels, bypass the gravels, be attracted by the gravels, or be terminated by the gravels (Zhu et al., 2005). Figure 12 shows conceptual interaction types between hydraulic fractures and gravels. The underlying mechanisms are complicated due to the heterogeneity in rock mechanical parameters.
To understand the mechanisms behind the interaction types presented above, in Figure 13, the effect of compression on stress reorientation in the rock sample is studied based on numerical modeling. The compression effect here resembles the tensile failure caused by the net pressure in hydraulic fracture. To calculate the reorientation of the principal stress in the horizontal plane presented in Figure 13, Eq. 10 is used as below: where θ is the change in the orientation; τ xy is the shear stress; σ x and σ y are stresses in two directions.
Since the rock sample is under a conventional tri-axial compression state (σ 1 > σ 2 σ 3 ) with uniform confining pressure and an eternal axial load, the initial direction of S Hmax is strictly radial. However, due to the heterogeneity introduced by gravels, it is noted that at various X-Y planes, S Hmax orientations are no longer strictly radial, and highly non-uniform reorientations are observed. In general, the S Hmax directions along the X axis and the Y axis are still largely radial, which is closer to the initial condition. At the circular interfaces between gravels and the matrix, S Hmax directions FIGURE 11 | Distribution of axial stresses at X 0, 3.125, and 6.25 mm in the Y-direction in planes of Z 10, 20, 30, and 40 mm in the base case. Vertical lines denote the gravel-matrix interfaces.
FIGURE 12 | Types of interactions between gravels and hydraulic fractures in tight sandy conglomerates.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 735208 8 are observed to be circled around the gravels. This is due to the nonuniform changes in the normal stresses and shear stresses around the gravels. The most prominent shear stress changes are in sandy matrix areas near the gravels, and this is the major reason for stress reorientations and stress reversals.
Based on the non-uniform reorientation of S Hmax in Figure 13, the different types of interactions between gravels and hydraulic fracturing propagations in the conglomerate can be explained. A hydraulic fracture initially propagating in the X-direction is used in the analysis as in Figure 14.
Firstly, for the penetration type, when the hydraulic fracture propagates in the direction of S Hmax since the directions of S Hmax in the gravel and in the matrix along the path do not change, the gravel-matrix interface does not change the fracture propagation direction. Although the hydraulic fracture propagates orthogonally with the direction of S Hmax in the gravel, the totally reversed S Hmax usually does not change the fracture path. This phenomenon of the effect of totally reversed S Hmax on fracture path is also verified in Safari et al. (2017). Secondly, for the bypassing type, the intersection between the gravel-matrix boundary and the hydraulic fracture has a S Hmax orientation that is nearly tangential with regard to the gravel-matrix boundary. Thus, the initial fracture propagation is diverted in the directions tangential to the gravel-matrix boundary. Except for the yellow fracture path in the plot, the fracture can also be captured by the gravel-matrix boundary and then be diverted by the reoriented S Hmax and move past the gravel. Thirdly, for the attraction type, the initial hydraulic fracture is governed by the S Hmax pointing in the radial direction. Note that the S Hmax pointing in the radial direction between the outer gravels is not affected by the existence of gravels. After the initial fracture propagation contacts the gravel-matrix boundary, the tangential S Hmax on the boundary diverts the fracture along the boundary. This explains the scenarios where a fracture propagating near a gravel is attracted to the gravel. Fourthly, it is also possible that a hydraulic fracture is terminated by the gravelmatrix boundary. This cannot solely be explained by the reorientation of S Hmax , as the different strengths between gravels and the sandy matrix can also affect the propagation of fractures. Since gravels usually have higher strengths, it is possible that the fracture can penetrate through the sandy matrix while the net pressure within the fracture cannot provide enough tensile traction to induce failures in the gravels. Also, modeling results indicate that there are stress concentrations in gravels, which makes it more difficult for the hydraulic fracture to propagate within the gravels with higher stresses. Thus, the propagation path of the hydraulic fracture is terminated by the interface between the gravel and the matrix.

Effects of Size and Elastic Modulus of Gravels on Heterogeneous Mechanical Responses
In the base case, a gravel elastic modulus of 27 GPa and a gravel radius of 2 mm are employed. To investigate the effects of gravel elastic modulus and gravel size on the deformation of the conglomerate sample, two sets of sensitivity analyses are carried out to show how the changes in elastic modulus and size of the spherical gravels affect the stress distributions in the conglomerate. In the analysis of elastic modulus of gravels, elastic modulus values of 35 and 19 GPa are used. In the analysis of gravel size, radii of 1 and 3 mm are used.  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 735208 In Figure 15, axial stress distributions in the Y-direction for X 0 mm at Z 40 mm are plotted. The effect of gravel elastic modulus and gravel size is exhibited. Based on the results with various gravel elastic moduli on the left, it is noted that the interfaces between gravels and the sandy matrix still denote the sharps changes in axial stresses. In general, a higher gravel elastic modulus leads to higher axial stress magnitudes in the gravels. However, when the elastic modulus is decreased to 19 GPa which is very close to the elastic modulus of the sandy matrix, axial stresses in gravels are no longer elevated. In contrast, increases in axial stresses in the sandy matrix between neighboring are observed. Therefore, when the elastic modulus of gravels is decreased from 35 to 19 GPa, the stress concentrations move from the gravels to the matrix between gravels. When the difference in elastic moduli between gravels and the matrix is relatively large (e.g., 35 GPa for gravels and 17 GPa for the matrix), it is harder for gravels to deform under tri-axial compression. When the difference is small (e.g., 19 GPa for gravels and 17 GPa for the matrix), it is easier for gravels to deform under tri-axial compression. However, in general, the average axial stresses in all elastic modulus simulation cases are very close. This indicates that although gravels and the matrix have their own moduli, the overall mechanical responses in the conglomerate to tri-axial compression do not change much. This is because the tri-axial compression parameters of axial loading and confining pressure are the same in all sensitivity cases.  When gravel radii of 1 mm, 2 mm, and 3 mm are simulated, it is noted that the peak axial stresses for a radius of 1 mm and a radius of 3 mm are at very similar levels of around 98 MPa in gravels. When the radius changes, the corresponding axial stress distribution drastically changes. In general, axial stresses in gravels are higher and axial stresses in the matrix are lower. The profile shapes are directly governed by sizes of spherical gravels. Regardless of the gravel size, sharp decreases in stresses are observed at the interfaces between gravels and matrices. However, there are no clear correlations between gravel radius and the axial stresses in either gravels or the matrix.
Since the stress reorientation around gravels is a key parameter to determine the interaction between hydraulic fracturing and gravels, Figure 16 plots the comparison of S Hmax reversals in the conglomerate at Z 40 mm for all the elastic modulus cases and gravel size cases. Results with elastic moduli of 35 GPa, 27 GPa, and 19 GPa indicate that the gravel modulus does not significantly alter the reorientation patterns of the maximum principal stress in the horizontal plane. Generally, S Hmax directions are circled around the gravel-matrix interfaces. For results with gravel radii of 1 mm, 2 mm, and 3 mm, the effect of gravel radius on reorientation of S Hmax directions is more pronounced. For the case of radius 3 mm, S Hmax in gravels generally points to the radial direction which is the direction caused by the loading of confining pressure. As the radius decreases from 3 to 1 mm, the radius of the circled S Hmax reorientation around gravels also decreases. Since S Hmax reorientations are circled around the gravels, changing the size of the gravels alters the reorientation patterns of S Hmax . Based on the observations in Figure 13, changing the size or modulus of gravels consequently alters the hydraulic fracturing patterns in the conglomerate.

CONCLUSION
In this study, the heterogeneous rock mechanical responses in a tight sandy conglomerate sample cored from a Lower Triassic Formation under tri-axial compression are experimentally and numerically investigated. The stress-strain relationships are first obtained from the lab. The lab data are then used to calibrate the finite element model for the simulation of the tri-axial compression process. Then, detailed numerical results for axial strain, axial stress, and S Hmax reorientation in the sandy matrix and gravels are presented and discussed.
In Conclusion 1) Before reaching peak strengths, the deformation of the tight sandy conglomerate sample exhibits linear elasticity. Consequently, the constitutive relationship in the numerical model is also determined as linear elastic.
2) Since the tri-axial compression test in the lab can only demonstrate the overall response of the conglomerate sample and cannot quantify the spatial distribution of stress and strain in the sample, it is meaningful to develop a three-dimensional finite element model with a reasonable and lab-calibrated parameterization for the quantification of the stress and strain evolutions in the heterogeneous sample. 3) In general, axial strains in the gravels are lower than those in the matrix, while axial stresses in the gravels are higher than those in the matrix. This phenomenon of stress concentration is explained by higher elastic moduli in the gravels, which makes it more difficult to deform. 4) The elastic modulus and size of gravels both affect the distribution of stresses and strains. Increased elastic moduli in gravels lead to more heterogeneous distributions of stresses and strains among the conglomerate, while their effects on S Hmax reorientation are not significant. The pattern of S Hmax reorientation is primarily governed by the size of gravels. 5) The nonuniformly reoriented S Hmax in the conglomerate is one of the reasons why hydraulic fracturing can have various interaction types with gravels. The interaction types include penetrating the gravels, bypassing the gravels, being attracted by the gravels, and being terminated by the gravels

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
BC, JJ, JL, HC, and XW contributed to the conception, design, and data interpretation for this study. XG and JL conducted the experimental and numerical study. XG wrote the draft of the manuscript. JJ, WY, and JL modified the manuscript.