# Structural Damage Evolution of Mesoscale Representative Elementary Areas of Mudstones

Authors: Qijun Hu; Junsen Zeng; Leping He; Chengyi Huang; Shiqing Feng
Journal: Geofluids
 (2023)

---

## Abstract

The mechanical behavior of mudstones is significantly affected by their mesostructure. This study investigates the damage evolution of red mudstone mesostructures under uniaxial compression through U-Net image segmentation, mesoscopic representative elementary area (mREA), meso-element equivalent method, and linearly superimposed model. Results show that the mREAs can be simplified as a binary structure comprising nonclayey minerals and a porous matrix. The inclusion-matrix interfaces constituted by the dominant inclusions, such as large inclusions and a series of inclusions with specific arrangements, are the weakest regions at the postpeak stage. The shear failure of these interfaces leads to general shear failure with a penetrating principal fracture or local shear with several nonpenetrating local fractures in the mREAs. The damage evolution and failure mode are sensitive to the inclusion geometric features, e.g., inclination, grain size, and spacing. However, the uniaxial compressive strength (UCS) is mainly determined by the included angle between the dominant inclusion inclination and the loading direction. With consistent composition, an increase in the included angle from 0 to 90° results in an increase in the UCS of the muddy facies and a decrease in the sandy facies, with the UCS of the sandy facies being more sensitive to changes in the included angle. Our findings indicate that the structural failure of red mudstones from mesoscopic to macroscopic is driven by the accumulation and expansion of shear failure at multiple scales due to the presence of hard–soft interfaces. This study enhances our understanding of the structural failure mechanism of mudstones and serves as a practical reference for the mesomechanical testing of clay rocks.


---

## Body

## 1. Introduction

Mudstones are fine-grained sedimentary rocks that are widely investigated in nuclear waste disposal, CO2 geological sequestration, and geotechnical engineering [1–3]. Mudstone includes tectosilicate minerals (e.g., feldspar and quartz), clay minerals (e.g., montmorillonite, kaolinite, and illite), and other minerals. The mudstone mechanical behavior is significantly affected by its meso-/microstructure at the submillimeter scale [4]. The mudstone mesostructures can be defined as a matrix mixed with inclusions [5]. For a stiff-based structure, the matrix is the nonclayey mineral grains, and the inclusions are the porous substance, which mainly consists of small nonclayey mineral grains, layered clay minerals, and pores. In contrast, the nonclayey mineral grains become the inclusions when their content is significantly lower than that of the porous substance. The stress–strain relationship of mudstones under uniaxial or triaxial compression can be divided into the elastic stage and the nonlinear strengthening, softening, and residual stages [6]. It is well accepted that damage is responsible for these nonlinear phenomena, but little is known about the specific damage mechanisms of mudstones [4].Historically, the phenomenological plastic damage models of rocks were derived in the framework of thermodynamics, continuum damage mechanics, and geomaterial constitutive models [7–9]. Recently, Shen et al. and Zhao et al. derived the multiscale plastic incremental model of porous geomaterials [10–12]. For the experimental approaches, the backscattered electron (BSE) mode of scanning electron microscopy (SEM) is commonly employed to observe the mineral composition and spatial distribution of mesostructures [13]. However, the corresponding BSE images are only tens of nanometers to tens of micrometers in size, which cannot illustrate the mesostructural heterogeneity and the representativeness of the region of interest (ROI). Therefore, a representative elementary area (REA) with a characteristic length LREA was proposed. It was determined by estimating the convergence of constituent geometric parameters in a domain that gradually expanded in the captured digital structure [5, 14–17]. The LREA of clay rocks was approximately hundreds of microns at an image resolution of dozens of nanometers [18–20]. Regarding the mesoscopic mechanical tests, Wang et al. obtained three types of deformation fields of mudstone microstructures under uniaxial compression by in situ loading tests withenvironmental scanning electron microscopy (ESEM) and digital image correlation (DIC) technology [4]. Sui et al. studied the microstructural deformation of a tight sandstone by the in situ loading/scanning experiment and deformation simulation method [21]. Shi et al. investigated the three-dimensional nonuniform strain fields and cracking patterns of claystone by the uniaxial compression creep test, digital volume correlation (DVC) method, and X-ray microtomography [22, 23].The use of numerical simulation has enabled investigating the mechanical responses inside rock mesostructures under loads. Discrete element methods (DEMs) have advanced in revealing the crack propagation and failure mode of mudstones [24–26]. Regrettably, most of the previous studies concentrated on macroscopic specimens in which the element size was much larger than the grain size of minerals. In addition, the material parameters were difficult to calibrate, even with mechanical tests and DIC technology [27]. Finite element methods (FEMs) are popular in the mesomechanical investigation of composite materials [28–30]. For example, to realize structural and material heterogeneity, Du et al. proposed the meso-element equivalent method (MEEM), which defined the unique material properties for each element based on the content ratio, Voigt parallel connection model, and trilinear model [31, 32]. The mesomechanical properties of mudstones are a combination of the mineral grains, interfaces, and pores/cracks, which are difficult to measure, even with nanomechanical instruments [33, 34]. Recently, Li et al. studied the mesodamage of shale REA by the concrete damaged plasticity (CDP) model [35], in which the mesomechanical properties were calibrated by macroscopic mechanical tests. However, their work primarily examined the pore influences but did not account for the potential interaction between inclusions and the matrix. In summary, although these studies investigated the relationships between the mudstone meso-/microstructure and mechanical properties, the representativeness of the region of interest and the influences of the mesostructure on the damage evolution remain unclear.The main objective of this study is to characterize the damage evolution of the mesoscopic representative elementary area (mREA) of mudstones. This paper is organized as follows: first, the mREAs of four red mudstone mesostructures from the Sichuan Basin in China are obtained by the proposed mREA estimated procedure [20]. Next, a series of uniaxial compression testing are simulated with the MEEM and CDP models, including the natural mREAs and artificially designed mesostructures. The results are presented and discussed in terms of the probability density distribution (PDD) of elastic parameters, stress–strain relations, strain fields, damage evolution, and failure mode. The emphasis is placed on the analysis of mesostructural damage characteristics and mechanism in relation to the composition, geometric features, and spatial distribution of constituents.


## 2. Materials and Methods

### 2.1. Material and Sample Preparation

The tested materials were four moderately weathered red mudstone samples cored at foundation pits in Chengdu (S1-S3, muddy-facies, and Guankou formation) and Chongqing (S4, sandy-facies, and Shaximiao formation), China. The mean uniaxial compressive strengths (UCS) of the core specimens (φ50×100 mm3) of samples S2-S4 were 2.73, 7.29, and 7.96 MPa. The core specimen of sample S1 was unsuccessfully prepared due to its weak water resistance. Table 1 lists the major chemical composition of the four red mudstone samples measured by an X-ray fluorescence analyzer (X’Pert PRO, PANalytical B.V.).Table 1
Major chemical composition of the red mudstone samples.
Oxides (wt%)Sample S1Sample S2Sample S3Sample S4SiO252.5754.6166.1163.68Al2O318.6019.5418.6418.58Fe2O38.867.166.305.93K2O4.194.283.313.79MgO3.353.452.182.98CaO10.438.591.062.37TiO20.890.910.980.87Na2O0.670.910.821.26P2O50.160.200.300.25To prepare the specimens for BSE imaging, the four samples were cut into cubes of approximately2×2×0.8 cm3. Then, the observing surface was polished with 800-10000 graded SiC abrasive papers. Finally, argon-ion polishing (Leica EM ACE200 and TIC3X) was used to polish the observing surface for 4-5 h. The BSE images were obtained by FESEM (FEI Quanta 650 FEG) with a 30 kV voltage and 9-12 mm working distance. A 5000x magnification with an approximately 54 nm·pixel-1 resolution allowed clear identification of the mesostructures of samples S1-S3, whereas 1300x magnification with an approximately 208 nm·pixel-1 resolution was used for sample S4. The numbers of images captured of the four samples were 100 (10×10), 30 (10×3), 100 (10×10), and 25 (5×5). After cropping the overlapping areas and edge areas, the 100 stitched images provided an image domain of approximately 691×320 μm2, as opposed to approximately 1331×832μm2 for the 25 stitched images of sample S4.


### 2.2. Estimation of Mesoscopic Representative Elementary Area

A three-step workflow was used to segment the BSE images and estimate the mREA [20]. First, the U-Net was employed to segment the mesostructure. According to the mineralogical characteristics and their mechanical properties, the mesostructure was classified into metallic minerals, nonclayey mineral aggregations (hereinafter referred to as the nonclayey minerals), the matrix comprising small nonclayey mineral grains and clay minerals, the transition parts between the nonclayey minerals and matrix, and pores/cracks, as shown in Figure 1. Next, based on the segmented stitched image, the characteristic length (LREA, i.e., the side length of the mREA) was estimated by calculating the geometric parameter convergence as the domain size increased, including the area fraction, mean area, and the pore inclination and circularity. Eventually, based on the mean structural similarity (MSSIM), the mREAs were determined from the square domains that had the side length of the LREA in the stitched image. The most representative region (MRR) had the highest score of the MSSIM, which could be considered the “average face” in these square domains. In contrast, the most special region (MSR) had the lowest score of the MSSIM.Figure 1
Schematic diagram of the image segmentation and the meso-element equivalent method.


### 2.3. Meso-Element Equivalent Method

In the original paper on the MEEM, the mesostructure was assumed as a matrix mixed with randomly distributed spheres [31, 32]. Then, the mesostructure was segmented into domains that had a consistent size with the elements. In this study, the segmented mREA was further simplified to a binary structure, as shown in Figure 1. The captured pores were ignored because their number and size were too small compared to the LREA and other constituents. The elements were defined with different material properties according to the inclusion area ratio. The explicit method with a mass scale factor (104) and loading time (0.01) was used when the implicit calculation was aborted due to the excessive nonlinearity of the material properties between elements. Although the UCS and residual strength under this optimization was still approximately 20% lower than that of the implicit method, their damage origination and evolution were consistent.


### 2.4. Linearly Superimposed Model

