Unconventional resource development has become an exercise in determining where horizontal wells should be placed to maximize hydro-fracture efficiency. A large part of optimizing how the resource is developed requires understanding the reservoir from both lithologic and rock properties perspectives and correlating these with production data to infer frac efficiency. Given the variety of unconventional plays in North America, it is important to understand the differences in rock properties so as to effectively map them seismically. Figure 1 shows backscattered electron images displaying the diversity in microstructure beyond the differences in lithology.

Fig. 1
Figure 1. Backscattered electron image of 9 different North American unconventional plays (Curtis et. al 2010).

Understanding the variations in terms of mineralogical composition and microstructure requires an in depth model. With such a model, templates can be generated for interpretation of seismic attributes.

The following is an examination of two distinct rock physics models and how they can be employed to generate seismic interpretation templates. These templates can then be used in conjunction with inverted volumes of seismic attributes to map lithologic and rock property variations. In addition, from the rock physics trends, AVO models are synthesized to evaluate potential offset dependent amplitude responses.

Two different rock physics models

The process begins with a comparison of two distinct rock physics models. Each model will be reviewed, referencing previous important works followed by suggesting potential application to unconventional reservoirs. The first is a granular model referred to as the modified Hertz-Mindlin Hashin-Shtrikman (HMHS) while the second, referred to as the Non-Interacting Approximation (NIA) is an inclusion based model which describes the number of pores and their shapes as opposed to the grain to grain contacts.

The modified Hertz-Mindlin Hashin-Shtrikman (HMHS)

This model is a granular representation of a rock in that it assumes a random packing of spheres under hydrostatic compression. Contact between spheres is used to estimate bulk and shear moduli in terms of normal and tangential stiffness. The expressions for the dry bulk and shear moduli, known as the Hertz-Mindlin estimates, are given by (Dvorkin and Nur, 1996)

Formula 1
Formula 2

where n is the coordination number (number of grain contacts), φ is the amount of porosity, ξ is a grain frictional coefficient (more on this later), and R is the grain radius. The normal and tangential stiffnesses are

Formula 3
Formula 4

where G and v are the shear modulus and Poisson’s ratio of the sphere. The variable a is the radius of the contact area between the two spheres and an estimate is given by Hertz,1882 (Duffaut et. al, 2010),

Formula 5

where FN is the confining force acting between the two spheres, given by

Formula 6

and P is the effective hydrostatic confining pressure.

The scalar value of ξ varies between 0 and 1 where 0 implies no friction between grains and 1 implies infinite friction. There are models that describe the variation of the parameter ξ. Duffaut et. al (2010) and Bachrach and Avseth (2008) studied the impact of the parameter and proposed non-linear and linear expressions, respectively, to model frictional variations. Bachrach and Avseth (2008) also investigate the effects of grain shape (round versus angular) with respect to the contact radius and its influence on the bulk and shear moduli.

The Hertz-Mindlin estimate of K and G is used to describe a rock at the critical porosity end member. This end member is assumed to represent the starting point of rock formation. Any porosity greater than the critical porosity will result in grains in suspension; not frame supported (Nur et. al, 1998). The Hashin- Shtrikman bounds are used to connect the critical porosity and the solid, zero porosity end member to create upper and lower limits for porosity values between the end members. This tren is schematically shown in figure 1, along with the associated LMR representation. These upper and lower bounds have been used to describe cementing and sorting trends (Avseth et. al, 2010).

Fig. 2
Figure 2. Schematic representation of Hashin-Shtrikman bound with critical porosity and mineral end members. From the bulk and shear modulus (not shown) expressions, LMR, or any other seismic attribute, templates can be constructed.

The Hashin-Shtrikman bounds describe lower and upper limits to n-phase mixtures. The Hashin-Shtrikman bounds have the form (Mavko et. al, 2009)

Formula 7
Formula 8


Formula 9

where Kn and Gn are the bulk and shear moduli of each phase, fn is the fractional quantity of each phase and G is the maximum (upper bound) or minimum (lower bound) shear modulus of all the phases.

In the modified HMHS, which incorporates the critical porosity concept, the bounds consist of two phases (n=2): the Hertz- Mindlin defined critical porosity and the mineral, zero porosity end member.

The Non-Interacting Approximation (NIA)

