Quantitative integration of seismic and GPR reflections to derive unique estimates for water saturation and porosity in subsoil

Pełen tekst

(1)GEOPHYSICAL RESEARCH LETTERS, VOL. 33, L05404, doi:10.1029/2005GL025376, 2006. Quantitative integration of seismic and GPR reflections to derive unique estimates for water saturation and porosity in subsoil R. Ghose1 and E. C. Slob1 Received 19 January 2006; revised 19 January 2006; accepted 24 January 2006; published 15 March 2006.. [1] For shallow subsoil, the estimates of in-situ porosity and water saturation are important, but until now it has been difficult to estimate them reliably. We relate seismic and GPR reflection coefficients to porosity and water saturation using a shared earth model. Using this model, we propose a method to integrate quantitatively seismic and GPR angledependent reflection coefficients. The new approach has been tested through numerical simulations. The results clearly show that from either seismic or GPR data alone it is impossible to obtain unique estimates for porosity and water saturation; however, a correct integration of those two data types leads to unique and stable estimates at a subsoil layer boundary. Potential applications of this approach exist in hydrogeology and environmental, agricultural and geotechnical engineering. Citation: Ghose, R., and E. C. Slob (2006), Quantitative integration of seismic and GPR reflections to derive unique estimates for water saturation and porosity in subsoil, Geophys. Res. Lett., 33, L05404, doi:10.1029/ 2005GL025376.. 1. Introduction [2] With improved hardware, acquisition and processing techniques, in recent years the resolution of shallow seismic methods has increased greatly, to an extent as to be able to image the very shallow depth range which is commonly addressed by GPR [e.g., Steeples et al., 1995; Ghose et al., 1998; Baker et al., 1999; Deidda and Balia, 2001]. Decimeter-scale seismic resolution in the top soil layers has become possible under favorable conditions. The seismic wavelength, thus, approaches the wavelength of GPR. Furthermore, the amplitude fidelity of the shallow seismic signals has improved; this has a positive influence on the results of amplitude-versus-offset (AVO) analysis of shallow seismic reflection data [Bradford et al., 1997; Ghose, 2003]. [3] For GPR data, AVO analysis has been attempted lately by many people [e.g., Zeng, 1998; Baker, 1998; Reppert et al., 2000; Bradford and Deeds, 2006]. Multioffset GPR data has been used to estimate interval velocity (Vem) of the high-frequency electromagnetic waves in the shallow subsurface. From this velocity, with some assumptions, it is possible to estimate for a partially saturated soil the volumetric water content, which is the product of porosity, f and water saturation, Sw [e.g., Greaves et al., 1996; Grote et al., 2002; Hubbard et al., 2002; Bradford and Harper, 2005; Lambot et al., 2004]. However, from. GPR data alone it is impossible to estimate f and Sw individually. [4] Empirical correlations between seismic and GPR (electromagnetic) velocities have also been recently made in order to deduce petrophysical properties [Koesoemadinata and McMechan, 2003]. Correlation of shallow seismic and GPR reflections has been attempted for improved interpretation [Cardimona et al., 1998; Bachrach and Rickett, 1999; Baker et al., 2001; Sloan et al., 2005], but these integrations are qualitative in nature. [5] Here we present the result of a quantitative integration of angle-dependent seismic and GPR reflection data, and suggest a means to estimate uniquely f and Sw in the shallow subsoil. Estimates of f and Sw and their lateral variation in the very shallow soil layers are crucial in many hydrogeological, environmental and geotechnical projects. A lack of reliable estimates for these two parameters is responsible for large risks and inefficiency in geoengineering design and constructions. First, we propose here a shared earth model to derive both elastic and electromagnetic bulk properties in terms of f and Sw, and the known values of the properties of the constituent solid, water and air phases. We present then a new approach for quantitative integration of seismic and GPR data. The approach is tested on synthetic data.. 2. Seismic-GPR Shared Earth Model [6] The seismic reflection coefficients at an interface between two layers are determined by the Lame´ constants (l, m) and the bulk density (r). The parameter that determines the electromagnetic (GPR) reflection coefficients is the electric permittivity, e. We assume magnetic permeability (mm) to be the free-space permeability (m0) and electric conductivity (s) to be constant; these are valid assumptions in most cases. Frequency-independence of s is generally a fair, first-order approximation in the frequency range of GPR [e.g., Davis and Annan, 1989]. [7] Therefore, l, m, r and e describe the bulk properties that determine the seismic and GPR reflection coefficients. Various mixing formulae have been theoretically and/or empirically developed so far for both seismic and dielectric responses of heterogeneous mixtures such as water-saturated soils. To estimate Lame´ constants for a 3-phase system (air-water-solid grains) we use the Bruggeman mixing rule, which is obtained by a theoretical effective medium calculation [Sihvola, 1999]. We express the bulk Lame´ constants in terms of the Lame´ constants of the constituent phases:. 1. Department of Geotechnology, Delft University of Technology, Delft, Netherlands.. ð1 fÞ. Copyright 2006 by the American Geophysical Union. 0094-8276/06/2005GL025376$05.00. L05404. Xs Xb Xw Xb Xa Xb þ Sw f þ ð1 Sw Þf ¼ 0; Xs þ 2Xb Xw þ 2Xb Xa þ 2Xb ð1Þ. 1 of 4.

