Discrete Element Study on the Mechanical Response of Soft Rock Considering Water-Induced Softening Effect

: Soft rocks are prone to softening upon contact with water, and their rapid deterioration in mechanical properties is a significant cause of instability and failure soft rock masses. Besides, the macroscopic mechanical response of rocks is closely related to the mineral composition and microstructure. The purpose of this research is to consider the heterogeneity factors and softening effects, and systematically investigate the influence of confining pressure and softening time on the damage and failure characteristics of soft rocks. The Voronoi polygons generated using a built-in Voronoi diagram algorithm and contact elements (the substances with cementing capacity) of UDEC discrete element method are employed to represent the clastic grains and interfacial cemented bonding (ICB) structures in soft rock. Based on the Voronoi probabilistic method, the grain-based discrete element model (GB-DEM) considering the softening effect is established by introducing a meso-scale softening damage factor, along with a detailed calibration method for meso-scale parameters. The damage parameters such as the crack initiation threshold, the crack damage threshold, the damage degree, and the tensile and shear crack ratio are then analyzed. The study results indicate that the simulated strengths of the heterogeneous models under different water immersion time are in good agreement with the experimental results. The thresholds for crack initiation and damage, the proportions of tensile and shear cracks, and the degree of damage are positively correlated with the confining pressure. The attenuation patterns of the crack initiation threshold and damage threshold in the heterogeneous models with water immersion time are highly consistent with the meso-scale softening damage factor. The damage parameters show a trend of increasing first and then decreasing with the extension of water immersion time. The cement–cement contact elements are the main locations for crack initiation and propagation. The research outcomes have significant theoretical and practical implications for understanding and predicting the mechanical behavior of soft rocks under a water–rock interaction.


