Micromechanical elastic properties solver

Abstract

This disclosure describes a novel method for predicting continuous wellbore mechanical properties such as static elastic stiffness and failure strength, where the properties solutions are deterministic and based on mechanical theory. It has at least three immediate applications: (a) continuous plots of mechanical properties vs. depth, (b) conceptual testing of the effect of changing constituent volume fractions, (c) ternary plots.

Claims

1. A method of deriving continuous mechanical properties, said method comprising: obtaining mineralogy data from drill cuttings from a wellbore using x-ray diffraction (XRD), wherein said mineralogy data comprises composition and texture parameters by depth, constituent type, volume fraction of constituents, porosity, shape factor, and elastic parameters including Young's modulus, Poisson's ratio, bulk modulus, and Lame's constants; compiling said mineralogy data by volume fraction and designating at least two of said elastic parameters; inverting said elastic parameters into a compliance tensor; converting said shape factor into Eshelby shape tensors; performing continuum micromechanical elastic properties modeling iteratively on the Eshelby shape tensors using Eshelby's inclusion method for heterogeneous composite materials; calculating an overall micro-mechanically averaged compliance tensor based on the continuum micromechanical elastic properties modeling; and displaying results as at least one of (1) continuous plots of the overall micro-mechanically averaged compliance vs. depth, (2) conceptual testing of an effect of changing the volume fraction of constituents, or (3) ternary plots of volume fraction and elastic properties, wherein the continuum micromechanical elastic properties modeling is performed with the following equation: $\stackrel{\_}{C}=C+\underset{=1}{\overset{N}{.Math.}}{f}^{}\left(C^{}-C\right){\left(\stackrel{\_}{C}-C^{}\right)}^{-1}{\stackrel{\_}{C}\left[{\left(\stackrel{\_}{C}-C^{}\right)}^{-1}\stackrel{\_}{C}-{\stackrel{\_}{S}}^{}\right]}^{-1}$ using an initial assumption of C and S.sup. to incrementally solve the equation by
C.sup.n+1=Cf.sup.(C.sup.C)[(C.sup.nC.sup.)C][(C.sup.nC.sup.).sup.1C.sup.nS.sup.].sup.1, wherein C is the compliance tensors, C is an overall micromechanically averaged compliance tensor, f.sup. is a volume fraction of an inclusion phase, C.sup. is an elasticity of the inclusion phase, and S.sup. is an elasticity of the Eshelby inclusion plane.

2. The method of claim 1, further comprising sorting said mineralogy data according to the volume fraction.

3. The method of claim 1, further comprising calibrating said compliance tensor using at least one of the following steps: i) substituting said elastic parameters with a field-tested value and/or experimentally derived values from literature; ii) grouping mechanically similar constituents in said mineralogy data; iii) separating mechanically dissimilar constituents in said mineralogy data; iv) comparing the overall micro-mechanically averaged compliance tensor with laboratory data; v) plotting curves according to the overall micro-mechanically averaged compliance tensor and iteratively adjusting unknown parameters to match predicted curves to one or more control points; and vi) interpreting discrepancies between the overall micro-mechanically averaged compliance tensor and one or more control points.

4. The method of claim 1, wherein the Eshelby's tensor is represented by the following matrix: $\left[\begin{array}{cccccc}{\mathrm{Qa}}_1^2{I}_{11}+{\mathrm{RI}}_1& {\mathrm{Qa}}_2^2{I}_{12}-{\mathrm{RI}}_1& {\mathrm{Qa}}_3^2{I}_{13}-{\mathrm{RI}}_1& 0& 0& 0\\ {\mathrm{Qa}}_1^2{I}_{21}-{\mathrm{RI}}_2& {\mathrm{Qa}}_2^2{I}_{22}+{\mathrm{RI}}_2& {\mathrm{Qa}}_3^2{I}_{23}-{\mathrm{RI}}_2& 0& 0& 0\\ {\mathrm{Qa}}_1^2{I}_{31}-{\mathrm{RI}}_3& {\mathrm{Qa}}_2^2{I}_{32}-{\mathrm{RI}}_3& {\mathrm{Qa}}_3^2{I}_{33}+{\mathrm{RI}}_3& 0& 0& 0\\ 0& 0& 0& \frac{Q}{2}\left(a_1^2+a_2^2\right)I_{12}+\frac{R}{2}\left(I_1+I_2\right)& 0& 0\\ 0& 0& 0& 0& \frac{Q}{2}\left(a_2^2+a_3^2\right)I_{23}+\frac{R}{2}\left(I_2+I_3\right)& 0\\ 0& 0& 0& 0& 0& \frac{Q}{2}\left(a_3^2+a_1^2\right)I_{31}+\frac{R}{2}\left(I_3+I_1\right)\end{array}\right]$ wherein a1, a2, a3 are semi-axes and a1>a2>a3; Q and R are constant where Q=(1) and R=(12)/8(1); the I terms are $I_1=\frac{4a_1{a}_2{a}_3}{\left(a_1^2-a_2^2\right)\sqrt{a_1^2-a_3^2}}\left[{}_0^{}\frac{1}{\sqrt{1-k^2{\sin }^2}}d-{}_0^{}\sqrt{1-k^2{\sin }^2}d\right]$ $I_3=\frac{4a_1{a}_2{a}_3}{\left(a_2^2-a_3^2\right)\sqrt{a_1^2-a_3^2}}\left[\frac{a_2\sqrt{a_1^2-a_3^2}}{a_1{a}_3}-{}_0^{}\sqrt{1-k^2{\sin }^2}d\right]$ $I_2=4-I_1-I_3$ $I_{\mathrm{ij}}=\frac{I_i-I_j}{3\left(a_i^2-a_j^2\right)},ij$ $I_{\mathrm{ii}}=\frac{4}{3a_i^2}-I_{\mathrm{ij}}-I_{\mathrm{ik}},ijk$ wherein the parameter and k are represented as $={\sin }^{-1}\sqrt{1-\frac{a_3^2}{a_1^2}},k=\sqrt{\frac{a_1^2-a_2^2}{a_1^2-a_3^2}}.$

