Detection of Permeable Bodies: From Laboratory Measurements to Seismic Measurements

1. Introduction

The seismic reflection method has the advantage of providing a picture of the subsurface in three dimensions (3D) with a regular grid. In high resolution seismic surveys, the size of the grid cell is about tens of meters for horizontal distances, and of several meters for vertical distances. The classical approach of seismic processing leads to obtain migrated seismic sections in 2D or migrated seismic blocks in 3D. It is indispensable to have a good velocity model to carry out the migration process. The migrated sections can then be transformed into acoustic impedance sections if well data, mainly acoustic and density logs, are available. The procedures used to obtain acoustic impedance sections are often referred to as Model-based seismic inversions which require an a priori impedance model (obtained from well data) which is iteratively refined so as to give a synthetic seismic section to match the seismic section to be inverted. The impedance model is the fruit of a processing sequence which is usually done in time. 3D seismic impedance blocks and logging data are often used for geological model building in time. The geological model must be then converted from time to depth, thanks to the computation of a time to depth conversion model.

The final impedance model can be converted into porosity by using an empirical relationship between porosity and acoustic impedance established at well locations. To model porosities, an other option is to use porosity at wells location and interpolate between the wells by means of kriging. Partly due to the small number of wells, this outcome is really smooth and usually does not seem geologically consistent. More dense information can be integrated in order to improve the estimation of porosity. As porosity is linked to acoustic impedance, it is relevant to use dense seismic acoustic impedance information. So a collocated cokriging of porosity integrating the seismic information is performed by using the normalized acoustic impedance as the secondary variable [1]. 3D cube makes possible to provide 3D imaging of the connectivity of the porous bodies [2]. Core analysis is usually carried out to establish porosity vs. permeability laws [3].

It has been shown that it is possible to extract new attributes from seismic sections, leading to a better understanding of the distribution of the porous and permeable bodies [4]. The attributes are also used to detect the impermeable layers. The methodology is based on laboratory experiments which have shown that a formation permeability indicator can be obtained via the computation of four input data: P-wave frequency, attenuation, porosity and specific surface. The procedure has been firstly conducted in acoustic logging to estimate permeability of porous layers and to detect water inflows [5]. In seismics, the processing is performed in order to measure these parameters. The analytic signal is used to compute the instantaneous frequency and attenuation (Q factor). The porosity and specific surface are computed from seismic impedances obtained by acoustic inversion of the migrated seismic sections. The input parameters are used to compute a new index named Ik-Seis factor. (Indicator (I) of permeability (k) from acoustic or seismic (Seis) data).The chapter is written mainly from two journal articles: “Detection of porous and permeable formations: from laboratory measurements to seismic measurements” [4] and “Characterization of geological formations by physical parameters obtained through full waveform acoustic logging” [5].

In a first step we show how the methodology has been transposed from laboratory measurements to geophysical data. Then we illustrate the potential of the proposed procedure via field examples: acoustic logging, 2D seismic line with acoustic impedance inversion after migration, 3D seismic section with elastic impedance inversion after migration and depth conversion.

2. From laboratory experiments to geophysical data

We present a short review of the laboratory measurements conducted by P. Morlier and J.P. Sarda [6]. The possibilities of using laboratory results for field geophysical applications are then discussed.

2.1. Laboratory measurements

Laboratory experiments [6] have shown that the attenuation of a clean formation can be expressed in terms of three structural parameters: porosity, permeability and specific surface. Both theoretical and experimental studies have identified the relationship between acoustic attenuation and petrophysical parameters:

$$\delta =\left(C.S/\phi \right).{\left(2\pi .k.f.{\rho }_f/\mu \right)}^{1/3}$$E1

with:

δ : attenuation (dB/cm), f : frequency (Hz), ρf : fluid density, μ : fluid viscosity (centipoise)

$\phi$: porosity, S : Specific surface (cm²/cm³), C : calibration coefficient, k : permeability (mD).

Figure 1 is an example of laboratory measurements on sandstone core plugs. The upper part of the figure shows the results obtained on cores with constant specific surface, the lower part on cores with variable specific surface, the specific surface being estimated on the basis of the average pore radius measurement. It is necessary for computing the permeability from eq 1, to measure the attenuation of the formation and calculate the effective specific surface of the formation.

Fabricius et al. [7] have found that the specific surface with respect to grain volume (Sg) is apparently independent from porosity. In an attempt to remove the porosity effect on Vp/Vs and mimic a reflected φ vs log (Sg) trend, they propose to use the following relationship between porosity φ, Vp/Vs and Sg

$$\log \left(Sg.m\right)=a.\phi +b.\left(Vp/Vs\right)+c\text{with}Sg=S/\left(1-\phi \right)$$E2

