A Three-Level Graded Characterization Framework for the Mesoscopic Mechanical Behavior of Ultra-Deep Heterogeneous Dolomite Based on the Combined Finite-Discrete Element Method

1. Introduction

Ultra-deep dolomite formations with burial depths greater than 6000 m have become important targets for petroleum scientific drilling and production engineering [1]. Such formations are typically characterized by high in situ stress, high pore pressure and strong rock heterogeneity [2,3,4]. Current understanding of rock mechanical behavior at these depths remains limited, and the mechanical parameters and failure descriptors required for drilling and stimulation design are still scarce. This deficit contributes to low rock-breaking efficiency, frequent wellbore instability and insufficient stimulation pressure, which collectively restrict the safe and efficient development of ultra-deep oil and gas resources [5,6,7].

The mechanical evaluation of ultra-deep dolomite is complicated by the multi-source heterogeneity of the host rock. The coexistence of varying mineral composition, multiscale pores, natural fracture networks and complex cementation textures gives rise to pronounced anisotropy in both elastic response and failure behavior [8,9,10]. For the Cambrian subsalt dolomite of the Tarim Basin, Yin et al. [11] demonstrated that integrating mechanical and petrophysical indicators into a reservoir-quality classification framework provides a practical pathway for differentiating reservoir-section behavior and guiding acid-fracturing design. Such field-scale practice illustrates the value of organizing heterogeneous mechanical descriptors into a graded framework, yet the upstream basis at the mesoscopic scale remains underdeveloped.

Numerical simulation has become a powerful tool for investigating the mechanical behavior of heterogeneous rocks [12,13,14]. The combined finite-discrete element method (FDEM), which integrates the advantages of continuum and discontinuum approaches, has been increasingly applied to capture crack initiation, propagation and coalescence in rocks [15,16,17,18]. To incorporate mesoscopic heterogeneity into FDEM, the Voronoi-based grain model captures polygonal mineral particles but suffers from boundary identification difficulty and limited engineering-scale applicability [19]. The Weibull-based random parameter assignment scheme, advanced by Deng et al. [20,21], offers a statistically efficient alternative. Recent FDEM heterogeneity studies, however, have largely addressed individual factors in isolation, such as the sensitivity to fracture energies [22], the influence of heterogeneity on the fracture process zone [23], or the effect of a single mineral phase [24]. A unified framework that systematically organizes mineralogical, structural and scale-related heterogeneity into a transferable grading template, with quantitative behavioral thresholds for each level, is still lacking and constitutes the central motivation of the present study.

The main novelty of this work is the reorganization of mineralogical, structural and scale-related heterogeneity into a single hierarchical framework with explicit, figure-traceable thresholds, rather than treating these factors as parallel, isolated effects. The resulting grading template, transferable to other mineralogically heterogeneous deep rocks, constitutes the principal methodological contribution of the present study.

2. Methodology

2.1. Combined Finite-Discrete Element Method

The combined finite-discrete element method (FDEM), originally proposed by Munjiza [15,25] and implemented in the Y-Geo code with the Y-GUI preprocessor by Mahabadi et al. [16,26], integrates the continuum description of intact material with the discontinuum description of fracturing. Initially thicknessless quadrilateral cohesive elements are inserted between adjacent triangular finite elements (Figure 1), and the transition from continuum to discontinuum is realized through the deformation and failure of these cohesive elements. Three failure modes are distinguished according to the local normal stress and the relative displacement of the crack edges, namely tensile failure (Mode I), shear failure (Mode II) and mixed tensile-shear failure (Mode I-II), each governed by a strain-hardening stage before peak and a strain-softening stage thereafter. At each stage the normal stress transmitted by a cohesive element under the three failure modes is given by Equation (1), and the corresponding shear stress by Equation (2), where the post-peak shear stress in the tension-shear and compression-shear states is treated with a continuous pre-peak to post-peak transition and a residual shear stress. The complete derivation of the cohesive constitutive model follows Munjiza [15], Mahabadi et al. [16] and Lisjak et al. [27]. Each cohesive-element failure is additionally recorded as one acoustic emission (AE) event, providing a discrete record of microcracking activity during loading.

(1)

(2)

2.2. Heterogeneous Properties of Ultra-deep Dolomite

