- Email: [email protected]
3D seismic data for shallow aquifers characterisation
Journal of Applied Geophysics 68 (2009) 394–403
Contents lists available at ScienceDirect
Journal of Applied Geophysics j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j a p p g e o
3D seismic data for shallow aquifers characterisation Michela Giustiniani ⁎, Flavio Accaino, Stefano Picotti, Umberta Tinivella Istituto Nazionale di Oceanografia e Geofisica Sperimentale — OGS, Borgo Grotta Gigante 42c, 34010 Sgonico, Trieste, Italy
a r t i c l e
i n f o
Article history: Received 25 June 2008 Accepted 12 March 2009 Keywords: Aquifer Amplitude versus offset 3D pre-stack depth migration Depth modelling
a b s t r a c t We present the results obtained from conventional and non-conventional analysis of 3D high-resolution seismic data acquired nearby the water spring line, which separates the upper from the lower Friuli–Venezia Giulia plain (Italy), in order to characterise an important multilayered confined aquifer. The main targets of this study were two shallow aquifers located at about 30 m and 180 m depth, respectively. The aquifer structures were reconstructed by adopting a technique consisting an iterative updating procedure, for refining and improving an initial model in depth. The method includes pre-stack depth migration, residual move-out analysis and seismic reflection tomography. In the final 3D migrated cubes, two high velocity layers were identified at about 270 m and 480 m respectively, which correspond to unknown deep aquifers, as confirmed by recent well data (stratigraphies and down-hole velocity measurements). Travel-time tomography and Amplitude Versus Offset analysis evidence that seasonal variation in the seismic response of the aquifers are not detectable. However, in this case, aquifers are well detectable by lithological changes. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Fresh water is becoming one of the most important natural resource and many nations have carried out long term sustainable policies for the protection and public use of groundwater. The technological development has permitted the application of the seismic method to shallow hydrogeological studies starting from the late 1980s. In the last 20 years, in fact, the seismic reflection technique has been applied to many hydrogeological problems purposes with varying degrees of success (Steeples and Miller, 1990; Bachrach and Nur, 1998; Bradford, 1998; Whiteley et al., 1998; Cardimona et al., 1998; Bradford and Sawyer, 2002). Bachrach and Nur (1998) have shown clearly the relation between a deep water table and its seismic image, as well as the seismic response of the subsurface under different wetting and draining conditions. Whiteley et al. (1998) and Cardimona et al. (1998) have shown good examples of mapping shallow aquifers by reflection seismic data. Bradford and Sawyer (2002) have illustrated that pre-stack depth migration can significantly improve image quality and accuracy of hydrogeologic data. In this paper, we present the results of conventional and nonconventional processing of 3D high-resolution seismic data in order to characterise a multilayered confined aquifer system, located near the spring water line of the Friuli–Venezia Giulia plain (North-East Italy; Fig. 1). The main targets of this study were the identification of two shallow aquifers located at about 30 m and 180 m depth, respectively. The site is particularly suitable for the experiment ⁎ Corresponding author. Tel.: +0039 40 2140259; fax: +0039 40 327521. E-mail address: [email protected] (M. Giustiniani). 0926-9851/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2009.03.005
described in the following because of the presence of highly permeable sedimentary layers, which host an important aquifer system. The Friuli Plain is divided in two parts by the water spring line (Fig. 1): the upper and lower Friuli Plain. The upper part is characterised mainly by the presence of gravel deposits, which host unconfined aquifers. Instead, alternating layers of sand, gravel and clay are present in the lower part hosting important multilayered confined aquifers. In the study area, the Plio-quaternary sediments lie on the miocenic sediments, which are directly in contact with mesozoic carbonates. The latter were deposited in a shallow marine environment during an extensional phase. Cenozoic compressional phase caused the deposition of turbiditic sediments (Nicolich et al., 2004). During the winter 2005, three 2D seismic lines were acquired to define a preliminary subsurface geometry (Giustiniani et al., 2008), to analyze the seismic response of the aquifers and to optimise a subsequent 3D survey, which is the object of this paper (Fig. 1). Four 3D seismic datasets were acquired in March 2006 and July 2006, providing information about the 3D geometry of the aquifers and about possible variations between different seasons. The data processing was focused on the signal-to-noise ratio increase and vertical resolution enhancement. A ‘preserving-amplitude’ processing was adopted in order to allow a successive Amplitude Versus Offset (AVO) analysis (Yilmaz, 2001). An iterative updating procedure involving pre-stack depth migration, residual move-out analysis and seismic reflection tomography to obtain the velocity model and a realistic depth imaging of the investigated area, was applied. Finally, AVO analysis was performed to evaluate the Poisson's ratio contrast (Accaino et al., 2005), which provides useful information
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
395
Fig. 1. Location map of the 2D and 3D seismic surveys.
about the petro-physical properties of the shallow layers, such as the fluid content. 2. Seismic data acquisition Seismic prospecting is an effective tool for underground exploration, but the costs limited so far the use of this technology for
hydrological applications. Nowadays, however, the increasing water request and the technological improvements with the decreasing costs of the seismic acquisition are making the seismic prospecting suitable in environmental problems. In the studied area, the 2D geophysical investigations and well data showed strong lateral variations in petro-physical underground characteristics. In particular, the results obtained from a previous 2D
Fig. 2. 3D velocity model, obtained from the interpolation of three 2D velocity fields.
396
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
seismic data analysis (Giustiniani et al., 2008) have pointed out the presence of lateral variation of the seismic velocity that could be related to lithological and/or to pressure variation of the water filling the pores. The catchment wells drilled by the Acquedotto Basso Livenza reached the shallow aquifers at 30 m and 180 m. In particular, hydrogeological investigations performed in the Southern part of the studied area pointed out the absence of the aquifer located at 180 m, as confirmed by the previous 2D seismic investigation (Fig. 2; Giustiniani et al., 2008). The strong lateral velocity variation is evident in the 3D velocity field obtained from the interpolation of the 2D sections (Fig. 2), which is adopted in the present work as starting model for the updating procedure described in the following sections. The results obtained with the previous 2D seismic profiles were utilized to plan the subsequent 3D seismic acquisitions in two different areas called A (about 0.4 km2) and B (about 0.6 km2), as shown in Fig. 1. In each area, we carried out two surveys: the first one in March 2006 acquiring two cubes called respectively CubeA-spring and CubeB-spring, and the second in July 2006 obtaining the CubeAsummer and CubeB-summer. It is important to underline that the two cubes show differences in the coverage because the shot and receiver displacements are different. In fact, environmental problems strongly influenced the seismic acquisition. In conclusion, two 4D prospections surveyed the areas A and B. Preliminary field tests were performed to define the best sweep lengths and frequency windows of the vibroseis and to analyze the strong ground-roll and the environmental noise. We used the DMT Summit recording instruments with 260 active channels. The MiniVib IVI T-2500 was the seismic source and 10 Hz single geophones were deployed. Both the receiver interval and the receiver lines spacing (in-line) were 18 m. The source interval and the shot line spacing (cross-line) were 18 m and 36, respectively. Considering the receiver and shot intervals and the receiver and the shot line spacing, we adopted a binning equal to 9 × 9 m. The signal sweep of 15 s length and a frequency band of 20–200 Hz were used to provide 2 s length records after cross-correlation of recorded data with the pilot trace with a time sampling interval of 1 ms.
2.1. Processing of seismic data The quality of the raw data was satisfactory, despite of the groundroll (see Fig. 3). The reflections of interest are evident from 0.1 s to about 1 s, while the refractions are not always detectable because of velocity inversions caused by the alternation of impermeable and water filled permeable layers (Bourbié et al., 1987). The processing was focused to enhance the signal-to-noise ratio and to increase the vertical resolution. A ‘preserving-amplitude’ approach was adopted in order to maintain the control of the amplitudes of the signals. This is required for a successful AVO analysis. Firstly, we zeroed the traces or part of them showing a persistent level of noise. Then, we applied a band-pass filter 20–160 Hz to remove the low frequencies mainly related to the ground-roll. After applying a spherical divergence correction using the seismic velocity field obtained from the pre-stack depth migration, a predictive deconvolution was applied with a 120 ms operator length and a lag of 5 ms. A white noise level of 1% was added to stabilize the procedure. In the studied area, the lateral variation of lithological composition in the shallow subsurface layers causes strong shallow lateral velocity variations. Thus, we performed a surface consistent residual static correction using the velocity field obtained from the iterative updating procedure (described in the following sections). The effects of residual static application (window between 160 and 260 ms) is evident in Fig. 4. The obtained velocity fields, converted from interval to RMS velocity, were used to stack the data, after the application of a muting function to remove the direct arrivals. We applied a time-variant filter and F-X deconvolution to the post-stack data to increase the signal-tonoise ratio and the lateral continuity of the signals of the different interfaces. 3. Depth modeling of 3D seismic data The seismic image of an earth model in depth is usually described by two sets of parameters: the seismic velocities and the interfaces
Fig. 3. Example of raw data.
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
397
Fig. 4. Stack section of a 3D seismic line without application of residual statics (on the left) and with application of residual statics (right).
geometry. The estimation of these parameters with a required level of accuracy make the earth model definition a challenging task for the seismic exploration. Subsurface seismic imaging is finally realised by pre-stack depth migration, which, converts the seismic travel-time data into the depth section using the velocity field. The more reliable is the velocity field we supplied to the migration, the more realistic will be the imaging.
The modelling flow-chart adopted in this work is shown in Fig. 5a; it consists in an iterative updating procedure for refining and improving an initial model in depth, involving pre-stack depth migration, residual move-out analysis and seismic reflection tomography. The initial model generally consists in a horizontally layered model with laterally invariant velocities, or in a partially refined model obtained from other techniques and information.
Fig. 5. Flow-chart of a typical modelling in depth procedure (a) and the concept of residual move-out on a Common Image Gather (CIG) (b).
398
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
Residual move-out analysis is a velocity analysis performed after applying an initial velocity function to the data, in order to find the residual errors in the velocity field. If the initial velocity field has been estimated with sufficient accuracy, then the common image gathers (CIGs) derived from pre-stack depth migration using this model should exhibit a flat sequence of events, as shown in Fig. 5b. Any error in layer velocities and/or reflector geometries, on the other hand, should give rise to a residual move-out along those distribution of events on the CIGs, which are no-flat. In other words,
the degree of no-flatness of the reflection events on the CIGs is a measurement of the error in the model, and residual move-out analysis identifies the correction required to flatten the reflection events (Yilmaz, 2001). The adopted inversion technique can be described step by step in the followings points (Yilmaz, 2001). 1. Generate a set of CIGs from pre-stack depth migration using the velocity field derived from the initial model.
Fig. 6. 3D pre-stack depth migration of the CubeA-spring (a) and the CubeA-summer (b) using the final tomographic velocity model (c).
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
2. If the CIGs exhibit flat events, then the model is good and the procedure stops giving the final velocity model and the imaging performed by the pre-stack depth migration. If the events are not flat, the model is not sufficiently accurate yet, and the procedure goes on. 3. Compute the horizon-consistent residual move-out for all offsets along events on CIGs that corresponds to layer boundaries included in the model, obtaining the semblance sections. 4. Pick on the residual move-out profiles for all the horizons by tracking the semblance peaks. 5. Build the travel-time error vector Dt using the picked residual move-out in depth.
399
6. Define the initial model by a set of slowness and depth parameters and construct the coefficient matrix L in the following equation: Δt = LΔp
ð1Þ
where Δp is the tomographic update column vector. 7. Estimate the change in parameters vector Dp, by way of the Generalised Linear Inversion (GLI) solution given by −1 T T L Δt ð2Þ Δp = L L 8. Update the parameter vector p + Dp. By combining the update interval velocity profiles with the new depth horizons, we obtain the updated model.
Fig. 7. 3D pre-stack depth migration of the CubeB-spring (a) and the CubeB-summer (b) using the final tomographic velocity model (c).
400
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
9. Extract the velocity field from the new model and perform the pre-stack depth migration. 10. Return to the point 2 if the result is not satisfactory. At each iteration, both velocity and reflector geometries are updated, until the two set of parameters reach a good degree of stability. In our case, the initial model is the 3D interpolated velocity field obtained from the tomographic inversion of the 2D seismic lines described in Giustiniani et al. (2008; Fig. 2). The initial result was not satisfactory, because the seismic energy is not well focused and the CIGs exhibit non-flat sequence of events. We refined the model and performed the 3D pre-stack depth migration using the commercial package GeoDepth©. During each iteration, we carried out a residual move-out analysis, for all the horizons, by picking the coherences on the semblance profiles. The residual maps and the interval velocity sections extracted from the model for each depth horizon, are the input data for the tomographic update. The result of modelling updating needs to be verified for consistency with the input seismic data and examined for any remaining residual move-outs to decide whether or not to continue with the iterations of tomographic update (Yilmaz, 2001). According with the modelling flow-chart in Fig. 5a, a pre-stack depth migration is performed, using the velocity field derived from the new velocity model. The main consistency checks are the following: 1. Overlay the depth horizons from the updated model onto the image section derived from pre-stack depth migration and note if they
coincide with the reflectors associated with the layer boundaries included in the model. 2. Compute the residual move-out semblance spectra from the image gathers at selected locations after the model update and compare them with those before the update. The residual move-out analysis of image gathers and the update of velocity–depth model should be performed iteratively until the velocity–depth model and the depth image are consistent. Consistency may be achieved with only a few iterations when the residual move-outs are small. Fortunately, this is the case of our seismic survey, because we obtained a good consistency after only three iterations. After the third tomographic update, the events on the CIGs looks nearly flat and the depth horizons associated with the updated model coincide with the reflectors that correspond to the layer boundaries when superimposed on the image section: the model reached a good accuracy. Fig. 6c shows the final 3D velocity model of CubeA. The vertical resolution is satisfactory and we notice that two high velocity layers are evident at about 270 m and 460 m respectively, which could be associated to deeper aquifers. Fig. 6a shows the final migrated CubeA-spring, while Fig. 6b shows the migrated CubeA-summer. As it can be seen, below 250 m depth, two strong reflections appear in both cubes and can be associated to the top of two aquifers (see the high velocity in Fig. 6c). The cubes are sectioned to show the good lateral continuity of the reflectors. In Fig. 6c only one velocity cube is shown because the difference in the
Fig. 8. Correlation of the well stratigraphy with a transversal section of CubeA-summer. The velocity section is overlapped to the seismic migrated section. The dashed violet lines identify the picked horizons on the depth seismic section.
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403 Fig. 9. Stack section of 3D dataset acquired in March 2006 by using the traces with an offset ranging between 9 m and 149 m (panel a), between 150 m and 299 m (panel b) and greater than 300 m (panel c). The xline 45 is shown.
401
402
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
P-wave velocity values between the CubeA-spring and CubeAsummer velocity models is less that the estimated velocity error. Assuming for the first two acquifers an average picking error of about ±6 ms, the resulting average velocity deviation inside the two aquifers is about 2%. For the deeper acquifers, the velocity error increase and could be estimated to be about 5% for this inversion algorithm (Tinivella et al., 2002). In Fig. 7 the final velocity CubeB and the resulting migrated cubes are shown. We notice that the resolution and lateral continuity is better than the CubeA because of the higher coverage. Note also that the two reflections identified in the two CubeA at about 270 m and 480 m are here better visualized. Moreover, other two layers were identified and picked between 30 and 180 m. In the chaircut of the migrated CubeA-spring (Fig. 7a) a lateral amplitude inversion in the shallow layers is evidenced (with a circle). In the migrated CubeAsummer (Fig. 7b) it is not evident, but the coverage in this zone is lower with respect to the CubeA-spring. This lateral amplitude inversion may be due to lateral changes in the lithology or to a fault. 4. Correlation between seismic and well data On the basis of the seismic data results, a 510 m deep well was drilled, which allowed us to check the real consistency of our 3D models. Fig. 8 shows a transversal depth section cutting the CubeAsummer, overlapped to the corresponding interval velocity section, together with the well stratigraphy and the down-hole velocity log. There is an excellent correspondence between the horizons picked on the migrated section and the main discontinuities reported in the well stratigraphy, except for some slight discrepancies, which could be due to the distance between the well and the survey area (about 200 m, see Fig. 1). The two high velocity layers at about 270 m and 480 m correspond in the well stratigraphy to two deep aquifers. Those aquifers, because the catching wells were limited to depths of about 200 m, were unknown. Moreover, Fig. 8 shows the good correlation between the down-hole velocity log and the interval velocity section for aquifers 2 and 3. The four aquifers indicated on the well stratigraphy correspond to high interval velocity values, which is mainly related to the lithological variation. In fact, the velocity is affected mainly by three factors: 1) the matrix, 2) the pore fluid and 3) the pore pressure. In our case, as shown in Fig. 8, aquifer layers are composed by sand and gravel, while the adjacent layers are mainly composed by water saturated sandy clays. Note that, because of compaction effects, clay velocity increases with depth. It is well known that, for unconsolidated sediments, sandstones are characterised by higher velocities than those of shales (Schön, 1996). Moreover, the pore fluid is in overpressure conditions, slightly reducing the seismic velocity (Giustiniani et al., 2008). So, the combined effect is the increase of the seismic velocity in correspondence of the aquifers. 5. Amplitude versus offset The 3D pre-stack data analysis underlines the presence of AVO effects associated either with lithological changes and the presence of fluids at high pressure. These effects are evident, for example, in the three in-line stacked sections of the 3D dataset obtained by using different offset ranges: 0–149 m, 150–299 m, and higher than 300 m (Fig. 9). This figure shows that the amplitudes at the two main shallow aquifers change with offset. In the far-offset section, the amplitude is stronger compared to the others: this is evident for the deeper reflectors, at about 270 m and 480 m. As expected, the near-offset stacked section provides more information about the shallowest targets. To better discriminate between lithological and fluid changes at each interface, we performed the AVO analysis by using commercial software (Vern and Hilterman, 1995). We linearized Zoeppritz's
equations using the Aki–Richards method (Aki and Richards, 1980). The adopted approach is described in details in Giustiniani et al. (2008). The AVO inversion was applied to the 3D seismic dataset. The reliability of the result depends on the seismic data quality and on the coverage. As well as travel-times tomography (see previous section), also AVO analysis is not able to detect any seismic response difference between the two seasons. As already observed in the analysis of 2D seismic data (Giustiniani et al., 2008), the comparison between P- and S-wave reflectivities confirms that the top of each aquifer is well evident, confirming that the seismic interface is mainly generated by lithological changes and only slightly by the fluid condition changes. 6. Discussions and conclusions In this work we studied some strategies in order to characterise an important multilayered aquifer in an area of the NE of Italy by using the seismic method, namely seismic tomography, imaging and AVO. The 3D seismic acquisition was based on the results obtained from a previous 2D seismic survey and performed in two different seasons (March and July, 2006). We successfully applied a modelling in depth procedure, obtaining detailed interval velocity and migrated cubes for both seasons. We also applied the AVO analysis in order to obtain petro-physical properties of the multilayered aquifer. In this study, we obtained two significant results. The first is that, even if the imaging in depth of the aquifer system is satisfactory, travel-time tomography cannot detect any difference in the seismic response between the two seasons, because the velocity differences are less that the estimated velocity error (2% for the shallow aquifers and 5% for the deeper). This is also confirmed by the results obtained from the AVO analysis, which cannot detect overpressure variations between the two seasons. In aquifer studies, the seismic analysis (both tomographic and AVO inversion) detects mainly lithological variations, because the fluid pressure conditions in the aquifer slightly influences the seismic properties compared to the lithological effects. In fact, the stratigraphies of several wells located in the study area evidence that the aquifers are host in gravel layers confined between clay layers. Therefore, the 4D approach does not give other information with respect to the 3D seismic method, because travel-time tomography and AVO analysis are not enough sensitive to detect the overpressure variations inside the aquifers. The second result is that, because of the strong lithological variation in the study area, we obtained a detailed geometry and interval velocity model that allowed the discovery of two unknown aquifers, located at about 270 m and 480 m. In fact, a new 510 m deep well was drilled at about 200 m from the investigated area A (Fig. 1), confirming the existence of two aquifers identified by the seismic data analysis. Concluding, the 3D seismic method can be successfully applied in aquifer exploration. 3D cubes allow extending the well information, mainly lithological, to the surrounding area. This information can be used to optimise the location of new wells, or to detect deeper sources. Acknowledgments We are very grateful to Elvio Del Negro for his contribution and technical support in collecting the field data. We acknowledge the European community that has supported the project “Water-bearing characterisation with integrated methodologies – CAMI” (LIFE Project Number – LIFE04 ENV/IT/000500) and the leader of the project Daniel Gustavo Nieto Yabar. We wish to thank the Acquedotto Basso Livenza, and in particular Enrico Marin, for the logistical support and the all information provided (hydrogeological data and the stratigraphies of the catchment wells). We are grateful to the reviewers (Rinaldo Nicolich and Adele Manzella) for useful suggestions.
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
References Accaino, F., Tinivella, U., Rossi, G., Nicolich, R., 2005. Geofluid evidence from analysis of deep crustal seismic data (Southern Tuscany, Italy). Journal of Volcanology and Geothermal Research 148, 46–59. Aki, K., Richards, P.G., 1980. Quantitative Seismology: Theory and Methods. W.H. Freeman, San Francisco. Bachrach, R., Nur, A., 1998. High-resolution shallow–seismic experiments in sand, part I: water table, fluid flow, and saturation. Geophysics 63, 1225–1233. Bradford, J.H., 1998. Characterizing shallow aquifers with seismic reflection, part II — prestack depth migration and field examples. Geophysics 67, 98–109. Bradford, J.H., Sawyer, D.S., 2002. Characterizing shallow aquifers with seismic reflection, part II — prestack depth migration and field examples. Geophysics 67, 98–109. Bourbié, T., Coussy, O., Zinzner, B., 1987. Acoustics of Porous Media. Gulf Publishing, Houston, Texas. Cardimona, S.J., Clement, W.P., Kadinsky-Cade, K., 1998. Seismic reflection and groundpenetrating radar imaging of a shallow aquifer. Geophysics 63, 1310–1317. Giustiniani, M., Accaino, F., Picotti, S., Tinivella, U., 2008. Characterisation of shallow aquifers by high-resolution seismic data. Geophysical Prospecting 56, 655–666.
403
Nicolich, R., Della Vedova, B., Giustiniani, M., Fantoni, R., 2004. Carta del sottosuolo della Pianura Friulana. Litografia Cartografica, Firenze. Steeples, D.W., Miller, R.D., 1990. Seismic reflection methods applied to engineering and groundwater problems. In: Ward, S.W. (Ed.), Geothechnical and Environmental Geophysics 1, pp. 1–30. Schön, S., 1996. Physical properties of rocks—fundamentals and principles of petrophysics. Handbook of Geophysical Exploration: Seismic Exploration, vol. 18. Pergamon, Oxford/Tarrytown. Tinivella, U., Accaino, F., Camerlenghi, A., 2002. Gas hydrate and free gas distribution from inversion of seismic data on the South Shetland margin (Antarctica). Marine Geophysical Research 23, 109–123. Vern, R., Hilterman, F., 1995. Lithology color-coded seismic sections: the calibration of AVO crossplotting of rock properties. The Leading Edge 14, 847. Whiteley, R.J., Hunter, J.A., Pullan, S.E., Nutalaya, P., 1998. Optimum offset” seismic reflection mapping of shallow aquifers near Bangkok, Thailand. Geophysics 63, 1385–1394. Yilmaz, O., 2001. Seismic Data Analysis: Processing, Inversion and Interpretation of Seismic Data. SEG Series: Investigation in Geophysics, Tulsa.
Contents lists available at ScienceDirect
Journal of Applied Geophysics j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j a p p g e o
3D seismic data for shallow aquifers characterisation Michela Giustiniani ⁎, Flavio Accaino, Stefano Picotti, Umberta Tinivella Istituto Nazionale di Oceanografia e Geofisica Sperimentale — OGS, Borgo Grotta Gigante 42c, 34010 Sgonico, Trieste, Italy
a r t i c l e
i n f o
Article history: Received 25 June 2008 Accepted 12 March 2009 Keywords: Aquifer Amplitude versus offset 3D pre-stack depth migration Depth modelling
a b s t r a c t We present the results obtained from conventional and non-conventional analysis of 3D high-resolution seismic data acquired nearby the water spring line, which separates the upper from the lower Friuli–Venezia Giulia plain (Italy), in order to characterise an important multilayered confined aquifer. The main targets of this study were two shallow aquifers located at about 30 m and 180 m depth, respectively. The aquifer structures were reconstructed by adopting a technique consisting an iterative updating procedure, for refining and improving an initial model in depth. The method includes pre-stack depth migration, residual move-out analysis and seismic reflection tomography. In the final 3D migrated cubes, two high velocity layers were identified at about 270 m and 480 m respectively, which correspond to unknown deep aquifers, as confirmed by recent well data (stratigraphies and down-hole velocity measurements). Travel-time tomography and Amplitude Versus Offset analysis evidence that seasonal variation in the seismic response of the aquifers are not detectable. However, in this case, aquifers are well detectable by lithological changes. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Fresh water is becoming one of the most important natural resource and many nations have carried out long term sustainable policies for the protection and public use of groundwater. The technological development has permitted the application of the seismic method to shallow hydrogeological studies starting from the late 1980s. In the last 20 years, in fact, the seismic reflection technique has been applied to many hydrogeological problems purposes with varying degrees of success (Steeples and Miller, 1990; Bachrach and Nur, 1998; Bradford, 1998; Whiteley et al., 1998; Cardimona et al., 1998; Bradford and Sawyer, 2002). Bachrach and Nur (1998) have shown clearly the relation between a deep water table and its seismic image, as well as the seismic response of the subsurface under different wetting and draining conditions. Whiteley et al. (1998) and Cardimona et al. (1998) have shown good examples of mapping shallow aquifers by reflection seismic data. Bradford and Sawyer (2002) have illustrated that pre-stack depth migration can significantly improve image quality and accuracy of hydrogeologic data. In this paper, we present the results of conventional and nonconventional processing of 3D high-resolution seismic data in order to characterise a multilayered confined aquifer system, located near the spring water line of the Friuli–Venezia Giulia plain (North-East Italy; Fig. 1). The main targets of this study were the identification of two shallow aquifers located at about 30 m and 180 m depth, respectively. The site is particularly suitable for the experiment ⁎ Corresponding author. Tel.: +0039 40 2140259; fax: +0039 40 327521. E-mail address: [email protected] (M. Giustiniani). 0926-9851/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2009.03.005
described in the following because of the presence of highly permeable sedimentary layers, which host an important aquifer system. The Friuli Plain is divided in two parts by the water spring line (Fig. 1): the upper and lower Friuli Plain. The upper part is characterised mainly by the presence of gravel deposits, which host unconfined aquifers. Instead, alternating layers of sand, gravel and clay are present in the lower part hosting important multilayered confined aquifers. In the study area, the Plio-quaternary sediments lie on the miocenic sediments, which are directly in contact with mesozoic carbonates. The latter were deposited in a shallow marine environment during an extensional phase. Cenozoic compressional phase caused the deposition of turbiditic sediments (Nicolich et al., 2004). During the winter 2005, three 2D seismic lines were acquired to define a preliminary subsurface geometry (Giustiniani et al., 2008), to analyze the seismic response of the aquifers and to optimise a subsequent 3D survey, which is the object of this paper (Fig. 1). Four 3D seismic datasets were acquired in March 2006 and July 2006, providing information about the 3D geometry of the aquifers and about possible variations between different seasons. The data processing was focused on the signal-to-noise ratio increase and vertical resolution enhancement. A ‘preserving-amplitude’ processing was adopted in order to allow a successive Amplitude Versus Offset (AVO) analysis (Yilmaz, 2001). An iterative updating procedure involving pre-stack depth migration, residual move-out analysis and seismic reflection tomography to obtain the velocity model and a realistic depth imaging of the investigated area, was applied. Finally, AVO analysis was performed to evaluate the Poisson's ratio contrast (Accaino et al., 2005), which provides useful information
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
395
Fig. 1. Location map of the 2D and 3D seismic surveys.
about the petro-physical properties of the shallow layers, such as the fluid content. 2. Seismic data acquisition Seismic prospecting is an effective tool for underground exploration, but the costs limited so far the use of this technology for
hydrological applications. Nowadays, however, the increasing water request and the technological improvements with the decreasing costs of the seismic acquisition are making the seismic prospecting suitable in environmental problems. In the studied area, the 2D geophysical investigations and well data showed strong lateral variations in petro-physical underground characteristics. In particular, the results obtained from a previous 2D
Fig. 2. 3D velocity model, obtained from the interpolation of three 2D velocity fields.
396
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
seismic data analysis (Giustiniani et al., 2008) have pointed out the presence of lateral variation of the seismic velocity that could be related to lithological and/or to pressure variation of the water filling the pores. The catchment wells drilled by the Acquedotto Basso Livenza reached the shallow aquifers at 30 m and 180 m. In particular, hydrogeological investigations performed in the Southern part of the studied area pointed out the absence of the aquifer located at 180 m, as confirmed by the previous 2D seismic investigation (Fig. 2; Giustiniani et al., 2008). The strong lateral velocity variation is evident in the 3D velocity field obtained from the interpolation of the 2D sections (Fig. 2), which is adopted in the present work as starting model for the updating procedure described in the following sections. The results obtained with the previous 2D seismic profiles were utilized to plan the subsequent 3D seismic acquisitions in two different areas called A (about 0.4 km2) and B (about 0.6 km2), as shown in Fig. 1. In each area, we carried out two surveys: the first one in March 2006 acquiring two cubes called respectively CubeA-spring and CubeB-spring, and the second in July 2006 obtaining the CubeAsummer and CubeB-summer. It is important to underline that the two cubes show differences in the coverage because the shot and receiver displacements are different. In fact, environmental problems strongly influenced the seismic acquisition. In conclusion, two 4D prospections surveyed the areas A and B. Preliminary field tests were performed to define the best sweep lengths and frequency windows of the vibroseis and to analyze the strong ground-roll and the environmental noise. We used the DMT Summit recording instruments with 260 active channels. The MiniVib IVI T-2500 was the seismic source and 10 Hz single geophones were deployed. Both the receiver interval and the receiver lines spacing (in-line) were 18 m. The source interval and the shot line spacing (cross-line) were 18 m and 36, respectively. Considering the receiver and shot intervals and the receiver and the shot line spacing, we adopted a binning equal to 9 × 9 m. The signal sweep of 15 s length and a frequency band of 20–200 Hz were used to provide 2 s length records after cross-correlation of recorded data with the pilot trace with a time sampling interval of 1 ms.
2.1. Processing of seismic data The quality of the raw data was satisfactory, despite of the groundroll (see Fig. 3). The reflections of interest are evident from 0.1 s to about 1 s, while the refractions are not always detectable because of velocity inversions caused by the alternation of impermeable and water filled permeable layers (Bourbié et al., 1987). The processing was focused to enhance the signal-to-noise ratio and to increase the vertical resolution. A ‘preserving-amplitude’ approach was adopted in order to maintain the control of the amplitudes of the signals. This is required for a successful AVO analysis. Firstly, we zeroed the traces or part of them showing a persistent level of noise. Then, we applied a band-pass filter 20–160 Hz to remove the low frequencies mainly related to the ground-roll. After applying a spherical divergence correction using the seismic velocity field obtained from the pre-stack depth migration, a predictive deconvolution was applied with a 120 ms operator length and a lag of 5 ms. A white noise level of 1% was added to stabilize the procedure. In the studied area, the lateral variation of lithological composition in the shallow subsurface layers causes strong shallow lateral velocity variations. Thus, we performed a surface consistent residual static correction using the velocity field obtained from the iterative updating procedure (described in the following sections). The effects of residual static application (window between 160 and 260 ms) is evident in Fig. 4. The obtained velocity fields, converted from interval to RMS velocity, were used to stack the data, after the application of a muting function to remove the direct arrivals. We applied a time-variant filter and F-X deconvolution to the post-stack data to increase the signal-tonoise ratio and the lateral continuity of the signals of the different interfaces. 3. Depth modeling of 3D seismic data The seismic image of an earth model in depth is usually described by two sets of parameters: the seismic velocities and the interfaces
Fig. 3. Example of raw data.
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
397
Fig. 4. Stack section of a 3D seismic line without application of residual statics (on the left) and with application of residual statics (right).
geometry. The estimation of these parameters with a required level of accuracy make the earth model definition a challenging task for the seismic exploration. Subsurface seismic imaging is finally realised by pre-stack depth migration, which, converts the seismic travel-time data into the depth section using the velocity field. The more reliable is the velocity field we supplied to the migration, the more realistic will be the imaging.
The modelling flow-chart adopted in this work is shown in Fig. 5a; it consists in an iterative updating procedure for refining and improving an initial model in depth, involving pre-stack depth migration, residual move-out analysis and seismic reflection tomography. The initial model generally consists in a horizontally layered model with laterally invariant velocities, or in a partially refined model obtained from other techniques and information.
Fig. 5. Flow-chart of a typical modelling in depth procedure (a) and the concept of residual move-out on a Common Image Gather (CIG) (b).
398
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
Residual move-out analysis is a velocity analysis performed after applying an initial velocity function to the data, in order to find the residual errors in the velocity field. If the initial velocity field has been estimated with sufficient accuracy, then the common image gathers (CIGs) derived from pre-stack depth migration using this model should exhibit a flat sequence of events, as shown in Fig. 5b. Any error in layer velocities and/or reflector geometries, on the other hand, should give rise to a residual move-out along those distribution of events on the CIGs, which are no-flat. In other words,
the degree of no-flatness of the reflection events on the CIGs is a measurement of the error in the model, and residual move-out analysis identifies the correction required to flatten the reflection events (Yilmaz, 2001). The adopted inversion technique can be described step by step in the followings points (Yilmaz, 2001). 1. Generate a set of CIGs from pre-stack depth migration using the velocity field derived from the initial model.
Fig. 6. 3D pre-stack depth migration of the CubeA-spring (a) and the CubeA-summer (b) using the final tomographic velocity model (c).
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
2. If the CIGs exhibit flat events, then the model is good and the procedure stops giving the final velocity model and the imaging performed by the pre-stack depth migration. If the events are not flat, the model is not sufficiently accurate yet, and the procedure goes on. 3. Compute the horizon-consistent residual move-out for all offsets along events on CIGs that corresponds to layer boundaries included in the model, obtaining the semblance sections. 4. Pick on the residual move-out profiles for all the horizons by tracking the semblance peaks. 5. Build the travel-time error vector Dt using the picked residual move-out in depth.
399
6. Define the initial model by a set of slowness and depth parameters and construct the coefficient matrix L in the following equation: Δt = LΔp
ð1Þ
where Δp is the tomographic update column vector. 7. Estimate the change in parameters vector Dp, by way of the Generalised Linear Inversion (GLI) solution given by −1 T T L Δt ð2Þ Δp = L L 8. Update the parameter vector p + Dp. By combining the update interval velocity profiles with the new depth horizons, we obtain the updated model.
Fig. 7. 3D pre-stack depth migration of the CubeB-spring (a) and the CubeB-summer (b) using the final tomographic velocity model (c).
400
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
9. Extract the velocity field from the new model and perform the pre-stack depth migration. 10. Return to the point 2 if the result is not satisfactory. At each iteration, both velocity and reflector geometries are updated, until the two set of parameters reach a good degree of stability. In our case, the initial model is the 3D interpolated velocity field obtained from the tomographic inversion of the 2D seismic lines described in Giustiniani et al. (2008; Fig. 2). The initial result was not satisfactory, because the seismic energy is not well focused and the CIGs exhibit non-flat sequence of events. We refined the model and performed the 3D pre-stack depth migration using the commercial package GeoDepth©. During each iteration, we carried out a residual move-out analysis, for all the horizons, by picking the coherences on the semblance profiles. The residual maps and the interval velocity sections extracted from the model for each depth horizon, are the input data for the tomographic update. The result of modelling updating needs to be verified for consistency with the input seismic data and examined for any remaining residual move-outs to decide whether or not to continue with the iterations of tomographic update (Yilmaz, 2001). According with the modelling flow-chart in Fig. 5a, a pre-stack depth migration is performed, using the velocity field derived from the new velocity model. The main consistency checks are the following: 1. Overlay the depth horizons from the updated model onto the image section derived from pre-stack depth migration and note if they
coincide with the reflectors associated with the layer boundaries included in the model. 2. Compute the residual move-out semblance spectra from the image gathers at selected locations after the model update and compare them with those before the update. The residual move-out analysis of image gathers and the update of velocity–depth model should be performed iteratively until the velocity–depth model and the depth image are consistent. Consistency may be achieved with only a few iterations when the residual move-outs are small. Fortunately, this is the case of our seismic survey, because we obtained a good consistency after only three iterations. After the third tomographic update, the events on the CIGs looks nearly flat and the depth horizons associated with the updated model coincide with the reflectors that correspond to the layer boundaries when superimposed on the image section: the model reached a good accuracy. Fig. 6c shows the final 3D velocity model of CubeA. The vertical resolution is satisfactory and we notice that two high velocity layers are evident at about 270 m and 460 m respectively, which could be associated to deeper aquifers. Fig. 6a shows the final migrated CubeA-spring, while Fig. 6b shows the migrated CubeA-summer. As it can be seen, below 250 m depth, two strong reflections appear in both cubes and can be associated to the top of two aquifers (see the high velocity in Fig. 6c). The cubes are sectioned to show the good lateral continuity of the reflectors. In Fig. 6c only one velocity cube is shown because the difference in the
Fig. 8. Correlation of the well stratigraphy with a transversal section of CubeA-summer. The velocity section is overlapped to the seismic migrated section. The dashed violet lines identify the picked horizons on the depth seismic section.
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403 Fig. 9. Stack section of 3D dataset acquired in March 2006 by using the traces with an offset ranging between 9 m and 149 m (panel a), between 150 m and 299 m (panel b) and greater than 300 m (panel c). The xline 45 is shown.
401
402
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
P-wave velocity values between the CubeA-spring and CubeAsummer velocity models is less that the estimated velocity error. Assuming for the first two acquifers an average picking error of about ±6 ms, the resulting average velocity deviation inside the two aquifers is about 2%. For the deeper acquifers, the velocity error increase and could be estimated to be about 5% for this inversion algorithm (Tinivella et al., 2002). In Fig. 7 the final velocity CubeB and the resulting migrated cubes are shown. We notice that the resolution and lateral continuity is better than the CubeA because of the higher coverage. Note also that the two reflections identified in the two CubeA at about 270 m and 480 m are here better visualized. Moreover, other two layers were identified and picked between 30 and 180 m. In the chaircut of the migrated CubeA-spring (Fig. 7a) a lateral amplitude inversion in the shallow layers is evidenced (with a circle). In the migrated CubeAsummer (Fig. 7b) it is not evident, but the coverage in this zone is lower with respect to the CubeA-spring. This lateral amplitude inversion may be due to lateral changes in the lithology or to a fault. 4. Correlation between seismic and well data On the basis of the seismic data results, a 510 m deep well was drilled, which allowed us to check the real consistency of our 3D models. Fig. 8 shows a transversal depth section cutting the CubeAsummer, overlapped to the corresponding interval velocity section, together with the well stratigraphy and the down-hole velocity log. There is an excellent correspondence between the horizons picked on the migrated section and the main discontinuities reported in the well stratigraphy, except for some slight discrepancies, which could be due to the distance between the well and the survey area (about 200 m, see Fig. 1). The two high velocity layers at about 270 m and 480 m correspond in the well stratigraphy to two deep aquifers. Those aquifers, because the catching wells were limited to depths of about 200 m, were unknown. Moreover, Fig. 8 shows the good correlation between the down-hole velocity log and the interval velocity section for aquifers 2 and 3. The four aquifers indicated on the well stratigraphy correspond to high interval velocity values, which is mainly related to the lithological variation. In fact, the velocity is affected mainly by three factors: 1) the matrix, 2) the pore fluid and 3) the pore pressure. In our case, as shown in Fig. 8, aquifer layers are composed by sand and gravel, while the adjacent layers are mainly composed by water saturated sandy clays. Note that, because of compaction effects, clay velocity increases with depth. It is well known that, for unconsolidated sediments, sandstones are characterised by higher velocities than those of shales (Schön, 1996). Moreover, the pore fluid is in overpressure conditions, slightly reducing the seismic velocity (Giustiniani et al., 2008). So, the combined effect is the increase of the seismic velocity in correspondence of the aquifers. 5. Amplitude versus offset The 3D pre-stack data analysis underlines the presence of AVO effects associated either with lithological changes and the presence of fluids at high pressure. These effects are evident, for example, in the three in-line stacked sections of the 3D dataset obtained by using different offset ranges: 0–149 m, 150–299 m, and higher than 300 m (Fig. 9). This figure shows that the amplitudes at the two main shallow aquifers change with offset. In the far-offset section, the amplitude is stronger compared to the others: this is evident for the deeper reflectors, at about 270 m and 480 m. As expected, the near-offset stacked section provides more information about the shallowest targets. To better discriminate between lithological and fluid changes at each interface, we performed the AVO analysis by using commercial software (Vern and Hilterman, 1995). We linearized Zoeppritz's
equations using the Aki–Richards method (Aki and Richards, 1980). The adopted approach is described in details in Giustiniani et al. (2008). The AVO inversion was applied to the 3D seismic dataset. The reliability of the result depends on the seismic data quality and on the coverage. As well as travel-times tomography (see previous section), also AVO analysis is not able to detect any seismic response difference between the two seasons. As already observed in the analysis of 2D seismic data (Giustiniani et al., 2008), the comparison between P- and S-wave reflectivities confirms that the top of each aquifer is well evident, confirming that the seismic interface is mainly generated by lithological changes and only slightly by the fluid condition changes. 6. Discussions and conclusions In this work we studied some strategies in order to characterise an important multilayered aquifer in an area of the NE of Italy by using the seismic method, namely seismic tomography, imaging and AVO. The 3D seismic acquisition was based on the results obtained from a previous 2D seismic survey and performed in two different seasons (March and July, 2006). We successfully applied a modelling in depth procedure, obtaining detailed interval velocity and migrated cubes for both seasons. We also applied the AVO analysis in order to obtain petro-physical properties of the multilayered aquifer. In this study, we obtained two significant results. The first is that, even if the imaging in depth of the aquifer system is satisfactory, travel-time tomography cannot detect any difference in the seismic response between the two seasons, because the velocity differences are less that the estimated velocity error (2% for the shallow aquifers and 5% for the deeper). This is also confirmed by the results obtained from the AVO analysis, which cannot detect overpressure variations between the two seasons. In aquifer studies, the seismic analysis (both tomographic and AVO inversion) detects mainly lithological variations, because the fluid pressure conditions in the aquifer slightly influences the seismic properties compared to the lithological effects. In fact, the stratigraphies of several wells located in the study area evidence that the aquifers are host in gravel layers confined between clay layers. Therefore, the 4D approach does not give other information with respect to the 3D seismic method, because travel-time tomography and AVO analysis are not enough sensitive to detect the overpressure variations inside the aquifers. The second result is that, because of the strong lithological variation in the study area, we obtained a detailed geometry and interval velocity model that allowed the discovery of two unknown aquifers, located at about 270 m and 480 m. In fact, a new 510 m deep well was drilled at about 200 m from the investigated area A (Fig. 1), confirming the existence of two aquifers identified by the seismic data analysis. Concluding, the 3D seismic method can be successfully applied in aquifer exploration. 3D cubes allow extending the well information, mainly lithological, to the surrounding area. This information can be used to optimise the location of new wells, or to detect deeper sources. Acknowledgments We are very grateful to Elvio Del Negro for his contribution and technical support in collecting the field data. We acknowledge the European community that has supported the project “Water-bearing characterisation with integrated methodologies – CAMI” (LIFE Project Number – LIFE04 ENV/IT/000500) and the leader of the project Daniel Gustavo Nieto Yabar. We wish to thank the Acquedotto Basso Livenza, and in particular Enrico Marin, for the logistical support and the all information provided (hydrogeological data and the stratigraphies of the catchment wells). We are grateful to the reviewers (Rinaldo Nicolich and Adele Manzella) for useful suggestions.
M. Giustiniani et al. / Journal of Applied Geophysics 68 (2009) 394–403
References Accaino, F., Tinivella, U., Rossi, G., Nicolich, R., 2005. Geofluid evidence from analysis of deep crustal seismic data (Southern Tuscany, Italy). Journal of Volcanology and Geothermal Research 148, 46–59. Aki, K., Richards, P.G., 1980. Quantitative Seismology: Theory and Methods. W.H. Freeman, San Francisco. Bachrach, R., Nur, A., 1998. High-resolution shallow–seismic experiments in sand, part I: water table, fluid flow, and saturation. Geophysics 63, 1225–1233. Bradford, J.H., 1998. Characterizing shallow aquifers with seismic reflection, part II — prestack depth migration and field examples. Geophysics 67, 98–109. Bradford, J.H., Sawyer, D.S., 2002. Characterizing shallow aquifers with seismic reflection, part II — prestack depth migration and field examples. Geophysics 67, 98–109. Bourbié, T., Coussy, O., Zinzner, B., 1987. Acoustics of Porous Media. Gulf Publishing, Houston, Texas. Cardimona, S.J., Clement, W.P., Kadinsky-Cade, K., 1998. Seismic reflection and groundpenetrating radar imaging of a shallow aquifer. Geophysics 63, 1310–1317. Giustiniani, M., Accaino, F., Picotti, S., Tinivella, U., 2008. Characterisation of shallow aquifers by high-resolution seismic data. Geophysical Prospecting 56, 655–666.
403
Nicolich, R., Della Vedova, B., Giustiniani, M., Fantoni, R., 2004. Carta del sottosuolo della Pianura Friulana. Litografia Cartografica, Firenze. Steeples, D.W., Miller, R.D., 1990. Seismic reflection methods applied to engineering and groundwater problems. In: Ward, S.W. (Ed.), Geothechnical and Environmental Geophysics 1, pp. 1–30. Schön, S., 1996. Physical properties of rocks—fundamentals and principles of petrophysics. Handbook of Geophysical Exploration: Seismic Exploration, vol. 18. Pergamon, Oxford/Tarrytown. Tinivella, U., Accaino, F., Camerlenghi, A., 2002. Gas hydrate and free gas distribution from inversion of seismic data on the South Shetland margin (Antarctica). Marine Geophysical Research 23, 109–123. Vern, R., Hilterman, F., 1995. Lithology color-coded seismic sections: the calibration of AVO crossplotting of rock properties. The Leading Edge 14, 847. Whiteley, R.J., Hunter, J.A., Pullan, S.E., Nutalaya, P., 1998. Optimum offset” seismic reflection mapping of shallow aquifers near Bangkok, Thailand. Geophysics 63, 1385–1394. Yilmaz, O., 2001. Seismic Data Analysis: Processing, Inversion and Interpretation of Seismic Data. SEG Series: Investigation in Geophysics, Tulsa.