The second model, referred to as the non-interacting approximation (NIA), has been explored (Kachanov 1992, Kachanov et. al, 1994 and Shafiro and Kachanov, 1996) in great detail. This approach constructs the elastic potential of a solid with cavities (pores of various shapes including cracks) for both isotropic and anisotropic scenarios. The basic form of the NIA expression is a sum of individual compliances of components comprising the rock (Vernik and Kachanov 2010) and is illustrated in figure 3.

Formula 10

Note that the expression is linear in compliance and differs from Husdon’s approach in that Hudson deals with stiffness and has a linearization approximation.

Fig. 3
Figure 3. Schematic representation of NIA model. Composite material is a function of solid and various cavities.

The elastic potential is

Formula 11

where σ is stress and e is the strain. Hooke’s Law is invoked to describe the strain on the object as

Formula 12

where Ssolid is the compliance of the solid material and Δe is the additional strain due to inclusions defined by

Formula 13

where Sincl is the so called cavity compliance tensor. For each cavity of interest, the cavity compliance tensor Sincl must be calculated. The total compliance of the material is the sum of the solid plus the cavity compliance for each pore shape being described. The elastic potential can then be expressed as

Formula 14

Once the equation has been set up with appropriate cavity compliance tensors for each pore shape, the moduli of interest can be determined from the stress strain relation. The compliance matrix Sijkl can be of any form and is not exclusively appropriate for isotropic materials. In this instance, the analysis will be restricted to the isotropic scenarios (a subsequent paper will delve into the anisotropic models).

To account for the interaction between cavities, Kachanov proposes using Mori-Tanaka’s scheme (1973) where the average stress of the effective material is

Formula 15

The expression for σavg now replaces σ in equation 13.

As per Vernik and Kachanov (2010), a spherical pores and crack density expression would be

Formula 16
Formula 17

where v is the Poisson’s ratio of the solid, φ is porosity and n is the crack density. Note that:

  1. there is no stress dependence in this expression
  2. porosity affects cracks, but cracks do not affect porosity.

Model comparison

The following will be an examination of the two different rock physics models introduced. Specifically, the impact of pore shape as described by the NIA model is investigated in an attempt to rationalize differences between the granular and inclusion models. For the trends that follow, the rock is composed of 50% quartz, 25% clay and 25% limestone, representative of unconventional reservoirs.

Figure 4 shows the differences between the granular model and inclusion model with spherical cavities. It is immediately apparent that the NIA model generates much larger values of bulk modulus (K) and shear modulus (G) for various amounts of porosity.

Fig. 4
Figure 4. (Left) Bulk modulus (upper) and Shear modulus (lower) as a function of porosity for granular (HMHS) and inclusion (NIA) models. The black and brown trends have an effective stress of 1 MPa and 10 MPa respectively. (Right) LMR representation of rock physics trends.

The blue dotted line overlies the solid blue line in LMR space. This means that the Vp/Vs ratio of the trend does not change by accounting for pore interactions; it is merely a change in porosity scale.

To rationalize differences between the models, an attempt to match the trends is made by adding pores of different shapes in the inclusion (NIA) model, shown in figure 5. The first attempt is the inclusion of infinitely thin cracks. The impact of cracks is to decrease K and G with increase in crack density, n. As in the previous figure, in addition to the K and G versus porosity / crack density plots, an LMR crossplot is included to demonstrate the expected response in a seismic attribute space. Note that at large porosities the effects of cracks is minimal.

Fig. 5
Figure 5. (Left) Bulk modulus (upper) and Shear modulus (lower) as a function of porosity for inclusion (NIA) model with porosity (blue) and crack density (red). Each blue curve in an increase in crack density (seen vertically) while the red curve is an increase in porosity.

The next pore shape shown in figure 6, is that of penny-shaped cracks. The impact is larger than in the infinitely thin cracks.

Fig. 6
Figure 6. (Left) Bulk modulus (upper) and Shear modulus (lower) as a function of porosity for inclusion (NIA) model with porosity (blue) and penny-shaped crack density (red). Each blue curve in an increase in crack density (seen vertically) while the red curve is an increase in porosity.

At first glance, it is observed that the appropriate pore shape and concentration could be constructed to mimic the granular trends. Though not demonstrated in the figures, it is postulated that the granular trends could be a representation of many different pore shapes ranging from spherical pores, penny-shaped cracks and thin cracks at high porosities, to just thin cracks at low porosities.

