Construction of a novel brittleness index equation and analysis of anisotropic brittleness characteristics for unconventional shale formations

The brittleness prediction of shale formations is of interest to researchers nowadays. Conventional methods of brittleness prediction are usually based on isotropic models while shale is anisotropic. In order to obtain a better prediction of shale brittleness, our study firstly proposed a novel brittleness index equation based on the Voigt–Reuss–Hill average, which combines two classical isotropic methods. The proposed method introduces upper and lower brittleness bounds, which take the uncertainty of brittleness prediction into consideration. In addition, this method can give us acceptable predictions by using limited input values. Secondly, an anisotropic rock physics model was constructed. Two parameters were introduced into our model, which can be used to simulate the lamination of clay minerals and the dip angle of formation. In addition, rock physics templates have been built to analyze the sensitivity of brittleness parameters. Finally, the effects of kerogen, pore structure, clay lamination and shale formation dip have been investigated in terms of anisotropy. The prediction shows that the vertical/horizontal Young’s modulus is always below one while the vertical/horizontal Poisson’s ratio (PR) can be either greater or less than 1. Our study finds different degrees of shale lamination may be the explanation for the random distribution of Vani (the ratio of vertical PR to horizontal PR).


Introduction
Shale reservoir exploration has become an important research area in recent years. The industrial scale production of shale oil/gas has led to a revolution in hydrocarbon exploration. Unlike conventional sandy reservoirs, shale formations are typically low porosity and low permeability, which means that conventional methods are no longer effective (Johnston and Christensen 2012). Deep learning and artificial intelligence have become more and more popular in hydrocarbon exploration (Lin et al. 2018), while the development of engineering technology brings more commercial potential into unconventional exploration. For instance, to optimize shale pore space, hydraulic fracturing has been widely implemented (Rickman et al. 2008). Through injecting fluids into the borehole, artificial fractures can be formed in the shale formations, and then, proppant is used to keep the fractures open, all of which improves the porosity and permeability of the target area.
Brittleness, a significant parameter of a shale "sweet spot," can play a critical role in hydraulic fracturing and horizontal well trajectory design (Jarvie et al. 2007). Brittleness is a research topic in many fields, such as petroleum engineering (drilling, fracturing, etc.), geology, geochemistry and geophysics. Different fields have their own definitions of brittleness (Rickman et al. 2008). The literature Edited by Jie Hao and Xiu-Qiu Peng shows that there are dozens of equations that can be used to define brittleness. In a variety of brittleness evaluation systems, laboratory measurements result in the most accurate brittleness values. However, laboratory measurements cannot be widely used owing to high coring costs and limited sampling. In order to estimate brittleness more economically and obtain 2D/3D brittleness distributions, using geophysical data (seismic data or well logging data) in brittleness prediction seems to be a straightforward choice (Goodway et al. 1999(Goodway et al. , 2010. Unfortunately, brittleness cannot be measured directly from either well log data or seismic data (Li et al. 2012). Therefore, geophysicists usually use elastic data and the volume fraction of minerals to characterize the brittleness indirectly with the help of a brittleness index equation. The advantage of geophysical data is that it can provide 2D/3D elastic information (such as velocity, density, etc.) on the target area, which can help us to model 2D/3D brittleness distributions. Furthermore, different elastic parameters show different sensitivity to brittleness, and the brittleness index (BI) can be used to reconstruct different elastic parameters and make them more sensitive to brittleness. This means that construction of a reasonable brittleness index is a significant factor in brittleness prediction (Huang et al. 2015). If the elements of the BI are elastic parameters, we call the brittleness prediction based on this category of BI-"the elastic parameter-based method." "The mineralogy-based method" is another category of BI, which obviously means the elements in BI equations are volume fractions of different minerals. The mineralogy method uses the volume percentage of brittle minerals to characterize the rock brittleness, since the brittleness of rock will be proportional to the volume fraction of brittle minerals (Liu and Sun 2015).
The existing brittleness index formulas are subdivided into the mineralogy method or the elastic parameter method; very few methods consider both two categories simultaneously. Both methods have their own advantages and disadvantages. On the one hand, the elastic parameter method seems to give a more accurate result, since many elastic parameters (i.e., Young's modulus, Lamé constants, etc.) can describe brittleness and these elastic parameters contain the information on both mineralogy and pore fluid. However, the accuracy of predictions based on this method is subject to the accuracy of the input elastic parameters. On the other hand, the accuracy of predictions based on the mineralogy method is not so good by comparison, but this method is still widely used in oil industry since the input parameters are easy to obtain. Due to the high cost of shale exploration wells, the required data is scarce. The accuracy and reliability of brittleness predictions of shale formations need to be improved. Hence, in this paper, we attempt to combine these two types of brittleness formulas and strive to obtain a novel BI equation with higher accuracy brittleness predictions while simultaneously using less input data.
Admittedly, the brittleness index can provide us an effective way to predict brittleness, but brittleness indexes are mostly based on the assumption of isotropy (Huang et al. 2015). It is widely known that shale is highly anisotropic. We believe therefore that we should consider the properties of shale brittleness in anisotropic terms. As shown in Fig. 1, the well trajectory of a horizontal well in a shale formation usually passes through anisotropic formations (e.g., vertical, tilted or horizontal transverse isotropic formations, VTI, TTI or HTI). Ignoring shale anisotropy will cause errors in brittleness estimation. In the oil industry, however, the common way to evaluate shale brittleness is based on using isotropic parameters (i.e., Young's modulus and Poisson's ratio), since these can be easily calculated from measured elastic wave velocities from geophysical data. The rock physics model has been proven as a useful tool for reservoir characterization, which can be used to estimate the brittleness of shale reservoirs (Hornby et al. 1994;Vanorio et al. 2008). Anisotropic rock physics modeling offers us a new way to estimate the brittleness of shale in an anisotropic way. In order to build a reasonable anisotropic shale model, we need to consider the possible origins of shale anisotropy. Based on published research, several factors contribute the anisotropy of shale. The causes of shale anisotropy can be summarized as follows: (1) A preferred orientation of clay and kerogen particles can be used to explain the VTI quality (Vernik and Nur 1992;Vernik and Landis 1996;Vernik and Liu 1997;Johansen et al. 2004); (2) SEM (scanning electron microscopy) and XRD (X-ray diffraction) observations show that sometimes shale pores are typically flattened with associated microfractures, which may further enhances shale anisotropy; (3) Clay minerals have their own intrinsic anisotropy (Sayers 1994(Sayers , 2005Sondergeld and Rai 2011;Guo and Li 2015). In our study, we intend to construct a rock physics model, which can take the different causes of shale anisotropy into consideration.
Rock physics templates (RPTs) are an effective tool to link elastic properties with reservoir properties in terms of a rock physics model. Guo et al. (2012a, b) constructed RPTs to analyze geo-mechanical properties in terms of seismic attributes. Avseth and Carcione (2015) built different RPTs to study the effects of kerogen content, pore fluid saturation and in situ pressure.
Therefore, in order to analyze the anisotropic properties of shale brittleness, we firstly create an anisotropic shale rock physics model to build a bridge between anisotropic elastic parameters and physical properties. Following on from this, we construct rock physics templates based on the model in order to analyze the sensitivity of brittleness parameters. Finally, we discuss anisotropic brittleness in terms of physical properties (lamination, pore structure, etc.).

The classic brittleness index equation
In order to describe the brittleness of a shale reservoir, various brittleness equations have been proposed (Li et al. 2012). In this chapter, we pay attention to their respective characteristics and try to reshape them to allow for extraction of more information on brittleness.
As described above, the commonly used brittleness equations are mainly based on the weight content of minerals and elastic parameters, as shown in Eqs. (1)-(9). BI_1 to BI_3 and BI_6 to BI_8 belong to the group of elastic parameters-based brittleness indexes, while BI_4 and BI_5 belong to mineralogy-based brittleness indexes. Rickman et al. (2008) propose an average brittleness equation based on the analysis of the Barnett shale: where YM is Young's modulus (GPa), PR represents Poisson's ratio (dimensionless). ΔYM and ΔPR represent the normalized Young's modulus and normalized Poisson's ratio, respectively. The real data range of Barnett shale Young's modulus and Poisson's ratio determines the constants in the equations.
The general form of Rickman's equation is that.
where YM 3 and PR 3 correspond to ΔYM and ΔPR, respectively. YM 1 and YM 2 represent the maximum and minimum value of Young's modulus for the target formation, respectively. PR 1 and PR 2 represent the maximum and minimum values of Poisson's ratio for the target formation, respectively. Guo et al. (2012a) constructed a brittleness index based on the opinion that the brittleness index is proportional to Young's modulus and inversely proportional to Poisson's ratio. Liu and Sun (2015) defined the brittleness index by the ratio of normalized Young's modulus and normalized Poisson's ratio.
where YM 3 and PR 3 represent the normalized Young's modulus and normalized Poisson's ratio, respectively.
The following Eqs. BI_4 and BI_5 are the brittleness index based on mineralogy-based equations: (1) (2) where Quartz, Total, Calcite represent the weight content of quartz, total minerals and calcite, respectively. Guo et al. (2012a) defined the brittleness index using Lamé parameters: Chen et al. (2014) indicated that the ratio of Young's modulus and Lamé parameters can represent the highly brittle strata: Huang et al. (2015) proposed a more sensitive brittleness index: where K is bulk modulus.
As discussed above, the various brittleness indexes proposed have different definitions and have been derived for different purposes. We should choose the most appropriate brittleness index to arrive at better brittleness evaluations.
BI_1 to BI_3 are all brittleness indexes defined by Young's modulus and Poisson's ratio. They all show that brittleness index is proportional to Young's modulus and inversely proportional to Poisson's ratio. However, the weighted average in BI_1 does not consider the difference between Young's modulus and Poisson's ratio. BI_1 should be further verified by more shale data to prove its feasibility in different scenarios. BI_2 and BI_3 define brittleness index using the ratio of Young's modulus and Poisson's ratio. The normalizations in BI_3 avoid outliers and improve the stability of calculation. However, the normalizations can also lead to a decrease in sensitivity and zero values of Poisson's ratio, which is the reason why a singular value occurred.
Besides, BI_1 to BI_3 and BI_6 to BI_8 have different magnitudes so we cannot compare these brittleness indexes using a uniform reference system.
BI_4 and BI_5 are all mineralogy-based brittleness indexes. BI_4 regards quartz as the brittle mineral, while BI_5 regards quartz and calcite as brittle minerals. BI_4 and BI_5 are simple to calculate, but the given brittleness indexes are relatively rough and inaccurate when compared to elastic parameter based brittleness indexes.

An improved brittleness equation based on the Voigt-Reuss-Hill average theory
As shown in Fig. 2, Fig. 2a, b shows two hypothetical rocks. Sample (a) contains 40% quartz and 10% calcite, and sample (b) consists of 40% calcite and 10% quartz. If the brittleness value of these two rocks is found to be the same using the formula BI_5, the mineral composition of the two rocks is different and therefore the brittleness is obviously different. This means that the formula considering the mineral percentages has some limitations and fails to predict the accurate brittleness under some circumstances. As mentioned above, the advantage of the mineral composition method is the requirement for fewer input parameters, but the accuracy of predictions is low. While the advantage of the elastic parameter method is the high prediction accuracy, but more input parameters are required.
The proposed hybrid algorithm in this paper combines the two methods. Based on the mineral composition method, the precision of our new method approaches to elastic parameter-based method while the input parameters are only the mineral composition of the shale.
Shale is usually composed of a variety of minerals, each having different brittleness. Minerals such as quartz and calcite have usually been considered as brittle minerals, and clay and kerogens are considered as non-brittle minerals. If we use the mineralogy-based method, the selection of the brittle mineral will have an impact on the final brittleness prediction, for example, BI_4 considers that the brittleness of the rock is only related to the quartz content, while BI_5 considers that the brittleness of the rock is related to both quartz and calcite content. In practice, all the various mineral components have an impact on the brittleness of the rock. Hence, we use elasticity information such as the velocity and density of pure minerals (as shown in Table 1) to calculate the brittleness of each pure mineral component (as shown in Table 2). Since the quartz content is usually known as the most important indicator affecting the brittleness index of a rock, the brittleness index of quartz is assumed to one, and the "pure brittleness index" of other minerals is calculated relative to quartz. Finally, the brittleness index is been obtained by bringing the pure brittleness index into the VRH brittleness equation by using the percentage content of each mineral.
When the available information is limited, in order to obtain more realistic prediction results, petrophysicists usually use the approach of determining the upper and lower limits of the target parameters by simulating the extreme situations (Mavko et al. 2009). This leads to the prediction of the reasonable ranges of the target parameter's upper and lower limits. In the absence of geometric details relating to the mutual coupling of each mineral, the Voigt-Reuss-Hill theory (Mavko et al. 2009) is often been used to predict the equivalent elastic modulus of a mixture of minerals and pores.
When the input information is only the "pure brittleness index" and volume fraction of each mineral component, the specific geometric details of the mineral coupling are missing. The upper and lower limits of the equivalent brittleness are analyzed by referring to the idea of VRH upper and lower limits, and the coupling state of each mineral is assumed to be the most extreme series-parallel connection. In addition, using the Hill average, the final brittleness prediction result is the average of the upper and lower brittleness. Figure 3 shows the Voigt model and Reuss model schematically; the figure shows two hypothetical minerals A and B in gray and white. Under the same vertical pressure, the Voigt model simulates parallel connection of the two minerals, while the Reuss model simulates serial connection of the two minerals. When the minerals A and B were mixed as in the Voigt model, the equivalent elasticity of the mixed mineral calculated by the Voigt formula is the upper limit of the modulus. Correspondingly, when the minerals A and B are mixed in series as in the Reuss model, the equivalent result is usually the lower limit. Therefore, this study attempts to use the brittleness index as a "modulus." By assuming that the coupling geometry of each mineral component is in either series or parallel as shown in Fig. 3, considering that mineral coupling details are missing and the input parameters are limited, the construction of brittleness upper and lower limits allows the prediction of the reasonable range of rock brittleness.
The VRH brittleness formula is been shown below: In the formula, i refers to the different mineral composition, f i is the volume content of the mineral, PBI i is the relative brittleness index of the mineral, BI_V and BI_R are the upper and lower limits of the brittleness, respectively. BI_VRH is the Hill average of the brittleness.
The PBI i value is shown in Fig. 4, and we calculate the relative brittleness indices of six common minerals in shale reservoirs. Firstly, the elastic modulus and shear modulus of each mineral have been calculated using the velocity and density of each mineral as shown in Table 1. Then we use formula 2 and formula 3 to calculate the corresponding pure mineral brittleness index. Finally, the brittleness of the quartz is assumed to be one, and the brittleness of other minerals is normalized to this to obtain the final relative brittleness index (as shown in Table 2).
It can be seen from Fig. 4 that BI_2 (Blue line, see Eq. 3) is more sensitive to mineral composition than BI_1 (Red line, see Eq. 2). Therefore, this paper uses BI_2 to calculate the brittleness factor. The corresponding brittleness factors after normalization are brought into formulas 10-12 to obtain the corresponding brittleness indexes. It can be seen from formulas 10-12 that the BI_V has a similar form to the upper limit of the Voigt model and describes the brittle exponential representation of each component in the case of equal strain. Therefore, this paper tries to follow the Voigt-Reuss-Hill concept, assuming that the brittleness index of BI_V is upper limit of brittleness, under the same stress condition, the brittleness index of BI_R is lower limit of the brittleness and the brittleness index of BI_VRH is the average value of these two limits.

Case study
In order to verify the applicability of our equation, we apply it to real data from China to predict brittleness. We also compare the results of our equation with those of conventional equations to assess the superiority of our method. Figure 5 shows well logs from a well in Sichuan province, China. The target layer, in the blue box, is located in the Longmaxi and Wufeng Formations, whose depth is between 2045 and 2065 m. The logs in Fig. 5, from left to right are GR, Porosity, Vp, Vs, TOC, Density, Vp/Vs and Poisson's ratio. In the left column, the red curve is GR and the blue is GR without uranium. The kerogen's radioactivity comes mainly from uranium. The gap between the red and the blue curves can be used to identify shale reservoirs. The bigger the gap is, the more organic content there will be. In our case, other well logs also show the obvious characteristics of shale reservoirs. Since the shale layer contains more clay and kerogen, and the elastic modulus (bulk modulus and shear modulus) of clay and kerogen is relatively low, the Vp, Vs and density values are often lower than in the surrounding rock. Figure 6 shows the mineral composition of the target. From Fig. 6, we know quartz and calcite content is higher in the target, which indicates the target has a high brittleness; it is therefore classified as a high-quality reservoir interval.
We calculated brittleness indexes based on brittleness Eqs. 1-12. The results are shown in Fig. 7. In Fig. 7, columns 3-5 are log interpretation results, and columns 6-10 are distributions of brittleness indexes calculated from different equations. BI_2 and BI_3 in column 6 are brittleness indexes based on Young's modulus and Poisson's ratio. Though these two curves show differences in the Wufeng Formation, the overall trend is consistent; this confirms that these two equations are applicable. BI_4 and BI_5 in column 7 are brittleness indexes based on mineral components. We can see that these two curves are less correlated, and this indicates that the equation based on mineral components is less accurate than the equation based on elasticity parameters.
Column 8 shows the comparison between brittleness indexes BI_VRH and BI_3. The red curve is brittleness index BI_VRH; the black curves show the upper and lower limits of brittleness index BI_V and BI_R.
Column 9 shows the comparison between brittleness indexes BI_VRH and BI_4. It is worth noting the trend of brittleness indexes based on mineral percentage is consistent with the trend of the upper limits BI_V curve. BI_6, BI_7 and BI_8 shown in column 10 are calculated based on Lamé parameters λ and μ. From the figure, we can see that these three curves have similar trends.
Brittleness indexes shown in column 8 and 9 are calculated based on our VRH equation. BI_VRH and BI_3 show a better match, both these curves lie between the upper and lower limits of brittleness index. This confirms the applicability of brittleness index BI_VRH. Across the majority of the target depth interval denoted by the blue box, the accuracy of our new brittleness index BI_VRH (red curve) is similar to brittleness index curve BI_3 based on elasticity parameters (blue curve). When calculating BI_VRH, required input parameters are mineral percentage and brittleness index minerals. Therefore, using our new VRH brittleness equation, we can obtain an accurate brittleness index without requiring elastic parameter data.
It is worth mentioning that brittleness indexes based on mineral percentages are consistent with those related to upper limits BI_V. This is because brittleness indexes based on mineral percentages ignore fluids in pores. Pore It can be found that at a depth of 2055 m, our results are not accurate. According to log interpretation results, at this depth, rock porosity is high and pore fluid effects are complex. It is also known, based on other test reports, that anisotropy at this depth is high. These may be the reasons for the increase in error. Brittleness index BI_VRH considers mineral fusion, but ignores rock porosity and fluid effects, indicating that our new equation is most suitable for application in intervals with low porosity. When porosity is high and fluid composition and pore structure is complex, the applicability of our equation needs further study.

Anisotropic characteristics of shale brittleness
Shale usually shows high anisotropy. The isotropic brittleness index is insufficient to describe the anisotropic characteristics of shale formations. Young's modulus and Poisson's ratio are two typical parameters, which are sensitive to brittleness. In order to calculate the anisotropic YM and PR, a 6*6 elastic constant matrix is needed for calculation based on a rock physics model. Therefore, in this section, we initially construct an anisotropic rock physics model. Following this, we build a rock physics template in order to discover the parameters sensitive to brittleness. Finally, we derive the relationship between the physical properties, (kerogen content, lamination, pore structure and layer dip) and anisotropic brittleness.  (2013) performed an extensive analysis of published shale data. He found that the vertical YM is typically lower than the horizontal YM, and vertical PR can be either higher or lower than horizontal PR. Based on our shale rock physics model analysis, our result indicates this observed phenomenon may due to the physical properties of shale. In the following section, this phenomenon was explained using our anisotropic brittleness analysis.

Construction of an anisotropic shale rock physics model
In order to discover the anisotropic behavior of shale brittleness, firstly, we need to construct an anisotropic model. The model we use is based on our previous work (Qian et al. 2014(Qian et al. , 2016. The main steps can be summarized as follows: (1) Kerogen-and clay-related pores are added into pure clay successively to build the clay-kerogen block using SCA + DEM (SCA is short for self-consistent approximation and DEM is short for differential effective medium. For more details of isotropic SCA/DEM and anisotropic SCA/DEM please see the paper by Hornby et al. 1994) (2) Many identical blocks are rotated to different angles; the range of angles is normally distributed. Bond transform is used to calculate the elastic stiffness of all the rotated blocks, which are then combined using a VRH average to form the clay-kerogen background with preferred orientation. The angles of blocks are decided by ODF. (ODF is short for orientation density function.) (3) The elastic stiffness of other minerals present in rock is modeled using the VRH average. Then, brittle mineral-related porosity is added using the isotropic SCA method. (4) Considering the preferentially orientated clay-kerogen mixture as background, we use anisotropic DEM to add the isotropic brittle mixture into the background, forming the final effective mixture Our previous work has shown that our organic-rich model is suitable of modeling the organic-rich shale properly. However, in real oil fields, a tough problem we need to face is that the mineralogy of the target well may vary significantly from shallow zone to deep zone, which means the organic-rich shale model may not be suitable for all the target zones. The rock physics workflow is decided based on the characteristics of the target rock type. Moreover, the effective theories that we used in the workflow need to satisfy several assumptions. Hence, our basic idea is that we can build a variable rock physics workflow that can choose the most reasonable work path according to the type of rock. Different work paths perform different effective theories, which give the best simulation of different types of rock.
In order to enhance the generalization ability of our model, we develop our organic-rich shale model into a general edition (shown in Fig. 8). The core feature of our new workflow is that our model can select a reasonable path according to the mineralogy of the formation. Two pink diamonds in Fig. 8 are two crossroads.
We build three paths here: (a) As the volume fraction of clay content is below 50%, the workflow chooses path one. In this scenario, the brittle minerals are the majority material of rock. Therefore, we perform anisotropic DEM to simulate an effective mixture with brittle mineral background and clay inclusion; (b) as volume fraction of clay content is above 50% and that of kerogen content below 3%, this scenario is some kind of mudstone blocking layer. The clay is treated as background material and the brittle minerals are treated as inclusions by using anisotropic DEM; (c) as volume fractions of clay and kerogen are both relatively high, the workflow chooses path three, which indicates the scenario of organic-rich shale and this type of path shows a similar workflow to our previous work (Qian et al. 2014(Qian et al. , 2016. In order to check the feasibility of our workflow, we perform our modified workflow to our target well and predict its elastic wave velocity, by performing the inversion procedure in our previous paper (Qian et al. 2016). Our predicted P and S wave velocities are shown in Fig. 9. Figure 9a is the velocity predictions based on our previously constructed model which is only suitable for organic-rich shale. Meanwhile, Fig. 9b shows the predictions based on our new modified workflow. Blue and red curves indicate measured and predicted velocities.
We divide our target zones into two main types and use pink and green arrows to notice them. The zones highlighted by a pink arrow are a typical kerogen-rich formation, which has relatively high porosity, high TOC and high velocity. Moreover, the green arrows indicate the zones with high brittle mineral content, which have relatively low porosity, low TOC and low velocity.
It is notable that both models are suitable for organicrich formation (pink arrow), while for brittle mineral-rich formation (green arrow), our new workflow shows better predictions, since we consider different workflow paths for different type of rock. Hence, our generalized shale rock physics workflow shows good feasibility and applicability, especially for formations with significant lithology change.

Construction of the rock physics template
In this section, rock physics templates are built based on our model to illustrate the influence of lithology and porosity on shale brittleness. Figure 10 shows the crossplots of λρ-μρ, Young's modulus (YM) versus Poisson ratio (PR) and Vp/Vs versus Ip (P wave impedance). The mesh in the figure is controlled by two parameters: clay content and porosity. The dashed lines indicate the effect of mineralogy, while the solid lines represent the effect of porosity. The dashed lines show that the clay content is gradually substituted by quartz, and other minerals remain unchanged. The clay content increases from 10% to 40% with an interval of 10%, and the quartz decreases from 40% to 10%. Meanwhile, the solid lines show that the porosity varies from 0.08 to 0.02 with step equals to 0.02. Figure 10a indicates that λρ is more sensitive to brittle mineral content than μρ especially for the low porosity situation. Figure 10b shows that high quartz content and low porosity correspond to high YM and low PR, which indicates that brittleness is affected not only by mineralogy but also by porosity. YM appears to be  Figure 10c shows that Vp/Vs decreases as quartz content increases.
As described in Sect. 3.1, our rock physics model first builds a clay-kerogen block and then rotates the block to model the inclination observed in the formation. Figure 11a shows the influence of the alignment angle of the clay-kerogen block on the anisotropic parameters. The angle varies from 0 to 90°. Figure 11a indicates that the Thomsen anisotropic parameters ε and γ decrease as angle increases, while δ reaches its peak at around 40°. ε and γ are always above zero, while δ can be either positive or negative, which agrees with Sayers's study. In addition, ε and γ have a linear relationship while δ seems to be independent. Figure 11b indicates the relationship between inclination angle and YM and PR. High YM and low PR usually represent high brittleness. The peak value of brittleness appears when the orientation angle around 45°.

Effect of pore geometry and kerogen
Our previous work has proven that pore geometry and kerogen are two significant factors that influence elastic properties. Based on our model, we link the physical properties of shale, i.e., pore geometry (pore aspect ratio) and kerogen content, to the anisotropic sensitive brittleness parameters (YM and PR). Figure 12 shows the predicted anisotropic YM and PR based on our model resulting from varying pore aspect ratios (from 0.1 to 1) and the volume fraction of kerogen (from 0 to 1). The elastic properties we use for kerogen and "shale" (the mixture of the other minerals) are quoted from Sayers (2013) as below. For kerogen C 11 = C 33 = 9.8GPa, C 44 = C 66 = 3.2GPa and C 13 = 3.4GPa. Moreover, for shale C 11 = 85.6GPa, C 33 = 65.5GPa, C 44 = 24.6GPa, C 66 = 29.7GPa and C 13 = 21.1GPa. Figure 12a, b indicates the ratio between vertical and horizontal YM and PR; Fig. 12a shows the anisotropic YM: E ani (the ratio of vertical YM to horizontal YM) decreases as pore aspect ratio decreases, since low pore aspect ratio leads to high shale anisotropy. As volume fraction of kerogen increases to around 50%, E ani becomes more sensitive, which may be due to the "interaction" between the shale phase and kerogen phase reaching a maximum as both phases have similar volume fractions. Figure 12c-e indicates the calculated brittleness based on Rickman's equation (Rickman et al. 2008). Modeling results show BI_horizontal (horizontal brittleness) is more sensitive to a low pore aspect ratio, while BI_vertical, (vertical brittleness) is more sensitive to a high pore aspect ratio. The BI_vertical/BI_horizontal ratio shows similar characteristics to the YM and PR ratio.

Effect of lamination
In our model, we use bond transforms and orientation density functions to represent the lamination of shale. Here, we vary the extent of lamination and analyze the anisotropic YM and PR. As shown in Fig. 13, as the lamination index increases, the extent of lamination decreases, the shale becomes more randomly distributed. Hence, the E ani value approaches one, which indicates a decrease in shale anisotropy. It is notable that, while varying the lamination index, the E ani value remains below one, while the V ani value can be greater than or less than one, which coincides with the observed phenomenon from Sayers (2013). Hence, we believe the lamination of shale can be used to explain the V ani phenomenon.

Effect of formation dip
As mentioned above, our model uses orientation density functions (ODF) to represent the lamination of the shale. The ODFs are been controlled by two main parameters. One is "expectation," and the other is "variance." "Expectation" controls the angle of the highest density of the clay-kerogen block, which represents the dip of formation. "Variance" controls the gap between the highest and the lowest densities, which can be used to model the lamination of shale. In the previous section, we set "expectation" equal to zero and vary the value of "variance" to model the effect of lamination. In this section, we vary "expectation" and set "variance" as a constant to model the dip angle of shale layers. Figure 14 shows the calculated brittleness obtained by varying the pore aspect ratio and the dip angle of shale layers. For a high pore aspect ratio, the anisotropy of shale is weak, and BI_vertical/BI_horizontal varies around 1. As pore aspect ratio decreases, anisotropy increases. The values of BI_horizontal and BI_vertical trade places with each other, since the dip angle of the shale formation varies from 0 to 90 which indicates the transformation of a VTI formation into a HTI formation.
Comparing the modeling results with the diagram of a real well trajectory (Fig. 1), different dip angles indicate different phases of the well trajectory. Considering the anisotropic information of the shale layers can help drilling engineers locate the target area more accurately.

Conclusions
In this paper, we attempt to link two brittleness equations (the mineralogy-based and elastic parameter-based methods) together and introduce the concept of "pure brittleness index," the mineral components were treated equally considering the common influence of each mineral component on the brittleness of rocks. Subsequently, based on the idea of the Voigt-Reuss-Hill average, a new brittleness equation has been introduced. The PBI has been reconstructed to predict the upper limit, lower limit and average of brittleness. The new equation has been applied to real logging data from southwest China. The predicted results show that the new equation can generate good predictions, comparable to those derived from the elastic parameter method, with only limited input information.
However, the accuracy of our equation is lower than the predicted results calculated from the elastic parameter method at some depths, which is because there is no input information about pore fluid in areas of relatively high porosity. It is clear that the applicability of the new brittleness equation in areas of relatively high porosity needs further research.
In terms of anisotropic analysis, our studies first introduce a new rock physics model for shale, based on our anisotropic shale model. The relationship between physical properties (like pore geometry, TOC, shale lamination and dip angle of shale formation) and elastic parameters (like anisotropic YM and PR) is then established. Our modeling results coincide with the real data points derived from field examples of shale. The prediction shows the vertical/horizontal YM always remains below one while the vertical/horizontal PR can be greater or less than one. We find different extents of shale lamination may be the explanation for the random distribution of V ani . Additionally, the brittleness of shale can 14 Predicted anisotropic brittleness based on our model resulting from varying pore aspect ratio (from 0.1 to 1) and dip angle of the shale formation (from 0 to 90 degree) vary significantly while drilling due to the strong anisotropy of the shale formations.