Introduction
Rock masses serve as the medium for the storage and flow of geological fluids, and often interact with water in engineering construction.These interactions occur in various projects, including dam foundations [1], reservoir slopes [2], and deep buried water diversion tunnels [3], where the rock masses are susceptible to different levels of softening due to groundwater or surface water [4].Soft rocks with argillaceous cementation, as a complex rock with intense water-rock interaction, exhibit significant water-induced softening characteristics.The significant reduction in mechanical properties due to water softening is an important cause of geological disasters in many engineering soft rock masses.Therefore, studying the softening effect of water on soft rocks has significant practical importance.
In recent years, numerous studies have investigated the effects of water on the mechanical responses of rocks.The water-induced weakening analysis mainly focuses on the influence of water content, immersion time, and wetting-drying cycles on rock strength.Many studies have shown that the mechanical parameters (compressive, tensile, shear strength, and Young's modulus) of rocks weaken with increasing water content, especially mudstone [5][6][7].In addition, the water content can reduce the crack toughness [8].Compared with water content, the relationship between the immersion time and the actual field situation is closer.The immersion time can significantly affect rock strength since the deep surrounding rock masses are immersed in water for months [9].For wetting-drying cycles (WDC), it was found that the first several WDC significantly affect the strength of rocks, but there was no change in the following cycles [10,11].In essence, natural rocks are heterogeneous, and the macroscopic failure behavior and mechanical properties of rock materials under water-rock interactions are closely related to the types of minerals (e.g., minerals with different shapes, sizes, and mechanical properties) and microstructures (e.g., pores, mineral particle skeleton, etc.) [12][13][14].The softening mechanism of water on rocks is complex, and many scholars have studied different issues using various techniques.Jiang et al. [15] investigated the evolution law of micro-fractures in mudstone during the waterrock interaction process based on CT scanning.Zhang et al. [16] studied the evolution characteristics of microstructures of red bed mudstone and sandstone under rainwater infiltration based on XRD and SEM, and found that the contact relationship of clay mineral particles changed from face-to-face to edge-to-edge and edge-to-face contact, leading to the loosening and porosity of the micro-structure.Ciantia et al. [17,18] proposed the concepts of diagenetic bonding (DG) and depositional bonding (DP) that support the macroscopic stability of rocks in order to study the impact of water on the macroscopic mechanical properties based on the polarizing microscope thin sections.However, obtaining and preparing rock samples for experiments and conducting experiments are redundant, and the microstructure of each rock sample is different, so the results of laboratory tests are almost irreproducible.
In recent years, with the rapid development of computational technology, numerical methods have emerged as a promising experimental alternative tool to study the mechanical behavior of rocks at different scales.However, numerical methods based on continuous media (e.g., FEM, FDM, BEM) are difficult to capture the initiation and propagation of cracks during the fracture process of rocks, and cannot characterize the rock heterogeneity.Discontinuous simulation methods [19,20] and continuum-discontinuum simulation methods [21,22] overcome the drawbacks of continuous medium-based methods, and have obvious advantages in characterizing the special morphology and spatial distribution of material components and in studying the impact of micro-scale heterogeneity on the macroscopic deformation and failure process of rocks.For the former, block discrete element and particle flow methods have been widely used to study the mechanical properties of soft rocks.Xie et al. [23] used the discrete element method (DEM) for the numerical study of uniaxial compression tests on layered rock, and discussed the influence of different thicknesses and dips of soft rock interlayers on the deformation and strength characteristics of rock specimens.Dong et al. [24] used 3DEC to study the toppling deformation and failure mode of soft-hard interbedded rock masses.Yang et al. [25] used a method based on the two-dimensional particle flow code to explore the mechanical behavior of mudstone specimens with different damage degrees under multi-stage triaxial compression.Zhao et al. [26] simulated the weakening effect of water on rock mechanical properties by reducing the friction coefficient and bonding strength.Potyondy [27] proposed a parallel bond stress corrosion model to simulate the stress corrosion effect of stress or water solution on rock cementation.For the latter, Indraratna et al. [28] explored the deformation characteristics of a single stone pillar and its surrounding clay based on the finite-differencediscrete-element coupling method (FDM-DEM).Wang et al. [29] established an FDM-DEM coupling model for coal, mudstone, and sandstone and explored the co-evolutionary laws of stress distribution and crack extension.Huang et al. [30] explored the failure mecha-nism of surrounding rock in soft rock tunnel excavation based on the FDM-DEM.Lisjak et al. [31] proposed an anisotropic strength model implemented in the combined finiteelement -discrete-element method (FDEM) to explore the impact of bedding planes on the mechanical response characteristics and failure characteristics of Opalinus claystone.However, the above studies rarely consider both the heterogeneity and the softening effect of water on microstructures simultaneously.
Based on the above issues, this paper aims to comprehensively consider the heterogeneity factors and softening effects, and to systematically explore the influence of confining pressure and softening time on the damage and failure characteristics of soft rocks.

Voronoi Diagram
The author [32] believes that the skeleton of soft rocks such as terrigenous clastic rocks is composed of clastic grains, cement (non-clastic minerals), and interfacial cemented bonding (ICB) structures formed by the substances with cementing capability.For numerical models, three geometric shapes are commonly used to discretize the study area, namely circular elements [33,34], Voronoi polygons [35], and triangular elements [36].According to the micro-scale image of silty mudstone under the polarizing microscope, clastic grains are often polygonal structures (Figure 1), so Voronoi polygons and contact elements are used to represent clastic grains with different properties and the contact types corresponding to ICB structures.The Voronoi diagram is a method that can generate any polygon or polyhedron in a specific area.Mathematically, given a series of seed points S i on the domain D, D ∈ R 2 , R 3 , each seed point is assigned a Voronoi cell, which is as follows: where d(•, •) represents the Euclidean distance.Voronoi diagram is a continuous polygonal structure formed by connecting the perpendicular lines obtained from the adjacent seed points randomly selected in the domain following uniform distribution Therefore, the area filled by the Voronoi polygons has no overlaps or gaps.Based on the above description, the discrete element method (DEM) is chosen in this study, which has great advantage in establishing Voronoi diagram and simulating fracture between them.

Mohr-Slip Joint Model of Discrete Element Method
DEM was proposed by Cundall [37], which regards the research object as a collection of multiple rigid or deformable units.In this study, the Mohr-slip joint model is adopted, which assumes a linear relationship between stress and displacement, as shown in Figure 2. Once the stress of the sub-contact element exceeds the strength threshold, the sub-contact element will undergo tensile or shear failure.The specific principle is as follows.In the normal direction of the sub-contact element, it is assumed that the normal stressdisplacement obeys a linear relationship, which can be expressed as: where ∆σ n is the normal stress increment, ∆u n is the normal displacement increment, and k n is the normal stiffness.In addition, the threshold of normal tensile stress is the tensile strength f t .When the normal tensile stress exceeds the tensile strength, σ n = 0.
Appl.Sci.2024, 14, x FOR PEER REVIEW 2 of 14  Similarly, the shear stress-displacement relationship in the shear process is controlled by the shear stiffness k s , which can be expressed as: where ∆τ s is the shear stress increment, and ∆u e s is the elastic part of the shear displacement increment.Then, the maximum shear stress τ max on the sub-contact element can be expressed as: where c and φ are the cohesion and internal friction angle of the sub-contact element, respectively.When the shear stress τ s ≥ τ max , the residual strength of the sub-contact element can be expressed as: where φ r is the residual internal friction angle of the sub-contact element.
According to the XRD test, the silty mudstone studied in this paper is mainly composed of quartz, feldspar, and cement (e.g., non-clastic minerals), with contents of 35%, 33%, and 32%, respectively.Therefore, the grain-based heterogeneous model should contain three components: quartz, feldspar, and cement, as well as six types of contacts: quartzquartz contact, quartz-feldspar contact, quartz-cement contact, feldspar-feldspar contact, feldspar-cement contact, and cement-cement contact.The FISH language built into Itasca software is utilized to construct an algorithm to generate the grain-based discrete element (GB-DEM) model.The distribution of material components within the model follows a uniform distribution, and the generation of Voronoi polygons uses the built-in algorithm in the UDEC software.The contents of each material in the numerical model are consistent with the experimental results, as shown in Figure 3.The size of the numerical model is consistent with the test conditions, with the standard cylindrical model having a height of 100 mm and a diameter of 50 mm; the Brazilian disk has a diameter of 50 mm and a thickness of 25 mm, according to the International Society for Rock Mechanics Suggested Methods.The range of loading rates used for uniaxial compression and Brazilian disk tests in the existing literature [35,36,38,39] is 0.01 m/s to 0.1 m/s.Since the default time step size in the UDEC software is about 10 −8 to 10 −7 s/step, a loading rate of 0.1 m/s is equivalent to loading 10 −6 to 10 −5 mm per time step.It means that loading 1 mm requires 10 5 to 10 6 time steps, indicating that the specimen is under quasi-static loading conditions.Therefore, considering computational accuracy and efficiency, the loading rate is taken as 0.01 m/s.

Calibration Procedure
In fact, meso-scale parameters in numerical models are difficult to obtain directly through experiments.Therefore, it is necessary to efficiently calibrate a set of reasonable meso-scale mechanical parameters through a standard calibration process.Some of the literature describes how to calibrate the meso-scale parameters of rock [38,[40][41][42].In summary, the calibration steps are as follows: (1) directly use the elastic modulus and Poisson's ratio obtained from experiments as the elastic constants of the block elements; (2) calibrate the elastic modulus based on the normal stiffness, and calibrate Poisson's ratio by changing the ratio of normal stiffness to shear stiffness; (3) adjust the contact cohesion, contact friction angle, and contact tensile strength to calibrate the cohesion, friction angle, and tensile strength of rock samples.

Implementation of Softening Effect
In this study, the mudstone samples were completely immersed in the water tank.At different immersion time (1 day, 4 days, 8 days, and 15 days), the samples were taken out to conduct the uniaxial compression and Brazilian splitting tests.From the perspective of the diagenetic process of soft rock, the cementation process is one where the microstructural bonding strength increases.The softening process during water immersion means that the bonding capability of the microstructure is weakened.Therefore, for numerical models, the softening process corresponds to the attenuation of contact cohesion (jcoh) and contact tensile strength (jten) of contact elements.Unlike the existing literature that uses macroscopic elastic damage factors to indirectly characterize the damage of microstructures [43], a damage variable to describe the water-induced weakening effect of ICB structure directly from the mesoscale is provided by the author as follows [32]: where d(t) is the meso-scale softening damage factor, D is the diffusion coefficient, and t denotes the time of water-rock interaction.The innovative grain-based discrete element model considering the heterogeneity and the softening effect is illustrated in Figure 4. Therefore, the water-induced softening effect under water-rock interaction can be expressed as follows: Appl.Sci.2024, 14, x FOR PEER REVIEW 4 of 14

Parameter Calibration
According to Section 2, the meso-scale parameters of the cement can be directly based on the mechanical parameters obtained in the experiments.For minerals, several literature and rock mineral monographs [39,[44][45][46] provide the basic mechanical parameters.The mechanical parameters of the three material components are summarized in Table 1.Regarding the meso-scale parameters for the six types of contacts, the contact parameters between the same components were first calibrated.In this paper, the numerical model is regarded as a virtual model containing only quartz or feldspar or cement.The contact parameters for cement-cement contact are equivalent to those of the homogeneous soft rock model.The calibration of the meso-scale contact parameters for the homogeneous model in its natural state are targeted at the average values of the mechanical parameters obtained from experiments.The uniaxial compression, Brazilian splitting, and biaxial compression numerical tests are then conducted.For quartz-quartz and feldspar-feldspar contact, only uniaxial compression numerical tests are conducted to calibrate the uniaxial compressive strength and elastic modulus.The meso-scale parameters of contact types between the same components after trialand-error process are shown in Table 2, and the simulation results are shown in Figure 5 and Table 3.According to Figure 5a,b, the failure modes under uniaxial compression and Brazilian splitting conditions are primarily dominated by single-plane shear failure and splitting tensile failure, respectively, which are in good agreement with the failure modes under experimental conditions.Due to the inherent pores and fractures in the soft rock, it often undergoes a stage of pore and fracture compaction under loading conditions.Therefore, the initial slope of the axial stress-strain curve under experimental conditions is relatively small and gradually increases, which is different from the numerical simulation results.As shown in Table 3, except for the tensile strength with an error of −9.36%, which is relatively large, the simulation errors for the other parameters are at most 2.78%, indicating that the meso-scale parameters of cement-cement contact are reasonable.Figure 5c shows that the axial peak stresses under different confining pressures are 29.99MPa, 35.5 MPa, 41.3 MPa, 53.6 MPa, and 80.9 MPa, respectively.Compared to the experimental results, the simulation errors are 2.81%, 3.71%, −7.17%, 3.45%, and 0.52%, respectively.Except for the error at a confining pressure of 2 MPa, which is greater than 5%, the rest of the errors are relatively small.Moreover, according to the simulation results of the virtual mineral model shown in Figure 5d, the uniaxial compressive strength and elastic modulus of quartz and feldspar are 224 MPa and 80.2 GPa, 181 MPa and 56 GPa, respectively, with errors of 0.90% and −0.99%, 0.56%, and 0.00%.For the contact parameters between different components, the average value of the contact parameters of two same components is taken in this paper, but the simulated uniaxial compressive strength and tensile strength are much higher than the mechanical parameters of silty mudstone under experimental conditions.Therefore, this paper simultaneously increases or decreases the same proportion factor for the contact parameters of the six contact types, and the calibrated contact parameters are shown in Table 4.
Biaxial compression tests under various confining pressures and Brazilian splitting tests were simulated based on Table 4, with the results shown in Figure 6 and Table 5.The simulation results indicate that the failure modes simulated under uniaxial compression and Brazilian splitting conditions are highly consistent with those observed under experimental conditions.As can be seen from Figure 6c, the axial peak stresses under different confining pressures are 29.77MPa, 35.91 MPa, 41.70 MPa, 53.46 MPa, and 79.56 MPa, respectively.The errors compared to the experimental results are 2.02%, 4.91%, −6.27%, 3.18%, and −1.14%, respectively.According to Table 5, the simulation error of elastic modulus is the largest, which is only −2.59%.The simulation errors for the other parameters are all less than 2%.Therefore, the mechanical parameters in Table 4 are reasonable.The crack initiation threshold (σ ci ), the crack damage threshold (σ cd ), the damage degree (D), and the tensile and shear crack ratio are used to investigate the damage characteristics of soft rock under different confining pressures and water immersion durations.Based on previous studies [47][48][49], this paper defines the stress corresponding to the first occurrence of tensile cracks as the crack initiation threshold, and the inflection point of the shear crack curve where it significantly increases as the crack damage threshold.The damage degree refers to the ratio of the length of broken contact elements to the total length of contact elements.The proportion of tensile (shear) cracks is the ratio of the number of tensile (shear) crack elements to the total number of contact elements.The number of cracks and the crack types are both identified and recorded using the built-in FISH language in UDEC software.Since the failure modes of heterogeneous models with different water immersion time are relatively similar, the influence of confining pressure on the damage and failure characteristics is investigated by taking the heterogeneous model in the natural state as an example.
Figure 7 shows the relationship curves between stress, the number of tensile and shear cracks, damage degree, and strain of heterogeneous models under different confining pressures.The results indicate that the heterogeneous models first generate tensile cracks at the initial stage of loading, with the transverse strain linearly decreasing and the volumetric strain linearly increasing, indicating some degree of damage in the rock during the early stage of loading.Moreover, the relationship curves between the number of tensile and shear cracks and the damage degree with strain fluctuate significantly due to the significant divergence in mechanical properties between different components.As loading continues, the shear crack curves will all exhibit a clear increase in slope, the number of tensile cracks continues to increase, and the rock also exhibits characteristics of nonlinear deformation.As the confining pressure continues to increase, the relationship curves between the deviatoric stress and axial strain, transverse strain, and volumetric strain become more prominent in the post-peak stage, indicating that the greater the confining pressure, the greater the inhibiting effect on the failure of rock specimen.At the approach to the peak strength, the shear crack curves gradually tend to stabilize.Although the tensile cracks continue to increase, the rate of increase slightly slows down, indicating that the main tensile and shear cracks are generated before the strength.The failure mode under low confining pressure conditions (0 MPa and 1 MPa) is mainly characterized by "X-type conjugate shear failure".As the confining pressure continues to increase, the main failure mode is single-plane shear failure, with local features observable as "X-type conjugate shear failure".According to the fitting results shown in Figure 8, the crack initiation threshold, crack damage threshold, and damage degree all show a linear increase trend with the confining pressure.The tensile and shear crack ratio also have linear positive correlation with confining pressure.

Mechanical Response Considering Softening Effect 4.2.1. Strength Characteristics
The softening effect considered in numerical simulation is essentially to attenuate the contact cohesion and contact tensile strength according to Equations ( 7) and (8).In this paper, soft rock specimens were immersed in water for 1 day, 4 days, 8 days, and 15 days, respectively.Since the attenuation degree of internal friction angles in the test is small, the contact internal friction angles of the heterogeneous models with different immersion time are consistent with the values in Table 4.The parameter calibration process described in Section 3 is used.
Figure 9 shows the stress-strain (displacement) simulation curves for heterogeneous models with different immersion time.The simulation results indicate that the post-peak characteristics of the stress-strain curves of uniaxial compression are quite evident, which is due to the presence of quartz and feldspar with very good mechanical properties in the heterogeneous model.During the loading process, the mineral composition of the model increases its resistance to external loads, resulting in a gradual decrease in the stress-strain curve after reaching the peak stress.According to Figure 10, the uniaxial compressive strengths for different immersion times of 0 days, 1 day, 4 days, 8 days, and 15 days are 29.65 MPa, 4.40 MPa, 2.18 MPa, 1.88 MPa, and 1.21 MPa, respectively, with errors compared to the average experimental values of 1.61%, −5.38%, −8.40%, −4.44%, and −5.47%, respectively.The tensile strengths for the same immersion time are 2.656 MPa, 0.332 MPa, 0.152 MPa, 0.114 MPa, and 0.085 MPa, respectively, with errors compared to the average experimental values of −0.11%, −7.78%, −5.59%, 3.64%, and 4.94%, respectively.In summary, the small errors observed in both uniaxial compressive strengths and tensile strengths support the validity of the meso-scale softening damage factor to some extent.These errors suggest that the meso-scale softening damage factor has the potential to capture the softening behavior of soft rocks reasonably well.However, the presence of discrepancies between simulated and experimental results indicates that there may be areas where the model can be further refined or enhanced to improve its accuracy in representing the mechanical response of materials.It is crucial to consider the limitations and assumptions associated with the meso-scale softening damage factor and validate its applicability across different materials, loading conditions, and degradation mechanisms.Overall, while the errors are relatively small, they serve as indicators for potential areas of improvement in the meso-scale softening damage factor.Further calibration, validation against additional experimental data, and refinement of the modeling approach could help enhance the prediction accuracy and reliability of the model, leading to a better understanding of the softening behavior of materials and improved applicability in practical engineering scenarios.

Damage Characteristics
Figure 11 illustrates the relationship between stress, damage degree, and the number of tensile and shear cracks with strain under uniaxial compression conditions for different durations of water immersion.The results indicate that the overall trends of variation in both tensile and shear cracks in the heterogeneous models are similar.The fluctuation of the shear crack curves is evident, especially at 4 and 8 days of water immersion.The tensile and shear cracks of the heterogeneous model tend to be stable near the peak strength only when immersed in water for 0 days and 1 day, which is because tensile and shear cracks may continue to occur even if the stress is small when immersed in water for a long time.In addition, Figure 12 shows the heterogeneous models with different durations of water immersion are primarily dominated by "X-type conjugate shear failure".
In order to further explore the relationship between the crack threshold and meso-scale softening damage factor shown in Equation ( 6), a factor D w was introduced to characterize the attenuation degree of the crack initiation threshold and crack damage threshold: where D w is the crack threshold damage factor in a numerical simulation, TS 0 is the crack threshold of the numerical model in its natural state, and TS t is the crack threshold of the numerical model with t days water immersion.According to Equation (9), the attenuation curves of the crack initiation threshold and damage threshold and the relationship curves of damage parameters with immersion time of the heterogeneous model can be obtained, as shown in Figure 13a.The results indicate that the curve of the meso-scale softening damage factor (theoretical damage factor) with immersion time is very similar to the evolution law of the crack threshold damage factor.Figure 13b shows that both the tensile and shear crack ratio and the damage degree at 1 day of immersion are greater than those without immersion, but then all three show a gradually decreasing trend.The possible reason is that although both contact cohesion and contact tensile strength attenuate nonlinearly, the absolute values of their attenuation are different.In addition, the location of crack initiation and propagation also has a certain randomness due to the heterogeneity; thus, there is a large difference in the number of cracks when the numerical model is broken.

Failure Characteristics
The heterogeneity of minerals and meso-structures within rocks has a significant impact on the initiation and propagation of cracks, as well as the macroscopic failure patterns.Since the failure modes of the numerical models under different immersion times exhibit a high degree of similarity, the heterogeneous model in its natural state is taken as an example to investigate the influence of heterogeneity on the failure characteristics of soft rock.The crack types and failure modes under uniaxial compression and Brazilian splitting conditions are shown in Figure 14.Generally, the failure of the numerical model is the result of the combined action of tensile cracks (T) and shear cracks (S), both of which can be comprehensively determined by the direction of displacement between block units and the magnitude of displacement on either side of the crack.According to Figure 14a, the failure mode can be observed as "Xtype conjugate shear failure" and the tensile cracks can be observed locally within the model.According to Figure 14b, it can be seen that the model is primarily tensile fracture.It is worth mentioning that under uniaxial compression, tensile and shear cracks can occur at six types of contact: quartz-quartz, quartz-feldspar, quartz-cement, feldspar-feldspar, feldsparcement, and cement-cement, among which the three types related to the cement are the main locations where cracks are generated.In the Brazilian splitting numerical simulation, cracks mainly occur at three types of contacts: quartz-cement, feldspar-cement, and cement-cement.The fundamental reason for the above phenomena is that the mechanical properties of cement and the strength of such contact type are significantly lower than those of quartz and feldspar.Therefore, the direction of crack propagation is closely related to the heterogeneity of material composition and contact types.

Conclusions
This paper established a time-dependent strength degradation grain-based discrete element method model by introducing a meso-scale softening damage factor, systematically exploring the influence of confining pressure and softening effect on the damage and failure characteristics of soft rock.The main conclusions are as follows: (1) The softening process of soft rock is numerically characterized by the attenuation process of contact cohesion and contact tensile strength with d(t).The simulated strengths of the heterogeneous models under different immersion times are in good agreement with the experimental results, thereby verifying the validity of the mesoscale softening damage factor.(2) The simulation accuracy of the heterogeneous models under different confining pressures is above 85%.The crack initiation threshold and damage threshold, the tensile and shear crack ratio, and the damage degree are linearly positively correlated with confining pressure.(3) The attenuation trends of the crack initiation threshold and damage threshold with the immersion time exhibit robust agreement with the meso-scale softening damage factor.During the loading process, the fluctuation degree of the tensile and shear crack curves and the damage degree curve is relatively large.The tensile and shear cracks do not completely tend to be stable near the peak strength.Each damage parameter shows a trend of increasing first and then decreasing with the increase of immersion time.(4) The location of crack initiation and propagation paths is closely related to the heterogeneity of material composition and contact properties.Contact elements related to cement are important locations for crack initiation and propagation.

Figure 1 .
Figure 1.Mesoscopic image of silty mudstone under a polarizing microscope.Figure 1. Mesoscopic image of silty mudstone under a polarizing microscope.

Figure 1 .
Figure 1.Mesoscopic image of silty mudstone under a polarizing microscope.Figure 1. Mesoscopic image of silty mudstone under a polarizing microscope.

Figure 5 .
Figure 5. Simulation results based on calibration parameters: (a) uniaxial compression simulati curve and failure mode; (b) Brazilian fracture simulation curve and failure model; (c) simulati error of the peak axial stresses under different confining pressures; and (d) stress-strain simulati curves of the quartz and feldspar models.

Figure 5 .
Figure 5. Simulation results based on calibration parameters: (a) uniaxial compression simulation curve and failure mode; (b) Brazilian fracture simulation curve and failure model; (c) simulation error of the peak axial stresses under different confining pressures; and (d) stress-strain simulation curves of the quartz and feldspar models.

Figure 6 .Figure 6 .
Figure 6.Simulation results of the heterogeneous model based on calibration parameters: (a) uni ial compression simulation curve and failure mode; (b) Brazilian fracture simulation curve and f ure model; and (c) simulation error of the peak axial stresses under different confining pressures

Figure 8 .
Figure 8. Relationship between damage parameters and confining pressure in heterogeneous models: (a) crack threshold; and (b) tensile and shear crack ratio and damage degree.

Figure 10 .Figure 10 .
Figure 10.Strength simulation results and errors of heterogeneous models with different immersion time: (a) uniaxial compression condition; and (b) Brazilian splitting condition.

Figure 12 .Figure 13 .
Figure 12.Failure modes of heterogeneous models with different immersion time.Figure 12. Failure modes of heterogeneous models with different immersion time.

Figure 13 .
Figure13.The relationship between damage parameters and immersion time in heterogeneous models: (a) crack threshold and its damage factor; and (b) tensile and shear crack ratio and damage degree.

Table 1 .
Basic mechanical parameters of the material components.

Table 2 .
Contact parameters between the same components.

Table 3 .
Simulation results and errors of mechanical parameters of the homogeneous model.

Table 4 .
Finalized meso-scale contact parameters of the six contact types.

Table 5 .
Simulation results and errors of mechanical parameters of the heterogeneous model.