5. The method of claim 1, wherein the elastic properties in the continuum micromechanical elastic properties modeling are imported from experimentally proven values.

6. The method of claim 1, wherein the initial assumption C.sup.0=C(,v) and S.sup.0=S (v), where , v are elastic constants.

7. The method of claim 1, wherein performing continuum micromechanical elastic properties modeling is terminated when norm(C.sup.n+1C.sup.n) is smaller than a predetermined critical value.

8. The method of claim 1, wherein performing continuum micromechanical elastic properties modeling is performed for more than 3 iterations and fewer than 8 iterations.

9. The method of claim 1, wherein the Eshelby's tensor is represented by the following matrix: $\left[\begin{array}{cccccc}{\mathrm{Qa}}_1^2{I}_{11}+{\mathrm{RI}}_1& {\mathrm{Qa}}_2^2{I}_{12}-{\mathrm{RI}}_1& {\mathrm{Qa}}_3^2{I}_{13}-{\mathrm{RI}}_1& 0& 0& 0\\ {\mathrm{Qa}}_1^2{I}_{21}-{\mathrm{RI}}_2& {\mathrm{Qa}}_2^2{I}_{22}+{\mathrm{RI}}_2& {\mathrm{Qa}}_3^2{I}_{23}-{\mathrm{RI}}_2& 0& 0& 0\\ {\mathrm{Qa}}_1^2{I}_{31}-{\mathrm{RI}}_3& {\mathrm{Qa}}_2^2{I}_{32}-{\mathrm{RI}}_3& {\mathrm{Qa}}_3^2{I}_{33}+{\mathrm{RI}}_3& 0& 0& 0\\ 0& 0& 0& \frac{Q}{2}\left(a_1^2+a_2^2\right)I_{12}+\frac{R}{2}\left(I_1+I_2\right)& 0& 0\\ 0& 0& 0& 0& \frac{Q}{2}\left(a_2^2+a_3^2\right)I_{23}+\frac{R}{2}\left(I_2+I_3\right)& 0\\ 0& 0& 0& 0& 0& \frac{Q}{2}\left(a_3^2+a_1^2\right)I_{31}+\frac{R}{2}\left(I_3+I_1\right)\end{array}\right]$ wherein a1, a2, a3 are semi-axes and a1>a2>a3; Q and R are constant where Q=(1) and R=(12)/8(1); the I terms are $I_1=\frac{4a_1{a}_2{a}_3}{\left(a_1^2-a_2^2\right)\sqrt{a_1^2-a_3^2}}\left[{}_0^{}\frac{1}{\sqrt{1-k^2{\sin }^2}}d-{}_0^{}\sqrt{1-k^2{\sin }^2}d\right]$ $I_3=\frac{4a_1{a}_2{a}_3}{\left(a_2^2-a_3^2\right)\sqrt{a_1^2-a_3^2}}\left[\frac{a_2\sqrt{a_1^2-a_3^2}}{a_1{a}_3}-{}_0^{}\sqrt{1-k^2{\sin }^2}d\right]$ $I_2=4-I_1-I_3$ $I_{\mathrm{ij}}=\frac{I_i-I_j}{3\left(a_i^2-a_j^2\right)},ij$ $I_{\mathrm{ii}}=\frac{4}{3a_i^2}-I_{\mathrm{ij}}-I_{\mathrm{ik}},ijk$ wherein the parameter and k are represented as $={\sin }^{-1}\sqrt{1-\frac{a_3^2}{a_1^2}},k=\sqrt{\frac{a_1^2-a_2^2}{a_1^2-a_3^2}}.$

10. The method of claim 1, wherein one or more drilling operations can be adjusted in view of overall micro-mechanically averaged compliance tensor.

11. A method of deriving continuous mechanical properties, said method comprising: obtaining mineralogy data from drill cuttings from a wellbore using x-ray diffraction (XRD), wherein said mineralogy data comprises composition and texture parameters by depth, constituent type, volume fraction of constituents, porosity, shape factor, and elastic parameters including Young's modulus, Poisson's ratio, bulk modulus, and Lame's constants; compiling said mineralogy data by volume fraction and designating at least two of said elastic parameters; inverting said elastic parameters into a compliance tensor; converting said shape factor into Eshelby shape tensors; performing continuum micromechanical elastic properties modeling iteratively on the Eshelby shape tensors using Eshelby's inclusion method for heterogeneous composite materials; calculating an overall micro-mechanically averaged compliance tensor based on the continuum micromechanical elastic properties modeling using at least one of the following steps: i) substituting said elastic parameters in step a) with field-tested value or experimentally derived values from literature; ii) grouping mechanically similar constituents in said mineralogy data; iii) separating mechanically dissimilar constituents in said mineralogy data; iv) comparing the overall micro-mechanically averaged compliance tensor with laboratory data; v) plotting curves according to the overall micro-mechanically averaged compliance tensor and iteratively adjusting unknown parameters to match predicted curves to one or more control points; and vi) interpreting discrepancies between the overall micro-mechanically averaged compliance tensor and one or more control points, displaying results as at least one of (1) continuous plots of the overall micro-mechanically averaged compliance vs. depth, (2) conceptual testing of an effect of changing the volume fraction of constituents, or (3) ternary plots of volume fraction and elastic properties, wherein the continuum micromechanical elastic properties modeling is performed with the following equation: $\stackrel{\_}{C}=C+\underset{=1}{\overset{N}{.Math.}}{f}^{}\left(C^{}-C\right){\left(\stackrel{\_}{C}-C^{}\right)}^{-1}{\stackrel{\_}{C}\left[{\left(\stackrel{\_}{C}-C^{}\right)}^{-1}\stackrel{\_}{C}-{\stackrel{\_}{S}}^{}\right]}^{-1}$ using an initial assumption of C and S.sup. to incrementally solve the equation by
C.sup.n1Cf.sup.(C.sup.C) [(C.sup.nC.sup.)C][(C.sup.nC.sup.).sup.1C.sup.nS.sup.].sup.1 wherein C is the compliance tensors, C is an overall micromechanically averaged compliance tensor, f.sup. is a volume fraction of an inclusion phase, C.sup. is an elasticity of the inclusion phase, and S.sup. is an elasticity of the Eshelby inclusion plane.

