# Thermo-hydrodynamic modelling in silicoclastic reservoirs: case study of the Albian geothermal reservoir at Saclay, Paris Basin, France

Codjo T.F. Essou, Benjamin Brigaud, Miklos Antics, Pierre Ungemach, Perrine Mas, Remy Deschamp, Eric Lasseur, Yara El Bayssari, Gillian Bethune

<figure><img src="image_1.png" /></figure>

## To cite this version:

Codjo T.F. Essou, Benjamin Brigaud, Miklos Antics, Pierre Ungemach, Perrine Mas, et al. Thermo-hydrodynamic modelling in silicoclastic reservoirs: case study of the Albian geothermal reservoir at Saclay, Paris Basin, France. 48th Workshop on Geothermal Reservoir Engineering, Université de Stanford, Feb 2023, Stanford, United States. [hal-05297690](https://hal.science/hal-05297690v1)

HAL Id: hal-05297690

https://hal.science/hal-05297690v1

Submitted on 4 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Thermo-hydrodynamic modelling in silicoclastic reservoirs: case study of the Albian geothermal reservoir at Saclay, Paris Basin, France

Codjo T.F. Essou<fnref n="1" />,², Benjamin Brigaud<fnref n="1" />, Miklos Antics², Pierre Ungemach², Perrine Mas<fnref n="1" />, Remy Deschamp<fnref n="3" />, Eric Lasseur⁴, Yara El Bayssari<fnref n="1" /> and Gillian Bethune²,<fnref n="3" />

codjo-thomas-florent.essou@universite-paris-saclay.fr

**Keywords:** Clastic, Geothermal, Reservoirs, Models, Flow simulations

## ABSTRACT

Conceptual geological models aim at producing a coherent image of the investigated porous media. Such reservoir models, generally calibrated on the production data histories, facilitate predictions addressing the development of the geothermal resources. However, in order to reduce uncertainties and to improve prediction of interference between geothermal wells or early thermal breakthrough using numerical flow simulators, a custom designed approach is suggested. In the present study, the geological models were based on careful examination of historical core descriptions and well-log analysis using PETREL to obtain 3D models. These geomodels are used to simulate the mass and heat transfers via TOUGH3, ECLIPSE300 and PUMAFLOW softwares. Several calibration simulations of the temperature, pressure and flow patterns were performed, based on the last three years (2019-2022) geothermal production histories of the Saclay geothermal development site, located 20 km southwest of Paris. We show that using a high-resolution 3D grid simulation workflow constrained by sedimentary facies, it was possible to operate the simulations from different software now available. These software allow us to solve the flow and heat transfers in a structurally complex and heterogeneous multilayered geothermal reservoirs. We also compare simulation of water drawdown on different grids. TOUGH3, ECLIPSE300 and PUMAFLOW codes are used to predict the future flow and temperature evolutions, suggesting that all codes are adapted to simulate flow and thermal breakthrough. Ultimately, we were able to predict the preferential flow paths related to the heterogeneity of the targeted Albian siliciclastic reservoir in the Paris Basin. Preferential paths are recognized in the upper part of the reservoir (clean shoreface sand) and locally at the base of the reservoir where coarse sand facies are present. In the future, we recommend to produce only clean shoreface sands (*Sables de Frécambault Formation*) present in the upper part of the reservoir and the coarse channel sands (*Sables Verts*) at the base of the reservoir present locally in some wells. The rest of the reservoir contains too many clays and need to be isolated from production, thus limiting the risk of well plugging during reinjection.

## 1. INTRODUCTION

Since the Paris climate agreement in 2015, more and more countries are interested in developing geothermal energy for industrial and domestic uses. In France, the enthalpy potential of the Paris Basin has already been proven by numerous geothermal resource development projects, in particular in the deep geothermal aquifers of the Dogger Limestones (1500 m deep) and the Albian Sands (ca 600 m deep). The Albian sands and clayey deposits constitute a strategic reservoir throughout the Paris Basin for drinking water, since it is present throughout the greater Paris metropole (12 millions of inhabitants, Vernoux, 1997).

In a context where energy sobriety is sought, the Public Institution of Territory Development, called *Etablissement Public Aménagement (EPA) Paris-Saclay* in the vicinity of the Saclay high (*Plateau de Saclay*), 20 km south of Paris, has identified the development of low-carbon energy as a priority. *EPA Paris-Saclay* decided to build a heating and cooling network supplied by two geothermal doublets (four wells) producing the Albian aquifer. The initial heat production of each doublet was estimated at 4.9 MW. Each planned doublet should supply a heat plant and contribute to 60% of the renewable energy objective in this area.

Testing of these four geothermal wells shows an excellent productivity index of approximately 150m<fnref n="3" />/h/bar and almost zero skin. However, after a few weeks of production, the injectivity dropped drastically, leading to a total shut down of geothermal production at the two heat plants. Several hypotheses were debated among which the formation of internal cake induced or entrained suspended particles of micrometer size in the injection wells. The heterogeneity of the aquifer, alternating clays and sands, is responsible of the possible reinjection of fine particles (clays) in the geothermal wells. Actually, fine particles would pass through the filtration system downstream of the heat exchanger. It is therefore necessary to better characterize these reservoirs and improve numerical simulations and modelling in highly heterogeneous geothermal reservoirs (very rapid lateral variations between sands and clays). The objective is to compare different configurations of the reservoir, by proposing 5 different representations from very simple representation to more complex architecture depending on depositional facies. We also compare the results of hydro thermodynamic evolution in the reservoir using three codes (TOUGH3, PUMAFLOW and ECLIPSE300). This work is a prerequisite for any improvements in (hydro-dynamic and thermal) flow prediction within these heterogenous sand-clay reservoirs.

<footnote n="1">¹ Université Paris-Saclay CNRS, GEOPS, 91405, Orsay, France ; ² GEOFLUID, 165 Rue de la belle étoile, 95700 Roissy CDG, France ;</footnote>

<footnote n="3">³ IFP-Energie Nouvelle, 1-4 Avenue de Bois-Préau, 92852 Rueil-Malmaison, France and ⁴ BRGM, 3 avenue Claude Guillemin, BP 36009, 45060 Orléans, France</footnote>
## 2. GLOBAL FRAMEWORK

### 2.1 Geological Setting

The Paris Basin is an intracratonic sedimentary basin that extends over 110,000 km² in the northern half of France (Figure 1). This basin is bounded at the west by the Armorican Massif, the Vosges Massif at the east, by the Central Massif at the south and by the Ardennes Paleozoic Massif at the north-east. The Aptian and Albian periods are characterized by marine silicoclastic deposits resulting from the erosion of Paris basin paleo-borders (i.e. mainly from Armorican, Central and Ardennes Massifs). These deposits are mainly composed of terrigenous detrital sands and clays. The sedimentary deposits were deposited in shallow marine environments from the end of Aptian to Late Albian (Jacquin et al., 1998). The Albian sedimentary deposits are composed of seven Formations: (1) the *Sables Verts*, (2) *Argiles de l'Armance*, (3) *Sables des Drillons*, (4) *Argiles des Drillons*, (5) *Sables de Frécambault* topped by (6) the *Argiles de Gault* and finally by (7) the *Marnes de Brienne*. These sedimentary Formations are well-described in the study area by Lemoine et al. (1939) in the Orsay core-drill (Figure 2). The *Sables Verts*, *Sables des Drillons* and *Sables de Frécambault* consist of sands and form 3 reservoir units (Figure 2). Three clayey units are interlayered between these 3 sand reservoir layers.

<figure><img src="image_3.png" /><figcaption>Figure 1: A. Geological map and sedimentary section of the Paris Basin show the location (black stars) of the Orsay wells investigated in this study. OR—Orsay well (CR12). B. The west-east geological section of the Paris Basin modified after Gély and Hanot (2014)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Sedimentary section described adapted from Lemoine et al. (1939) in Orsay well.</figcaption></figure>
## 2.2 Review of the Paris Saclay geothermal reservoir

The geothermal field has been developed in an area called *Plateau de Saclay* (Saclay high), where two major universities are located (1) the *University Paris-Saclay* and (2) the *Institut Polytechnique de Paris*, at 20 km southwest of Paris. Two doublets have been completed each including 2 production and 2 injection wells. These two doublets, targeted the Albian sands, were drilled as open hole in 2017-2018. (Figure 3). Two heatplants were built (1) the Moulon heatplant and (2) the Ecole Polytechnique heatplant. Production started in 2019 but in early 2020, decrease of injectivity has been noticed. All wells were shut down in 2020. In 2020-2021, further to cleaning well test operation in the 2 injection wells Moulon-2 and Ecole-Polytechnique-2 have been made. In 2021, the 2 doublets are re-started

<figure><img src="image_5.png" /><figcaption>Figure 3: A. Overall; B. Paris Saclay geothermal well locations. Production wells are conventionally coded in red and injection wells in blue</figcaption></figure>

<table><caption>Table 1: Table 1: Reservoir hydraulic, hydrologic, petrophysical and thermal parameters. *PLT : Production Logging Tool (Flow/Temperature profiling)</caption><thead><tr><th>Parameter's</th><th>Moulon-1 Well</th><th>Moulon-2 Well</th><th>Ecole Polytechnique-1 Well</th><th>Ecole Polytechnique-2 Well</th></tr></thead><tbody><tr><td>Hydrostatic level on November 2017 in meter (m) from sea level</td><td>49.72</td><td>51.94</td><td>47.87</td><td>49.89</td></tr><tr><td>Temperature in °C</td><td>31.5</td><td>32.8</td><td>32.6</td><td>33.9</td></tr><tr><td>Hydraulic transmissivity (kh) in m²/s</td><td>4.3*10<sup>-3</sup></td><td>5.2*10<sup>-3</sup></td><td>4.8*10<sup>-3</sup></td><td>8.6*10<sup>-3</sup></td></tr><tr><td>Intrinsic transmissivity in D.m</td><td>330</td><td>311</td><td>364</td><td>626</td></tr><tr><td>Net cumulative producing height (h) in m</td><td>14m for 12 producing layer (Interpreted from PLT* log)</td><td>17.5 m for 7 producing layer (Interpreted from PLT* log)</td><td>18 m for 7 producing layer (Interpreted from PLT* log)</td><td>22.1 m for 10 producing layer (Interpreted from PLT* log)</td></tr><tr><td>Intrinsic permeability (k) in Darcy</td><td>23.6</td><td>17.8</td><td>20.2</td><td>28.3</td></tr><tr><td>Effective porosity (%)</td><td>39</td><td>34.</td><td>33.4</td><td>36</td></tr><tr><td>Wellbore storage coefficient (S)</td><td>1.6*10<sup>-4</sup><S<1.9*10<sup>-4</sup></td><td>1.6*10<sup>-4</sup><S<1.9*10<sup>-4</sup></td><td>1.1*10<sup>-4</sup><S<3.3*10<sup>-4</sup></td><td>1.1*10<sup>-4</sup><S<3.3*10<sup>-4</sup></td></tr></tbody></table>
## 3. MODELLING METHODOLOGY

To benchmark the best fit model, five reservoir simulations were run using PETREL (Schlumberger)/SKUAGOCAD (Aspentech Emerson) geomodelling software. Simulations compared to wells-interference data were made in order to select the accurate reservoir representation used to benchmark ECLIPSE300 (Schlumberger), THOUGH3 (University of Berkeley-LBNL) and PUMAFLOW (Beicip Franlab) softwares.

### 3.1 Reservoir models

The 2D grid boundary is squared as 10km*10km. It has been obtained by logarithmic discretization (Figure 4). Cell dimensions around the wells stand at 25m*25m, layering with 46 vertical layers. Thus, the total cells in the 3D grid amount to 419520. The caprock and bedrock are physically set to maintain constant pressure and temperature during simulation while using TOUGH3.

<figure><img src="image_6.png" /><figcaption>Figure 4: A. 2D grid → ON NE VOIT PAS LE NOM DES PUITS. B. 3D grid around the Albian conceptual model.</figcaption></figure>

#### 3.1.1 Model Case 1

In this simple case it is assumed that whole the reservoir thickness, porosity and permeability correspond to the average of the four wells data listed in Table 1. Hence, the thickness of this simple model stand at 17.9 m (average thickness in the four wells) filled with 35% porosity (average porosity in the four wells) and 22.7 Darcy permeability (average thickness in the four wells), (Figure 5A).

#### 3.1.2 Model Case 2

This conceptual, so called “sandwich model”, is derived from the work of Antics et al., (2011) who introduced the model concept described below in (1), (2) and (3). Here the methodology consists of extracting the total productive thickness, as well as the cumulative non-productive thickness to achieve the following reservoir structure:

(1) a first layer of thickness equal to half of the total productive thickness,

(2) an impermeable layer separating the reservoir in two parts,

(3) a second productive layer symmetrical to the first.

Petrophysical properties are averaged per well and propagated throughout the model. Thus each reservoir layer has a thickness equal to half of 17.9 m (i.e. 8.95 m), 35% porosity and 11.385 Darcy permeability (Figure 5B).

#### 3.1.3 Model Case 3

Here the model consists of three distinct reservoir layers from bottom to top (1) the Sables Verts, (2) the Sables de Drillons and (3) the Sables de Frecambault, separated by two hydraulically impervious barriers. The properties listed in Table 1 are averaged and are distributed according to the weight of each layer thickness (Figure 5C).

#### 3.1.4 Model Case 4

The reservoir was constructed using each well data issued from Production Logging Tool (flowmetry).

For all these wells, permeabilities of each productive layer have been derived from individual transmissivities using the following formula:

$$K_i h_i = K_h \times Q_i / Q,$$
where Kᵢ is the permeability of each productive layer individualised from flowmetry (in m²), hi is the thickness of each productive layer, Kh the total reservoir transmissivity in the well (in m²/s), Qᵢ the flow of each productive layer (m) and Q the well total flow.

A variogram is used to estimate the permeability trend recalculated in each productive layer and in each well (Figure 5D). Random connection is simulated and the grid is populated with the average porosities listed in Table 1.

<figure><img src="image_7.png" /><figcaption>Figure 5: Pseudo 3D views of reservoir models. A: Case 1; B: Case 2; C: Case 3. D: Case 4</figcaption></figure>

### 3.1.5 Model Case 5

The model includes well correlation, 2nd and 3rd depositional sequences and facies interpretation. Depositional environment have been identified in each well, then “upscaled” using the “Most of” algorithm. Facies were propagated in the 3D grid using the “Truncated Gaussian with trends” an algorithm which allows us to strongly control the distribution of facies, zone wise.

The neutron porosity log “NPHI” has been processed to obtain the effective porosity by shifting the clay impact from the recorded signal with the following relationship $Phie=NPHI(1-V_{shale})$. The obtained Phie logs were further “upscaled” by the arithmetic method, then propagated through the grid by the “Gaussian Random Function Simulation” algorithm and constrained by the interpreted depositional environments.

To obtain permeability, we used classical porosity-permeability relationships available in the literature, in our case used the relationship given by Zinszner et al. (2007) and Al Saadi et al. (2017):

$$K = 10^{[6,5+4,6*Log(Phie)]},$$

where K is the permeability in mD and Phie the corrected porosity obtained by NPHI well-logs from the 4 wells. The estimated permeability logs in the wells are further “upscaled” as geometric mean. Subsequently, permeability values were assigned to the model cells using the “Gaussian Random Function Simulation”. The variogram is a relatively layered association displaying equal and constant anisotropies. The distribution of the output data range is derived as absolute, and the distribution in each depositional facies defined from the upscaled permeability logs. Ultimately, a control of the 3D distribution of the constrained permeability, with respect to the three-dimensional (3D) distribution of the porosity was carried out by co-kriging in order to benefit from a supplementary constraint in the permeability distribution/process.
**3.2. Simulation assumptions:**

During simulations we considered the following assumptions:

a. Both thermal conduction and convection are simulated in the producing levels and the of mass and heat transfers in porous media are time dependent (non-steady state).

b. Equations of state consider

i. density effects (temperature-dependent density))

