Using a Poroelastic Theory to Reconstruct Subsurface Properties: Numerical Investigation

value the ﬂuid viscosity The Biot propagative, the displacement P-wave (PP), wave the P direct wave in the (PS1), S-wave (PS2) identiﬁed. two S-waves (PS1 and PS2) offsets. very small differences between both modelling The importance of seismic wave research lies not only in our ability to understand and predict earthquakes and tsunamis, it also reveals information on the Earth's composition and features in much the same way as it led to the discovery of Mohorovicic's discontinuity. As our theoretical understanding of the physics behind seismic waves has grown, physical and numerical modeling have greatly advanced and now augment applied seismology for better prediction and engineering practices. This has led to some novel applications such as using artificially-induced shocks for exploration of the Earth's subsurface and seismic stimulation for increasing the productivity of oil wells. This book demonstrates the latest techniques and advances in seismic wave analysis from theoretical approach, data acquisition and interpretation, to analyses and numerical simulations, as well as research applications. A review process was conducted in cooperation with sincere support by Drs. and


Introduction
The quantitative imaging of the Earth subsurface is a major challenge in geophysics. In oil and gas exploration and production, aquifer management and other applications such as the underground storage of CO 2 , seismic imaging techniques are implemented to provide as much information as possible on fluid-filled reservoir rocks. Biot theory (Biot, 1956) and its extensions provide a convenient framework to connect the various parameters characterizing a porous medium to the wave properties, namely, their amplitudes, velocities and frequency contents. The poroelastic model involves more parameters than the elastodynamic theory, but on the other hand, the wave attenuation and dispersion characteristics at the macroscopic scale are determined by the intrinsic properties of the medium without having to resort to empirical relationships. Attenuation mechanisms at microscopic and mesoscopic scales, which are not considered in the original Biot theory, can be introduced into alternative poroelastic theories (see e.g. Pride et al., 2004).
The inverse problem, that is, the retrieval of poroelastic parameters from the seismic waveforms, is much more challenging. Porosity, permeability and fluid saturation are the most important parameters for reservoir engineers. The estimation of poroelastic properties of reservoir rocks from seismic waves is however still in its infancy. The classical way of estimating these is to first solve the elastic problem and then interpret the velocities in terms of poroelastic parameters by using deterministic or stochastic rock physics modelling. However, unlike Full Waveform Inversion (FWI), these methods do not make full use of the seismograms.
In the poroelastic case, eight model parameters are required to describe the medium, compared with only one or two in the acoustic case, and three in the elastic case if wave attenuation is not considered. The advantages of using a poroelastic theory in FWI are (1) to directly relate seismic wave characteristics to porous media properties; (2) to use information that cannot be described by viscoelasticity or elasticity with the Gassmann (1951) formulae and (3) to open the possibility to use fluid displacement and force to determine permeability and fluid properties.
As an example of geological target, figure 1 presents the data recorded by a seismic survey on a seashore in the South of France. As water is pumped inland, saltwater is intruding into the coastal aquifer. This can affect the ground water and lead to severe problems with water supplies in the area. The monitoring of this phenomenon requires knowledge of the soil characteristics, including the permeability and porosity, and the properties of the fluid. A simple elastic approach cannot fully solve this problem, however the poroelastic theory may offer an alternative solution. In this example, the medium comprises alternating layers of sand, silt and clay with varying levels of compaction and a wide range of porosity and permeability. This layering produces strong reflected waves as shown in figure 1. The aim of the paper is to investigate how a poroelastic theory can be used to monitor water flow, and identify preferential pathways, using reflected waves. However, at this early stage, this chapter will only focus on numerical tests. Distance (m) Fig. 1. Example of data for which a poroelastic-based interpretation may be useful. Data are recorded on a seashore in the South of France. Source is a hammer shot, and the 24 receivers are equally spaced between 10 and 55 m from the source. The medium is very soft and water saturated, with a direct P wave velocity of c. 1600m/s and a S wave velocity lower than 200 m/s. We will investigate the gain and the feasibility of using a poroelastic approach, rather than the classical elastic one, in full waveform methods. The forward modelling is solved using different algorithms: a reflectivity approach, a 3D finite difference scheme and a 2D discontinuous Galerkin method. The comparison of synthetic data computed in the elastic and poroelastic cases shows that poroelastic modelling leads to some typical patterns that cannot be explained by elastic theory. This proves that the use of poroelastic theories may bring more insight to the model reconstruction, particularly, in relation to the fluid properties. Moreover, mesoscopic attenuation can be introduced in the poroelastic laws for double porosity medium, adding extra changes in the waveforms. This demonstrates the utility of using such theories to correctly reproduce measured seismic data.
Analytical formulae are then derived to compute the first-order effects produced by plane inhomogeneities on the point source seismic response of a fluid-filled stratified porous medium. The derivation is achieved by a perturbation analysis of the poroelastic wave equations in the plane-wave domain using the Born approximation. The sensitivity of the wavefields to the different model parameters can be investigated: the porosity, consolidation parameter, solid density, and mineral shear modulus emerge as the most sensitive parameters in the forward and inverse modelling problems. However, the amplitude-versus-angle response of a thin layer shows strong coupling effects between several model parameters.
The inverse problem is then tackled using a generalized least-squares, quasi-Newton approach to determine the parameters of the porous medium. Simple models consisting of plane-layered, fluid-saturated and poro-elastic media are considered to demonstrate the concept and evaluate the performance of such a full waveform inversion scheme. Numerical experiments show that, when applied to synthetic data, the inversion procedure can accurately reconstruct the vertical distribution of a single model parameter, if all other parameters are perfectly known. However, the coupling between some of the model parameters does not permit the reconstruction of several model parameters at the same time.
To get around this problem, we consider composite parameters defined from the original model properties and from ap r i o r iinformation, such as the fluid saturation rate or the lithology, to reduce the number of unknowns. We then apply this inversion algorithm to time-lapse surveys carried out for fluid substitution problems, such as CO 2 injection, since in this case only a few parameters may vary as a function of time. A two-step differential inversion approach allows the reconstruction of the fluid saturation in reservoir layers, even though the medium properties are mainly unknown.