12. The method of claim 11, wherein the elastic properties in the continuum micromechanical elastic properties modeling are imported from experimentally proven values.

13. The method of claim 11, wherein the initial assumption C.sup.0=C(,v) and S.sup.0=S(v), where , v are elastic constants.

14. The method of claim 11, wherein performing continuum micromechanical elastic properties modeling is terminated when norm(C.sup.n+1C.sup.n) is smaller than a predetermined critical value.

15. A method of deriving continuous mechanical properties, said method comprising: obtaining mineralogy data from drill cuttings from a wellbore using x-ray diffraction (XRD), wherein said mineralogy data comprises composition and texture parameters by depth, constituent type, volume fraction of constituents, porosity, shape factor, and elastic parameters including Young's modulus, Poisson's ratio, bulk modulus, and Lame's constants; compiling said mineralogy data by volume fraction and designating at least two of said elastic parameters; inverting said elastic parameters into a compliance tensor; converting said shape factor into Eshelby shape tensors; performing continuum micromechanical elastic properties modeling iteratively on the Eshelby shape tensors using Eshelby's inclusion method for heterogeneous composite materials; calculating an overall micro-mechanically averaged compliance tensor based on the continuum micromechanical elastic properties modeling; and displaying results as at least one of (1) continuous plots of the overall micro-mechanically averaged compliance vs. depth, (2) conceptual testing of an effect of changing the volume fraction of constituents, or (3) ternary plots of volume fraction and elastic properties, wherein the Eshelby's tensor is represented by the following matrix: $\left[\begin{array}{cccccc}{\mathrm{Qa}}_1^2{I}_{11}+{\mathrm{RI}}_1& {\mathrm{Qa}}_2^2{I}_{12}-{\mathrm{RI}}_1& {\mathrm{Qa}}_3^2{I}_{13}-{\mathrm{RI}}_1& 0& 0& 0\\ {\mathrm{Qa}}_1^2{I}_{21}-{\mathrm{RI}}_2& {\mathrm{Qa}}_2^2{I}_{22}+{\mathrm{RI}}_2& {\mathrm{Qa}}_3^2{I}_{23}-{\mathrm{RI}}_2& 0& 0& 0\\ {\mathrm{Qa}}_1^2{I}_{31}-{\mathrm{RI}}_3& {\mathrm{Qa}}_2^2{I}_{32}-{\mathrm{RI}}_3& {\mathrm{Qa}}_3^2{I}_{33}+{\mathrm{RI}}_3& 0& 0& 0\\ 0& 0& 0& \frac{Q}{2}\left(a_1^2+a_2^2\right)I_{12}+\frac{R}{2}\left(I_1+I_2\right)& 0& 0\\ 0& 0& 0& 0& \frac{Q}{2}\left(a_2^2+a_3^2\right)I_{23}+\frac{R}{2}\left(I_2+I_3\right)& 0\\ 0& 0& 0& 0& 0& \frac{Q}{2}\left(a_3^2+a_1^2\right)I_{31}+\frac{R}{2}\left(I_3+I_1\right)\end{array}\right]$ wherein a1, a2, a3 are semi-axes and a1>a2>a3; Q and R are constant where Q=(1) and R=(12)/8(1); the I terms are $I_1=\frac{4a_1{a}_2{a}_3}{\left(a_1^2-a_2^2\right)\sqrt{a_1^2-a_3^2}}\left[{}_0^{}\frac{1}{\sqrt{1-k^2{\sin }^2}}d-{}_0^{}\sqrt{1-k^2{\sin }^2}d\right]$ $I_3=\frac{4a_1{a}_2{a}_3}{\left(a_2^2-a_3^2\right)\sqrt{a_1^2-a_3^2}}\left[\frac{a_2\sqrt{a_1^2-a_3^2}}{a_1{a}_3}-{}_0^{}\sqrt{1-k^2{\sin }^2}d\right]$ $I_2=4-I_1-I_3$ $I_{\mathrm{ij}}=\frac{I_i-I_j}{3\left(a_i^2-a_j^2\right)},ij$ $I_{\mathrm{ii}}=\frac{4}{3a_i^2}-I_{\mathrm{ij}}-I_{\mathrm{ik}},ijk$ wherein the parameter and k are represented as $={\sin }^{-1}\sqrt{1-\frac{a_3^2}{a_1^2}},k=\sqrt{\frac{a_1^2-a_2^2}{a_1^2-a_3^2}}.$

16. The method of claim 15, further comprising sorting said mineralogy data according to the volume fraction.

17. The method of claim 15, further comprising calibrating said compliance tensor using at least one of the following steps: i) substituting said elastic parameters with a field-tested value and/or experimentally derived values from literature; ii) grouping mechanically similar constituents in said mineralogy data; iii) separating mechanically dissimilar constituents in said mineralogy data; iv) comparing the overall micro-mechanically averaged compliance tensor with laboratory data; v) plotting curves according to the overall micro-mechanically averaged compliance tensor and iteratively adjusting unknown parameters to match predicted curves to one or more control points; and vi) interpreting discrepancies between the overall micro-mechanically averaged compliance tensor and one or more control points.

