Multi-Scale Characterization of an Heterogeneous Aquifer Through the Integration of Geological, Geophysical and Flow Data: a Case Study

This paper presents the integrated modelling study of an experimental hydrogeological site recently developed near Poitiers city. The objective assigned to this site is the design of reliable methodologies to simulate flow and transport in highly-heterogeneous reservoirs. Around 30 wells have been drilled on this low-depth captive aquifer, according to a square and centred pattern, with an average well spacing of 70 meters. All the observations and measurements acquired on this Site lead to a characterization of the spatial distribution of the aquifer heterogeneities, and to the recognition of their probable origin. Firstly, pumping tests and production logs revealed the heterogeneous flow behaviour at the Site scale. A static characterization of the reservoir was then undertaken from geophysical surveys and especially from the joint geological analysis of cores and logs, while taking into account flow measurements. Most information were found to be consistent with a paleo-karstic origin of this aquifer, an assumption that was recently corroborated by wellbore images. The first step of this project ends up with the construction of a geostatistical model that illustrates the spatial organization of the main drains and their link with the background facies. This model, once enriched by recent data, will constitute the support model for simulating and interpreting the pumping and interference tests. This second step, i.e. the dynamic simulation of the Site, will be indispensable to constrain the location and size of the drains that ensure the productivity of each well and their hydraulic connections.


INTRODUCTION
Many underground aquifers were developed as experimental sites during the past decade. These sites are designed for insitu measurements and calibration of flow, transport and/or reactions in underground reservoirs that are heterogeneous by nature. Some of them are considered as potential repositories for nuclear waste, such as the Yucca Mountain site in the USA (Liu and Sonnenthal, 2003) or for geothermal energy production, such as the Soultz-sous-Forêts project in the Rhine Graben, France (Bruel, 2002). However, at the present time, they all remain in situ research laboratories where experiments are performed to design and test various innovative techniques and methodologies for future industrial applications. Among such applications, one may quote the development of new geophysical methods for the detection of conductive paths in aquifers (Hunt and Worthington, 2000;Hubbard et al., 2001), or remedial methods for contaminated near-surface aquifers taking advantage of the underground microbial natural activity (Nielsen et al., 2006).
The Experimental Hydrogeological Site (EHS) studied in this paper was also developed as a research laboratory aiming at simulating flow and transport in highly-heterogeneous, i.e. fractured and/or karstic, aquifers and also on the design of water quality indicators. The present publication concerns the design of a modelling methodology of flow heterogeneities that starts with the integration of available direct and indirect information from geosciences at various scales. This methodology is to be compared to other research teams approaches, whereby the measured flow data in time and space are directly interpreted with conceptual flow models, that take into account a fractal behaviour of the reservoir for instance (Delay et al., 2004).
The most straightforward industrial applications of the EHS research studies concern the survey and monitoring of heterogeneous drinkable aquifers, and also "field-scale" remedial methods for near-surface aquifer decontamination, from nitrates for instance.
The site under consideration is located on the borderline, named "the Poitou threshold", between the Paris and the Aquitaine sedimentary basins. The structural setting consists in a low-depth platform overlying a granitic basement.
This aquifer consists of about 100 meters of tight diagenetized carbonates of Middle Jurassic (Dogger) age, including Aalenian, Bajocian, Bathonian and Callovian stages. It is separated from the granitic basement by 45 m of Liassic marls and lies 20 to 30 meters below a superficial impermeable layer made up of argillaceous and organic-rich weathered deposits. The studied area, denoted as the Experimental Hydrogeological Site (EHS), covers an area of 210 m per 210 m. Well location is shown in Figure 1.
In spite of a minimum well spacing of only 50 meters, the long-duration pumping tests performed after drilling revealed surprisingly-contrasted flow behaviours from one well to another. Indeed, because the carbonate sediments of this water-bearing reservoir were deposited in the shallow marine environment of a platform margin, one does not expect any significant variation of facies or deposits thickness at the Site scale. However, as will be unveiled throughout this paper, this relatively uniform porous carbonate series was transformed into a complex karstified multiple-porosity reservoir. Fractures do not seem to constitute the main direct origin of flow heterogeneity, but probably influenced the development of diagenetic overprints that modified the texture and nature of deposits.
In the following, the heterogeneous distribution of fluid flows over the Site will firstly be illustrated both at the wellbore scale and the Site scale. This flow heterogeneity justified a multi-scale characterization of the reservoir to be undertaken. The corresponding geophysical and geological studies constitute the main section of this paper. Eventually, the resulting set of interpreted data is integrated into a preliminary geostatistical model illustrated in the third section.

FLOW OBSERVATIONS AS EXPRESSIONS OF THE RESERVOIR HETEROGENEITY
Evidences of the heterogeneous distribution of flows within the Site are found both at the Site scale from pumping tests and at the local wellbore scale from flow profiles.

348
and logs, while taking into account flow measurements. Most information were found to be consistent with a paleo-karstic origin of this aquifer, an assumption that was recently corroborated by wellbore images. The first step of this project ends up with the construction of a geostatistical model that illustrates the spatial organization of the main drains and their link with the background facies. This model, once enriched by recent data, will constitute the support model for simulating and interpreting the pumping and interference tests. This second step, i.e. the dynamic simulation of the Site, will be indispensable to constrain the location and size of the drains that ensure the productivity of each well and their hydraulic connections. The Experimental Hydrogeological Site -Location of wells: cored wells C1 (reference central well) and C2 are shown in red; the red plain line is a correlation line between wells studied herein; the black dotted line represents the in-line profile of the seismic survey described in the next section of this paper.

Figure 2
Permeability map of the EHS: a thick black line delimits the drilled zone of interest within the Site; wells are shown as white circles; well permeabilities range from around 100 millidarcys to nearly 11 Darcys.

Pumping Tests -Overall Permeability Distribution
After drilling, each well was pumped while measuring the pressure interferences in all other wells. Such field-scale measurements constitute a very rich flow data set that will be used to calibrate the hydrodynamical model of the Site. As a very preliminary analysis, we estimated an average reservoir permeability from each well pumping results, considering the aquifer as homogeneous and infinite. Large productivity/permeability differences, in a ratio exceeding 10, are observed between wells. They are illustrated in Figure 2 showing the iso-permeability regions of the Site. Note the existence of a low-permeability region in the South-Eastern region of the Site, with well pumping permeabilities of one to a few hundred millidarcys, close to 100 times less than in the Western region. Obviously, the assumption of homogeneity is not valid at all for this reservoir. Furthermore, the heterogeneity of permeability shown in Figure 2 represents a smoothed image of the local distribution of flows which may still be much more contrasted.