Fluid Effects

For fluid-filled pore spaces, Gassmann’s theory for fluid substitution (Smith, 2003) can be used to estimate the bulk and shear moduli for various amounts of porosity and fluids. Gassmann’s fluid substitution is suitable for the NIA and HSHM models (Grechka, 2007). The special case of fluid-filled cracks is shown to have a dependence on the crack aspect ratio as well as the fluid itself (Shafiro and Kachanov, 1997). Figure 7 above shows Gassmann’s results for gas and oil filled trends.

Fig. 7
Figure 7. (Left) LMR crossplot of gas filled rock physics trends. (Right) LMR crossplot of oil filled rock physics trends.

Note that the oil filled low pressure HMHS trends almost completely overly each other.

Figure 8 shows the inclusion based model (NIA) incorporating fluids and displaying some interesting results. Figure 8 shows the differences between gas and oil filled penny-shaped cracks (aspect ratio 0.1).

Fig. 8
Figure 8. (Left) LMR crossplot of gas filled penny-shaped cracks trends. (Right) LMR crossplot of oil filled penny-shaped cracks trends.

For thin cracks the response in LMR space is shown in figure 9.

The left figure shows the gas response while the right shows the oil response. Note that there is a slight increase in λρ as the cracks become filled with incompressible fluids. The presence of fluid strengthens the rock as seen from the parameter λ. The combination of figure 8 and 9 show the impact of crack aspect ratio in conjunction with the presence of fluids.

Fig. 9
Figure 9. (Left) LMR crossplot of gas filled thin cracks trends. (Right) LMR crossplot of oil filled thin cracks trends.
Fig. 10
Figure 10. LMR crossplot of inclusion based model (NIA) for three dual mineral trends.

Multi-mineral models

Potential unconventional templates can be modelled as a mixture of quartz, clay, and carbonate material, each with varying quantities. How to visualize the available permutations with the number of variables (mineralogy, fluid, porosity, pore shape, etc.) becomes a large part of the interpretation process. Comparing the previous LMR crossplots, there is overlap amongst the trends and thus ambiguity in the interpretation. One way to visualize the various lithologic permutations is to construct a two mineral trend demonstrating the variations in mineralogical composition and porosity. Figure 10 illustrates three such dual mineral trends: quartz-clay, quartz-limestone and limestone-clay.

Another method could be to select three minerals and determine a distribution that would encompass the majority of potential outcomes. Figure 11 demonstrates the effects of a three mineral combination; quartz, clay and limestone.

Fig. 11
Figure 11. LMR crossplot of inclusion based model (NIA) for different combination of three mineral trends.


Examples of well logs from three different unconventional resource plays in western Canada are presented with potential rock physics trends that can be used for interpretation. Each of the three examples will have a unique template in addition to a more generic multi-play template. The generic multi-play template is added to demonstrate the variability in each play with respect to a fixed template.

The importance of data integration cannot be overemphasized as using the templates in isolation would leave an interpreter with a myriad of possibilities that would minimize the predictive power of the trends. Without additional constraints, the use of modelled trends is not optimized.

Horn River

Figure 12 shows LMR and Gamma Ray (GR) logs along side two LMR crossplots each superimposed with two different rock physics interpretation templates.

The template on the left assumes a single mineralogical mix of 50% quartz, 25% limestone and 25% clay. The trends are generated using the NIA model. Each blue trend corresponds to increasing amount spherical porosity for a fixed penny-shaped crack density. The red trends correspond to increasing crack density for a fixed amount of spherical porosity. In this case, the LMR reservoir variations are attributed entirely to the concentration of penny-shaped cracks. Decreasing LR values are purely a function of increased crack density. Conversely, the crossplot on the right is the template introduced in figure 9. These NIA trends assume spherical pores and varying amounts of minerals. In this case the low LR values correspond to increasing quartz content. The “correct” interpretation template is likely a combination of these templates as well as other pore shape and mineral combinations.

Fig. 12
Figure 12. LMR and GR logs (left) and crossplots (right). Two interpretation scenarios plotted on the crossplots – left: gas filled penny-shaped cracks with a constant mineral content of 50% Quartz, 25% Limestone and 25% Clay. Right – three NIA dual mineral trends; quartz-clay (black), quartz-limestone (green) and limestone-clay (brown).