Wave propagation in stratified porous media
The governing equations for the poroelastodynamic theory were first derived by Biot (1956), and are thus often referred to as "Biot's theory". The main hypothesis behind these equations is that the seismic wavelengths are longer than the pore size; the medium can then be described by homogeneised laws. Poroelastic theories have since been derived and improved by many authors (e.g. Auriault et al., 1985;Geertsma & Smith, 1961;Johnson et al., 1987).

Governing equations
Assuming a e −iωt dependence, Pride et al. (1992) rewrote Biot's (1956) equations of poro-elasticity in the frequency domain in the form where u and w respectively denote the average solid displacement and the relative fluid-to-solid displacement, ω is the angular frequency, I the identity tensor, ∇∇ the gradient of the divergence operator and ∇ 2 the Laplacian operator. The other quantities appearing in equations (1) are properties of the medium. The bulk density of the porous medium ρ is related to the fluid density ρ f , solid density ρ s and porosity φ: K U is the undrained bulk modulus and G is the shear modulus. M (fluid storage coefficient) and C (C-modulus) are mechanical parameters. In the quasi-static limit, at low frequencies, these parameters are real, frequency-independent and can be expressed in terms of the drained bulk modulus K D , porosity φ, mineral bulk modulus K s and fluid bulk modulus K f (Gassmann, 1951): It is also possible to link the frame properties K D and G to the porosity and constitutive mineral properties (Korringa et al., 1979;Pride, 2005): where G s is the shear modulus of the grains. The consolidation parameter c s appearing in these expressions is not necessarily the same for K D and G (Korringa et al., 1979). However, to minimize the number of model parameters, and following the recommendation of Pride (2005), we consider only one consolidation parameter to describe the frame properties. c s typically varies between 2 to 20 in a consolidated medium, but can be much greater than 20 in asoftsoil.
Finally, the wave attenuation is explained by a generalized Darcy's law which uses a complex, frequency-dependent dynamic permeability k(ω) defined via the relationship (Johnson et al., 1994):ρ In equation (5), η is the viscosity of the fluid and k 0 the hydraulic permeability. Parameter n J is considered constant and equal to 8 to simplify the equations.
The relaxation frequency ω c = η/(ρ f Fk 0 ),withF the electrical formation factor, separates the low frequency regime where viscous losses are dominant from the high frequency regime where inertial effects prevail. We refer the reader to the work of Pride (2005) for more information on the parameters used in this study.
The solution of equation (1) leads to classical fast P-and S-waves, and to an additional slow P-wave (often called Biot wave) . The fast P-wave has fluid and solid motion in phase, while the Biot wave has out-of-phase motions. At low frequency, the Biot wave has a diffusive pattern and can be seen as a fluid pressure diffusion wave. At high frequency, the inertial effects are predominant. This wave becomes propagative and can be seen in data, giving an experimental justification to the dynamic poroelasticity theory Plona (1980).

