Effective moduli of rocks predicted by the Kuster–Toksöz and Mori–Tanaka models

Natural rocks are polymineral composites with complex microstructures. Such strong heterogeneities significantly affect the estimation of effective moduli by some theoretical models. First, we have compared the effective moduli of isotropic rocks predicted by the Kuster–Toksöz (KT) model and the Mori–Tanaka (MT) model. The widely used KT model only has finite precision in many cases because of its assumption that is restricted to the first-order scattering approximation. However, the MT model based on the Eshelby tensor in mesomechanics has the advantage of predicting effective moduli of rocks, especially when the volume fraction of embedded inclusions is sufficiently large. In addition, the MT model can be used to predict the effective modulus of anisotropic rocks, but the KT model cannot. For a certain kind of shale or tight sandstones, which are viewed as isotropic composites, both the models work well. For the medium containing spherical pores, both the models produce the same results, whereas for ellipsoidal pores the MT model is more accurate than the KT model, validated by the finite element simulations. In what follows, the applicable ranges of simplified formulas for pores with needle, coin and disk shapes, widely used in engineering, are quantitatively given based on the comparison with the results according to the reduced ellipsoidal formulas of the MT and KT models. These findings provide a comprehensive understanding of the two models in calculating the effective modulus of rocks, which are beneficial to such areas as petroleum exploration and exploitation, civil engineering, and geophysics.


