Rock-physics modeling guided by depositional and burial history in low-to-intermediate-porosity sandstonesRock-physics of well-consolidated sandstones | Geophysics

ABSTRACT

We present a new rock-physics modeling approach to describe the elastic properties of low-to-intermediate-porosity sandstones that incorporates the depositional and burial history of the rock. The studied rocks have been exposed to complex burial and diagenetic history and show great variability in rock texture and reservoir properties. Our approach combines granular medium contact theory with inclusion-based models to build rock-physics templates that take into account the complex burial history of the rock. These models are used to describe well log data from tight gas sandstone reservoirs in Saudi Arabia, and successfully explain the pore fluid, rock porosity, and pore shape trends in these complex reservoirs.

modeling, rock-physics, porosity, Middle East, log analysis

INTRODUCTION

Understanding the seismic properties of low-to-intermediate-porosity gas-saturated sandstones is essential for improved quantitative analysis and interpretation of the reservoir quality of tight reservoirs. Tight sandstones have generally undergone severe mechanical and geochemical compaction processes that cause a large reduction in their original porosity from the time of their deposition. To improve the reservoir characterization of tight sands, we use a rock-physics modeling scheme that honors the burial history of the rocks. This new methodology applies models that include the effects of the initial and pervasive cementation, and subsequent dissolution and generation of microporosity and cracks.

For unconsolidated to cemented sandstones with relatively high porosity (larger than 20%), granular medium models are often applied (Avseth et al., 2005). For this type of modeling, the amount and localization of the cement within the grain structure are the main geometrical input, along with the average number of contact points per grain (coordination number). The elastic properties of low-porosity rocks depend more on the specific pore structure than on the actual porosity, and, in general, the bulk modulus within the crack-like porosity dominates the overall fluid response (Johansen et al., 2002). Dvorkin et al. (1999) combines the contact cement model with effective medium theory to capture the transition from 0% to 100% cement concentration. Ruiz and Dvorkin (2010) extend this study by incorporating intragranular porosity, as they apply their model to carbonates.

We expand upon these studies, and present a slightly different strategy that is more in agreement with our observed diagenetic processes and associated changes in rock texture, including cement and cracks. The overall strategy of the modeling is to combine contact cement theory (CCT) as described by Dvorkin and Nur (1996), with inclusion-based modeling, in this case the differential effective medium (DEM) model as first described by Berryman (1992). The two different modeling schemes are linked using the upper Hashin-Shtrikman bound (Hashin and Shtrikman, 1963). In this way, our modeling approach captures the transition from granular dominated textures to pore dominated textures. Pore fluid effects are usually modeled as per Gassmann (1951). However, for the low porosities, we use DEM to obtain dry rock properties to be input to the Gassmann model when the pore system is considered open (low-frequency behavior). Alternatively, the DEM model may also be used to predict the fluid effects in case pores are considered isolated (high-frequency behavior). The results of the modeling are demonstrated in rock-physics templates (RPTs) (Ødegaard and Avseth, 2004; Avseth, 2005), where we compare our model trends with selected well log data from a tight gas sandstone reservoir in Saudi Arabia.

GEOLOGIC MOTIVATION

The data used in this study are comprised of a suite of well logs obtained from the Late Ordovician Sarah Formation in Saudi Arabia. The sandstones, located at approximately 3000 m depth, represent reworked glacial moraines deposited by alluvial fans and braided rivers in paleovalleys created by glaciers. Three main depositional facies have been defined for this formation (Clark-Lowes, 2005), as shown in Figure 1 (left). Facies F1 consists of medium- to coarse-grained, poorly sorted sandstones and diamictites. These are predominantly debris flows deposited at the base of the paleovalleys. Facies F2 consists of medium- to fine-grained, moderately well-sorted sandstones deposited in a braided fluvial environment. These sandstones represent the major part of the central axis of the paleovalley and overlay the F1 facies. Finally, facies F3 consists of fine-grained sandstones, siltstones, and minor shales that were mainly deposited on the margins of the paleovalley, and are sometimes found interbedded with F1 in marginal alluvial fan depositional settings. The fine-grained sedimentary rocks of F3 were deposited in quiet water pools at the valley margins, or as overbank flooding events splayed from the channel axis of the central braided river systems. The F3 sediments represent the later stages of valley infill when the braided fluvial systems evolved more into meandering river systems.