Production Logging
As a matter of fact, the flow profiles that were measured along most wellbores revealed that water production was not uniformly distributed along the wellbores but concentrated at very few given depth intervals, generally one or two, more rarely three. The 6 production profiles shown in Figure 3 clearly demonstrate how these flow drains are unpredictable: -the three wells M2, M21 and M22 are neighbouring wells distant by only 50 or 70 m from one another, however they produce water from drains located at various depths, around 65 m for M2, 85 m for M21, 115 m and 80 m for M22; -on the contrary, wells M3 and M5 show very similar flow profiles with a unique water entry point located between 110 and 115 m, although those two wells are far from one another, 210 m apart; -well M16, located at an intermediate distance between wells M5 and M2, produces via three drains whereas M5 and M2 produce via a single drain. Therefore, both the location and number of drains feeding wellbores can change considerably within distances less than the minimum well spacing, equal to 50 meters for this Site. The vertical extent of a given drain also seems to be variable, ranging from less than 5 m to 15-20 m.
Considering the 16 wells where a production log was measured, water entry points are encountered at various depths ranging from 115 m to 50 m, with the 80-90 depth interval being the most productive level. Figure 4 shows the normalized flow profiles along two wellbores, M7 and M3, with a superimposition of the caliper log. The production of these two wells comes from a single major drain found respectively at around 85 m and 110 m. The normalized caliper logs show wellbore enlargements that are due to the presence of large vugs, even caves, at some levels of this limestone reservoir. The well diameter, that is about 22 cm, can increase up to 40 cm at those levels and the vertical extent of such caves, generally decimetric, can reach one meter for the largest ones. It is worthwhile pointing out that a drain is most often associated with the presence of a vuggy/karstified rock but that the reciprocal assertion is not true. For instance, contrary to well M7, the diameter of well M3 shows many enlargements distributed along most of the reservoir interval while water production occurs at only one depth interval, between 115 and 105 m. These observations are valuable indices that will sustain the conceptual model proposed later on, i.e. selective drains crossing the karstic diagenetized layers or bodies found in this limestone reservoir.
The role played by highly-conductive drains is not only suspected from wellbore observations but also from largescale flow measurements over fairly long time ranges, as briefly illustrated hereafter.

Interferences
Unexpected interferences were observed in some wells while testing a given well. The water level can actually decrease rapidly and considerably in some wells located far away from the pumped well, contrary to other nearby wells showing delayed and lower responses. For instance, Figure 5 shows that, after pumping in well M7 during 57 hours, a drawdown of less than 2 meters is observed in well M8 located at a distance of 50 m from M7 whereas a higher drawdown, of more than 3 m, is recorded in M5, a well three times as far from M7 as M8.
A preliminary diagnosis of this behaviour can be attempted through a qualitative analysis of transient well responses.

Transient Flow Behaviour while Pumping
Considering again the same three wells, M7, M5 and M8, the drawdown evolution with logarithmic time (Fig. 6)  Production logging of various wells of the Site: a few water entry points ensure wells deliverability. Production and caliper logs of Wells M7 and M3.
of two water-feeding media with contrasted properties. However, the slope of these two segments differs contrary to well-differentiated dual-medium systems where only one medium ensures fluid flow at the large scale: in the present case, the late-time period shows a faster drawdown increase than would be observed for a dual medium. This observation can be attributed to the fact that the overall permeability of the draining system decreases with time, i.e. when farther regions from the pumped wellbore are sollicitated. -Interference wells responses differ very much during the initial period of pumping: the interference well M5 shows a quasi-immediate drawdown whereas a significant drawdown is observed much later in well M8 although located much closer from M7 than M5: this observation can be attributed to the presence or absence of highly-conductive flow paths ensuring the hydraulic communication of interference wells with the pumped well. -Two important observations concern the late-time flow period: -first, the drawdown does not show any stabilization with time, even after nearly three days of pumping at a significant rate, that is, a permanent regime of the wells is not yet established, apparently because the pumping goes on sollicitating farther and farther regions from the wellbore; -however, at that time, all wells show parallel evolutions of drawdown with logarithmic time: this may result from the fact that all the wells are at that time included in a unique large-scale draining system, which includes the Site but is still expanding with time.

Conclusion
Available flow information confirm the heterogeneous distribution of flows both at the Site scale and at the local scale along wellbores. An extensive geological characterization of the Site was thus undertaken, combining preliminary geophysical surveys, a detailed core-log analysis of reservoir facies in reference wells and a study of correlations between wells.

MULTI-SCALE MULTI-DISCIPLINARY RESERVOIR CHARACTERIZATION
The study of heterogeneity distribution over the Site, still under way, involves a multi-disciplinary approach confronting observations and measurements of different resolution scales. Wellbore data are the main source of information, however they do not constitute sufficient 3D modelling constraints because of the high and hardly-predictable variability in reservoir lithology and properties from one well to another. That is, the inter-well space definitely has to be explored, thus justifying the special effort made to implement the geophysical methods described hereafter.

Multi-Scale Geophysical Surveys
The seismic and acoustic methods, based on wave propagation, can be used to describe a geological structure, a reservoir zone or an aquifer at different scales. The surface seismic method (reflection survey) has a poor vertical resolution but a huge lateral investigation in comparison with the well- bore seismic methods. The wellbore methods, i.e. Offset or Vertical Seismic Profiles (VSP), have a restricted lateral investigation but a much finer vertical resolution than the surface seismic method. These differences are still enhanced for the acoustic logging that has a very high vertical resolution but a very poor lateral investigation.
The vertical resolution is usually evaluated as the quarter of the dominant wavelength and the lateral resolution as the diameter of the first Fresnel zone. In practice, the full processing sequence, as described hereafter, improves the lateral resolution up to the vertical resolution. As a guide for the application of seismic methods to this low-depth aquifer, the following table summarizes these points: The acquisition of usable seismic data is particularly difficult for low-depth reservoirs underlying a thick weathered zone, such as the aquifer studied herein. Different surface seismic surveys were however attempted with different acquisition schemes that were implemented in the vicinity of central reference well C1. A refraction survey has been recorded to image in depth the limit between the weathered zone and the unweathered carbonate formation, close to the aquifer upper limit. Some seismic shots have also been recorded to assess reflection seismic as a suitable method to detect and image reservoir internal markers. However, the geological or petrophysical identification of such markers requires higher resolution methods, i.e. well seismic surveys and acoustic logging. These two calibration methods were implemented in well C1 and a few neighbouring wells.

Surface Seismic -Refraction Survey
Two profiles were recorded, a first profile, called the in-line profile, coarsely oriented SW -NE (Fig. 1), and a second cross-line profile, orthogonal to the first one. The two profiles are crossing in their middle, in the vicinity of well C1. A detonating impulse source was used. The recording spread consisted of 48 single geophones distant by 5 m from one another in the in-line profile, and 96 single geophones distant by 2.5 m in the cross-line profile. Whatever the profile, two shots (a direct shot and a reverse shot) were fired at the extremities of the recording spread. The time sampling interval was 0.25 ms and the recording length 0.5 s.
The picked times of the refracted arrivals along each profile (Fig. 7, left and middle parts) enable to get an image in depth of the low-velocity superficial zone contrasting with the underlying high-velocity carbonates. This rough refraction surface (Fig. 7, right part) is a clearly-detectable marker but would not correspond to the aquifer top surface. Actually, the latter was invariably encountered at a depth of -30 m according to the three neighbouring piezometers, PZ1, PZ2 and PZ3, drilled along the in-line profile, whereas the refractor fluctuates within a large depth range, between -10 m and -30 m (Fig. 7). One assumption, yet to be confirmed, is that the measured refraction surface would correspond to the bottom limit of a residual paleo-karst locally infiltrating the Results of the refraction survey: evidence of an irregular marker above the aquifer. limestone aquifer formation via caves or lineaments. During drilling operations, such caves or lineaments were indeed encountered and found to be filled in by clayey surface erosion material.

