© J.-L. Mari et al., published by IFP Energies nouvelles, 2021

1 Introduction

The need for images of the shallow subsurface at high resolution is claimed by various applications including geotechnics, hydrology, reservoir engineering, etc. Geophysics has the capability to render such images. For example, seismic data obtained from shot points with a single geophone per trace and small distances between geophones (1–5 m) allows for the simultaneous recording of different waves (refracted wave, surface wave, reflected wave) which can then be processed independently or in combination.

Some specific processing methods have been developed such as: travel time tomography to obtain P-wave velocity models [1, 2], surface wave analysis to obtain shear velocity models [3, 4], deterministic and statistical seismic inversion [5], and Full Waveform Inversion (FWI) for P-wave velocity estimation of near-surface seismic data [6]. Surface seismic methods, such as refraction tomography and reflected wave processing, can be combined to extend a P-wave velocity model from the surface to hundred meters over depth [6, 7]. FWI is an advanced seismic imaging method [8] with the objective to determine a model (velocity, density, and possibly anisotropy and attenuation) of the subsurface in which the synthetic shots gather the best fit of observed data [9]. Since the early 2000s, FWI has been developed from 2D to 3D, from acoustic to elastic and visco-elastic, in anisotropic contexts, from offshore to onshore datasets. FWI remains under development, especially for multi-parameter estimation (i.e., not only pressure velocity models from body waves, but also shear velocity models from surface waves, and anisotropy or attenuation parameters). The difficulty is to extract more than one parameter [10].

Through an integrated field case obtained on an experimental site, this work shows that a solution exists to improve our knowledge regarding the structure and the properties of the subsurface. This solution is clearly associated with investigations crossing and inter-comparing the results from diverse geological and geophysical techniques: 2D–3D Very High Resolution (VHR) seismic imaging including refraction tomography, acoustic logging, and flow measurements.