Mesoscopic attenuation and more complex theories
Although the slow P-wave does not appear on the seismograms at low frequency, it plays an important role in the attenuation process, as it produces loss of energy by wave-induced fluid-flow. However, the attenuation as described in the Biot theory is not strong enough to model the attenuation in geological media, especially at low (i.e. seismic) frequencies.
Attenuation processes can actually be separated into 3 different spatial scales, namely microscopic, mesoscopic and macroscopic (Pride et al., 2004). Within this classification, the Biot mechanism of attenuation takes place at macroscopic scale (on the order of the seismic wavelength). The microscopic attenuation is due to mechanisms that occur at the grain size, such as the squirt flow mechanism (Mavko & Jizba, 1991). This mechanism leads to wave attenuation mainly at high frequencies. The attenuation mechanism that prevails at low frequency comes from the mesoscopic scale (Pride et al., 2004), and it is due to fluid flow that occurs at boundaries between any medium heterogeneities whose sizes are between the grain sizes and the seismic wavelengths. This is particularly true for layered media (Gurevich et al., 1997;Pride et al., 2002) or when the medium contains 1) inclusions of different materials such as composite medium or double porosity medium (Berryman & Wang, 2000;Pride et al., 2004;Santos et al., 2006), or 2) different fluids (Santos et al., 1990) or patches of different saturation (Johnson, 2001).
Double porosity medium (DP in the text) refers to a porous medium which contains inclusions with different porosity and permeability (Pride et al., 2004;Santos et al., 2006). Assuming that the most compressible phase (patches, phase 2) is embedded into the least compressible one (host rock, phase1), the fluid flow inside the phase 2 can be eliminated from the homogenised equations. This assumption allowed Pride & Berryman (2003a, ,2003b to write the DP equations under the form of the classical Biot theory (eq. 1). This involves the use of complex frequency dependent moduli K U , C and M. Particularly, these parameters are functions of the respective volume of each phase and of the size of the patches (a denotes here the average radius of the inclusions in the host rock). Finally, as the patches are assumed to be spherical, the shear modulus of the medium is still real and not frequency dependent, and can be approximated by the geometrical mean of the modulus inside each phase. For the derivation and the detailed expressions of the parameters, please refer to the work of Pride et al. (2004) and Pride (2005).
To show how the mesoscopic attenuation due to DP media impacts the seismic properties, we look for the changes in the P-wave velocity and attenuation. The medium is composed of little patches of high permeable and high porosity in a less permeable host rock. Following the work of Liu et al. (2009), we consider a sandstone with 3% sand inclusions. The complex moduli K U , C and M are computed using the DP effective theory of Pride et al. (2004), leading to the P-wave velocity and attenuation with respect to the frequency. The results are compared to the seismic properties of each single phase and to the response using an average single porosity medium, where moduli are computed by geometrical averages. Figure 2 shows the P-wave velocity and attenuation (via the inverse of the quality factor) for the double porosity medium (for the inclusion radius a equal to 1, 5 and 10 cm), for the average single porosity medium and for both single phase media.
The P-wave velocity is much more dispersive for the double porosity medium than for the equivalent single porosity medium. At high frequency, it is much higher in the double porosity medium than in the equivalent single porosity medium. It shows two main changes:

137
Using a Poroelastic Theory to Reconstruct Subsurface Properties: Numerical Investigation www.intechopen.com  1) at the relaxation frequency of the host rock medium, which are at the transition between the diffusive and the propagative regime of the Biot wave, and 2) at low frequencies, around the relaxation frequency of the patches, with a strong dependance on a. It is worth noting that the size of inclusions, a, has an strong effect on the low frequency behaviour. We observe similar patterns for the P-wave attenuation (inverse of the quality factor). In the double porosity medium, attenuation shows two main peaks, associated with the two phases, while the equivalent single poroelastic medium produces only one single peak at high frequency. At low frequency (seismic frequencies), the attenuation is very high for the double porosity medium, leading to quality factors that are in good agreement with quality factors in geological materials. This pattern strongly depends on the size of the inclusions. On the other hand, attenuation produced by single phase medium is too low to be realistic. This means that the fluid flow at the boundaries between heterogeneities plays a fundamental role in the attenuation process, that cannot be neglected. These P-wave characteristics will have a strong influence on the seismic waveforms (see section 3.2).
Using a poroelastic theory is much more complex than an equivalent visco-elastic theory. However, modelling seismic waves with poroelastic theories take into account the attenuation induced by fluid equilibration at layer interfaces or heterogeneity boundaries, whereas a viscoelastic approach neglects this attenuation process. As shown by Pride et al. (2004), this is the most important attenuation process at low frequency. As the shallow subsurface has strong lateral and vertical heterogeneities, one should solve the full poroelastic theory to deal with attenuation.