The dolomite studied here was sampled from the Cambrian section at 8715 to 8720 m in the Tarim Basin, where the wellbore rock mass is poorly stable and complex drilling accidents occur frequently. X-ray diffraction on the cores gives a mineral composition of 45.74% dolomite, 37.97% calcite, 14.54% quartz, 1.32% clay minerals, 0.41% feldspar and 0.02% pyrite, so dolomite, calcite and quartz are the three dominant minerals (Figure 2b). Pronounced differences in density and acoustic velocity among the recovered blocks confirm that heterogeneity is an intrinsic characteristic of this rock. Scanning electron microscopy with energy-dispersive spectroscopy (Figure 2c,d) reveals a densely packed yet heterogeneous microstructure. Closely embedded dolomite crystals build a stable stress-transfer skeleton, whereas microcracks, pores and localized secondary mineral filling are unevenly distributed between regions. This dense but heterogeneous fabric provides an efficient load-bearing system and implies a complex mechanical response.

3. Mesoscopic Heterogeneous Dolomite Model

3.1. Model Construction

The mesoscopic model is built in five steps (Figure 3). The dominant mineral composition is first obtained from the mineral-composition experiments, and the elastic modulus and Poisson ratio of the dominant minerals are measured by nanoindentation. A two-dimensional geometric model of the specimen is then constructed in Y-GUI, and the basic properties of the triangular finite elements are randomly assigned from the experimental parameters. An element-property optimization algorithm removes malformed grains so that the mineral distribution approaches the mineral-identification image. To represent heterogeneity within a single mineral, the elastic modulus is assigned following a two-parameter Weibull distribution. The cumulative distribution function used to draw the modulus value of each element is given by Equation (3), in which the expected modulus and the homogeneity index m appear, and a larger m corresponds to a more homogeneous specimen. Two independent indices are used, because heterogeneity arises from two distinct sources. The matrix index m₁ governs the dispersion of matrix-particle properties such as elastic modulus and Poisson ratio, whereas the cementation index m₂ governs the dispersion of the cohesive-element strength parameters such as tensile strength, cohesion and Mode I and Mode II fracture energy. Because cementation between mineral grains is weak and randomly distributed, and deep dolomite exhibits complex physical and mechanical characteristics under in situ conditions [28], this dual-Weibull scheme, following Deng et al. [20,21], captures the resulting mesoscopic heterogeneity.

(3)

3.2. Model Parameters and Macroscopic Response under Uniaxial Compression

The model parameters are summarized in Table 1, with the elastic properties of the three dominant minerals determined experimentally by nanoindentation. The mineral-identification map, the randomly distributed model and the optimized model are compared in Figure 4. The optimized specimen is highly similar to the identification map and satisfies the randomness requirements of mineral spatial distribution, size and shape, qualifying as a numerical specimen of ultra-deep heterogeneous dolomite. Figure 5 compares the stress-strain curves and AE responses of homogeneous and heterogeneous specimens under uniaxial compression. The homogeneous specimen shows a linear pre-peak slope, an essentially constant elastic modulus, a rapid linear post-peak drop and AE concentrated near the peak, which does not reproduce the pronounced pre-peak nonlinearity commonly observed in ultra-deep rocks. The heterogeneous specimen develops pronounced pre-peak nonlinearity, because the higher dispersion of cohesive-element properties promotes local microcracking; its peak strength, peak strain and elastic modulus are lower, while AE initiates at about 30% of the peak load and surges at the post-peak stress drop. Figure 6 shows the damage evolution of the heterogeneous specimen. In the initial stage only a few discretely distributed tensile cracks appear; as loading increases the number and length of tensile cracks grow while shear and tensile-shear cracks emerge; near the peak the tensile-dominated damage forms several band-like clusters; and after the peak the cracks coalesce into steep through-going macroscopic fractures. Macroscopic failure is thus reproduced as a progressive accumulation of non-uniform microcracks, indicating that the heterogeneous model captures the mechanical behavior of ultra-deep dolomite. A quantitative comparison of the simulated and experimental responses is presented separately as the triaxial validation in Section 5.

Building on this model, the following section systematically examines how mineral composition, matrix and cementation heterogeneity, and the representative observation scale govern the macroscopic mechanical behavior, and organizes them into a three-level graded characterization framework.

