## Abstract

Extracting information about reservoir quality from seismic data is a key challenge in exploration, appraisal and production of hydrocarbons. We demonstrate how to perform quantitative reservoir characterization by using inverse rock physics modelling on seismic inversion data. This allows us to evaluate the non-uniqueness of our predictions. We demonstrate our methodology on a gas–condensate Norwegian Sea field under appraisal and production, and perform reservoir quality predictions along a selected seismic cross-section where we have well control. Even though such a seismic dataset is more uncertain than well log data, which have been used previously in similar analysis, we still achieve reasonable and consistent predictions of reservoir quality.

## Introduction

Robust and flexible methods for estimation of reservoir properties from seismic data are essential for quantitative interpretation in reservoir characterization and monitoring. Several methods are commonly used to assess reservoir conditions from seismic reflection amplitudes. The conventional way is to perform qualitative seismic interpretation to outline geological structures and reservoir architecture from seismic reflection events, their geometry and character. Going further, amplitude v. offset (AVO) screening (e.g. Rutherford & Williams 1989) is often used to highlight lithology and fluid information from angle-dependent seismic reflection amplitudes in the seismic pre-stack domain. However, for better quantitative interpretations of seismic data, more sophisticated methods are required.

Rock physics addresses a direct link between parameters describing reservoir properties (e.g. porosity, lithology and fluid saturation) and its effective (or upscaled) rock properties (e.g. compressional and shear impedances and velocities) that may be inferred from seismic observations at a scale characterized by the seismic wavelength. Hence, a stronger integration of rock physics in reservoir characterization is a reasonable approach to achieve a more quantitative interpretation. Forward rock physics modelling of effective rock properties from reservoir properties can be challenging if input model parameters are poorly constrained, but is commonly feasible. However, reservoir characterization based on seismic data is an inverse problem that includes several obstacles:

Similar rocks can exhibit very different effective rock properties: for example, unconsolidated and consolidated rocks with equivalent composition give different seismic responses. Hence, a library of models is required to describe the various types of rocks that exist in a reservoir setting.

Rock properties continuously change throughout burial history, which may yield different models for the same rock unit.

Various combinations of reservoir properties may have equivalent seismic properties: for example, high-porosity clean sands may have the same seismic response as low-porosity shaly sands.

The scale of geophysical measurements must also be addressed according to the scale of rock heterogeneities relative to the applied signal wavelength: for example, ultrasonic velocity measurements of a small rock sample may differ from seismically derived velocities that probe macroscopic volumes of the same rock for reasons that are the result of physics (frequency-dependent properties) and statistics (different sample volumes).

Consequently, the inverse problem is complicated to solve and has non-unique solutions (Bosch *et al.* 2010; Dvorkin *et al.* 2014). Yet, both deterministic and probabilistic methods have been formulated to achieve a seismic-to-reservoir transform. For example, one may use statistical methods in forward rock physics modelling to assess the probabilities of various reservoir scenarios (Avseth *et al.* 2005; Doyen 2007). Furthermore, a Bayesian seismic inversion workflow for lithology–fluid classification can be used when *a priori* information about anticipated lithology–fluid scenarios is available (e.g. Buland *et al.* 2008; Grana & Rossa 2010; Rimstad *et al.* 2012). A more simple and, possibly, more accessible method was introduced by Ødegaard & Avseth (2004), and is referred to as the rock physics template (RPT) approach. Here, a so-called template for a rock physics model is projected on top of a cross-plot derived from data, allowing us to make an interpretation of reservoir properties. Among others, Chi & Han (2009) and Avseth *et al.* (2014) demonstrated applications of RPT on seismic data. A further review and discussion of various methods can be found, for instance, in Helgesen *et al.* (2000), Doyen (2007) and Bosch *et al.* (2010).

In our study, we perform a quantitative seismic interpretation based on seismic inversion data, which serve as input into the inverse rock physics modelling (IRPM) strategy described by Johansen *et al.* (2004, 2013). The approach is useful for evaluating various rock physics models (Moyano *et al.* 2011) to inspect and constrain the non-uniqueness of solutions. Bredesen *et al.* (2013) have evaluated the method using synthetic data. In this study, using data from the Norwegian Sea, we demonstrate its application using real seismic data as input into the IRPM for predicting reservoir properties away from well locations.

## Forward and inverse rock physics modelling

