Detection of Porous and Permeable Formations: From Laboratory Measurements to Seismic Measurements

Résumé — Détection des formations poreuses et perméables : des mesures de laboratoire aux mesures sismiques — Dans le but d’avoir une meilleure compréhension de la distribution des corps poreux et perméables d’une formation géologique, nous montrons que de nouveaux attributs peuvent être extraits des données sismiques. Les attributs peuvent être également utilisés pour détecter les niveaux imperméables. La méthodologie est basée sur des mesures expérimentales effectuées en laboratoire, qui ont montré qu’un indicateur de perméabilité peut être obtenu à partir de quatre grandeurs: fréquence et atténuation des ondes de compression, porosité et surface spécifique. La procédure a d’abord été appliquée à des données de diagraphie acoustique dans le but d’estimer la perméabilité d’horizons poreux et de détecter des venues d’eau [Mari et al. (2011) Phys. Chem. Earth 36 , 17, 1438-1449]. En sismique, le traitement est réalisé dans le but d’estimer ces paramètres. Le signal analytique est utilisé pour calculer la fréquence instantanée et l’atténuation (facteur de qualité Q ). Porosité et surface spécifique sont estimées à partir des impédances sismiques obtenues par inversion des sections migrées. Les paramètres Abstract — Detection of Porous and Permeable Formations: From Laboratory Measurements to Seismic Measurements — We present a seismic processing method which shows 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. 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 4 input data: P-wave frequency and 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 [Mari et al. (2011) Phys. Chem. Earth 36 , 17, 1438-1449]. In seismic, 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 potential of the proposed procedure is demonstrated via a field case, both in full waveform acoustic and in seismic The example shows that the Ik-Seis factor can be used to map both the distribution of the permeable bodies in the carbonate formations and the non permeable shaly layers associated with the Callovo-Oxfordian claystone.


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 to seismic processing can be summarized in two main steps. The first step includes pre-processing of the data and the application of static corrections. The purpose of pre-processing is to extract reflected waves from individual shots, by filtering out the parasitic events created by direct and refracted arrivals, surface waves, converted waves, multiples and noise. It is intended to improve resolution, compensate for amplitude losses related to propagation, and harmonize records by taking into account source efficiency variations and eventual disparities between receivers. If the data are sorted in common midpoint gathers, the second processing step is the conversion of the seismic data into a migrated seismic section after stack. This second step includes the determination of the velocity model, using velocity analyses, the application of Normal Move-Out (NMO) corrections, stacking and migration. If the data are sorted in constant offset sections, a Pre-Stack Time (PSTM) or Depth (PSDM) Migration procedure which simultaneously performs dip correction, NMO correction, common mid-point stack and migration after stack, is applied. It is indispensable to have a good velocity model to carry out the migration process. The migrated section can then be transformed into acoustic impedance sections if well data (such as acoustic 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 final impedance model can be converted into porosity by using an empirical relationship between porosity and acoustic impedance established at well locations. 3D cube makes it possible to provide 3D imaging of the connectivity of the porous bodies . Core analysis are usually carried out to establish porosity versus permeability laws (Timur, 1968;Zinszner and Pellerin, 2007).
In this paper, we present a short review of the laboratory measurements conducted by Morlier and Sarda (1971). The possibilities of using laboratory results for field geophysical applications are then discussed. The methodology has been firstly conducted in acoustic logging to estimate permeability of porous layers and to detect water inflows . We present a seismic processing method which makes possible the extraction of attributes (P-wave frequency and attenuation, porosity and pseudo specific surface) from seismic sections and leads to a better understanding of the distribution of the porous and permeable bodies. Then, we show, with an example, how the methodology has been transposed to field data (acoustic logging and 2D seismic line). Morlier and Sarda (1971) have looked at ultra-sonic data (P-wave and S-wave velocities, frequencies and attenuations) and petrophysical data (porosity, permeability, specific surface) of numerous core plugs of different rock types (sandstone, limestone, carbonate). Their laboratory experiments have led them to the following results: -when there is only one saturating fluid, the attenuation is an increasing function of frequency f and of the reverse of the kinematic viscosity (ρf/μ with ρf: fluid density, μ: fluid viscosity (centipoise)); -the attenuation δ depends on the structure of the rock (i.e. pore geometry); -the attenuation δ can be expressed in terms of three structural parameters: porosity, permeability and specific surface.

Figure 1
Relationship between attenuation and petro-physical parameters (after Morlier and Sarda, 1971). Laboratory measurements on cores with a) constant specific surface and b) variable specific surface.
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 Equation (1) to measure the attenuation of the formation and to calculate the effective specific surface of the formation. Theoretically, the effective specific surface S can be calculated from the porosity ϕ and the Klinkenberg permeability k (given in m 2 in Eq. 2 but typically reported in mD) using Kozeny's equation (Kozeny, 1927). (2) (3) with ϕ: porosity, S: specific surface, Sg: specific surface with respect to grain volume, C k : Kozeny's factor. The Kozeny's factor can be calculated from the porosity via a simple model of linear 3D interpenetrating tubes (Mortensen et al., 1998). The specific surface Sg with respect to the bulk volume is given in 1/m in Equations (2) and (3) but typically reported in m 2 /cm 3 .
Fabricius et al. (2007) have found that the specific surface with respect to grain volume (Sg) apparently does not depend porosity. In an attempt to remove the porosity effect on Vp/Vs and mimic a reflected ϕ versus log (Sg) trend, they propose to use the following relationship between porosity ϕ, Vp/Vs and Sg: (4) where it should be observed that Sg is multiplied by m to make Sg dimensionless.
To establish Equation (4), Fabricius et al. (2007) have looked at ultra-sonic data, porosity and the permeability of 114 carbonate core plugs.