18. The method of claim 15, wherein the elastic properties in the continuum micromechanical elastic properties modeling are imported from experimentally proven values.

19. The method of claim 15, wherein performing continuum micromechanical elastic properties modeling is terminated when norm(C.sup.n+1C.sup.n) is smaller than a predetermined critical value.

20. The method of claim 15, wherein performing continuum micromechanical elastic properties modeling is performed for more than 3 iterations and fewer than 8 iterations.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a flow chart illustrating the method of this disclosure.

(2) FIG. 2A-B shows the modeling results of Young's modulus comparing Finite Element Method (FEM) and MMEPS solution for a 2-phase and 3-phase composite.

(3) FIG. 3A-B shows the continuous modeling for Young's Modulus along two wellbores.

(4) FIG. 4A-C shows the modeling results of the disclosed method.

(5) FIG. 4A is a ternary plot of claycarbonate+pyritequartz+feldspar.

(6) FIG. 4B is the grouping of data points based on similar parameters.

(7) FIG. 4C shows the linking between elastic rock properties and rock composition and texture.

(8) FIG. 5A-B illustrates the basic assumptions of inhomogeneities and same volume fraction for the phases in the inclusion scheme.

DETAILED DESCRIPTION

(9) Mechanical rock properties are fundamental input for reservoir stress modeling. Current means to determine mechanical properties from logs is not accurate in shale gas systems. Therefore if a quantitative link between depositional stratigraphy and mechanical stratigraphy can be established to derive continuous mechanical properties, well-drilling can be more efficiently accomplished.

(10) Non-conventional reservoirs refers to reservoirs that are not produced through conventional oil production techniques and may include oil sand, shale, tight-gas sands, coalbed methane, heavy oil, tar sands, etc. These reservoirs generally cannot be simulated through conventional techniques and data. Therefore, the strategy for non-conventional reservoirs could be to constrain the uncertainty in dynamic stress measurements by understanding the relationship between observed deformation behavior and lithology.

(11) The disclosed method therefore focuses on implementing the micromechanical solutions for the overall elastic moduli and its application to mineralogical data from the well, especially as applied to wells located in or near non-conventional reservoirs.

(12) In one embodiment, this disclosure describes a method of deriving continuous mechanical properties in a subterranean formation having at least one wellbore therein, said method comprising: a) obtaining mineralogy data around said wellbore, wherein said mineralogy data comprises composition and texture parameters by depth, constituent type, volume fraction of constituents, porosity, shape factor, and elastic parameters including Young's modulus, Poisson's ratio, bulk modulus, and Lame's constants; b) compiling said mineralogy data by volume fraction and designating at least two of said constituent elasticity parameters; c) inverting said elastic parameters into a compliance tensor [C]; d) converting said shape factor parameters into Eshelby shape tensors; e) performing continuum micromechanical elastic properties modeling iteratively using Eshelby's inclusion method for heterogeneous composite materials; f) calculating overall micro-mechanically averaged compliance tensors; and g) printing or displaying results as at least one of (1) continuous plots of mechanical properties vs. depth, (2) conceptual testing of the effect of changing constituent volume fractions, or (3) ternary plots of volume fraction and elastic properties.

(13) The results of the compliance tensor can be further calibrated by using at least one of the following steps: i) substituting said elasticity parameters in step a) with field-tested value or experimentally derived values from literature; ii) grouping mechanically similar constituents in said mineralogy data; iii) separating mechanically dissimilar constituents in said mineralogy data; iv) comparing the calculated results with laboratory data; v) plotting curves according to the calculated results and iteratively adjusting unknown parameters to match predicted curves to control points; and vi) interpreting discrepancies between the calculated results and control points.

(14) The MicroMechanical Elastic Properties Solver (MMEPS) tool has been developed to perform the approach. MMEPS reads in column-based mineralogical data, typically each data point will include a depth and the volume fraction of each constituent. The user then has the option to modify a data file for each potential constituent that includes its parameters, or for some common constituents, default parameters from the literature can be used.

(15) The program then creates a data structure linking the input mineralogy data to the parameter definitions. Each data point (either sorted by depth or sample name) has its constituents sorted by volume fraction and the software builds compiled data structures depending on if a matrix constituent (maximum volume %) needs to be defined, as is the case for the self-consistent scheme.

(16) A separate micromechanical Matlab script is modified and implemented to invert the elastic input parameters into a compliance tensor [C] and to convert the shape factors into the Eshelby shape tensor [S]. These are also inputs to choose one of two multi-mineral micromechanical solutions (self-consistent or a multi-component e.g., polycrystalline scheme, after Nemat-Nasser and Hori, 1993). Other schemes may be added if necessary. MMEPS assembles all the data structures and then runs a script that solves the multi-mineral micromechanical equations:
C=C+.sub.=1.sup.Nf.sup.(C.sup.C)(CC.sup.).sup.1C[(CC.sup.).sup.1CS.sup.].sup.1

(17) The above equation is solved using an iterative numerical scheme using an initial guess of C and S.sup., for example C.sup.0=C(,v) and S.sup.0=S(v), where ,v are elastic constants. The incremental form of the micromechanical equation is solved by
C.sup.n+1=Cf.sup.(C.sup.C)[(C.sup.nC.sup.)C][(C.sup.nC.sup.).sup.1C.sup.nS.sup.].sup.1

(18) From the new matrix C.sup.n+1 the updated elastic parameter .sup.n+1 is determined and used to update S as S.sup.n+1(.sup.n+1). The iteration ends when norm(C.sup.n+1C.sup.n)<critical value, e.g. 10-.sup.6. The scheme generally converges in about 47 iterations with an error of <10%.

(19) The program then outputs the overall micromechanically averaged compliance tensor [C] for each data point. It has options to export compliance or to solve for other stiffness moduli such as Young's modulus, shear modulus, etc.