Forward modelling solution
The model properties m, which are the material parameters introduced in the previous section, are nonlinearly related to the seismic data d via an operator f , i.e., d = f (m). The forward problem has been solved by many authors, using different methods. Analytical solutions have been derived for a homogeneous medium (Boutin et al., 1987;Philippacopoulos, 1997). The response of porous layered medium has been computed using reflectivity methods in the frequency-wavenumber domain, such as the Kennett (1983) approach (De Barros & Dietrich, 2008;Pride et al., 2002). This method was also used to solve the coupling between seismic and electromagnetic waves (Garambois & Dietrich, 2002;Haartsen & Pride, 1997). The poroelastic equations have been solved in 2D and 3D cases, mainly using finite difference schemes (Carcione, 1998;Dai et al., 1995;Masson & Pride, 2010;O'Brien, 2010) in the time-space domain. For discretisation issues, the equations (1) should be decomposed in propagative and diffusive parts, which have to be solved independently (Carcione, 1998). Other time domain numerical schemes have been used, such as finite elements (Morency & Tromp, 2008) or finite volume (de la Puente et al., 2008). Finally, Dupuy et al. (2011) solved this problem in the frequency domain using a discontinuous Galerkin approach. For a complete and precise review of the numerical modelling used to solve the poroelastic problem, we refer the reader to Carcione et al. (2010).
In this paper, we will use three different techniques to illustrate our points: www.intechopen.com by using the generalized reflection and transmission method of Kennett (1983). The synthetic seismograms are finally transformed into the time-distance domain by using the 3D axisymmetric discrete wavenumber integration technique of Bouchon (1981). • a Discontinuous Galerkin Method (DGM; Dupuy et al., 2011): For 2 dimension medium, the discrete linear system for the Biot theory has been deduced in the frequency domain for a discontinuous finite-element method, known as the nodal discontinuous Galerkin method. Solving this system in the frequency domain allows accurate modelling of the wave propagation for all frequencies.
The last two approaches are in the frequency domain, which has several advantages: 1) there is no need to decompose the problem into diffusive and propagative parts; 2) all frequencies, i.e., in the low-and high-frequency regimes, can be accuarately treated; 3) solving more complex theories, such as double porosity or poroviscoelastic theories is straighforward and does not require any modifications of the solver; and 4) frequency domain has been shown to be the most efficient way to solve the Full waveform inverse problem, as the solution has to be calculated only for a few frequencies (Pratt et al., 1998).