The linearly superimposed model, based on the Voigt parallel model, was employed to describe the stress–strain relationship of elements. Since theLREA was much smaller than the size of the macroscopic core specimen, the mesomechanical properties of the two basic elements, i.e., the stiff and soft elements, were assumed to meet the macromechanical properties of two typical red mudstones. It was divided into the sandy facies with a high content and large grain size of nonclayey minerals (mainly from the Badong formation) and the muddy facies with a relatively high content of the matrix (mainly from the Guankou formation). The data were collected from the uniaxial compression test and the triaxial compression tests conducted with low confining pressures, as listed in Table 2 [36–40].Table 2
Compressive stress–strain parameters of the stiff and soft elements.
ParametersStiff elementSoft elementPoisson’s ratio0.100.30Peak stressfc/MPa248Elastic stress ratioλ10.500.50Residual stress ratioλ20.500.50Elastic strainε0/%0.400.50Peak strainε1/%0.701.20Residual strainε2/%0.851.60Failure strainεcu/%2.002.00Based on the area ratio and Voigt parallel model, the effective elastic modulusE∗and Poisson’s ratio v∗ of an equivalent element were given as
(1)E∗=CagEag+CmaEma,v∗=E∗2G∗−1=CagEag+CmaEmaCagEag/1+vag+CmaEma/1+vma−1,where ag represents the nonclayey minerals, ma represents the matrix, G∗ is the effective shear modulus, E is the elastic modulus, C is the area fraction, and v is the Poisson’s ratio.The compressive stress–strain relationship for each equivalent element was defined as(2)σ=12Cagfcagε0ag+Cmafcmaε0maε,ε≤ε0ag,12Cagfcag1+ε−ε0agε1ag−ε0ag+12Cmafcmaεε0ma,ε0ag<ε≤ε0ma,12Cagfcag1+ε−ε0agε1ag−ε0ag+12Cmafcma1+ε−ε0maε1ma−ε0ma,ε0ma<ε≤ε1agCagfcag1−ε−ε1ag2ε2ag−ε1ag+12Cmafcma1+ε−ε0maε1ma−ε0ma,ε1ag<ε≤ε2ag12Cagfcag+12Cmafcma1+ε−ε0maε1ma−ε0ma,ε2ag<ε≤ε1ma,12Cagfcag+Cmafcma1−ε−ε1ma2ε2ma−ε1ma,ε1ma<ε≤ε2ma,12Cagfcag+12Cmafcma,ε2ma<ε≤εcu,,where.fc is the peak stress, ε0 is the elastic limit strain, ε1 is the peak strain, ε2 is the softening strain, and εcu is the failure strain and was assumed to be 0.02.This representation, although somewhat simplified, captured the main features of stress–strain relationship of the red mudstones under uniaxial compression. The compressive and tensile stress–strain relationship was implemented using the concrete damaged plasticity (CDP) model. To meet the calculation and convergence of the CDP model, the compressive stress–strain relationship was further simplified as(3)E0eq=σeqε=ε1agε1ag,σε=ε2ag=ε1ma−ε2agε1ma−ε1agσε=ε1ag+ε2ag−ε1agε1ma−ε1agσε=ε1ma,σε=ε1ag≤σε=ε1ma,σε=ε1ma=ε2ma−ε1maε2ma−ε2agσε=ε2ag+ε1ma−ε2agε2ma−ε2agσε=ε2ma,σε=ε1ag>σε=ε1ma,σε=εcu=0.412Cagfcag+12Cmafcma.There is a relatively high dispersion in the tensile parameters of the clay rocks (including the red mudstone) according to their Brazilian testing results [35, 41, 42]. In addition, according to the preparameter test, the tensile stiffening and postfailure stress–strain relationships had little influence on the mREA mechanical behaviors under the uniaxial compression. Thus, after comparing the relationships between the compressive and tensile parameters of these rocks, the CDP tensile parameters were defined as
(4)σt1=0.1σc1,εt1=0.006,σt2=0.1σt1,εt1=0.007,σtcu=σt2,εcu=0.02,whereσt1 is the peak tensile strength, σc1 is the UCS, σt2 is the residual tensile strength, σtcu is the failure tensile strength, εt1 is the peak tensile strain, εt2 is the residual tensile strain, and εtcu is the failure tensile strain.The Sidoroff method was introduced to describe the compressive damagedt and tensile damage dc during the uniaxial compression [43]:
(5)d=1−σεE∗.


## 2.1. Material and Sample Preparation

The tested materials were four moderately weathered red mudstone samples cored at foundation pits in Chengdu (S1-S3, muddy-facies, and Guankou formation) and Chongqing (S4, sandy-facies, and Shaximiao formation), China. The mean uniaxial compressive strengths (UCS) of the core specimens (φ50×100 mm3) of samples S2-S4 were 2.73, 7.29, and 7.96 MPa. The core specimen of sample S1 was unsuccessfully prepared due to its weak water resistance. Table 1 lists the major chemical composition of the four red mudstone samples measured by an X-ray fluorescence analyzer (X’Pert PRO, PANalytical B.V.).Table 1
Major chemical composition of the red mudstone samples.
Oxides (wt%)Sample S1Sample S2Sample S3Sample S4SiO252.5754.6166.1163.68Al2O318.6019.5418.6418.58Fe2O38.867.166.305.93K2O4.194.283.313.79MgO3.353.452.182.98CaO10.438.591.062.37TiO20.890.910.980.87Na2O0.670.910.821.26P2O50.160.200.300.25To prepare the specimens for BSE imaging, the four samples were cut into cubes of approximately2×2×0.8 cm3. Then, the observing surface was polished with 800-10000 graded SiC abrasive papers. Finally, argon-ion polishing (Leica EM ACE200 and TIC3X) was used to polish the observing surface for 4-5 h. The BSE images were obtained by FESEM (FEI Quanta 650 FEG) with a 30 kV voltage and 9-12 mm working distance. A 5000x magnification with an approximately 54 nm·pixel-1 resolution allowed clear identification of the mesostructures of samples S1-S3, whereas 1300x magnification with an approximately 208 nm·pixel-1 resolution was used for sample S4. The numbers of images captured of the four samples were 100 (10×10), 30 (10×3), 100 (10×10), and 25 (5×5). After cropping the overlapping areas and edge areas, the 100 stitched images provided an image domain of approximately 691×320 μm2, as opposed to approximately 1331×832μm2 for the 25 stitched images of sample S4.


## 2.2. Estimation of Mesoscopic Representative Elementary Area

A three-step workflow was used to segment the BSE images and estimate the mREA [20]. First, the U-Net was employed to segment the mesostructure. According to the mineralogical characteristics and their mechanical properties, the mesostructure was classified into metallic minerals, nonclayey mineral aggregations (hereinafter referred to as the nonclayey minerals), the matrix comprising small nonclayey mineral grains and clay minerals, the transition parts between the nonclayey minerals and matrix, and pores/cracks, as shown in Figure 1. Next, based on the segmented stitched image, the characteristic length (LREA, i.e., the side length of the mREA) was estimated by calculating the geometric parameter convergence as the domain size increased, including the area fraction, mean area, and the pore inclination and circularity. Eventually, based on the mean structural similarity (MSSIM), the mREAs were determined from the square domains that had the side length of the LREA in the stitched image. The most representative region (MRR) had the highest score of the MSSIM, which could be considered the “average face” in these square domains. In contrast, the most special region (MSR) had the lowest score of the MSSIM.Figure 1
Schematic diagram of the image segmentation and the meso-element equivalent method.


## 2.3. Meso-Element Equivalent Method

In the original paper on the MEEM, the mesostructure was assumed as a matrix mixed with randomly distributed spheres [31, 32]. Then, the mesostructure was segmented into domains that had a consistent size with the elements. In this study, the segmented mREA was further simplified to a binary structure, as shown in Figure 1. The captured pores were ignored because their number and size were too small compared to the LREA and other constituents. The elements were defined with different material properties according to the inclusion area ratio. The explicit method with a mass scale factor (104) and loading time (0.01) was used when the implicit calculation was aborted due to the excessive nonlinearity of the material properties between elements. Although the UCS and residual strength under this optimization was still approximately 20% lower than that of the implicit method, their damage origination and evolution were consistent.


## 2.4. Linearly Superimposed Model

The linearly superimposed model, based on the Voigt parallel model, was employed to describe the stress–strain relationship of elements. Since theLREA was much smaller than the size of the macroscopic core specimen, the mesomechanical properties of the two basic elements, i.e., the stiff and soft elements, were assumed to meet the macromechanical properties of two typical red mudstones. It was divided into the sandy facies with a high content and large grain size of nonclayey minerals (mainly from the Badong formation) and the muddy facies with a relatively high content of the matrix (mainly from the Guankou formation). The data were collected from the uniaxial compression test and the triaxial compression tests conducted with low confining pressures, as listed in Table 2 [36–40].Table 2
Compressive stress–strain parameters of the stiff and soft elements.
ParametersStiff elementSoft elementPoisson’s ratio0.100.30Peak stressfc/MPa248Elastic stress ratioλ10.500.50Residual stress ratioλ20.500.50Elastic strainε0/%0.400.50Peak strainε1/%0.701.20Residual strainε2/%0.851.60Failure strainεcu/%2.002.00Based on the area ratio and Voigt parallel model, the effective elastic modulusE∗and Poisson’s ratio v∗ of an equivalent element were given as
(1)E∗=CagEag+CmaEma,v∗=E∗2G∗−1=CagEag+CmaEmaCagEag/1+vag+CmaEma/1+vma−1,where ag represents the nonclayey minerals, ma represents the matrix, G∗ is the effective shear modulus, E is the elastic modulus, C is the area fraction, and v is the Poisson’s ratio.The compressive stress–strain relationship for each equivalent element was defined as(2)σ=12Cagfcagε0ag+Cmafcmaε0maε,ε≤ε0ag,12Cagfcag1+ε−ε0agε1ag−ε0ag+12Cmafcmaεε0ma,ε0ag<ε≤ε0ma,12Cagfcag1+ε−ε0agε1ag−ε0ag+12Cmafcma1+ε−ε0maε1ma−ε0ma,ε0ma<ε≤ε1agCagfcag1−ε−ε1ag2ε2ag−ε1ag+12Cmafcma1+ε−ε0maε1ma−ε0ma,ε1ag<ε≤ε2ag12Cagfcag+12Cmafcma1+ε−ε0maε1ma−ε0ma,ε2ag<ε≤ε1ma,12Cagfcag+Cmafcma1−ε−ε1ma2ε2ma−ε1ma,ε1ma<ε≤ε2ma,12Cagfcag+12Cmafcma,ε2ma<ε≤εcu,,where.fc is the peak stress, ε0 is the elastic limit strain, ε1 is the peak strain, ε2 is the softening strain, and εcu is the failure strain and was assumed to be 0.02.This representation, although somewhat simplified, captured the main features of stress–strain relationship of the red mudstones under uniaxial compression. The compressive and tensile stress–strain relationship was implemented using the concrete damaged plasticity (CDP) model. To meet the calculation and convergence of the CDP model, the compressive stress–strain relationship was further simplified as(3)E0eq=σeqε=ε1agε1ag,σε=ε2ag=ε1ma−ε2agε1ma−ε1agσε=ε1ag+ε2ag−ε1agε1ma−ε1agσε=ε1ma,σε=ε1ag≤σε=ε1ma,σε=ε1ma=ε2ma−ε1maε2ma−ε2agσε=ε2ag+ε1ma−ε2agε2ma−ε2agσε=ε2ma,σε=ε1ag>σε=ε1ma,σε=εcu=0.412Cagfcag+12Cmafcma.There is a relatively high dispersion in the tensile parameters of the clay rocks (including the red mudstone) according to their Brazilian testing results [35, 41, 42]. In addition, according to the preparameter test, the tensile stiffening and postfailure stress–strain relationships had little influence on the mREA mechanical behaviors under the uniaxial compression. Thus, after comparing the relationships between the compressive and tensile parameters of these rocks, the CDP tensile parameters were defined as
(4)σt1=0.1σc1,εt1=0.006,σt2=0.1σt1,εt1=0.007,σtcu=σt2,εcu=0.02,whereσt1 is the peak tensile strength, σc1 is the UCS, σt2 is the residual tensile strength, σtcu is the failure tensile strength, εt1 is the peak tensile strain, εt2 is the residual tensile strain, and εtcu is the failure tensile strain.The Sidoroff method was introduced to describe the compressive damagedt and tensile damage dc during the uniaxial compression [43]:
(5)d=1−σεE∗.