ii. viscosity effects (viscosity depends on the reinjected fluid temperature, inducing coupled hydraulic and thermal transfers).

**4. RESULTS AND INTERPRETATION**

**4.1 Sedimentary architecture: Model case 5**

<figure><img src="image_8.png" /></figure>

**Figure 6:** Wells cross-section and depositional facies portraying the architecture of the Aptian-Albian reservoir in the study area. In all wells, columns represent (1) the depth (m from Sea Surface Through Vertical Depth), (2) Gamma Ray (GR) log, (3) facies interpretation, (3) Production Logging Tool (PLT) and (4) the facies model. The second column in the Orsay well represents the 2^nd^ and 3^rd^ depositional sequences.

Facies are mainly based from Gamma Ray responses, which are very low in shoreface sands (Sables de Frecambault), suggesting the high quality of this reservoir (Figure 6 and 7). The Sables des Drillons Formation consists mainly of heterolithic sand bars facies and clays from tidal flat, whose contribution in the production of geothermal water is relatively low (as shown by PLT in Figure 6). The Sables Verts Formation consists of coarse sands and gravels, sands from tidal flat and heterolithic sand bars facies, and has locally good reservoir quality especially in its coarse sands and gravels facies (Figure 6 and 7). The clays from the lower offshore partition the reservoir, from the permeability barriers.
<figure><img src="image_9.png" /><figcaption>Figure 7: 3D view of reservoir architecture in case 5 geomodel</figcaption></figure>