(20) A plotting script is run to create e.g., one or more graphs, such as (1) depth vs. elastic property (in the case of log-based data); (2) volume fraction vs. elastic property (in the case of conceptual modeling of compositional mixing effects); (3) a ternary plot comprised of a color contoured equilateral triangle, where the color is the elastic property and the points of the triangle are volume fractions for sedimentalogically sorted end-members (e.g., total clay-carbonate-siliciclastics).

(21) The solution results can be calibrated by at least one of the following: (1) inputting the best-known constituent parameters using field-tested values or experimentally derived values from the literature; (2) grouping mechanically similar constituents in the original mineralogy data file to reduce the number of independent variables that must be calibrated; (3) separating mechanically dissimilar constituents in the original data file if necessary (e.g., separating quartz grains from quartz overgrowth cement if thin section data is also available); (4) comparing the solver predictions with laboratory data for selected intervals to lock-in constituent parameters and reduce the number of unknowns; (5) iteratively adjusting remaining unknown parameters to match the continuous predicted curves to control points; and (6) using geologic experience to interpret remaining discrepancies (e.g., a certain zone has undergone a unique cementation process or thin micro-laminations are introducing anisotropic effects).

(22) The technique offers an alternative way to derive continuous mechanical elastic properties in the wellbore. Previously, the only approach readily available to derive continuous mechanical properties was based on dynamic moduli measured with sonic velocity-based logging tools or using empirical transforms for gamma ray and density logs. The sonic-approach requires conversion using a dynamic-static transform, which may be difficult to attain in horizontal wells. The gamma ray and density log approach is typically poorly calibrated and yields extreme non-unique solutions.

(23) The developed micromechanical approach in this disclosure does not require a dynamic-static transform, and is instead based on geologic description of the composition and texture, and when applied to drilling cuttings data (e.g., XRD), can be used continuously along horizontal wells, where there are often difficulties running the logging tools.

(24) Referring to FIG. 1, which is a flow chart of the disclosed method, as exemplified in a script for MATLAB. However, the same concept can be applied to any software that process geomechanical information.

(25) In step 101, the first step of the method is to obtain and import the mineralogy data around a particular wellbore of interest. The mineralogy data is preferably column-based and typically contains information such as the depth of the location where the sample was taken, and the volume fraction of each constituent. For example, for a sample that was taken at 2,500 feet below surface and having quartz, calcite, kerogen with 50% porosity, the data would read 2500-0.25-0.5-0.2-0.5.

(26) In step 103, the user may modify and compile the data file for each potential constituent that includes certain useful parameters. In some embodiments, default parameters from literature or other verified sources may be used for certain common constituents. For example, the parameters for a quartz grain may read: [QuartzYoung's modulus=30 GPa; Poisson's ratio=0.15; shape factor a=b=c=1 (i.e. sphere)]. The program then creates a data structure linking the input mineralogy data in step 101 to the parameter definitions. The data file can be sorted however the user wishes, and each data point has its constituents sorted by volume fraction. The software thus builds a compiled data structure.

(27) In step 105, with this modified data file, the software then inverts the elasticity parameters into an initial compliance tensor [C]. The typical elastic inputs are Young's modulus E and Poisson's ratio v. For the isotropic case there are only two parameters and the elastic tensor can be written as
C.sub.ijkl=.sub.ij.sub.kl+(.sub.ij.sub.kl+.sub.il.sub.jl)

(28) Where and are Lame constants, and =K2/3, and K is the Bulk modulus, K=E/3(12), so that for example, component

(29) $C_{11}=\frac{E}{1+v}+\frac{1-v}{1-2v}.$

(30) In step 107, the shape factor is converted to Eshelby 66 shape tensor [S]. This is required to apply the Eshelby inclusion as discussed above. For the most general cases, the form of Eshelby's shape tensor for an ellipsoidal inclusion with semi-axes a.sub.1>a.sub.2>a.sub.3 has the form as

(31) $\left[\begin{array}{cccccc}{\mathrm{Qa}}_1^2{I}_{11}+{\mathrm{RI}}_1& {\mathrm{Qa}}_2^2{I}_{12}-{\mathrm{RI}}_1& {\mathrm{Qa}}_3^2{I}_{13}-{\mathrm{RI}}_1& 0& 0& 0\\ {\mathrm{Qa}}_1^2{I}_{21}-{\mathrm{RI}}_2& {\mathrm{Qa}}_2^2{I}_{22}+{\mathrm{RI}}_2& {\mathrm{Qa}}_3^2{I}_{23}-{\mathrm{RI}}_2& 0& 0& 0\\ {\mathrm{Qa}}_1^2{I}_{31}-{\mathrm{RI}}_3& {\mathrm{Qa}}_2^2{I}_{32}-{\mathrm{RI}}_3& {\mathrm{Qa}}_3^2{I}_{33}+{\mathrm{RI}}_3& 0& 0& 0\\ 0& 0& 0& \frac{Q}{2}\left(a_1^2+a_2^2\right)I_{12}+\frac{R}{2}\left(I_1+I_2\right)& 0& 0\\ 0& 0& 0& 0& \frac{Q}{2}\left(a_2^2+a_3^2\right)I_{23}+\frac{R}{2}\left(I_2+I_3\right)& 0\\ 0& 0& 0& 0& 0& \frac{Q}{2}\left(a_3^2+a_1^2\right)I_{31}+\frac{R}{2}\left(I_3+I_1\right)\end{array}\right]$
In which Q and R are constants as Q=3/8(1), and R=(12)/8(1). The I terms are given as follows