The subsequent burial history affected the three depositional facies differently (Figure 1, right). Due to the large differences in sorting and porosity, the cleaner sands of F2 would be more resistant to mechanical compaction than F1 and F3. Furthermore, the presence of a nearby shallow marine environment, and possibly repeated marine inundations of the paleovalleys, resulted in an abundance of carbonaceous ooze and shells that would dissolve and become a source of early calcite cement in the sandstones. In particular, in certain areas, F2 has extensive calcite cement (Clark-Lowes, 2005). Petrographic studies of the Sarah Formation show evidence of secondary dissolution and preservation of primary porosity in these sandstones (Figure 2). For the F2 facies, we observed sandstones where the calcite cement is still widely preserved, sandstones where the calcite has been completely dissolved and primary porosity is retained, and where the sandstones are pervasively cemented with quartz cement. Figure 2 also shows that the microstructure of the Sarah Formation sandstones is highly variable between the different facies. The combination of broad depositional variability and different diagenetic scenarios within the Sarah Formation requires a new rock-physics modeling approach that can honor the great geologic variability we observe in these rocks.

ROCK-PHYSICS MODELING

For the rock-physics modeling, we assume the porous rocks to be isotropic, which implies that the elastic stiffness tensor can be completely defined using the two elastic parameters: bulk modulus

K

and shear modulus

μ

⁠. In our data examples to follow, P-velocity (⁠

V_{p}

⁠), S-velocity (⁠

V_{s}

⁠), and density (⁠

ρ

⁠) data are available, enabling a simple transform to measured dynamic values for

K

and

μ

⁠, via

K = ρ (V_{p}^{2} - \frac{4}{3} V_{s}^{2}),

(1)

and

μ = ρ V_{s}^{2} .

(2)

The effective density is accurately determined from

ρ = (1 - ϕ) \sum_{i = 1}^{N_{s}} V_{i} ρ_{i} + \sum_{j = 1}^{N_{f}} S_{j} ρ_{j},

(3)

where

V_{i}

and

ρ_{i}

denote mineral volume fraction and density of mineral of index

i

among a total number of

N_{s}

minerals. Similarly,

S_{j}

and

ρ_{j}

are fluid saturation and density of fluid of index

j

of a total number of

N_{f}

fluids.

A further strategy for the modeling of the elastic parameters is schematically shown in Figure 3. The porosity range is divided into three intervals. The first is a low-porosity interval valid for porosities within $0 \leq ϕ \leq ϕ_{\max, DEM}$ ⁠, and where the DEM model is used to predict the elastic properties. Here, $ϕ_{\max, DEM}$ defines a preset maximum porosity for which the DEM model is used. This value has to be evaluated in correspondence with the particular pore model being used. A high-porosity $ϕ_{CCT}$ interval is then defined, where the elastic moduli are found using a cemented granular medium model, in our case the CCT model. Finally, the elastic properties in the middle porosity interval $(ϕ_{\max, DEM}, ϕ_{CCT})$ are found using the Hashin-Shtrikman upper bound (⁠ ${HS}^{+}$ ⁠), where the stiff component is defined by the elastic moduli obtained from DEM modeling of the low-porosity interval, and those of the soft component are defined at the maximum porosity point using the CCT model. Furthermore, in our approach, the stiff component becomes gradually more porous as the porosity decreases from the maximum porosity point (⁠ $ϕ_{CCT}$ ⁠) to the highest porosity point of the lowest porosity interval (⁠ $ϕ_{\max, DEM}$ ⁠). More details of the three modeling steps are given below.

Step 1.

Computation of elastic parameters in the low-porosity interval,

0 \leq ϕ \leq ϕ_{\max, DEM}

⁠, using the DEM model. The effective elastic parameters of the two-phase composite, consisting of a solid mineral with effective elastic moduli

K_{S}

and

μ_{S}