Surface Seismic -2D Reflection Survey
Data were recorded on the in-line profile used for the refraction study, with the same recording spread, composed of 48 single geophones. Two sources, either a detonating impulse source or a mini vibrator system, were tested but the latter was discarded because of an insufficient frequency spectrum. Figure 8 shows a shot point obtained with the impulse source. Data are presented both as frequencies versus wavenumbers and as times versus distances. Undesired energetic wave fields (Fig. 8, left section), composed of low and high apparent velocity surface waves and of direct and refracted body waves, mask possibly-recorded waves emanating from reservoir markers. In such a situation, waves had to be separated sequentially thanks to a dedicated procedure involving the selection of suitable filters (Mari, 2006). Although addressed to seismic treatment specialists, this procedure is summarized hereafter, to underline the difficulties of near-surface seismic surveys.
The elimination of an undesired wavefield is performed with the following procedure: -the extraction of the considered wavefield by a specific filtering method, either in the frequency-wavenumber domain (F-K filtering), or by Singular Value Decomposition coupled with Independent Component Analysis (SVD-ICA filtering); -the subtraction of the extracted wavefield from the input section to obtain a residual section; -the residual section becomes the input section for the following step.
For the near-surface objective of our survey, the above separation procedure had to be implemented four times sequentially. The first step consisted in eliminating by a F-K filter the direct wave and the slow (low apparent velocity) near-surface Rayleigh wave. In a second step, the same procedure was applied with a F-K filter to eliminate the fast (high apparent velocity) Rayleigh wave. The residual section shows negative apparent velocity events that are converted refracted waves, which were eliminated by F-K filtering in a third step. The last step consisted in extracting the refracted wave by SVD-ICA filtering. Independent Component Analysis improves the separation of waves obtained by the conventional Single Value Decomposition technique. Details on the SVD-ICA technique can be found in Vrabie et al. (2004), with synthetic data set examples.
The residual section thus obtained clearly shows infinite apparent-velocity events which are associated with reflected waves (Fig. 8, middle section). Corrections for source-detector offset (Normal Move-Out) and resolution enhancement have finally been applied to this residual section. The seismic section thus obtained is shown in Figure 8 (right part). It shows a few reflected waves in the [0-100] ms time range that investigates the entire reservoir section, in particular a reflected wave at 75 ms.
Then, these reflected waves have to be identified and depth-calibrated by comparison with well data, namely Vertical Seismic Profiles (VSP) and acoustic logs. Results of the reflection survey: hardly-detectable reservoir markers as the result of a specific treatment procedure.

Vertical Seismic Profile (VSP)
Vertical and slightly-Offset Seismic Profiles were acquired in the reference well C1 in order to identify major reservoir markers and calibrate the surface seismic sections. The weathered zone obstacle was overcome thanks to buried sources. The frequency bandwidth of received signals was large and reached 800 Hz. Figure 9 (left part) shows the deconvolved upgoing reflected waves of the VSP, both in the entire 30-840 Hz bandwidth and in the limited 20-140 Hz bandwidth of a conventional surface seismic survey. Obviously, the latter delivers an impoverished signature of possible reservoir markers. Figure 9 (right part) shows the VSP data inserted in the surface seismic section of Figure 8 (right-hand section). VSP data allow the identification of seismic reflectors on the surface seismic section, in particular the event observed at 75 ms. The latter turns out to be a reservoir marker located at a depth of about 110 m, and would then correspond to the Aalenian-Bajocian transition, close to the aquifer bottom limit. Other internal reservoir markers are also observed above, a priori within the Bajocian-Bathonian unit.
That is, this well seismic survey validated the surface seismic results and pointed out the interest of carrying out a 3D surface seismic survey to obtain an image of the aquifer. For that purpose, an experimental spread, described hereafter, was tested.

Feasibility of a Full (3D) Seismic Coverage of the Site
Inland 3D seismic surveying generally involves a cross spread, i.e. shot lines that are perpendicular to the receiver lines. In our case, a simple experimental spread was tested. The same recording spread, 48 geophones, as for the refraction study was used on the same in-line profile. A single shot line was also implemented. It was located on the cross line profile used for the refraction study. A single shot has been recorded and processed. For this feasibility study, the source was a light vibrating system. Because the latter generated a strong air wave, the distance between the source point and the receiver line, of 110 m, was chosen in Results of the well seismic survey: calibration of the reflection survey. order to avoid the interference between this air wave and the useful body waves. The processing of the 3D shot provided a seismic section parallel with the 2D seismic line described before, but laterally shifted by half the sourcereceiver line distance, i.e. 55 m. Fig. 10 shows a composite section, comprising: -the 2D section obtained before with the impulse source for distance abscissas ranging from 0 to 50 m, -the parallel section obtained with the vibrator source for distance abscissas ranging from 50 to 300 m. Both sets of data correlate well together in spite of the shift between those 2D sections. That is, a 3D seismic reflection survey looks as a feasible alternative to get a 3D image of the aquifer over the Site area, provided however that a low spacing between recording lines is adopted.
To end with, preliminary surface seismic and VSP surveys demonstrated the possibility to identify the aquifer bottom limits as well as the main markers of the aquifer. A surface seismic survey extended to the whole Site is therefore considered in order to correlate seismic markers of major rock heterogeneities possibly playing the role of drains or related to these drains. Assigning such a challenging objective to a 3D surface seismic survey assumes that markers are carefully identified from VSP data, or from acoustic logging data as described hereafter.