4. Three-Level Graded Characterization Framework

The sensitivity behavior of the established model is organized below into three hierarchical levels. The compositional level (Section 4.1) sets the strength magnitude, the structural level (Section 4.2) sets the stiffness, nonlinearity and failure mode of the response, and the scale level (Section 4.3) sets the numerical resolution at which the first two are represented. Their integration and underlying mechanisms are presented in Section 4.4.

4.1. Compositional Level

The compositional level constitutes the base of the proposed grading framework. It captures how the volume fractions of the constituent minerals govern the macroscopic strength magnitude of the heterogeneous dolomite specimen. The mineral identification in Section 2.2, together with the model-construction procedure in Section 3.1, provides the foundation for this level.

Four representative compositions with a dolomite content ranging from 40% to 60%, bracketing the measured value, were investigated. Figure 7 presents the vertical stress fields and failure morphologies under uniaxial compression for the four compositions. The load is preferentially carried by the stiffer dolomite and quartz grains, which have higher elastic moduli, while the softer calcite carries a smaller share of the load. Figure 8a further shows that the peak strength increases monotonically as dolomite content rises from 40% to 60%. The compositional control of macroscopic strength magnitude is therefore robust over the investigated range.

A complementary test examined whether mineral spatial layout, at fixed composition, perturbs the macroscopic response. Three independent realizations of the reference composition were generated by varying the random seed of the spatial-assignment algorithm. As shown in Figure 9 and Figure 8b, the three realizations yield nearly coincident stress-strain curves and comparable failure morphologies, with peak strengths stabilizing around 58 MPa. The macroscopic response is therefore statistically insensitive to local spatial layout, indicating that the compositional level can be reliably parameterized by the bulk mineral content alone.

In summary, the compositional level reflects the volume-fraction control of the load-bearing skeleton. Dolomite content acts as the primary compositional indicator, and the macroscopic strength magnitude scales with this indicator robustly against random spatial perturbations. The mechanistic basis is elaborated in Section 4.4.3.

4.2. Structural Level

The structural level builds upon the compositional base and captures how the statistical distribution of mesoscopic mechanical properties, at a given composition, governs the dispersion and nonlinearity of the macroscopic response. Following the Weibull-based discretization in Section 3.1, two independent shape parameters are assigned: m₁ controls the distribution of matrix-particle properties such as elastic modulus and Poisson ratio, and m₂ controls the distribution of cohesive-element strength parameters such as tensile strength, cohesion and Mode I and Mode II fracture energy.

A two-dimensional sensitivity scan was performed with m₁ and m₂ each taking six values, namely 1.5, 2.0, 3.0, 4.0, 8.0 and 12.0, generating a six-by-six matrix of simulations. Figure 10 shows the stress-strain curves and AE responses for fixed m₁ with varying m₂: as m₂ decreases, the pre-peak nonlinearity becomes more pronounced, the post-peak stress drop more gradual, and the AE onset shifts earlier. Figure 11 shows the complementary scan with fixed m₂ and varying m₁: m₁ exerts a modest influence on peak strength and strain but a comparatively weaker influence on the overall curve shape.

Figure 12 summarizes the evolution of the three macroscopic parameters across the full m₁- m₂ matrix, and three observations are central. First, the influence of each parameter saturates: once m₁ or m₂ exceeds about 4, further increase produces only marginal change in peak strength, peak strain and elastic modulus, so the macroscopic response becomes weakly dependent on the shape parameter beyond this value, with the saturation consistently located near m = 4. Second, the elastic modulus is governed almost entirely by m₁ and is essentially insensitive to m₂, since the modulus curves obtained at different m₂ collapse onto a single m₁-controlled trend (Figure 12e,f); matrix-particle heterogeneity therefore controls the overall stiffness, whereas the pre-peak nonlinearity, the AE onset, the failure mode and the magnitude and dispersion of peak strength and peak strain are governed primarily by m₂ (Figure 12a,c), so that cementation-strength heterogeneity controls microcrack initiation, the load-bearing magnitude and the local instability of the skeleton. Third, m₁ and m₂ act jointly rather than independently: the influence of m₁ on peak strength and peak strain is conditional, remaining weak at small m₂ and becoming appreciable only once m₂ is large (Figure 12b,d); when the cementation is strongly heterogeneous the response is governed by the weakest cohesive elements and is largely insensitive to the matrix, whereas for a more uniform cementation the matrix heterogeneity controlled by m₁ begins to govern the peak strength and strain.