An experimental site (Fig. 1), located in the Cher region (central part of France), has been developed for the training of both students from universities and professional practitioners. The site belongs to Géocentre Forsol (https://www.geocentre-forsol.fr), a company involved in geotechnical studies and drilling. It is also used for experimental studies in near subsurface geophysics. The training in geophysics covers the acquisition and processing of surface-monitored seismic data in 2D or 3D. In addition, two boreholes are available (Fig. 2). These boreholes allow for the acquisition of well seismic data such as Vertical Seismic Profiles (VSP) and loggings such as full wave form acoustic data and flow measurements.

After a short description of the geological context and a review of the seismic imaging carried out at the site, this work shows how acoustic logging and refraction tomography can be merged to form a high-resolution continuous velocity model from the surface up to the terminal depth of the boreholes.

The benefit of combining hybrid seismic methods (reflection seismic processing, refraction tomography, VSP processing), and acoustic logging is also emphasized with the idea to extend laterally the velocity model obtained at a borehole.

A seismic porosity distribution is identified, which informs on areas where preferential flow would occur. Because the resolution of the pseudo-porosity sections is not sufficient to detect the presence of fractures, criss-cross events revealed by acoustic sections can be used to build an index log for detecting fractures.

2 Geological context

The site extends straight at the transition between Triassic and Jurassic geological layers, south to the Paris basin (Fig. 1a). A simplified geological cross-section oriented North–West South–East has been established (Fig. 1b).

Two boreholes are available at the site (Fig. 2), one (B1) hardly investigable as it is steel cased over its whole depth (0–90 m), the other (B2) being more accessible to measurements, as it is only steel cased in its upper part (0–78 m), and then slotted PVC cased in its lower part (78–200 m). Borehole B2 has been cored in the 78–200 m depth interval. A core description has been done by Anthony Hardy (Géocentre Forsol):

• 0–78 m: no core, drilling in a destructive mode,

• 78–104 m: dark grey/brown dolomitic limestone,

• 104–112.2 m: bioclastic grey/beige limestone,

• 112.2–122.4 m: dark grey to light grey bioclastic limestone,

• 122.4–138 m: alternating of flaky greenish argillite and dolomitic gregarious limestone banks,

• 138–144 m: grey/orange coarse sandstone with greenish pasts,

• 144–152 m: slightly sandy reddish/greenish argillite with sandy-gravel pasts,

• 152–158 m: grey/orange coarse sandstone with greenish pasts,

• 158–160.3 m: reddish argillite,

• 160.3–162 m: grey/orange coarse sandstone with greenish pasts,

• 162–164.5 m: alternation argillite/greyish coarse sandstone,

• 164.5–177 m: grey coarse sandstone with reddish and purple pasts,

• 177–179.5 m: reddish undid argillite,

• 179.5–183 m: dolomite and red/greenish argillite,

• 183–187 m: coarse grey sandstone (little recovery),

• 187–191 m: dolomite and red/greenish argillite,

• 191–200 m: coarse sandstone and then medium-grain grey to past wine/purple binds.

The core analysis also highlighted fractured zones mainly between 70 and 87 m depth, between 120 and 130 m, and between 150 and 155 m. Fractures are clearly visible at: 122, 126 and 152 m (Fig. 1c, unpublished results from A. Hardy, Géocentre Forsol).

A simplified description of the geological series can be summarized as follows: Recent superficial deposits overlay a sedimentary formation with a thickness of approximately 200 m. The sedimentary formation is mainly composed of Liassic limestone down to 122 m depth, the alternation of argillite and dolomitic, sandy limestone banks between 122 and 138 m (probably the intermediary Rhaetic facies), and Triassic sandstones with some argillite and dolomite intercalations between 138 and 200 m depth. The Permian substratum could be reached at 230–250 m depth.

3 Seismic imaging

Usually, 3D seismic imaging requires a great amount of data and powerful memory and CPU computer. A 2D profile has been recorded. A 3D survey has been designed to obtain a 3D cube by recording a low amount of data [11]. The example in [12] reports on the optimum choice of the seismic spread design and on the diverse signal processing steps allowing for transforming the recorded seismic data into a 2D vertical section in multiple folds (up to 48) and 3D cube in multiple folds (up to 22).

The field case presents a refraction-reflection imaging strategy with the capability to evaluate reflectivity information from the acquisition at the ground surface [7]. The procedure developed to obtain a composite section over depth was already applied successfully on these 2D–3D datasets [12].

3.1 2D seismic imaging

For the 2D seismic profile (Fig. 2a), the receiver spread is fixed. It is composed of 48 geophones, 2 m apart. The source is a lightweight dropper (Fig. 2a) which is moved and fired between 2 adjacent geophones. Consequently, the distance between 2 adjacent Common Mid Points (CMP) is 1 m. The listening time is limited to 250 ms, the sampling time interval is 0.5 ms.

The composite sections obtained by hybrid seismic methods are processed via a procedure in four steps [13]:

• Construction of a velocity-over-depth model from first arrival times, constructed iteratively by tomographic inversion. The Simultaneous Iterative Reconstructive Technique (SIRT) is an Analysis Step and hybrid algorithm used to invert the first arrival travel times (refracted waves) from seismic reflection data for building a near-surface velocity model [1]. This process works for both 2D and 3D datasets. The depth velocity model obtained by refraction tomography is used to build a depth-to-time conversion model (variable in distance) and consequently a velocity model in two-way times. The velocity model over time is converted in a reflectivity section up to 20 ms. The results of the tomographic inversion are shown in Figure 3a.

  Fig. 3 2D Seismic imaging (hybrid seismic processing). (a) Velocity model obtained by tomography in depth and time. (b) 2D seismic section in time. (c) 2D composite seismic section in time. (d) 2D composite seismic section in depth.

• Construction of a time reflectivity section by classical reflection seismic processing. The processing is carried out with the SPW (Seismic Processing Workshop) software developed by Parallel Geoscience [14]. The processing sequence of each shot includes amplitude recovery, deconvolution in the 15–150 Hz frequency bandwidth, tail mute, static corrections. A deconvolution is performed to increase the resolution and attenuate the surface waves. A tail mute is used to cancel out the air waves and the surface waves. Static corrections are made to compensate the effects of the weathering zone. In the reported example, the static corrections are very weak. A stacking velocity model performed by velocity analysis is used to obtain stacked sections. Surface consistent residual statics are computed to enhance the signal to noise ratio and preserve the high resolution of the data in the stack procedure. The residual statics improve the quality of the stacked sections. To increase the lateral coherency of the seismic section, a 3-terms moving average filter has been applied. The 2D section is presented in Figure 3b. Theoretically, the 2D section is composed of 96 CMP, 1 m apart. Practically, the tail mute used to cancel out the surface waves reduces the length of the profile to the 88 central CMP, which gives a trapezoidal shape to the seismic section.

• Extension up to the surface of the time reflectivity section. The operation is carried out by converting the shallowest depth velocity model to time reflectivity, associated with velocity contrasts in the subsurface. The time reflectivity sections require the application of a scaling factor before being merged in a final time reflectivity section [13]. The 2D section after extension up to the surface is shown in Figure 3c.

• Depth conversion. A Vertical Seismic Profile (VSP) [15, 16] was recorded in borehole B1. The measurement of the arrival time of the first down-going waves propagating at close-to-normal incidence is used to provide both a time versus depth law and a velocity distribution in the subsurface. After processing, the VSP provides a seismic trace (VSP corridor stack trace) directly comparable to the surface seismic section passing in the vicinity of the borehole. VSP is used to predict events (reflectors or heterogeneities) located beneath the borehole. The law of VSP time versus depth, t = f(z), can be used to convert the time seismic sections (Fig. 3c) into seismic sections in depth scale (Fig. 3d).

Figure 3 shows the four steps of the procedure used to obtain the 2D composite section.

3.2 VSP and 3D seismic imaging

Figure 4a shows the VSP recorded in borehole B1, with the horizontal axis representing the different depths of the borehole geophone (Fig. 4b), and the vertical axis being the listening time. In this example, the depth of the sensor varies between 25 and 90 m, the surface source (Fig. 2a) is slightly offset (5 m) from the wellhead. The distance between two successive positions of the borehole geophone in the borehole is 5 m. The sampling time interval is 0.25 ms for a listening time of 250 ms. In Figure 4, the representation of the listening time was limited to 100 ms.

After editing and calibrating with a reference geophone located at the surface, the first step of a VSP processing sequence is that of picking the first arrival times of the down-going wave at each depth position. The calibration compensates for the variations of the source energy and of the Time Breaks (TB). The second processing step is the separation of body waves (down-going and up-going waves) in the wave field. As the number of VSP traces is limited to 13, the wave separation is done via a SVD filtering (singular value decomposition) after amplitude compensation. The wave separation is performed according to the following steps [13]:

1. Time shift of the VSP section to flatten the down-going wave;

2. Extraction of the down-going wave by SVD filter (first eigensection);

3. Computation of a residual section which is the difference between the VSP section and the estimated down-going wave;

4. Time shift of the residual section to flatten the up-going wave;

5. Extraction of the up-going wave by SVD filter (first eigensection).

The extracted down-going and up-going waves are then relocated at their initial time positions.

The third processing step applied to data is the deconvolution and the computation of a corridor stack trace. The deconvolution of up-going waves by the down-going waves is performed trace-by-trace to increase the vertical resolution of the VSP and to obtain a zero-phase VSP section. A corridor is designed to select and to stack a part of the deconvoluted and flattened up-going waves. The result is a corridor stack trace directly comparable to a seismic trace extracted from the processed surface seismic data. The VSP corridor stack trace shows seismic reflections up to 200 Hz. Figure 5 shows the different steps of the VSP processing.

The 3D seismic spread (Fig. 2b) is composed of a receiver spread and a source spread. The receiver spread, displayed in blue in Figure 2b, is composed of two receiver lines, indicated by the geophones numbered 1–24 and 25–48. The receiver line direction is called the in-line direction. The distance between the receiver lines is 4 m, each line counts 24 geophones with a lag distance of 2 m between two neighbors. The source spread, displayed in yellow in Figure 2b, is composed of 11 source lines, noted T_n,m (n: line number and m: shot number) oriented perpendicularly to the receiver lines, with 11 shots fired per line and a distance between shots of 2 m. The distance between source lines is 4 m. The distance between the receiver spread and the source spread is 4 m, with no overlap between spreads. The CMP domain (white rectangle in Fig. 2b) is composed of 11 CMP lines (1 m apart) in the receiver line direction with 44 CMP (1 m apart) per line [12]. The azimuth of the 3D spread is given by a dotted line displayed in white in Figure 2. Notably, this line was also that used for the implementation of the 2D seismic profile. The 3D recording parameters, i.e. the listening time and the sampling time interval are similar to the ones employed in the 2D survey.

The processing sequence, defined for the 2D profile, has been applied to the 3D data set to obtain a composite seismic block extended up to the surface. The tail mute used to cancel out the surface waves gives a trapezoidal shape to each 3D in-line seismic section. Figure 6 shows an example of an in-line (line #3) and a cross-line (line #15) composite seismic sections extracted from the 3D block, both in time (Fig. 6a) and in depth (Fig. 6b).

The CMP point #27 of the in-line seismic section #3 is located approximately 40 m away from the borehole B1 in which the VSP has been recorded. The VSP stacked trace, duplicated five times, is inserted in the in-line section #3 at the CMP location #27 (Fig. 7). The results are shown both in time (Fig. 7a) and in depth (Fig. 7b).

The correlation coefficient between the seismic trace and the VSP stacked trace at this CMP point is greater than 0.75, thus showing a good fit of the depth evidenced by the seismic horizons. It is worth noting that the VSP renders a calibration of the surface seismic in the range of depths 25–90 m where the measurements were done, but also beneath the borehole.

The VSP stacked trace highlights reflectors that appear at depths greater than the borehole depth. This example illustrates the predictive role of the VSP, as an identification of features beneath the borehole. On all sections (Figs. 6 and 7), we observe strong reflectors at 40 m, 80–90 m, 120 m, and 160 m as predicted by the VSP. We can also observe a strong reflected event at 245 m (Fig. 7b), which could correspond to the top of the Permian.

The 2D composite section (Figs. 3c and 3d) has a better resolution and a better signal to noise ratio than a 3D composite section (Figs. 7a and 7b), these differences being the consequence of the diverse folds employed and differences in the distribution of offsets:

• Fold: The maximum fold in the central part is 48 for the 2D line and 22 for the 3D block. Consequently, we have a better signal to noise ratio in the 2D composite section than in a 3D composite section.

• Range of offsets: The offset distribution is more regular in 2D than it is in 3D. There is a lack of small offsets in 3D to highlight the reflectors between the surface down to 40 m depth.

4 Acoustic logging

As indicated earlier, two boreholes are available at the site. They are marked by green and red crosses in Figure 2. Borehole B1 was drilled by Geocentre in 2006 but is now fully steel cased and cemented. Borehole B2 was drilled by Géocentre Forsol in two phases between September 2019 and September 2020. During the drilling phases, some parameters such as the Rate Of Penetration (ROP) and torque were continuously recorded. The first drilling phase from the surface up to 78 m depth resulted in a steel cased but not cemented borehole. Re-handling the borehole within the second drilling phase allowed to reach the depth of 200 m with a borehole completely cored between 78 and 200 m, then equipped with a slotted PVC casing.

To evaluate the borehole conditions, a full wave form acoustic tool was run into the two boreholes. The acoustic tool [17] is a monopole-type flexible tool with a small diameter of 50 mm. The transmitter is magnetostrictive (transmission frequencies: 17–22 kHz). It can be equipped with two pairs of receivers (near receivers at 1 and 1.25 m beneath the source, and far receivers 3 and 3.25 m beneath the source).

In a vertical well, monopole tools enable the recording of different wave fields [18]:

• Body waves such as refracted waves. Arrival times of refracted P-waves are used to measure the P-wave velocity of the formation.

• Dispersive borehole guided waves, such as Stoneley waves.

In addition, for a constant offset (distance between emitter and receiver) in acoustic sections, coherent slanted events, also referred to as criss-cross events, can be observed. They correspond to refracted waves reflected at the edges of geological discontinuities (acoustic impedance discontinuities), such as fractures. The amplitudes of criss-cross events are weak in comparison with the amplitudes of the refracted waves.

The acoustic tool is used here with its far offsets’ configuration. The sampling interval over depth is 2 cm for borehole B1 and 4 cm for borehole B2. The sampling interval over time is 5 μs and the listening time is 5 ms. Figure 8a shows the 3 m-offset acoustic sections (R1 section) recorded in boreholes B1 and B2.

The acoustic section recorded in borehole B2 shows strong resonances in the time interval 0.6–3 ms, between 30 m and 78 m. The resonances with extraordinarily strong amplitude free of any attenuation clearly indicate that there is no coupling between the steel casing and the geological formation. B2 is an uncemented cased hole up to 78 m depth. Consequently, this part of the borehole cannot be used to conduct any type of logging for geological studies.

Since boreholes B1 and B2 are located close to one another, acoustic data recorded in borehole B1 along the 30–78 m depth interval and in borehole B2 along the 78–192 m depth interval were merged together to obtain a continuous acoustic information over the 30–192 m depth interval with a depth sampling rate of 4 cm. Figure 8b shows the composite acoustic section obtained by merging acoustic data recorded both in borehole B1 (steel cased hole) in the 30–78 m depth interval and in borehole B2 (slotted PVC cased hole) in the 78–192 m depth interval. Composite acoustic sections were processed to obtain acoustic logs. Figure 8b shows the Cemented Bond Log (CBL) which can be correlated with the resonances observed on the acoustic section. The CBL can also be used as a noise level indicator in the 40–70 m and 100–192 m depth intervals. Figure 8b also shows the P-wave velocity log, obtained by picking the arrival times of the refracted P-wave observed at the 2 offset sections (3 m–3.25 m). The picked times are used to flatten out (at time 0) the refracted P-wave on both sections (sampled at 3 and 3.25 m beneath the source). Then, for a prescribed position z of the probe, the amplitudes of the signal at 3 m and 3.25 m are sampled within a short time window to establish the correlation between observations at 3 and 3.25 m. A high value of correlation coefficients indicates a good velocity measurement. More than 86% of P-wave velocity values have a correlation coefficient greater than 0.7, and more than 77% of P-wave velocity values have a correlation coefficient greater than 0.8.

4.1 Calibration of seismic section by full waveform acoustic logging and VSP

The velocity log obtained by a Full Waveform Acoustic Logging (FWAL) is a measurement carried out continuously over the logged depth interval, sampled at regular space steps (4 cm in the present case). It can be used to identify the velocity contrasts associated with seismic horizons observed on the seismic section. As an example, the seismic horizon at 120 m depth, predicted by the VSP stack trace, corresponds to a strong decrease of velocity as indicated by the acoustic velocity log.

Figure 9 shows the comparison between the 3D composite seismic section, the acoustic velocity log (inserted in the seismic section at CMP 15#, as a projection of borehole B2), and the VSP stacked trace. The figure highlights the difference of vertical resolution between full waveform acoustic logging (several tens of centimeters) and seismic imaging (several meters).

4.2 Fracture detection

The full-bore Formation Micro-Imager (FMI) provides images of micro-resistivity or formation images. It is currently used to reveal sedimentary forms, determining dip in water-based mud, and understanding structural and stratigraphic features, including fractures [19, 20].

Here we show how an FWAL log can be used to detect fractures from the criss-cross events. The composite acoustic section (Fig. 8b) was processed to separate, on the one hand, refracted waves and Stoneley waves and, on the other hand, criss-cross events (Fig. 10). After separation, 2 sets of criss-cross events are evidenced (Fig. 10a): slanted events with a positive apparent velocity, and slanted events with negative apparent velocity. The slanted events with a positive apparent velocity are shifted over depth to become vertical events at various depths. The same procedure is applied to slanted events with negative apparent velocity. Then, the amplitudes of vertical events are stacked in a short time window (1-ms duration) to obtain an energy log as a function of depth. The energy log is normalized by the highest energy value to obtain a log ranging between 0 and 1.

The normalized energy log computed from the distribution of criss-cross events observed on constant-offset section is used to estimate the criss-cross density index over depth, via a criss-cross index log (Icriss). Criss-cross events are refracted-reflected waves at the limits of acoustic impedance contrasts. Figure 10a exemplifies these criss-cross events and the associated criss-cross index logs. The automatic detection of criss-cross events allows for detecting the levels associated with strong acoustic impedance contrasts, some of them pointing out the presence of fractures. Low amplitude criss-cross events can be seen in the depth intervals of 30–65 m, 95–110 m, 170–192 m. The peak of criss-cross index between 92 and 95 m is an artifact due to the presence of a piece of casing. Many high amplitude criss-cross events are clearly visible in the following depth intervals: 65–90 m, 110–135 m, 145–160 m. A high density over depth of criss-cross events can often be related with the presence of fractures. The fracture zones are also revealed by a decrease in the P-wave velocity. A fracture can be detected by both a high value of criss-cross index (Icriss(z)) and a low value of P-wave velocity (VP(z)). A fracturation index is defined as follows:

$$\mathrm{IFrac}\left(z\right)=\mathrm{ICriss}\left(z\right)\left(1-\frac{\mathrm{VP}\left(z\right)}{{\mathrm{VP}}_{\max }}\right).$$(1)

The fracturation index highlights fractured zones mainly between 70 and 87 m, between 120 and 130 m, and between 150 and 155 m (Fig. 10b). Fractures are clearly visible at: 122, 126 and 152 m. The analysis of cores (Fig. 11) confirms the results given by the fracturation index log.

To confirm the interpretation given to the acoustic logs, additional measurements were carried out in borehole B2 (slotted PVC cased hole) by means of flow logs.

5 Flow-log measurements

Flow-log measurements were performed with a downward moving probe at constant speed corresponding to a constant rotation of the propeller of the flowmeter, the rotation speed being on the order of 3.3 rotations/second. In addition to the flow logs, a gamma ray (GN) was recorded. It is of weak amplitude with locally some peaks, associated with low values of velocity, which indicate the presence of thin shaly layers. The velocity log is used to compute a porosity log, by using the empirical equation of the sonic porosity ∅_RH established by Raymer et al. [21], adapted to carbonate formation [22]:

$${\varnothing }_{\mathrm{RH}}=C\times \left(\frac{∆t -{∆t}_{\mathrm{ma}}}{∆t}\right),$$(2)where ∆t is the measured sonic “travel time” (inverse of a velocity), ∆t_ma is the theoretical sonic travel time for the rock matrix, and C is a calibration coefficient. The ∆t_ma and C values are usually calibrated on core data via laboratory measurements. In our case, as we do not have laboratory measurements, we used C = 0.72 and ∆t_ma = 212.1 μs/m as the values given in [22].

Figure 12 shows from left to right: P-wave velocity, Gamma ray, porosity, temperature, conductivity, and flow. The temperature trend (red curve) indicates a gradient of 0.026 °C/m.

We remind that:

• In a full cased well, the propeller rotation value corresponds to the descent speed of the logging tool (here, 3.3 rotations/second).

• A constant rotation corresponds to zero vertical water speed in the borehole.

• A decrease in rotation corresponds to a downward flow of water.

• An increase in rotation corresponds to an upward flow of water in the borehole.

• An upward or downward flow of water in a borehole is related to the hydraulic head between two levels.

• Because of the very high sensitivity of the propeller, only significant variations of propeller rotation can be interpreted in terms of flow. The weak variations can be associated with a change in the borehole diameter or even to the slight variations of diameter at the junction between elements of the slotted PVC casings.

The interpretation of the flow (rotation speed) and temperature profiles (Fig. 12) leads to the following results.

• Between 30 and 83 m depth, over a full casing, the constant rotation corresponds only to the downward speed of the tool, therefore null water flow in the borehole.

• At 83 m, a sudden decrease in rotation is observed.

• Between 83 and 143 m, the rotation is low (1.6 rotations/second). The decrease in rotation indicates a downward flow between 83 (arrival) and 143 (exit) m depth. The observation is corroborated by an almost constant temperature (14.2 °C).

• The disturbance observed at 132 m depth could be interpreted, at least, as a decrease of the downward velocity of flow into the well. This feature could be the consequence of a leak from the well toward the formation, but the slight increase (hardly visible in Fig. 12), at the same depth, of water temperature would also indicate an inflow from the formation to the well. Since a partial inflow in the well (for a constant prescribed flux moving from 83 to 143 m depth in the well) would also very locally diminish the overall downward velocity in the well, the 132 m depth horizon is probably a very local inflow from formation to well.

• The disturbance observed at approximately 137 m is probably linked to a fracture but is not visible on the temperature profile.

• At 143 m depth, a very sudden increase in temperature occurs.

• Between 143 and 152 m, the rotation of the propeller is 3.3 rotations/second, which corresponds to the descent speed of the tool. Therefore, between 143 and 152 m depth, no flow occurs in the borehole.

• Between 152 and 156 m, the rotation is low (on the order of 1.8 rotations/second). It is the same type of behavior as between 83 and 143 m, with here a very local entry at 152 and an exit at 156 m depth. The constant temperature between 153 and 156 m would also tend to confirm this hypothesis.

• Between 159 and 181 m, the rotation is low (on the order of 1.8 rotations/second), witnessing a downward flow between 159 m (entry) and 181 m (exit). This observation is corroborated by a constant temperature (15.8 °C) between these depths.

• Beyond 183 m, the rotation corresponds to the lowering speed of the tool.

The flows, observed at different depth intervals (83–143 m, 152–156 m and 159–181 m), are characterized by both a sharp decrease in the rotation of the propeller of the flowmeter and a constant temperature of the water between each entry/exit. Those are associated with fractures detected by acoustic logging at 84 m, between 85 and 87 m, between 120 and 130 m, and at 152 m. Flows are evidenced by low velocity levels corresponding to high porosity levels. We can also notice that the flow areas are completely independent (not related) to each other.

6 Discussion

After calibration, and under the assumption that the slower the wave velocity, the higher the permeability – porosity, the lines of the seismic block inform on preferential areas where flow should occur. At the scale of an open wellbore, acoustic loggings that measure wave velocities over a short distance within the well, also inform on open features crosscut by the well.

Seismic imaging of the site is of limited extension over space. For the 3D, the extension of the block is 44 m in the in-line direction and 12 m in the cross-line direction. For the 2D, the extension of the line is 88 m, but limited to 40 m (CMP 30 to CMP 70) for a depth investigation from the surface up to a depth of 200 m (Fig. 3d). The velocity log obtained from the composite acoustic section has been extrapolated up to the surface by relying upon both the velocity distribution given by the surface tomography (Fig. 3a) and the torque to ROP ratio [23]. The extended acoustic velocity log is shown in Figure 13a. Each CMP stacked trace being assumed as an estimate of the reflectivity function R(z) at the CMP location, the reflectivity can be transformed into a pseudo-velocity function V(z) using the following equation:

$$V\left(z+\Delta z\right)=V\left(z\right)\frac{1+R\left(z\right)}{1-R\left(z\right)}.$$(3)

In this way, the 2D depth section has been transformed into a 2D pseudo-velocity section. A calibration operator has been designed to fit the pseudo-velocity trace of the 2D section, located at CMP 50, with the extended acoustic velocity log, and resampled at the seismic sampling rate. To extend the 2D velocity distribution over depth, the operator has been applied to the 2D pseudo-velocity section, between CMP 30 and CMP 70. Figure 13a shows the seismic velocity distribution compared with the extended velocity log. The velocity distribution has been converted in porosity using the Raymer equation (2). The results are shown within the 30–190 m depth interval in Figure 13b. The high-porosity layers appear in red.

The flowmeter measurements evidence that water flows from the formation into the well at 83–85 m depth, then flows downward into the well and goes back into the formation at 143 m. The same geometry of flow occurs between 152 and 156 m, and then between 159 m and 181–183 m. The water “loops” are indicated by blue arrows on the flow log, the arrows also pointing out that the Stoneley waves (between 2 and 3 ms on the FWAL, in Fig. 13) are strongly attenuated at the depths of water loops.

This field investigation exemplifies how both seismic and acoustic data might become well-suited tools to image porous layers within the subsurface. That being said, the spatial resolution is not sufficient to detect the presence of fractures (revealed by acoustic logging) on the pseudo-porosity sections.

As another example [24] justifying the worthiness of the method discussed above, a summary of an extensive study previously undertaken at a larger scale and over a karstic limestone aquifer, is presented hereafter.

7 Example of a near surface karstic reservoir

This study consisted of a 3D seismic survey of the near-surface karstic reservoir located at the Hydrogeological Experimental Site (HES) of the University of Poitiers (France). The detailed information on the study is available in [24, 25].

The processing of the 3D data leads to a 3D velocity block. The velocity block was then converted into pseudo-porosity values. The resulting 3D seismic pseudo-porosity block reveals three high-porosity, presumably-water-productive, discontinuous layers, at depths of 30–40, 85–87 and 110–115 m (Fig. 14).

The HES field case shows how Full Wave Acoustic Logging (FWAL) can be used to validate the results obtained from the 3D seismic if the karstic body has a lateral extension over several seismic cells. If karstic bodies have a small extension, FWAL in open hole can be fruitfully used to:

• Detect highly permeable bodies, thanks to measurements of acoustic energy and attenuation;

• Detect the presence of karstic bodies characterized by a very strong attenuation of the different wave trains and a loss of continuity of acoustic sections;

• Confirm the results obtained by Vertical Seismic Profile (VSP) data.

The field example also shows that acoustic attenuation of the total wave field as well as conversion of downward-going P-wave in Stoneley waves observed on VSP data are strongly correlated with the presence of flow. The HES site encloses more than 35 boreholes, but only results from the well C1, located in the middle of the site, are selected for the presentation below. Several other wells at the HES showed results similar to that in C1 [24, 25].

Figure 15 exemplifies via the well C1 the type of VSP recorded at the HES. For the VSP acquisition, the seismic source is a lightweight dropper (Fig. 2a), and the borehole sensor is a hydrophone. The sampling step is 2.5 m for an acquisition performed over the 22.5–117.5 m depth interval.

The VSP data are highly corrupted by Stoneley waves. A conversion of downward propagating P-waves into downward and upward propagating Stoneley waves was observed at the level of the karstic bodies. This phenomenon occurs in highly permeable formations [15, 24, 26] and would indicate the presence of flow to be confirmed by Production Logging Tool (PLT) data [24, 25].

We can notice, for example, the conversion from P- to Stoneley waves at 55 m depth (Fig. 15a). The first arrival which is the downward moving P-wave is strongly attenuated at 55 m depth. The P-wave being partly converted to a downward moving Stoneley wave which is then reflected at the bottom of the well. The VSP data were processed to extract the downward-going and upward-going Stoneley waves. The upward-going Stoneley waves are shown in Figure 15b. At 55 m depth, we can clearly see an upward-going Stoneley wave, created by the conversion of downward-going P-waves. The Hilbert transform applied to the different wave fields allows for estimating their amplitude (instantaneous envelope). The instantaneous amplitudes of the upward-going Stoneley waves were stacked in a small corridor located after the arrival time of the downward-going P-wave, with the aim to infer a P-wave to Stoneley wave conversion factor which indicates a karstic level between 50 and 55 m depth. Assuming that the conversion factor is linearly proportional to water fluxes, it can be integrated over depth from bottom to top to mimic a flowmeter and compared with a PLT log. The integrated conversion factor is normalized by its highest amplitude value to obtain a VSP flow log ranging between 0 and 1 (Fig. 15c, center).

The normalized VSP flow (Fig. 15c) has a low vertical resolution, the sampling interval over depth being 2.5 m. Figure 15c (left part) shows a comparison between PLT and VSP flow logs. Despite the difference in vertical resolutions between PLT and VSP flow logs (1 cm for the PLT, 2.5 m for the VSP), both profiles are in good agreement and confirm that there exists an active flowing karstic body at 52 m depth as detected by the total acoustic signal attenuation log [27].

8 Conclusion

The hybrid seismic imaging tool showed that seismic data obtained from shot points with a single geophone per trace and small distances between geophones (1–5 m) is a valuable tool rendering information regarding the reflectivity for targets located in the near and/or very near subsurface. Based on a four-step procedure, the processing of refracted and reflected waves provided two sections, which after merging produced, first, an extended time reflectivity section starting from the surface and, second, a section over depth after calibration with VSP and acoustic data.

The overall procedure for carrying out FWAL experiments is relatively simple and affordable, but the scale investigated does not exceed the close vicinity of the probed well. Acoustic wave monitoring calls for a rigorous signal processing because of the multiple type of waves recorded within a wellbore by receivers. The FWAL and hydrophone VSP techniques locally complement 3D surface seismic in detecting fractured bodies or small karstic bodies not revealed by large-scale seismic wave propagation. It also renders by-products such as “synthetic” flow logs that can be compared to actual ones, and information that allows for constraining classical loggings when the question is to identify flowing and non-flowing objects crosscut by a wellbore. In that sense, FWAL might also improve the conditioning of 3D seismic raw data on borehole logs to better image fractured or karstic bodies at the large scale.

Coupling 3D seismic and FWAL is probably not a suited combination to image fractured or karstic aquifers at the regional scale, to fix the ideas, over territories of several km of horizontal extension. But as shown with the presented study, the coupling produces images of sufficient resolution for the accurate delineation of fractured or karstic bodies at a scale compatible with the security perimeter usually surrounding subsurface catchments for water supply. In view of the short transit times associated with preferential flow paths in karstified aquifers, and the subsequent risks of rapid contamination of catchments, it makes sense to envision 3D seismic and FWAL as valuable tools for assessing the effectiveness of such perimeters.

Geophysical investigations of the near subsurface are increasingly appealing to the Hydrological community, mainly for their ability to image under-sampled systems usually only visible via wells. The feature is here exemplified by high-resolution seismic data able to probe local variations in compression-wave propagation velocities. If these velocities are processed to transform them into pseudo-porosity values, seismic data image the widespread but discontinuous water bearing bodies of a fractured carbonate formation or a karstified limestone aquifer.

Acknowledgments

We thank APEC (Association Pédagogique et Expérimentale du Cher) and IFP-School for granting us their permission to use the seismic data. We thank APEC for granting us the permission to use full waveform acoustic data. Special thanks to Anthony Hardy (Géocentre Forsol) for the description and analysis of the cores. The processing of seismic data (2D 3D seismic imaging and refraction tomography) has been done with the SPW software (Parallel Geoscience). We thank Dan Herold for his advice in the processing of the data. The time picking of acoustic sections has been done with Earth-Quick software.

References

All Figures