Acoustic Logging
Acoustic logs constitute a second and higher-resolution method to calibrate surface seismic surveys. Full waveform acoustic profiles were recorded along wellbore C1 and four neighbouring wells, MP6, MP5, M8 and M9. Full waveform acoustic profiles of well C1 and neighbouring wells.  The acoustic tool was a monopole tool designed and manufactured by S.E.M.M. service company. It includes one transducer and two receivers. The distance between the source and the nearest receiver is 3 m and the distance between the two receivers is 0.25 m. Source and receivers are multidirectional. The transducer generates in the fluid a compressional wavefield which is conveyed in the formation as a compression wave (P wave) and a shear wave (S wave) at the refraction limit angles. More precisely, in a vertical well, such a tool permits the recording of five propagation modes: the refracted P wave, the refracted S wave (only in fast formations -definition below), the fluid wave, and 2 dispersive guided modes (the pseudo-Rayleigh and the Stoneley waves). S waves can be generated only if the formation S-wave velocity is higher than the P-wave velocity in the mud; the formation is then called a fast formation (contrary to the socalled slow formations). Figure 11 shows a representation versus depth of the recorded waves that are at first the refracted compression wave along the wellbore wall, then the refracted shear wave, and finally the wellbore-conveyed Stoneley waves. The amplitudes of the different waves are coded in different colours.
Such a representation allows to classify the geological entities in acoustic facies. For instance, considering the reference well C1, three main acoustic units could be defined by the following depth intervals: the 32-60 m unit including a discontinuity at 50-52 m, and two heterogeneous units at 60-85 m and 85-110 m that show a fairly similar acoustic facies and are hardly distinguished from one another. However, the sharp velocity change seen at 85 m on the P-wave velocity log also shown in Figure 11 confirms a major acoustic discontinuity at that depth.
The acoustic profile also gives a qualitative information about the internal heterogeneity of identified units. Whereas wells M8 and M9 show fairly unaltered waves, wells MP5 and MP6 show perturbed and discontinuous acoustic profiles. As a matter of fact, wells MP5 and MP6 reveal themselves as prolific wells contrary to wells M8 and M9 that have a very poor water deliverability.
Furthermore, the Stoneley wave response is strongly related to the state of continuity of the well wall along which it travels. In practice, the disappearance of the Stoneley wave at a given level along the wellbore signs the presence of an open fracture or joint, or of a cavernous layer, as it can be seen here at a depth of about 50-52 m for well C1, and apparently also at depths of 50-53 m and 75 m for well MP6.
Quantitatively, the picking of the arrival times of the different waves gives access to velocity logs. The time shift between the two receivers, distant by 0.25 m from one another, was measured during this acoustic logging. Figure

Figure 12
Acoustic velocity profiles of well C1 and neighbouring wells (note that the sharp discontinuity visible in all 5 wells, at a depth varying between 21 m for well C1 and 38 m for well MP5, corresponds to the bottom limit of the steel casing set in the upper section of each borehole).
To conclude, acoustic logs give a local image of the major reservoir units and of their internal heterogeneity. Moreover, they confirm the presence of major discontinuities, that are in our case related to the presence of vugs and caves. However, such discontinuities may or not reveal themselves as effective water entry points for the well, depending on their connection to the major flow network far from wellbore. In that respect, the production log of well MP6 showed a major water arrival at a depth comprised between 74 and 57 m, which would exclude the 50-53 m discontinuity as a water feeding point for this well, whereas well C1 flow profile confirmed the 50-52 m discontinuity as the single major water entry point for that well.
The main conclusion of this geophysical study is that a preliminary surface seismic survey calibrated from a well seismic survey and logging demonstrated the possibility to identify major markers within the aquifer or limiting it. A 3D surface seismic survey extended to the whole Site is therefore considered along with its calibration from well VSPs and acoustic logs. Expected information is the correlation over the entire Site of major discontinuities or rock heterogeneities possibly playing the role of drains. Due to its limited resolution, hardly less than ten or fifteen meters, one cannot expect that the flow paths be so accurately identified from the surface seismic method as from production logs. Then, the minimum objective assigned to a 3D seismic survey is the detection and 3D visualization of the major caved or vuggy bodies. The assessment of the size and connectivity of such presumably-conductive bodies is essential to build a geostatistical model of the Site, as described later on. In this respect, 3D seismic-derived information may help in defining the variogram parameters and the main reservoir regions reflecting the distribution of those conductive bodies and/or their background facies.

Geological Characterization of the Site
The reservoir facies were characterized from core observations of two wells, from a limited suite of logs and also drilling reports. Objectives were first to characterize the reservoir facies and understand their origin and genesis, then to establish correlations between wells for the subsequent construction of a 3D model of the Site. In addition to the recognition and correlation of primary sedimentary features, a special attention was devoted to the post-deposition phenomena that altered the lithology, texture and properties of deposits.

Sedimentary Facies and Depositional Sequences
Reservoir facies as well as depositional sequences, as proposed by Cross (1988) and Van Wagoner et al. (1990), were characterized thanks to a thorough analysis of the cores from two wells C1 and C2, distant by 300 m and located respec-tively at the centre of the Site area and in the South-East. Several stratigraphic markers of Maximum Flooding Surface (MFS) and of sequence boundary were identified on the basis of macro/microfacies observations, with a specific attention to the lithology, the sediment texture and granulometry, the biogenic or detrital origin of figured elements, and also the nature of sedimentary surfaces, either bioturbated or erosive. Except for ammonites found in well-C1 and well-C2 cores at the top of Liassic shales that give a Toarcian age after E. Cariou (personal communication), the proposed ages and sequences come from the calibration of cores to regional literature data (Gonnin et al., 1992;Gonnin et al., 1994;Gabilly et al., 1997).
The following vertical succession of sedimentary facies and depositional sequences was thus identified by combining both wells observations (Fig. 13): -The granitic basement was reached in well C2 alone. It was emerged and deeply eroded from the Late Paleozoic time; -The first sedimentary deposition overlapping this basement consists of 45 meters of Liassic marls including: -a 8 meter-thick alternation of azoic shales and dolostones, dated Hettangian-Sinemurian in the Vienne district area (after Gabilly et al., 1997), and ending with a strongly-bioturbated surface; -mudstones and marls with pelagic and benthic fauna revealing a clear deepening of the depositional setting. Bioclastic marls with common bivalves and belemnites, observed between 150 m and 135 m in well C2, are supposed to be Pliensbachian outer-platform deposits. A Toarcian age of deeper-deposition marls was confirmed and correlated between wells C1 and C2, thanks to ammonites found between 130 m and 127 m in well C1 and between 130 m and 122 m in well C2. -Overlying Liassic marls, two deposition sequences attributed to the Aalenian formation: -a 10-meter thick shallowing-up sequence with a clear lithological change from shale to carbonate that is associated with the Aalenian climate change (after Gonnin et al., 1992). The bottom limit of the studied aquifer is assumed to lie within this sequence, i.e. at a depth of 120 to 125 m. Sedimentary facies and depositional sequences from core observations of wells C1 and C2 (the Fracture Intensity log is also reported).
highly-bioturbated bioclastic wackstone and packstone rich in skeletals (bivalves, crinoïds, sponges, corals, bryozoans), while regressive cycles are underlined by coarser and cleaner facies (packstone to grainstone) that are very rich in non-skeletal grains (peloïds, ooïds and oncoïds). 3 oncoïd-rich levels (no. 2, 3 and 4) make possible a reliable correlation of the 3 sequence boundaries between the 2 cored wells. -Between 46 m and 25 m (C1) or 26 m (C2), a fining-up unit attributed to Bathonian age. It starts with a nonskeletal, oncoïd-rich grainstone and rudstone that changes vertically to a packstone and then to a skeletal wackstone that reveals a clear deepening event at the end of the deposition unit. A noticeable point however is that the Bathonian unit observed in this well lumps together the five condensed or eroded depositional sequences, named Bt 1 to Bt 5, that are usually reported in adjacent areas for the Bathonian stage (Gonnin et al., 1994). This observation supports the assumption of a karstic origin of this aquifer discussed later on. -The carbonate column ends with a 20 m-thick, fine-grain, chalky limestone unit assumed to be of Late Bathonian to Callovian age. These well-known limestones around Poitiers city have a high microporosity and the fine texture of micro-packstones. -The last 5-6 meters of Tertiary-Quaternary shaly deposits and soils infiltrate caves that may be locally abundant in the underlying aquifer formations. In summary, the vertical succession of sedimentary facies and the thickness of deposits are similar in the two cored wells, C1 and C2. Deposition sequences can be identified and correlated, especially for the Aalenian-Bajocian time interval. That is, the origin of the heterogeneous hydrodynamic behaviour of this aquifer has to be searched for elsewhere, namely in post-depositional phenomena that deeply altered primary reservoir facies. Further investigation of this heterogeneity follows with the observation of structural features and of diagenetic fingerprints.