The saturation threshold m = 4 is therefore adopted as the behavioral threshold of the structural level. Its statistical interpretation, rooted in the skewness transition of the Weibull distribution near m = 3.6, is presented in Section 4.4.3. The differentiated yet jointly acting control of m₁ and m₂ reflects the complex internal heterogeneity of ultra-deep dolomite and is treated as an intrinsic feature of the structural level rather than a separately parameterized interaction term.

4.3. Scale Level

The scale level constitutes the representational prerequisite of the framework. It captures how the chosen mesh size influences the model capacity to resolve the compositional and structural heterogeneity established at the first two levels. Unlike the previous two levels, which describe intrinsic rock properties, the scale level concerns the numerical resolution at which those properties are sampled.

A series of uniaxial-compression simulations was performed with mesh sizes of 0.5, 0.6, 0.8, 1.0, 1.2, 1.4, 1.5, 1.6, 1.8 and 2.0 mm, all other parameters identical to those of the reference case in Table 1. Figure 13 compares the microcrack distributions and failure morphologies at the ten mesh sizes. For a mesh size not exceeding 1.2 mm, a well-defined penetrating main failure surface develops with concentrated microcracks. As the mesh size increases beyond 1.2 mm, the failure morphology becomes progressively dispersed, the main failure surface ceases to be clearly identifiable, and microcracks distribute over a larger volume in a weakly localized manner; the coarser meshes thus underestimate stress concentrations at high-to-low-modulus interfaces and fail to reproduce the dominant failure mechanism.

Figure 14 confirms this quantitatively: both peak strength and elastic modulus stabilize for a mesh size not exceeding 1.2 mm, with relative variations below 5%, while further refinement yields only marginal change at sharply rising computational cost. The behavioral threshold of the scale level is therefore set at a mesh size not exceeding 1.2 mm. This threshold ensures the mesh resolves the mineral-domain topology at the mesoscopic scale, and acts as the prerequisite for the compositional and structural levels to be reliably captured: any sensitivity analysis at coarser resolution would be confounded by resolution-related artifacts.

4.4. Integrated Grading Framework and Mechanisms

4.4.1. Hierarchical Relationship among the Three Levels

The three levels established in Section 4.1 to 4.3 are not parallel sensitivity factors but possess an intrinsic hierarchy, which is the central methodological contribution of the proposed framework (Figure 15). The compositional level forms the base of the hierarchy: mineral composition is the material substrate and determines the magnitude of the load-bearing capacity; without a defined composition, neither the heterogeneity of mesoscopic properties nor the resolution of their representation can be meaningfully discussed. The structural level acts as a modifier on the compositional base: at a fixed composition and fixed mean properties, the Weibull-distributed heterogeneity of matrix and cementation governs the elastic modulus, the pre-peak nonlinearity, the failure mode and the dispersion of the macroscopic response, and, through the weak-link behavior of the more heterogeneous cohesive elements, also modulates the peak strength and strain about the level set by the compositional baseline. The scale level serves as the representational prerequisite for the other two: the mesh size governs whether the lower levels are resolved at sufficient numerical fidelity. This three-tier architecture, comprising a material base, a statistical modifier and a representational prerequisite, distinguishes the present framework from the conventional treatment of mesoscopic heterogeneity as a set of parallel sensitivity factors.

4.4.2. Quantitative Grading Indicators and Thresholds