## 4.2 Porosities and permeabilities and temperature

The different representations of porosities and permeabilities in two geomodels are expressed in **Figure 8**. The reservoir shows high heterogeneity in vertical and horizontal directions. In clean shoreface sands of the *Sables de Frécambault* Formation, porosity and permeabilities value stand at 35% and 10 Darcy respectively, and the temperature from the top to the bottom of the reservoir varies between 32 and 33°C.

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8: 3D views of porosity models in case 5; B: Case 5; C: Case 4; D: Case 4; E. Temperature model</figcaption></figure>

4.3 Interference tests

The wells of the Saclay plateau were evaluated prior to the geothermal production. Prior to the geothermal production, the production wells have been submitted to pumping tests in order to detect interferences and inter-doublet connectivities. At the western part of the geothermal field, the Moulon-2 was flowing at a stabilized 200m³/h water flow rate. A pressure gauge was placed at bottomhole of well Moulon-1 (initially shut down) to measure the drawdown induced by the producing well Moulon-2. Similarly at the eastern part of the geothermal field, the Ecole-Polytechnique-2 well was produced at a stabilized 205m³/h water flow rate. A gauge was placed at the bottomhole of the Ecole-Polytechnique-1 well (initially shut down) to measure the drawdown induced by the producing well Ecole-Polytechnique-2.