## 3. Results

### 3.1. Mesoscopic Representative Elementary Area

Figure2 shows the stitched BSE image, MRR, MSR, and simplified structures of four samples. According to the mREA estimation, the LREA was approximately 194, 108, 280, and 582 μm for four samples, respectively. Broadly, the MRR described a relatively intact region in the stitched BSE image rather than the most intact region with the highest content of the nonclayey mineral grains/aggregations or the regions mainly comprising the porous matrix. To be specific, sample S1 mesostructure mainly consisted of stacked clay minerals and pores, and it hardly contained nonclayey minerals. Thus, the S1-mREAs were directly estimated by the BSE images and MSSIM, and its simplified structure only contained the matrix (Figures 2(a1)–2(a5)). For sample S2, the nonclayey mineral grains were scattered in the porous matrix, and the transition parts were hardly found (Figures 2(b1)–2(b5)). Samples S3 and S4 were stiff-based mesostructures (Figures 2(c1)–2(d5)), in which sample S4 exhibited the highest nonclayey mineral content and largest grains size. For their simplified structures, the loose transition parts of sample S3 were contained in the matrix (Figures 2(c3) and 2(c5)). In contrast, the dense transition parts of sample S4 were contained in the nonclayey minerals (Figures 2(d2) and 2(d5)). Eventually, in the simplified structure of the four MRRs, the area fraction of the nonclayey minerals (φag) was 0, 36.9, 60.2, and 90.2%, respectively. In the simplified structure of the four MSRs, the φag was 0, 34.1, 53.4, and 69.4%. The LREA accuracy was affected by the stitched image area and the constituent geometric features, e.g., grain size. Comparing in the four MRRs, samples S1, S2, and S4 were limited by the mineral composition, stitched image area, and the ratio of the stitched image area to grain size, respectively. Their MRRs showed representativeness to the mesostructures, but the S3-MRR exhibited the highest representativeness. Therefore, the following sections mainly discussed the mechanical behavior of the S3-MRR, and the other three mREAs were used for comparison. To improve the structural accuracy, the boundary of the nonclayey mineral grains in the simplified S3-MRR was further delineated artificially, as shown in Figure 2(c3). As a result, its φag decreased from 60.2 to 51.6%.Figure 2
Backscattered electron image, MRR, MSR, and simplified structures of four red mudstone samples. (a), (b), (c), and (d) represent samples S1-S4, respectively. In each row, labels 1-5 represent the stitched BSE image, the segmented MRR image, the simplified MRR for simulations, the segmented MRR image, and the simplified MRR for simulations. The MRR and MSR are highlighted in the stitched BSE image.
(a)(b)(c)(d)


### 3.2. Elastic Parameter Distribution of the S3-MRR

Figure3 presents the probability density distributions (PDDs) of the elastic modulus and Poisson’s ratio of the S3-MRR at different mesh sizes. The PDD of the elastic modulus could be fitted by Gaussian normal distribution with an R2 higher than 0.9 when the mesh size decreased from 70 to 14 μm, as shown in Figure 3(a). In addition, the maximum elastic modulus shifted to a larger value and the randomness of the values increased. The PDD of the elastic modulus could not be fitted by the Gaussian normal distribution at a mesh size of 10 μm, but it could be fitted by the Weibull distribution (R2=0.15). At mesh sizes smaller than 10 μm, the PDD of the elastic modulus gradually polarized, and the polarization enhanced as the mesh size decreased. However, the PDD of the elastic modulus cannot be fitted by common statistical functions. For Poisson’s ratio, the maximum value continually decreased as the mesh size decreased, as shown in Figure 3(b). The PDD of the Poisson’s ratio at all mesh sizes could not be fitted by common statistical functions. Since the element would only contain the nonclayey minerals or matrix when the mesh size was smaller than 10 μm, the following modeling was conducted with a 10 μm mesh size or an approximate ratio of mesh size to LREA.Figure 3
Hot maps of the probability density distributions (PDD) of the (a) elastic modulus and (b) Poisson’s ratio of the S3-MRR at different mesh sizes.
(a)(b)


### 3.3. Stress–Strain Relationships and Damage Evolution of the S3-MRR

The stress–strain relationships of the S3-MRR under two loading directions are shown in Figure4. Broadly, the two S3-MRRs presented a strain-softening behavior. The model under loading in the y-axis direction possessed higher UCS than that under loading in the x-axis direction. In contrast, the model under loading in the x-axis direction had higher stress than that under loading in the y-axis direction from points B (softening stage) to C (residual stage). Figure 4 also presents the minimum principal strain (MPS) at point O (elastic stage) under loading in the y-axis direction. The nonclayey minerals deformed only slightly, while large strains were found mostly in the matrix. The MPS of the interfaces between the nonclayey minerals and the matrix (hereinafter referred to as the inclusion-matrix interfaces) was between those of the matrix and nonclayey minerals.Figure 4
Stress–strain relationship of the S3-MRR under two loading directions.The compressive and tensile damages at points A (peak strain), B (softening stage), and C (residual stage) are illustrated in Figure5. For the S3-MRR under loading in the y-axis direction, the compressive damage originated from the matrix elements of the inclusion-matrix interfaces that had an included angle of approximately 45° with the loading direction (seen in Figure 5(a1)). The compressive damage propagated along the inclusion-matrix interfaces as the strain increased and then expanded into the connected matrix regions (Figure 5(a2)). Notably, these compressive damages were more easily found in the matrix regions that had an included angle smaller than 90° with the loading direction. Finally, the S3-MRR showed a general shear failure mode with a principal fracture in the direction of the model diagonal, i.e., the compressive damaged elements with high compressive damage degrees of 0.625-0.75. Additionally, there was a short fracture perpendicular to the principal fracture. These two fractures constituted a V-shaped principal fracture (Figure 5(a3)). The tensile damage appeared later than the compressive damage (Figure 5(b1)). At the postpeak stage, the highest damage degree of the tensile damage (0.95) was higher than that of the compressive damage (0.75), as presented in Figures 5(b2) and 5(b3)). Most of the tensile damaged elements with high damage degrees were in the inclusion-matrix interfaces that had an included angle of 45° with the loading direction. In addition, these tensile damaged elements propagated along the loading direction. Nevertheless, the number of tensile damaged elements was much lower than that of the compressive damaged elements. Both the compressive and tensile damage characteristics of the model under loading in two directions were almost consistent (Figures 5(c1)–5(d3)). However, the V-shaped principal fracture under loading in the x-axis direction appeared earlier than that under loading in the y-axis direction (Figures 5(c1)–5(c3)). In addition, the inflection point of the V-shaped principal fracture was different. These results indicated that the inclusion-matrix interfaces played a significant role in the damage evolution, particularly their inclination and connectivity. Combined with the loading direction, it provided the channels and specific orientation for fracture propagation. In addition, due to the large deformation of the matrix under loads, the damage evolution was also correlated with the area and location of the matrix regions. Eventually, the mesostructural failure was mainly driven by the shear failure of the connected inclusion-matrix interfaces which were dominant in the mesostructure.Figure 5
Compressive and tensile damage evolution under two loading directions: (a) compressive and (b) tensile damages under loading in they-axis direction and (c) compressive and (d) tensile damages under the loading in the x-axis direction. For each row, the subimages with labels 1-3 represent the damage cloud charts at points A (peak strain), B (softening stage), and C (residual stage) in Figure 4.
(a)(b)(c)(d)


### 3.4. Mesomechanical Behaviors of the Simplified S3-MRRs and Random Models