Structural Features -Fracturing
The Jurassic sediments constituting the reservoir were deposited in a low-depth platform environment located at the threshold of two sedimentary basins, named the Poitou threshold. At that location, sedimentary deposits are thinner than in the Paris and Aquitaine Basins and they overlie a Hercynian granitic basement linking the Vendée and Limousin domains.
The structure of this basement is influenced by N120°to N140°and N060 to N070°tectonic lineaments which have been reactivated during Jurassic, as suggested by the thickness and lithology variations of sediments (see for instance Mourier and Gabilly, 1985).
Following the sedimentation on the post-Hercynian peneplain, several tectonic events, related to distant Pyrenean and Alpine orogenesis of Tertiary age, induced a fracturing of the carbonate deposits. Because of their possible impact on the aquifer flow behaviour, a characterization of these fractures was undertaken on analogue outcrops located near the hydrogeological site, and from the observation of available cores from wells C1 and C2. -Analogue outcrops study. Outcrops show that Bajocian and Bathonian limestone formations are mainly crossed by sub-vertical swarms with a high vertical extension rather than bed-controlled diffuse fractures. These fracture swarms are oriented in three directions, N120°, NS and N50° (Fig. 14), with the N120°direction probably influenced by the pre-existing Hercynian lineaments. The spacing of these fracture sets is metric to decametric. The Callovian outcrops show a higher density of fractures, but with clay filling.
-EHS cores study. Regarding the Site, fairly few fractures were observed within the Bajocian-Bathonian sections of the two vertical wells C1 and C2 (Fig. 13). More fractures were found on the Callovian upper sections, but often filled in by a clayey or organic material due to the proximity to the surface. Typical tectonic fractures, i.e. high-dip tension and shear fractures, are rare and most often cemented. The few observed fractures, of decimetric size, are often located at the interfaces between cherts and rock matrix or within cherts. Dominant micro-structures are sub-horizontal stylolithic surfaces and compaction bands, apparently associated with burial processes. Some of them are not cemented and may constitute flow conduits. Fluid circulation seems to have influenced the fracturing process and the evolution of fracture properties, as shown by the presence of clay smearing (Fig. 15).
Although distant by only 200 meters, the two wells revealed a high fracture density contrast: -in well C1, few decimetric-scale fractures are found in the host rock and cherts and are cemented; -on the contrary, well C2 shows a high fracture intensity, particularly in the 85-95 m and the 55-25 m depth intervals, with more or less connected fractures, associated with vugs and molds. Organic matter and shale are sometimes filling those fractures. This difference is not yet fully understood. In the absence of any identified large-scale structure near by these wells, such a contrast may be related to the distribution of the diagenetic fingerprints within the reservoir. The greater presence of siliceous nodules in well C1 than in well C2 may be underlined to that respect. To conclude, typical diffuse tectonic fractures do not seem to constitute the main direct origin of flow heterogeneity. However, horizontal stylolite-associated discontinuities at bed boundaries may constitute preferential flow paths if remained open. In addition, outcrops revealed the existence of fracture corridors at a larger scale. Such discontinuities might influence the flow connectivity between wells at the reservoir scale, either as preferential flow paths or as flow barriers depending on their conductivity. They may also have played a role in the development of the diagenetic features described hereafter.

Diagenetic Fingerprints
Particular attention was given to the detection of the following diagenetic phenomena that deeply altered the porous limestones of this aquifer: -dolomitization: -local dolomitization of vuggy carbonate layers (sucrosic dolomite), in association with stylolithes that may have played the role of conduits for Mg-rich solutions; -diffuse dolomitization of the muddy limestone that mainly affected the bioturbation framework; Figure 15 Cores from C2 well at about 97 m depth: clay material is filling the upper "damaged" zone (karst residue?); a vertical fracture is observed between two marly (dark) layers in the lower core.

Figure 16
Correlation between wells C1 and C2: evidence of the high spatial variability of diagenetic fingerprints.

Figure 17
Reference well C1 data integration (fracture observations reported in Figs. 13 and 16). Texture codes (underlined): Mudstones, Wackstones, (μ(micro-))Packstones, Grainstones Grain sizes (underlined): Very fine, fine, Medium, Coarse, Very coarse, Pebbles Presence of cherts shown as black triangles on the lithological column -silica precipitation: silica is present in the form of cherts grown around calcite biological structures such as burrows or bioclasts, or in the form of cherty bands replacing calcite within formerly-permeable layers; -argillaceous deposits and clay filling in previously-diagenetized or -dissolved layers and fractures; -fracturing features; -and finally also breccia and caves, as a result of extensive dissolution phenomena. Those diagenetic fingerprints, in particular the presence of silica and dolomite, were identified and characterized in relation with the primary sedimentary features. Hereafter, the core observations from wells C1 and C2 give clear evidences of their spatial variability (Fig. 16): -Dolomitization: C1 cores reveal the presence of a vuggy re-crystallized dolomitic sub-unit between 91 and 87 m. A similar dolomitized layer is also observed on C2 cores but 1 meter thinner and embedded in a tighter limestone than in well C1. That is, dolomitization does not uniformly affect the aquifer. -Silicification: well C1 shows cherty limestones over three depth intervals, 120-112 m, 100-98 m and 56-22 m, whereas C2 reveals the presence of only one chert interval between 112 and 108 m, i.e. with a reduced thickness compared to the deepest chert interval in well C1. This high variability of chert distribution at the Site scale may contribute to the heterogeneous flow behaviour of the aquifer. -Stylolithes: they are observed in both wells, more frequently below 60 m. Stylolithes are the result of physico-chemical processes linked to sediment compaction, related to a burial of probably several hundreds of meters for this reservoir. -Fractures: whereas fracturing is intense over the entire aquifer section of well C2, fractures are only found in the first 30 meters of well C1. Some of these fractures are filled with clays, which were identified as transported clays (El Albani et al., 2005): their assumed origin would be vertisols of Tertiary age and thus, they could constitute the remaining evidence of a paleo-karstic nature of this aquifer.
To conclude, the core observations of wells C1 and C2 reveal that diagenetic fingerprints completely modify the lithology over short distances. Cherty limestones and fractured/karstic units may vanish or pinch out within a distance in the order of the well spacing. Therefore, the main origin of the contrasted well responses at the Site scale may lie in the great variability of diagenetic fingerprints that affected this limestone formation. The confrontation of static and dynamic observations made in reference well C1 further confirms this assumption.