**Figure 9** emphasizes the interference impact, predicted by different conceptual models without any additional mathematical adjustment whatsoever of the history calibration matching. Models 1, 2 and 5 reflect the pace of the observed drawdown. The history matching curve of model 5 stands close to observed data. The conceptual model 5, whose petrophysical properties have been estimated via a single probabilistic approach, seems to present the best representation to predict the observed drawdown.

The modelling approach to obtain the conceptual model case 5 is a best practice to improve the simulation of history and prediction of geothermal reservoir performance. Even more times consuming, this high-resolution modelling approach should contribute to reduce the risks of early thermal breakthroughs, improve predictive location of reservoirs and de-risk future operations planned in the Paris Basin Albian reservoir.

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /><figcaption>Figure 9: Results of the simulated interference tests at initial state. A. Water drawdown at Ecole-Polytechnique-1; B. Drawdown at Moulon-1</figcaption></figure>

## 4.4 Simulation for prediction of thermal breakthrough with ECLIPSE300, PUMAFLOW and TOUGH3

The benchmark between ECLIPSE300 and TOUGH2 was exercised by Pham et al., (2019). For our study, we coupled thermo hydrodynamic simulations on conceptual model case 5 using ECLIPSE300, TOUGH3 AND PUMAFLOW. The production control rate of the geothermal water stands at 200 m3/h (representing the initial average rate of each doublet of the Paris Saclay project) and the injection temperature at 10° C over 30 years.