with porosity

ϕ

and fluid bulk modulus

K_{f}

⁠, are obtained solving a coupled system of ordinary differential equations (Berryman, 1992),

(1 - ϕ) \frac{d}{d ϕ} [K_{DEM} (ϕ)] = [K_{f} - K_{DEM} (ϕ)] P^{*} (ϕ)

(4)

and

(1 - ϕ) \frac{d}{d y} [μ_{DEM} (ϕ)] = - μ_{DEM} (ϕ) Q^{*} (ϕ),

(5)

where

K_{DEM} (ϕ)

and

μ_{DEM} (ϕ)

are the effective bulk and shear moduli, respectively, and starting from initial conditions

K_{DEM} (ϕ = 0) = K_{s}

and

μ_{DEM} (ϕ = 0) = μ_{s}

⁠, which in our case is set to the elastic moduli of the mineral.

P *

and

Q *

are geometrical factors associated with the aspect ratios of the pores. Norris (1985) proved that DEM results are always within Hashin-Shtrikman bounds (Hashin and Shtrikman, 1963), but their result depends on the elastic properties of the constituents, the shape of the inclusions and the integration paths (see also Mavko et al., 2009).

An alternative way of modeling fluid effects is by using the model of Gassmann (1951) where the dry rock properties

K_{d}

and

μ_{d}

are estimated considering dry pores using DEM. Hence,

K = K_{s} \frac{ϕ K_{d} - (1 - ϕ) K_{f} K_{d} / K_{s} + K_{f}}{(1 - ϕ) K_{f} + ϕ K_{s} - K_{f} K_{d} / K_{s}},

(6)

where the shear modulus is assumed unaffected by the pore fluid, i.e.,

μ = μ_{d}

⁠.

The DEM model used relies on the formulation of Kuster and Toksøz (1974), which assumes a noninteraction between inclusions. This criterion means that if all pores have the same aspect ratio, say $α$ ⁠, then the maximum porosity to be considered is restricted by $ϕ \leq α$ ⁠. This implies that a high aspect ratio model of spherical pores can only constitute a highly porous material with noninteracting pores, whereas a pore model of low-aspect ratios has a much lower maximum porosity limit. In the DEM approach, we consider that the strict noninteraction criterion given by Kuster and Toksöz (1974) can be moderately exceeded. Although the DEM model assumes isolated pores, we can obtain connected pores by incrementally adding more pores to the background medium (i.e., the background medium with isolated pores serves as an updated background medium for the next iteration where we add new pores). In practice, this implies that the choice of $ϕ_{\max, DEM}$ can be higher than the aspect ratio, but it will still be a low-porosity value that increases slightly with increasing aspect ratio. There are several alternative inclusion models to the DEM model, for instance, the self-consistent approach of Berryman (1980a, 1980b), and the T-matrix formulation of Jakobsen et al. (2003), which take into account other microscopic features of the pore space.

Step 2.

Computation of the elastic properties at the maximum porosity point

ϕ_{CCT}

using the CCT model. The shear modulus obtained from the CCT model is reduced by multiplication of a shear relaxation factor due to slip and/or heterogeneous stress at the grain contacts, as discussed by Bachrach and Avseth (2008). The dry rock elastic moduli

K_{d}

and

μ_{d}

are obtained from

K_{CCT} = \frac{n [1 - ϕ_{CCT} / (1 - S_{C})]}{6} M_{C} {\hat{S}}_{n},

(7)

and

μ_{CCT} = \frac{1}{5} K_{CCT} + \frac{3 n [1 - ϕ_{CCT} / (1 - S_{C})]}{20} μ_{C} S_{τ} .

(8)