Figure 13 shows LMR and GR logs for a Duvernay well, again with two LMR crossplot with two interpretation templates. The carbonate rich section (low GR data approaching 100% limestone trend) corresponds to overlying section.

Fig. 13
Figure 13. LMR and GR logs (left) and crossplots (right). Two interpretation scenarios plotted on the crossplots – left: NIA gas filled thin cracks with a constant mineral content of 50% Quartz, 25% Limestone and 25% Clay (blue lines) and gas-filled spherical pores with thin crack density of 0.1 and varying mineralogy (black lines). Right – three NIA dual mineral trends; quartz-clay (black), quartz-limestone (green) and limestone-clay (brown).

It is noted that the Duvernay interval displays a more linear trend compared to the Horn River data. The log data also manifests itself in a range that can be suitably modeled by the quartz rich templates of figure 11. With that in mind, the blue trend lines correspond to an NIA model with a constant mineral mix of 50% quartz, 25% clay and 25% limestone, with increasing spherical porosity and a constant amount of gas-filled thin cracks. The black trend lines correspond to an NIA quartz-clay mix with increasing spherical porosity and a constant amount thin crack density. The crossplot on the right are the same trends shown in figure 11. Given the narrow range of Vp/Vs in the crossplot, one can make the assertion that this well contains very little lithologic variation (or has other minerals with similar bulk and shear moduli as quartz, clay or limestone) and any variations seen could be a function of pore shape.


Figure 14 shows LMR and GR logs for a Montney well as well as set of two interpretation templates. Again, the right crossplot is the dual mineral mix introduced in figure 10. The crossplot on the left however, has NIA trends composed of 60% quartz, 30% dolomite and 10% clay with spherical porosity and various amounts of penny-shaped cracks. The red trends are gas-filled inclusions while the green trends are fluid-filled. Like the Duvernay, Montney rocks display a linear trend in LMR space, indicative of small Vp/Vs changes within the interval.

Fig. 14
Figure 14. LMR and GR logs (left) and crossplots (right). Two interpretation scenarios plotted on the crossplots – left: gas filled penny-shaped cracks with a constant mineral content of 60% Quartz, 30% Dolomite and 10% Clay (red lines) and thin oil filled cracks (green lines). Right – three NIA dual mineral trends; quartz-clay (black), quartz-limestone (green) and limestone-clay (brown).

AVO Template Generation

The rock physics trends can also be used to generate a range of potential AVO responses that can be used to interpret seismic data without inverting the data. A large number of different scenarios can be modelled. Note that the four examples shown follow the expected behaviour as outlined in Perez 2010. Figure 15 shows the AVO responses expected from a clay rich layer overlying quartz, clay and limestone dominated layer. This may represent the transition from overlying non-reservoir shale to a potential unconventional reservoir. The ability to differentiate between quartz, clay or limestone dominated unconventional targets helps in understanding frac efficiency.

Fig. 15
Figure 15. LMR crossplot (left) and AVO plot (right) showing simple interface representation of potential seismic expressions for various mineral combinations. These scenarios could be representative of reflection coming from non-reservoir shale to reservoir boundary. Open circle represents overlying layer, closed circle underlying layer.

Figure 16 shows AVO and LMR plots of what could be representative of a transition from clay dominated to quartz or limestone dominated reservoir. As demonstrated by figure 10, the largest impact in both LMR and AVO signature is a function of quartz content. This could be representative of intra-Montney reflections and demonstrates the subtlety in mapping lithologic variations with AVO.

Fig. 16
Figure 16. LMR crossplot (left) and AVO plot (right) showing simple interface representation of potential seismic expressions for various mineral combinations. These scenarios could be representative of inter-reservoir reflections coming from clay dominated to quartz (red) and limestone (blue) dominated reservoirs. Open circle represents overlying layer, closed circle underlying layer.

Figure 17 illustrates the effect of porosity and the ability of LMR and AVO signatures to detect such porosity variations. This shows a decreasing near offset amplitude and decreasing gradient with increasing porosity.

Fig. 17
Figure 17. LMR crossplot (left) and AVO plot (right) showing simple interface representation of potential seismic expressions for various mineral combinations. These scenarios could be representative of reflections coming from non-reservoir shale to reservoir with varying amounts of porosity. Open circle represents overlying layer, closed circle underlying layer.