Reference Well C1 Data Integration in Relation with Fluid Flow Behaviour
The central well of the Site, C1, disposes of the most complete dataset among all the Site wells, with detailed core descriptions from the surface down to 129 m in depth, a large set of wireline logs, including gamma-ray, resistivity, acoustic and borehole image logs, petrophysical data, and a flow profile. All these information sources were confronted together to derive possible indicators/drivers of lithological changes and of flow properties along the wellbore. Figure 17 shows the confrontation of all information sources as well as the proposed units, labelled A to M. Core observations of grain size and rock texture, but also the acoustic, gamma-ray and electrical logs responses, reveal the heterogeneity of Bajocian units. The upper and lower limits of this formation seem to coincide with major resistivity changes. On the contrary, the Bathonian facies look more uniform in terms of texture, except for the presence of numerous cherts.
Regarding fluid flows, the major water arrival observed on the production log is found at a depth of 50 to 52 m and is responsible for about 75% of the well deliverability. This water entry point lies at the base of a cherty limestone unit. As mentioned before, the presence of this drain at that depth is consistent with the full-wave acoustic log showing the disappearance of Stoneley waves in relation with the presence of a fluid-filled discontinuity. According to the recorded flow profile, an entirely-dolomitized vuggy layer at about 90 m and another cherty section, 110 to 120 m deep, may also contribute to fluid flow, however this contribution is not so obvious as for the water-entry point at 50 m.
Petrophysical measurements on well-C1 cores. The matrix porosities and permeabilities were measured on 31 core fragments distributed along the wellbore. Porosity values are ranging from 7 to 25% with an average of 15%, and permeability values from less than 0.01 md to some 200 md. Most Bajocian permeability values are comprised between 0.1 and 1 md, around one order of magnitude less than Bathonian-Callovian values. However, high matrix permeability contrasts, of up to three orders of magnitude, are observed locally within each unit.
The porous structure of these carbonates was characterized from Nuclear Magnetic Resonance measurements. It varies significantly between samples, from unimodal to bimodal pore-size distributions, as illustrated below by the time-T2-NMR-responses of the samples distributed along the wellbore (Fig. 18).
Conclusions regarding the local hydrodynamic behaviour are that: -water is feeding well C1 via a single entry point, that is a vuggy unit (D) located in a cherty limestone background, -other vuggy/dolomitized and/or cherty units (L,H) may also contribute but to a lesser extent to wellbore water deliverability. However, considering the fairly-high porosities and significant permeabilities measured on some cores, the contribution of the matrix medium to the refilling of the main drains cannot be discarded although it cannot be captured from the observation of the wellbore flow profile.
More generally, the main drains ensuring the well deliverability appear to be embedded in diagenetized units, or found as locally-dolomitized levels within clean limestones: they would then effectively result from post-deposition diagenesis episodes that strongly altered the lithology and petrophysical properties of the original sediments.

Study of Correlations between Wells
Following the advent of the sequence stratigraphy methods in the late seventies and eighties (Vail et al., 1977;Cross, 1988;Galloway, 1989;Van Wagoner et al., 1990), well-to-well correlations were usually based on the wireline log signatures of given time surfaces of stratigraphic sequences, such as Maximum Flooding Surfaces (MFS) and/or Sequence Boundaries (SB). In this study, although wells are very close from one another (50 or 70 m), the correlation of stratigraphic surfaces is unworkable as strong diagenetic fingerprints completely and abruptly modify the lithology, and in consequence the wireline log signatures as well.
At the initial stage of this ongoing project, the only available data for well-to-well correlations were essentially: -the drilling reports of all wells: they provide useful but only-qualitative information regarding the lithologic changes and water arrivals in wellbore, -the gamma-ray logs available for all wells.
In addition, acoustic and production logs were also acquired in a limited number of wells. Hereafter, we analyze the well-to-well correlations between the 5 wells located near the Site centre and disposing of acoustic logs, that are MP5, MP6, C1, M9 and M8. This line of wells is drawn in Figure 1 and the results of a tentative correlation between them is shown in Figure 19. For geological interpretation purposes, this correlation study takes also advantage of recentlyacquired wellbore images -for C1, M8 and M9 -and of production logs -for C1, MP5, MP6 -as those pieces of information give a more detailed and reliable insight of the reservoir structure and behaviour. Therefore, the lithology log was built from wellbore images when available, from cores for the reference central well C1, and from the drilling reports in the absence of any core-log information. It involves a limited number of lithotypes that integrates both the sedimentary facies and the diagenetic fingerprints described in the previous section. This lithotype definition is close to the one adopted for the geostatistical model described in the last section.
A few stratigraphic correlations consisting in two Sequence Boundaries could be established between the 5 wells. They constitute the minimum evidence that depositional facies are of similar nature and thickness at the EHS scale. Otherwise, the nature and the sequence of lithotypes are found to change considerably from one well to its neighbour. In particular, cherty and fractured/karstic limestone units may vanish or pinch out within a few tens of meters, as can be observed between wells C1, M9 and M8 that dispose of the most reliable information, derived from cores or borehole images.
The location of drains is reported from available production logs, i.e. for wells C1, MP5 and MP6. Obviously, these water drains do neither follow the stratigraphy, nor the lithotype distribution. The only information regarding their distribution is that they are embedded in diagenetized limestone bodies, i.e. fractured or cherty limestones, or found as locally-dolomitized levels within clean limestones. In addition, the absence of such drains can be reasonably assumed for M8 and M9 that revealed a poor deliverability from pumping tests.
Eventually, the correlation of the main drains from one well to another turns out to be hazardous as these drains seem to result almost exclusively from local diagenetic phenomena that are not correlated from one well to another. Furthermore, the apparently-random distribution of flow paths within diagenetized units is in favour of a karstic origin of the aquifer drain system. A recent extensive analysis of borehole images did corroborate a paleokarst origin of the aquifer. Qualitative pore size distribution of well C1 cores. Time T2 is directly related to the pore size and the signal amplitude is the approximate contribution of a given pore size to the sample porosity: unimodal to bimodal pore size distributions are observed with a very variable spreading from one sample to another.

Figure 20
Bore-hole images as evidences of the paleo-karstic origin of the EHS heterogeneities.