Introduction
Unconventional reservoir rocks, such as tight sandstones and shales, are characteristic of multi-mineral compositions and complex microstructures. The resulting strong heterogeneities and complicated mechanical properties (Prasad et al. 2009;Guo et al. 2013;Jian et al. 2020) challenge the evaluation of effective moduli. In particular, the directionally distributed microcracks in tight sandstones or shales from deep reservoirs due to deviatoric stress fields often make reservoir rocks transversely isotropic, which significantly increases the difficulty in estimating effective elastic moduli by conventional theoretical models (Shapiro & Kaselow 2005;Gui et al. 2018;Zhang et al. 2020a).
There are two well-known theoretical models that have been widely used to evaluate the effective bulk and shear moduli of porous rocks. The Kuster-Toksöz (KT) model (Kuster & Toksöz 1974) is formulated based on first-order wave scattering approximation for the effect of elasticity. The Mori-Tanaka (MT) model (Mori & Tanaka 1974) is an equivalent method based on the Eshelby inclusion theory in mesomechanics, where the interaction among inclusions has been fully considered. Both models consider the influence of inclusion shape on effective modulus, the pores with a special shape such as cracks or air voids also played important roles (Zhang et al. 2020b). In the KT model the effect of fluid flows in rocks is ignored and, therefore, it only applies to the case of high frequency (Mavko et al. 1998). In addition, in the KT model it is assumed that the inclusions inside a rock are normally sparse and randomly oriented, which, however, is not satisfied in practice (Xu & White 1995). To handle these problems with the KT model, the Gassmann equation (Gassmann 1951) is usually used to transform the results to the case of low frequency. Moreover, the KT model is usually used in combination with the differential effective medium method so as to overcome the drawback that the model only applies to a case where the inclusions have a small volume content (Xu & White 1995;Liu & Sun 2015;Mollajan et al. 2018;Zhang et al. 2018;Zhang et al. 2020c). To meet the condition of sparse inclusions in the KT model, the inclusions are added to the rock matrix by iteratively replacing the background medium, which may greatly reduce the computational efficiency. It should be noted that the effect of directionally distributed inclusions on the elastic properties of rocks has not been taken into account in the KT model. In general, inclusions with small modulus in reservoirs, such as kerogens (Vernik & Landis 1996;Ma & Chen 2014;Liu et al. 2017) and microcracks (Sayers 1994), tend to have an oriented distribution due to deviatoric stress.
To consider the anisotropy of rocks, many models have been proposed, but they are still not perfectly developed. Among these models, the Backus model (Vernik & Nur 1992) only aims at the anisotropy caused by the multilayer arrangement of the reservoir. Another is the Hudson model (Hudson 1981(Hudson , 1986, where the average wavefield scattering theory is used to explore anisotropic rocks, but only thin coin-shaped inclusions are considered. Based on the Eshelby effective inclusion theory (Eshelby 1957), Cheng (1993) calculated the effective modulus of a transversely isotropic rock with cracks (the Eshelby-Cheng model), which has been widely adopted in the field of rock physics (e.g. Liu et al. 2006;Zhao et al. 2016;Figueiredo et al. 2018).
In practice, the theory of mesomechanics has been extensively studied for composite materials (Berryman 1995). Although the rock can be regarded as a composite material with complex structures and compositions, the mesomechanics has not been widely extended. One of the effective and efficient methods in mesomechanics is the MT model, where the influence of inclusions on the average stress of matrix is considered, and thus it can better describe the interaction among inclusions. Previous results show that the MT model has been applied in many fields (e.g. Shen et al. 2013;Wang et al. 2017;Park et al. 2020), but there are very few reports on the application to rock physics. Zimmerman (1991) compares the MT model with other two classical models, but the analysis is only limited to spherical inclusions. Deng et al. (2015a) use the MT model to calculate the effective moduli of a dry sandstone containing randomly distributed spherical inclusions. Therefore, in the present work we make a comprehensive study on the effective moduli based on the MT and KT models, with an attempt to clarify their applicable ranges.
The article is organised as follows. In Section 2, we outline the fundamentals of the KT and MT models that are applied to the evaluation of effective moduli for Longmaxi formation shales. The results are compared with the experimental results. To illustrate the difference between these two models, in Section 3, we model a two-phase sandstone containing spherical and ellipsoidal pores, respectively. Both of the results calculated by the KT and MT models are validated by finite element (FE) simulations. In Section 4, we extend the application to rocks containing coins, disks and needles, with an attempt to quantitatively explore the applicable ranges of the MT and KT models for convenience of engineering usage.

Fundamentals of the MT and KT models
According to the theory of mesomechanics for evaluation of effective moduli of a composite material, we select the phase with the largest volume fraction as the matrix. For example, in a simple model for a porous sandstone we would consider solid phases as the matrix with pores as inclusions. Based on the fundamental spirit of mesomechanics, a representative volume element (RVE) can be defined to be big enough relative to the microstructures inside the volume, but small enough compared to the macroscopic material. Such an RVE has its stress or displacement boundaries, and it provides an efficient way to calculate the effective moduli of composite materials. As a result, both the mesoscopic and macroscopic RVEs have the same response to external loads. That is, the average of loading-induced strains from individual phases of multi-mineral composites is equivalent to the effective strain by taking the multi-mineral composites as a homogeneous RVE .
We assume the material of an RVE satisfies linear elasticity with small deformations. Based on the generalised Hooke's law, the constitutive relationship defined for this RVE can be expressed as where i , M i and i are the stress tensor, elastic stiffness tensor and strain tensor of the i-th phase, respectively, in which i = 0 represents the matrix phase. For any physical quantity 540 , its average valuēin the RVE is defined as where V is the volume of the RVE. Consequently, the effective moduli of the RVE with multiple components can be calculated bȳ=M wherēand̄represent the average stress and strain of the RVE, respectively, andM is the effective elastic stiffness tensor to be determined. Among the various models for the effective moduli of composite materials, the interaction between the matrix and inclusions can be taken into account in the MT model, and this fact facilitates the analysis for the material with inclusions whose volume fraction is larger. For the RVE with n components that can be expressed as the matrix plus n-1 kinds of inclusions, its effective moduli can be computed by (Weng 1984 where c i is the volume fraction of each phase, satisfying ∑ n−1 i=0 c i = 1. The strain concentration tensor A i is where M 0 is the elastic stiffness tensor of the matrix, and the stress concentration tensor B i is where I represents the identity tensor, and S i is the Eshelby tensor of the i-th phase inclusion (Zhao & Weng 1990;Shen et al. 2013). In Appendix A we list different expressions of S i for the i-th phase inclusion with specific shapes. As one typical example for the two-phase composite material (n = 2), substituting equations (5) and (6) into (4), one can simplify the expression ofM as follows Due to the symmetry of stress, strain, and elastic stiffness, the second-order tensor is a column vector with six elements, whereas the fourth-order tensor is a 6 × 6 matrix. The corresponding three tensors can be formulated as One can see that the fourth-order tensor is actually a symmetric matrix because of the Voigt symmetry. The corresponding equivalent tensor I is also a 6 × 6 identity matrix. Solving inverse fourth-order tensors for effective moduli is detailed in Appendix B, which can be implemented by the self-developed code via MATLAB.
The KT model has been widely used to predict the effective moduli of rocks. In the model it is assumed that the inclusions must be randomly distributed in the matrix to make the equivalent material be isotropic. The KT model is generally expressed as where 0 = (9K 0 + 8G 0 )G 0 /[6(K 0 + 2G 0 )],K andḠ are respectively the effective bulk and shear moduli, K 0 and G 0 are, respectively, those of the matrix, and K i and G i are the bulk and shear moduli of the i-th inclusion, respectively, with P 0,i and Q 0,i its shape coefficients given in Appendix C. As a typical example, we use the experimental data of a shale by Deng et al. (2015b) to compare the results on effective moduli from the KT in equation (11), the MT in equation (4) and the Hill-average method (Mavko et al. 1998). As a simple and convenient method, the Hill-average strategy averages the Voigt upper and Reuss lower bounds. The mineral components with relevant elastic parameters are shown in Table 1  The effective moduli of the shale predicted by these three methods are shown in Table 2 for TOC-and water-saturated pores, respectively. One can see that the predicted shale moduli by the three methods are close to the experimental results for TOC-saturated pores. For water-saturated pores, the resulting shale moduli by the KT and MT models are similar to those of the shale with TOC-saturated pores, but the Hill-average shear modulus is obviously smaller because of the nearly zero modulus of water included. For both the TOC-and water-saturated pores, the KT-predicted shear moduli of the shale are slightly bigger than the experimental results because of its limited accuracy, while the MT model gives more accurate predictions. In addition, both the KT-and MT-predicted bulk moduli of the shale are slightly bigger than the experimental results possibly because not all inclusions are spherical in nature. This indicates that the shape of inclusions significantly affects both of the values ofK andḠ based on the KT and MT models.

Comparison of the KT and MT models
We compare the KT and MT models for estimation of effective moduli of a two-phase sandstone. Based on the previous work ( Jian et al. 2020), the bulk and shear moduli of the matrix (quartz) are selected as K 0 = 37 and G 0 = 44 GPa, respectively. We assume the sandstone contains randomly dis-tributed pores (as inclusions) that are saturated by water with K 1 = 4.07 and G 1 = 0.001 GPa. Both spherical and ellipsoidal pores are used to validate the KT and MT models. For the former, the resulting RVEs are distributed homogeneously to make the sandstone be isotropic. For the latter, the resulting RVEs are aligned to yield a directional distribution of inclusions. In the Cartesian coordinate system O-xyz, the ellipsoid is given as where the aspect ratio = a 1 /a 2 = a 1 /a 3 with a 1 , a 2 and a 3 the principal radii in the x, y and z directions, respectively. The rotation axis of the ellipsoid is along the x axis. Especially, when = 1, the ellipsoidal inclusion is reduced to the spherical one.
The MT and KT models are compared to calculate the effective bulk and shear moduli for a two-phase sandstone with spherical and ellipsoidal water-saturated pores, which are shown in figure 1. One can see that the values ofK andḠ are the same for spherical inclusions. As indicated by Mavko et al. (1998), both the KT and MT models for isotropic spherical inclusions reduce to the Hashin-Strikman upper bound,K For aligned ellipsoidal inclusions, the MT model gives more accurate predictions ofK andḠ than the KT model, but both tend to be similar with increasing aspect ratios. As is well known, the effective moduli of spherical inclusions (SI) are the upper bounds of ellipsoidal inclusions (EI) (Mavko et al. 1998), which can be verified by equation (13). From figure 1c and d, one can see that the results based on the KT model have large errors at low aspect ratios, even yielding an unreasonable value smaller than K 1 = 4.07 GPa, while this case does not occur in the MT model.
When the effective modulus of the sandstone amounts to the averaging modulus of all inclusions (e.g. K 1 ), the whole sandstone reduces to pores, with the volume fraction of pores being 100% in theory. However, in the KT model, the maximum volume fraction of pores only reaches a limited value smaller than 100%, which is defined as the critical volume fraction c 1, cr . Figure 2 shows the dependence of the critical volume fraction c 1, cr of inclusions on the aspect ratio for both of the two models. One can see that for = 1 (spheroids), the KT model reaches its maximum value (the level of the MT model) where the shadow area is defined as the applicable range of the KT model. In addition, the 542 Figure 1. Comparison of the KT and MT models to estimate the effective bulk and shear moduli for a two-phase sandstone with spherical and ellipsoidal water-saturated pores, respectively.  KT curve in figure 2 is not monotonic, with the watershed at = 1. The applicable range of < 1 (flat ellipsoids) is much smaller than that of > 1 (elongated ellipsoids). Previous work (Kuster & Toksöz 1974;Zhang et al. 2018) shows that the KT model only applies when c 1 / ≪ 1, but without the critical value of c 1 / quantitatively detailed. The condi-tion of c 1 / ≪ 1 becomes quite different between < 1 and > 1.

Validation with experimental data by FE simulation
In this section, we apply the KT and MT models to the experimental data of a tight sandstone (Yan et al. 2011  bulk and shear moduli of the matrix are 39.89 and 30.01 GPa, respectively. The bulk moduli of dry and saturated inclusions are 0.001 and 2.505 GPa, respectively. Considering some soft or flexible pores and cracks existing in the rock, we take the aspect ratios as 1, 0.1 and 0.01. Figure 3 compares the experimental data and the effective moduli predicted by the MT (part a) and KT (part b) models for pores with different aspect ratios. One can see the MT predictions of effective moduli with various aspect ratios basically cover the data distribution of experiments, which is much better than the KT predictions, especially for those with small aspect ratios. It is difficult to validate theoretical models directly using the experimental data of natural rocks because of complex inclusions with varying sizes and shapes. As an efficient alternative, numerical experiments can facilitate the examination of theoretical models. Herein, we use the FE software COM-SOL to model porous RVEs with complex inclusions in sizes, shapes, numbers and arrangements. In the COMSOL simulation, periodic boundary conditions are applied to the RVE, and six sets of average strains are used for loading, which are expressed as follows The corresponding average stress matrices can be obtained by FE simulations. The effective stiffness matrix [M] can be evaluated based on equation (3), with its inverse matrix expressed as 544 Figure 4. Schematics of a RVE with three types of inclusions of arbitrary orientation (a), in-plane arbitrary orientation (b) and single orientation (c).
whereĒ 1 ,Ē 2 andĒ 3 are the effective elastic moduli along the x, y and z directions, respectively, andḠ 23 ,Ḡ 31 andḠ 12 are the effective shear moduli in the y-z, z-x and x-y planes, respectively. The Poisson's ratiov ij takes the negative value of the ratio of the j-direction contraction to the i-direction extension when the normal stress is applied only in the i-direction. The symbols i, j = 1, 2, 3 represent the directions of x, y and z, respectively. Figure 4 shows a RVE schematics with three types of inclusion of arbitrary orientation, in-plane (the x-y plane) arbitrary orientation and single orientation along the x axis (Vernik & Landis 1996;Ma & Chen 2014;Liu et al. 2017), yielding increasing anisotropic characteristics. These inclusions are selected to be ellipsoidal configurations with = 0.5. In figure 5 we compare the effective elastic moduli by the KT and MT models and the FE numerical method for randomly distributed SI, EI arbitrarily oriented in space, EI arbitrarily oriented in plane and EI single oriented along a coordinate axis, respectively. For randomly distributed SI, the KT and MT models give the same predictions as expected from equation (13). In the FE simulation for this case, the RVE consists of one inclusion, eight inclusions and 64 inclusions, respectively. One can see that the KT and MT predictions agree well with the FE simulations for randomly distributed SI (figure 5a).
For EI arbitrarily oriented in space, the resulting RVE can be regarded as an isotropic material. The FE simulations are conducted for an RVE including 27 and 64 EI, respectively. The resulting matrix and inclusions are all divided into tetra-hedral meshes. For instance, there are 364 182 tetrahedral elements in the RVE with 64 ellipsoid inclusions. The convergence of FE simulations has been verified by setting meshes with different precisions. As shown in figure 5b, both the analytical KT and MT predictions match the FE simulations well. Compared to the KT prediction, the MT prediction is closer to the FE simulation, but the error between them is not so much, implying that both the KT and MT models can be applicable for isotropic materials. We should stress that the applicability is limited to the EI of = 0.5. For the other values of , more calculations should be further performed.
For EI arbitrarily oriented in plane (e.g. the x-y plane), the averaged RVE can be viewed as a transversely isotropic material, with each inclusion having the rotary axis z. As shown in figure 5c, the MT predictions of effective elastic moduli along the x axis (Ē 1 ) and the z axis (Ē 3 ) agree well with the FE simulations for the RVE including 8 and 64 EI. This again demonstrates the reliability of the MT model. Previous results also show that the KT model only applies in the case where the inclusions have arbitrary orientations (i.e. random distribution), rather than in plane. This shows the MT model has a wider range of applications including in-plane arbitrary orientations.
For single directional EI along a coordinate axis, the resulting RVE is definitely a transversely isotropic material. As shown in figure 5d, the MT predictions of effective elastic moduli along the x axis (Ē 1 ) and the z axis (Ē 3 ) are consistent with the FE simulations for the RVE including 1 and 64 EI. The maximum error ofĒ 3 between the MT prediction and the FE simulation is only 4.305% for c 1 = 0.2. This again indicates that the MT model can predict the mechanical behaviors of anisotropic materials.

Typical inclusions for convenience of engineering usage: coin, disk and needle
In engineering, as the formula of ellipsoidal inclusion is often very complicated when calculating the effective modulus, three simplified formulas for the inclusions with extreme shapes are often provided, including the coin inclusion ( needle inclusion (NI) and disk inclusion (DI). From the viewpoint of geometric configuration, the CI, NI and DI correspond to EI when << 1, → ∞ and → 0, respectively. Although these simplified formulas have been widely used in engineering for convenience, their comparison with the formulas of EI, which has a strong theoretical background, has rarely been explored. In the KT model, the three kinds of inclusion with extreme shapes actually correspond to the simplified formulas of three different shape factors P 0, i and Q 0, i respectively, see Appendix C. As a result, the symbols Eq-KT-CI, Eq-KT-DI and Eq-KT-NI are used to denote corresponding simplified formula with respect to CI, DI and NI based on the KT model, respectively. Similarly, different extreme inclusion shapes correspond to different S tensors in the MT model, and Eq-MT-CI, Eq-MT-DI and Eq-MT-NI represent the corresponding simplified formulas of CI, DI and NI based on the MT model, respectively. For convenience of description, the classical ellipsoid formulas in the KT and MT models, i.e. the reduced formulas, are also designated as Eq-KT-EI and Eq-MT-EI, respectively.

Comparison between the reduced KT model and simplified formulas
Previous results (Liu & Sun 2015;Zhang et al. 2020c) show that the Eq-KT-CI only applies when the parameter is very small. When the value of is sufficiently big, the Eq-KT-CI has no meaning and cannot be used any more. However, the real range of when the Eq-KT-CI can be properly used has not been quantitatively given, and this is beneficial to reasonable selection of this model. Therefore, the calculation results of effective modulus of sandstone based on the Eq-KT-EI are compared with those with Eq-KT-CI, Eq-KT-NI and Eq-KT-DI, which are shown in figure 6a, b, c and d.
As shown in figure 6a and b, for the predicted effective modulus, the difference between the results calculated by the Eq-KT-CI and that by the Eq-KT-EI increases with the increase of , which can be seen from the differences between lines 1-2, 3-4, 5-6 and 7-8. In particular, when = 0.05 and c 1 = 0.25, the error ofḠ calculated by the two models is 2.375 GPa, which cannot be ignored. Similarly, the difference related withK andḠ between the Eq-KT-NI (line 4) and Eq-KT-EI (lines 5 and 6) both increase slightly with the decrease 546 of , as displayed in figure 6c and d. By comparing the calculation results on the effective modulus, it is suggested that the Eq-KT-CI and Eq-KT-NI can be used when < 0.05 and > 5, respectively. Furthermore, only when < 0.01, such as when = 0.001 in figure 5c and d, the calculation result based on the Eq-KT-EI (line 2) is close to that of the Eq-KT-DI (line 1), as shown in figure 6c and d. During this situation, the values ofK andḠ become negative when the volume fraction is very small, which can be seen from the lines 1-4 in figure 6c and d, and this again indicates that there are flaws in the KT model. This shows that the Eq-KT-DI has a very narrow range of applicability and few practical situations conform to this configuration.

Comparison between the reduced MT model and simplified formulas
As mentioned above, the MT model can be used to analyse the composite material where the embedded inclusions have three typical orientations. Figure 7 shows the comparison of results between the Eq-MT-EI and the three simplified models, where three orientations of inclusions in space are considered. In practice, the real range of when the Eq-MT-CI can be used has not been quantitatively given. Therefore, the calculation results on effective moduli of sandstone with Eq-MT-EI are compared with those with Eq-MT-CI, Eq-MT-NI and Eq-MT-DI, which are shown in figure 7, figure 8 and figure 9, respectively. figure 7a and b, for the predicted effective modulus of sandstone with inclusions of arbitrary orientations, the difference between the results calculated by the Eq-MT-CI and that by the Eq-MT-EI increases with the increase of , which can be seen from the difference between lines 1-2, 3-4 and 5-6. Although the difference is smaller than that using the KT model, the value ofK predicted by Eq-CI-MT is obviously wrong when = 0.5, as it even exceeds the HS upper bound. Similar to the KT model, for both the values ofK andḠ, the difference between the Eq-MT-NI (line 4) and the Eq-MT-EI (lines 5 and 6) increases slightly with the decrease of , as shown in figure 7c and d. By comparing the calculation results on the effective modulus, it is suggested that the Eq-MT-CI and Eq-MT-NI can be used when < 0.1 and > 5, respectively. The difference onḠ between the Eq-MT-DI (line 1) and Eq-MT-EI (lines 2 and 3) is obviously larger than that onK, which is shown in figure 7c and d. This indicates that the Eq-MT-DI is more rigorous in predictingḠ thanK. Combining the results ofK andḠ, it is suggested that the applicable condition of Eq-MT-DI in this case is that < 0.01.

Inclusions with arbitrary orientations in the x-y plane.
The effective moduli of the sandstone with inclusions of arbitrary orientations in the x-y plane are shown in figure 8, where the results based on Eq-MT-CI, Eq-MT-NI and Eq-MT-DI are compared with Eq-MT-EI. The difference between the results calculated by Eq-MT-CI and Eq-MT-EI increase with the increase of , which can be found from the difference between lines 1-2, 3-4 and 5-6 in figure 8a, b, c and d. Among the results of four effective moduli, the difference between two the formulas onĒ 3 is the biggest, and the monotonic laws betweenĒ 3 and calculated by the two models are even opposite, which can be observed from lines 1-6 in figure 8b.
Combining the prediction results of the four effective moduli, when the parameter < 0.05, means that taking the Eq-MT-CI is reasonable. As can be seen from figure 8e, f, g and h, the difference related with four effective moduli based on the Eq-MT-NI (line 4) and Eq-MT-EI (lines 5 and 6) increases slightly with the decrease of . As the value of decreases to 1, the effective modulus calculated by Eq-MT-NI tends to be that calculated by the spherical formula. The results show that the transverse isotropy of RVE caused by in-plane arbitrarily oriented NIs is very small. Therefore, Eq-MT-NI is recommended to be used when > 5. For the effective moduli on Eq-MT-DI (line 1) and Eq-MT-EI (lines 2 and 3) in figure 8e, f, g and h, the difference onĒ 1 ,Ḡ 23 andḠ 12 is significant with the increase of , whileĒ 3 calculated by Eq-MT-EI decrease slightly with the increase of . It is found that only when < 0.005, the shape of inclusion can be considered to satisfy the disk condition. On the other hand, it can be clearly seen in figure 8 that the anisotropy of RVE caused by flat typed inclusions ( << 1) is significantly bigger than that caused by slender type inclusions ( >> 1) in the case of arbitrary orientation in the x-y plane. In the calculation results of flat typed inclusions ( << 1), the anisotropy of elastic modulus is bigger than that of the shear modulus.

Inclusions with one single orientation.
For a more anisotropic case, the effective modulus of the sandstone with one single orientation of inclusions along the x direction is shown in figure 9, where the Eq-MT-CI, Eq-MT-NI and Eq-MT-DI are compared with Eq-MT-EI. The difference between the results calculated by Eq-MT-CI (lines 2, 4 and 6) and Eq-MT-EI (lines 1, 3 and 5) increases with the increase of , which can be noticed from figure 9a, b, c and d. It is reasonable that the laws of difference between the two formulas are similar to those when the inclusions have arbitrary orientations in plane. However, by comparing figure 9c, d and figure 8c, d, it can be clearly found that the anisotropy of elastic and shear modulus caused by the single oriented inclusions is obviously bigger than that caused by the arbitrarily oriented inclusions in the plane. Especially, the shapes of lines 1-6 between figure 8c and figure 9c are different. This phenomenon also appears in other cases. According to the calculation results for each effective modulus, Eq-MT-CI can be used when < 0.05.
The comparison of results on Eq-MT-DI (line 1) and Eq-MT-EI (lines 2 and 3) are shown in figure 9e, f, g and h. It can be seen that when = 0.01, the difference between Eq-MT-EI and Eq-MT-DI in calculatingĒ 1 ,Ē 3 andḠ 23 are very small, as shown in figure 9e, f and h. Although there is still a big gap between the two formulas in calculatingḠ 12 , the result calculated by Eq-MT-DI drops rapidly to nearly 0.001 GPa with a very small c 1 , and this is almost impossible in actual reservoir rocks. Therefore, Eq-MT-DI can be applied when < 0.01. In the comparison of Eq-MT-NI (line 4) and Eq-MT-EI (lines 5 and 6) in figure 9e, f, g and h, the difference inĒ 1 , E 3 andḠ 23 is bigger than that in figure 8e, f, g and h. However, the difference in effective moduli between Eq-MT-NI and Eq-MT-EI with = 5 is acceptable, so when > 5 Eq-MT-NI can still be used.
Altogether, based on the analyses from figure 6 to figure 9, the range of corresponding to three typical inclusions according to the KT and MT models can be summarised in Table 3, which can provide some guidance for the rational selection of a simplified formula to treat mechanical behaviors of rocks.

Conclusions
We compared the KT and the MT models for the calculation of the effective moduli of sandstone, and the results are discussed in detail. Accounting for six components, either the KT or the MT model predicts experimental results for a shale well. Then the sandstone composed of quartz and pore is comprehensively studied. For SI, the results predicted by the KT and MT models are almost the same. However, for the EI, the KT model has some defects as the critical volume fraction of inclusions is normally smaller than one. Compared with the KT model, the MT model can predict the effective moduli of anisotropic sandstone with inclusions of different orientations in space, which has also been verified by the FEM simulation. The applicable scopes of several simplified formulas, such as Eq-KT-CI, Eq-KT-DI, Eq-KT-NI, Eq-MT-CI, Eq-MT-DI and Eq-MT-NI are discussed, which are grouped into one table feasible in engineering. The orientation distribution of inclusions will reduce the scope of application of the simplified formula. This basic research can deepen the understanding of the two models and provide guidance for the application of the two models in rock physics and other engineering fields.

Appendix B.
In the MT model, in the calculation of the effective modulus of the pore inclusions in the arbitrary orientation in the plane and the arbitrary orientation in the space, the angle conversion formula is used to convert the stiffness matrix of the local coordinate system. The fourth-order tensor corresponds to 6 × 6 and the angle conversion formula of the matrix can be written as D = TD ′ T −1 , (B1) where D is the 6 × 6 matrix in the global coordinate system, and D ′ is the 6 × 6 matrix in the local coordinate system. T can be expressed as where Q is the angle conversion matrix of the vector, which is Q = ⎡ ⎢ ⎢ ⎣ cos sin cos sin sin −sin cos cos cos sin 0 − sin cos where is a polar angle, and is an azimuthal angle. Therefore, equation (B1) can be expressed as Then, the mean value of tensor with arbitrary orientation in space can be expressed as The mean value of tensor with arbitrary orientation in x-y plane can be expressed as (2) Ellipsoidal inclusion (rotation axis is x, a 2 = a 3 = a, aspect ratio = a 1 /a) where )] , ] , In equation (C3), A, B and R are, respectively, The value of g is the same as equation (A3), and f = 2 1 − 2 (3g − 2) .