Seismic waveforms in poroelastic medium
As already stated, the main improvement of using poroelastic versus elastic theories is in the description of the attenuation from intrinsic medium parameters. Figure 3 gives an example of poroelastic and elastic data computed by equivalent finite difference codes (FD, O'Brien, 2010). The medium is a complex 200m thick reservoir embedded in a homogeneous half-space. The reservoir, modified from Manzocchi et al. (2008) is composed of 7 different facies, with different mineral and frame properties (see Fig. 3, top). In order to mimic a time lapse survey for CO 2 geological storage in a saline aquifer, a baseline is first computed for fully brine saturated medium. CO 2 is then injected in the center of the reservoir and spread out according to the permeability, leading to areas containing gas of roughly 500 m and 2000 mdiameter .
Figure 3 presents the differential data (data with CO 2 in the reservoir minus baseline) for both gas extensions. The same models are run using equivalent elastic properties (computed through the Gassmann formulation). As the velocities are equal, arrival times of the elastic and poroelastic waves are the same. However, changes appear in amplitude, mainly in the multiple reflected waves and coda. The amplitude differences are up to 40% of the data. It stresses the importance of using attenuation in the forward modelling and inversion processes. Even if the Biot theory cannot explain the full range of attenuation, it leads to significant changes in the seismic waveforms.
Poroelastic attenuation is due to the fluid movement. In the poroelastic theories, evidence for this lies in the relative fluid-to-solid motion and in the existence of the slow P-wave. Figure  4 reproduces the experiment performed by Plona (1980). The models are made using the Discontinuous Galerking method (DGM) of Dupuy et al. (2011) and are checked against the reflectivity approach (SKB) used by De Barros & Dietrich (2008). The SKB synthetic data have been corrected from the 3D effects, using an infinite line of sources, in order to be directly comparable to the DGM solutions. A P-wave is generated with an explosive source (central frequency of 200 Hz) in an quasi-elastic layer (the Biot wave is entirely diffusive in this layer) and is transmitted and converted into S-wave and Biot wave at the interface with a porous layer. Receivers are set into the second layer and record the three waves.  Fig. 3. Example of poroelastic and elastic numerical modelling to mimic differential data for CO 2 storage. Top) Resevoir models (left: Facies, right: Permeability) used in the modelling and modified from Manzocchi et al. (2008). Elastic properties are computed using the Gassmann relationships from the porous parameters. Middle) Left: Differential seismograms (model with CO 2 minus baseline) for elastic (red) and poroelastic data (black), and right: differences between the elastic and poroelastic differential data. Bottom) As the middle panels for a larger extension of CO 2 . Receivers and source are located on the free surface.
The explosive source has a 8Hz Ricker wavelet signature. The SKB solution is indicated by a continuous line and the DGM by crosses; dashed-dotted lines indicate the differences between the two solutions (multiplied by a factor of 5). PP, PBiot and PS2 stand for the transmitted P and the converted Biot and S waves, respectively. PS1 stands for the conical wave associated with the direct P wave in the first layer. The explosive source is in a quasi-elastic layer, while receivers stand in a porous half space.
Two cases are studied: 1) a low frequency case, where the source frequency is smaller than the cut-off frequency ( f c = 6400 Hz). The Biot wave, in this case, is not propagative and cannot be seen; 2) a high frequency case, the source frequency is higher than the cut-off frequency ( f c = 0.64 Hz). This is obtained by decreasing the value of the fluid viscosity by 10000. The Biot wave becomes propagative, and can be observed, mainly in the fluid displacement data. In the seismograms of figure 4, the transmitted P-wave (PP), the conical wave associated with the P direct wave in the first layer (PS1), and the transmitted S-wave (PS2) are identified. These two S-waves (PS1 and PS2) will be uncoupled at further offsets. Finally, the very small differences between both modelling methods prove their accuracy.
A similar case is studied in a double porosity medium. Taking the same source/receiver layout with an explosive source in the first layer comprising a quasi-elastic sandstone and a line of receivers in the second layer comprising the double porosity medium described in the part 2.2 (sandstone with 3% of 1 cm spherical sand inclusions), we compute the seismograms and compare the double porosity results with the effective single porosity results. The seismograms of solid and relative fluid/solid displacements are given in figure  5. The influence of the double porosity homogeneization (via complex frequency dependent mechanic moduli) is clearly visible on the transmitted P-waves. Particularly, the waveforms are strongly distorted as the attenuation and dispersion are higher (see figure 2). As we are in the low frequency domain, the Biot waves are not visible in the seismograms, but they are responsible for the loss of seismic energy. As predicted by the theory, the converted S-waves are not impacted by the double porosity approach.
In these examples, we have demonstrated the importance of taking into account complex poroelastic theories in order to understand and reproduce real seismic signals whose waveforms are strongly impacted by the presence of fluid and medium heterogeneities.

