Characterizing ﬂ ow in the ﬁ rst hundred-meter depth of a fractured aquifer using hybrid seismic methods, acoustic logging, and ﬂ ow-log measurements

. Understanding subsurface ﬂ ow, especially in fractured rocks only housing water through a few preferential pathways, is still challenging. The point is mainly associated with the poor accessibility of the subsurface and the lack of accurate representations for both heterogeneity and spatial distribution of water bearing bodies. This notwithstanding, highly-resolved geophysical investigations bring new images of the subsurface. This is exempli ﬁ ed over a fractured limestone aquifer at the site scale (for example, that of the radius of in ﬂ u-ence of an extraction well). On an experimental site, situated in the Cher region (France), two boreholes have been drilled for ﬁ eld experiments. Full Waveform Acoustic Logging (FWAL) and seismic experiments were con-ducted. Hybrid seismic imaging, which consists in combining refraction and re ﬂ ection seismic results, has been carried out. Based on a four-step procedure, the processing of refracted and re ﬂ ected waves provided two sections. After assemblage, these sections produced in a ﬁ rst step an extended time re ﬂ ectivity section starting from the surface and, in a second step, a section over depth after calibration with Vertical Seismic Pro ﬁ le (VSP) and acoustic data. However, even the Very High Resolution (VHR) seismic methods do not have a suf- ﬁ cient vertical resolution to describe accurately the geological formation. The acoustic sections were processed to separate the different wave ﬁ elds, to extract the criss-cross events and to build a criss-cross index log. A log of fracturation index, based on both criss-cross index and P-wave velocity measurements, was computed to detect the presence of fractures. After calibration, and under the assumption that the slower the P-wave velocity, the higher the permeability – porosity, a 3D seismic block of re ﬂ ection can inform on preferential areas where ﬂ ow 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. Finally, ﬂ ow log measurements con ﬁ rm the occurrence of ﬂ owing horizons that were previously marked by both seismic and acoustic data. Seismic and acoustic data are therefore suited to image contrasted hydraulic properties over fractured subsurface systems usually poorly documented.


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.

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

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].

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. 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. 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]:

VSP and 3D seismic imaging
(i) Time shift of the VSP section to flatten the downgoing wave; (ii) Extraction of the down-going wave by SVD filter (first eigensection); (iii) Computation of a residual section which is the difference between the VSP section and the estimated down-going wave; (iv) Time shift of the residual section to flatten the up-going wave; (v) 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.

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 ls 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.

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).

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: 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.

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]: where Dt is the measured sonic "travel time" (inverse of a velocity), Dt ma is the theoretical sonic travel time for the rock matrix, and C is a calibration coefficient. The Dt 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 Dt ma = 212.1 ls/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. 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.

Discussion
After calibration, and under the assumption that the slower the wave velocity, the higher the permeabilityporosity, 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: 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 pseudoporosity 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.

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 downwardgoing 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 upwardgoing 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 downwardgoing 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].

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 pseudoporosity values, seismic data image the widespread but discontinuous water bearing bodies of a fractured carbonate formation or a karstified limestone aquifer.