(32) $I_1=\frac{4a_1{a}_2{a}_3}{\left(a_1^2-a_2^2\right)\sqrt{a_1^2-a_3^2}}\left[\underset{0}{\stackrel{}{}}\frac{1}{\sqrt{1-k^2{\sin }^2}}d-\underset{0}{\stackrel{}{}}\sqrt{1-k^2{\sin }^2}d\right]$$I_3=\frac{4a_1{a}_2{a}_3}{\left(a_2^2-a_3^2\right)\sqrt{a_1^2-a_3^2}}\left[\frac{a_2\sqrt{a_1^2-a_3^2}}{a_1{a}_3}-\underset{0}{\stackrel{}{}}\sqrt{1-k^2{\sin }^2}d\right]$$I_2=4-I_1-I_3$$I_{\mathrm{ij}}=\frac{I_i-I_j}{3\left(a_i^2-a_j^2\right)},ij$$I_{\mathrm{ii}}=\frac{4}{3a_i^2}-I_{\mathrm{ij}}-I_{\mathrm{ik}},ijk$
In which the parameter and k have the form as

(33) $={\sin }^{-1}\sqrt{1-\frac{a_3^2}{a_1^2}},k=\sqrt{\frac{a_1^2-a_2^2}{a_1^2-a_3^2}}$
For example, one special geometric case is an oblate spheroid with axes (a.sub.1=a.sub.2>a.sub.3)

(34) $I_1=I_2=\frac{2a_1^2{a}_3}{{\left(a_1^2-a_3^2\right)}^{3/2}}\left[{\cos }^{-1}\frac{a_3}{a_1}-\frac{a_3}{a_1}{\left(1-\frac{a_3^2}{a_1^2}\right)}^{1/2}\right]$$I_{12}{I}_{21}=/3a_1^2-I_{13}/4$

(35) In step 109, the user elects either the self-consistent or multi-component schemes, which means the elastic properties of the dominant constituents are either predetermined or iteratively averaged and corrected. Either scheme provides effective medium approximations based on Eshelby's elasticity solution for an inhomogeneity embedded in an infinite medium. The multiphase composite method (Nemat-Nasser and Hori, 1993) is a generalization of the double-inclusion method and assumes an ellipsoidal volume with inclusions and its overall properties, all embedded in an infinite elastic domain with a known elasticity. Interactions are approximated and it is not necessary to have a matrix phase. If a matrix phase is required, then the infinite elasticity can be set in one of two ways. If the infinite elasticity is set to the matrix constituent C, then the result is the Mori-Tanaka model. If the infinite elasticity C of the medium is selected as equal to the overall material properties of the composite, the self-consistent method is modeled and numerical iteration is required to determine the converged overall stiffness from an initial estimate.

(36) The basic assumptions are inhomogeneities and same volume fraction for the phases, as illustrated in FIG. 5A-B.

(37) The user has the option to import experimentally proven values for the elastic properties to be used as the inclusion properties. Each constituent is assumed to be an ellipsoidal elastic inclusion with a shape factor determined by the axes (a.sub.1, a.sub.2, a.sub.3) that are input into the equation for Eshelby's 66 shape tensor [S]. The shape factor can be determined from thin section analysis or other methods.

(38) The internal stresses and strains are solved for within the heterogeneous volume depending on whether boundary strain or stress is applied. The convention of Nemat-Nasser and Hori (1999) is implemented in the software, however other equivalent forms of the homogenization can be used. For example, consider the self-consistent scheme and that a composite is subjected to the displacement boundary condition:
u.sup.0=.sup.0x
Each phase is considered as a single inhomogeneity embedded in the effect medium:

(39) $A_r={\stackrel{\_}{T}}_r\mathrm{and}{A}_0=\frac{1}{c_0}\left(I-\underset{r=1}{\overset{N}{.Math.}}{c}_r{A}_r\right)$
The localization tensor is then obtained from the effective elastic properties of the medium:
T.sub.r=[I+SM(L.sub.rL)].sup.1
where L=.sub.r=0.sup.Nc.sub.rL.sub.rA.sub.r
Consider the composite being subjected to prescribed tractions:
t.sup.0=.sup.0n
At the boundary, the strain is expressed using the Hooke's law:
.sup.0=M.sup.0.fwdarw..sup.0=M
T.sub.r.sup.0=T.sub.rM.fwdarw..sub.r=T.sub.rM
By applying the Hooke's law again:
.sub.r=L.sub.rT.sub.rM
The stress concentration tensor is therefore written as:

(40) $B_r=L_r{\stackrel{\_}{T}}_r\stackrel{\_}{M}\mathrm{and}\text{:}{B}_0=\frac{1}{c_0}\left(I-\underset{r=1}{\overset{N}{.Math.}}{c}_r{B}_r\right)$
Finally, the effective stiffness tensor is obtained from the expression of the concentration tensors because:

(41) $\stackrel{\_}{M}=\underset{r=0}{\overset{N}{.Math.}}{c}_r{M}_r{B}_r$$\mathrm{and}\stackrel{\_}{M}\stackrel{\_}{L}=I$

(42) In summary, for self-consistent scheme, average strain of the composite is:
.sub.Ref=.sup.0= where L.sub.Ref=L

(43) Instead of dealing with phases of the above-mentioned heterogeneous solid, a technique is used to approximate an equivalent homogeneous solid with uniform properties. The difference between the properties of the inclusions and matrix are addressed by introducing the micromechanical concept of eigenstrains or eigenstresses due to either far-field stresses or strains, respectively. The eigenstrain * is the suitable strain field introduced in the domain such that the equivalent homogeneous solid has the same strain and stress fields as the actual heterogeneous solid under applied tractions or displacements, whichever may be the case. The strain and stress fields over the matrix M and inclusion are thus given as a function of the uniform (e.g., .sup.0) and disturbed (.sup.d) states:

(44) 0$\mathrm{.Math.}\left(x\right)={\mathrm{.Math.}}^0+{\mathrm{.Math.}}^d\left(x\right)\mathrm{and}\left(x\right)={}^0+{}^d\left(x\right)$$\left(x\right)=C\text{:}(\mathrm{.Math.}\left(x\right)-{\mathrm{.Math.}}^{*}\left(x\right)=\{\begin{array}{cc}C\text{:}\left({\mathrm{.Math.}}^0+{\mathrm{.Math.}}^d\left(x\right)\right)& \mathrm{in}M\\ C\text{:}\left({\mathrm{.Math.}}^0+{\mathrm{.Math.}}^d\left(x\right)-{\mathrm{.Math.}}^{*}\left(x\right)\right)& \mathrm{in}\end{array}$

(45) The average eigenstrains over the matrix and inclusions are substituted back into exact solutions for the average strain and stress over the volume. The elastic tensors are derivations of the stress and strain fields.

(46) In step 111, the overall elastic moduli are then calculated according to predetermined equations and the overall micromechanically averaged compliance tensor [C] is obtained. The combined effect of all the inclusions among the domain is determined by their respective volume fractions. The constituents with higher volume fractions will be assigned higher weight in the calculation. Here the inclusions are assumed to be evenly distributed in the rock so that the solution is not spatially variant. In other words, the Representative Volume Element or RVE is maintained. The micromechanical solution scheme calculates the stress-strain relationship inside and outside of the inclusion, including at the boundaries, and then the overall elastic moduli are solved for.

(47) In one embodiment, equations are given for the multi-phase composite, Mori-Tanaka, or self-consistent methods, depending on how the matrix stiffness is prescribed, and according the form:
C=C+.sub.=1.sup.Nf.sup.(C.sup.C)(CC.sup.).sup.1C[(CC.sup.).sup.1CS.sup.].sup.1

(48) However, other suitable solution methods may also be used interchangeably or integrated herein. For example, a Hashin-Shtrikman definition may be added so that the results are compatible with other two-component rock physics tools.

(49) Hashin-Shtrikman Bounds.

(50) The Hashin-Shtrikman bounds are the tightest bounds possible from range of composite moduli for a two-phase material. Specifying the volume fraction of the constituent moduli allows the calculation of rigorous upper and lower bounds for the elastic moduli of any composite material. The so-called Hashin-Shtrikman bounds for the bulk, K, and shear moduli is given by:

(51) $K_{\mathrm{HS}}^{}=K_2+\frac{}{{\left(K_1-K_2\right)}^{-1}+\frac{1-}{K_2+4{}_2/3}}$${}_{\mathrm{HS}}^{}={}_2+\frac{}{{\left({}_1-{}_2\right)}^{-1}+\frac{2\left(1-\right)\left(K_2+2{}_2\right)}{5{}_2\left(K_2+4{}_2/3\right)}}$

(52) The upper bound is computed when K.sub.2>K.sub.1. The lower bound is computed by interchanging the indices in the equations. For the case of a solid-fluid mixture, K.sub.2 is K.sub.S, the bulk modulus of the solid component, and K.sub.1 is K.sub.f, the bulk modulus of the fluid component.

(53) In step 113, the user has the option to plot the results as depth vs. elastic properties, volume fraction vs. elastic properties, or as a ternary plot, depending on user's need. Depth vs. elastic property plots is typically used for log-based data; the volume fraction vs. elastic property plot is typically used for conceptual modeling of compositional mixing effects; the ternary plot usually comprises a color-contoured equilateral triangle, where the color is the elastic property and the points of the triangle are volume fractions for sedimentalogically sorted end-members, such as total clay-carbonate-siliciclastics as the three points.

(54) Alternatively or in addition, the final results can be recorded and stored and/or sent to another reservoir modeling program. Ultimately, the results are used to make and implement decisions on an optimal development plan for a reservoir and/or to predict production levels and the like.

(55) Steps 115-125 are several ways to calibrate the modeling results in steps 111. These steps may be used alone or two or more of these steps may be used together to obtain more accurate results.

(56) In step 115, the inputs for parameters of the constituents are replaced with field-tested values or experimentally derived values. The results are then compared with the results from step 111 to see if further adjustment to the method is necessary.

(57) In step 117, the inputs for parameters of the constituents are grouped based on mechanical similarities. The similar ones are grouped together to reduce the number of independent variables that need to be calibrated. By similar what is meant is the type of the constituent is the similar, such as grouping two types of clay, and that grouping the parameters results in no more than 30% variation, preferably no more than 25%, and most preferably no more than 20% variation from a solution that does not group similar constituents.

(58) In step 119, the inputs for parameters of the constituents are separated if they are mechanically dissimilar. For example, thin section petrographic characterization may be available to distinguish quartz grains from quartz overgrowth cement. By dissimilar what is meant is the type of the constituent is believed to be mechanically distinct, the volume fraction of each constituent may be different, and differentiating the components yields results with values are at least 30% different from leaving the components grouped. For example, matrix supporting clay should be distinguished from grain-coating cement occurring from diagenetic alteration.

(59) In step 121, the modeling results are compared with laboratory data for selected intervals (i.e. those with same parameters), so as to lock-in the constituent parameters and reduce the number of unknowns. For example, a laboratory sample with 90% quartz volume of and 10% porosity can be used to set the constituent parameters for quartz by assuming the porous volume has zero stiffness.

(60) In step 123, the modeling results from step 111 are iteratively adjusted for the remaining unknown parameters to match the continuous prediction curves to control points. For example, if the user can use step 121 to determine the parameters for quartz and porosity phases, then a third sample containing 80% quartz, 10% porosity, and 10% clay could be used to solve for the unknown stiffness of the clay phase.

(61) In step 125, already existing geologic experience can be helpful in interpreting the remaining discrepancies between the calculated results and the actual observation. For example, for certain zones that undergo a unique cementation process or thin micro-lamination could introduce anisotropic effects to the formation. Such experience can be added to the software so that the discrepancies can be accounted for and corrected.

(62) The final results can be displayed in any suitable manner, e.g., on a monitor, projected as a hologram, or printed for use. Alternatively, the results can be stored in memory, and/or forwarded to another program for further use, e.g., in reservoir modeling, planning, and the like.

(63) The following descriptions are intended to be illustrative only, and not unduly limit the scope of the appended claims.

Test 1

(64) The hypothetical variations of a two-phase rock was tested for validating the model. Under the FE (finite element) Model, a triaxial loading of heterogeneous plug was assumed, where the Young's moduli (E) varied for each element. For this test, the results of the finite element method are assumed to represent the actual overall stiffness of the rock as a function of compositional variation. The overall elasticity of the same rock formation was also modeled by using Micromechanical Elastic Properties Solver (MMEPS) of this disclosure. The results are shown in FIG. 2A-B.

(65) Results show that Young's modulus decreases non-linearly with addition of the less stiff volume constituent. For the two-phase model (FIG. 2A) the MMEPS solution using the self-consistent scheme (E_ SC) closely matches the finite element result, whereas the multiphase-component model (E_MPC) has greater error. The MMEPS implementation allowed for very rapid testing of an appropriate micromechanical solution for this configuration, with generation of the input file taking only a minute and computational run-time of only a few seconds. In contrast, the finite element solution set took several hours to build, run, and compile results.

(66) Similarly, FIG. 2B shows the robustness of the MMEPS model to add a third volume phase, porosity, shown for 5% and 10% volume fractions, and using the self-consistent model. Setting up and running the model using MMEPS took only a few seconds, whereas adding porosity to the same finite element would have increased model setup time and run time significantly.

Test 2

(67) The MMEPS modeling approach was further tested with two wells in the Muskwa formation, Horn River Basin Canada. The method was as described above. FIG. 3A-B shows the depth vs. Young's Modulus chart, where the actual core data taken at different depth was compared with different volume fractions of clay, calcite and quartz. FIG. 3A-B shows that the modeling results closely follow the actual core data.

Test 3

(68) The third test shows the ability to link potentially proprietary MMEPS mechanical solutions to other supporting data, for example micro-facies characterization from published papers. An exemplary ternary plot based on the MMEPS approach is shown in FIG. 4A-C, where the mineralogy data was taken from the Lewis Shale in Wyoming.

(69) FIG. 4A is a ternary plot showing volume fraction total claycarbonate+pyritequartz+feldspar. The ternary plot shows data from a ConocoPhillips well overlaid on digitized results from a publication by Almon et al. (2005), who provide the mineral categorization scheme. All the mineralogy data is overlaid on an MMEPS elastic properties solution that is calibrated to unpublished laboratory data. This is an easy way to quickly identify the elastic properties of the formation.

(70) FIG. 4B is an exemplary grouping of data points based on similar properties for the data (Almon et al. 2005). For example, Group 1 is likely to show similar properties as massive organic mudstones; Group 2 is likely to show similar properties as organic laminated shales; Group 3 is likely to show similar properties as calcareous laminated shales; Group 4 is likely to show similar properties as organic bioturbated shales; and Group 5 is likely to show similar properties as massive calcareous.

(71) In this grouping, a centroid is iteratively calculated by averaging the values of all included data points, and the variation between each data point and the centroid is then calculated. The overlapped grouping will adopt smaller variation to remove the overlap. This grouping further simplifies the modeling process, because instead of 100 data points, now only 5 groups of similar properties are considered, and which can be linked back to the ternary plot. Accuracy may be sacrificed for simplification, but an optimal balance can be achieved by adjusting the grouping criteria.

(72) FIG. 4C shows the linking between elastic rock properties and rock composition and texture (Almon et al. 2005). This compiled data is not only useful for current wellbore, but can serve as additional references for future wellbores that have similar properties in terms of composition, depth and elastic properties.

(73) The minimum computer requirements are programming language software that can incorporate the algorithms. Because the equations include a numerical iterative scheme, software capable of programming and running iterative loops must be used. A robust graphical plotting tool is necessary to display multiple solutions, such as continuous logs or the properties ternary plot. The present implementation was performed in MATLAB, though other programming environments such as FORTRAN could be used.

(74) The following references are incorporated by reference in their entirety for all purposes:

(75) WO2009108432 Rock physics model for simulating seismic response in layered fractured rocks.

(76) US20140373616 Mechanical characterization of core samples.

(77) US20140352949 Integrating rock ductility with fracture propagation mechanics for hydraulic fracture design.

(78) Voigt, W. (1889), Uber die Beziehung zwischen den beiden Elastizitatskonstanten isotroper Korper, Wied. Ann. Physik, Vol. 38, 573-587.

(79) Budianksy, B. (1965), On the elastic moduli of some heterogeneous materials, J. Mech. Phys. Solids, Vol. 13, 223-227.

(80) Hashin, Z. (1962), The elastic moduli of heterogeneous materials, ASME J. Appl. Mech., Vol. 29, 143-150.

(81) Hashin, Z. and Shtrikman, S. (1963), A Variational approach to the theory of the elastic behavior of multiphase materials, J. Mech. Phys. Solids, Vol. 11, 127-140.

(82) Hashin, Z. (1968), Assessment of the self-consistent scheme approximation, J. Compos. Mater., Vol. 2, 284-300.

(83) Hill, R. (1965), A self-consistent mechanics of composite materials, J. Mech. Phys. Solids, Vol. 13, 213-222.

(84) Kroner, E. (1958), Berechnung der elastischen Konstanten des Vielkristalls aus den Konstanten des Einkristalls, Z. Phys. Vol. 151, 504-518.

(85) Mori, T., and Tanaka, K. (1973), Average stress in matrix and average elastic energy of materials with misfitting inclusions, Acta Met., Vol. 21, 571-574.

(86) Nemat-Nasser, S., and Hori, M. (1999). Micromechanics: overall properties of heterogeneous materials. Elsevier Science B.V., Amsterdam, The Netherlands, Second edition, 786 pages.