(2) L05404. GHOSE AND SLOB: SEISMIC-GPR QUANTITATIVE INTEGRATION. L05404. where X stands for both l and m, the subscripts b, s, w and a denote the effective bulk property, and the properties of the constituent solid, water and air phases, respectively. The Lame´ constants for the solid, water and air phases are known [e.g., Landolt-Bo¨rnstein, 1982], or can be realistically assigned. The effective bulk values are obtained by solving equation (1) with an appropriate root finding routine. [8] The effective bulk density (rb) and the bulk dielectric permittivity (eb) are computed using a power law: Yba ¼ ð1 fÞYsa þ Sw fYwa þ ð1 Sw ÞfYaa ;. ð2Þ. where Y stands for e with a = 0.5 for sandy soils with little clay content (Complex Refractive Index Method or CRIM [Birchak et al., 1974]) and a = 0.65 for clayey soils [Wang and Schmugge, 1980], and Y stands for r with a = 1 (which gives the volumetric mean of the constituent densities); the subscripts have the same meaning as in equation 1. For eb, a number of models have been proposed before [e.g., Dobson et al., 1985; Jackson and O’Neill, 1986]. We have used CRIM as it provides accurate values for eb for many soil types. [9] In the absence of information about the geometrical distribution of the constituents, Hashin-Shtrikman (HS) bounds [Hashin and Shtrikman, 1962] offer the best (with narrowest range) upper and lower bounds for the predicted effective moduli. We have checked that the values of the effective bulk properties estimated from the mixing rules satisfy the HS bounds well. At low values of f, our estimates are close to the upper HS bound (determined by the solid phase), whereas as f increases the estimated values approach the HS lower bound (determined by the fluid phase). For Sw, the nature of variation is determined by the fluid components. The mixing models used in this study, therefore, give reasonable results. Figure 1 shows the effective seismic and electromagnetic velocities as functions of f and Sw, estimated using mixing rules for a 3-phase system. The velocity variations seen here are realistic, which indicates goodness of the models. [10] Note that in equations (1) and (2), the only free variables are f and Sw. Integration is achieved by using these two free variables for both elastic and electromagnetic properties.. 3. Integration Methodology [11] We assume that 3-component seismic and 2-component GPR data are acquired with identical acquisition geometry at a given location, the wavelengths are similar in the two data sets, and a shallow reflection event present in both data sets represents a given interface or depth region in the subsurface. CMP (common midpoint) measurements are needed. We also assume that the data have been preprocessed to minimize all surface-related effects, source-coherent and other noises, and then decomposed into SH, SV and P modes for seismic data and TE and TM modes for GPR data. [12] Using the effective bulk properties obtained by mixing rules as discussed above, the angle-dependent reflection coefficients, for both seismic and GPR, can be expressed as functions of f and Sw. We choose SH+TE and. Figure 1. Velocities from the mixing rules as a function of (a) porosity (f) and (b) water saturation (Sw). SV + TM mode combinations, because along a given measurement profile, seismic SH and GPR TE modes are identically polarized and they sense the subsurface heterogeneities in a similar way; the same is true for seismic SV and GPR TM modes. Assuming negligible influence of losses, seismic reflection coefficients are frequency-independent. The theoretical expressions for angle-dependent reflection coefficients are given, e.g., by Aki and Richards [1980] for seismic waves and by Kong [1986] for electromagnetic waves. [13] Finally, the integration is accomplished by combining the seismic and GPR reflection coefficients in a cost functional: X q;f. Cb ¼. X þ. q;f. b jRSH ðq; f; Sw Þ RSH data ðq; f Þj X b 2 jRSH data ðq; f Þj q;f b jRTE ðq; f; Sw ; f Þ RTE data ðq; f Þj X b 2 jRTE data ðq; f Þj q;f. !1b ;. ð3Þ. where q is the angle of incidence, f is frequency, RSH and RTE denote respectively the angle-dependent reflection coefficients of seismic SH and GPR TE wave modes estimated as functions of f and Sw using the mixing rules, TE RSH data and Rdata are the observed reflection coefficients for the same wave modes, and b = 2 for L2-norm (global rootmean-square error) and b = 1 for L1-norm (absolute error). The cost functional Cb is minimized to obtain estimates for in-situ f and Sw. If the data is noisy and the noise has roughly a zero mean then L1-norm should be used. The models for the reflection coefficients are adequate filters. Similar minimization can be performed also for SV + TM mode combination. Ideally, the minimization should be performed using all data (P + SV + SH + TE + TM), which will help with noise and erroneous amplitudes. For this purpose, equation (3) can be adapted. Here we use noisefree data, and hence, SH + TE and SV + TM combinations are examined separately. We invert for f and Sw using full, nonlinear Zoeppritz and Fresnel reflection coefficients. [14] Assuming that the bulk properties (or velocity and density) of the layer above are known from other processing. 2 of 4.