The Proposed Conceptuel Model of the EHS Aquifer: a Paleo-Karst System
Paleo-karst systems result from a complex history, involving a succession of near-surface karst generation phases and of subsequent sedimentation and burial phases. More precisely, Loucks (1999) proposed the following genesis model: -first, a long-term exposure of carbonate sediments to the chemically-aggressive conditions of a fluctuating vadose zone leads to the development of a karst made up of superposed near-surface cave systems; -then, the onset of a new sedimentation phase leads to the burial of this karst and to the collapsing of walls and roofs of superimposed caves, with a concomitant development of breccia and fractures. The result is a paleokarst, i.e. a coalesced system of partly-collapsed caves of complex geometry with an apparently-random distribution in space. Paleo-karst facies are of peculiar nature and distribution. Loucks et al. (2004) proposed the following classification: -undisturbed host rock facies, -breccia facies, mainly resulting from the collapse of cave walls, -cave-filling sediments, generally of silici-clastic (clayey) nature. Regarding our experimental site, recent borehole images of the EHS wells revealed the presence of a complex partlycollapsed cave system. The presence and vertical location of this system change very much and cannot be correlated from one well to another. The various types of facies described above are clearly identified on those images, as shown from left to right in Figure 20: -the host rock facies, that consist of vuggy limestone, of limestone-dolostone alternations or of cherty limestone with cherts distributed as bands or nodules. Obviously, this host rock underwent several phases of diagenesis. -breccia facies resulting from the partial collapse of the paleo-caves: they are represented either by a stack of unsorted coarse clasts leaving millimetric to decimetric open vugs or caves between them, or by finer clasts embedded in a shaly matrix. -cave-filling shales, possibly laminated, that may appear in different colours on borehole images, probably in relation with different episodes of filling. Whereas the coarse brecchia facies has an autochton origin by cave walls collapsing, the cave-filling facies results from the filling of the karst conduits by allochton finer sediments of single or multiple origins (El Albani et al., 2005).
Eventually, borehole images bring a direct evidence of the suspected paleo-karstic origin of the major flow network of this aquifer.
Although some questions remain regarding the sequence and dating of the karst generation, burial and re-juvenation phases, a paleo-karstic conceptual model looks consistent with the paleogeographical evolution of the region. Indeed, the following observations support the assumption that a first karst may have developed during the Bathonian age: -the shallow marine Bathonian deposits are drastically thinner (between 20 and 40 m) on the Poitou threshold of the Paris Basin than in the central part of this Basin (about 200 m); -furthermore, the stratigraphic correlations between the central part of Paris Basin and its Southwestern Poitou margin would reveal missing sequences and erosion surfaces during the Middle Bathonian age (P. Houel, personal communication). This initial karst system would then have been buried by the thick marine sediments of Upper Jurassic age, leading to a collapsed paleo-karst. This paleo-karst may then have been rejuvenated during the Lower Cretaceous period following an erosion of the overlying Upper Jurassic cover, before another phase of burial during the Upper Cretaceous period, and other episodes of rejuvenation during the Tertiary-Quaternary eras, finally resulting in the complex filled-in paleo-cave system observed today.
A similar Bathonian paleokarst is known in the Grands Causses area, along the southern margin of the Central Massif (Charcosset et al., 2000).
Beside the geological understanding of the Site, borehole images point out some limitations regarding the detection of the water feeding points of a wellbore. Indeed, image logs turn out to be more reliable indicators of the productive reservoir sections than cores, as core sampling is often incomplete in vuggy or caved wellbore sections. Production logs must then confirm the effective contribution to flow of these potentially-productive sections.
The following section presents a preliminary geostatistical model that attempts to integrate in a consistent way the characterization data described, cross-checked and correlated all through the present section.

GEOSTATISTICAL MODELLING OF THE SITE
This model was built at an early stage of the project. At that time, wells disposed of very few pieces of information, essentially a gamma-ray profile and a drilling report giving the succession of contrasting rock units as well as water arrival events. Very few flow profiles and no other borehole images than well-C1 image log were available. This information was used to reconstruct a facies log for each well of the Site. Well-to-well correlations were attempted combining both stratigraphic surfaces and diagenetic fingerprints. Well-C1 core observations/measurements were used as the reference information to calibrate the correlated units in terms of sedimentary, lithological and petrophysical properties.
The resulting set of data and constraints was deemed sufficient to build a preliminary geostatistical model of the site.
The complex distribution of facies heterogeneities was modelled using the truncated gaussian algorithm, implemented in the 5.2 version of Heresim 3D software package. The background of this facies modelling workflow and methodology has been presented in several papers (e.g. Matheron et al., 1987;Galli et al., 1993;Doligez et al., 1999;Ravenne et al., 2002). In particular, the method offers a large choice of options to constrain modelling from external geological or seismic-derived information.

Modelling Features
For modelling purposes, the various facies defined from core observations were lumped into 8 groups, called lithotypes.
These lithotypes take into account both the primary and diagenetic features. The three main reservoir lithotypes are the limestone facies and two diagenetized expressions of this facies that are fractured limestone and cherty limestone.
Four other lithotypes represent pure argillaceous or clayey facies: clayey limestone, marl, clay and dark clay. Note that the latter facies, dark clay, can be considered as a diagenetic facies as it mostly constitutes a fill-in sediment of caves and fractures.
The eighth facies, denoted as "Water", represents the conductive facies levels, as detected in wells during the drilling operations or from the few flow profiles available at the time of model construction. This "Water" facies is only defined on Vertical Proportion Curves (VPC) of lithotypes in the different regions of the geostatistical model (distance is given from a reference level that is the soil surface). the basis of its contribution to flows and not according to its lithology or its petrophysical properties. It is most often found within cherty or fractured limestone bodies.
Wellbore facies data revealed that the presence or absence of cherts was the sole major factor of spatial variability of lithotypes that could be quantified at the Site scale. The vertical distribution of lithotypes was then determined for two groups of wells showing either abundant or scarce cherty levels. Figure 21 shows the Vertical Proportion Curves (VPC) for those two groups of wells and for all the wells together. The probability of presence of the conductive "Water" lithotype at a given depth is clearly correlated to the presence of diagenetic cherty or fractured limestones, themselves embedded in the virgin limestone facies.
The Site was then mapped according to the silicification criterion, with three regions that were "cherty", "noncherty", or "undefined" in the absence of any wellbore information. Of course, once other logs will be acquired and new wells be drilled, other heterogeneity drivers than the presence of silica might be integrated into the geostatistical model. For the present model, one can state that the delimitation of the Site regions was based on a diagenetic factor, that is the intensity of silicification phenomenon. The latter may have taken place along the main paleo-karst conduits and thus constitute a possible indicator of the present flow paths, but not exclusively as other indicators, such as the presence of a leached or dolomitized facies or more directly the occurrence of caves or fractures, may also explain the flow distribution.
Variograms were built to analyze and quantify the spatial continuity of the previously-defined lithotypes. Vertical variograms are shown in Figure 22. The variogram ranges of the carbonate and/or clayey facies having a primary sedimentary origin cannot be defined as they exceed the hectometric scale of observation. On the opposite, the variogram ranges of the diagenetized facies are much lower, in the order of one or a few decameters, about 10 meters for the "Water" lithotype and around 30 meters for the cherty limestone. These variogram characteristics are consistent with wellbore observations, that show a very limited vertical extent of conductive bodies, and cherts distributed over pluri-decametric sections in most wellbores.
No meaningful horizontal variogram could be defined. Actually, the frequent discontinuity of diagenetized facies between neighbouring wellbores emphasizes the nugget effect. To account for this random continuity or discontinuity of diagenetized facies from one wellbore to another, the range of the horizontal variogram was estimated to be in the order of the well spacing. However, as productivity contrasts between neighbouring wellbores are less often observed in the North-South direction than in the West-East direction (Fig. 2), a higher range value was taken in the North-South direction, i.e. 100 m instead of 50 m the East-West direction. These values were used to parameterize a gaussian variogram model that reflected essentially the spatial distribution of the conductive "Water" lithotype.
Using the regionalized Vertical Proportion Curves of lithotypes and the gaussian variogram described above, the distribution of facies was simulated within a 3D grid of the Site. A random truncated gaussian function was used with probability thresholds derived from the proportions of each lithotype at the reservoir location under consideration. These lithotype proportions are read on the Vertical Proportion Curves. In addition, the statistical realization honours the wellbore observations of facies and the calibrated variogram model.
One last remark on this preliminary geostatistical modelling concerns the adoption of a single gaussian function although several geological processes, i.e. primary deposition and multiple subsequent diagenesis phases, determine the nature and distribution of lithotypes and would justify the use of at least two gaussian functions as suggested in the following. The reason for such a choice is the absence of sufficient and reliable data to constrain the spatial distribution of the different lithotypes. Despite such limitations in modelling, the availability of a significant amount of wellbore data and the attention awarded to the "Water" lithotype distribution made possible a realistic modelling of the conductive bodies.
The attribution of petrophysical properties to this model was performed using the statistical distributions of porosity and permeability of each lithotype. The parameters of those distributions were derived from petrophysical measurements on well-C1 core fragments. However, the permeability of the eighth lithotype, "Water", was given for each productive well directly from the pumping test results, taking into account the valid approximation that this facies mostly contributes to the well deliverability. Experimental vertical variograms.

