Prediction of the Anisotropic Effective Moduli of Shales based on the Mori–Tanaka Model and the Digital Core Technique

Natural rocks belong to the class of polymineral composite materials with complex microstructures. Such strong heterogeneity of rocks makes it difficult to estimate the effective moduli by traditional models in theory. In the present study, a Mori–Tanaka (MT) model considering the shape and orientation of inclusion minerals obtained by micro-computed tomography (micro-CT) is established, and then it is applied to evaluate the anisotropic parameters of shales. In the MT model, the principal radii and Eulerian angles of the ellipsoidal inclusion are obtained by solving its inertia matrix through micro-CT. According to this inclusion information, we determined the statistics on the ratio of average principal radii and the distribution of Eulerian angles of inclusions with different minerals. In what follows, the effective elastic stiffness matrix of shale samples is predicted by the MT model, and the corresponding digital core is input for finite element method (FEM) analysis to verify the accuracy of the theoretical results. It is shown that the anisotropy of the elastic stiffness matrix predicted by the MT model and FEM is consistent under two sizes of representative volume elements. These findings highlight the potential of this approach for applications in rock mechanics, civil engineering, and oil exploitation, among others.


Introduction
Unconventional reservoir rocks, such as tight sandstone and shale, are characterized by multimineral composition and complex microstructures. The resulting strong heterogeneity and complex mechanical properties can make it challenging to predict the overall elastic properties of reservoir rocks (Prasad et al., 2009;Jian et al., 2020). In particular, the directionally distributed pores or other minerals in tight sandstones or shales from deep reservoirs as deviatoric stress fields often make reservoir rocks anisotropic, and this fact significantly increases the difficulty in theoretically estimating effective moduli by conventional models (Shapiro & Kaselow, 2005;Gui et al., 2018;Zhang et al., 2020).
Many models have been proposed to consider the anisotropy of reservoir rocks (Vernik & Nur, 1992;Hudson, 1986;Cheng, 1993), but the consideration of the shape and orientation of minerals is missing. In practice, the theory of mesomechanics starting from the Eshelby tensor (Eshelby, 1957) has been extensively applied in composite materials (Berryman, 1995;Raju et al., 2018), among which the Mori-Tanaka (MT) model (Mori & Tanaka, 1974) is one of the most representative methods (Park et al., 2020) and it has been put to use in many fields (Shen & Wang, 2012;Wang et al., 2017). Zimmerman and David (1991, 2011a, 2011b have systematically analyzed many effective medium models including the MT model, and they deduced the pore compressibility and shear compliance, which has greatly promoted the application of the MT model in rock physics. In succession, Deng et al. (2015) used the MT model to calculate the effective moduli of a dry sandstone containing randomly distributed spherical inclusions. Goodarzi et al. (2016) compared the applicability of MT, self-consistent and generalized self-consistent models in evaluating the effective moduli of shales, and they found that the MT and generalized self-consistent models could accurately predict the bulk modulus. Similarly, Zhao et al. (2018) compared the shale moduli predicted by the MT and self-consistent models based on the micro-indentation results, and reported that the MT model could better predict the anisotropy of shales.
Although the MT model has been introduced into the evaluation of the elastic properties of shales, the shapes and orientations of inclusions in the MT model have been assumed theoretically in previous studies. With the development of micro-computed tomography (micro-CT) and digital core technologies (Zhu et al., 2019), it is now possible to introduce the real information for inclusions in shales into the MT and other theoretical models. The micro-CT has been widely used to predict the effective moduli of composites (Sahimi & Tahmasebi, 2021), and the digital core technique is effectively used to describe the porous media in the analysis of rock conductivity, heat conduction, and seepage (Zhao et al., 2020). In terms of the elastic properties of rocks, is the most convenient tool for analysis of the digital core is the finite element method (FEM), which is accurate but time-consuming. However, there are few investigations on the application of digital core methods to effective medium models, including the MT model, to evaluate the elastic properties of shales.
For the combination of the MT model and the digital core technique, the key point is how to extract the geometric parameters of inclusions required in the Eshelby tensor. Drach et al. (2011Drach et al. ( , 2013 developed a method to extract individual pore inclusions from the micro-CT data and classified the irregular characteristic groups using the moment of inertia of any irregularly shaped pores. Park et al. (2020) obtained the distribution laws of pore shape and orientation in air-entrained (AE) cement by the pore's inertia matrix and predicted the effective moduli of anisotropic AE cement by the MT model. Considering that inclusions with small stiffness, such as pores and clay, in the shale follow similar laws, Giraud et al. (2007) proposed a theoretical model to describe the orientation distribution of inclusions. However, there are few studies investigating the application of the microscale (10 -5 *10 -6 m) inclusion's inertia matrix obtained by micro-CT to the evaluation of shale elastic properties.
As a consequence, the direct goal of this study is to obtain the inertia matrix of micron-scale inclusions in shale samples through micro-CT (X-ray computed tomography, X-CT), in order to analyze the shape and orientation laws of different kinds of inclusions. This information on the inclusions is combined with the MT model to evaluate the anisotropic properties of shale samples and the contribution of different kinds of inclusions.
The remainder of this article is organized as follows. In Sect. 2, the principle of the MT model for predicting the anisotropy of shale's elastic properties is introduced, in which the definition and introduction of Eulerian angles are formulated. In Sect. 3, through the reconstruction of CT images of the shale sample in the Longmaxi formation, information regarding the shape and orientation of inclusions with different components is obtained and analyzed. In Sect. 4, the elastic stiffness matrix of the shale sample is quantified by the MT model, and the effect of different kinds of inclusions on the shale's anisotropy is discussed. Moreover, a comparison of the results obtained by the MT model and FEM analysis is presented to verify the accuracy of the MT model in predicting the anisotropic properties of shales.

Fundamentals of the MT Model Considering the Orientation of Inclusions
We first provide an illustration of the MT model, which is one of the most widely used methods for predicting the effective properties of composite materials in various fields. The interaction among inclusions can be considered by distinguishing the far-field strain and the matrix strain in the MT model, so it is suitable for the analysis of representative volume element (RVE) including large volume fraction inclusions. For the shale with multiple and morphological inclusions, the prediction formulas of the effective modulus can be expressed as (Shen et al., 2013) where C is the effective elastic stiffness tensor, C 0 represents the elastic stiffness tensor of the matrix, and N is the number of components in the shale including the matrix and (N-1) kinds of inclusions. The subscript r is the component number, where 0 2982 Z. Wang et al. Pure Appl. Geophys. corresponds to the matrix. The parameter f r is the volume fraction of each component, satisfying P NÀ1 r¼0 f r ¼1. The quantity C r is the elastic stiffness tensor of the rth component, and S r is the Eshelby tensor of the rth component, satisfying S r = P r C 0 . The expression of the Hill tensor P r of general ellipsoidal inclusions is listed in Appendix A.
In the prediction of the effective modulus using the MT model, each inclusion can be regarded as an ellipsoid, and the three principal radii of an ellipsoidal inclusion satisfy a 1 [a 2 [a 3 , i.e. the principal inertia of the ellipsoid satisfies M 11 \ M 22 \ M 33 . The orientation of the ellipsoid with principal axes Ox 1 x 2 x 3 with respect to a set of axes O-XYZ fixed in the shale may be specified by the three Eulerian angles w, h, and u. The principal coordinate O-x 1 x 2 x 3 of any orientation in space can be obtained by three rotations of the axes O-XYZ and the specific process is schematized in Fig. 1.
Therefore, for any fourth-order tensor expressed as A 0 in the axes O-x 1 x 2 x 3 , its form A in the axes O-XYZ can be derived by the transformation equation, which is where T is a fourth order transformation tensor, and its corresponding 696 matrix form is The symbol R kl represents the component of a second-order tensor R in its 393 matrix form, and the corresponding matrix reads where R w , R h , and R u correspond to the second-order tensors of w, h and u to transform the vector from axes O-x 1 x 2 x 3 to axes O-XYZ, which are written as Then the matrix form of R is derived as When considering the orientation of inclusions, Eq. (1) can be rewritten as where n r represents the total number of inclusions of the rth component, the subscript i is the inclusion number, the subscript (r,i) represents the ith inclusion in the rth component, and T r,i , C r,i , and S r,i are the corresponding transformation tensor, elastic stiffness tensor and Eshelby tensor, respectively. On the premise that geometric parameters are obtained, the effective elastic stiffness tensor of the shale can be predicted by the MT model according to Eq. (7). The shape and orientation parameters of all inclusions will be determined in Sect. 3.

Digital Core of Longmaxi Formation Shale
As mentioned earlier, this study focuses on the influence of micron-scale inclusions on anisotropy in shales. Therefore, considering that the sample size is in the range of 10 -2 *10 -3 m and that the voxel size of the image is in the range of 10 -5 *10 -6 m obtained by the X-CT method, the following two factors affecting the anisotropy in the shale sample are ignored: (1) The rock physics method such as the MT model usually focuses on the shape and orientation of inclusions, and the mechanical stratigraphy of shales cannot be considered. Therefore, the stratigraphy caused by microscale layers which can be observed in the core with a length being 50 mm and a diameter being 25 mm in the work of Ramos et al. (2019), will be ignored in the following analysis. As a result, a smaller core with a diameter and a height of less than 5 mm is selected to describe the geometric information of inclusions more accurately, and the impact of microscale layers in the model can be ignored.
(2) Although recent studies (Vernik & Milovac, 2011;Sayers & den Boer, 2018) show the importance of accounting for the elastic anisotropy of clay particles, this property of clay cannot be effectively combined with the CT image. Locally regularly arranged clay platelets observed by scanning electron microscopy have a length of 700 Å (1 Å = 10 -10 m) (Sayers & den Boer, 2019), and they are difficult to identify at the micro-CT scale. Asaka et al. (2021) simplified this kind of clay domain into anisotropic microscopic inclusions, but the information for each platelet was still difficult to distinguish in the micro-CT images. Therefore, the anisotropic effect of clay is not considered in this analysis.
Next, the core sample of Longmaxi-formation shale with a diameter of about 4 mm is used to generate the so-called digital core. The core is processed via the tomograph in the Shanghai Synchrotron Radiation CT Facility of the Shanghai Institute of Applied Physics, Chinese Academy of Sciences. Based on the different absorption capacities of different substances to light, a set of 1397 highresolution core slice images are obtained, which represent different components by different gray levels. Figure 2 shows the first and last horizontal slice images of this core, and each image contains 200091400 pixels. Different components of the shale can be clearly distinguished according to different gray levels of the image, such as pores and TOC (the darkest part, shown in yellow circles), pyrite (shown in white and red circles), clay, quartz, feldspar and other minerals (gray-black to gray-white, shown in blue and green boxes). It can be seen that the Longmaxi shale demonstrates very strong heterogeneity, which has a major impact on the elastic properties and anisotropy of the rock. Specifically, the directional arrangement of pores and clays in the shale may lead to obvious anisotropy (Sayers, 1994).
For the subsequent analysis of the shale's effective modulus, it is necessary to batch-preprocess the images and divide the core into six components according to the gray level, as follows: (1) pore, (2) TOC, (3) clay, (4) quartz, (5) feldspar, calcite, and (6) pyrite, siderite (Jian et al., 2020). The preprocessed image is divided into six relatively fixed gray values, and the first and last preprocessed images are shown in Fig. 3. The different components are clearly distinguished, and the difficulty of 3D reconstruction is reduced.
Then, all the slice images are imported for 3D reconstruction. Based on the results of computed tomography, we determine that the sizes of a pixel are 2.841 lm, 2.841 lm, and 1.665 lm along the X, Y, and Z-axes, respectively. The reconstructed core of the scanning section is shown in Fig. 4a. However,  the cylindrical core cannot meet the boundary symmetry conditions required for the RVE analysis, so the cuboid containing 900990091200 pixels is intercepted from the core center and denoted as RVE total , as shown in Fig. 4b.

Geometric Information of Each Component
Binarization is carried out to obtain the corresponding components in RVE total by setting different thresholds, and the binarization results of pores are shown in Fig. 5a. In fact, the RVE total in Fig. 5a contains more than 10 5 individual pore inclusions, which is a huge task for mesh generation. In addition, Figure 4 Three-dimensional reconstruction of shale digital core from Longmaxi formation: a sample of scanning section, b RVE total for analysis. we take a single pore inclusion from the RVE total , which is displayed in Fig. 5b-d. It can be seen that the pores obtained only by threshold segmentation have irregular and sharp boundaries, as shown in Fig. 5c, d. These defects and bulges have little effect on the anisotropy of the shale caused by pores, but they will greatly affect the calculation efficiency of the MT model and FEM. Therefore, we have carried out the digital image processing on the core, and the process includes removing small holes and closing, opening, and removing small spots. In brief, these operations increase the rationality of the MT model and reduce the calculation amount in the MT model and FEM. The processed results of pores in the RVE are demonstrated in Fig. 5e, h. It can be seen from Fig. 5e that the small-volume pores are deleted and the number of pores decreases significantly. On the premise of maintaining the overall characteristics of pores, the opening and closing operation can reduce the defects and bulges of pores, which can be obtained by comparing Fig. 5c, d and g, h.
By adjusting the relevant parameters in the image processing, the volume fraction of each component before and after processing is maintained at a similar value. The volume fraction of each component after processing is shown in Table 1. It can be seen from Table 1 that the component with the largest volume fraction in the RVE is quartz, which exceeds 50%, so it should be used as the matrix phase when using the MT model. Accordingly, other components are considered as inclusions. In addition, the bulk and shear moduli of each component are listed in Table 1. In particular, the anisotropy caused by the clay lamination usually requires the observation of largerscale ([10 mm) cores (Ramos et al., 2019). For the convenience of analysis, all types of inclusions, including clay, are considered as isotropic material.

Ellipsoidal Approximation of Inclusions
Inclusions of different components in the shale usually have irregular shapes. For example, the pores in Fig. 5h show a shape similar to a flat ellipsoid. However, the ellipsoid is usually adopted as a representative element to characterize the contribution of inclusions to the mechanical properties of the RVE in mesomechanics (Drach et al., 2011), and so is the MT model. In a previous study, Drach et al. (2013) used the moment of inertia of the pore to classify the irregular pores of carbon/carbon composites, and it was verified that the inertia matrix of the inclusion can significantly affect the anisotropy of the RVE. As a consequence, the principal radius and principal direction of the approximate ellipsoidal inclusion used in the MT model are obtained through the inertia matrix of the irregular inclusion.
The inertia matrix of an irregularly shaped inclusion in the coordinate system O-XYZ is where I is the inertia tensor of the irregularly shaped inclusion, [I] is the corresponding inertia matrix, I XX , I YY , and I ZZ are the moments of inertia of the inclusion along the X, Y, and Z directions, respectively, and I XY , I YZ , and I ZX are the inertia products of the inclusion in the X-Y, Y-Z, and Z-X planes, respectively, m is the mass of the inclusion, and M 2X , M 2Y , M 2Z , M 2XY , M 2YZ , M 2ZX are detailed in Appendix B.
After establishing the inertia matrix of the inclusion, the principal moment of inertia and the directions of the principal axes can be obtained. Then, the principal radii of the approximate ellipsoidal inclusions can be obtained according to Eq. (8), and the Eulerian angles corresponding to the principal axes of each inclusion can be obtained using the dcm2angle function in MATLAB via the relations where I 1 , I 2 and I 3 are the first, second, anIn order to evaluate the contributiond third principal moments of inertia of the inclusion, respectively.

Statistics of Shape and Orientation of Inclusions
Previous results show that the shape and orientation of inclusions will significantly affect the elastic properties and anisotropy of the shale. Generally, when the effective medium model is used for calculation, the inclusion is considered a rotating ellipsoid, that is, a 1 = a 2 or a 2 = a 3 . The aspect ratio a is used to describe the shape of the rotating ellipsoidal inclusion. Since the size relationship between the three principal radii is defined in this study, there are two cases for the value of a: (1) when a 1 = a 2 , a = a 3 /a 1 \1; (2) when a 2 = a 3 , a = a 1 /a 3 [1.
Previous studies have generally discussed the influence of the aspect ratio of the inclusion on the elastic properties of the rock, but the selection of the aspect ratio of the inclusion is typically assumed in advance. Therefore, extracting the shape information for components in the shale by the digital core technology and making reasonable ellipsoid approximations are of clear significance for the selection of the aspect ratio.
The average principal radii of all inclusions in different components are counted to determine the representative shapes of inclusions in different components, which can be obtained by where a 1;r , a 2;r and a 3;r are the first, second, and third average principal radii of the rth component, respectively, a 1,i , a 2,i and a 3,i are the first, second, and third principal radii of the ith inclusion, V i is the volume of the ith inclusion, and V r is the total volume of the rth component.
Since the Hill tensor of an ellipsoidal inclusion is only related to the ratio between the principal radii and is independent of its size, the two ratios (a 2 =a 1 and a 3 =a 1 ) of each inclusion are the shape parameters affecting the calculation results of the MT model. The values a 2 =a 1 and a 3 =a 1 of five components are shown in Fig. 6. The two ratios of inclusions can be divided into four cases: (1) when a 2 =a 1 , a 3 =a 1 ! 1, the inclusions tend to take the spherical shape; (2) when a 2 =a 1 ! 1, a 3 =a 1 ( 1, the inclusions tend to take the disk shape; (3) when a 2 =a 1 ,a 3 =a 1 ( 1, the inclusions tend to take the needle shape; (4) when 0 ( a 2 =a 1 ( a 3 =a 1 ( 1, the inclusions are general ellipsoids without obvious characteristics. It can be seen from Fig. 6 that in the shale samples, the clay, feldspar, and pyrite tend to be needle-shaped, while the pores and TOC show general ellipsoids without obvious characteristics. Similarly, the Eulerian angles of inclusions with different components are statistically analyzed, where we take the pores as an example. The Eulerian angle distributions of the pore inclusions are shown in Fig. 7. It can be seen that there is a large pore with volume of 1.315910 8 lm 3 in the digital core, which accounts for 30.661% of the total pore volume. If the influence of the maximum pore on the law is ignored, it can be seen from the orientation distribution that w, h and u tend to the maximum at 0 degrees. From the orientation distribution of the three Eulerian angles, w has obvious symmetry, while the distributions of h and u are not regular. Interestingly, the distribution law of w is similar to the theoretical model given by Giraud et al. (2007) to describe the orientation distribution function of inclusions in rocks, which can be expressed as where W(w) is the orientation distribution function of w, and r is the parameter for the degree of the preferred alignment. There are two typical distribution orientation cases: (1) when r ! 0, the distribution of w of inclusions is uniform; (2) when r ! 1, the inclusions of the rth component have the same w. In other words, the effect of inclusions on the anisotropy of materials increases with the increase in r. Therefore, the parameter r can be used to describe the degree of anisotropy of orientation distribution. The orientation distribution functions of the five kinds of components can be fitted by Eq. (8), and the fitting results are shown in Fig. 8. The value of r of the five component inclusions in the sample core are TOC, clay, pore, pyrite, and feldspar from large to small. This means that, when considering only the orientation distribution, the TOC, clay, and pores have more obvious effects on the anisotropy of the sample core, followed by the pyrite, and the feldspar has no contribution to the anisotropy of the sample core because its r = 0.

Effective Moduli of the Shale Sample Predicted by the MT Model
In this section, the anisotropic effective moduli of the shale sample (RVE total ) are predicted using Eq. (7) and the inclusion information in Sect. 3.3. The volume fraction of each component inclusion (f r ) is recorded in Table 1, and the volume fraction of matrix/quartz (f 0 ) is f 0 ¼ 1 À P NÀ1 r¼1 f r . When analyzing the effect of each component on the anisotropy of RVE total , the RVE total is considered as a two-phase material, and thus f 0 =1f r .
In order to evaluate the contribution of different mineral components to anisotropy, five kinds of twophase samples with a single type of inclusion and quartz matrix are predicted by the MT model. The components of the effective elastic stiffness matrix of the two-phase RVE inclusion-matrix predicted by the MT model are shown in Table 2. From the prediction results of C 11 , C 22 and C 33 , it can be seen that the pores, TOC and clay contribute more to the anisotropy of the shale sample than the feldspar and pyrite. Moreover, the effects of the pores, TOC and clay on the anisotropy are the same, that is, C 11 [C 22 [C 33 and C 44 \C 55 \C 66 , and their differences (C 11 -C 33 ) are 2.466 GPa, 2.860 GPa, and 2.844GPa, respectively. However, the results of the feldspar and pyrite are C 22 [C 11 [C 33 , C 44 \C 66 \C 55 and C 11 [C 33 [C 22 , C 44 \C 66 \C 55 , and their differences between C 11 and C 33 are only 0.014 GPa and 0.069 GPa, respectively. Different components with different moduli will take different distribution forms under the action of in situ stress, which directly reflects their contributions to the anisotropy of the shale. For example, the moduli of pores, TOC, and clay are generally smaller than those of the quartz/matrix, so their contributions are similar. The moduli of the feldspar are similar to that of the quartz/matrix, and the moduli of the pyrite are significantly larger than that of the quartz/matrix, so their contributions are different.
The components of the elastic stiffness matrix of RVE total with five kinds of inclusions and quartz matrix predicted by the MT model are shown in Table 2. The anisotropy laws of this RVE total are  C 11 [C 22 [C 33 and C 44 \C 55 \C 66 , which are mainly affected by the internal pores, TOC, and clay. The difference between C 11 and C 33 is 3.299 GPa, and the Thomsen anisotropy coefficients (Thomsen, 1986) e and c are 0.0259 and 0.0175, respectively, which shows that the anisotropy of this sample is still relatively small. However, as a preliminary attempt, it is important to analyze the rationality of the MT model in the prediction of the shale's effective moduli.

Validation by FEM
Although the influence of the shape and orientation of inclusions on the effective modulus of materials is considered in the MT model, there are still two factors which may cause prediction errors. The first is that the shape of the inclusion is simplified as an ellipsoid, which may cause variation in the stress distribution of the matrix around the inclusion and affect the prediction results. The other is that, in the effective medium models including the MT model, the location of inclusions in the RVE cannot be considered, and the inclusions with an obvious regular arrangement may also affect the anisotropy of materials. As a result, the accuracy of the MT model in predicting the effective moduli and anisotropy of RVE is verified by FEM using Abaqus software.
It should be noted that the geometry of the RVE obtained in Sect. 3.2 needs to be meshed, which means that the calculation using FEM is much more complex than that using the MT model. There are tens of thousands of inclusions in the RVE total , and thus it is unrealistic to analyze directly by FEM. Therefore, the RVE total is divided into 36 RVEs containing 30093009300 pixels and numbered as RVE 1 to RVE 36 for subsequent analyses.
In order to facilitate comparison, the analysis object is taken as the quartz-pore material. The moduli of the quartz and saturated pores are shown in Table 1. The deformation of the generated quartz/matrix mesh is consistent with that of the pore/inclusion mesh, that is, the boundaries of the two kinds of meshes are the consistent pair. In addition, both the quartz/matrix and pores/inclusions are reconstructed by tetrahedral grids. By adding the Micromechanics Plug-Ins modulus into Abaqus, the following Lagrangian strains (Moon et al., 2015) are applied to the RVE: According to the stress field results, the average stress is calculated, and the effective stiffness matrix is inversely calculated.
To verify the rationality of ellipsoidal approximation without considering the influence caused by the second factor mentioned above, an RVE sm containing only three pores is intercepted in the RVE 1 for analysis. The three models on the RVE sm with different grid sizes are shown in Fig. 9. Table 3 lists the components of the elastic stiffness matrix predicted by FEM and the MT model. It can be seen that the prediction results of FEM and the MT model are very similar, with the relations C 33 [C 11 [C 22 and C 66 \C 44 , C 55 . The difference in the prediction results between the MT model and FEM decreases with the increase in grid density. Fig. 10 shows the Mises stress of the RVE sm under e 1 (tension along X-axis), where half of the length of the RVE sm along the X-axis is truncated to observe the stress distribution around the pores. It can be seen from Fig. 10 that the stress in the pore is obviously less than that in the matrix, and there is a stress disturbance caused by the pores in the matrix around the pore. It can be seen from the local figures that this stress disturbance becomes obvious with the increase in the grid density and tends to be more practical. With the improvement of the grid quality, the prediction results of effective moduli become smaller, which shows that the variation range of the low stress area of the matrix is higher than that of the high stress area. In the MT model, this Vol. 180, (2023) Prediction of the Anisotropic Effective Moduli of Shales based on the Mori-Tanaka Model 2991 disturbance is averaged to the whole matrix region, so the effective moduli predicted by the MT model are more suitable for the case of inclusions with large volume fraction than the sparse model or the Kuster-Toksöz model. By comparing the results between the MT model and FEM in Table 3, it can be seen that the averaging treatment of stress disturbance in the MT model is reasonable in predicting the effective moduli. Further, the three kinds of RVE 1 with different grid densities are analyzed by FEM, of which the grid division is shown in Fig. 11a, b, and c. There are 8599, 49,565, and 585,803 elements representing the matrix and 3855, 9686, and 37,346 elements representing the pores in the three models, respectively. Figure 11d, e and f show the Mises stress of RVE 1 under e 1 (tension along the X-axis), where we cut half of RVE 1 to observe the stress distribution. Similar to the analysis results for RVE sm , the stress disturbance area of the matrix around the pores gradually increases with the densification of the grid, and the effective moduli of RVE 1 generally decrease slightly with the densification of the grid, which can be seen from Table 4. It can be seen that the anisotropy laws of RVE 1 predicted by FEM are in agreement with those predicted by the MT model, i.e., C 11 [ C 22 [ C 33 and C 44 \C 55 \C 66 . Therefore, in the RVE 1 , the second factor mentioned above has little effect on the anisotropy.
However, comparing the results between Tables 3 and 4, the difference between the results of RVE 1 predicted by FEM and the MT model is larger than that of the RVE sm . The main reason is that the volume ratio between the pore and element in RVE 1 is smaller than that in RVE sm , and there are relatively large errors in the matrix stress disturbance near the pores. However, further refinement of the grid will greatly affect the computational efficiency of the FEM and thus impose high memory requirements for the computer. Mesh division of the RVE sm : a low mesh quality, b medium mesh quality, and c high mesh quality. Table 3 Comparison of the effective stiffness matrix of the RVE sm predicted by the MT model and FEM (GPa)

Effect of the Volume Fraction of Inclusions
In what follows, the effective moduli of the RVE 1 to RVE 36 containing the quartz/matrix and pores/ inclusions are predicted by the MT model and FEM for analyzing the effect of the inclusion's volume fraction. One of the reasons that the MT model can be widely used is that it can be applied to the case dealing with a large volume fraction of inclusions (Raju et al., 2018). The effective moduli of the RVE 1 to RVE 36 predicted by the MT model and FEM are shown in Fig. 12, and the average values are shown in Table 5.
The prediction results of the two methods are close, and the laws of anisotropy are the same: C 11 [ C 22 [ C 33 and C 44 \ C 55 \ C 66 . Moreover, the effective moduli predicted by the two methods decrease linearly with the increase in the volume fraction of inclusions. In fact, it can be seen from previous studies that when the shape of the inclusions is not a disk shape (a 2 =a 1 ! 1, a 3 =a 1 ( 1), the relationship between the effective moduli predicted by the MT model and the volume fraction of the inclusions basically changes linearly (Wang et al., 2021). Interestingly, the result predicted by the MT model is more discrete than that of FEM, especially when the inclusion volume fraction is large, which indicates that the MT model may overestimate the influence of inclusion's orientation on the anisotropy of the RVE. However, when using the digital core for modeling in FEM, we have to calculate a large number of elements, and efficiency is relatively low. Therefore, the combination of the MT model and digital core proves a practical and efficient method for predicting the anisotropy of shales and other kinds of rock.

Conclusions
In the present study, we establish and verify a prediction method for anisotropic effective moduli of Longmaxi formation shales based on the micron digital core, the MT model and the finite element method (FEM). We confirm that the moment of inertia of inclusions in shales can be used to obtain the principal radius and Eulerian angles of its representative ellipsoid. By analyzing the shale Figure 10 Mises stress of the RVE sm under tension along the X-axis Vol. 180, (2023) Prediction of the Anisotropic Effective Moduli of Shales based on the Mori-Tanaka Model 2993 components extracted by the micron-scale digital core, it is found that the orientation distribution of different components in the shale is similar to the theoretical function given by Giraud. At the micron scale, inclusions of pores, TOC, and clay show a more anisotropic orientation distribution, which leads to these three kinds of inclusions contributing greatly to the anisotropy of the shale. In addition, the geometry of inclusions in micron shale cores obtained by micro-CT imaging cannot be accurately simplified into rotating ellipsoids typically used in traditional rock physics models. The anisotropy of elastic properties of the shale sample was predicted by the MT model, and the five groups of two-phase material models, where it was determined that pores, TOC, and clay contribute greatly to the anisotropy of shale, while feldspar and pyrite have almost no effect on the anisotropy of the shale samples. The effective elastic stiffness matrixes of multiple groups of quartz-pore RVE obtained by the digital Figure 11 Mesh division and Mises stress of the RVE, where a, b, and c are the low, medium, and high mesh quality, and d, e and f are the corresponding Mises stress under tension along the X-axis. core are predicted by the MT model and FEM, and the anisotropy of RVE predicted by the two methods is consistent. The mesh accuracy of FEM has a great influence on the prediction results. The effective moduli change linearly with the inclusion volume fraction, and the discrete prediction result of MT is greater than that of FEM. These findings demonstrate potential for application in fields such as rock mechanics, civil engineering, and oil exploitation.
Author contributions For author contributions, we declare that ZW, as the first author of this manuscript, is responsible for the theoretical and numerical calculations for all the models in the paper. GC has helped with the model and results analyses. JL has helped to modify the paper and all the result  analyses. LF has discussed the details of the paper with the other authors.

Funding
This study was funded by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA14010303), the National Natural Science Foundation of China (11972375), and the Science and Technology Project in the West Coast New Area of Qingdao (2020-81).

Data availability
The supporting information involved in this study is available within this article.

Declarations
Conflict of interest The authors declare that they have no conflicts of interest to report regarding the content of this article.
Ethics approval and consent to participate Written informed consent for publication of this paper was obtained from the China University of Petroleum and all authors.