(3) GHOSE AND SLOB: SEISMIC-GPR QUANTITATIVE INTEGRATION. L05404. Table 1. Layer Parameters and Ranges in the Solution Space VP, m/s. VS, m/s. r, gm/cc. Vem, cm/ns. s, S/m. f, %. Sw, %. Layer 1 1100 110 1.8 10.4 0 30 50 (sand) Layer 2 1803 83 1.3 5.7 0.5 45 80 (clay) Layer 2 1044 – 2089 46 – 127 0.7 – 1.7 4.5 – 8.9 0.5 0 – 60 0 – 100 (sol. sp.). operations, then the integrated cost functional is convex and a unique solution is obtained for f and Sw of the layer below. We study this through numerical experiments.. 4. Numerical Simulations: Results [15] We have computed SH, SV and P-SV seismic reflection coefficients and TE and TM mode electromagnetic reflection coefficients for multiple source-receiver offsets using the formulations mentioned above. The synthetic data is free from all surface-related effects. The CMPsorted data is converted to incidence angle gathers. There are 13 traces in the 0 –45 incidence angle range. The model we have used is a 2-layered model representing an unsaturated loose sand layer overlying an almost saturated clay layer. For the upper layer (layer 1) we assume the soil properties to be known, while for the lower layer (layer 2) they are unknown. The used values of f and Sw for layer 2 are respectively 45% and 80%. In the solution space, the VP, VS, Vem and r for layer 2 are estimated using the mixing rules for a range of values for f and Sw. The layer parameters and the ranges used for the solution space are listed in Table 1. As the data is noise-free, here we have used only the dominant frequency; however in case of noise, using more frequencies will contribute to robustness. [16] Figure 2 shows the cost functional as a function of f and Sw when either seismic SH (Figure 2a) or GPR TE (Figure 2b) data are used (i.e., when only one of the two terms in equation (3) is present). The lines of local minima are shown in deep blue. Clearly, many combinations of f and Sw produce almost the same minimum in the cost functionals, and there is no clear global minimum. Angledependent reflection information does not help eliminate this nonuniqueness. Thus, from SH or TE reflection data. Figure 2. L2 norm cost functional for SH and TE modes as a function of f and Sw.. L05404. alone it is not possible to determine f and Sw uniquely and reliably. The same result is obtained when we use either seismic SV or GPR TM data alone. [17] The result of integration of seismic and GPR reflection coefficients is shown in Figure 3. The integrated cost functional (equation (3)), as a function of f and Sw, is plotted. Figure 3a shows the SH + TE mode integration, Figure 3b shows the SV + TM mode integration. We find a well-defined, single minimum in both cases. Figures 3a and 3b are very similar, indicating that the solution spaces for SV and TE modes separately resemble very much those for SH and TE modes, respectively (Figure 2). Because the elastic and electromagnetic properties behave in different ways as functions of f and Sw (Figure 1), the integrated cost functional is convex in the solution space, and this leads to stable and unique estimates of f and Sw. In this integrated model inversion, the local minima lines shrink to a sharp global minimum. No global search method is needed. The global minimum can be found using a local optimization routine; we have used the Matlab routine ‘fminsearch’ for this purpose. In Figure 3, the 45% f and 80% Sw for the second layer can be estimated correctly. We have checked with a maximum 5% error in the values of the solid constituent properties (which are least known) and found that the estimated f and Sw values remain unchanged. This indicates the stability of the inversion approach. [18] In this research we have considered only two variables, f and Sw; the reconstruction methodology will still work if more complexity is introduced by increasing the number of free variables, possibly at the cost of introducing local minima in the 3- or more-dimensional space. Further, one could also consider estimating f and Sw for both upper and lower layers. These issues are open for future research. It is important to have as many source and receiver components (3-C seismic, 2-C GPR) as possible. Full decomposition of the acquired multi-component data into up- and downgoing waves of various modes (P, SV, SH, TE, TM) is mandatory. To implement this approach on field TM. Figure 3. Integrated cost functional for SH + TE and SV + TM mode integration.. 3 of 4.