It was laborious to delineate the boundary of the nonclayey mineral grains. Therefore, the failure modes of three different simplified S3-MRRs under loading in they-axis direction were investigated, as shown in Figure 6. The first model (i.e., the segmented model) was the binarized structure directly generated by the three-step mREA estimating procedure, as shown in Figure 1. The second model (i.e., the delineated model) was the delineated structure presented in Figure 2(c3). In the third model (i.e., the ellipse model), the nonclayey mineral grains were further simplified to ellipses based on the delineated model. Eventually, the area fraction of the mineral aggregations in the three models was 60.2, 51.6, and 53.2%, respectively.Figure 6
Stress–strain relationships of the (a) simplified S3-mREAs and (b) random models under loading in they-axis direction.
(a)(b)Owing to the highest nonclayey mineral content, the segmented model possessed the highest elastic modulus and UCS and the smallest peak strain (Figure6(a)). It was also the fastest to enter the residual stage and had the highest residual strength. The other two simplified S3-MRRs had almost consistent stress–strain relationship, in which the ellipse model had higher residual strength than the delineated model. Figure 7(a) shows the compressive damage of the three models at point C in Figure 6. The three models had a similar general shear failure mode with a V-shaped principal fracture. However, the damaged area of the segmented model (Figure 7(a1)) was smaller than that of the other two models (Figures 7(a2) and 7(a3)). In addition, in the ellipse model, the middle part in the principal fracture shifted from the inclusion-matrix interface to the large matrix region (Figure 7(a3)). These results indicated that the peak and residual strengths were positively correlated with the nonclayey mineral content (φag) for the structures with consistent mesostructural characteristics. However, the damage evolution and failure mode were still determined by the dominant inclusion-matrix interface.Figure 7
Compressive damage of the (a) simplified S3-mREAs and (b) random models at point C in Figure6 under loading in the y-axis direction.
(a)(b)To further investigate the damage evolution and failure mode of the mesostructures with consistentφag, three random models were generated based on the PSD of the φag in the S3-mREA (Figures 6(b) and 7(b)). The grayscale images of the three random models were enlarged based on bilinear interpolation. In these images, the blurred regions between white and black regions represent the inclusion-matrix interfaces.The UCS was negatively correlated with theφag, as shown in Figure 6(b). The first and second random models with φag=0.504 and φag=0.511 had the highest and second highest UCS. However, the first two random models had consistent softening and residual stages. The stress–strain relationship of the third model with φag=0.513 was consistent with that of the delineated model. In addition, the third random model was the fastest to enter the softening and residual stages. Although the third random model had the lowest residual strength among the three random models, its residual strength was slightly higher than that of the delineated model, which was almost consistent with the first two random models. According to the grayscale images and compressive damage at point C (Figures 7(b1) and 7(b2)), the first two random models had more local compressive damaged elements than the third random model. It was defined as a local shear failure mode with several nonpenetrative cracks. In contrast, the third random model presented a general shear failure mode with a penetrating V-shaped principal fracture (Figure 7(b3)). Compared to the delineated model, the length and width of its principal fracture were shorter and narrower. Notably, it propagated along the dominant inclusion-matrix interfaces comprising the large nonclayey minerals and the wide matrix regions between nonclayey minerals. These results illustrated that the geometric features (e.g., grain size, inclination, and spatial distribution) of the dominant inclusions significantly affected the mechanical behavior of the structures with similar mineral composition.


### 3.5. Mesomechanical Behaviors of Different Red Mudstone MRRs and MSRs

The stress–strain relationships of the MRRs and MSRs under two loading directions are presented in Figure8. Notably, the uniaxial compression of the S1-mREA was conducted by the explicit method. Broadly, the UCS was positively correlated with the φag, whereas the peak and residual strains were negatively correlated with the φag. For the MRRs, both the S2-MRR and S3-MRR presented similar stress–strain relationships under two loading directions, in which the UCS under loading in the y-axis direction was higher than that under loading in the x-axis direction (seen in Figure 8(a)). Nevertheless, the residual strengths were almost consistent under two loading directions. The S4-MRR under loading in the x-axis direction exhibited a lower UCS but a much longer softening stage with higher residual strength. For the MSRs, the MSR-UCS under loading in the y-axis direction was higher than that under loading in the x-axis direction, as depicted in Figure 8(b). Notably, the difference in UCS between the two loading directions increased as the φag increased.Figure 8
Stress–strain relationships of (a) the most representative regions (MRRs) and (b) the most special regions (MSRs) under two loading directions.
(a)(b)To further illustrate the difference in the UCS, the relationships between the meso-/macro-UCSs andφag are presented in Figure 9. On the one hand, both the UCS difference between the MRRs and the UCS difference between the MSRs increased as the φag increased (Figure 9(a)). In addition, the mean UCS difference between the MRRs and the MSRs increased as the φag increased. On the other hand, the UCS of the MRRs and MSRs can be fitted by the quadratic function with a high R2 of 0.99 (Figure 9(b)). Both the meso-UCS (the MRRs and MSRs) and the macro-UCS were positively correlated with the φag. However, with consistent φag, the meso-UCS was higher than the macro-UCS.Figure 9
Uniaxial compressive strength (UCS) of (a) the most representative regions (MRRs) and (b) the most special regions (MSRs) under two loading directions.
(a)(b)Figure10 presents the compressive damage of the MRRs and MSRs of samples S2-S4 at the beginning of the residual stage. The MRRs fractured in different but similar paths under loading in two directions. The S1-mREA presented an axisymmetric shear failure mode with X-shaped principal fractures (Figure 8(a)). It was the evolution of the V-shaped principal fracture as the strain increased. The compressive damage was mainly distributed inside the top and bottom wedge-shaped regions. The S2-MRR under loading in two directions presented the general shear failure mode with a principal fracture in the direction of the model diagonal. The fracture path was affected by the largest nonclayey mineral aggregation (Figures 10(a1) and 10(a2)). Consistent with the S3-MRR under loading in the y-axis direction (Figure 10(b1)), the S3-MRR showed a V-shaped principal fracture under loading in the x-axis direction, and the compressive damage expanded to the matrix regions with small nonclayey mineral grains (Figure 10(b2)). For the S4-MRR under loading in the y-axis direction, the fractures originated from the narrow matrix regions and propagated along its major axis (Figure 10(c1)). Although these fractures further propagated into the nonclayey minerals as the strain increased, no penetrating fractures were found. For the S4-MRR under loading in the x-axis direction, only the matrix regions at the left bottom were crushed (Figure 10(c2)). However, an approximate form of the V-shaped principal fracture was identified in the nonclayey mineral aggregations of the S4-MRRs under two loading directions. Compared to the MRRs, in the soft-based structures, the principal fracture of the MSRs propagated in the channels between the large nonclayey mineral grains (e.g., the S2-MSRs in Figures 10(a3) and 10(a4)). In stiff-based structures (e.g., the S3-MSRs and S4-MSRs, seen in Figures 10(b3), 10(b4), 10(c3), and 10(c4)), the principal fracture of the MSRs propagated along the dominant inclusion-matrix interfaces. In addition, the principal fracture was affected by the large matrix regions and loading direction, e.g., the S3-MSRs and S4-MSRs.Figure 10
Compressive damage of the most representative regions (MRRs) and the most special regions (MSRs) at the beginning of the residual stage under two loading directions.
(a)(b)(c)These results indicated that the red mudstone mREAs (the MRRs and MSRs) generally exhibited the general shear failure with a penetrating principal fracture or the local shear with several nonpenetrating local fractures. The geometric features of the dominant inclusions, e.g., large inclusions and a series of inclusions with specific arrangements, determined the failure mode. In addition, the stress–strain relationship, damage evolution, and failure mode of the mesostructure with high contents of nonclayey minerals were more sensitive to the loading direction.


### 3.6. Influences of Inclusion Geometric Features on Mesomechanical Behaviors