where it should be observed that Sg is multiplied by m to make Sg dimensionless. To establish eq 2, Fabricius et al. [7] have looked at ultra sonic data, porosity, and permeability of 114 carbonate core plugs.

Figure 1.

Relationship between attenuation and petro-physical parameters [6]. Laboratory measurements on cores with constant specific surface (top) and with variable specific surface (bottom).

An ideal porous medium made up of spheres (porosity 30%) of any size from cm (pisolites), mm (oolites) to a fraction of mm has an inverse trending specific surface that is the surface of solid exposed to the intergranular void. Fine grained well rounded sands or sandstones have a higher specific surface because the number of grains per unit volume is high. Clay content either as a coating of grains or as lumps or lenses enhances the specific surface up to a limit of effectiveness as long as pore connections prevail. The mineralogical nature of the grains, either silica (quartz or any siliceous rock debris) or carbonate is of little effect on the specific surface as long as secondary changes for instance diagenesis (for instance dolomitization), dissolution (karstification) or fracturation did not occur. Any process that simplifies the mineral surfaces exposed to the pore space tends to decrease the specific surface, whereas fractures and fissures, often occurring as networks or swarms, create new specific surfaces but very little added porosity.

Such processes remain minimal in tectonically quiet basins (Paris basin) where and when connate waters remain in equilibrium with the reservoir minerals. In subsurface exploration when no exposures more existing wells (logs, core and cutting information's) are available, it is difficult to assess the degree of transformation that a granular sedimentary reservoir formation has undergone through the elapsed geologic time.

2.2. To geophysical measurements

In practice, the parameter Ik-Seis (Indicator (I) of permeability (k) from acoustic or seismic (Seis) data) computed from equation 1 is proportional to permeability k.

$$\text{Ik-Seis = ( φ}{\text{.δ/S)}}^{\text{3}}{\text{/f = ( φ/SQ)}}^{\text{3}}\text{/f}$$E3

with f : P-wave frequency, Q quality factor, δ : attenuation,

S : specific surface, φ : porosity

The natural rocks are never perfectly elastic. The viscoelastic media always exhibit a wave amplitude decline as a function of time, independently of geometrical effects. From a wave propagation point of view, this results in attenuation on the one hand, and on the other hand in a dispersion of the propagation velocity, i.e. velocity depends on wave frequency. In order to conduct calculations and measurements on these parameters, it is necessary to build a model, the best known and used being the constant Q model. Q, called quality factor, is inversely proportional to the attenuation often designated by δ or α. In the case where Q is constant with frequency, which is true in the domain of frequencies used with seismic, we have: α = πf/Q, f being the frequency.

Wave propagation modelling in a saturated porous medium reveals the existence of three propagation waves: two P-waves and one S-wave. The two P-waves (P₁-wave and P₂-wave) have two different propagation velocities V_p1 and V_p2 (V_p1 > V_p2)and two different characteristic particle motions corresponding to a movement in which overall and fluid displacements are in phase for the P₁-wave (fast P-wave) and out of phase for the P₂-wave (slow P-wave) [8]. The laboratory experiments conducted by Plona [9] have confirmed the Biot 's theory [8]. Frequency is very important. Two major domains, separated by a critical frequency must be distinguished. Above the critical frequency, it is possible to estimate a permeability knowing that the calculated permeability is only an approach of hydraulic permeability. As shown in figure 2, above 2 kHz, permeability has an influence on velocities and attenuations [10]. The attenuation may reach a maximum for frequency in the order of 10 kHz, which is the domain of full waveform acoustic logging. It is the reason why several authors [2, 5-7, 10-12] attempted to predict permeability from acoustic data. The historical focus has been on predicting permeability from P-wave velocity and attenuation.

Figure 2.

Theoretical sensitivities of velocity and attenuation as a function of permeability, for a 19%-porosity sandstone (Schmitt, 1986). Curves 1 to 7 were calculated for a permeability of 2, 32, 200, 500, 1000, 1500 and 5000 milliDarcys. P₁-wave and P₂-wave are respectively the fast and the slow P-waves.

The transmission of an acoustic wave through geological formations is used for formation characterisation, in the acoustic frequency domain (ranging between 1 and 25 kHz). Acoustic logging allows the measurement of the propagation velocities and frequencies of the different waves (P-wave, S-wave and Stoneley wave) which are recorded by an acoustic tool. The analysis of the acoustic waves recorded simultaneously on both receivers of the acoustic tool is used to compute additional logs, defined as acoustic attributes, useful for the characterization of the formation, such as amplitude, shape index, wavelength and attenuation logs. Singular value decomposition (SVD) filtering method can be used to attenuate the noise, to measure the attenuation and to extract the acoustic wavelets [4].

Above the critical frequency, in the domain of acoustic frequencies, the Ik-Seis factor can be seen as a pseudo- permeability log after calibration on core data or hydraulic tests. Below the critical frequency (low frequency approximation), in the domain of seismic frequencies, the Ik-Seis factor can only be seen as a relative indicator which varies from 0 for less porous and permeable bodies to 1 for more porous and permeable bodies. In the same way, the specific surface is a pseudo specific surface which varies from 0 for less shaly bodies to 1 for more shaly bodies.

The seismic data must be inverted in order to obtain seismic impedance sections. If an elastic inversion is done, it is possible to obtain the elastic impedances Ip and Is. At well location, it is usually possible to obtain cross plots between acoustic impedance and porosity φ and to define a law between the two. Usually a linear or polynomial law can be extracted. The Ip, Is and φ quantities are used to compute the seismic specific surface (equation 2). If an acoustic inversion is done, well logs must be used to define an experimental law between Ip and Is. The analytic signal [13] is computed in order to extract, from the migrated seismic section, the variation of the seismic frequency and the Q factor versus time. The instantaneous frequency gives the frequency variation versus time and the envelop decrease leads to an estimation of the Q-factor. However a high signal to noise ratio is required. For that purpose, a singular value decomposition (SVD) filtering method [13] is used to enhance the coherent reflections and to attenuate the noise. Whatever the geophysical method (acoustic logging or reflection seismic surveying), equation 3 is used to obtain the Ik-Seis factor. The Ik-Seis factor can be used to detect in acoustic logging or on seismic sections permeable and impermeable bodies. For that purpose, we need to compute four quantities : P-wave frequency f, Q factor or attenuation δ, specific surface S and porosity φ. More information concerning the data processing and analysis is given in [4] and in the field examples.

3. Fields examples

2D and 3D seismic data were recorded in France at the boundary of the Meuse and Haute-Marne departments in the vicinity of the Andra Center (National radioactive waste management Agency). The acoustic data were recorded in well "EST 431" located on the Ribeaucourt township, in the national forest of Montiers-sur-Saulx, 8 km North-North-West of the Andra Center. Figure 3 indicates the location of 2D and 3D seismic surveys and of the selected seismic lines: line "07EST09" from the 2D survey and the in line 405 from the 3D survey. Figure 3 also indicates the location of the well " EST 431".

In a first step, we present and discuss the results obtained in acoustic logging. In a second step we describe those obtained with seismic data (2D then 3D).

3.1. Geological context

One of the drilling platforms, located in the centre of the studied zone, was used to study formations ranging from Oxfordian to Trias. The analysis presented here concerns borehole "EST 431" and covers the Oxfordian formation. The objective of this borehole is to complement the geological and hydrogeological knowledge of this formation. This formation consists essentially of limestone deposited in a vast sedimentary platform. The limestone facies, which vary from one borehole to another, are generally bio-detritic with reef constructions. The base of the Kimmeridgian shale was observed at –258.3 m (100 m ASL) and the base of the Oxfordian limestone at - 544.3 m (-186 m ASL).In this formation, porosity ranges between 5 and 20% and "porous horizons" of kilometric extension have been identified. As far as hydrogeology is concerned the observed water inflows are usually located in high porosity zones [14]. During the drilling, water inflows were detected at – 368 m and - 440 m. At the end of the drilling, the well was left in its natural water.

Figure 3.

Location map Well EST431 and seismic surveys. The selected seismic lines are the line 07EST09 from the 2D survey and the in line 405 from the 3D survey.

3.2. Acoustic logging

The acoustic tool used for the field experiment described in this paper is a flexible monopole tool with two pairs of receivers: a pair of near receivers (1 and 1.25 m offsets) and a pair of far receivers (3 and 3.25 m offsets). The data have been recorded through the far offset configuration. The sampling depth interval is 10 cm. The sampling time interval is 5 microseconds. The length of recording is 5 ms.

The acoustic log has been run in the Oxfordian carbonate formation, in the 333 – 510 m depth interval. Figure 4 (top left) shows the 3 m constant offset section, opposite the geological description. On the acoustic section, the refracted P-waves appear in the 0.6 – 1.2 ms time interval, the converted refracted shear waves in the 1.2 – 2 ms time interval, and the Stoneley wave in the 2 – 2.4 ms time interval.

On the acoustic section, we can differentiate: an event at 345 m showing a very strong attenuation of all the waves; an interval showing a very strong slowing down of the P and S waves (363 – 375 m); a relatively homogeneous mid-level interval (375 – 400 m); a level which stands out because of its strong variations in P, S and Stoneley velocities (400 – 455 m); a very homogeneous zone below 455 m with easily identifiable P and S waves, and an image of alteration between 501 and 507 m.

Figure 4.

Permeability estimation from acoustic logs [5] Top : acoustic section and P-velocity and frequency logs Bottom from left to right: porosity and attenuation logs; specific surface and predicted permeability (Ik-Seis) logs.

The processing of the acoustic data has been described in detail in [5]. Figure 4 is a display of acoustic logs: P-wave velocity (VP), P-wave frequency, acoustic porosity, P-wave attenuation, Sg specific surface and acoustic permeability. For the carbonate formation, a "V_P - V_S " cross plot led us to define a linear relationship between the two logs and to compute a shear velocity model. The experimental linear relationship computed as a regression line between "V_P and V_S" is V_S = 0.37 V_P + 879. The correlation coefficient between the V_S measurements and the V_S values given by the linear equation is equal to 0.93. For a clean formation, if the matrix and fluid velocities are known, an acoustic porosity log can be computed from the acoustic V_p velocities by using the formula given by Wyllie [15] expressed in velocities. As the maximum P-wave velocity is 6150 m/s at a depth of 464.5 m, the matrix velocity value has been chosen at 6300 m/s, and the fluid velocity at 1500 m/s, since the formation fluid is water. The acoustic porosity log, valid only in the clean part of the formation, shows a strong correlation (correlation coefficient : 0.86) with a NMR porosity log (not displayed here) recorded in the well.

The attenuation (expressed in dB/m) of the formation is computed from the first eigensection (obtained by SVD) of the refracted P-wave acoustic signal recorded by the two adjacent receivers of the acoustic tool. The results obtained by the SVD processing procedure are shown in figure 5.

Figure 4 (bottom right) shows the Sg specific surface log and the acoustic permeability log (Ik-Seis) calculated from eq. 2 and 3. The fluid viscosity μ and density ρ_f have been assumed to be constant (μ = 1 centipoise, ρ_f = 1 g/cm³). The Ik-Seis log detects three permeable zones at 368 m, between 400 and 440 m, and 506 m. The permeable zone located at 506 m corresponds to a high value of conductivity and is characterized by a low porosity (6 %), a 10 dB/m attenuation, but a significant decrease of the P-wave frequency and of the specific surface. During drilling, water inflows have been detected at – 368 m and - 440 m. At the end of the drilling, the well was left in its natural water. The hydraulic tests and conductivity measurements conducted later on did not confirm the inflow at 368 m seen during the drilling, but have validated the 400 – 440 m and 506 m permeable zones detected by the acoustic logging.

A short pumping test was conducted between the 4th and 7th of March, 2008. This test, associated with four geochemical logs, highlighted five productive zones between - 297 and - 507 m. The overall productivity is weak with a 15 L/min flow below 23.5 m drawdown. In general, there is a good correlation between the zones identified through the analysis of the geochemical logs, the natural gamma-ray log and the NMR porosity log. All the identified production zones correspond to low clay content zones.

The mean overall transmissivity of the Oxfordian formation below the EST431 drilling pad is 7.2 10^-06 m²/s, ranging between 5 10^-06 and 1 10^-05 m²/s. A total of five inflows have been identified between – 297.5 and -506 m. The most productive inflow (- 506 m) ranges between 2.5 10^-06 and 5 10^-06 m²/s. These inflows are associated with pore porosity and correspond in part to the porous horizons (HP) described in the Oxfordian formation below the Underground Laboratory:

• -297.5 to -301 m inflow, alternation of bioclastic and oolites packstone/mudstone, which corresponds to HP7,

• -328.5 m inflow, coral reefs packstone/grainstone, which corresponds to HP6,

• -413 m inflow, carbonated sand with oolites and oncolites, and numerous coral polyps, which corresponds to HP4, -439 m inflow, oolites and coral polyps grainstone-packstone interface, which corresponds to HP3,

• -506 m inflow, oolites and coral polyps grainstone-packstone interface, no correspondence with the porous horizons. The acoustic porosity and free fluid NMR porosity do not exceed 6%. The inflow, detected by the acoustic permeability log, is characterized by a 10 dB/m attenuation, a significant decrease of the P-wave frequency and of the specific surface.

The logs reveal a strong acoustic discontinuity at a depth of 345 m, clearly visible on the attenuation. The acoustic discontinuity is also revealed by a strong increase of the specific surface, a significant decrease of the acoustic velocities. The acoustic discontinuity is due to the presence of a thin shaly layer in the Oxfordian carbonate formation. It is confirmed by a change in the borehole diameter and a high value of the gamma ray log (not displayed here). Since at that depth, the acoustic porosity log has high values (larger than 25%, figure 4 bottom left), the shaly layer is probably water saturated.

The analysis of the acoustic waves recorded simultaneously on both receivers of the acoustic tool is used to compute additional logs defined as acoustic attributes useful for the characterization of the formation, such as amplitude, shape index and attenuation logs. The results obtained are optimum if the studied wave is extracted from the records and if the signal to noise ratio is high. We show the benefit of using Singular Value Decomposition (SVD) for that purpose [13]. The SVD processing is done on the 2 constant offset sections independently, in a 5 traces (N=5) depth running window. After flattening of each constant offset section with the picked times of the refracted wave, the refracted wave signal space is given by the first eigensection obtained by SVD:

$${\underset{\_}{\text{<munder accentunder="true">r <mo>\_</mo></munder>}}}^{\text{sig}}{\text{=λ}}_{\text{1}}{\text{<munder accentunder="true">u <mo>\_</mo></munder>}}_{\text{1}}{\text{<munder accentunder="true">v <mo>\_</mo></munder>}}_{\text{1}}{}^{\text{T}}$$E4

v₁ is the first singular vector giving the time dependence, hence named normalized wavelet, u₁ is the first singular vector giving the amplitude in depth, therefore called propagation vector and λ₁ the associated eigenvalue. The amplitude variation of the refracted wavelet over the depth interval is λ₁u_1.

Figure 5 (top) shows the normalized wavelet (v₁) and the associated amplitude (λ₁u₁) log versus depth, for the two constant offset sections associated with the two receivers (R1 and R2) of the acoustic tool. Figure 5 (bottom left) also shows the refracted wave signal space versus depth for the two receivers. The amplitude logs have been used to compute the attenuation log (figure 4, bottom left) expressed in dB/m. The correlation coefficient (figure 5, bottom right) between the two normalized wavelets has been computed at each depth. We can notice some anomalies at local depth (358, 390, 460, 492 and 503 m) and a significant decrease of the correlation coefficient in the 400 – 440 m depth interval. The interval corresponds to the porous and permeable zone detected by the Ik-Seis factor (figure 5, bottom right). It is therefore suggested that changes in phase or distortion of the acoustic signal is linked to propagation through a porous and permeable zone. The distortions can be measured by a shape index attribute. To measure the shape variation, an acoustic attribute, named Ic, independent of the energy of the source, has been introduced [10]. The Ic parameter is given by the following equation:

$$\text{Ic = }{\left(\left({\text{ A}}_{\text{2 }}{\text{+A}}_{\text{3}}\right){\text{/ A}}_{\text{1}}\right)}^{\text{n}}$$E5

where A₁, A₂ and A₃ are the amplitudes of the first three arches, respectively, of the studied signal and n an exponent.

The shape index is computed from formula 5 with an exponent value of 3 (n = 3). Figure 5 (top) allows the comparison between the two shape index logs (Ic-R1 and Ic-R2).

Figure 5.

Analysis of refracted waves in acoustic logging by SVD [4]. Top : SVD analysis (Wavelet, amplitude, shape index) on receivers R1 and R2. Bottom from left to right: signal space (R1, R2); predicted permeability (Ik-Seis), shape index and correlation coefficient logs.

The shape index logs highlight anomalic zones, in the 350 – 380 m depth interval with a maximum value at 358 m, and in the 400 – 440 m depth interval with a peak at 415 m which corresponds to the porous layer HP4. In order to reduce the noise to extract the common component of the two shape index logs (Ic-R1 and Ic-R2), the geometric mean of the two has been computed. The resulting shape index log Ic is compared with the Ik-Seis log and the correlation coefficient log (Figure 5, bottom right). The comparison shows a good coherence between the Ic log and the correlation coefficient log. The permeable and porous layers in the 400-440 m depth interval and the inflow at 506 m are seen both by the shape index log and by the Ik-Seis log. Shape index and correlation coefficient logs can be used as a quick look method to detect permeable bodies and inflows which must be confirmed by the Ik-Seis log.

3.3. Seismic surveying: 2D seismic line with acoustic impedance inversion after migration

The 2D seismic line was recorded in 2007 (see location map, figure 3). The 2D design is a split dip spread composed of 240 traces. The distance between 2 traces is 25 m. The source is a vibroseis source generating a signal in the 14 - 140 Hz frequency bandwidth. The bin size is 12.5 m. The nominal fold is 120.

3.3.1. Seismic procedure

A conventional seismic sequence was applied to the data set. It includes amplitude recovery, deconvolution and wave separation, static corrections, velocity analysis, CMP-stacking and time migration.

However, the amplitude recovery has been done in two main steps: spherical divergence compensation by using the Newman’s law [16], and attenuation compensation. Besides amplitude changes at each interface, the signal undergoes a general decay as a function of the propagation distance from the source (decaying as 1/r) and with its transmission through the interface. The application of a gain scalar that is a function of t can compensate for these effects. The commonly used gain function is:

$${\text{G}}_{\text{s}}\left(\text{ t}\right)\text{ = t}{\text{. V}}^{\text{2}}{}_{\text{rms}}\left(\text{t}\right){\text{/ V}}_{\text{1}}$$E6

where t is propagation time, V₁ is the propagation velocity of sound in the first layer and V_rms(t) is the rms velocity at time t.

In order to be able to apply the Newman’s law, the velocity model given by the rms velocity model must be estimated by velocity analysis. An a priori gain scalar t.V is applied to each shot point. The deconvolution is then performed by spectrum equalization in the 10-130 Hz frequency band. The wave separation by frequency – wavenumber filter (f-k filter) is then done to cancel the direct, refracted and surface waves and to enhance the reflected waves. The data are sorted in CMP gathers and the velocity analysis is performed. The knowledge of the velocity model allows the computation of the gain function G_s(t). The initial gain function is retrieved and replaced by the function G_s(t). After such a processing, a residual decay of the amplitude of the stacked traces has been noticed. The residual decay observed on the envelope of the stacked trace has been used to extract a residual compensation law G_r(t), after a strong smoothing in time of the envelope. The amplitude recovery residual law has been then approximated locally by an exponential law e^αt, where t is propagation time and α attenuation factor. Figure 6 (top, left) shows the migrated trace at CMP 600 after amplitude recovery by Newman’s law, the amplitude recovery function which is used to compensate the attenuation, and the seismic trace after compensation of attenuation.

The migrated section has been filtered by SVD in order to enhance the signal to noise ratio before computing seismic attributes (instantaneous frequency and envelope). Before SVD filtering, the migrated section has been shifted in time to flatten a reference seismic horizon, noted S1. In that case, the SVD processing is done on the migrated section, in a 5 traces (N=5) CMP running window and the signal space is composed of the two first eigensections in order to take into account the local dips which can been present on the time migrated section. Figure 6 (top, right) shows the amplitude recovery function G_r(t), the instantaneous frequency trace f(t) and the associated Q function. The Q factor at time t is obtained from the following equation:

$$\text{Q}\left(\text{t}\right)\text{ =}\pi \text{f}\left(\text{t}\right)\text{/}\alpha \text{ with }\alpha \text{= Ln }\left({\text{G}}_{\text{r}}\left(\text{ t}\right)\right)\text{/t}$$E7

After migration, a model-based stratigraphic inversion (a priori impedance model obtained from well data) provides a 2D impedance model section. The 2D impedance model section has also been shifted in time to flatten the seismic horizon S1. At well locations, the logs of velocity V_p, density ρ and porosity φ have been used to define laws between porosity and acoustic impedance (φ vs Ip) and between V_p and I_p. The porosity vs impedance cross-plot, displayed in figure 6 (bottom right) was used to define a linear law between the two. The porosity law obtained in the Oxfordian limestone and the relation between V_s and V_p (V_s = 0.37 V_p + 879) defined at well EST 431 allow the computation of both the pseudo specific surface and Ik-Seis functions, thanks to equations 2 and 3 (figure 6, bottom left). The processing is done on each migrated CMP independently in order to build the Ik-Seis section. The results are shown in the vicinity of CMP 600 in figure 7 and for the total line in figure 8.

Figure 6.

Seismic analysis at CMP 600. [4] Top from left to right: seismic trace before and after amplitude recovery, amplitude recovery function, instantaneous frequency trace and Q seismic trace. Bottom from left to right: acoustic impedance trace, Specific surface and Ik-Seis traces, porosity – acoustic impedance relationship.

3.3.2. Seismic analysis at CMP 600

Figure 6 (bottom left) displays, at location of CMP 600, the acoustic impedance trace, the pseudo specific surface seismic trace with its associated Ik-Seis trace. The specific surface seismic trace clearly shows the main geological units: the Kimmeridgian marls (between 0.1 and 0.2 s), the Oxfordian carbonates (between 0.2 and 0.33 s), the Callovo Oxfordian claystone (between 0.33 and 0.42 s), the Dogger carbonates (between 0.42 and 0.53 s) and the Toarcian claystone after 0.53 s. The shaly units have a high pseudo specific surface and a low value of the Ik-Seis factor. Figure 7 shows the migrated section and the associated Ik-Seis section in the vicinity of the CMP 600. As far as hydrogeology is concerned the observed water inflows are usually located in high porosity zones located in the lower part of the Oxfordian limestone, as it can be seen on the Ik-Seis section between 0.28 and 0.33 s (Figure 7, left). The Ik-Seis section also shows the distribution of the permeable bodies in the Dogger formation (Figure 7, right).

Figure 7.

Migrated and Ik-Seis sections at the vicinity of the CMP 600. [4] Left: 150 – 400 ms time interval (Oxfordian carbonate formation), Right: 350 – 600 ms time interval (Dogger carbonate formation).

3.3.3. Seismic line analysis

Figure 8 shows the distribution of the permeable bodies in the carbonate formations via the Ik-seis factor.

Several carbonated shelves have been emplaced during the middle and upper Jurassic, according to second order stratigraphic sequences. The Callovo-Oxfordian claystone is deposited over the Bathonian platform and is in turn overlaid by the Oxfordian platform.

Figure 8.

Ik-Seis section [4] 1: base of the Callovo-Oxfordian claystone; 2: carbonated layer used for flattening (S1); 3: top of the Callovo-Oxfordian claystone; Interval 3-4: Oxfordian limestone

The contact between the Dogger formation and the Callovo-Oxfordian claystone is sharp, associated to retrogradation condensed facies. The Dogger carbonate formation displays discontinuous thin porous layers in the Figure 7, but the correlation with the two known thin porous layers cannot be assumed due to the seismic vertical resolution.

The Callovo-Oxfordian claystone displays very low values of the Ik-Seis factor, with no significant lateral variation in good agreement with the known lithology. A carbonated layer, corresponding to a sequence boundary, has been used for flattening. This layer is a known marker-bed in the Paris basin (RIO: repère inférieur oolithique).

The value of the factor in the lower part of the Oxfordian limestone (around 0.30 s) is low in the interval CMP 850 - 1600 according to the transition from high energy inner ramp carbonated facies toward outer ramp marly facies in the west direction. The porosity decreases drastically in the western facies.

3.4. Seismic surveying: 3D seismic section with elastic impedance inversion after migration and depth conversion

The 3D seismic survey was recorded in 2010. The 3D design is a cross spread. The active spread is composed of 12 receiver lines with 120 stations each. The source lines are perpendicular to the receiver lines. The receiver and source line spacing’s are respectively 80 m and 120 m. The receiver and source point spacing’s are 20 m. The source is a vibroseis source generating a signal in the 14 - 140 Hz frequency bandwidth. The bin size is 10 x 10 m². The nominal fold is 60.

A conventional seismic sequence was applied to the data set. It includes amplitude recovery, deconvolution and wave separation, static corrections, velocity analysis and pre-stack time migration. The pre-stack time migration allows a model-based inversion (a priori impedance model obtained from well data) which provides a 3D elastic impedance block. No well being located in the 3D area, additional 2D lines have been recorded to calibrate 3D data on 3 wells situated outside of the 3D area. After migration, the seismic block is depth converted using consistent velocity model thanks to Bayesian kriging [17]. The velocity model is also used for estimating a Ik-Seis model in depth.

3.4.1. Time to depth conversion

Geostatistical approach for time-to-depth conversion of seismic horizons is often used in many geo-modelling projects. From a geostatistical point of view, the time-to-depth conversion of seimic horizons is a classical estimation problem involving one or more secondary variables.The converted depth and associated uncertainty can be estimated using a kriging method which can be constrained by the well markers, velocity model and interpreted horizons. For the multilayer case, the kriging estimator should take into account all the relationships between horizons determined by the velocity model associated to each layer.

The more appropriate kriging method for this problem is the Bayesian Kriging (BK) [17 -19]. Bayesian approach provides an excellent estimator which is more general than the traditional kriging with external drift(s) and fits very well to the needs for time-to-depth conversion of seismic horizons. The advantage of BK as estimator comparing to the others consists in the fact that we can manage simultaneously the uncertainty on the trend velocity model and the local uncertainty defined by the uncertainty of interpreted time maps and local fluctuations of interval velocities.

The input information's for BK are:

• Two-way-time (TWT) maps for interpreted horizons

• Well markers for each horizon

• Prior velocity model and associated uncertainty for each layer

• Local uncertainty definition for each time map (picking uncertainty, and spatial variogram definition)

• Local uncertainty definition of interval velocity for each layer (local velocity fluctuations around the trend model, and spatial variogram definition)

As any Kriging based estimator, the Bayesian Kriging provides as results:

• The estimated variable (estimated depth for each horizon)

• Variance of estimation (associated uncertainty of estimated depth).

The use of BK in depth conversion has the advantage to combine the prior knowledge of the velocity model with a certain degree of uncertainty and the well data. All sources of uncertainty (velocity and time) are integrated in a consistent way in a unique probabilistic model used for estimation or simulation.

For each selected horizon, the Bayesian Kriging provides its estimated depth Z associated with its time t. The “ Z versus t “ data set is interpolated in the whole space (3D block) at the time sampling rate (1 ms) in order to obtain a time to depth conversion model.

The time to depth conversion procedure is illustrated via the In line 405 extracted from the 3D seismic block (see location map, figure 3).

10 seismic horizons numbered from 1 to 10 have been picked in time and depth converted. The 10 seismic horizons are:

1. Top of Kimmeridgian White Limestones

2. Top of Porous Horizon HP4

3. Top of Lower Oxfordian (Top of target interval)

4. Top of Upper Callovian (RIO)

5. Top of Carbonated Dogger (Base of target interval)

6. Base of Argillaceous limestone and Longwy marls

7. Base of Carbonated Dogger

8. Top of Domerian

9. Base of Lias (base of Gryphees limestone)

10. Top of Beaumont dolomite

Figure 9 shows the acoustic impedance section (I_p) in time and the picked times of the 10 seismic horizons, the time- to – depth conversion model and the depth conversion of the 10 horizons. The time-to-depth conversion model can be used to convert in depth any type of seismic sections (amplitude, velocity, acoustic impedance, Ik-Seis factor).

3.4.2. Ik-Seis section in depth

The Ik-Seis section has been computed in time following the procedure described for the 2D line (figure 8) and then depth converted using the time to depth conversion model shown in figure 9.

Figure 10 shows the Ik-Seis sections in time (a) and in depth (b). The Ik-Seis sections allow the identification of a porous and permeable horizon (HP4) in the Oxfordian limestone (depth : 500 m, CMP : 100). They show the distribution of porous layers in the Dogger formation (between 800 m and 1000 m depth). The contact between the Dogger carbo nate formation and Callovo-Oxfordian argillite is clearly marked. The Callovo-Oxfordian claystone displays very low values of the Ik-Seis factor. The results obtained confirm the observations done on the 2D line (figure 8).

Figure 9.

Time to depth conversion. Top : acoustic impedance section in time and the picked times of the 10 seismic horizons, Bottom: time- to – depth conversion model and the depth conversion of the 10 horizons.

Figure 10.

Ik-Seis sections a : in time, b: in depth.

4. Conclusion

Knowledge about porosity and permeability is essential to evaluate fluid content and to detect fluid flow. We have presented a procedure which allows to use geophysical data (i.e. full waveform acoustic data and reflection seismic data) for a better understanding of the distribution of the porous and permeable bodies. The methodology is based on laboratory experiments which have shown that a formation permeability indicator, named Ik-Seis factor, can be obtained via the computation of 4 quantities: P-wave frequency, attenuation, porosity and specific surface.

The ability to use laboratory results for field geophysical applications has been discussed. The frequency content is important. Regarding frequencies above 2 kHz, permeability has an influence on velocities and attenuations. The attenuation may reach a maximum for frequency in the order of 10 kHz, this being the domain of full waveform acoustic logs.

Consequently, the procedure has been firstly conducted in acoustic logging to estimate the permeability of porous layers and to detect water inflows. Full waveform acoustic data were recorded in an Oxfordian carbonate formation. The Ik-Seis factor computed in the acoustic frequency domain (ranging between 10 and 25 kHz) has detected permeable zones, both associated with high porosity (20 %) but also with low porosity (6 %). The hydraulic tests and conductivity measurements conducted later on have validated the permeable zones detected by acoustic logging. The benefit of using the SVD method to evaluate the signal to noise ratio, to compute the attenuation log and to extract the acoustic wavelets from the acoustic data has been shown. It has also been observed that the correlation coefficient computed between acoustic wavelets recorded by two adjacent receivers of an acoustic tool significantly decreases in porous and permeable zones. It is therefore suggested that changes in phase or distortion of the acoustic signal is linked to propagation through a porous and permeable zone. The distortions can be measured by a shape index attribute. After calibration on core data or hydraulic tests, the Ik-Seis could be seen as a pseudo acoustic permeability log.

In seismic, after signal to noise ratio enhancement by the SVD method, processing is carried out in order to measure the needed parameters (frequency, attenuation, impedance) to compute the Ik-Seis factor. The analytic signal is used to compute the instantaneous frequency and attenuation (Q factor). The porosity and specific surface are computed from seismic impedances obtained by acoustic inversion of the migrated seismic sections. The Ik-Seis factor should only be used as a relative indicator which varies from 0 for less porous and permeable bodies to 1 for more porous and permeable bodies.

Our results suggest that it is possible to extract a significant Ik-Seis factor from seismic sections. This factor leads to a better understanding of the distribution of the porous and permeable bodies. The potential of the proposed procedure has been demonstrated via a 2D seismic profile and a seismic in-line extracted from a 3D block. We have demonstrated the benefit of combining time–to-depth conversion of seismic horizons by Bayesian kriging, consistent seismic velocity model and acoustic impedance in time for building a 3D geological model in depth. The field data example illustrates the potential of the proposed depth conversion procedure for estimating a Ik-Seis model in depth.

Acknowledgments

We thank Andra for permission to use the data presented in the field examples. We thank Jean Dellenbach, Jean Paul Sarda and Daniel Guillemot for their valuable help and advices. We thank Arben Shtuka (Seisquare) for very useful discussions on various occasions, specifically for his experience in the use of geostatistical methods for time to depth conversion.