Figure 18 models the ability to detect potential carbonate rich frac barriers. In general, because the frac barrier is a “stronger” material, seismic expression is unsurprisingly a peak. The main control of the zero offset response is the high P-impedance while the lithology variations of the overlying layer control the gradient.

Fig. 18
Figure 18. LMR crossplot (left) and AVO plot (right) showing simple interface representation of potential seismic expressions for various mineral combinations. These scenarios could be representative of reflection coming from reservoir shale to potential high carbonate volume frac barriers. Open circle represents overlying layer, closed circle underlying layer.

Two different rock physics models are introduced and used to interpret unconventional reservoirs from a seismic attribute perspective. The differences between the models are rationalized and potential interpretation templates are provided. Some AVO models are also generated demonstrating potential amplitude variations and the subtle nature of the seismic reflection in unconventional settings. These tools allow for a more in depth interpretation of seismic data.


About the Author(s)

Marco Perez received his B.Sc. at McGill University before completing an M.Sc. in Geophysics at the University of Calgary. He started working at PanCanadian, later EnCana, focusing on AVO, inversion and LMR analysis. After moving to Apache in 2007 Marco has continued to work with advanced geophysical techniques within the Exploration & Production Technology group where he is currently a Senior Staff Geophysicist.


Avseth, P., Mukerji, T., Mavko, G., and Dvorkin, J. (2010). ”Rock-physics diagnostics of depositional texture, diagenetic alterations, and reservoir heterogeneity in high-porosity siliciclastic sediments and rocks — A review of selected models and suggested work flows.” GEOPHYSICS, 75(5), 75A31–75A47.

Bachrach, R. and Avseth, P. (2008). ”Rock physics modeling of unconsolidated sands: Accounting for nonuniform contacts and heterogeneous stress fields in the effective media approximation with applications to hydrocarbon exploration.” GEOPHYSICS, 73(6), E197–E209.

Curtis, M.E., Ambrose, R.J., Sondergeld, C.H. and Rai, C.S., 2010, “Structural Characterization of Gas Shales on the Micro- and Nano-Scales, CUSG/SPE 137693.

Dvorkin, J. and Nur, A. (1996). ”Elasticity of high��porosity sandstones: Theory for two North Sea data sets.” GEOPHYSICS, 61(5), 1363–1370.

Duffaut, K., Landrø, M., and Sollie, R. (2010). ”Using Mindlin theory to model friction- dependent shear modulus in granular media.” GEOPHYSICS, 75(3), E143–E152.

Grechka, V., 1997, Fluid Substitution in porous and fractured solids: the non-interaction approximation and gassmann theory, Int J. Fract., 148, 103-107.

Kachanov, M., 1992, “Effective elastic properties of cracked solids: critical review of some basic concepts”, Appl Mech Rev, 45, 8, 304-335.

Kachanov, M. Tsukrov, I. ad Shafiro, B. 1994, “Effective moduli of solids with cavities of various shapes”, Appl. Mech. Rev. 47, 1, S151-S174.

Mavko, G., T. Mukerji, and J. Dvorkin, 2009, The rock physics handbook, 2nd ed.: Cambridge University Press.

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

Nur, A., Mavko, G., Dvorkin, J., and Galmudi, D. (1998). ”Critical porosity: A key to relating physical properties to porosity in rocks.” The Leading Edge, 17(3), 357–362.

Perez M., Beyond isotropy in AVO and LMR: CSEG RECORDER, 35(7), 37-41

Shafiro, B. and Kachanov, M. 1996, “Materials with fluid-filled pores of various shapes: effective elastic properties and fluid pressure polarization”, Int. J. Solids Structures, 34, 27, 3514-3540.

Smith, T., Sondergeld, C., and Rai, C. (2003). ”Gassmann fluid substitutions: A tutorial.” GEOPHYSICS, 68(2), 430–440.

Vernik, L. and Kachanov, M. (2010). ”Modeling elastic properties of siliciclastic rocks.” GEOPHYSICS, 75(6), E171–E182.


Join the Conversation

Interested in starting, or contributing to a conversation about an article or issue of the RECORDER? Join our CSEG LinkedIn Group.

Share This Article