The cooling kinetics plotted in Figure 10 indicate a 0.2°C moderate drop at production well Moulon-1. The Orsay well, producing geothermal waters for heating the city swimming pool, is not impacted by any thermal breakthrough (Figure 10). At the above mentioned production conditions, the Ecole-Polytechnique-1 production well would have reached a thermal breakthrough of 0.9°C after 30 years (Figure 10).

The three simulators ECLIPSE300, PUMAFLOW and TOUGH3 yield similar results. However, computing times vary significantly between simulators depending on the internal solver architecture. This parameter will determine the performance of the different tools.

The current work highlights the benefit of the state-of-the-art flow simulator used in the oil and gas industry for the development of geothermal energy, especially when it addresses deep and very deep low and high enthalpy geothermal targets.

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /><figcaption>Figure 10: Cooling kinetics; A. Benchmark between ECLIPSE300, PUMAFLOW and TOUGH3 software; B. Reservoir temperature final state at year 2052</figcaption></figure>

## 6. CONCLUSION

The conceptual geological models of cases 1, 2 3 and 4 failed to predict water drawdown during interference tests. This study shows that upon data availability, the modeling workflow should integrate the sedimentary facies in order to better express the vertical meter-scale heterogeneity of the reservoir. PETREL, ECLIPSE300, PUMAFLOW and TOUGH3 have been successfully tested. These tools make it possible to accurately simulate the thermo-hydrodynamic transfer using complex grids architectures. The predicted thermal breakthrough is very limited, which is very promising regarding the sustainable use of the Albian deep geothermal aquifer in supplying heat to the Saclay plateau. Preferential paths were identified in the upper part of the reservoir (clean shoreface sand) and, locally, at the base of the reservoir where coarse sand facies are present. To limit the well plugging during reinjection, it is recommended to produce only the clean sands (Sables de Frécambault Formation) present in the upper part of the reservoir (20m thick) and the coarse channel sands (Sables Verts) at the base of the reservoir present locally in the western part of the geothermal field.