Result: the Initial EHS Geostatistical Model
The resulting 3D model including the drilled Experimental area is shown in Fig. 23, with clear evidence of the drains, denoted as the "Water" lithotype, in blue colour. Two horizontal cross-cuts of this model at the Site scale (Figs. 24,25) show clearly the distribution of drains at the 53 and 85 m depths often observed as the major productive levels. West-East and North-South vertical cross-sections (Figs. 26,27) illustrate the vertical connectivity of lithotypes, in particular the conductive "Water" lithotype. Figures 28 and 29 show the same cross-sections in terms of permeability. Those permeability cross-sections actually reflect the distribution of the "Water" lithotype, that is the sole conductive lithotype.

Figure 24
Facies distribution at a depth of 53 m: the blue conductive facies forms a complex network.

Figure 25
Facies distribution at a depth of 85 m: the blue conductive facies forms a loose network.
According to this geostatistical model, the EHS aquifer appears as a limestone reservoir within which diagenetized (post-deposition) facies, i.e. the cherty limestone and fractured limestone lithotypes, are distributed as complex 3D volume elements that may or not be connected both horizontally and vertically. The conductive "Water" lithotype itself is found as isolated or connected 3D bodies of variable size within the previous diagenetized limestone volume elements. The horizontal connectivity of these conductive bodies looks hazardous, even at the two main productive depths around 50-55 m and 85 m. A vertical connection between the conductive channels found at these two depths may exist in some places but looks still much less effective than the areal connectivity.
Eventually, this distribution pattern of the conductive lithotype explains why neighbouring wells have similar or contrasted deliverabilities whether or not they cross conductive bodies and whether or not these bodies are connected away from wellbore.

Future Improvement of this Model
The wellbore images and production logs recently acquired in about half of the site wells provide a refined description of the lithologic column along wellbores, including a more precise location of the water entry points. Borehole images are particularly adapted to detect the presence of the partly-collapsed cave systems with clay fill-ins. In conjunction with those images, production logs reveal that the flow paths form a subdomain of the diagenetized limestone bodies without any apparent link with the original sedimentary units. That is, the assumption of a paleo-karst system, probably resulting from several episodes of burial and emersion, is being corroborated. These images and logs constitute valuable pieces of information to further constrain the geostatistical simulation of facies (or lithotypes) at wellbore locations. In addition, they will give the possibility to build a more realistic model by taking into account two mechanisms driving the distribution of facies instead of one in the preliminary model. These two mechanisms include respectively: -primary sediment deposition processes, that drive the distribution of all primary facies, essentially limestone, clayey limestone, marl and clay, -diagenetic processes that control the location of fractured and cherty limestones, of the flow paths ("water" facies), and also of dark clays that mainly constitute a filling sediment of fractures and caves.  Figure 26 West-East facies cross-section of the Site: a poor connectivity of the conductive facies is expected.

Figure 28
West-East permeability cross-section of the Site.

Figure 29
North-South permeability cross-section of the Site.

Figure 27
North-South facies cross-section of the Site: the conductive lithotype looks here more areally-continuous than in the West-East direction, but no vertical connection is observed between the 50 m and 85 m productive levels.
Separate geostatistical parameters will be determined for these two sets of facies, in particular the range of their respective variograms. A long variogram range will characterize the primary-deposition lithotypes, whereas a short range, probably less than the well spacing, will control the geostatistical simulation of embedded diagenetized bodies.
Eventually, the distribution of lithotypes will be simulated using the bi-gaussian distribution qualitatively schematized on the right-hand side of Figure 30, by comparison with the initial distribution shown on the left-hand side. This way, the diagenetized facies, including the conductive one, should clearly form nested bodies or channels within the primary layered deposits, thus yielding a realistic image of flowimpacting heterogeneities in the resulting model.
Nevertheless, the exact location and the detailed size, geometry and connectivity of these heterogeneities are expected to remain incompletely defined as their lateral extension is often less than the well spacing.

CONCLUSION -PERSPECTIVES
This hydrogeological Site constitutes a very instructive multi-disciplinary data base to set up consistent methodologies for characterizing the hardly-tractable heterogeneity of complex diagenetized reservoirs and for modelling their hydrodynamic behaviour.
The static characteristics and the dynamic behaviour of the aquifer were identified and cross-checked from data sources of various nature at different scales: -flow measurements at the Site scale and at the wellbore scale revealed the heterogeneous dynamic behaviour of this aquifer at all investigated scales; -preliminary seismic surveys indicated that meaningful markers above and inside the aquifer could be imaged over the Site, following their identification from well seismic or acoustic logging data; -core observations unveiled the existence of different diagenetic fingerprints strongly modifying the lithology, texture and properties of facies, thus explaining the complexity of most log responses; -the study of well-to-well correlations confirmed the variability of diagenetic fingerprints over short distances, less than the well spacing. The geostatistical model resulting from the integration of available data shows a distribution of facies that is consistent with the conceptual model of paleo-karst assumed for this aquifer.
Only the first step of data integration has been undertaken yet. The geostatistical model presented herein will soon be updated thanks to newly-available data, among which wellbore images and flow profiles.
Regarding new data acquisition, the detailed pattern of heterogeneities and flow paths will be further characterized thanks to dedicated seismic surveys and directional drillings that will explore the inter-well space. Next, the model resulting from the integration of this extensive database will be used as a support for the simulation of well pumping tests and of pressure interferences within the Site. This subsequent stage will be determinant to constrain the remaining uncertainties of the aquifer model that concern in particular the connectivity of conductive bodies or channels, and maybe also the role of the matrix background as a fluid-feeding source for the long-term deliverability of wells. To that end, history matching techniques that preserve the ascertained geological features will be applied. At a more advanced stage, once the connected flow network has been identified from well dynamic data, the realization of tracer tests is considered to further investigate this network in terms of flow and transport capabilities.