The quantitative grading indicators and behavioral thresholds derived from the three sensitivity analyses are consolidated in Table 2. Each row lists one controlling factor, its investigated range or behavioral threshold, the macroscopic indicator it governs, the underlying physical basis, and the source figures. Two features are worth emphasizing. First, all behavioral thresholds are directly traceable to the source figures: the compositional range of 40% to 60% dolomite reflects the actual scan in Figure 7 and Figure 8a; the saturation threshold m = 4 corresponds to the converged region in Figure 12; and the convergence threshold of 1.2 mm is read directly from Figure 14. No threshold is extrapolated beyond the simulated range. Second, the physical-basis column anchors each threshold to a material property, a statistical structure or a numerical-resolution requirement, so the thresholds acquire mechanistic meaning beyond their numerical values, as elaborated in Section 4.4.3. The Weibull-skewness mechanism underlying the structural-level threshold is plotted in Figure 16, which marks the symmetry transition at m = 3.6 and locates the simulation-based threshold m = 4 within the statistically symmetric regime. Through the joint use of Table 2, Figure 15 and Figure 16, the framework offers a transparent and reproducible grading template for the mesoscopic mechanical evaluation of ultra-deep heterogeneous dolomite.

4.4.3. Mechanistic Interpretation of the Grading Thresholds

The hierarchical structure reflects three independent physical mechanisms. The compositional level expresses the volume-fraction control of the load-bearing skeleton: as listed in Table 1, dolomite and quartz possess elastic moduli of 96.7 GPa and 99.3 GPa, whereas calcite is markedly softer at 69 GPa, so an increase in the combined fraction of the two stiffer phases enlarges the specimen volume that carries the principal compressive load, consistent with the monotonic strength increase in Figure 8a; the vertical stress fields in Figure 7 further show load concentrating within the connected network of high-modulus grains.

The structural-level threshold m = 4 has a clear statistical interpretation: for a two-parameter Weibull distribution the skewness transitions from positive, that is right-skewed and dominated by sub-mean values, to nearly zero as m approaches 3.6. Below this transition, a substantial fraction of cohesive elements carry sub-mean strength and act as preferential nucleation sites for microcracks, producing the pronounced pre-peak nonlinearity and early AE onset in Figure 10, consistent with the microcrack-controlled stress-strain behavior of brittle materials [29]. Once m exceeds about 4, the distribution is sufficiently symmetric that the cohesive-element population is statistically homogeneous at the mesoscale, microcrack initiation approaches spatial uniformity, and further increase of m yields diminishing influence. The simulation-based threshold m = 4 in Figure 12 is therefore consistent with the underlying statistical structure of the Weibull-distributed properties rather than an empirical artifact.

The scale-level convergence threshold of 1.2 mm is governed by the requirement that the mesh resolve the compositional and structural heterogeneity at a representative mesoscale; the chosen size is small relative to the characteristic mineral-grain dimension, providing enough elements to discretize the mineral-domain topology. Coarser meshes blur the mineral-boundary topology and underestimate local stress concentrations, as manifested by the dispersed, weakly localized failure for mesh sizes above 1.2 mm in Figure 13; the convergence of peak strength and modulus in Figure 14 corroborates the threshold, in line with the principle that the representative element size should scale with, and not below, the characteristic heterogeneity length of the material.

4.4.4. Generalizability of the Framework

The three-tier architecture, comprising a compositional base, a structural modifier and a scale prerequisite, is designed to be transferable beyond the present case. Although established for Cambrian ultra-deep dolomite, the architecture itself is expected to remain valid for a broad class of mineralogically heterogeneous deep rocks, including other carbonates, tight sandstones and igneous rocks. Application to a new lithology requires only re-determination of the quantitative thresholds of each level against the corresponding mineral composition, heterogeneity statistics and grain scale, whereas the hierarchical logic, namely composition setting the baseline strength magnitude, structure governing the stiffness, nonlinearity and failure mode while modulating the strength and strain, and scale setting the resolution, remains unchanged.

5. Validation under Triaxial Loading

This section validates the mesoscopic model under triaxial loading and verifies that the three-level grading established under uniaxial compression remains consistent in the high-confining-pressure regime relevant to ultra-deep formations. Numerical simulations and laboratory experiments were conducted at six confining pressures, namely 0, 60, 100, 140, 180 and 220 MPa. The model parameters established in Section 3 (Table 1) were adopted, with other inputs set per Tatone and Grasselli [30]. Figure 17 compares the stress-strain curves from simulation and experiment across the six confining pressures, showing good agreement over the full range. Figure 18 compares the elastic modulus, peak strength and residual strength versus confining pressure: all three increase nonlinearly with confining pressure, with the strengthening most pronounced between 0 and 100 MPa and tending to saturate above about 100 MPa. The model established under uniaxial compression therefore reproduces the experimental response across the full 220 MPa range, indicating that the three-level grading remains applicable under varying confining pressure. The three-level grading remains internally consistent across the triaxial range: the compositional-level ordering of peak strength is preserved at all pressures; the structural-level saturation threshold m = 4 continues to delimit the regime within which heterogeneity governs pre-peak nonlinearity and failure mode; and the scale-level convergence threshold, set by the mineral-domain topology rather than the loading condition, is preserved in the triaxial simulations. The framework is thus applicable across the loading conditions relevant to ultra-deep dolomite formations.