Here, $n$ denotes the coordination number (usually defined by a correlation to porosity using Murphy, 1982), $M_{c}$ and $μ_{c}$ are the P-wave bulk modulus and the shear modulus of the cementing material, respectively, $S_{c}$ denotes the fraction of the cemented porosity of the unconsolidated granular medium, and $S_{n}$ and ${\hat{S}}_{τ}$ are normal and shear stiffness at the grain contacts, which are functions of the elastic parameters of the grains and cement, the geometrical details of how the cement is distributed (grain contact dominated versus grain coated), and the porosity before and after cementation. For more details, see Avseth et al. (2005). The cement properties are defined by the elastic moduli of the grain mineral. Due to slip and/or torsion at the grain contact interfaces, the shear modulus as predicted by CCT is usually considered to be too high (Dvorkin and Nur, 1996; Bachrach and Avseth, 2008); therefore, the tangential shear stiffness in the CCT model is normally multiplied by a shear relaxation factor of between zero and one. For loose sands, this factor will give a shear stiffness of somewhere between a Walton Rough and a Walton Smooth contact theory (associated with no slip and no friction, respectively), and it can be estimated from the dry Poisson’s ratio (Bachrach and Avseth, 2008). For cemented sandstones, it is not fully understood why the CCT overpredicts shear stiffness, and the choice of shear relaxation factor becomes a calibration procedure. Finally, the fluid saturated elastic moduli are found using Gassmann equation 6, $K_{d} = K_{CCT}$ and $μ_{d} = μ_{CCT}$ ⁠.

Step 3.

Computation of the elastic properties for the intermediate porosities,

ϕ_{\max, DEM} < ϕ < ϕ_{CCT}

⁠, using the Hashin-Shtrikman upper bound (⁠

{HS}^{+}

⁠). Here, we assume that the elastic moduli at the low-porosity point

ϕ_{\max, DEM}

are higher than those at the maximum porosity point

ϕ_{CCT}

⁠. In our modeling, the elastic properties of the infill cement at

ϕ_{CCT}

are those of the nonporous cement, which gradually alters as pores are embedded, so that the porosity of the cement and matrix material are equivalent at

ϕ_{\max, DEM}

⁠. The effective elastic moduli between these two porosities are found using

{HS}^{+}

⁠, where the stiff components are the actual modeled cement properties. Hence,