Sensitivity analysis
In the next sections, we use the reflectivity algorithm SKB (De Barros & Dietrich, 2008) and focus on backscattered energy, i.e., we consider reflected seismic waves as in a seismic reflection experiment. We further assume that, whenever they exist, waves generated in the near surface (direct and head waves, surface and guided waves) are filtered out of the seismograms prior to the analysis. The assumption of plane-layered media is admittedly too simple to correctly describe the structural features of geological media, but it is nevertheless useful to explore the feasibility of an inversion process accounting for the rheology of porous media.
The sensitivity of the seismic waveforms to the model parameters is investigated for layered medium by computing the first-order derivatives of the seismic displacements with respect to the relevant poroelastic parameters. These operators, which are often referred to as the Fréchet derivatives, are expressed via semi-analytical formulae by using the Born approximation (De Barros & Dietrich, 2008). They can be readily and efficiently evaluated numerically because they are only functions of the Green's functions of the unperturbed medium. In each layer, we consider the eight following quantities as model parameters: 1) the porosity φ,2 ) the mineral bulk modulus K s , 3) the mineral density ρ s , 4) the mineral shear modulus G s ,5 ) the consolidation parameter c s , 6) the fluid bulk modulus K f , 7) the fluid density ρ f and 8) the permeability k 0 . This parameter set allows us to distinguish the parameters characterizing the solid phase from those describing the fluid phase. The fluid viscosity η is one of the input parameter but it is not considered in the inversion tests as its sensitivity is similar to the permeability. than the porosity) on the wave amplitudes. The inversion for the parameters with such a low influence on the seismic waves will therefore be very delicate if other parameters are imperfectly known.
Figure 6 (right) shows the sensitivity of the P waves to the the fluid modulus K f and the mineral solid modulus K s . We note that the fluid modulus K f has a stronger influence than the solid modulus K s if the medium is poorly consolidated. The inverse is true for a consolidated medium. Similar patterns can be observed with the porosity: the higher the porosity, the stronger the influence of the fluid on the seismic waves. This means that it will be easier to determine fluid properties for an unconsolidated medium. For example, fluid substitution due to CO 2 injection leads to clear bright spot in Sleipner area (Norway, Arts et al., 2004), where the medium is poorly consolidated. The same set-up in stiff rocks does not produce such clear images, like in Weyburn field (Canada, White, 2009) Fig. 6. Left) Sensitivity of the reflected waves to the porous parameters for the P (downward triangles) and S (circles) waves, i.e., maximum of the energy reflected back by a perturbation of the porous parameters. Right) Sensitivity of the reflected waves to the fluid K f and the mineral bulk modulus K s as a function of the consolidation parameter cs. Note that higher consolidation parameter cs corresponds to softer materials.
To evaluate the coupling between parameters, we look at the Amplitude Versus Angle (AVA) curves in figure 7 for the PP (left) and SS (right) reflected waves due to a small and localized perturbation of a model parameter. We note that for some parameters, the model perturbations lead to similar modifications of the seismic response. For example, perturbations in densities and permeability show identical AVA responses. The same is true for the bulk moduli. This strong coupling between parameters will prevent simultaneous reconstruction of these parameters in an inversion process. Morency et al. (2009) also investigated the sensitivity of the seismic waves in porous media. They determined finite-frequency kernels based upon adjoint methods and investigated different parameter sets, in order to find the set that leads to the minimal coupling between parameters. They concluded that decomposing the input parameters into seismic velocities is the most stable approach in an inversion code.