A series of structures were designed to further investigate the influences of the inclusion geometric features on the mesomechanical behaviors. TheLREA and the inclusion area fraction were 280 μm and 30%. The inclusions were simplified as the ellipses with 80 μm major axis length, 40 μm minor axis length, and 90 μm spacing. According to the results of the natural mREAs, this section considered the included angle between the inclusion and loading direction, inclusion type, inclusion spacing, inclusion grain size, and the large inclusion location and inclination. The uniaxial compression of these samples was conducted with the explicit method due to the convergence issues.Figure11 presents the compressive damage of the structures with consistent inclination and grain size but different included angles and spacings at the residual stage (ε=0.01). Generally, the fractures propagated along the inclusion-matrix interfaces and exhibited an included angle of ±45° with the loading direction. For the soft-based structures, the longest principal fracture more easily appeared at the lower boundaries of the stiff inclusions that were in the direction of the model diagonal. The number of the principal fractures was correlated with the included angle between the stiff inclusions and loading direction. The structures with the included angles of 22.5-45 and 135-157.5° were found to present three principal fractures (Figures 11(a2) and 11(a3) and 11(b2) and 11(b3). In contrast, the structures with the included angles of 67.5 and 112.5° exhibited a single principal fracture (Figures 11(a4) and 11(b4)). The principal fractures between the structures under loading in two directions were different because the fixed constraint at the lower left corner restrained the shear failure of the principal fractures. Notably, the a1 structure presented the general shear failure mode with a V-shaped principal fracture (Figure 11(a1)). In contrast, the b1 structure presented the local shear failure mode (Figure 11(b1)), in which the local fractures propagated along the inclusion-matrix interfaces.Figure 11
Compressive damage at the beginning of the residual stage (ε=0.01) of the structures with consistent inclusion inclination and grain size but different included angles and inclusion spacings.
(a)(b)(c)(d)(e)For the stiff-based structures, the included angle showed a more significant influence on the principal fracture propagation than that of the soft-based structures. There were several local fractures in the structures when the included angle was 0, 22.5, 90, and 157.5° (Figures11(c1), 11(c2), 11(d1), and 11(d2)). Nevertheless, these local fractures gradually transformed to a single penetrating principal fracture in the direction of the model diagonal when the included angle increased from 45 to 67.5° or from 112.5 to 135° (Figures 11(c3) and 11(c4) and 11(d3) and 11(d4). Notably, the compressive damage degree of the stiff elements between the soft elements in the fracturing direction was higher than that of these soft elements.In addition to the included angle, the inclusion spacing also affected the principal fracture propagation. For the soft-based structures, the three principal fractures changed to a single and wider principal fracture when the inclusion spacing changed from 90 to 70μm. It propagated along the lower boundaries of the stiff inclusions in the direction of the model diagonal (Figures 11(e1) and 11(e2)). Particularly, the e1 structure presented a V-shaped principal fracture, which was consistent with the S3-MRR. In contrast, for the stiff-based structures, the single principal fracture changed to three principal fractures when the inclusion spacing changed from 90 to 70 μm (Figures 11(e3) and 11(e4)).The principal fracture in the stiff/soft structures exhibited a symmetry relationship about the included angle. However, there was an opposite rule about the principal fracture form in the soft-based and stiff-based structures as the change of the included angle. To further illustrate the symmetry relationship, Figure12 presents the UCS of the structures in Figure 11. The UCS distribution was symmetric at 90° (Figure 12(a)). For the soft-based structures, the UCS was almost unchanged as the included angle increased from 0 to 45°, and then, it increased by 27.2% as the included angle increased from 45 to 90°. In contrast, the stiff-based structure UCS decreased as the included angle increased from 0 to 67.5°. It first slowly decreased as the included angle increased from 0 to 22.5°. Next, the UCS sharply decreased by 27.5% as the included angle increased from 22.5 to 45°. Then, the UCS slowly decreased as the included angle increased from 45 to 67.5°. Finally, the UCS reincreased to 9.66 MPa as the included angle increased from 67.5 to 90°, but it was still lower by 21.1% than that of the structure with the included angle of 0°. For the structures with different inclusion spacings (Figure 12(b)), the accumulation of the stiff inclusions hardly changed the UCS of the soft-based structures, whereas the accumulation of the matrix regions was conducive to increasing the UCS of the stiff-based structures.Figure 12
Uniaxial compressive strength (UCS) of the structures presented in Figure11 with different (a) inclinations and (b) spacings (Ds). The label “soft/stiff” represents the structure mainly comprising the porous matrix or nonclayey minerals. The labels on the data correspond to the subfigure labels in Figure 11.
(a)(b)Figure13 presents the compressive damage of the structures with consistent inclination but different inclusion grain sizes and included angles at the beginning of the residual stage (ε=0.01). For both the soft-based and stiff-based structures, the principal fracture was still mainly affected by the included angle. For the soft-based structures, the location and inclination of the large stiff inclusion significantly affected the principal fracture propagation (Figures 13(a1)–13(b4)). As a current example, there was a wide principal fracture that propagated along the lower boundary of the large stiff inclusion when the major axis of the large stiff inclusion was in the model diagonal (Figures 13(a2) and 13(b4)). In contrast, there were two approximate parallel principal fractures when the major axis of the large stiff inclusion was parallel to but not in the model diagonal (Figures 13(a4) and 13(b2)). One fracture was in the lower boundary of the large stiff inclusion, and the other fracture propagated along the lower boundaries of the two ordinary stiff inclusions at the lower half of the structure.Figure 13
Compressive damage at the beginning of the residual stage (ε=0.01) of the structures with consistent inclusion inclination but different grain sizes and included angles.
(a)(b)(c)(d)(e)Similar phenomena were found in the stiff-based structures (Figures13(c1)–13(d4)). The principal fracture in the direction of the model diagonal propagated straight through the large soft inclusion when its major axis was in the model diagonal. In contrast, there was a nonpenetrating fracture that propagated along the short boundary of the large soft inclusion when the major axis of the large soft inclusion was parallel to but not in the model diagonal (Figures 13(c4) and 13(d2)). Compared to the above four conditions, the principal fractures of the other four stiff-based structures were more complex. Notably, the large soft inclusion in the two structures with an included angle of 90° (Figures 13(c3) and 13(d3)) was more easily damaged than that with an included angle of 0° (Figures 13(c1) and 13(d1)). In addition, consistent with the stiff-based structures in Figure 11, for all the stiff-based structures, the compressive damage degree of the stiff elements between the soft elements in the fracturing direction was higher than that of these soft elements. For the structures where the major axis of the large soft inclusion was vertical to that of the ordinary soft inclusions, there was a V-shaped principal fracture that propagated along the stiff/soft inclusion lower boundaries (Figures 13(e1)–13(e4)).Consistent with the structures in Figure11, the symmetric relationship about the included angle was found in the principal fractures of the structures in Figure 13. To further illustrate the symmetric relationships, their UCSs are presented in Figure 14. The symmetric relationship about the included angle was also found in the UCS distribution. The location of the large inclusion hardly changed the UCS. The soft-based structure UCS increased by 37% and the stiff-based structure UCS decreased by 28% as the included angle increased from 0 to 90° (Figure 14(a)). The reincrease of the UCS at the included angle of 90° was not found in the stiff-based structures in Figure 14. For the symmetric stiff/soft structures, the major axis location of the large inclusion also hardly affected the UCS (Figure 14(b)).Figure 14
Uniaxial compressive strength (UCS) of the structures presented in Figure13 with different (a) inclinations and (b) major axis positional relationship. The label “soft/stiff” represents the structure mainly comprising the porous matrix or nonclayey minerals. The label “vertical/parallel” represents the major axis positional relationship between the large and ordinary inclusions. The labels on the data correspond to the subfigure labels in Figure 13.
(a)(b)These data indicated that the damage evolution and failure mode were sensitive to the inclusion geometric features. However, compared to the inclusion spacing and grain size and the large inclusion inclination and location, the stress–strain relationship of the red mudstone mesostructures was mainly determined by the nonclayey mineral content and the included angle between the dominant inclusions and loading direction. In addition, the included angle and inclusion spacing played more significant roles in controlling the mechanical behaviors of the structures with high contents of nonclayey minerals.


## 3.1. Mesoscopic Representative Elementary Area

Figure2 shows the stitched BSE image, MRR, MSR, and simplified structures of four samples. According to the mREA estimation, the LREA was approximately 194, 108, 280, and 582 μm for four samples, respectively. Broadly, the MRR described a relatively intact region in the stitched BSE image rather than the most intact region with the highest content of the nonclayey mineral grains/aggregations or the regions mainly comprising the porous matrix. To be specific, sample S1 mesostructure mainly consisted of stacked clay minerals and pores, and it hardly contained nonclayey minerals. Thus, the S1-mREAs were directly estimated by the BSE images and MSSIM, and its simplified structure only contained the matrix (Figures 2(a1)–2(a5)). For sample S2, the nonclayey mineral grains were scattered in the porous matrix, and the transition parts were hardly found (Figures 2(b1)–2(b5)). Samples S3 and S4 were stiff-based mesostructures (Figures 2(c1)–2(d5)), in which sample S4 exhibited the highest nonclayey mineral content and largest grains size. For their simplified structures, the loose transition parts of sample S3 were contained in the matrix (Figures 2(c3) and 2(c5)). In contrast, the dense transition parts of sample S4 were contained in the nonclayey minerals (Figures 2(d2) and 2(d5)). Eventually, in the simplified structure of the four MRRs, the area fraction of the nonclayey minerals (φag) was 0, 36.9, 60.2, and 90.2%, respectively. In the simplified structure of the four MSRs, the φag was 0, 34.1, 53.4, and 69.4%. The LREA accuracy was affected by the stitched image area and the constituent geometric features, e.g., grain size. Comparing in the four MRRs, samples S1, S2, and S4 were limited by the mineral composition, stitched image area, and the ratio of the stitched image area to grain size, respectively. Their MRRs showed representativeness to the mesostructures, but the S3-MRR exhibited the highest representativeness. Therefore, the following sections mainly discussed the mechanical behavior of the S3-MRR, and the other three mREAs were used for comparison. To improve the structural accuracy, the boundary of the nonclayey mineral grains in the simplified S3-MRR was further delineated artificially, as shown in Figure 2(c3). As a result, its φag decreased from 60.2 to 51.6%.Figure 2
Backscattered electron image, MRR, MSR, and simplified structures of four red mudstone samples. (a), (b), (c), and (d) represent samples S1-S4, respectively. In each row, labels 1-5 represent the stitched BSE image, the segmented MRR image, the simplified MRR for simulations, the segmented MRR image, and the simplified MRR for simulations. The MRR and MSR are highlighted in the stitched BSE image.
(a)(b)(c)(d)


## 3.2. Elastic Parameter Distribution of the S3-MRR

Figure3 presents the probability density distributions (PDDs) of the elastic modulus and Poisson’s ratio of the S3-MRR at different mesh sizes. The PDD of the elastic modulus could be fitted by Gaussian normal distribution with an R2 higher than 0.9 when the mesh size decreased from 70 to 14 μm, as shown in Figure 3(a). In addition, the maximum elastic modulus shifted to a larger value and the randomness of the values increased. The PDD of the elastic modulus could not be fitted by the Gaussian normal distribution at a mesh size of 10 μm, but it could be fitted by the Weibull distribution (R2=0.15). At mesh sizes smaller than 10 μm, the PDD of the elastic modulus gradually polarized, and the polarization enhanced as the mesh size decreased. However, the PDD of the elastic modulus cannot be fitted by common statistical functions. For Poisson’s ratio, the maximum value continually decreased as the mesh size decreased, as shown in Figure 3(b). The PDD of the Poisson’s ratio at all mesh sizes could not be fitted by common statistical functions. Since the element would only contain the nonclayey minerals or matrix when the mesh size was smaller than 10 μm, the following modeling was conducted with a 10 μm mesh size or an approximate ratio of mesh size to LREA.Figure 3
Hot maps of the probability density distributions (PDD) of the (a) elastic modulus and (b) Poisson’s ratio of the S3-MRR at different mesh sizes.
(a)(b)


## 3.3. Stress–Strain Relationships and Damage Evolution of the S3-MRR

The stress–strain relationships of the S3-MRR under two loading directions are shown in Figure4. Broadly, the two S3-MRRs presented a strain-softening behavior. The model under loading in the y-axis direction possessed higher UCS than that under loading in the x-axis direction. In contrast, the model under loading in the x-axis direction had higher stress than that under loading in the y-axis direction from points B (softening stage) to C (residual stage). Figure 4 also presents the minimum principal strain (MPS) at point O (elastic stage) under loading in the y-axis direction. The nonclayey minerals deformed only slightly, while large strains were found mostly in the matrix. The MPS of the interfaces between the nonclayey minerals and the matrix (hereinafter referred to as the inclusion-matrix interfaces) was between those of the matrix and nonclayey minerals.Figure 4
Stress–strain relationship of the S3-MRR under two loading directions.The compressive and tensile damages at points A (peak strain), B (softening stage), and C (residual stage) are illustrated in Figure5. For the S3-MRR under loading in the y-axis direction, the compressive damage originated from the matrix elements of the inclusion-matrix interfaces that had an included angle of approximately 45° with the loading direction (seen in Figure 5(a1)). The compressive damage propagated along the inclusion-matrix interfaces as the strain increased and then expanded into the connected matrix regions (Figure 5(a2)). Notably, these compressive damages were more easily found in the matrix regions that had an included angle smaller than 90° with the loading direction. Finally, the S3-MRR showed a general shear failure mode with a principal fracture in the direction of the model diagonal, i.e., the compressive damaged elements with high compressive damage degrees of 0.625-0.75. Additionally, there was a short fracture perpendicular to the principal fracture. These two fractures constituted a V-shaped principal fracture (Figure 5(a3)). The tensile damage appeared later than the compressive damage (Figure 5(b1)). At the postpeak stage, the highest damage degree of the tensile damage (0.95) was higher than that of the compressive damage (0.75), as presented in Figures 5(b2) and 5(b3)). Most of the tensile damaged elements with high damage degrees were in the inclusion-matrix interfaces that had an included angle of 45° with the loading direction. In addition, these tensile damaged elements propagated along the loading direction. Nevertheless, the number of tensile damaged elements was much lower than that of the compressive damaged elements. Both the compressive and tensile damage characteristics of the model under loading in two directions were almost consistent (Figures 5(c1)–5(d3)). However, the V-shaped principal fracture under loading in the x-axis direction appeared earlier than that under loading in the y-axis direction (Figures 5(c1)–5(c3)). In addition, the inflection point of the V-shaped principal fracture was different. These results indicated that the inclusion-matrix interfaces played a significant role in the damage evolution, particularly their inclination and connectivity. Combined with the loading direction, it provided the channels and specific orientation for fracture propagation. In addition, due to the large deformation of the matrix under loads, the damage evolution was also correlated with the area and location of the matrix regions. Eventually, the mesostructural failure was mainly driven by the shear failure of the connected inclusion-matrix interfaces which were dominant in the mesostructure.Figure 5
Compressive and tensile damage evolution under two loading directions: (a) compressive and (b) tensile damages under loading in they-axis direction and (c) compressive and (d) tensile damages under the loading in the x-axis direction. For each row, the subimages with labels 1-3 represent the damage cloud charts at points A (peak strain), B (softening stage), and C (residual stage) in Figure 4.
(a)(b)(c)(d)


## 3.4. Mesomechanical Behaviors of the Simplified S3-MRRs and Random Models

It was laborious to delineate the boundary of the nonclayey mineral grains. Therefore, the failure modes of three different simplified S3-MRRs under loading in they-axis direction were investigated, as shown in Figure 6. The first model (i.e., the segmented model) was the binarized structure directly generated by the three-step mREA estimating procedure, as shown in Figure 1. The second model (i.e., the delineated model) was the delineated structure presented in Figure 2(c3). In the third model (i.e., the ellipse model), the nonclayey mineral grains were further simplified to ellipses based on the delineated model. Eventually, the area fraction of the mineral aggregations in the three models was 60.2, 51.6, and 53.2%, respectively.Figure 6
Stress–strain relationships of the (a) simplified S3-mREAs and (b) random models under loading in they-axis direction.
(a)(b)Owing to the highest nonclayey mineral content, the segmented model possessed the highest elastic modulus and UCS and the smallest peak strain (Figure6(a)). It was also the fastest to enter the residual stage and had the highest residual strength. The other two simplified S3-MRRs had almost consistent stress–strain relationship, in which the ellipse model had higher residual strength than the delineated model. Figure 7(a) shows the compressive damage of the three models at point C in Figure 6. The three models had a similar general shear failure mode with a V-shaped principal fracture. However, the damaged area of the segmented model (Figure 7(a1)) was smaller than that of the other two models (Figures 7(a2) and 7(a3)). In addition, in the ellipse model, the middle part in the principal fracture shifted from the inclusion-matrix interface to the large matrix region (Figure 7(a3)). These results indicated that the peak and residual strengths were positively correlated with the nonclayey mineral content (φag) for the structures with consistent mesostructural characteristics. However, the damage evolution and failure mode were still determined by the dominant inclusion-matrix interface.Figure 7
Compressive damage of the (a) simplified S3-mREAs and (b) random models at point C in Figure6 under loading in the y-axis direction.
(a)(b)To further investigate the damage evolution and failure mode of the mesostructures with consistentφag, three random models were generated based on the PSD of the φag in the S3-mREA (Figures 6(b) and 7(b)). The grayscale images of the three random models were enlarged based on bilinear interpolation. In these images, the blurred regions between white and black regions represent the inclusion-matrix interfaces.The UCS was negatively correlated with theφag, as shown in Figure 6(b). The first and second random models with φag=0.504 and φag=0.511 had the highest and second highest UCS. However, the first two random models had consistent softening and residual stages. The stress–strain relationship of the third model with φag=0.513 was consistent with that of the delineated model. In addition, the third random model was the fastest to enter the softening and residual stages. Although the third random model had the lowest residual strength among the three random models, its residual strength was slightly higher than that of the delineated model, which was almost consistent with the first two random models. According to the grayscale images and compressive damage at point C (Figures 7(b1) and 7(b2)), the first two random models had more local compressive damaged elements than the third random model. It was defined as a local shear failure mode with several nonpenetrative cracks. In contrast, the third random model presented a general shear failure mode with a penetrating V-shaped principal fracture (Figure 7(b3)). Compared to the delineated model, the length and width of its principal fracture were shorter and narrower. Notably, it propagated along the dominant inclusion-matrix interfaces comprising the large nonclayey minerals and the wide matrix regions between nonclayey minerals. These results illustrated that the geometric features (e.g., grain size, inclination, and spatial distribution) of the dominant inclusions significantly affected the mechanical behavior of the structures with similar mineral composition.


## 3.5. Mesomechanical Behaviors of Different Red Mudstone MRRs and MSRs

The stress–strain relationships of the MRRs and MSRs under two loading directions are presented in Figure8. Notably, the uniaxial compression of the S1-mREA was conducted by the explicit method. Broadly, the UCS was positively correlated with the φag, whereas the peak and residual strains were negatively correlated with the φag. For the MRRs, both the S2-MRR and S3-MRR presented similar stress–strain relationships under two loading directions, in which the UCS under loading in the y-axis direction was higher than that under loading in the x-axis direction (seen in Figure 8(a)). Nevertheless, the residual strengths were almost consistent under two loading directions. The S4-MRR under loading in the x-axis direction exhibited a lower UCS but a much longer softening stage with higher residual strength. For the MSRs, the MSR-UCS under loading in the y-axis direction was higher than that under loading in the x-axis direction, as depicted in Figure 8(b). Notably, the difference in UCS between the two loading directions increased as the φag increased.Figure 8
Stress–strain relationships of (a) the most representative regions (MRRs) and (b) the most special regions (MSRs) under two loading directions.
(a)(b)To further illustrate the difference in the UCS, the relationships between the meso-/macro-UCSs andφag are presented in Figure 9. On the one hand, both the UCS difference between the MRRs and the UCS difference between the MSRs increased as the φag increased (Figure 9(a)). In addition, the mean UCS difference between the MRRs and the MSRs increased as the φag increased. On the other hand, the UCS of the MRRs and MSRs can be fitted by the quadratic function with a high R2 of 0.99 (Figure 9(b)). Both the meso-UCS (the MRRs and MSRs) and the macro-UCS were positively correlated with the φag. However, with consistent φag, the meso-UCS was higher than the macro-UCS.Figure 9
Uniaxial compressive strength (UCS) of (a) the most representative regions (MRRs) and (b) the most special regions (MSRs) under two loading directions.
(a)(b)Figure10 presents the compressive damage of the MRRs and MSRs of samples S2-S4 at the beginning of the residual stage. The MRRs fractured in different but similar paths under loading in two directions. The S1-mREA presented an axisymmetric shear failure mode with X-shaped principal fractures (Figure 8(a)). It was the evolution of the V-shaped principal fracture as the strain increased. The compressive damage was mainly distributed inside the top and bottom wedge-shaped regions. The S2-MRR under loading in two directions presented the general shear failure mode with a principal fracture in the direction of the model diagonal. The fracture path was affected by the largest nonclayey mineral aggregation (Figures 10(a1) and 10(a2)). Consistent with the S3-MRR under loading in the y-axis direction (Figure 10(b1)), the S3-MRR showed a V-shaped principal fracture under loading in the x-axis direction, and the compressive damage expanded to the matrix regions with small nonclayey mineral grains (Figure 10(b2)). For the S4-MRR under loading in the y-axis direction, the fractures originated from the narrow matrix regions and propagated along its major axis (Figure 10(c1)). Although these fractures further propagated into the nonclayey minerals as the strain increased, no penetrating fractures were found. For the S4-MRR under loading in the x-axis direction, only the matrix regions at the left bottom were crushed (Figure 10(c2)). However, an approximate form of the V-shaped principal fracture was identified in the nonclayey mineral aggregations of the S4-MRRs under two loading directions. Compared to the MRRs, in the soft-based structures, the principal fracture of the MSRs propagated in the channels between the large nonclayey mineral grains (e.g., the S2-MSRs in Figures 10(a3) and 10(a4)). In stiff-based structures (e.g., the S3-MSRs and S4-MSRs, seen in Figures 10(b3), 10(b4), 10(c3), and 10(c4)), the principal fracture of the MSRs propagated along the dominant inclusion-matrix interfaces. In addition, the principal fracture was affected by the large matrix regions and loading direction, e.g., the S3-MSRs and S4-MSRs.Figure 10
Compressive damage of the most representative regions (MRRs) and the most special regions (MSRs) at the beginning of the residual stage under two loading directions.
(a)(b)(c)These results indicated that the red mudstone mREAs (the MRRs and MSRs) generally exhibited the general shear failure with a penetrating principal fracture or the local shear with several nonpenetrating local fractures. The geometric features of the dominant inclusions, e.g., large inclusions and a series of inclusions with specific arrangements, determined the failure mode. In addition, the stress–strain relationship, damage evolution, and failure mode of the mesostructure with high contents of nonclayey minerals were more sensitive to the loading direction.


## 3.6. Influences of Inclusion Geometric Features on Mesomechanical Behaviors

A series of structures were designed to further investigate the influences of the inclusion geometric features on the mesomechanical behaviors. TheLREA and the inclusion area fraction were 280 μm and 30%. The inclusions were simplified as the ellipses with 80 μm major axis length, 40 μm minor axis length, and 90 μm spacing. According to the results of the natural mREAs, this section considered the included angle between the inclusion and loading direction, inclusion type, inclusion spacing, inclusion grain size, and the large inclusion location and inclination. The uniaxial compression of these samples was conducted with the explicit method due to the convergence issues.Figure11 presents the compressive damage of the structures with consistent inclination and grain size but different included angles and spacings at the residual stage (ε=0.01). Generally, the fractures propagated along the inclusion-matrix interfaces and exhibited an included angle of ±45° with the loading direction. For the soft-based structures, the longest principal fracture more easily appeared at the lower boundaries of the stiff inclusions that were in the direction of the model diagonal. The number of the principal fractures was correlated with the included angle between the stiff inclusions and loading direction. The structures with the included angles of 22.5-45 and 135-157.5° were found to present three principal fractures (Figures 11(a2) and 11(a3) and 11(b2) and 11(b3). In contrast, the structures with the included angles of 67.5 and 112.5° exhibited a single principal fracture (Figures 11(a4) and 11(b4)). The principal fractures between the structures under loading in two directions were different because the fixed constraint at the lower left corner restrained the shear failure of the principal fractures. Notably, the a1 structure presented the general shear failure mode with a V-shaped principal fracture (Figure 11(a1)). In contrast, the b1 structure presented the local shear failure mode (Figure 11(b1)), in which the local fractures propagated along the inclusion-matrix interfaces.Figure 11
Compressive damage at the beginning of the residual stage (ε=0.01) of the structures with consistent inclusion inclination and grain size but different included angles and inclusion spacings.
(a)(b)(c)(d)(e)For the stiff-based structures, the included angle showed a more significant influence on the principal fracture propagation than that of the soft-based structures. There were several local fractures in the structures when the included angle was 0, 22.5, 90, and 157.5° (Figures11(c1), 11(c2), 11(d1), and 11(d2)). Nevertheless, these local fractures gradually transformed to a single penetrating principal fracture in the direction of the model diagonal when the included angle increased from 45 to 67.5° or from 112.5 to 135° (Figures 11(c3) and 11(c4) and 11(d3) and 11(d4). Notably, the compressive damage degree of the stiff elements between the soft elements in the fracturing direction was higher than that of these soft elements.In addition to the included angle, the inclusion spacing also affected the principal fracture propagation. For the soft-based structures, the three principal fractures changed to a single and wider principal fracture when the inclusion spacing changed from 90 to 70μm. It propagated along the lower boundaries of the stiff inclusions in the direction of the model diagonal (Figures 11(e1) and 11(e2)). Particularly, the e1 structure presented a V-shaped principal fracture, which was consistent with the S3-MRR. In contrast, for the stiff-based structures, the single principal fracture changed to three principal fractures when the inclusion spacing changed from 90 to 70 μm (Figures 11(e3) and 11(e4)).The principal fracture in the stiff/soft structures exhibited a symmetry relationship about the included angle. However, there was an opposite rule about the principal fracture form in the soft-based and stiff-based structures as the change of the included angle. To further illustrate the symmetry relationship, Figure12 presents the UCS of the structures in Figure 11. The UCS distribution was symmetric at 90° (Figure 12(a)). For the soft-based structures, the UCS was almost unchanged as the included angle increased from 0 to 45°, and then, it increased by 27.2% as the included angle increased from 45 to 90°. In contrast, the stiff-based structure UCS decreased as the included angle increased from 0 to 67.5°. It first slowly decreased as the included angle increased from 0 to 22.5°. Next, the UCS sharply decreased by 27.5% as the included angle increased from 22.5 to 45°. Then, the UCS slowly decreased as the included angle increased from 45 to 67.5°. Finally, the UCS reincreased to 9.66 MPa as the included angle increased from 67.5 to 90°, but it was still lower by 21.1% than that of the structure with the included angle of 0°. For the structures with different inclusion spacings (Figure 12(b)), the accumulation of the stiff inclusions hardly changed the UCS of the soft-based structures, whereas the accumulation of the matrix regions was conducive to increasing the UCS of the stiff-based structures.Figure 12
Uniaxial compressive strength (UCS) of the structures presented in Figure11 with different (a) inclinations and (b) spacings (Ds). The label “soft/stiff” represents the structure mainly comprising the porous matrix or nonclayey minerals. The labels on the data correspond to the subfigure labels in Figure 11.
(a)(b)Figure13 presents the compressive damage of the structures with consistent inclination but different inclusion grain sizes and included angles at the beginning of the residual stage (ε=0.01). For both the soft-based and stiff-based structures, the principal fracture was still mainly affected by the included angle. For the soft-based structures, the location and inclination of the large stiff inclusion significantly affected the principal fracture propagation (Figures 13(a1)–13(b4)). As a current example, there was a wide principal fracture that propagated along the lower boundary of the large stiff inclusion when the major axis of the large stiff inclusion was in the model diagonal (Figures 13(a2) and 13(b4)). In contrast, there were two approximate parallel principal fractures when the major axis of the large stiff inclusion was parallel to but not in the model diagonal (Figures 13(a4) and 13(b2)). One fracture was in the lower boundary of the large stiff inclusion, and the other fracture propagated along the lower boundaries of the two ordinary stiff inclusions at the lower half of the structure.Figure 13
Compressive damage at the beginning of the residual stage (ε=0.01) of the structures with consistent inclusion inclination but different grain sizes and included angles.
(a)(b)(c)(d)(e)Similar phenomena were found in the stiff-based structures (Figures13(c1)–13(d4)). The principal fracture in the direction of the model diagonal propagated straight through the large soft inclusion when its major axis was in the model diagonal. In contrast, there was a nonpenetrating fracture that propagated along the short boundary of the large soft inclusion when the major axis of the large soft inclusion was parallel to but not in the model diagonal (Figures 13(c4) and 13(d2)). Compared to the above four conditions, the principal fractures of the other four stiff-based structures were more complex. Notably, the large soft inclusion in the two structures with an included angle of 90° (Figures 13(c3) and 13(d3)) was more easily damaged than that with an included angle of 0° (Figures 13(c1) and 13(d1)). In addition, consistent with the stiff-based structures in Figure 11, for all the stiff-based structures, the compressive damage degree of the stiff elements between the soft elements in the fracturing direction was higher than that of these soft elements. For the structures where the major axis of the large soft inclusion was vertical to that of the ordinary soft inclusions, there was a V-shaped principal fracture that propagated along the stiff/soft inclusion lower boundaries (Figures 13(e1)–13(e4)).Consistent with the structures in Figure11, the symmetric relationship about the included angle was found in the principal fractures of the structures in Figure 13. To further illustrate the symmetric relationships, their UCSs are presented in Figure 14. The symmetric relationship about the included angle was also found in the UCS distribution. The location of the large inclusion hardly changed the UCS. The soft-based structure UCS increased by 37% and the stiff-based structure UCS decreased by 28% as the included angle increased from 0 to 90° (Figure 14(a)). The reincrease of the UCS at the included angle of 90° was not found in the stiff-based structures in Figure 14. For the symmetric stiff/soft structures, the major axis location of the large inclusion also hardly affected the UCS (Figure 14(b)).Figure 14
Uniaxial compressive strength (UCS) of the structures presented in Figure13 with different (a) inclinations and (b) major axis positional relationship. The label “soft/stiff” represents the structure mainly comprising the porous matrix or nonclayey minerals. The label “vertical/parallel” represents the major axis positional relationship between the large and ordinary inclusions. The labels on the data correspond to the subfigure labels in Figure 13.
(a)(b)These data indicated that the damage evolution and failure mode were sensitive to the inclusion geometric features. However, compared to the inclusion spacing and grain size and the large inclusion inclination and location, the stress–strain relationship of the red mudstone mesostructures was mainly determined by the nonclayey mineral content and the included angle between the dominant inclusions and loading direction. In addition, the included angle and inclusion spacing played more significant roles in controlling the mechanical behaviors of the structures with high contents of nonclayey minerals.


## 4. Discussion

### 4.1. Structural Damage Mechanism

The mREA and linearly superimposed model, i.e., the mREA mechanical model, revealed the damage mechanism of the red mudstone mesostructures under uniaxial compression. In the prepeak stage, the strain difference between constituents was essentially caused by the contrast in stiffness [4]. Our study indicates that the inclusion-matrix interfaces are the weakest regions in the postpeak stage, particularly the interfaces constituted by the dominant inclusions, e.g., large inclusions and a series of inclusions with specific arrangements. As the strain increases, the shear failure of these interfaces causes the fractures to appear on the matrix side and then propagate along the interfaces. This finding is consistent with the literature [22] that the strain concentration zones (SCZs) are located in the matrix regions around the stiff inclusions. Our findings further illustrate that the damage evolution is sensitive to the inclusion geometric features, e.g., spacing, grain size, and the large inclusion inclination and location. However, the UCS of the structures with consistent mineral composition is mainly determined by the included angle between the dominant inclusion inclination and loading direction. This finding extends the conclusion in the literature [23] suggesting that the failure strength is strongly dependent on the included angle between the loading direction and bedding planes. The inclusion-matrix interfaces constituted by the dominant inclusions are the mesoscopic structural planes. Driven by this damage mechanism, the red mudstone mesostructures present two shear failure modes, i.e., the general shear failure with a penetrating principal fracture and the local shear with several nonpenetrating local fractures. However, for the stiff-based structures, the second failure mode is established based on ignoring the grain boundaries in nonclayey mineral aggregations. At present, the impacts of interfaces between the nonclayey minerals on the fracture propagation of stiff-based mREAs are not clear. In the future, the adhesion difference between the nonclayey minerals and porous matrix based on atomic force microscopy may provide some suggestions. Notwithstanding its limitation, this study suggests that the macrostructural failure of red mudstones is driven by the accumulation and expansion of the shear failure of a series of the hard–soft interface at multiple scales.


### 4.2. Compatibility between Mesoscopic and Macroscopic Structural Failure

The mREA mechanical model builds a rational bridge to predict meso-to-macrostructural failure. On the one hand, although the weight of the MRRs and MSRs on the mechanical behaviors of the red mudstone mesostructures was not examined, the damage characteristics found in the mREAs could be applicable to multiscale structures, e.g., the different damage evolution of the S3-MSR under two loading directions. The MRR provides a compromised understanding of the damage evolution of the red mudstone mesostructure. However, for example, the results of the S3-MSR under loading in thex-axis direction indicate how the large matrix regions influence the principal fracture propagation. It can be speculated that the MSRs are the significant nodes of the principal fractures in a larger mesostructure under specific loading directions, while the bedding planes and natural fractures will exhibit their influences on the principal fracture propagation of the structure at a much larger scale. On the other hand, both the macro-UCS and mREA-UCS are positively related to the nonclayey mineral content, in which the mREA-UCS was higher than its macro-UCS. However, it is mainly because the segmented model has a higher content of nonclayey minerals than the realistic mesostructure. In addition, the transition parts were defined as having the same material properties as the nonclayey minerals. The mechanical properties of the constituents should be different in various red mudstones. Therefore, there is insufficient evidence to show that the UCS compatibility between the macroscopic and mesoscopic is positively correlated with mesostructural representativeness. Despite this, considering the structural planes at larger scales, the meso-UCS should be higher than the macro-UCS.


### 4.3. Potentiality of the mREA Mechanical Model

Mesomechanical experiments are still the most suitable approach to reveal the mesoscopic damage evolution of clay rocks. The mREA mechanical model will be a convenient and rational reference for the mesostructural mechanical testing of clay rocks and rock-like materials. On the one hand, it provides a prior judgment of the ROI representativeness to the mesostructure. Based on the damage mechanism in this study, the failure mode can be also rapidly estimated according to the captured mesostructure. On the other hand, the mREA mechanical model can simulate the damage evolution and failure mode. For instance, destructive experiments were generally prohibited in ESEM because the owners of these instruments need to prevent the fragments from contaminating the lens. It means that the mechanical experiments would be limited to the elastic or nonlinear strengthening stage. Since the mREA provides a relatively large area that can be found again, the mechanical test can be conducted outside the ESEM apparatus, while the BSE images of the original and destroyed mesostructures can be recorded inside the ESEM apparatus [21]. The damage evolution can be rapidly adjusted to match the realistic situation based on the high computational efficiency which was a few minutes to ten minutes with an i5-10600KF CPU and 16 GB memory in this study.


### 4.4. Pore Influences on Damage Evolution and Seepage

The pore influence on the damage evolution and seepage should be discussed. On the one hand, this gap restrains the failure mode rationality of the mREAs with high contents of nonclayey minerals. Most of the pores captured are the intergranular pores between nonclayey minerals and between the mineral grains in the matrix. The former geometric features and spatial distribution are determined by that of the nonclayey minerals, while the latter is commonly too small in the mREA. Literature [22] indicated that the SCZs were attributed to a large number of nano- to microscale pores in the clay matrix surrounding the stiff inclusion zones. Although the damage evolution of the mesostructure is mainly correlated with the dominant inclusions, the influences of the two types of pores on the damage mechanism should be identified. The pores and cracks possibly induce partial grain boundaries possibly to fracture and slide earlier than the inclusion-matrix interfaces.On the other hand, the weak water resistance causes the strength reduction, swelling, and disintegration of the red mudstone. It troubles the structural safety of foundations, slopes, and tunnels, inducing building inclination, landslides, and floor heave. At the mesoscale captured, the pores between nonclayey minerals and the pores in the matrix are one of the major seepage channels for the water in the mesostructure, as the pores of four samples shown in Figures15(a1)–15(d2). In addition to the swelling of clay minerals, the small particles in the matrix can be immigrated under the seepage according to the water–rock–silt flow model [44–46], particularly for the loose matrix such as sample S2 (Figures 15(b1) and 15(b2)). The fracture propagation under loads will aggravate the particle immigration under seepage. Meanwhile, the seepage will accelerate fracture propagation. However, the mREA estimated by the constituent geometric features may not be directly used to study the representativeness of the seepage. The mREA of the seepage is more likely to be correlated with large porous matrix regions.Figure 15
Pores in the four red mudstone samples: (a) sample S1, (b) sample S2, (c) sample S3, and (d) sample S4.
(a)(b)(c)(d)


## 4.1. Structural Damage Mechanism

The mREA and linearly superimposed model, i.e., the mREA mechanical model, revealed the damage mechanism of the red mudstone mesostructures under uniaxial compression. In the prepeak stage, the strain difference between constituents was essentially caused by the contrast in stiffness [4]. Our study indicates that the inclusion-matrix interfaces are the weakest regions in the postpeak stage, particularly the interfaces constituted by the dominant inclusions, e.g., large inclusions and a series of inclusions with specific arrangements. As the strain increases, the shear failure of these interfaces causes the fractures to appear on the matrix side and then propagate along the interfaces. This finding is consistent with the literature [22] that the strain concentration zones (SCZs) are located in the matrix regions around the stiff inclusions. Our findings further illustrate that the damage evolution is sensitive to the inclusion geometric features, e.g., spacing, grain size, and the large inclusion inclination and location. However, the UCS of the structures with consistent mineral composition is mainly determined by the included angle between the dominant inclusion inclination and loading direction. This finding extends the conclusion in the literature [23] suggesting that the failure strength is strongly dependent on the included angle between the loading direction and bedding planes. The inclusion-matrix interfaces constituted by the dominant inclusions are the mesoscopic structural planes. Driven by this damage mechanism, the red mudstone mesostructures present two shear failure modes, i.e., the general shear failure with a penetrating principal fracture and the local shear with several nonpenetrating local fractures. However, for the stiff-based structures, the second failure mode is established based on ignoring the grain boundaries in nonclayey mineral aggregations. At present, the impacts of interfaces between the nonclayey minerals on the fracture propagation of stiff-based mREAs are not clear. In the future, the adhesion difference between the nonclayey minerals and porous matrix based on atomic force microscopy may provide some suggestions. Notwithstanding its limitation, this study suggests that the macrostructural failure of red mudstones is driven by the accumulation and expansion of the shear failure of a series of the hard–soft interface at multiple scales.


## 4.2. Compatibility between Mesoscopic and Macroscopic Structural Failure

The mREA mechanical model builds a rational bridge to predict meso-to-macrostructural failure. On the one hand, although the weight of the MRRs and MSRs on the mechanical behaviors of the red mudstone mesostructures was not examined, the damage characteristics found in the mREAs could be applicable to multiscale structures, e.g., the different damage evolution of the S3-MSR under two loading directions. The MRR provides a compromised understanding of the damage evolution of the red mudstone mesostructure. However, for example, the results of the S3-MSR under loading in thex-axis direction indicate how the large matrix regions influence the principal fracture propagation. It can be speculated that the MSRs are the significant nodes of the principal fractures in a larger mesostructure under specific loading directions, while the bedding planes and natural fractures will exhibit their influences on the principal fracture propagation of the structure at a much larger scale. On the other hand, both the macro-UCS and mREA-UCS are positively related to the nonclayey mineral content, in which the mREA-UCS was higher than its macro-UCS. However, it is mainly because the segmented model has a higher content of nonclayey minerals than the realistic mesostructure. In addition, the transition parts were defined as having the same material properties as the nonclayey minerals. The mechanical properties of the constituents should be different in various red mudstones. Therefore, there is insufficient evidence to show that the UCS compatibility between the macroscopic and mesoscopic is positively correlated with mesostructural representativeness. Despite this, considering the structural planes at larger scales, the meso-UCS should be higher than the macro-UCS.


## 4.3. Potentiality of the mREA Mechanical Model

Mesomechanical experiments are still the most suitable approach to reveal the mesoscopic damage evolution of clay rocks. The mREA mechanical model will be a convenient and rational reference for the mesostructural mechanical testing of clay rocks and rock-like materials. On the one hand, it provides a prior judgment of the ROI representativeness to the mesostructure. Based on the damage mechanism in this study, the failure mode can be also rapidly estimated according to the captured mesostructure. On the other hand, the mREA mechanical model can simulate the damage evolution and failure mode. For instance, destructive experiments were generally prohibited in ESEM because the owners of these instruments need to prevent the fragments from contaminating the lens. It means that the mechanical experiments would be limited to the elastic or nonlinear strengthening stage. Since the mREA provides a relatively large area that can be found again, the mechanical test can be conducted outside the ESEM apparatus, while the BSE images of the original and destroyed mesostructures can be recorded inside the ESEM apparatus [21]. The damage evolution can be rapidly adjusted to match the realistic situation based on the high computational efficiency which was a few minutes to ten minutes with an i5-10600KF CPU and 16 GB memory in this study.


## 4.4. Pore Influences on Damage Evolution and Seepage

The pore influence on the damage evolution and seepage should be discussed. On the one hand, this gap restrains the failure mode rationality of the mREAs with high contents of nonclayey minerals. Most of the pores captured are the intergranular pores between nonclayey minerals and between the mineral grains in the matrix. The former geometric features and spatial distribution are determined by that of the nonclayey minerals, while the latter is commonly too small in the mREA. Literature [22] indicated that the SCZs were attributed to a large number of nano- to microscale pores in the clay matrix surrounding the stiff inclusion zones. Although the damage evolution of the mesostructure is mainly correlated with the dominant inclusions, the influences of the two types of pores on the damage mechanism should be identified. The pores and cracks possibly induce partial grain boundaries possibly to fracture and slide earlier than the inclusion-matrix interfaces.On the other hand, the weak water resistance causes the strength reduction, swelling, and disintegration of the red mudstone. It troubles the structural safety of foundations, slopes, and tunnels, inducing building inclination, landslides, and floor heave. At the mesoscale captured, the pores between nonclayey minerals and the pores in the matrix are one of the major seepage channels for the water in the mesostructure, as the pores of four samples shown in Figures15(a1)–15(d2). In addition to the swelling of clay minerals, the small particles in the matrix can be immigrated under the seepage according to the water–rock–silt flow model [44–46], particularly for the loose matrix such as sample S2 (Figures 15(b1) and 15(b2)). The fracture propagation under loads will aggravate the particle immigration under seepage. Meanwhile, the seepage will accelerate fracture propagation. However, the mREA estimated by the constituent geometric features may not be directly used to study the representativeness of the seepage. The mREA of the seepage is more likely to be correlated with large porous matrix regions.Figure 15
Pores in the four red mudstone samples: (a) sample S1, (b) sample S2, (c) sample S3, and (d) sample S4.
(a)(b)(c)(d)


## 5. Conclusions

The main objective of this study was to characterize the influence of the composition, geometric features, and spatial distribution of the constituents on the mesostructural mechanical behavior of the red mudstones through the mREA estimating procedure, MEEM, linearly superimposed model, and CDP model. The probability density distribution of the elastic parameters, stress–strain relationships, strain field, damage evolution, and failure mode were investigated and discussed. The basic conclusions were drawn as follows:(1)
The inclusion-matrix interfaces constituted by the dominant inclusions are the weakest regions in the postpeak stage, e.g., large inclusions and a series of inclusions with specific arrangements. As the strain increases, the shear failure of these interfaces causes the fractures to appear on the matrix side and then propagate along the interfaces. Driven by this damage mechanism, the mREAs generally present two shear failure modes, i.e., the general shear failure with a penetrating principal fracture and the local shear with several nonpenetrating local fractures.(2)
The damage evolution and failure mode are sensitive to the inclusion geometric features, e.g., inclination, spacing, grain size, and the location of the large inclusion. In contrast, the UCS of the structures with consistent mineral composition is mainly determined by the included angle between the dominant inclusion inclination and loading direction. With consistent composition, an increase in the included angle from 0 to 90° results in an increase in the UCS of the muddy facies and a decrease in the sandy facies. In addition, compared to the muddy structures, the included angle shows a more significant influence on the damage evolution and the UCS of the sandy structures.(3)
The mREA mechanical model builds a rational bridge to predict the mechanical behavior of rock and rock-like materials from the mesoscale to the macroscale. The inclusion-matrix interfaces comprising the dominant inclusions are the mesoscopic structural planes in the red mudstone structure. Considering the potential structural planes at larger scales, the structural failure of red mudstone is driven by the accumulation and expansion of the shear failure of a series of the hard–soft interface at multiple scales.This study improves our understanding of the structural failure mechanism of mudstones and will be a convenient and rational reference for the mesomechanical testing of clay rocks and rock-like materials.



---

*Source: 2106513-2023-03-04.xml*