FROM LABORATORY MEASUREMENTS 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: 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. This is due to the fact that a portion of the vibration energy is transformed into heat because of various causes the viscoelasticity theory does not attempt to identify. This results directly from the second law of thermodynamics. From a mechanical point of view, this means that the stresses (σ) and strains (e) are no longer related by a constant (Hooke's law) but rather depend on time-varying elastic pseudo modulus m(t). This medium therefore has a memory: , where * is a convolution operator. 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 three best-known models are the following: standard linear or Zener model, nearly constant Q model, constant Q model. In the scope of this paragraph, we cannot elaborate further on the theory development, except to say that 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.
While the visco-elastic theory considers a homogeneous, non-elastic equivalent medium, the poro-elastic theory, developed by Biot (1956), differentiates between the fluid and solid. When taking into account the specific vibrations of solid and fluid, as the wave passes by, there is evidence of fluid displacements relative to the solid, hence an access to permeability. This theory uses a generalized Darcy's law with a complex permeability K(f) depending on frequency f with 1/K(f) = H 1 (f) + i.H 2 (f). The limit for f = 0 represents the hydraulic permeability. Coupling solid and fluid vibrations depends on a viscosity term, H 1 (f), and an inertia term, H 2 (f). Other coefficients, which depend on pore shape, can be computed using these terms. 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 1 -wave and P 2 -wave) have two different propagation velocities V p1 and V p2 and two different characteristic movements. It can be shown (Biot, 1956) that one of the characteristic movements corresponds to a movement in which overall and fluid displacements are in phase (P 1 -wave), and the second to a movement in which the displacements are out of phase (P 2 -wave). The P 2 -wave is called the slow wave or the wave of the second kind. This terminology derives from the fact that the associated velocity V p2 is much lower than the velocity V p1 of the in-phase movement wave (P 1 -wave), called wave of the first kind or fast P-wave. P 1 -wave and P 2 -wave correspond to classic P-waves, with which they merge in the absence of fluid. The laboratory experiments confirming the theory were conducted by Plona (1980). It is therefore possible for a given type of reservoir, to model and to compute some velocities (V p1 , V p2 , V s ) and some attenuations (Schmitt, 1986). More information concerning wave propagation in saturated porous media is given in Bourbié et al. (1987).
Frequency is very important. Two major domains, separated by a critical frequency, introduced by Biot (1956), 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 (top), above 2 kHz, permeability has an influence on velocities and attenuations. Figure 2 (bottom) represents the variation of attenuation as a function of the wave frequency and water saturation (Murphy, 1982). The fluid complement to water is air, the effect would be lowered with a heavier gas. 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 (Morlier and Sarda, 1971;Dupuy et al., 1973;Morlier, 1983, 2010;Klimentos and MacCann, 1990;Yamamoto, 2003;Fabricius et al., 2007; attempted to predict permeability from acoustic data. The historical focus has Acoustics and porous -permeable media. Top: 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, 1 000, 1 500 and 5 000 milliDarcys. P 1 -wave and P 2 -wave are respectively the fast and the slow P-waves. Bottom: Schematic representation of attenuation variations as a function of frequency and saturation (Murphy, 1982).
been on predicting permeability from P-wave velocity and attenuation. 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 which are recorded by an acoustic tool. Velocities and frequencies are computed from the picked arrival times of the waves (P-wave, S-wave and Stoneley wave). For a clean formation, if the matrix and fluid velocities are known, an acoustic porosity log can be computed from the acoustic P-wave velocity log. 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. We will show that the Singular Value Decomposition (SVD) filtering method is used to attenuate the noise, to measure the attenuation and to extract the acoustic wavelets.
Below the critical frequency (low frequency approximation), in the domain of seismic frequencies, Gassmann's formula (1951) is often used and 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 Sg with Equation (4) and then S with Equation (3). If an acoustic inversion is done, well logs must be used to define an experimental law between Ip and Is.
The analytic signal 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 is used to enhance the coherent reflections and to attenuate the noise.
Whatever the geophysical method (acoustic logging or reflection seismic surveying), Equation (5) 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 will be given in the field example. More information about Singular Value Decomposition (SVD) is given in Appendix (Glangeaud and Mari, 2000).

FIELD EXAMPLE: ACOUSTIC LOGGING AND SEISMIC SURVEYING
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 "EST431" 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 the seismic line "07EST09" and of the well "EST431".

Geological Context
One of the drilling platforms, located in the center of the studied zone, was used to study formations ranging from Oxfordian to Trias. The analysis presented here concerns borehole "EST431" 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 (Ferry et al., 2007). The limestone facies, which vary from one borehole to another, are generally bio-detritic with reef constructions. 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 (Delay et al., 2007). 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). 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.

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 source is a magnetostrictive transducer. The receivers are independent and each receiver has its own integrated acquisition device. 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. Location map: Well EST431 and seismic line 07EST09. Figure 4 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.

Characterization of a Carbonate Formation by Acoustic Logging
The processing of the acoustic data has been described in detail in . Figure 5 is a display of acoustic logs: P-wave velocity (VP), P-wave frequency, S-wave velocity (VS), Poisson's ratio, acoustic porosity, P-wave attenuation, Sg specific surface and acoustic permeability. For the carbonate formation, a "VP -VS" 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 "VP and VS" is VS = 0.37 VP + 879. The correlation coefficient between the measured VS log and the VS log 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 et al. (1956) expressed in velocities. It is given by the following equation: with V ma the matrix velocity, V f the fluid velocity.
As the maximum P-wave velocity is 6 150 m/s at a depth of 464.5 m, the matrix velocity value has been chosen at 6 300 m/s and the fluid velocity at 1 500 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 (Nuclear Magnetic Resonance) 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. This point will be discussed in the next section (Sect. 3.2.2). Figure 5 (bottom right) shows the Sg specific surface log and the acoustic permeability log calculated from Equations (4) and (5). The fluid viscosity μ and density ρ f have been assumed to be constant (μ = 1 centipoise, ρ f = 1 g/cm 3 ). The acoustic permeability 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 Permeability estimation from acoustic logs (after .
Top from left to right P-velocity and frequency logs; S-velocity and Poisson's ratio logs.
Bottom from left to right porosity and attenuation logs; specific surface and predicted permeability (Ik-Seis) logs.
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 below the EST431 drilling pad is 7.2 × 10 -06 m 2 /s, ranging between 5 × 10 -06 and 1 × 10 -05 m 2 /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 2 /s. These inflows are associated with pore porosity and correspond in part to the porous horizons described in the Oxfordian below the Underground Laboratory: --297.5 to -301 m inflow, alternation of bioclastic and oolites packstone/mudstone, which corresponds to HP7 (L2c); --328.5 m inflow, coral reefs packstone/grainstone, which corresponds to HP6 (L2b); --413 m inflow, carbonated sand with oolites and oncolites, and numerous coral polyps, which corresponds to HP4 (bottom of L1b); --439 m inflow, oolites and coral polyps grainstonepackstone interface, which corresponds to HP3 (top of L1a); -506 m inflow, oolites and coral polyps grainstonepackstone interface, no correspondence with the porous horizons (C3b). 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 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%, Fig. 5 bottom left), the shaly layer is probably water saturated.

Analysis of Refracted Acoustic Waves by Singular
Value Decomposition (SVD) 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. 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 (NS = 1) obtained by SVD: r v 1 is the first singular vector giving the time dependence, hence named normalized wavelet, u 1 is the first singular vector giving the amplitude in depth, therefore called propagation vector and λ 1 the associated eigenvalue. The amplitude variation of the refracted wavelet over the depth interval is λ 1 u 1. Figure 6 (top left) shows the refracted wave signal space versus depth for the two constant offset sections associated with the two receivers (R1 and R2) of the acoustic tool. The noise space is the stack of the eigensections from rank 2 to rank 5 (N = 5). The signal to noise ratio can be evaluated by: (8) and expressed in dB.
The signal to noise ratio logs associated with the two constant offset sections (R1 and R2) are shown in Figure 6 (top right). We can notice a poor signal to noise in the shaly layer at 345 m and two anomalies near 390 m. Figure 6 (bottom left) and Figure 7 (top) show for the two constant offset sections (R1 and R2) the normalized wavelet (v 1 ) and the associated amplitude (λ 1 u 1 ) log versus depth. The amplitude logs have been used to compute the attenuation log (Fig. 5, bottom left) expressed in dB/m. The correlation coefficient between the two normalized wavelets has been computed at each depth and the correlation coefficient log is shown in Figure 7 (bottom left). 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 (Fig. 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. The study of the shape variation of the acoustic signals is not new. Previous works have shown that this variation can be used to measure the attenuation of the formation and to detect criss-cross patterns, fractures (Mari et al., 1996) and permeable zones (Lebreton and Morlier, 1983;Lebreton and Morlier, 2010). To measure the shape variation, an acoustic Analysis of refracted waves in acoustic logging by Singular Value Decomposition (SVD). Top from left to right: signal space and signal to noise ratio on the two receivers (R1 and R2) of the acoustic tool. Bottom from left to right: wavelet and amplitude obtained by SVD, shape index, amplitudes of the first three arches of the wavelet (receiver R2). Analysis of refracted waves in acoustic logging by Singular Value Decomposition (SVD). Top from left to right: signal spaces (wavelet, amplitude, shape index) on receivers R1 and R2. Bottom from left to right: wavelet R1, wavelet R2, correlation coefficient log, predicted permeability (Ik-Seis) and shape index logs. attribute, named Ic, independent of the energy of the source, has been introduced (Dupuy and Morlier, 1973;Lebreton and Morlier, 1983). The Ic parameter is given by the following equation: where A 1 , A 2 and A 3 are the amplitudes of the first three arches, respectively, of the studied signal and n an exponent. Figure 6 (bottom left) shows the shape index Ic-R2 log computed from the amplitude logs (A 1 , A 2 and A 3 , Fig. 6 bottom right) associated with the normalized wavelet R2. The shape index is computed from formula (9) with an exponent value of 3 (n = 3). Figure 7 (top) allows the comparison between the two shape index logs (Ic-R1 and Ic-R2). 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 pick 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 (Fig. 7, 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.

Seismic Surveying
The 2D seismic line was recorded in 2007 (see location map, Fig. 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.

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 (Common Mid Point)-stacking and time migration.
However, the amplitude recovery has been done in two main steps: spherical divergence compensation by using the Newman's rule (1973) 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. Newman's rule is commonly used. The decay rule proposed by Newman is: 1/G s (t) = V 1 /(t.V 2 rms (t)) and consequently the gain function is: where t is propagation time, V 1 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 tV 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 8 (top, left) shows the stacked 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 in 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 in the time migrated section. Figure 8 (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: 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 Vp, density ρ and porosity ϕ have been used to define laws between Seismic traces before and after amplitude recovery Seismic analysis at CMP 600.
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.
porosity and acoustic impedance (ϕ versus Ip) and between Vp and Ip. The porosity versus impedance cross-plot, displayed in Figure 8 (bottom right) was used to define a linear law between the two. The porosity law obtained in the Oxfordian limestone and the relation between Vs and Vp (Vs = 0.37 Vp + 879) defined at well EST431 allow the computation of both the pseudo specific surface and Ik-Seis functions, thanks to Equations (4) and (5) (Fig. 8, bottom left). Figures 9 and 10 show the procedure of residual amplitude recovery and Q factor estimation for the seismic line "07EST09". One can see the acoustic impedance section (Fig. 9, top left) and the variations of the residual amplitude recovery law G r (t, n) along the seismic line, where t is propagation time and n CMP number (Fig. 9, top right). As quality control, the effect of amplitude recovery for attenuation can be clearly seen on the envelops of the migrated section before and after application of the G r (t, n) law (Fig. 9, bottom). After computation of the instantaneous frequency section, Equation (11) is used to compute the Q factor section (Fig. 10). Finally Equations (4) and (5) allow the computation of Ik-Seis sections. The results are shown in the vicinity of CMP 600 in Figure 11 and for the total line in Figure 12 Seismic analysis: amplitude recovery. a) Acoustic impedance seismic section. b) Amplitude recovery function section. c) Envelop seismic sections before amplitude recovery. d) Envelop seismic sections after amplitude recovery.
(between 0.42 and 0.53 s) and the Toarcien claystone after 0.53 s. The shaly units have a high pseudo specific surface and a low value of the Ik-Seis factor. Figure 11 shows the migrated section and the associated Ik-Seis section in the vicinity of the CMP 600. As far as hydrogeology is con-cerned 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 (Fig. 11a). The Ik-Seis section also shows the distribution of the permeable bodies in the Dogger (Fig. 11b). JL Mari and D Guillemot / Detection of Porous and Permeable Formations: From    Figure 12 shows the distribution of the permeable bodies in the carbonate formations via the Ik-Seis factor.

Seismic Line Analysis
Several carbonated shelves have been emplaced during the middle and upper Jurassic, according to second order stratigraphic sequences. The Callovo-Oxfordian argillite is deposited over the Bathonian platform and is in turn overlaid by the Oxfordian platform.
The contact between the Dogger formation and the Callovo-Oxfordian argillite is sharp, associated to retrogradation condensed facies. The Dogger formation displays discontinuous thin porous layers in Figure 12 but the correlation with the two known thin porous layers cannot be assumed due to the seismic vertical resolution.
The Callovo-Oxfordian argillite 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-1 600 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.

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 in the paper, 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. The procedure will be extended to 3D data.

APPENDIX: SINGULAR VALUE DECOMPOSITION (SVD)
Seismic or acoustic signals, recorded on an array of scalar sensors (singular component sensor) are generally described as summation of several events related to the different sources (reflected, refracted, converted waves, surface waves) propagating through the media. The observed scalar signal depending on time t and space x (a seismic record) is described as: ( 1) where a i (t) is the wavelet of source i, s i (t,x) the propagation vector of the source and * the symbol for the convolution product. Signal r(t,x) can be described in dual domains associated to time and distance variables as: • S i (f,x) = FT t [s i (t,x)], the distance -frequency space representation; • S i (f,k) = FT x [S i (f,x)], the frequency -wavenumber representation (2D Fourier transforms on time and distance variables).
After propagation through the media, the received signal r i (t) on sensor i results from the superposition of NS waves [a 1 (t), ..., a NS (t)] via the transfer functions s i, p (t): ( 2) where b i (t) is a noise supposed to be Gaussian, white and centred and N c the number of sensors. With signals sampled in time, we write the received signals in a data matrix as: The Singular Value Decomposition of the time-space data matrix r = provides two orthogonal matrices u and v and one diagonal matrix Δ = made up of singular values. The initial data matrix is expressed as: where: • u = [u 1 , ..., u k , ..., u N ] is a N c × N orthogonal matrix made up of left singular vectors u k giving the amplitude in the real case (amplitude and phase in the complex case), therefore called propagation vectors; • v = [v 1 , ..., v k , ..., v N ] is a N t × N orthogonal matrix made up of right singular vectors v k giving the time dependence, hence named normalized wavelets; • Δ = = diag(λ 1 , ..., λ k , ..., λ N ) a N × N diagonal matrix with the diagonal entries ordered λ 1 > ... > λ k > ... > λ N > 0.
The product u k v k T is an N c × N t unitary rank matrix named the kth singular image of data matrix r = . Therefore, r = is given by the sum of all the kth singular images multiplied by their correspondent kth singular values λ. The rank of the matrix r = is the number of non-zero singular values in Δ = . In the noise free case, if the recorded signals are linearly dependent (for example if they are equal to within a scale factor; that means one wave with an infinite velocity) the matrix r = is of rank one and the perfect reconstruction requires only the first singular image. If the N c recorded signals are linearly independent, the matrix r = is full rank and the perfect reconstruction requires all singular images.
Using the SVD filter, separation between the signal and the noise subspace is given by: The signal subspace r = sig is characterized by the first NS higher singular images (associated to the first NS higher singular values). It gives roughly the waveform of the dominant wave, its energy and its amplitude repartition on the sensors. The reminder subspace r = noise contains the waves with a low degree of sensor-to-sensor correlation and the noise. In practice, before performing the SVD filtering, a flattening operation on the initial data is applied to obtain an infinite apparent velocity for the selected wave. For refracted waves, the flattening pre-processing is obtained by time shifting the data, the time shifts are derived from the picking of the first arrival times.