6. Limitations

This study is a two-dimensional mesoscopic modeling investigation, and several limitations should be acknowledged; each indicates a direction for future work.

(1) Two-dimensional modeling. The framework is established under plane two-dimensional conditions, which resolve the in-plane mineral-domain topology but do not represent out-of-plane microstructural heterogeneity. Three-dimensional extension may shift the scale-level convergence threshold and modify the relative roles of the structural-level parameters.

(2) Compositional coverage. The compositional level was characterized within a dolomite-content range of 40% to 60%, bracketing the measured composition. Whether the monotonic strength dependence holds beyond this range, and whether additional sub-levels are required for clay-rich, silica-rich or evaporitic end-members, requires a broader sample library.

(3) Dry, room-temperature conditions. All simulations and the corresponding experiments were performed under dry conditions at room temperature. The weakening effect of pore fluid, the temperature-dependent transition toward ductility, and chemo-mechanical interactions in deep-formation environments are therefore not represented, so the present results may underestimate these effects in situ.

(4) Resolution of the parameter scan. The structural-level thresholds were obtained from a six-point sensitivity scan of the two Weibull shape parameters, and the model was established for a single dolomite lithology. Finer scans near the m = 4 transition and extension to additional lithologies would further constrain the quantitative thresholds, although the underlying Weibull-skewness mechanism is expected to remain valid.

7. Conclusions

Based on the combined finite-discrete element method, a two-dimensional mesoscopic numerical model for Cambrian ultra-deep dolomite from the Tarim Basin was developed. By incorporating X-ray diffraction and scanning electron microscopy with energy-dispersive spectroscopy characterization and a Weibull-distribution discretization scheme for the heterogeneity of both matrix particles and inter-particle cementation, a three-level graded characterization framework was proposed and validated. The main conclusions are as follows.

(1) A dual Weibull discretization scheme for matrix and cementation heterogeneity was established. The two independent shape parameters m₁ and m₂ exert differentiated and coupled control on the macroscopic response: m₁ governs the elastic modulus, essentially independently of m₂, whereas m₂ governs the pre-peak nonlinearity, the AE onset, the failure mode and the magnitude and dispersion of peak strength and peak strain; the influence of m₁ on peak strength and strain is conditional, becoming appreciable only when m₂ is large. This differentiated yet coupled control reflects the complex internal heterogeneity of ultra-deep dolomite, which cannot be captured by a single-parameter discretization.

(2) A three-level graded characterization framework was proposed, comprising a compositional level that sets the baseline strength magnitude, a structural level that governs the elastic modulus, pre-peak nonlinearity and failure mode, and a scale level that controls the representativeness of the mesoscopic response. Quantitative grading indicators and behavioral thresholds were derived directly from the sensitivity analyses, namely the monotonic strength dependence on dolomite content within 40% to 60%, the saturation of structural-level influence at m = 4, and convergence at a mesh size not exceeding 1.2 mm. The structural-level threshold is supported by the statistical skewness transition of the two-parameter Weibull distribution.

(3) The framework was further validated under triaxial loading. Simulations and experiments at 0 to 220 MPa show good agreement in stress-strain response, with the elastic modulus, peak strength and residual strength all increasing nonlinearly and saturating above about 100 MPa. The three-level grading remains internally consistent throughout the loading range, confirming the applicability of the framework to ultra-deep formation conditions.

The proposed framework provides a transferable methodological reference for the mesoscopic mechanical evaluation of ultra-deep heterogeneous dolomite and similar mineralogically heterogeneous deep rocks. Re-determination of the thresholds for new lithologies and extension to three-dimensional and multi-physics conditions are envisaged as natural directions for future work.