(4) L05404. GHOSE AND SLOB: SEISMIC-GPR QUANTITATIVE INTEGRATION. data, the main challenges are acquisition of high-quality CMP data in identical condition for both seismic and GPR, and correct data processing to obtain reliable reflection amplitudes. For this purpose, recent developments in specialized field hardware appear promising. Works on development of amplitude-preserved processing schemes to minimize surface-related effects, source-coherent and other noises, and correct wavefield decomposition are now in progress. Use of multiple frequencies and many wave modes in the integrated cost functional will assist further in tackling noise and amplitude errors.. 5. Conclusions [19] A new idea for quantitative integration of shallow seismic and GPR data is presented. A shared earth model relates seismic and GPR angle-dependent reflection coefficients to f and Sw. This model is based on mixing rules to obtain bulk seismic and electromagnetic properties from the properties of the solid, water and air constituents in soil. Integrating seismic and GPR reflection coefficients, we defined a cost functional with a convex solution space. The convexity is caused by the elastic and electromagnetic properties, and not by the mixing rules used. The integrated cost functional was minimized using a local optimization routine to obtain robust estimates of f and Sw. Numerical studies on a 2-layer model clearly show that if we use seismic or GPR data alone, then it is impossible to obtain unique estimates for f and Sw at the boundary between the two layers. However, if we combine these two data types, then we can estimate f and Sw uniquely and reliably. Although additional work must be done to apply this methodology to field data, the proposed approach of integration is powerful as it can provide laterally continuous, in-situ estimates of f and Sw at a shallow subsoil boundary. This will be vital in various fields, e.g., environmental and geotechnical engineering, hydrogeology, and agriculture. [20] Acknowledgments. We acknowledge the support provided by the Dutch Top Research School, ISES. Constructive reviews by John Bradford and Jan van der Kruk have improved the manuscript. The authors acknowledge the contribution of Jorien Schaaf at an early stage of this research.. References Aki, K., and P. G. Richards (1980), Quantitative Seismology: Theory and Methods, W. H. Freeman, New York. Bachrach, R., and J. Rickett (1999), Ultra shallow seismic reflection in depth: Examples from 3D ultra shallow survey with application to joint seismic and GPR imaging, paper presented at 69th Annual Meeting, Soc. of Explor. Geophys., Houston, Texas. Baker, G. S. (1998), Applying AVO analysis to GPR data, Geophys. Res. Lett., 25(3), 397 – 400. Baker, G. S., C. Schmeissner, D. Steeples, and R. G. Plumb (1999), Seismic reflections from depths less than two meters, Geophys. Res. Lett., 26(2), 279 – 282. Baker, G. S., D. W. Steeples, C. Schmeissner, and R. Plumb (2001), Nearsurface imaging using coincident seismic and GPR data, Geophys. Res. Lett., 28(4), 627 – 630. Birchak, J. R., L. G. Gardner, J. W. Hipp, and J. M. Victor (1974), High electric constant microwave probes for sensing soil moisture, Proc. IEEE, 62(1), 93 – 98.. L05404. Bradford, J. H., and J. C. Deeds (2006), Theory and application of groundpenetrating radar thinbed offset-dependent reflectivity, Geophysics, in press. Bradford, J. H., and J. T. Harper (2005), Wave field migration as a tool for estimating spatially continuous radar velocity and water content in glaciers, Geophys. Res. Lett., 32, L08502, doi:10.1029/2004GL021770. Bradford, J. H., D. S. Sawyer, and C. A. Zelt (1997), AVO analysis of low velocity, shallow sands (<50 m), paper presented at 67th Annual Meeting, Soc. of Explor. Geophys., Dallas, Texas. Cardimona, S. J., W. P. Clement, and K. Kadinsky-Cade (1998), Seismic reflection and ground penetrating radar imaging of a shallow aquifer, Geophysics, 63(4), 1310 – 1317. Davis, J. L., and A. P. Annan (1989), Ground-penetrating radar for highresolution mapping of soil and rock stratigraphy, Geophys. Prospect., 37, 531 – 551. Deidda, G. P., and R. Balia (2001), An ultrashallow SH-wave seismic experiment on a subsurface ground model, Geophysics, 66(4), 1097 – 1104. Dobson, M. C., R. T. Ulaby, M. T. Hallikainen, and M. A. El-Rayes (1985), Microwave dielectric behavior of wet soil: Part II. Dielectric mixing models, IEEE Trans. Geosci. Remote Sens., 23, 35 – 46. Ghose, R. (2003), AVO analysis of shallow seismic horizons: Effect of accuracy and uniformity of the effective source wavelet, paper presented at 73rd Annual Meeting, Soc. of Explor. Geophys., Dallas, Texas. Ghose, R., V. Nijhof, J. Brouwer, Y. Matsubara, and Y. Kaida (1998), Shallow to very shallow, high-resolution reflection seismic using a portable vibrator system, Geophysics, 63(4), 1295 – 1309. Greaves, R. J., D. P. Lesmes, J. M. Lee, and M. N. Toksoz (1996), Velocity variations and water content estimated from multi-offset, ground-penetrating radar, Geophysics, 61(3), 683 – 695. Grote, K., S. Hubbard, and Y. Rubin (2002), GPR monitoring of volumetric water content in soils applied to highway construction and maintenance, The Leading Edge, 21(5), 482 – 485. Hashin, Z., and S. Shtrikman (1962), A variational approach to the theory of the effective magnetic permeability of multiphase materials, J. Appl. Phys., 33(10), 3125 – 3131. Hubbard, S., K. Grote, and Y. Rubin (2002), Mapping the volumetric soil water content of a California vineyard using high-frequency GPR ground wave data, The Leading Edge, 21(6), 552 – 559. Jackson, T. J., and P. E. O’Neill (1986), Microwave dielectric model for aggregated soils, IEEE Trans. Geosci. Remote Sens., 24, 920 – 929. Koesoemadinata, A. P., and G. A. McMechan (2003), Correlations between seismic parameters, EM parameters, and petrophysical/petrological properties for sandstone and carbonate at low water saturations, Geophysics, 68(3), 870 – 883. Kong, J. A. (1986), Electromagnetic Wave Theory, John Wiley, Hoboken, N. J. Lambot, S., M. Antoine, I. van den Bosch, E. C. Slob, and M. Vanclooster (2004), Electromagnetic inversion of GPR signals and subsequent hydrodynamic inversion to estimate effective vadose zone hydraulic properties, Vadose Zone J., 3, 1072 – 1081. Landolt-Bo¨rnstein (1982), Numerical Data and Functional Relationships in Science and Technology: New Ser., Group V, Springer, New York. Reppert, P. M., D. F. Morgan, and M. N. Toksoz (2000), Dielectric constant determination using ground-penetrating radar reflection coefficients, J. Appl. Geophys., (43), 189 – 197. Sihvola, A. (1999), Electromagnetic Mixing Formulas and Applications, 284 pp., Inst. of Electr. and Electron. Eng., New York. Sloan, S. D., P. D. Vincent, G. P. Tsoflias, and D. W. Steeples (2005), Combining near-surface seismic reflection and ground-penetrating-radar data in the depth domain, paper presented at 75th Annual Meeting, Soc. of Explor. Geophys., Houston, Texas. Steeples, D. W., B. Macy, C. Schmeissner, and R. Miller (1995), Contrasting near-surface and classical seismology, The Leading Edge, 14(4), 271 – 272. Wang, J. R., and T. J. Schmugge (1980), An empirical model for the complex dielectric permittivity of soils as a function of water content, IEEE Trans. Geosci. Remote Sens., 18, 288 – 295. Zeng, X. (1998), Numerical modeling of GPR wavefields using ray-based, Fourier, and finite-difference algorithms with applications to field data, Ph.D. thesis, 107 pp., Univ. of Tex. at Dallas, Richardson. . R. Ghose and E. C. Slob, Department of Geotechnology, Delft University of Technology, P.O. Box 5028, 2600 GA Delft, The Netherlands. (r.ghose@tudelft.nl; e.c.slob@tudelft.nl). 4 of 4.

(5)