## 7. ACKNOWLEDGMENTS

This work is issued from a collaborative project (No 2021-2034) between the University Paris-Saclay and GEOFLUID. This work is supported by the French Agency for Technological Research Grant n°2021/0773 and GEOFLUID for co-funding this doctorate project. The UPGEO "UPscaling and heat simulations for improving the efficiency of deep GEOthermal energy" program, funded by the French Agency for Research [ANR-19-CE05-0032-01], also supported this project. University Paris-Saclay benefits from the Schlumberger and the Beicip-Franlab Software Donation Programs. We are grateful to Schlumberger and Beicip-Franlab for the academic licenses for PETREL/ECLIPSE300 and PUMAFLOW respectively, and for their valuable support during their utilization.

## REFERENCES

* Al Saadi, F., Wolf K., Kruijsdijk C.V., (2017). Characterization of Fontainebleau Sandstone: Quartz Overgrowth and its Impact on Pore-Throat Framework. Journal of Petroleum & Environmental Biotechnology. 7, 328

* Antics M., Ungemach P., (2011). Hydraulic and Thermal Modelling practice for Geothermal Exploitation of the Dogger in the Paris Region. Rapport RP-59591-FR-Alforville. GEOFLUID, Roissy, pp. 105–106.

* Gély J.-P., Hanot F., (2014). Coupe Géologique du Bassin parisien et du Fossé hénan, in Le Bassin parisien, un Nouveau Regard sur la Géologie, Paris, Association des Géologues du Bassin de Paris, Bull. Inf. Géol. Bass. Paris, Mémoire hors-série, v. 9, p. 1 plate.

* Jacquin T., Dardeau G., Durlet C., de Graciansky C., Hantzpergue P., (1998). The North Sea Cycle: an Overview of 2nd-order Transgressive/Regressive Facies Cycles in Western Europe, in: de Graciansky, P.-C., Hardenbol, J., Jacquin, T., Vail, P.R. (Eds.), Mesozoic and Cenozoic Sequence Stratigraphy of European Basins. SEPM Special Publication, pp. 445–466.

* Lemoine P., Humery R., Soyer R., (1939). Les Forages Profonds du Bassin de Paris : la Nappe Artésienne des Sables Verts. Mémoires du Muséum National d'Histoire Naturelle, Nouvelle Série volume 11. Paris, 700 p.
Pham M., Sullera M., Williams M. J., Henneberger R., (2019). ECLIPSE Geothermal - A Next-Generation Geothermal Reservoir Simulator. European Geothermal Congress 2019. Den Haag, The Netherlands, 11-14 June 2019

Vernoux J.F., Maget P., Afzali H., Blanchin R., Donsimoni M., Vairon J., (1997). Synthèse Hydrogéologique du Crétacé Inférieur du Bassin de Paris. Rapport BRGM DSGR/IDF No. R39702. Bureau de Recherches Géologiques et Minières, Orléans. 93 p.

Zinszer B., Pellerin F., (2007). A Geoscientist's Guide to Petrophysics, Technip, IFP Publication, 384 p.