Full waveform inversion
Full waveform inversion has shown to be an efficient and accurate tool to study the subsurface in the acoustic and elastic wave theory (Brossier et al., 2009). Historically, most of the FWI methods (Lailly, 1983;Tarantola, 1984) have been implemented under the acoustic approximation, for 2D model reconstruction (e.g. Gauthier et al., 1986;Pratt et al., 1998) or 3D structures (for instance, Ben-Hadj- Ali et al., 2008;Sirgue et al., 2008). Applications to real data are even more recent (Hicks & Pratt, 2001;Operto et al., 2006;Pratt & Shipp, 1999). The elastic case is more challenging, as the coupling between P and S waves leads to ill-conditioned problems. Since the early works of Mora (1987) and Kormendi & Dietrich (1991), the elastic problem has been addressed several times over the last years with methodological developments (Brossier et al., 2009;Choi et al., 2008;Gélis et al., 2007).
Using a poroelastic theory makes the problem even more difficult, especially because it adds much more unknowns. To the best of our knowledge, the first attempt to solve this Fig. 7. Energy of plane waves reflected from perturbations in ρ s , ρ f , k 0 , φ, K s , K f , G s ,andcs, as a function of incidence angle. The eight upper panels and eight lower panels correspond to PP and SS reflections, respectively. The curves are normalized with respect to the maximum value indicated above each panel.
problem was made by De Barros & Dietrich (2008) and De Barros et al. (2010) for stratified media and by Morency et al. (2009) and Morency et al. (2011) in 3-dimensional media. In the following sections, we will describe the main results obtained by De Barros & Dietrich (2008) and De Barros et al. (2010).

Inversion algorithm
Our method to determine the intrinsic properties of porous media is based on a full waveform iterative inversion procedure. It is carried out with a gradient technique to infer an optimum model which minimizes a misfit function. The latter is defined by a sample-to-sample comparison of the observed data d obs with a synthetic wavefield d = f (m) in the time-space domain, and by an equivalent term describing the deviations of the current model m with respect to an apriorimodel m 0 , i.e., where the L2-norms || . || D and || . || M are defined in terms of a data covariance matrix C D and an ap r i o r imodel covariance matrix C M Tarantola (1987). The model m contains the description of one or several parameters in layers whose thicknesses are defined by the peak content of the data (Kormendi & Dietrich, 1991). The model is updated using a quasi-Newton algorithm Tarantola (1987), which involves the Fréchet derivatives obtained earlier. As this problem is strongly non-linear, several iterations are necessary to converge toward an optimum model m, i.e, a model whose response d satisfactorily fits the observed data d obs .

Numerical results
In order to determine the accuracy of the inversion procedure for the different model parameters considered, we first invert for a single parameter, in this case the mineral density ρ s , and keep the others constant. The true model to reconstruct and the initial model used to initialize the iterative inversion procedure (which is also the ap r i o r imodel) are displayed in figure 8. The other parameters are assumed to be perfectly known. Their vertical distributions consist of four 250 m thick homogeneous layers. Parameters φ, c s and k 0 decrease with depth while parameters ρ f , K s , K f and G s are kept strictly constant.
Vertical-component seismic data (labelled DATA, fig 9) are then computed from the true model for an array of 50 receivers spaced 20 metres apart at offsets ranging from 20 to 1000 metres from the source. The latter is a vertical point force whose signature is a perfectly known Ricker wavelet with a central frequency of 25 Hz. Source and receivers are located at the free surface. As mentioned previously, direct and surface waves are not included in our computations to avoid complications associated with their contributions. Figure 9 also shows the seismogram (labelled INIT) at the beginning of the inversion, i.e., the seismogram computed from the starting model. Figure 8 shows that the true model, which consists in 10 metre thick layers from the surface to 1000 metre depth, is very accurately reconstructed by inversion. As there are no major reflectors in the deeper part of the model, very little energy is reflected toward the surface, which leads to some minor reconstruction problems at depth.
In figure 9, we note that the final synthetic seismograms (SYNT) almost perfectly fit the input data (DATA) as shown by the data residuals (RES) which are very small. The inversions carried out for the φ, ρ f , K s , K f , G s and c s parameters (not shown) exhibit the same level of accuracy. However, as predicted by De Barros & Dietrich (2008) and Morency et al. (2009) with two different approaches, the weak sensitivity of the reflected waves to the permeability does not allow us to reconstruct the variations of this parameter. Being related to seismic wave attenuation and fluid flow, permeability appears as not only the most difficult parameter to estimate but also the one which would have the greatest benefits to the characterization of porous formations, notably in the oil industry (Pride et al., 2003). One possibility to estimate it is to measure the fluid motion, or, by reciprocity, to use fluid pressure sources.
As observed in the sensitivity study, parameters are strongly coupled. Multiparameter inversion is thus an ill-posed problem, which is, in most of the cases, not reliable, as errors on one parameter will map into the reconstruction of the other parameters. The use of analytical expressions for the sensitivity kernels allows an easy rearrangement of the parameter set, in order to invert for the most pertinent parameters. Using some aprioriinformation, it is then possible to efficiently decrease the number of unknown parameters. For example, there is no reason to invert for both fluid parameters ρ f and K f , if we know that pores are filled by either gas or water. It is much more efficient to invert only for the saturation rate.