K_{HS} (ϕ) = K_{DEM}^{'} (ϕ) + \frac{V_{CCT} (ϕ)}{{[K_{CCT} - K_{DEM}^{'} (ϕ)]}^{- 1} + V_{DEM} (ϕ) {[K_{DEM}^{'} (ϕ) + \frac{4}{3} μ_{CCT}]}^{- 1}},

(9)

and

μ_{HS} (ϕ) = μ_{DEM}^{'} (ϕ) + \frac{V_{CCT} (ϕ)}{{[μ_{CCT} - μ_{DEM}^{'} (ϕ)]}^{- 1} + 2 V_{DEM} (ϕ) (K_{DEM}^{'} + 2 μ_{DEM}^{'}) {[5 μ_{DEM}^{'} (ϕ) (K_{DEM}^{'} (ϕ) + \frac{4}{3} μ_{DEM}^{'} (ϕ))]}^{- 1}},

(10)

where

V_{DEM} (ϕ) = \frac{ϕ_{CCT} - ϕ}{ϕ_{CCT} - ϕ_{\max, DEM}}

and

V_{CCT} (ϕ) = 1 - V_{DEM} (ϕ)

⁠. The porosity

ϕ^{'}

of the stiff part increases from zero to

ϕ_{\max, DEM}

⁠, as we move from the CCT end-member to the DEM end-member:

ϕ^{'} = \frac{ϕ_{\max, DEM} (ϕ_{CCT} - ϕ)}{ϕ_{CCT} - ϕ_{\max, DEM}} .

(11)

Here,

K_{DEM}^{'} (ϕ)

and

μ_{DEM}^{'} (ϕ)

are corresponding elastic properties so that

K_{DEM}^{'} (ϕ) = K_{DEM} (ϕ^{'})

and

μ_{DEM}^{'} (ϕ) = μ_{DEM} (ϕ^{'})

⁠. Again, the fluid saturated elastic moduli are found using Gassmann equation ,

K_{d} = K_{HS}

and

μ_{d} = μ_{HS}

⁠.

Figure 4 illustrates the elastic and seismic response of the rock-physics modeling. Here, we model three pore geometries (P1–P3) describing the low-porosity interval, each containing one single pore aspect ratio (⁠ $α$ ⁠): P1 has spherical pores (⁠ $α = 1.0$ ⁠), P2 has $α = 0.1$ ⁠, and P3 has $α = 0.01$ ⁠. For the modeling of the contact cement point, the volume fraction of cement is 5% with an initial critical porosity of 40%; thus, $ϕ_{CCT} = 35$ ⁠. Because the CCT model is known to overpredict shear stiffness (Dvorkin and Nur, 1996), the shear modulus predicted from the CCT model is multiplied with a shear reduction factor of between zero and one (normally set to approximately 0.5). Furthermore, the maximum porosities of the low-porosity interval for the pore models P1, P2, and P3 are set to 0.1, 0.05, and 0.02, respectively. The mineral properties of the grains, cement, and matrix in the low-porosity interval are all equivalent and defined in Table 1, which also provides the necessary fluid parameters used in the modeling. Figure 4a shows the strong effect the pore geometry has within the low-porosity interval, with decreasing aspect ratios (⁠ $α$ ⁠) resulting in a strong decrease in the elastic moduli with porosity. The fluid effect is also strongly dependent on the pore structure (Figure 4b). Although the bulk modulus of a composite with spherical pores poorly discriminates dry and water saturated pore systems, the effect of a system of pores with low aspect ratios is large. Figure 4c displays an RPT (Avseth et al., 2005) showing the V_P/V_S-ratio versus acoustic impedance of the water and dry saturated models shown in Figure 4a and 4b. Again, we see a large span in possible signatures at low to intermediate porosities. In the subsequent section, we use similar templates to review a set of well data. To incorporate the effects of heterogeneous mineralogy, a set of effective mineral elastic moduli are modeled using the approach of Hill (1965), where the volume fractions and the elastic moduli of the various constituents are used as input.

DIAGNOSTICS OF SARAH FORMATION WELL LOG DATA

In the following, we use various RPTs, derived using the approach outlined in the previous sections, to analyze a set of data from two wells drilled in different fields. The Sarah sandstone in these two wells, hereafter referred to as Well A and Well B, is characterized by low-to-intermediate-porosity and contains either water or gas pore fluid.

Figure 5 shows a comparison of the data from well A after applying different RPTs. The thickness of the reservoir zone is approximately 30 m and it contains gas and water-saturated sands. Figure 5a and 5b shows the modeled shear and bulk moduli, respectively, using the pore geometry models P1–P3 described in the previous section. The reference mineral properties have been estimated from the interpretation of the composition logs and thin sections, and are a mixture of 60% quartz and 40% calcite. The physical properties used are shown in Table 1. The elastic properties at the contact cement point were defined using a value of 5% contact cement in sand and a porosity of 40%. A shear reduction factor of 0.5 was used. Figure 5a and 5b shows that the fluid properties and porosities of the modeled and observed data are quite consistent, although the shaly data are less well captured by the model. Figure 5c–5e shows various crossplots with the data color coded with respect to gas saturation, clay content, and porosity. We see that the data plots mostly within the bounds of the models. The gas-saturated data points plots with overlapping acoustic impedance values, but lower V_P/V_S values than the water saturated data (Figure 5c). Note that the low-porosity data points (Figure 5e) correspond well with the water-saturated data points, and these plot closer to the model with low aspect ratio.

Figure 6 displays the same type of crossplots as in Figure 5, but using the data from well B, where the reservoir thickness is approximately 100 m. The data and model trends from the RPTs again seem to match fairly consistently. The models nicely capture the variability in the data, and the varying pore aspect ratios at the low porosities explain why the gas saturated and the brine saturated sandstones can have overlapping V_P/V_S ratios.

DISCUSSION

A workflow for rock-physics modeling of well-consolidated sandstones varying from intermediate-to-low-porosities has been proposed. The main purpose of the workflow is to enable prediction of the elastic properties in the transition zone from one type of pore system, associated with higher porosity rocks of granular composition, to that of lower porosity rocks characterized by a suite of specific pore shapes. It is well known that pore shape plays a dominant role in the overall elastic behavior of rocks, and their response to the pore fluid.

The suggested modeling approach reflects the true geologic transition from poorly to well-consolidated sandstones associated with deep burial and high temperatures. Furthermore, the use of the DEM inclusion model with varying aspect ratios can account for microporosity resulting from dissolution of calcite cement. The Sarah Reservoir sandstones investigated in this study may also have been exposed to episodic uplift (not depicted in the schematic burial trends in Figure 1b), and these tectonic events may have resulted in faults and associated microcracks. Thus, we expect large variability in elastic stiffnesses at lower porosities, and this is accounted for by varying the aspect ratios used in the DEM modeling. This is apparent in the large variability in V_P/V_S ratios, which makes this particularly valuable for pore fluid detectability in these rocks. The presence of gas in these reservoirs may also be detectable using seismic amplitude versus offset (AVO) data analysis. We have not investigated any anisotropy effects and these could be significant in these rocks. So, before any AVO type analysis is performed, the effects of anisotropy on the elastic properties of these reservoirs should be modeled.

There is also some uncertainty as to whether the proposed models will be representative for the complete porosity range. The validations and data observations are only from a selected depth interval with a limited range of porosities (0%–25%). The fact that the CCT high-porosity end member (35%–35%) is modeled with an assumed initial contact cement volume of 5%, and an assumed shear relaxation factor of 0.5, without enough data coverage for good calibration means that this high-porosity end member is somewhat uncertain. This point is probably more realistic than a loose sand end member, because the studied sandstones have been exposed to early calcite cement. The contact cement model we used has been found to be quite good for high-porosity sandstones (e.g., Dvorkin and Nur, 1996; Avseth et al., 2005, 2010). Nevertheless, we suggest testing and validating the modeling approach presented in this paper on more data sets spanning a larger range of porosities and over a larger compaction window.

The mineral point also has some inherent uncertainty as we expect significant variability in mineral composition in the detrital grains as well as in the authigenic matrix. Nevertheless, we have found that most of these Sarah Reservoir sandstones are quite quartz rich, with abundant amounts of calcite cement. Some wells in the area, not included in this paper, showed lower calcite volume, and the data in these wells were better described by RPTs where 100% quartz was selected as the mineral point.

The modeling approach presented in this study may be quite useful for tight carbonate reservoirs and shale gas plays. The templates generated can be used to guide the seismic interpretation and to constrain seismic inversions in low-to-intermediate-porosity rocks. Note that the well log data in Well B (i.e., Figure 6) do not show clear trends, but are bounded by our models. The competing effects of porosity, pore aspect ratio, saturation, and clay volume will obscure any clear trends, but with the use of our models we can understand what is controlling the scatter in different directions. In other words, we demonstrate that the scatter in the data is not random or due to noise, but is a result of complex geology.

CONCLUSIONS

We have established a rock-physics modeling strategy that combines high-porosity granular models with low-porosity inclusion models. The modeling approach takes into account the burial history of a tight gas sandstone interval encountered in Saudi Arabia, and successfully predicts the rock-physics and seismic properties as a function of rock texture and pore fluid characteristics of these sandstones. These models can be used to better understand the seismic signatures of the studied Sarah sandstone reservoirs or other low-to-intermediate-porosity reservoirs elsewhere.

ACKNOWLEDGMENTS

We would like to thank Saudi-Aramco for financing this study and for permission to publish the results of this research.

Mineral or fluid	Bulk modulus (GPa)	Shear modulus (GPa)	Density (⁠ $g / ccm$ ⁠)
Quartz	36.0	44.0	2.65
Calcite	77	32.0	2.72
Water	3.58		1.17
Gas	0.10		0.28

Mineral or fluid	Bulk modulus (GPa)	Shear modulus (GPa)	Density (⁠ $g / ccm$ ⁠)
Quartz	36.0	44.0	2.65
Calcite	77	32.0	2.72
Water	3.58		1.17
Gas	0.10		0.28

Rock-physics modeling guided by depositional and burial history in low-to-intermediate-porosity sandstones

ABSTRACT

INTRODUCTION

GEOLOGIC MOTIVATION