We assume isotropic rocks in order that their seismic properties can be derived from the bulk modulus (*K*), shear modulus (*µ*) and density (*ρ*). The *K* and *µ* properties can be predicted from various rock physics models based on a set of reservoir parameters (see Appendix A for details). Three parameters are particularly important for reservoir characterization: porosity (*f*), the clay-to-sand fraction (*C*) (referred to as lithology) and the pore fluid saturation (*S*) (pore volume fraction of the pore fluids). These parameters and fluid permeability describe, in essence, the reservoir quality. Although predictions of fluid permeability are not considered in this study, in clastic formations, permeability typically increases with decreasing clay content and increasing porosity.

In the present study, we use inverse rock physics modelling (IRPM) (Johansen *et al.* 2004, 2013) implemented in a three-step procedure, comprising: (1) seismic inversion; (2) forward rock physics modelling; and (3) inverse rock physics modelling (Fig. 1):

**Step 1 – Seismic inversion:**in this study, simultaneous AVO inversion has been applied to derive acoustic impedance (AI), the P-wave to S-wave velocity ratio (*V*_{P}/*V*_{S}), and density (*ρ*) as input into the IRPM. This is a pre-step to the IRPM, and had already been performed by a third party in another project. As such, we were provided a set of seismic inversion data for our study, and thus our focus is on Step 2 and Step 3.**Step 2 – Forward rock physics modelling:**to define rock physics constraints relating reservoir parameters to seismic parameters, such as seismic velocities, acoustic impedances, elastic moduli and density.**Step 3 – Inverse rock physics modelling:**to find possible solution sets of reservoir parameters for a given set of seismic parameters consistent with suitable rock physics models.

An example of the results from Step 2 is shown in Figure 2 for a partly consolidated rock, where the *x*-, *y*- and *z*-axes denote porosity, lithology and gas saturation, respectively. One value of the bulk modulus corresponds to numerous possible reservoir parameter combinations, forming a surface in the reservoir parameter domain. We denote this surface as an isosurface because all points correspond to the same bulk modulus value. Hence, any combination of reservoir parameters corresponding to a point on such an isosurface is, in principle, a possible solution for the given value of the bulk modulus, which may derive from seismic inversion or well data. To constrain the number of solutions further, we correlate the bulk modulus isosurface with similar isosurfaces from other input data (e.g. shear modulus and density). Solutions are defined where the various isosurfaces intersect.

Figure 3a, b shows examples of solution sets of reservoir parameters estimated by combining density, and bulk and shear modulus. The solution sets are, of course, more constrained using three rather than two data parameters: however, we still observe non-uniqueness in both cases. High-quality, modern seismic data can enable inversion of two or three seismic parameters (e.g. AI, *V*_{P}/*V*_{S} and *ρ*), which, in turn, place tighter constraints on reservoir parameter solution sets (Avseth *et al.* 2005), as demonstrated in Figure 3. Some combinations of data parameters are better than others in constraining solutions and when considering solution stability (Jensen 2011). Nevertheless, our ability to investigate the non-uniqueness of estimated solutions is of vital importance in reservoir characterization. Thus, as we will demonstrate later, IRPM is a useful tool for testing the suitability of various rock physics models for a given dataset and to constrain corresponding solutions.

It is important to include uncertainties in reservoir parameter predictions. In our case, we have: (1) well and seismic data errors and limited signal frequency bandwidth; and (2) model errors due to shortcomings in the theoretical models and uncertainties in the input model parameters. One implementation of the IRPM solver identifies solutions only where all isosurfaces intersect. However, when dealing with the mentioned data uncertainties, the isosurfaces can be shifted, which implies a shift in the intersections and, hence, the predictions. In some cases, the isosurfaces might be shifted in such a way that previously intersecting isosurfaces now no longer intersect. To manage these implications, we use a so-called proximity-based solver that includes a tolerance factor (Moyano *et al.* 2011, 2015). Now, possible solution sets are identified as long as the isosurfaces are within a maximum distance from each other in the reservoir parameter domain. Hence, possible solution sets are defined not just to be at the intersection of the isosurfaces, but also within an expanded volume around them. Furthermore, it is not required that isosurfaces actually intersect for the solver to define a solution set.

## Application to Norwegian Sea dataset

### Regional setting

The inverse rock physics modelling (IRPM) methodology has been demonstrated on data from a Norwegian Sea gas–condensate field of Jurassic age (Fig. 4). The Jurassic reservoirs comprise two gas-bearing sandstone units named the Garn and Ile formations (Fig. 5a). These reservoir units are separated by a thin layer of shale, the Not Formation. The geological characteristics of these sandstone reservoirs are very similar, but the Ile Formation contains more thin layers of shale and is somewhat more consolidated. Both the Garn and the Ile formations were deposited in a shallow-marine environment. The reservoir quality is considered excellent, with porosities ranging from around 25 to 35%. Tectonic rifting and Permian salt remobilization, occurring during late Jurassic times, created fault blocks separated by narrow graben synclines. Late Jurassic erosion of the fault blocks may have caused the Garn Formation sands and Ile Formation sands to be redeposited into the graben lows, and they are then referred to as the Rogn Formation sands. The late Jurassic sequence also includes organic-rich shales of the Melke and Spekk formations, in which the Rogn Formation sands may be embedded.

### Step 1 – Seismic inversion

The seismic reflection data used in this study are three-dimensional (3D) broadband, pre-stack time-migrated (PSTM) data from a multi-client survey, and cover our study area as shown in Figure 5b. The data have been pre-conditioned for AVO inversion and the processing includes the following key steps: designature to zero phase; noise attenuation/filtering; partial Q-compensation (phase only); Kirchhoff pre-stack time migration; stacking velocity analysis; normal moveout correction; Radon multiple removal; estimation of angle gathers; wide-angle mute; and time alignment. No angle-dependent scaling was carried out, as this will be accounted for by angle-dependent wavelets in the simultaneous AVO inversion (see below). Figure 6 shows the data zoomed in at our focus area, together with inverted acoustic impedance (AI), the P-wave to S-wave velocity ratio (*V*_{P}/*V*_{S}) and density. The green and black dashed lines, and dots, respectively, denote interpreted horizons of target formations and the log data from well A in Figure 7 projected on top of the sections. The seismic section extends 1.3 km horizontally and 130 m vertically, covering the Garn and Ile formations (31 and 61 m in thickness, respectively) at the well location.

The seismic inversion data were obtained from a simultaneous seismic AVO inversion from 3D angle-stack volumes generated from a set of offset gathers. In a simultaneous AVO inversion (Ma 2002; Rasmussen *et al.* 2004; Russell *et al.* 2006), the reflectivities are converted to absolute elastic parameters, using a low-frequency, *a priori* starting model and iteratively updating this Earth model until forward modelled seismic data fit the real seismic data sufficiently well. The low-frequency model is made from smoothed well log data and seismic interval velocities, the latter combined with empirical rock physics relationships (the Greenberg–Castagna equation for the estimation of *V*_{S} and the Gardner equation for density estimation: see Mavko *et al.* 2009). The Aki & Richards (1980) linearized approximation to the Zoeppritz AVO equations is used to model angle-dependent reflectivities:

where , and , and is equivalent for the S-wave velocity and density. The angle-dependent reflectivities are then convolved with angle-dependent wavelets. These were estimated using spectral analysis of reflectivities from wells A and B, together with the seismic data along the wellbore.

The AI is derived from near-offset reflectivities, and tends to be relatively reliable and well constrained to other inverted parameters when well log data are available for calibration, and when seismic data quality is adequate. However, the Aki & Richards (1980) approximation assumes small contrasts in physical properties across boundaries, small incidence angles, plane-wave propagation, plane boundaries and isotropic rocks; assumptions of which may incur errors in the inverted parameters. The *V*_{P}/*V*_{S} ratio is estimated from angle-dependent reflectivities up to 30° when using a two-term AVO equation, but can be even more constrained via a three-term AVO inversion using angles of up to 45°. Nevertheless, *V*_{P}/*V*_{S} estimates are more sensitive than AI to errors (e.g. Escobar & Hansen 2010). Densities can also be obtained from a three-term AVO inversion, but this requires data with large incident angles unless reflections are dominated by density contrasts. Density estimation may be affected by anisotropy and refractions interfering with reflection amplitudes at larger offsets. Density suffers most from the ill-posed nature of the inversion problem, making its estimates more noise-sensitive (Downton & Lines 2004; Behura *et al.* 2007). Moreover, one of the key uncertainties related to estimation of both *V*_{P}/*V*_{S} and density is associated with the offset-to-angle transform, which requires a good velocity model for the ray-tracing algorithm to be reliable (Chopra & Castagna 2014). Finally, the low-frequency starting models for the different layer parameters and the angle-stack wavelets are also significant sources of uncertainty in the absolute values of inverted parameters.

Figure 8a shows the quantitative match between upscaled elastic logs (black lines) obtained from a simple average (Tiwary *et al.* 2009) and seismic inversion data (red lines) via estimated normalized root-mean-square error (NRMSE) (see Appendix B). As expected, the AI and *V*_{P}/*V*_{S} ratio generally show smaller errors than the density, and their NRMS errors are fairly similar. We also see a general higher error within the shallower Garn Formation than in the deeper Ile Formation. Figure 8b shows a synthetic trace generated from a zero-phase Ricker wavelet convolved with reflectivities obtained from the well log superimposing the seismic reflection amplitudes. The red peaks imply decreasing acoustic impedance that is consistent with top Garn and Ile reservoirs. A good match between the seismic inversion and the well log data is also demonstrated in Avseth *et al.* (2014), who studied the same dataset.

### Step 2 – Forward rock physics modelling

One major challenge in seismic reservoir characterization is to find and calibrate rock physics models that provide appropriate relationships between the well data (i.e. between the rock composition) and the seismic properties. The geological and geophysical characteristics of rocks change continuously from time of deposition through burial via mechanical and chemical processes. To describe the possible rock physics behaviour of this transition, we need a library of rock physics models to apply (Avseth *et al.* 2010).

We use well A containing P- and S-wave velocities and density logs (Fig. 7) to obtain the needed rock physics calibration of the reservoir formations. Rock physics templates (Avseth *et al.* 2005) were used in an initial screening phase to identify relevant model candidates (Fig. 9). These were, in turn, further evaluated and calibrated by using the IRPM. Figure 10a, b shows an example where the IRPM is used to compare the effects of homogeneous v. patchy fluid saturation. The black and coloured dots denote, respectively, petrophysical well observations and IRPM solutions, where the non-uniqueness of solutions is reflected by the width of solution ranges. The dry rock elastic moduli are obtained using a patchy cement model (Avseth & Skjei 2011; Avseth *et al.* 2012) combining contact theory (Walton 1987) and contact cement theory (Dvorkin & Nur 1996). Details of the final rock physics model are given in Appendix A.

Figure 11 shows the NRMSE between observations and the IRPM solutions in Figure 10a, b for a quantitative evaluation of the different fluid models. Error in predicted porosity and lithology for the two fluid mixing models are almost the same, whereas the saturation solutions have an average of 7.4% less error for patchy relative to homogenous saturation. Hence, a patchy saturation was selected to represent the Garn and Ile formations. However, the performance of our rock physics model varies with the different reservoir facies (Fig. 10b). For instance, better-constrained porosity solutions are obtained in the Garn Formation reservoir compared with the Ile Formation reservoir: the opposite trend is seen, however, for the saturation. The reservoir characterization might be improved by fine tuning the model for different reservoir facies: for example, it might be that Garn Formation reservoir is better described using a less pure patchy fluid mixing model. However, using a single model is robust enough for applications to low-frequency seismic data, which reflect an effective medium at the macroscopic scale.

Even though we have found a rock physics model that gives fairly low errors between predictions and observations at the well location, there is no guarantee that it will perform equally well away from the well location. This is particularly important when performing reservoir characterization in a prospect-hunting stage between wells in a frontier basin. Hence, we tested the robustness of our IRPM modelling by applying it to well B, which penetrates the same reservoir formations some distance away (Fig. 12). Although the Garn and Ile formations are buried slightly deeper at well B, the NRMSE (green lines) is close to zero for all reservoir parameters in the Garn Formation, whereas it is somewhat higher and varies more in the Ile Formation. These test results are encouraging, as we seek to use the rock physics model for reservoir characterization away from the well location.

### Step 3 – Inverse rock physics modelling

Figure 13 shows the estimated reservoir properties using the IRPM with AI, *V*_{P}/*V*_{S} and density as inputs. We have chosen to plot the mean values together with the corresponding standard deviations to examine the solutions and associated non-uniqueness, respectively. Again, the dashed black lines represent the interpreted horizons as in Figure 6 and the white patches represent zones where no solutions can be found. We see zones with a high porosity consistent with a low clay content in which we also expect to predict high gas saturation. However, the IRPM using AI, *V*_{P}/*V*_{S} and density for the Garn Formation yields inconsistent results.

As previously mentioned, density is often a problematic and uncertain parameter to invert from seismic data. Hence, we repeated the IRPM, but now using only AI and *V*_{P}/*V*_{S} as input data (see Fig. 14). The main difference between the two solution sets is that the high gas saturation is now more consistent with high porosities and low clay content when using only AI and *V*_{P}/*V*_{S}. In addition, we also see fewer and smaller white gaps representing missing solutions. This is because two isosurfaces are more likely to have a common intersect than three (see Fig. 3). Comparing the standard deviations for the two modelling scenarios confirms this: they have increased significantly compared to modelling with three parameters. However, the standard deviations within the Ile Formation are low relative to the surrounding rocks, as well as within the Garn Formation for the saturation.

To evaluate the solution accuracy, we matched the IRPM solutions along the well path with the petrophysical logs in well A. We performed a sonic-to-seismic upscaling using a simple average of the petrophysical well logs over a 7.5 m window to approach the seismic scale and then superimposed them on the solutions, as shown in Figure 15. We see, however, that the upscaling is insufficient to obtain a close match, as the averaged well logs still represent a finer-scale measurement than the IRPM solutions. The green lines are estimated NRMSE between the solutions and the upscaled logs, showing slightly higher errors on average than was seen when performing the IRPM directly on the well logs. The IRPM solutions follow a similar trend to that observed on the well logs, but generally underestimate the porosity, as reflected by the porosity error. The saturation error is smaller in the Ile reservoir than in that of the Garn, whereas the lithology error is fairly low within both reservoir units. Note that solutions outside the reservoir zones are expected to be less reliable, as the rock physics model is specifically designed to represent the reservoir rocks.

## Discussion

Quantitative interpretation of reservoir quality from seismic data is more challenging and uncertain than interpretation from well log data. This is also the case when using inverse rock physics modelling (IRPM) to perform reservoir characterization. Based on the results from our Norwegian Sea demonstration, in this section we will discuss some of the associated uncertainties and possible improvements to the IRPM for reservoir quality predictions away from well locations.

The misfit of some data points between our predicted porosity, lithology and gas saturation and the upscaled petrophysical well logs (Fig. 15) can have several causes. Issues with or uncertainties in the seismic inversion data are a possibility. The current seismic inversion data show a fairly good match to the elastic well logs (Fig. 8). However, the misfit is seemingly correlated with the errors in the seismic inversion data: for example, the opposite trends in predicted and observed porosities in the top of the Ile Formation (Fig. 15) correspond closely with high data errors (Fig. 8). Another possible explanation is that seismic and well log measurements are taken at two different scales, such that the corresponding measurements can deviate. The upscaling of the petrophysical well logs is intended to compensate for scale-dependency, but is not a rigorous approach and can introduce new uncertainties (Tiwary *et al.* 2009), especially if the different volumes of investigation of seismic and well log measurements encompass lateral heterogeneity. Furthermore, for vertical heterogeneity introduced by interbedding shales, Avseth *et al.* (2010) suggested performing a Backus average with varying net-to-gross (NTG) ratios to study how elastic moduli change from well log scale to seismic scale. However, the Garn and Ile reservoirs are not characterized by considerable heterogeneity and we have therefore not performed a Backus average. Our rock physics model remains calibrated to the original well log scale where model parameters are better fine-tuned (Fig. 10).

The white gaps within the Ile Formation reservoir unit in Figures 13 & 14 are zones where the IRPM does not find solutions. These are where the P-wave to S-wave velocity ratio (*V*_{P}/*V*_{S}) is very low, consistent with high gas saturation. The rock physics model shows strong sensitivity in reservoir quality to such low *V*_{P}/*V*_{S} values. Such sensitivity is illustrated in Figure 16, where there are no intersections of isosurfaces (dashed borders) and there is a neighbouring zone where solutions are found (solid borders). This shows that, even when the fractional changes in the three reservoir parameters are comparable, the acoustic impedance (AI) and density (*ρ*) isosurfaces vary only slightly relative to the large change in the *V*_{P}/*V*_{S} isosurface. Figure 16 also illustrates how uncertainties in general may propagate, and amplify, through to uncertainty in reservoir-quality predictions because of isosurface sensitivity. An uncritical increase in the tolerance factor in order to capture solutions where the isosurfaces do not intersect is not a good strategy for several reasons: it ignores the cause of the missing intersections, which could be due to noisy data or inappropriate rock physics modelling; it treats all data as being equally uncertain; and it ignores the differing sensitivities of the isosurfaces to uncertainty. In our particular case, it is reasonable to assume that the zones with missing solutions inside the Ile Formation reservoir are zones of high gas saturation.

Reservoir characterization guided by model-based seismic inversion data will become less certain as the distance away from wells increases. As a remedy, one can increase the level of uncertainty in the modelling, which for the IRPM means increasing the tolerance factor. Moyano *et al.* (2015) showed how increasing the tolerance factor typically results in decreased prediction accuracy: decreasing it, however, would have the reverse effect. Wider solution ranges do not guarantee more robust predictions and can, in fact, lead to misinterpretations when considering mean values. Because the tolerance factor is a simplified approach for including uncertainties with no direct relationship to data or model uncertainties, uncritical increase of it is not recommended. Instead, it should be specified as part of the model calibration process when using IRPM.

The implementation of the tolerance factor to handle uncertainties is certainly not ideal. An alternative approach would be to apply probability density functions (PDF) to the input data and model parameters, which, of course, can be based on their derived or assumed uncertainties. A Monte Carlo simulation can then be used to identify possible solution sets and associated probabilities (Houck 2002). If *a priori* information about anticipated lithology–fluid scenarios is also available, it can serve as input to a full Bayesian implementation of the IRPM workflow. In practice, this will provide a certain thickness to the isosurfaces determined by the PDFs, resulting in a cluster of solution sets around the intersections of various isosurfaces. This would definitely be an improvement to the IRPM as it gives direct control over the different uncertainties, and probabilities can be estimated for the various predictions.

## Conclusions

We have demonstrated inverse rock physics modelling for quantitative reservoir characterization based on real seismic data, covering a Norwegian Sea gas–condensate field under appraisal and production. The seismic data used had already been inverted to acoustic impedance (AI), P-wave to S-wave velocity ratio (*V*_{P}/*V*_{S}) and density.

We used petrophysical well observations intersecting the seismic section of study to evaluate a range of rock physics models for the gas sandstone reservoirs by using inverse rock physics modelling. Finally, we selected a patchy cementation model that gave well-constrained solutions. We then used this model in the inverse rock physics modelling on a seismic cross-section to predict porosity, clay-to-quartz volume fraction and gas-to-brine fluid saturation away from the well location. We compared solutions obtained using AI and *V*_{P}/*V*_{S} to solutions based additionally on density; neglecting the density parameter gave more consistent predictions of the reservoir properties, but also less constrained solutions with higher standard deviations. The method opens up for better handling and quantification of uncertainties, and for scrutinizing probabilities denoting various reservoir scenarios. However, our study demonstrates reasonable predictions consistent with rock physics theory, and with regional geological, geophysical and petrophysical information.

## Acknowledgements and Funding

We thank Tullow Oil for sharing data and useful input. We would also like to thank Tullow Oil, Christian Michelsen Research and Statoil for financing the PhD programme of Kenneth Bredesen.

## Appendix A: Rock physics modelling

We have used rock physics models to estimate the effective bulk modulus (*K*) and the effective shear modulus (*µ*) related to the seismically inverted acoustic impedance (AI) and the ratio of *P*-wave velocity (*V*_{P}) and *S*-wave velocity (*V*_{S}) (i.e. *V*_{P}/*V*_{S}), via:

The effective density (*ρ*) is expressed as:

where *f, C* and *S* denote porosity, lithology and fluid saturation, respectively. The indexes b, o and g denote brine, oil and gas, respectively, and *S*_{b} + *S*_{o} + *S*_{g} = 1.

### Hashin–Shtrikman–Walpole bounds

The Hashin–Shtrikman–Walpole (HSW) bounds (Walpole 1996*a*, *b*) represent theoretical upper and lower limits of an isotropic mixture. For a mixture of two constituents, the elastic moduli are given by:

where *V, K* and *µ* are the volume fraction, bulk modulus and shear modulus, respectively, and subscripts 1 and 2 denote the two constituents. When *K*_{m} and *µ*_{m} refers to the stiffest constituent, upper bounds are found; when *K*_{m} and *µ*_{m} refer to the softest constituent, however, lower bounds are found.

### Walton contact theory

The contact theory (CT) of Walton (1987) assumes an idealized package of spherical grains often used for unconsolidated sediments. The effective bulk modulus (*K*_{W}) is then given by:

where *f, C* and *P* are porosity, coordination number and pressure, respectively, and *B* is given by:

where *µ* and *λ* are the shear modulus and Lamé’s coefficient of the grain constituent, respectively. If we assume smooth grain surfaces, a lower effective shear modulus is given by:

If, however, rough grain surfaces are assumed, an upper effective shear modulus is given by:

where *v* is Poisson’s ratio of the grain constituent. For intermediate roughness, we can use a weighted average between *μ*_{W–} and *μ*_{W+}, given by:

where *f* is a friction factor between 0 and 1, representing smooth and rough grains, respectively.

### Contact cement theory

The contact cement theory (CCT) of Dvorkin & Nur (1996) provides the elastic moduli of cemented sediments (*M*_{c}) via:

where *f*_{0} is the critical porosity. *ρ*_{c}, *V*_{Pc} and *V*_{Sc} and *M*c are the density, and P- and S-wave velocities and compressional modulus of the cement constituent, respectively. The *S*_{n} and *S _{τ}* parameters are related to the normal and shear stiffness, respectively, and vary with the cement and grain properties, together with the amount of contact cement.

### Patchy cement model

The patchy cement model (Avseth & Skjei 2011; Avseth *et al.* 2012) is valid for granular rocks characterized by a mixture of uncemented and cemented grain contacts, as described by CT and CCT, respectively. One approach is to use the HSW bounds in equations A5 and A6 to mix CT and CCT for high-porosity patchy cemented sediments. Then, a lower HSW is used to model porosity-sorting effects down to the mineral point at zero porosity. Figure A1 shows an example of the effective bulk modulus as function of porosity when using a lower HSW to mix uncemented (blue line) and cemented (red line) high-porosity sands at various ratios (dashed lines). Here, the CCT model corresponds to a 7% porosity loss due to cementation (dashed brown line), whereas the CT is calculated at a non-cemented 40% critical porosity (small brown dot).

### The Gassmann fluid substitution recipe

The rock physics models presented so far predict the dry rock elastic moduli. If, however, rocks are saturated with a pore fluid of bulk modulus *K*_{f}, the shear modulus (*μ*_{sat}) equals the dry shear modulus, and the bulk modulus (*K*_{sat}) can be calculated from the Gassmann (1951) equation:

where superscripts d and s denote dry and solid bulk modulus, respectively, and *f* is porosity.

In case of a homogenous saturation of *N* fluids, *K*_{f} can be estimated using the equation from Wood (1955):

where *S _{i}* and

*K*denote the saturation and bulk modulus of the

_{i}*i*-th fluid. Likewise, a patchy saturation can be approximated using the equation by Voigt (1928):

### Rock physics modelling for the Norwegian Sea demonstration

Several rock physics models were tested in the demonstration of inverse rock physics modelling applied to a Norwegian Sea dataset. We finally selected a patchy cement model for the dry rock elastic moduli where a 70%:30% mixing between, respectively, CCT and CT was calibrated using the lower HSW bound. This implies that 70% of the grain contacts in the gas sandstone reservoirs are cemented, whereas the remaining 30% are still uncemented. The mineral and fluid constituent properties are shown in Table 1, whereas Table 2 summarizes the remaining modelling details.

## Appendix B: Normalized root-mean-square error (NRMSE)

The normalized root mean square error (NRMSE) is used to measure the difference between a solution (*X*_{sol}) and a reference value (*X*_{ref}) via:

where *X* can, for instance, be a reservoir or seismic parameter for a sample, *z. X*_{ref,max} and *X*_{ref,min} refers to the maximum and minimum reference value of all samples within a depth interval, respectively. *X*_{ref,z} and *X*_{sol,z} in Figure 8a represents the upscaled elastic logs (black trends) and seismic inversion data (red trends), respectively. In the remaining NRMSE calculations, *X*_{ref,z} refers to petrophysical well logs (i.e. the black dots in Figs 10, 12 & 15), whereas *X*_{sol,z} are mean values of porosity, lithology and fluid saturation solutions for each depth *z* (i.e. the coloured dots in Figs 10, 12 & 15).

- © 2015 The Author(s)