Differential inversion
To reduce the ambiguities of multiparameter inversion, a differential inversion has been considered and implemented (De Barros et al., 2010). Instead of dealing with the full complexity of the medium, we concentrate on small changes in the subsurface properties such as those occurring over time in underground fluid-filled reservoirs. This approach may be particularly useful for time-lapse studies to follow the extension of fluid plumes or to assess the fluid saturation as a function of time.
For example, the monitoring of CO 2 underground storage sites mainly aims at mapping the CO 2 extension. Time lapse studies performed over the Sleipner CO 2 injection site in the North Sea (see e.g. Arts et al., 2004) highlight the variations of fluid content as seen in the seismic data after imaging and inversion. In this fluid substitution issue, the parameter of interest is the carbon dioxide/saline water relative saturation. A differential inversion process will allow us to free ourselves from the unknown model parameters. This approach is valid for any type of fluid substitution monitoring problem, such as water-table variation, gas and oil extraction or hydrothermal activity.
The first step in this approach is to perform a base or reference survey to estimate the solid properties before the fluid substitution occurs. When performing a multiparameter inversion, the model properties are poorly reconstructed in general. However, the seismic data are reasonably well recovered. Thus, in spite of its defects, the reconstructed model respects the wave kinematics of the input data. In other words, the inverted model provides a description of the solid earth properties which can be used as a starting model for subsequent inversions. The latter would be used to estimate the fluid variations within the subsurface from a series of monitor surveys (second step).
To test this concept, we perform an inversion for two strongly coupled parameters, namely the porosity and the consolidation parameter. The resulting models are then used as starting models. We perturb the fluid properties of the true model to simulate a fluid variation over time. Two 30-metre thick layers located between 50 and 80 metres depth and between 110 and 140 metres depth are water depleted due to gas injection. The water saturation then varies between 60 and 80% in these two layers (figure 10). Our goal is to estimate the fluid properties by inverting the seismic data for the water saturation.
The model obtained is displayed in figure 10. We see that the location and extension of the gas-filled layers are correctly estimated. The magnitude of the water saturation curve, which defines the amount of gas as a function of depth, is somewhat underestimated in the top gas layer but is nevertheless reasonably well estimated. In the bottom gas layer, the inversion procedure only provides a qualitative estimate of the water saturation. These computations show that the differential inversion approach is capable of estimating, with reasonably good quality, the variations of fluid content in the subsurface without actually knowing all properties of the medium.

Conclusion
Using a poroelastic theory is much more complex than an elastic or a visco-elastic theory. However, poroelastic theories are an attempt to quantitatively describe the attenuation processes from the physical properties of the geological material. Furthermore, at seismic frequencies, attenuation is dominated by the mesoscopic scale mechanism, involving fluid flow at the boundaries of any heterogeneities. Poroelastic theories intrinsically take into account this loss of energy, while the equivalent visco-elastic approach neglects it. As the near surface media are stronly heterogeneous, with strong lateral and vertical contrasts, and different fluids involved, one has to deal with full poroelastic theories to accurately consider attenuation and fluid-solid motions.
The sensitivity to the different parameters varies hugely among parameters, and parameters are strongly coupled. Using a poroelastic theory to reconstruct model properties is in its nature an ill-conditioned problem. It shows however very promising possibility for differential inversion, and for certain issues where the problems can be reduced to the determination of only few parameters.
Poroelastic theories are, of course, not perfect yet, as they fail to give an universal law to explain seismic wave attenuation and propagation. They are however the direction to go if one wants to use full waveform inversion to make quantitative imagery of the rock physics and subsurface fluids. In particular, the permeability is a key parameter for exploration; the reconstruction of such a parameter from seismic waves will necessitate the use of complex poroelastic theories. Development in this direction still has to be continued. In particular, for imagery problems, data have to be improved to get around the problems of coupling between parameters. This can be done by using the information carried by fluid motions, which can give new insights into the permeability and the fluid properties, or by exploring in deeper details the coupling between seismic and electromagnetic waves.