Hadrien Thomas

Benjamin Brigaud

Thomas Blaise

Elodie Zordan

Hermann Zeyen

Maxime Catinat

Simon Andrieu

Emmanuel Mouche

Marc Fleury

Upscaling of geological properties in a world-class carbonate geothermal system in France: from core scale to 3D regional reservoir dimensions

Keywords: limestone, Paris Basin, Jurassic, 3D geomodeling, reservoir, facies, geothermics

The greater Paris area has some 12 million inhabitants and 48 heating network production units that exploit the heat capacity of a 1.5 km deep aquifer, the Bathonian limestone. This is one of the most productive aquifers in the world for district heating, with an annual output of about 1.7 TWh of energy. The current challenge for Paris is to triple the number of heating networks using geothermal energy so as to reduce dependence on fossil fuels (40% in 2020 in France). As no detailed geothermal reservoir model is available, drilling involves geological risk. A recent well, drilled at Grigny (20 km south of Paris), turned out to be a failure in transmissibility because the bed of permeable limestones it tapped into was very thin (meterscale). The main aim of this study is to create a digital database and a 3D geological model of this aquifer to minimize geological risks and optimize the location of future geothermal operations around Paris. By compiling data from 168 wells, a high-resolution 3D geological model of 360 km 3 size is constructed (about 40 km x 50 km x 0.2 km), made up of 12.2 million cells and displaying sedimentary facies, sequence stratigraphy, porosity (F) and permeability (k). About 20% of the oolitic and bioclastic facies are of good reservoir quality (F > 13% and k > 350 mD), especially in two targeted, high-quality reservoir sequences. These facies of interest probably correspond to giant dunes and a shoal/barrier prograding from east to west.

In these facies, permeable zones are generally 4 m thick and form patches of 1600 m x 1100 m, on average, elongated perpendicular to the depositional slope. 2D and 3D maps of temperature, salinity, porosity, transmissivity, and permeability allow us to understand the areas of interest for geothermal exploration, as demonstrated around Grigny. This model helps us to apprehend better the heterogeneous character of the reservoir for geothermal prospection and to reduce the risk of future doublets during well implantation. Detailed local models may be extracted to anticipate better the implantation of new doublets in areas with already densely spaced existing wells.

Introduction

The use of renewable energies has progressed in France over the last ten years, from 250 TWh in 2010 to 322 TWh in 2020, and they now make up 23% of the country's total energy production [START_REF] Beck | Chiffres clés de l'énergie édition[END_REF][START_REF] Phan | Chiffres clés des énergies renouvelables Edition 2021[END_REF]. Despite large aquifers in sedimentary basins, deep geothermal heat accounts for only 0.6% of renewable energy production in France, with 2 TWh in 2020 [START_REF] Beck | Chiffres clés de l'énergie édition[END_REF][START_REF] Phan | Chiffres clés des énergies renouvelables Edition 2021[END_REF]. Heating networks largely depend on fossil fuels (40% in 2020), and deep geothermal energy represents less than 6% of the energy mix of heating networks in France [START_REF] Phan | Chiffres clés des énergies renouvelables Edition 2021[END_REF]. France's energy-climate strategy is to produce 5 TWh of geothermal heat by 2028, compared to 1 TWh in 2010 and 2 TWh in 2020 [START_REF] Beck | Chiffres clés de l'énergie édition[END_REF][START_REF] Phan | Chiffres clés des énergies renouvelables Edition 2021[END_REF]. The greater Paris area (Fig. 1A) has exceptional deep geothermal energy resources, exploited in two main aquifers: (1) the Middle Jurassic limestones between 1.5 and 2 km depth with a temperature range of 55-80 C and (2) the Lower Cretaceous sandstones, between 0.5 km and 0.8 km depth and at 30-38 C [START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF][START_REF] Bel | Etude préalable à l'élaboration du schéma de développement de la géothermie en Île-de-France[END_REF]. This geothermal energy is mainly used for district heating [START_REF] Lund | Direct utilization of geothermal energy 2020 worldwide review[END_REF]. Seventy-six geothermal doublets (152 wells), nine triplets (27 wells) and six single geothermal wells have been drilled in the greater Paris area since 1931 (Orsay well). In 2022, 43 of these doublets, 9 of these triplets and 3 of the single wells are still active in these two main reservoirs. The greater Paris region is home to the highest concentration of low-energy geothermal heating network production units worldwide [START_REF] Lund | Direct utilization of geothermal energy 2020 worldwide review[END_REF]. While a single conventional production network unit in the Middle Jurassic aquifer provides heating for 4000 to 6000 housing units, all the production units together provide direct collective heating for 6-7% of the 12 million inhabitants of the greater Paris area. With 48 geothermal heating units currently producing geothermal water, the Middle Jurassic carbonate is the main aquifer exploited in the greater Paris area. This aquifer is one of the most productive in the world, supplying 1.7 TWh of energy per year. As some areas of this aquifer are intensively exploited, the objective of tripling geothermal energy production implies targeting new areas around Paris. A 3D high-resolution and predictive geological model could allow better understanding the reservoir heterogeneity. By transferring knowledge, 3D models allow more doublets to be drilled efficiently in an area where the density of geothermal wells is significant. There is a real risk that an installation may not reach a geothermal resource with sufficiently high transmissivity and temperature to ensure the project is cost-effective during its lifetime. Sites with poor transmissivity are challenging to predict and constitute real obstacles to the future development of geothermal energy in the greater Paris area. In the late 1980s, pioneering work by the French Geological Survey (BRGM) proposed maps of facies, porosity, and transmissivity of the Middle Jurassic carbonate reservoir based on 110 wells drilled between 1976 and 1986 [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF]. These maps have been used to guide the development of new operations since the 1990s, but they do not reliably predict low transmissivity areas because they fail to integrate small-scale variations in facies, porosity, and permeability.

Although 671 wells have been drilled for petroleum exploration or production since the 1950s, and 186 wells have been drilled for geothermal production, there is currently no detailed 3D reservoir model on regional scale compiling all the well data of the Middle Jurassic limestones [START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF]. As no 3D heterogenous geological and no recent geological study of the reservoir are available, a recently drilled well at Grigny (well GGR-2, at the same location as GGR-1 Fig. 1) turned out to be a failure regarding water flow as the aquifer turned out to be composed of very thin beds of porous limestones (meter-scale) with too low transmissivity. Geological models would help to reduce the risk involved in future geothermal operations around the Paris region for new doublet implantations [START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF]. Hydro-dynamic simulations are currently based on highly simplified synthetic models, assuming one homogeneous productive layer or multiple permeable layers (2 to 7 layers) with homogeneous hydraulic conductivity (typically 10 -6 m.s -1 ) or homogenous permeability (typically 5 Darcys), separated by impermeable strata [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF][START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF][START_REF] Le Brun | Hydraulic and thermal impact modelling at the scale of the geothermal heating doublet in the Paris Basin, France[END_REF][START_REF] Ungemach | Subhorizontal well architecture and geosteering navigation enhance well performance and reservoir evaluation A field validation[END_REF]. However, in Middle Jurassic limestones, with an emphasis on the Bathonian limestones, permeabilities exhibit a high degree of spatial heterogeneity across all scales, depending largely on the depositional environment and the subsequent diagenetic processes [START_REF] Brigaud | Acoustic Properties of Ancient Shallow-Marine Carbonates: Effects of Depositional Environments and Diagenetic Processes (Middle Jurassic, Paris Basin, France)[END_REF](Brigaud et al., , 2014a)). To achieve the energy-climate objective set out for France and for the greater Paris area it is essential to de-risk future operations that are often supported by public investors. This de-risking can be done by better predicting variations in the 3D carbonate reservoir and the reservoir quality of the Middle Jurassic limestones. The main problem is finding a way to spatially interpolate the porosity and permeability measured on rock samples or estimated by well logging from the centimeter or micrometer scale to the scale of the entire reservoir (hundreds of meters to kilometers). The objective of this work is to propose a fine-scaled geological model (XYZ dimension of the entire model of about 40 km x 50 km x 0.2 km and XYZ cell dimension of 100 m x 100 m x 3.73 m on average), based on upscaled porosity and permeability data measured from mini-cores and well logs and revealing small-scale heterogeneities. This objective was achieved through the following steps:

1. data compilation (well log, porosity and permeability) on a volume of 360 km 3 through the creation of a new, reliable and exploitable database by digitizing old well-logs available in .pdf or .tiff to ASCII format; 2. facies heterogeneity analysis by description of cores and cuttings and determination of facies with good reservoir qualities; 3. construction of the reservoir architecture by conducting high-resolution sequence stratigraphy correlation (third-order sequences) and determining the sequence(s) presenting the best reservoir quality; 4. definition of a method for upscaling facies distribution applicable to the Middle Jurassic carbonate platform; 5. definition of a method for upscaling porosity and permeability from mini-core scale to regional reservoir proportions; 6. creation of temperature maps of the geothermal water in the reservoir, temperature gradient, salinity, and transmissivity to better map the geothermal characteristics; 7. construction and prediction of reservoir heterogeneity on the regional scale in a 3D grid composed of 12.2 million cells of X, Y, Z dimensions of 100 m x 100 m x 3.73 m size (on average for Z) using Petrel© geomodeler.

The main objective of this study is to construct a 3D static model of petrophysical properties.

The resulting 3D model is one possible 3D realization showing the geological heterogeneities and complex sedimentary architecture. Based on the proposed workflow and existing data, different realizations of static models may be produced allowing quantification of the risk of failure. Th resulting 3D model will in addition be a prerequisite for thermo-hydro-dynamic simulations and for optimizing the location of future geothermal operations around Paris.

Geological, geothermal and structural settings

This study is based on the well log data available for a specific area of the Paris Basin, between Paris and Melun (Fig. 1A-B), located in the center of the greater Paris area.

Structural, palaeogeographic and stratigraphic settings

The Paris Basin is an intracratonic sedimentary basin, formed by three significant tectonic phases: (1) extension with crustal thinning from the Permian to the end of the Jurassic, marking a period of intense subsidence; (2) a period of slow subsidence (Early Cretaceous-Eocene) with an increased rate during the Late Cretaceous; and (3) a period of uplift since the Neogene [START_REF] Brunet | Subsidence of the Paris basin[END_REF][START_REF] Guillocheau | Meso-Cenozoic geodynamic evolution of the Paris Basin: 3D stratigraphic constraints[END_REF]Brigaud et al., 2014a). The basement of the basin comprises Permian-Carboniferous sedimentary deposits and Palaeozoic granitic and metamorphic rocks. It stretches from east to west for about 600 km, from the edge of the Armorican massif in the west to the Vosges massif in the east [START_REF] Guillocheau | Meso-Cenozoic geodynamic evolution of the Paris Basin: 3D stratigraphic constraints[END_REF]. The Paris Basin extends from the English Channel to the Massif Central from north to south. The structural context of the Paris Basin center has recently been reappraised from seismic data [START_REF] Beccaletto | Overview of the subsurface structural pattern of the Paris Basin (France): Insights from the reprocessing and interpretation of regional seismic lines[END_REF][START_REF] Dentzer | Thermal and seismic hints for chimney type crossstratal fluid flow in onshore basins[END_REF]. These studies reveal faults within the study area, such as the Beynes-Meudon fault-anticline system and the Valpuiseaux fault (Fig. 2B). The Beynes-Meudon fault complex originated from the N-S compressive Pyrenean phase that peaked in the Paris Basin during the Late Eocene [START_REF] Bergerat | La fracturation nivernaise: Influences bourguignonne et centralienne sur la structuration du Nivernais[END_REF][START_REF] Guillocheau | Meso-Cenozoic geodynamic evolution of the Paris Basin: 3D stratigraphic constraints[END_REF][START_REF] Lenoir | Uranium incorporation in fluorite and exploration of U-Pb dating[END_REF]. The Valpuiseaux fault is rooted in the basement and forms the south-eastern boundary of the Hurepoix block inherited from the Gondwana-Avalonia collision [START_REF] Cazes | Étude de la croûte terrestre par sismique profonde: profil nord de la France; structure hercynienne[END_REF][START_REF] Autran | Gravity and Magnetic Expression of Terranes in France and Their Correlation Beneath Overstep Sequences[END_REF][START_REF] Averbuch | Is there a remnant Variscan subducted slab in the mantle beneath the Paris basin? Implications for the late Variscan lithospheric delamination process and the Paris basin formation[END_REF][START_REF] Briais | Response of a low subsiding intracratonic basin to long wavelength deformations: the Palaeocene-early Eocene period in the Paris basin[END_REF]. (Berthuis-1), LJU-1D (Longjumeau-1), GGR-1 (Geothermy Grigny-1), LPQ-1 (Le-Petit-Quincy-1), VLG-6D (Vert-le-Grand-6D), , , VSD-1 (Vert-Saint-Denis-1), ETY-1D (Etrechy-1D), BR119 . The transect illustrated in figure 9 (green line from SW to NE) is located between ETY-1D in the SW and BTS001 in the NW. The transect illustrated in figure 10 (green line from NW to SE) is located between GCDN-1 in the NW and BRIE-119 in the SE.

Names in blue boxes are French administrative departments, their boundaries being indicated by black lines.

During the Middle Jurassic, the study area was a vast epicontinental sea located at subtropical latitudes of about 30° to 25°N. This marine domain was open to the Atlantic, Tethys, and Northern seas [START_REF] Contini | Evolution paléogéographique de la France au Jurassique moyen[END_REF][START_REF] Thierry | Middle Callovian, map 9[END_REF]. Deposits from the upper Bajocian to lower Callovian consist of alternating marls and carbonates. From bottom to top, the five sedimentary formations can be described as follows (Fig. 2B, [START_REF] Mégnien | Synthèse géologique du bassin de Paris[END_REF][START_REF] Brosse | Selection and Characterization of Geological Sites able to Host a Pilot-Scale CO2 Storage in the Paris Basin (GéoCarbone -PICOREF)[END_REF][START_REF] Delmas | Petrophysical Properties of the Middle Jurassic Carbonates in the PICOREF Sector (South Champagne, Paris Basin, France)[END_REF]:

(1) The Marnes à Ostrea acuminata consist of marls ranging from a few meters to nearly 50 m in thickness, deposited in lower to upper offshore environments and dated from the late Bajocian to the mid-Bathonian (zigzag ammonite zone).

(2) The Calcaires marneux à Pholadomyes, dated from the mid-Bathonian, are mainly composed of marls deposited in upper offshore environments. Ooid limestone patches are present locally. The thickness of this formation varies from 30 to 60 m.

(3) The Oolithe Blanche Formation consists of ooid limestones deposited on a high-energy shallow platform [START_REF] Vincent | Giant subaqueous carbonate dunes: a revised interpretation of large-scale oo-bioclastic clinoforms in the middle Jurassic of the Paris Basin and its implications[END_REF]. This formation is the main reservoir for geothermal energy in the Middle Jurassic limestones [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF], Lopez et al., 2010). The thickness of this formation varies from 50 to 100 m.

(4) The Calcaire de Comblanchien consists of micritic limestones and oncoidal-bioclastic limestones. This formation migrates from the south-eastern to the north-western part of the basin [START_REF] Gaumet | Middle Jurassic production rates and "patchy" architecture of the carbonate systems along the north western Thetysian margin (Paris basin to Subalpine basin)[END_REF] (Fig. 2B). Small production levels are found locally in this formation. This formation is absent west of Paris (Fig. 2A,B).

(5) The top of the Middle Jurassic limestones is composed of a thin (from a few meters to 20 m) limestone level of the early Callovian age called the Dalle Nacrée Formation. It is a wellknown oil reservoir in the Paris Basin [START_REF] Javaux | La plate-forme parisienne et bourguignonne au Bathonien terminal et au Callovien: dynamique sedimentaire, sequentielle et diagenetique: place et creation des reservoirs potentiels[END_REF][START_REF] Perrodon | M 51: Interior Cratonic Basins[END_REF]. Boreholes are cemented down to the base of the Dalle Nacrée Formation [START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF], to avoid oil contamination of geothermal fluids.

The total thickness of the upper Bajocian to lower Callovian formations is about 200 m in the study zone. [START_REF] Purser | Subdivision et interprétation des séquences carbonatées[END_REF] recognized three "lithoclines" in the Middle Jurassic limestones of the eastern Paris Basin, corresponding to biosedimentary packages and forming three regressive sequences. These three sequences form a progressive upward transition from clays and marls to shallow marine limestones, topped by a hardground surface [START_REF] Purser | Subdivision et interprétation des séquences carbonatées[END_REF]. The Bathonian lithocline is well-recognized in the study area from the Marnes à Ostrea acuminata at the base to the Calcaire de Comblanchien formation at the top. The cap rock consists of clayey and marly sedimentary formations dated from the early/mid Callovian to the early Oxfordian [START_REF] Vail | La stratigraphie sequentielle et son application aux correlations chronostratigraphiques dans le Jurassique du bassin de Paris[END_REF][START_REF] Perrodon | M 51: Interior Cratonic Basins[END_REF]. In this cap rock, a well-defined meterthick iron ooid layer appears in all the well-logs of the Paris Basin [START_REF] Vail | La stratigraphie sequentielle et son application aux correlations chronostratigraphiques dans le Jurassique du bassin de Paris[END_REF]. It is dated by ammonite fauna from the late Callovian (athleta and lamberti zones) to the early Oxfordian (cordatum zone) [START_REF] Vail | La stratigraphie sequentielle et son application aux correlations chronostratigraphiques dans le Jurassique du bassin de Paris[END_REF][START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF]Collin et al., 2005, Figs 2A, 3).

The Jurassic limestones of the Paris Basin were selected by [START_REF] Vail | La stratigraphie sequentielle et son application aux correlations chronostratigraphiques dans le Jurassique du bassin de Paris[END_REF] for the first application of the seismic stratigraphic model [START_REF] Vail | Seismic stratigraphy and global changes of sea level: Part 4. Global cycles of relative changes of sea level : Section 2. Application of seismic reflection configuration to stratigraphic interpretation[END_REF] in an intracratonic basin.

Following this first application, Purser's "lithoclines" were interpreted as second-order stratigraphic sequences [START_REF] Jacquin | Séquences de dépôt et cycles régressif/transgressifs en domaine marin carbonaté : exemple du Dogger du Bassin de Paris; Depositional sequences and regressive/transgressive facies cycles in marine carbonate settings : example of the Dogger of the Paris Basin[END_REF][START_REF] Jacquin | The North Sea cycle: an overview of 2nd-order transgressive/regressive facies cycles in western Europe[END_REF][START_REF] Vincent | Giant subaqueous carbonate dunes: a revised interpretation of large-scale oo-bioclastic clinoforms in the middle Jurassic of the Paris Basin and its implications[END_REF]. In the 1990s, biostratigraphic work resulted in cross-referencing the ammonite zonation with brachiopod zonations, particularly in sections located in the distal sedimentation area (in the Nièvre department) [START_REF] Garcia | Les associations de brachiopodes du Jurassique moyen du bassin de Paris; une echelle biochronologique ponctuee de niveaux-reperes pour la contrainte des correlations sequentielles a haute resolution[END_REF]. This work combines the brachiopod association zonation with the ammonite zonation, increasing the resolution [START_REF] Garcia | Les associations de brachiopodes du Jurassique moyen du bassin de Paris; une echelle biochronologique ponctuee de niveaux-reperes pour la contrainte des correlations sequentielles a haute resolution[END_REF][START_REF] Garcia | The validity of two biostratigraphic approaches in sequence stratigraphic correlations: brachiopod zones and marker-beds in the Jurassic[END_REF]. This zonation enabled the Bathonian and Callovian carbonates to be dated accurately in oil exploration boreholes in the Paris Basin [START_REF] Garcia | Les associations de brachiopodes du Jurassique moyen du bassin de Paris; une echelle biochronologique ponctuee de niveaux-reperes pour la contrainte des correlations sequentielles a haute resolution[END_REF][START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF].

On the carbonate platform, the brachiopod-rich levels correspond to maximum flooding surfaces [START_REF] Garcia | Les associations de brachiopodes du Jurassique moyen du bassin de Paris; une echelle biochronologique ponctuee de niveaux-reperes pour la contrainte des correlations sequentielles a haute resolution[END_REF][START_REF] Garcia | The validity of two biostratigraphic approaches in sequence stratigraphic correlations: brachiopod zones and marker-beds in the Jurassic[END_REF] detailing a third-order sequence stratigraphy framework, described in the Paris Basin boreholes by [START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF]. From the Bathonian base to the middle Callovian end, 14 third-order sequences have been recognized [START_REF] Garcia | Les associations de brachiopodes du Jurassique moyen du bassin de Paris; une echelle biochronologique ponctuee de niveaux-reperes pour la contrainte des correlations sequentielles a haute resolution[END_REF][START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF][START_REF] Jacquin | The North Sea cycle: an overview of 2nd-order transgressive/regressive facies cycles in western Europe[END_REF]. The maximum flooding surfaces of each sequence exhibit a specific brachiopod association, enabling us to date these sequences [START_REF] Garcia | Les associations de brachiopodes du Jurassique moyen du bassin de Paris; une echelle biochronologique ponctuee de niveaux-reperes pour la contrainte des correlations sequentielles a haute resolution[END_REF][START_REF] Thomas | Contribution of drone photogrammetry to 3D outcrop modeling of facies, porosity, and permeability heterogeneities in carbonate reservoirs (Paris Basin, Middle Jurassic)[END_REF][START_REF] Vincent | Giant subaqueous carbonate dunes: a revised interpretation of large-scale oo-bioclastic clinoforms in the middle Jurassic of the Paris Basin and its implications[END_REF] (Fig. 2B). [START_REF] Hendry | Geochemical trends and palaeohydrological significance of shallow burial calcite and ankerite cements in Middle Jurassic strata on the East Midlands Shelf (onshore UK)[END_REF] and [START_REF] Brigaud | The origin and timing of multiphase cementation in carbonates: Impact of regional scale geodynamic events on the Middle Jurassic Limestones diagenesis (Paris Basin, France)[END_REF]. The same color code for facies association types has been used in subsequent figures. The red rectangle in the figures marks the study area. B-Schematic lithostratigraphic section of the Middle Jurassic carbonates reconstructed from wells drilled in the Paris Basin modified from [START_REF] Nader | Diagenetic effects of compaction on reservoir properties: The case of early callovian "Dalle Nacrée" formation (Paris basin, France)[END_REF], flattened at the top to the ooid iron layer (Callovian/Oxfordian boundary).

Geothermal setting

The first geothermal development in the Middle Jurassic limestones of the Paris Basin was initiated in the late 1960s at Melun [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF]. The operation consisted of drilling two wells targeting the Bathonian limestones (Oolite Blanche and Calcaire de Comblanchien formations) from 1600 to 1700 m depth to form an open loop with one production well and one injection well, a technology known as a geothermal doublet [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF]Lopez et al., 2010). 78 doublets or triplets have produced or are producing geothermal water ranging from 55 to 80 °C, which is re-injected into the Middle Jurassic limestone aquifer at about 40 °C (Fig. 1B) [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF][START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF]. Of the 78 installations, 36 geothermal doublets and triplets lie within the study area and were drilled into or through the Middle Jurassic carbonate reservoir (Fig. 1B). Subsurface temperatures were compiled by [START_REF] Housse | Potentiel géothermique du bassin parisien (Rapport No. 74-7-0990[END_REF][START_REF] Le Nir | Ressources géothermiques du département de l'Essonne (91)[END_REF], [START_REF] Bonté | Subsurface temperature maps in French sedimentary basins: new data compilation and interpolation[END_REF] and [START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF]. The temperature of the aquifer is heterogeneous (Fig. 1B). The lowest temperatures are found in the northern part of the study area, below Paris and the Seine-Saint-Denis Department (Fig. 1B), reaching about 50 °C at a depth of 1650 m (thermal gradient about 23 °C/km). Temperatures in the study area are the highest south-east of Créteil (around 80 °C at a depth of 1660 m, thermal gradient of 40.9 °C/km) and south of Cachan (east of Fresnes, around 75 °C at a depth of 1590 m, thermal gradient of 39.6 °C/km). From Fresnes to the south of Paris, the temperature drops by 18 °C from 72 °C in Fresnes GFR1D well (1600 m below the surface, thermal gradient of 37.5 °C/km) to 54 °C south of the Paris boundary in the GBA1 well located at Bagneux (1510 m below the surface, thermal gradient of 27.8 °C/km, Fig. 1B). All the gradients are computed assuming a mean surface temperature of 12 °C. A possible explanation for temperatures dropping locally by nearly 20°C at the top of the reservoir may be the presence of pervasive fractures, leading to downward flow from the Cretaceous aquifers, but this explication is still debated in the community [START_REF] Dentzer | Thermal and seismic hints for chimney type crossstratal fluid flow in onshore basins[END_REF]. The temperature in the western part of the study zone is not well defined as no geothermal wells are available and historical exploration wells did not report reliable well temperatures [START_REF] Fabris | Les températures dans les forages -Gisements pétroliers du bassin Parisien et de l'Aquitaine[END_REF]. Additional information on groundwater flow and salinity can be found in the works of [START_REF] Fabris | Les températures dans les forages -Gisements pétroliers du bassin Parisien et de l'Aquitaine[END_REF], [START_REF] Aubertin | Determination experimentale de la vitesse naturelle d'ecoulement de la nappe geothermique du Dogger en region parisienne[END_REF], [START_REF] Wei | Regional modelling of groundwater flow and salt and environmental tracer transport in deep aquifers in the Paris Basin[END_REF], [START_REF] Menjoz | Flow of formation water in the Jurassic of the Paris Basin and its effects[END_REF] and [START_REF] Matray | Origin and age of deep waters of the Paris Basin[END_REF].

Material and methods

Cores, well-logs, petrophysical data, temperature, and salinity

The digital database was compiled from logs, reports and other technical data from 168 wells (94 oil exploration wells and 74 geothermal wells) available in the French Geological Survey subsurface database (BRGM, https://infoterre.brgm.fr/). The Sergines-1 and Berthuis-001 wells are considered in this study as reference wells but are not included in the final modeled zone (Fig. 1A). The well logs generally consisted of gamma-ray logs (GR, 151 logs), sonic logs (delta time, DT, 74 logs), resistivity logs (mainly laterolog deep, LLd and laterolog shallow, LLs, 62 logs), spontaneous potential logs (SP, 39 logs) and calliper logs (CAL, 5 logs). Porosity logs (neutron porosity, NPHI, neutron, NEUT, total combinable magnetic resonance, TCMR or related) were available for 76 wells, and density logs (Rhob) for 39 wells. An Nuclear Magnetic Resonance (NMR) logging tool (Combinable Magnetic Resonance -CMR©) was only available for the Cachan GCAH-2 well [START_REF] Ungemach | Subhorizontal well architecture and geosteering navigation enhance well performance and reservoir evaluation A field validation[END_REF][START_REF] Wielemaker | Delineating the geothermal structure and flow properties in a sub-horizontal well with the use of wireline and LWD data in a multiphysics approach[END_REF]. The NMR logging tool provides continuous porosity and permeability along a well [START_REF] Westphal | NMR Measurements in Carbonate Rocks: Problems and an Approach to a Solution[END_REF]. The permeability (kSDR) from the NMR tool is obtained with the Schlumberger Doll Research (SDR) default relationship [START_REF] Kenyon | A Three-Part Study of NMR Longitudinal Relaxation Properties of Water-Saturated Sandstones[END_REF]. Routine petrophysical analyses of porosity (F) and permeability (k) from mini-cores (plugs from cores) were compiled from wellreports, from 40 wells for F (4561 measurements in total) and 37 wells for k (3317 measurements in total). F-k data were tabulated in an Excel spreadsheet and converted into ".las" data files to be geolocated in 3D. 64 wells included water temperatures, and 66 wells reported salinity values in the aquifer. All the log data are available in the supplementary data sheet (Supplementary Table A). The well coordinates, trajectory, logs, temperature, salinity, and F-k data measured on mini-cores were imported into the Petrel© software for 3D processing. The well list is available in Supplementary Table A andF

-k data in Supplementary Table B.

Petrography and sequence stratigraphy

Vert-Saint-Denis 1 (VSD1) and Sergines 1 (SER1) well cores, stored at the BRGM core storage facility, were described. Thirty-four samples for these two core-wells as well as 123 cutting samples from nine other wells (Paris Aiport 1&2, Bagneux 1&2, Fresnes 3D, Cachan 1&2, Villejuif 1&2, Supplementary data A) were selected for thin-section manufacturing. These thin sections were impregnated with blue resin to highlight the porosity and easily differentiate it from the diagenetic calcite cements. Samples were extracted from the drilling mud to obtain thin sections from cuttings. Once dry and clean, they were sorted (clast diameter of more than 1 mm) to extract the coarse (1 mm) fraction of each sample. These 157 thin sections enabled us to determine facies, facies associations, pore space, and diagenetic cements in the reservoir. These thin sections were observed and described under an optical microscope (Leica DM750 P). A Micro-Computed Tomography scanner (µCT scan) was used to investigate the 3D pore space at high resolution in a sample in the Vert-Saint-Denis 1 well (VSD-1_1712,9m). The µCT scans were performed at IFPEN with a microtomograph X EasyTOM from RX Solution.

Following the definition of [START_REF] Embry | Two approaches to sequence stratigraphy[END_REF], sequence stratigraphic interpretations were performed, integrating data from 168 wells (Fig. 1). Maximum flooding surfaces (MFS) mark the boundary between deepening-upward and shallowing-upward trends. Subaerial unconformities or maximum regressive surfaces correspond to sequence boundaries (SB), labelled following the definition of recently compiled for the Bathonian of the Paris Basin. For the sequence stratigraphy and biostratigraphic framework, the wellknown Coulommes oilfield discovery in 1958 was taken as reference [START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF].

Brachiopod fauna in the Coulommes 001 well (CS001) provided a reference for sequence stratigraphy in the Middle Jurassic limestones [START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF]. Wells were correlated between Coulommes and Berthuis (BTS001) by [START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF] and throughout the basin [START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF][START_REF] Delmas | Petrophysical Properties of the Middle Jurassic Carbonates in the PICOREF Sector (South Champagne, Paris Basin, France)[END_REF][START_REF] Houel | Synthèse des données géologiques sur le bassin de Paris (Rapport IFP No. Rapport n° 63030[END_REF]. As brachiopod levels correspond to maximum flooding surfaces (MFS) in this platform, [START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF]

used 34

MFSs to define third-order sequences (Fig. 3). Maximum flooding surfaces are recognized in well-logs by Gamma-Ray maxima (Fig. 3). In this study, we also included the location of MRSs, identified, from bottom to top, by the onset of an increase of Gamma-Ray (Fig. 3). MRS corresponds to the top of a shallowing-upward trend, for example top of ooid-bioclast limestone units (Fig. 3). A total of 15 MRSs and 15 MFSs were correlated from BTS001 to Champigny-sur-Marne 1 well (CSM1D) (Fig. 3). These surfaces were propagated along the 168 wells of the study area. Following the labeling of [START_REF] Garcia | Les associations de brachiopodes du Jurassique moyen du bassin de Paris; une echelle biochronologique ponctuee de niveaux-reperes pour la contrainte des correlations sequentielles a haute resolution[END_REF], [START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF][START_REF] Hardenbol | Mesozoic and Cenozoic sequence chronostratigraphic framework of European basins[END_REF], two MRSs (Bj3 and Bj4) are defined in the Bajocian, six in the Bathonian (Bt0 to Bt5) and six in the Callovian from Ca0 to Ca4 (Fig. 3). including the well-log of BTS001, is defined by [START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF]. Maximum flooding surfaces (MFS) with biostratigraphy in Coulommes oilfield appear in blue continuous lines, without biostratigraphy in blue dotted lines. 3.3 Facies,stratigraphic,and petrophysical 3D modeling 3.3.1. 3D grid creation Eighteen stratigraphic surfaces were defined on well sections: MRSs from Bj3 to Ca4, the iron oolite layer (condensed MFSs Athleta and Lamberti), the Toarcian/Aalenian top, the Toarcian/Aalenian base and the basement top. Of these 18 stratigraphic surfaces, 14 (surfaces Bj4 to Ca4) were used to define 13 stratigraphic zones in the model (sequences MJVII, MJVIIIa, MJVIIIb, MJVIIIc, MJIV, MJXa, MJXb, MJXI, MJXIIa, MJXIIb, MJXIII, MJXIV, MJXV, Fig. 3). The other surfaces were not within the stratigraphic framework of the model. Four facies associations were encoded in 168 wells using well logs, well reports, cuttings, and core samples. Two cross-sections at a regional scale were produced, one oriented North-South and one oriented East-West. To ensure compatibility with other software, the corner point gridding algorithm (based on the pillar gridding algorithm) was used to create the 3D grid.

Faults were modeled using the structural maps provided by the BRGM [START_REF] Beccaletto | Overview of the subsurface structural pattern of the Paris Basin (France): Insights from the reprocessing and interpretation of regional seismic lines[END_REF][START_REF] Dentzer | Thermal and seismic hints for chimney type crossstratal fluid flow in onshore basins[END_REF]. They were then transformed as fault sticks in Petrel© and used as inputs in the structural and fault framework of the model. Cell size was set to 100 m in the I and J (horizontal) directions. A total of 52 layers were modeled with 3.76 m mean thickness in the K direction. This resulted in a grid of 12.2 million cells.

Property upscaling

The facies log was upscaled using the "most of" method without weighting. This method took the log values going through the cell and encoded the cell with the most commonly represented value. In total, the model has 6384 facies cells encoded using this method.

The well-interpreted and porosity-calibrated NPHI (28 logs), NMR (1 log), and NEUT (7 logs), as well as the min-core data (40 logs) were combined and upscaled to maximize cells containing porosity data (1686 cells). Sixty-nine porosity logs were considered for upscaling porosity in the model (Supplementary data A). Mini-cores were the calibrating tool for neutron-based data. Mini-cores and TCMR data were prioritized over NPHI-based data. The porosity logs were then upscaled using the arithmetic averaging method. This method takes the average value of all the data points of a well crossing a cell. Thirty-eight permeability logs were used for upscaling (37 permeabilities from mini-cores and one from NMR). The permeability data were upscaled using the geometric averaging method. The temperature and salinity properties were upscaled using the arithmetic averaging method. Temperature data is given in 65 wells at the top of the reservoir, measured at the beginning of the geothermal production. Only one well presents a continuous temperature log.

Facies, porosity, and permeability modeling

The "Truncated Gaussian with Trends" algorithm was used to model facies associations, stratigraphic zone by stratigraphic zone. This algorithm allowed us to define boundaries between facies associations. Facies were locally refined to better fit the stratigraphic well correlation diagrams (brush tool). This step was critical as the petrophysical properties are highly correlated to the facies [START_REF] Nurmi | The land of giants[END_REF][START_REF] Borgomano | Stratigraphic well correlations for 3-D static modeling of carbonate reservoirs[END_REF][START_REF] Tomassetti | Modeling lateral facies heterogeneity of an upper Oligocene carbonate ramp (Salento, southern Italy)[END_REF][START_REF] Trippetta | Carbonate-ramp reservoirs modelling best solutions: Insights from a dense shallow well database in Central Italy[END_REF].

Facies is commonly accepted as a factor controlling porosity in the Middle Jurassic carbonates of the Paris Basin (Brigaud et al., 2014b). For porosity, and depending on the facies, the transformation part of the Data Analysis with Petrel© was used to prepare the data and determine the distribution used in the model (i.e., input truncation and normal score). For the normal score transformation, the "define curve" option was used to match the porosity distribution with facies from the mini-core data in the wells and to apply it statistically to the entire volume of the model. The same approach was used for permeability by adding a logarithmic transformation to account for the porosity versus permeability logarithmic relationship.

To detect the anisotropy of the data (porosity and permeability), a variogram map based on the upscaled porosity was created. Variograms were conducted on the overall model for each facies association to populate porosity in the grid, and then zone by zone, facies by facies in 2 key zones of the reservoir with high data density. The same process was used for permeability.

Results

Sedimentary facies

Thirteen facies are identified for the entire study, grouped into four depositional environments (FA -Facies Associations) ranging from lower offshore to lagoon (Fig. 4, Fig. 5, Table 1). Facies associations have been defined for the eastern part of the Paris Basin in the synthesis of Brigaud et al. (2014a) on Jurassic formations. The distribution of facies associations is shown on a theoretical profile from the lower offshore in the outer ramp to lagoon in the inner ramp (Fig. 5). and6A, andTable 1). Facies F1a is composed of clays and silts (Table 1). Faunas are composed of ammonites, belemnites, brachiopods, and echinoderms. Facies F1a occurred mainly during the Callovian, and locally during the late Bajocian and early Bathonian (Figs 3 and4). Facies F1b is composed exclusively of clays (Figs 4, 5 and6A, Table 1). Although these two facies mainly appear in the Callovian deposits, forming the cap-rock of the aquifer, some local appearances are also detected in the late Bajocian deposits.

Interpretation -The absence of hydrodynamic sedimentary structures, the presence of ammonites and belemnites, bivalves or brachiopods in life position, and abundant bioturbation indicate very low energy conditions in an open distal marine environment, below the storm wave base, probably at a water depth of more than 40 m [START_REF] Andrieu | Disentangling the control of tectonics, eustasy, trophic conditions and climate on shallow-marine carbonate production during the Aalenian-Oxfordian interval: From the western France platform to the western Tethyan domain[END_REF] (Fig. 5). Considering the predominance of clays and marls in this facies, effective porosity and permeability are expected to be very low. This facies forms the cap rock of the reservoirs and locally impermeable layers in the Calcaires marneux à Pholadomyes. (VLG-4, 1525 m)). B-Core picture of marls (facies F2a, upper offshore) from the Vert-Le-Petit-4D well (VLP-4D, 1924 m). C-Oo-bioclastic grainstone (facies F3b and F3c, shoal) from the Vert-Saint-Denis-1 (VSD-1) core ( 1716

m). D-Intraclast grainstone (facies F4d, lagoon) from the Vert-Saint-Denis-1 (VSD-1) core (1694 m). E-Peloidal wackestone (facies F4a, lagoon from the Sergines-1 (SER-1) core (1656 m).

4.1.2 Alternating marls and mud-dominated limestones -Upper offshore environment (facies association FA2) Description -This facies association contains three facies: (1) alternating marls and peloidal wackestone to packstone (facies F2a), (2) siliciclastic bioturbated packstone (facies F2b), and

(3) marls with iron ooids (facies F2c), (Table 1, Figs 4, 5 and6B). Wackestone to packstone layers in Facies F2a are composed of peloids, quartz, clays, foraminifera, bivalves, echinoderm fragments, crinoids and gastropods. The main characteristic of facies F2a is the clear alternations of marl layers and packstone layers (Table 1; Figs 5 and6B). Facies F2b is roughly composed of the same bioclasts as the packstone of facies F2a but is more continuous vertically (Table 1;Fig 5). Facies F2c consists of marly deposits rich in iron ooids (Table 1). This facies also includes peloids and quartz grains. Facies F2c is present in all wells, forming a thin (metres-thick) and very well-defined marker across the Paris Basin. The ammonite fauna dated this layer from the late Callovian (lamberti Zone) to the early Oxfordian (cordatum Zone) [START_REF] Vail | La stratigraphie sequentielle et son application aux correlations chronostratigraphiques dans le Jurassique du bassin de Paris[END_REF][START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF][START_REF] Collin | Depositional environments and iron ooid formation in condensed sections (Callovian-Oxfordian, south-eastern Paris basin, France)[END_REF].

Interpretation -Silty, marly or mud-dominated facies mark a calm environment. Some sandy or bioclast-rich layers are slightly graded with undulating draping lamination, suggesting storm-graded beds. The depositional zone is situated between the storm wave base and the fair-weather wave base, probably at a water depth of 20 to 40 m [START_REF] Andrieu | Disentangling the control of tectonics, eustasy, trophic conditions and climate on shallow-marine carbonate production during the Aalenian-Oxfordian interval: From the western France platform to the western Tethyan domain[END_REF] and [START_REF] Brigaud | Sedimentary architecture, depositional facies and diagenetic response to intracratonic deformation and climate change inferred from outcrops for a pivotal period (Jurassic/Cretaceous boundary, Paris Basin, France)[END_REF], Table 1, Fig. 4). Considering the predominance of clays and marls in this facies, effective porosity is expected to be very low as is the effective permeability. This facies association mainly occurred in the Calcaires marneux à Pholadomyes and Dalle Nacrée 442 Formations (Fig. 4A) 443

Figure 7. A-Peloidal packstone (facies F3a) with bioclasts and ooids displaying brachiopods (Bra.), peloids (Pel.), echinoderm fragments (Ech.), coral fragments (Coral). Porosity appears in blue and is low in this sample. Some micro-porosity (Mic. Por.) and styloliths (Sty.) are present. B-Low to moderately cemented ooid grainstone (Facies F3b) with rare peloids (Pel.), bioclasts such as echinoderm fragments (Ech.) locally form nucleus of ooid, isopachous cements (Iso. Cem.) are well-developed.

Rhombohedral pore spaces are preserved as flat surfaces (fs) printing the dissolved dolomite rhomb faces or as points (p) printing the dissolved dolomite rhomb terminations (see arrows and zoom with possible location of dissolved dolomite rhombs in red dotted line). Dissolution (Dissol.). The porosity appears in blue due to the resin impregnation. C-Highly cemented oopeloidal (Pel.) grainstone (facies F3c) with echinoderm fragments (Ech.), crinoids (Cri.), bivalves (Biv.), formaminifera (For.). Micro-porosity (Mic. Por.) is present in ooids. Bc. : Blocky calcite. D-Rudstone to grainstone with bivalves (Biv.) and intraclasts (Int.), (facies F3d), oncoids (Onc.), echinoderm fragments (Ech.), brachiopods (Bra.), bryozoans (Bryo.), intraclasts and syntaxial overgrowth calcite cement around echinoderm fragments (Syn.). Porosity (in blue) is variable but can sometimes be very substantial. E-Peloidal (Pel.) wackestone/packstone with foraminifera (For.) and oncoids (Onc.), facies F4a. Blocky calcite (Bc) fills macroporosity. Micro-porosity is still present. F-Low to moderately cemented peloidal (Pel.) grainstone (facies F4b) with echinoderm fragments (Ech.), bivalves (Biv.), intraclasts (Int.) and oncoids (Onc.). Syntaxial overgrowth calcite cement around echinoderm fragments (Syn.) and isopachous cements (Iso. Cem.) reduce the porosity (in blue), while dissolution (Dissol.) enhances it. G-Oncolithic (Onc.) floatstone (facies F4c) with foraminifera (For.). H-Rudstone with intraclasts (Intra.), facies F4d, peloids (Pel.), echinoderm fragments (Ech.), crinoids (Cri.), brachiopods (Bra.) and bryozoans (Bryo.). Blocky calcite (Bc) is present, and micro-porosity (Mic. Por.) is current in the cortex of some grains. SER-1 corresponds to Sergines-1 well and VSD-1 to Vert-Saint-Denis-1 well.

Ooids or bioclast-dominated facies -Shoal environment (facies association FA3)

Description -This facies association encompasses four grain-supported facies: (1) peloidal packstone with ooids and bioclasts (facies F3a), (2) low to moderately cemented ooid grainstone with bioclasts (facies F3b), (3) highly cemented ooid-peloidal grainstone with bioclasts (facies F3c) and ( 4) ooid-oncoidal rudstone to grainstone with bioclasts and intraclasts (F3d) (Figs 4, 5, 6C-E and 7A-D). Facies F3a is a packstone composed mainly of peloids, bivalves, echinoderm fragments, ooids and rare oncoids, coral fragments, intraclasts, foraminifera (in ooids), bryozoans and gastropods (Fig. 7A). Grain sizes range from 100 µm to 2 mm (Fig. 7). In this facies association, it is the only facies with a bounding matrix phase and it has a packstone texture. The matrix is unevenly distributed in the facies, and interparticle porosity can sometimes be observed. Facies F3a displays compaction textures, reducing the effective porosity and permeability (Table 1, Fig. 7A). Facies F3b is a grainstone mainly composed of 200-500 µm-large ooids. Peloids are present, together with oncoids, intraclasts and bioclasts (bivalves, echinoderm fragments, foraminifera, bryozoans, and gastropods, Table 1, Figs 5, 6C, 7B). Large interparticle pores (100-500 µm) form the macroporosity considered as the effective porosity, contrary to the porosity under 10 µm considered as microporosity and non-effective porosity. Early isopachous cements, typical of early marine cement, are present. This early cement limits mechanical compaction, maintaining large intergranular pore space with sharp borders. Limited blocky calcite zones are also present (Fig. 7B). Locally, part of the ooids and blocky calcite have been dissolved, suggesting a late dissolution process. Rhombohedral pore spaces of about 500 μm in size are locally observed, labelled as mouldic pores (Fig. 7B). Facies F3c is mainly composed of ooids, with common content of peloids, oncoids, bivalves, and rare intraclasts, echinoderm fragments, foraminifera, bryozoans, gastropods, and coral fragments. This facies is very similar to F3b, but differs with blocky calcite cement filling intergranular space (Fig. 7C). Only microporosity in some ooids is present (Table 1, Fig. 7C, Fig. 7a). F3d is a grainstone to rudstone composed of common intraclasts, gastropods, bivalves, brachiopods, and echinoderm fragments, with rare ooids, peloids, oncoids, bryozoans and foraminifera (Fig. 7D). This facies was deposited in thin layers with millimetric bivalve and brachiopod shell fragments. Millimeter-sized blocky to syntaxial overgrowth calcite cement around echinoderm fragments fills the intergranular space. Dissolution forms 200 µm to 2 mm pores (Fig. 7D). CT scan 3D models of the porous network show that, despite a high porosity value, pores are not well connected in facies F3d (Fig. 8B). At the microscale, this facies displays moderate porosity and moderate permeability around 100 mD. Porosities extracted from the CT scan have similar variations in all three directions, ranging from 2 to 12% (Fig. 8). Horizontal permeability is anisotropic with a minimum of 70 mD and a maximum of 145 mD. Vertical permeability is close to the maximum horizontal value (Fig. 8). This facies presents a good effective porosity and permeability overall. Still, it is less common than F3b in this facies association (Table 1).

Interpretation -The presence of echinoderms, brachiopods, or bryozoans suggests welloxygenated water with normal salinity [START_REF] Brigaud | The origin and timing of multiphase cementation in carbonates: Impact of regional scale geodynamic events on the Middle Jurassic Limestones diagenesis (Paris Basin, France)[END_REF]). Common sigmoid cross-bedding in sandwaves and grainstone/rudstone textures indicates shallow and high-energy wavedominated conditions in an oolitic shoal complex (Fig. 4). Based on pore development observed in thin sections, two groups of facies have been distinguished. Facies F3a and F3c are considered to have low reservoir potential due to the presence of a mud-matrix, high compaction (F3a) and the development of highly blocky calcite (F3c). Facies F3b and F3d are the main reservoir facies in the shoal environment. Facies F3b is more common than F3d, which only appears locally in thin but porous and permeable layers. Rhombic pseudomorphs can be locally recognized by careful observations of the pore boundaries, either forming flat surfaces along ooids or forming points, printing the crystal terminations (Fig. 7B). Pores with rhomb shapes suggest dissolution of dolomite in facies F3b, as already observed by [START_REF] Purser | Dedolomite Porosity and Reservoir Properties of Middle Jurassic Carbonates in the Paris Basin, France[END_REF], (Fig. 7B). The dolomite is entirely dissolved, while blocky calcite and clasts have been only partially affected by dissolution. Dissolution appears to be an essential process of porosity development in the reservoir. This facies association mainly occurred in the Oolithe Blanche Formation (Fig. 4A).

Peloid to oncoidal facies -Lagoon environments (Facies association FA4)

Description -This facies association encompasses four grain-supported facies with peloid dominant non-bioclastic components: (1) peloidal wackestone to packstone with foraminifera and echinoderms (facies F4a), (2) low to moderately cemented peloidal grainstone (facies F4b), (3) oncoidal packstone to floatstone (facies F4c) and ( 4) rudstone with intraclasts (facies F4d) (Table 1, Figs 7E-H). Facies F4a is mainly composed of peloids, with common oncoids, foraminifera (presence of Meyendorfina bathonica), echinoderm fragments, bivalves and rare ooids, bryozoans, gastropods, brachiopods and coral fragments. Grain size ranges from 200 to 800 µm. This facies is mud-dominated, with a wackestone to packstone texture and it displays low effective porosity and permeability, mainly microporosity in the matrix between the grains (Table 1, Fig. 6E). Facies F4b is mainly composed of peloids but differs from facies F4a by the absence of a mud matrix. This facies has a grainstone texture and exhibits high effective porosity and permeability (Fig. 6F). Grain size varies from 100 to 400 µm, displaying 50 to 500 µm pore sizes with very high microporosity. Cementation is directly related to the echinoderm content, around which syntaxial cement systematically develops. As echinoderm fragments are rare, cementation is not well developed in facies F4b (Fig. 6F). Local dissolution of grains improves the porosity of this facies. Facies F3c is mainly composed of oncoids (grain size from 100 µm to 5 cm), with common peloids, and rare bivalves, foraminifera and gastropods. This facies is matrix-supported with a packstone to floatstone texture. No macro-pores are observed in facies F4c, resulting in very low effective porosity and permeability (Table 1, Fig. 7G). Facies F4d is composed mainly of intraclasts, with common oncoids, bivalves, brachiopods and rare foraminifera, gastropods and green algae. Pores are rare due to calcite overgrowths (Table 1, Fig. 7H).

Interpretation -The dominance of peloid-oncoid dominated facies and the muddy matrix indicates low-energy environments, probably within a lagoon. As for the shoal environment, two groups can be distinguished in this facies association. A first group is composed of the non-porous and impermeable facies F4a and F4c (permeability from 0.1 mD to a few mD in the Vert-Saint-Denis-1 well). A second group displays high effective porosity and permeability and contains F4b (permeability from 0.5 mD to >1500 mD in the Vert-Saint-Denis-1 well).

Facies F4d may belong to the first or second group depending on its degree of cementation. This facies association shows that, even in a sheltered environment like a lagoon, porous and permeable rocks are present and must be modeled. Petrographically, facies association FA4 has a lower reservoir value than facies association FA3, but some local non-cemented granular facies (F4b) may form productive layers. This facies association mainly occurs in the Calcaire de Comblanchien Formation (Fig. 4A), dated as Late Bathonian by Meyendorfina bathonica [START_REF] Thomas | Contribution of drone photogrammetry to 3D outcrop modeling of facies, porosity, and permeability heterogeneities in carbonate reservoirs (Paris Basin, Middle Jurassic)[END_REF][START_REF] Vincent | Giant subaqueous carbonate dunes: a revised interpretation of large-scale oo-bioclastic clinoforms in the middle Jurassic of the Paris Basin and its implications[END_REF]. The SW-NE stratigraphic cross-section from the Etrechy-1D to Berthuis-1 wells allows us to reconstruct vertical and lateral facies heterogeneity (Fig. 9A). The Berthuis-1 well is located outside the area modeled in this study. It allows to correlate the stratigraphic sections of the model with the facies and biostratigraphic attribution of the Coulommes oil field [START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF]. The Early Bajocian sequences beneath Bj4 are not modeled. They exhibit alternating marls, clays, and packstones (Fig. 9B). Between Bj4 and Bt2 (sequences VIIIa, VIIIb, VIIIc), the facies associations described on wells mainly correspond to marls and packstones of the upper offshore environment (facies association FA2), and locally oolitic deposits (facies association FA3) at the Croissy-Beaubourg-1D and Etrechy-1D wells. In the Vert-le-Grand-6D well, oolitic facies are absent, and clay layers (lower offshore) are predominant (Fig. 9). Between Bt4 and Bt2 (sequences MJIX and MJXa), the dominant facies along the cross-section are ooid limestones (facies association FA3). Marl-limestone alternations are still present at the base of each sequence, corresponding to the flooding of the platform (Fig. 9). In the western part of the cross-section (between the Etrechy-1D and Vert-le-Grand-6D wells), alternating marls and limestones (facies association FA2) dominate and form a trench (Fig. 9). These marllimestone deposits vastly reduce the reservoir properties of the area. Between Bt5 and Bt4

(sequence MJXb), the Calcaire de Comblanchien Formation (facies association FA4) is present between the Le Petit Quincy-1 well and the Croissy-Beaubourg-1D well between Bt5 and Bt4 surfaces (Fig. 9B). This formation is thicker towards the Berthuis-1 well, where the oolitic shoal facies (facies association FA3) passes laterally into lagoonal facies (from Bt2 to Bt5 surfaces).

West of the lagoon (west of the Le Petit Quincy-1 well), facies are dominated by oolitic deposits (at Vert-le-Grand-6D). In sequence MJXb, westwards, the dominant facies association is alternating marl/packstone (upper offshore environment: facies association FA2), corresponding to the distal part of the carbonate ramp. In sequence MJXb, oolitic limestones correspond to a shoal environment. Between Ca2 and Bt5 surfaces (sequences MJXc to MJXIII), the Dalle Nacrée formation consists mainly of alternating marl and limestone, with local ooid or bioclastic limestones at Berthuis between Bt5 and Ca1 or at Vert le Grand and Le Petit Quincy between Ca2 and Ca0. The cross-section shows the homogeneity of the cap rock (middle Callovian deposits) in the section (Fig. 9). From Ca4 to Ca2 surfaces, deposits are mainly composed of clays (lower offshore facies association FA1).

The second transect illustrates the facies architecture between Cachan in the northwest and the historic oil field of Chailly-en-Bière in the southeastern part of the study area (Fig. 1B, Fig. 10). and Grigny-1 wells (Fig. 10). A lagoon (facies association FA4) developed at the top of the sequence in the Grigny well. In the western and northern parts of the study area, sequence MJXb deposits are composed exclusively of lagoonal facies. In this sequence, the oolitic limestones form an oolitic shoal or barrier, protecting sheltered lagoons located around the Cachan-1, Grigny-1, and Brie-119 wells. This cross-section exhibits a similar architecture for the middle Callovian deposits as the E-W cross-section (Fig. 9), ranging from Ca0 to Ca4. An exception is found at the base of the lower Callovian deposits near Brie-119 and Villoison-1, with the oolitic facies corresponding to the reservoir of the Dalle Nacrée formation [START_REF] Javaux | La plate-forme parisienne et bourguignonne au Bathonien terminal et au Callovien: dynamique sedimentaire, sequentielle et diagenetique: place et creation des reservoirs potentiels[END_REF]. These deposits are usually framed between Bt5 and Ca1.

The third-order depositional sequences described displaying a regressive stacking pattern.

Facies pass upwards from (1) lower offshore marls or clays into (2) upper offshore marl/limestone alternations, (3) ooid limestones, and (4) lagoonal peloid to oncoidal facies (Figs 9 and 10). Sequences never contain the complete succession of facies, and most of them form vertical successions passing upwards from marl/limestone alternations into ooid limestones. Sequences MJIX, MJXa and MJXb are the most complete sequences with a vertical succession passing upwards from marl/limestone alternations into ooid limestones and peloid to oncoidal facies (Figs. 9 and 10). Using the two cross-sections and a view of the upscaled cells in each well, facies associations in each stratigraphic sequence have been computed using the "truncated Gaussian with trends" algorithm. This algorithm is well suited for modeling large-scale ordered facies progradation and retrogradation, defining transgressive and regressive systems tracts (Brigaud et al., 2014b). For each zone, this algorithm requires a workflow to be followed (Brigaud et al., 2014a): (1) choose the facies and check the facies order, (2) set up the depositional geometry, i.e., edit the transition lines between the facies (boundaries of the facies belts), and (3) decide on the variogram ranges (spatial extension of facies belts) and variances (degree of interfingering between facies belts). The "truncated Gaussian with trends" algorithm is used to model the evolution of the depositional system in time and space in each zone. Detailed cross-sections and palaeogeographic maps from hard data (e.g., facies

Facies modeling

in wells) at the bottom and base of each sequence are used to design facies trends and honor hard data from each well. To control the quality of facies modeling, we compared the facies association proportion in wells to ensure they are statistically compatible (Fig. 11) with the percentage of each facies association in the upscaled cells in the wells and in the entire grid.

In well logs, facies association content is very homogeneous: clays (facies association 1) = 16.9%; marl/limestone alternations (facies association 2) = 37.7%; ooid limestones (facies association 3) = 34.6% and peloid to oncoidal facies (facies association 4) = 10.8%. Facies contents (in percentage) in upscaled cells and in the entire model are similar to well logs, except for the marl/limestone alternations in the entire model, which seem to be overestimated (43.9%), at the expense of clays (facies association FA1). Since these rocks are not considered reservoirs, no impact is expected on the reservoir location. Facies evolution can be observed in 2D maps or in a 3D block-diagram along the model (Fig. 12). It appears that in the modeled area, the early-late Callovian deposits act as cap rock for the reservoir, with the partial exception of the Brie region, where some ooid limestones are present (Fig. 12). In the north-western area, the anticline structure of Beynes-Meudon forms a major structural feature of the geothermal reservoir (in the Cachan area). The ramp was deepening westward as attested by the progressive transition from the muddy lagoon and ooid shoals (facies associations FA3 and FA4) to upper offshore marl-limestone alternations (facies association FA2) (Fig. 12). The lagoonal facies association is mainly present in the eastern and northern parts of the model, and it progrades westward in sequence MJXb. Deposits of sequences MJIX and MJXa form a good reservoir in most of the study area. They are mainly composed of ooid limestones, with some marl and limestone alternations at the bottom of the sequences. The reservoir facies passes into marl-limestone alternations in the western part of the study area (Fig. 12 B,C). Below the Bt2 surface, the likelihood of encountering reservoir facies is low, except locally in the eastern part (Fig. 10), with more distal deposits ranging from marls to clays. data in shales [START_REF] Serra | Diagraphies différées -bases de l'interprétation[END_REF]. D-Mini-core and Schlumberger Doll Research (SDR) permeability upscaled data sorted by facies association. The line inside the box is the median, and the top and base of the box correspond to the 25 and 75% quartiles.

Petrophysical characteristics

The horizontal bars represent the 5% (for the bottom bar) and the 95% (for the upper bar) of the data range. The dots represent the data out of range.

The porosity measured on mini-cores ranges between 0% and 39% (Fig. 13A). The best porosities are found in ooid limestones (3rd quartile at 17%). However, they may also be particularly low (<5%, Fig. 13A). In ooid limestones, the median of mini-core data indicates 14% porosity confirming that these facies are of good reservoir quality (Fig. 13A). The lagoonal facies have lower median porosity (8.5%), but may locally display better porosity values, up to 25% (Fig. 13A). Marl-limestone alternations have a better median porosity than peloid to oncoidal limestones (median of 9.5%, Fig. 13A). Clays/marls have the lowest effective porosities (median of 4%, Fig. 13A). Note that few data are available for the clayey/marly facies (N = 18, Fig. 13A). The marl/limestone alternation and clay/marl facies present in the upper and lower offshore deposits are less porous overall and considered as non-reservoir or cap rocks.

The permeability measured on mini-cores ranges between 0.01 mD and 8783.5 mD (Fig. 13B).

Permeability values of lagoonal and oolitic limestones present very large ranges from 0.01 mD to 8711 mD for the lagoonal limestones and from 0.01 mD to 8783.5 mD for the oolitic limestones (Fig. 13B). The median value in oolitic limestone is 61 mD. In contrast, the median value of permeability in lagoonal facies is 24 mD, showing that these two sedimentary facies are of good reservoir quality (Fig. 13B). Marl-limestone alternations and clays have significantly lower permeabilities (median of 1.3 mD and 0.1 mD, Fig. 13B).

NPHI logs are known to be unrealistic in marls and clayey deposits if treated for limestone porosities [START_REF] Serra | Diagraphies différées -bases de l'interprétation[END_REF]. Mini-core porosity never exceeds 20%, whereas the 4 th quartile has an unrealistic value of 28% for marl-limestone alternation (Fig. 13C). Porosity larger than 10% in clayey deposits and larger than 15% in marl-limestone alternations has been filtered by data analysis during modeling (Fig. 13C). Upscaled porosity of ooid limestones in wells is lower (median at 11%, Fig. 13C) than mini-core porosity and porosity of the lagoonal facies has similar porosity in the upscaled wells as mini-core porosity (median value of 9%, Fig. 13C).

The upscaled permeabilities range between 0.1 mD and 4000 mD (Fig. 13D). The lagoonal and oolitic limestone permeabilities (mini-cores and kSDR) are similar (median at 18.5 mD and 19.1 mD respectively, Fig. 13D). Marl-limestone alternations and clays have significantly lower upscaled permeability values (median of 0.64 mD and 0.1 mD).

Horizontal and vertical porosity and permeability are positively correlated, except for a few samples that reveal larger vertical permeabilities than horizontal ones (Fig. 14). Porosity vs. logarithmic permeability relationship is presented in Figure 15 according to facies.

The data available for the clay/marl facies association (FA1) are not sufficient to determine a possible correlation between porosity and permeability (Fig. 15A). The marl/limestone alternation facies association (FA2) presents two datasets, (1) mini-core data and (2) porosity from TCMR versus kSDR data originating from GCAH-2 (Fig. 15B). Overall, the data in the FA2 is concentrated between 5 to 15% porosity, but with quite low permeability (0.1 to a few mD).

The ooid/bioclast limestone facies association (FA3) presents a cloud of data clustered in a zone of higher porosity and permeability (Fig. 15C). This cloud is present between 10 and 20% of porosity and from 10 to 1000 mD permeability (Fig. 15C). This trend is mainly observed on the mini-core data. In contrast, TCMR/kSDR tend to show a more homogenous distribution (Fig. 15C). In the peloid to oncoidal limestones, the sample distribution follows a linear regression between 8 and 25% of porosity (Fig. 15D). The lower porosity range, from 2 to 8%, follows a steeper linear regression line (Fig. 15D). Overall, there is a good overlap between mini-core data and data from the TCMR log (Fig. 15D). They follow a linear regression line (with permeability expressed as a logarithm) and display a symmetric dispersion for facies associations 2 (marl/limestone alternations), 3 (ooid/bioclast limestones) and 4 (peloid to oncoidal limestones) (Fig. 15). Of the 12.2 million cells of the entire grid, 6384 cells have upscaled facies, 1686 cells have upscaled porosity and 469 cells have upscaled permeability. Permeability data are only available in a small number of wells (34 wells) as it is not a property that could be easily logged in ancient wells. Porosity, on the other hand, was often logged using mini-core plugs or neutron-derived logs. All these cells with upscaled data represent cells where original data When calculating the anisotropy, the lag distance parameter is very important as it will determine the resolution and accuracy steps for the studied property. For the anisotropy analysis, as the horizontal lag distance is 50 m, the value is rounded to 50 m, and the uncertainty is +/-25 m (Fig. 17C, Table 2). In the vertical direction, the lag distance is set from 1.5 to 6 m, depending on the facies association and the thickness of the studied zone. The well-log data from the mini-cores are taken as a standard to calibrate the "Data Analysis" The Gaussian Random Function Simulation algorithm propagates porosity within the grid following the anisotropy of porosity for each facies association of the general model or following the anisotropy of porosity for specific facies associations in the MJXb and MJXa sequences (Table 2). The variograms used to populate porosity in ooid/bioclast limestones are mainly elongated perpendicular to the depositional slope, with a general relationship of 1300 m x 1050 m and a thickness of 6 m in all sequences (except for MJXa and MJXb), 1850 m

Porosity and permeability distribution in the grid

x 1150 m, and a thickness of 8 m for sequence MJXb, and 1150 m x 850 m and a thickness of 4 m for sequence MJXa (Table 2). The variogram for marl/limestone alternations also shows a correlation distance mainly elongated perpendicular to the depositional slope, with a general relationship of 1500 m x 1000 m and a thickness of 3 m in all sequences. Anisotropy is low for mud-dominated limestones with variograms giving elongated porous layers of 1000 m x 850 m and a thickness of 4 m for all sequences (Table 2). The porosity variogram for clay facies does not exhibit anisotropy, with a correlation distance size of 2150 m x 2150 m (Table 2).

Permeability is modeled similarly using distributions of permeability in well-logs, spatial variograms (Table 2), and the Gaussian Random Function Simulation algorithm (Brigaud et al., 2014b). As the permeability is correlated with the facies and porosity (Fig. 15), a co-kriging of permeability against porosity is used in the Gaussian Random Function Simulation.

Permeability and porosity are co-kriged using a varying local mean as a correlation and porosity as a secondary variable (Fig. 15). Each facies association is submitted to its respective variogram to introduce some heterogeneity into the reservoir (Fig. 16A-B). In clays, the same variograms are used to populate porosity and permeability. Ooid/bioclast limestones are mainly elongated perpendicular to the depositional slope, with a general relationship of 1600 m x 1100 m and a thickness of 4 m in all sequences (Table 2). The variogram for marl/limestone alternations also shows a correlation distance mainly elongated perpendicular to the depositional slope, with a general relationship of 1200 m x 900 m and a thickness of 3 m in all sequences. Anisotropy is low for mud-dominated limestones with variograms giving elongated porous layers of 850 m x 750 m and a thickness of 6 m for all sequences (Table 2).

The final porosity and permeability distributions are displayed in figure 16, showing the heterogeneity in the 3D grid along transects.

Histograms are used to compare the distribution of petrophysical properties between well logs, upscaled cells, and the total model classes after property modeling. These graphs help to perform quality control on the petrophysical properties (porosity and permeability) for each coded facies in the model to verify whether the data analysis and data transformation steps were correctly executed (Table 2; Fig. 18). Porosity and permeability can be extracted on a regional scale (Fig. 19) and averaged to quickly localize areas with high geothermal potential (Fig. 19). An area around Arpajon, from Palaiseau, Grigny, Vert-Le-Grand to Marolles-en-Hurepoix exhibits poor reservoir properties for MJVIIIb to MJXI (area 1 in Fig. 19A-B). On the other hand, the best reservoir quality is identified in the eastern part of the model, and south of Melun, with high-permeability patches reaching a few hundreds of mD (Area 3 in Fig. 19A-B). [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF] presented permeability maps indicating low global permeability in the western part of the studied zone with an overall upward trend. Between these two areas, intermediate porosities are present (area 2 in Fig. 19A). As the most productive reservoir facies are located in MJXa and MJXb sequences, average porosity and permeability in these two sequences allow us to refine the characterization of spatial heterogeneity of reservoir quality (Fig. 21C-F). In sequence MJXb, the area around Arpajon is unsuitable for geothermal exploitation (area 1 in Fig. 19D).

However, some porous and permeable patches are predicted. This corresponds to the westward progradation of ooids or bioclastic grainstone facies forming a shoal at the end of the Bathonian. An area from Creteil to the southern part of the model is characterized by the presence of many very permeable kilometer-scale patches (area 3 in Fig. 19D), from a few hundred mD (yellow) to 1 D (red, Fig. 19D). This area is the best suited for geothermal prospects in this sequence (Fig. 21C-D). The rest of the studied zone, north of Creteil and south-east of Melun (area 2 in Fig. 19D), is characterized by smaller, highly permeable patches (Fig. 19D). In sequence MJXa, a non-permeable area is clearly identified east of Palaiseau and Arpajon (area 1 in Fig. 19F). A triangular area, between Palaiseau/Arpajon and Evry (area 2 in Fig. 19F), presents highly heterogenous reservoir quality, with local permeable patches up to a few hundred mD (green to yellow, Fig. 19B), separated by more extensive areas of impermeable patches (blue, Fig. 19). It should be noted that the shape of the permeable lenses is robust, but their location is uncertain. They are not necessarily well located, perhaps 1 or 2 km away. In the rest of the modeled zone, the averaged petrophysical properties are very good and display numerous high-permeability patches (around one kilometer in size) separated by non-permeable kilometric patches (area 3 in Fig. 19F). A map has been constructed from the temperatures available in the well reports at the top of the reservoir using a convergent interpolation algorithm, helping to understand the temperature distribution (Fig. 20A). The south-eastern part of the study area does not have useful data, as temperature data in the area are issued from the drilling reports and are not usable as they are. Temperatures show high variability ranging from 58 °C in the northern part of the area to 82 °C in the western part (Fig. 20A). Also, in the southwest of the studied area, no trustworthy temperature could be extracted from the old well data available in this specific area. A transition from low temperatures of about 60 °C around Paris to high temperatures of 80 °C in the eastern part of the model is clearly visible (Fig20A). A geothermal gradient map using temperature and depth maps of the top of the reservoir has been extracted (Fig. 20B).

Temperature, salinity, transmissivity, and thermal gradient

The gradient resulting from this study at the top of sequence MJXb indicates some interesting variations concerning the raw temperature map (Fig. 20B). A corridor showing high gradients runs from Arpajon to the south of Torcy, with values around 40 °C/km, which is in the high range of geothermal gradients in the Paris Basin, and more generally in intracratonic basins [START_REF] Fabris | Les températures dans les forages -Gisements pétroliers du bassin Parisien et de l'Aquitaine[END_REF][START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF]. The southern part of the basin presents homogeneous values of around 36 °C/km (Fig. 20B). In the northern part of the model, however, the gradient drops drastically to 27 °C/km around Paris. Compared to temperature, a different pattern appears for salinity (Fig. 20C). A corridor of low salinity (less than 10 g/L) is present from Palaiseau through Evry to the south of the model, as already described by [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF] and [START_REF] Matray | Origin of the oil-field brines in the Paris basin[END_REF]. Higher values, reaching 26 g/L, are observed west of Creteil (Fig. 20C).

With this study database, we also calculated transmissivity distribution maps (Fig. 20D-F). We combined the permeability maps averaged for groups of sequences or individual sequences with the thickness maps to estimate the possible reservoir transmissivity (Fig. 20D-F). The transmissivity property is bound to the calculated permeability in the model, averaged between the cells of a sequence multiplied by the thickness of the sequence. In Bathonian limestones (sequences MJVIIIb to MJXI), the lowest transmissivity is located around Arpajon (area 1 in Fig. 20D) whereas the most transmissive regions are located in the eastern and southern parts of the modeled zone with transmissive patches of 50 D.m (area 3 in Fig. 20D).

This is because sequences MJVIIIc and MJIX are mainly composed of ooid/bioclast limestone facies in the eastern part of the model (Figs 10A and12C). In the west, sequences MJVIIIc and MJIX consist mainly of low transmissivity marl/limestone alternations (Figs 10A and 12C).

Intermediate transmissivity patches (about 10 D.m) are present around Créteil and north of Evry (area 2 in Fig. 20D), due to the smaller vertical thickness of ooid/bioclast limestones of sequences MJVIIIc and MJIX (Fig. 12C). In detail, the MJXb sequence presents very good transmissivity around Palaiseau, between Creteil and Torcy, and west of Melun (area 3 in Fig. 19E). Area 2 also presents some interesting zones, forming scattered patches around 10 D.m (green in Fig. 20E). Transmissivity is low in the south-west and south-east of the model, due to the occurrence of mud-dominated facies without a granular facies unit (Area 1 in Fig. 20E).

The MJXa sequence presents very low transmissivity in the western and eastern parts of the model (area 1 in Fig. 20F). An intermediate area (area 2 in Fig. 20F) presents local patches of tens of D.m and a high potential area (area 3 in Fig. 20F) presents short-wavelength patches of 30 D.m.

Discussion

5.1 Temperature, salinity, transmissivity, and thermal gradient heterogeneity at the regional scale

The high heterogeneity of petrophysical properties is known [START_REF] Housse | Potentiel géothermique du bassin parisien (Rapport No. 74-7-0990[END_REF][START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF][START_REF] Le Nir | Ressources géothermiques du département de l'Essonne (91)[END_REF][START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF] and is confirmed by the supplementary data added to our model. Namely, temperature, thermal gradient, salinity, and transmissivity have been showing that geothermal gradient is heterogeneous. Therefore, another mechanism must explain the drop in observed temperatures in the northwest of the studied area. This local geothermal gradient variation from 27 °C/km to 40 °C/km (Fig. 20B) is not well understood.

To explain the heterogeneous geothermal gradient in the reservoir, [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF] suggested (1) possible heat flux variation or, from pressure gradient, (2) fluid flow in the reservoir such as cold horizontal water flows from west of the greater Paris area to the north of Paris. This horizontal fluid flow would, however, be opposite to the rest of the formation water flows [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF]). The regional flow, reconstructed by [START_REF] Menjoz | Flow of formation water in the Jurassic of the Paris Basin and its effects[END_REF],

shows water flow from the southeast to the northwest, going from the Middle Jurassic limestone outcrops in Burgundy and Lorraine areas (Eastern Paris Basin) toward the English Channel. Regional flow modeling at the scale of the Paris Basin shows a possible effect of flow velocity increase to explain a positive anomaly of + 10 °C south of Paris [START_REF] Menjoz | Flow of formation water in the Jurassic of the Paris Basin and its effects[END_REF].

However, on re-interpreted seismic profiles, [START_REF] Dentzer | Thermal and seismic hints for chimney type crossstratal fluid flow in onshore basins[END_REF] identified reduced seismic amplitude along chimneys corresponding to pervasive fractures. Cold water drawdown flows, from the Late Jurassic and Early Cretaceous aquifers, along these chimneys, could partly explain the 20 °C fall in temperature between the northern and eastern parts of the model [START_REF] Dentzer | Thermal and seismic hints for chimney type crossstratal fluid flow in onshore basins[END_REF].

The transmissivities are in the same range as published earlier [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF], from less than 5 D.m west of Palaiseau to 50-75 D.m around Créteil. However, the transmissivity maps extracted from our model provide higher resolution than what was known in the literature [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF], showing locations of small transmissive patches unknown hitherto (Fig. 20D-F). Transmissivity at the doublet scale is highly uncertain where no well data is available.

Nonetheless, this study helps to better determine interesting areas of high transmissivity

Discussion -facies and impact on reservoir quality

3D heterogeneity at the core scale

In the high reservoir-quality facies association (facies association FA3; ooid/bioclast limestones), the pore network is relatively heterogeneous. Whereas pore sizes of facies F3b are relatively homogenous, with very well connected pores between 100 and 500 µm (Fig. 7B), facies F3c, which presents a grain type similar to the excellent reservoir facies F3b, displays intergranular spaces filled by large blocky calcite, reducing the porosity to intraparticle micropores (Fig. 7C). The lack of micrite and the presence of well-formed ooids (diameter of 400 to 500 µm; Fig. 7B) indicate a high-energy environment, corresponding either to oolitic shoal or giant dunes [START_REF] Thomas | Contribution of drone photogrammetry to 3D outcrop modeling of facies, porosity, and permeability heterogeneities in carbonate reservoirs (Paris Basin, Middle Jurassic)[END_REF][START_REF] Vincent | Giant subaqueous carbonate dunes: a revised interpretation of large-scale oo-bioclastic clinoforms in the middle Jurassic of the Paris Basin and its implications[END_REF]. Diagenesis plays an essential role in the reservoir quality of this facies association. Isopachous cements protected the rock against any decrease in its intergranular pore space during compaction [START_REF] Purser | Early diagenesis and the preservation of porosity in Jurassic limestones[END_REF][START_REF] Brigaud | Acoustic Properties of Ancient Shallow-Marine Carbonates: Effects of Depositional Environments and Diagenetic Processes (Middle Jurassic, Paris Basin, France)[END_REF]Brigaud et al., , 2014a)). Dedolomitization enhanced the effective porosity of facies F3b, and probably F3d (Fig. 7B-D), as already observed in the Coulommes oilfield 20 km from the study area [START_REF] Purser | Dedolomite Porosity and Reservoir Properties of Middle Jurassic Carbonates in the Paris Basin, France[END_REF], or in the Massangis outcrop analogue in Burgundy [START_REF] Thomas | Contribution of drone photogrammetry to 3D outcrop modeling of facies, porosity, and permeability heterogeneities in carbonate reservoirs (Paris Basin, Middle Jurassic)[END_REF].

Facies F3d exhibits larger pores than facies F3b, from 100 µm to 2 cm. However, the pores are less interconnected (Fig. 7D, Fig. 8). Large calcite crystals grew between bioclastic grains, especially along echinoderm fragments, filling pore space. Facies F3c is another good example of the diagenesis effects on reservoir properties. The heterogeneity of carbonate cementation is challenging to predict. The modeling of individual facies is complex due to the scarcity of core samples in the study area. The mud-dominated limestone facies association also has heterogeneous reservoir quality. Facies F4a and F4c were deposited in a calm environment where carbonate mud could settle, forming muddy facies with poor reservoir qualities. Facies

F4b and F4d present grainstone to rudstone textures, indicating the presence of local moderate-energy zones in the otherwise calm lagoon [START_REF] Rankey | Contrasts between wave-and tide-dominated oolitic systems: Holocene of Crooked-Acklins Platform, southern Bahamas[END_REF][START_REF] Harris | Mapping bathymetry and depositional facies on Great Bahama Bank[END_REF].

Facies F4b is the only facies in this association displaying high porosity and permeability. As in facies F3b, pores are well connected and from 50 µm to 500 µm in size. Facies F4b forms permeable patches in a mud-dominated environment. Unfortunately, it has not been possible

to model these individual facies due to the scarcity of core descriptions. Facies F4d, while similar in grain type and texture to facies F3d, underwent more intense carbonate cement development. The porosity only exists as intraparticle porosity, as well described around Saint-Martin-de-Bossenay [START_REF] Delmas | Petrophysical Properties of the Middle Jurassic Carbonates in the PICOREF Sector (South Champagne, Paris Basin, France)[END_REF] (Fig. 7H).

Permeabilities can also exhibit a high degree of spatial heterogeneity at all scales, largely depending on the depositional environment and the subsequent diagenetic processes. High space and local dissolution led to complex pore media at the centimetric scale.

3D heterogeneity at the doublet or multi-doublet scale

For the Middle Jurassic limestones of the greater Paris area, no published geological model capturing its fine-heterogeneity (facies, porosity and permeability) at the reservoir scale is currently available, which was highlighted a decade ago by [START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF]. Companies as GEOFLUID or the French Geological Survey (BRGM) maintain an up-to-date hydraulic and thermal model of the aquifer to predict the lifetimes of doublets and the occurrence of the thermal breakthrough, but no heterogeneous geological model capturing the fine heterogeneity of facies, porosity and permeability in the reservoir (XY dimension of the order of hundreds of m and Z dimension of the order of a meter), [START_REF] Ungemach | Sustainable Geothermal Reservoir Management[END_REF]Lopez et al., 2010;[START_REF] Hamm | Mise en oeuvre de la gestion de la ressource géothermique du Dogger de la région Ile-de-France dans le cadre du dispositif d'acquisition et de stockage des données[END_REF][START_REF] Hamm | Gestion du Dogger et corrélation entre niveaux producteurs[END_REF][START_REF] Papachristou | Geothermal resource management -a reservoir simulation approach -the Paris Basin case[END_REF]. The BRGM has built a welldocumented database on the "Geological Data Modeling" (GDM) tool [START_REF] Hamm | Gestion du Dogger et corrélation entre niveaux producteurs[END_REF]. An analysis of the flowmetry well-logs and sedimentological facies distributions along wells was carried out with two important information: (1) the depths of the producing levels with their contribution to the total flow and (2) the depths of the sedimentological facies (lagoonal, otolithic and marly facies, [START_REF] Hamm | Gestion du Dogger et corrélation entre niveaux producteurs[END_REF]. Simple correlations, without variogram analysis, of the producing levels between wells separated by less than 1 km were able to be carried out for the Chelles and Coulommiers doublets, indicating continuities between some producing levels but not for others [START_REF] Hamm | Gestion du Dogger et corrélation entre niveaux producteurs[END_REF]. Geostatistical distributions of available reservoir transmissivities published in the area were made more than 30 years ago and are shown as a map at the scale of 100 km x 100 km [START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF] We also propose very fine heterogeneity in a total of 12.2 million cells capturing heterogeneity of facies, porosity and permeability (XYZ cell dimension of 100 m x 100 m x 3.73 m on average).

As 3D heterogenous geological model and recent geological study at the scale of the reservoir were unavailable, geological risks for total or partial failure existed as mentioned above for the Grigny well [START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF]. Operations along the outer part of the oolitic ramp (facies Local porous patches, smaller in size (1000 m x 850 m x 6 m), are located in mud-dominated limestones, in granular facies (Fig. 21A-C). As permeability is correlated to porosity, permeable bodies also form slightly larger patches (1600 m x 1100 m x 4 m) disconnected from each other (Fig. 21D). It is important to note that this distribution is strongly controlled by facies association distribution in the model. Around Evry, the ooid/bioclast limestones form a wedge prograding westward (Fig. 21B). Total thickness of the ooid/bioclast limestones decreases from east to west, from 50 m to 28 m respectively at the reservoir scale (Fig. 21B). For example, the Grigny 2 well (GGR2) has 3 m of oolitic limestone and 6.5 m of mud-dominated limestone facies (which could also have some local high reservoir quality) in sequence MJXb.

In sequence MJXa, 19 m of oolitic limestones with marl/limestone alternations are found (Fig. 21B). A few kilometers further southeast, a good reservoir example can be observed in the Evry 1 well (GEV1) with 32 m thickness for sequence MJXb and MJXa (Fig. 21B). In the Grigny 2 well (GGR2), the high-quality ooid/bioclast limestones change laterally into marl/limestone alternations, which have low reservoir quality. As no model was available before the installation of the Grigny doublet, a geological risk has been taken to drill GGR2 in this area.

The GGR2 well turned out to be a failure in terms of transmissivity with very thin beds of porous limestones (meter-scale). The 3D zoom extracted from our model clearly shows the thickening of the marl/limestone alternations and thinning of the oolitic limestone reservoir from GGR-1 to GGR-2. The representation given here is very precise, predicting depths of the top of the reservoir and possible permeable levels, typically located 20 m below the caprock in this area. This indicates that models as the one presented here in this study could de-risk future operations and provide a tool to assist investors in deciding on the drilling of future geothermal wells by providing risk criteria for future projects (Fig. 21).

3D seismic is proving to be a very promising tool to strengthen the predictions of geothermal resources [START_REF] Casini | Results of a 3D seismic survey at the Travale (Italy) test site[END_REF][START_REF] Lüschen | 3D seismic survey for geothermal exploration at Unterhaching[END_REF][START_REF] Drouiller | 3D seismic for design and derisking of dual geothermal boreholes in sedimentary sequences and new prospects in the Paris Basin (Adapted methodology using petroleum industry techniques)[END_REF][START_REF] Salaun | High-resolution 3D seismic imaging and refined velocity model building improve the image of a deep geothermal reservoir in the Upper Rhine Graben[END_REF].

Use of 3D seismic data is increasing in the exploration and production of geothermal resources in European basins (Munich, Vienna, …). This kind of data is, however, not available in this study area, but it would be an excellent addition to de-risk the implementation of new doublets as it could help to localize areas of high potential in zones where well data is not present nearby and its general use should be seriously discussed also for the Paris area [START_REF] Drouiller | 3D seismic for design and derisking of dual geothermal boreholes in sedimentary sequences and new prospects in the Paris Basin (Adapted methodology using petroleum industry techniques)[END_REF]. Deep neural networks are also promising to estimate porosity and permeability in the Middle Jurassic geothermal reservoir from old 2D seismic lines acquired in the 1980s and well log data [START_REF] Allo | Characterization of a carbonate geothermal reservoir using rockphysics-guided deep neural networks[END_REF].

The facies property model is very deterministic. This is due to the "Truncated Gaussian with

Trends" algorithm, enabling the possibility to strongly influence the facies repartition both horizontally and vertically inside a zone (zone as a thickness layer in the model, here, depositional sequence) to constrain the property to mimic the "trends" observed in the facies data. However, the porosity and permeability properties are distributed randomly in the volume, such as to respect the variograms, histograms, and hard data issued from the upscaling process. A good way to reduce the uncertainty of porosity and permeability prediction as the distance of a cell increases from the closest cell containing hard data would be to use 3D seismic data, precise gravimetry and recalculated/synthetic well logs or reservoir analogue to populate the grid with usable trends [START_REF] Petrovic | Facies heterogeneities in a ramp carbonate reservoir analogue: a new high-resolution approach for 3D facies modelling[END_REF][START_REF] Arnold | Uncertainty Quantification in Reservoir Prediction: Part 1-Model Realism in History Matching Using Geological Prior Definitions[END_REF][START_REF] Feng | Lithofacies classification of a geothermal reservoir in Denmark and its facies-dependent porosity estimation from seismic inversion[END_REF]. This could enhance the predictability of properties in the volume.

Creating synthetic logs from other well logs could be a way to de-risk future operations, allowing for more upscaled cells in the model. For example, it is possible to create a permeability log by applying porosity-permeability relationships from the Neutron logs [START_REF] Lucia | Carbonate Reservoir Characterization[END_REF]. However, in this study, as porosity-permeability data is relatively heterogeneous along facies association (Fig. 15), only real data were used from well logs as synthetic logs could be a new vector of uncertainty. For new geothermal operations, the NMR logging tool is very promising, giving a continuous record of permeability [START_REF] Wielemaker | Delineating the geothermal structure and flow properties in a sub-horizontal well with the use of wireline and LWD data in a multiphysics approach[END_REF].

Conclusions

To guide the positioning of future geothermal doublets, this paper presents the first highresolution 3D reservoir model of the Middle Jurassic geothermal reservoir in the southeastern greater Paris area (about 40 km x 50 km x 0.2 km). This 3D model, compiling all the available well-log data for the area, details the geological architecture and specifies the 3D reservoir heterogeneity of the late Bajocian to early Callovian limestones. The resulting 3D model is one possible 3D realization showing the geological heterogeneities. The grid of the model including zones, facies, porosity and permeability propagated in the volume, as well as temperature and salinity with their representative maps, are available online as a .rar file. This realization should be used to inspire further modeling and doublet planning. Important information such as the top of the reservoir, facies heterogeneity or possible location of permeable layers will be available. The geostatistical results presented here could be taken into account in future geothermal developments in order to propose other possible realizations and to investigate, for example, the uncertainties. This 3D representation will be a prerequisite for thermo-hydro-dynamic simulations and could provide support to assist decision-making on potentially low-risk areas for future geothermal development in a carbonate reservoir. The primary outcomes of this study are:

1. The creation of a digital database compiled for 168 wells. This database includes 4561 porosity and 3317 permeability values from routine analysis of mini-cores from 40 wells, precisely located in 3D. It is the most complete published database of the Middle Jurassic geothermal reservoir in this part of the basin.

2. The late Bajocian to early Callovian deposits contain 13 facies grouped into four facies associations: (1) clays and marls from lower offshore, (2) marl/limestone alternations from upper offshore, (3) ooid/bioclast limestones from shoreface, and (4) muddominated lagoonal limestones. Many oolitic facies have good reservoir quality (25% of the porosities are higher than 17% and permeabilities higher than 340 mD). The granular facies located in the mud-dominated facies association locally has good reservoir units (25% of the porosities above 12.7% and permeabilities above 252 mD).

Figure 1 .

 1 Figure 1. A-Location of the study area in the Paris Basin on a simplified geological map of the northern part of France with the greater Paris boundary (grey outline: administrative Ile-de-France Region) and the model area (red outline). Note the location of two cores studied: Vert-Saint-Denis-1 (VSD-1) located in the model area near Melun and Sergine-1 (SER-1), out of the model area. B-Location of the wells on a temperature map at the top of the Bathonian aquifer. Blue ellipses: geothermal doublets in the aquifer used in the study; white dots: wells used in this study; black dots: other wells present in the region. The following abbreviations are used: GCDN-1 (Geothermal Cachan Doublet North-1), CBG-1D (Croissy-Beaubourg-1D), BTS001

Figure 2 .

 2 Figure 2. A-Palaeogeographic map of France and part of England during the Late Bathonian with the location of the study area, adapted from[START_REF] Hendry | Geochemical trends and palaeohydrological significance of shallow burial calcite and ankerite cements in Middle Jurassic strata on the East Midlands Shelf (onshore UK)[END_REF] and[START_REF] Brigaud | The origin and timing of multiphase cementation in carbonates: Impact of regional scale geodynamic events on the Middle Jurassic Limestones diagenesis (Paris Basin, France)[END_REF]. The same color code for facies association types has been used

Figure 3 .

 3 Figure 3. Biostratigraphic and well-log calibration between Champigny-sur-Marne (CSM 1D) well (see Fig. 2B for location in the study area) and the reference well of Berthuis 001 (BTS001) in the Coulommes oil field. The biostratigraphy of this oil field,

Figure 4A .

 4A Figure4A. Section of the Vert-Saint-Denis-1 core located in the Paris Basin (VSD-1 located in Fig.1). This core section displays the sedimentological log with the constituent allochems and texture, well-logs laterolog deep (LLD), neutron (NEUT) and

Figure 4B .

 4B Figure 4B. Legend of the core section presented in Fig. 4A.

Figure 5 .

 5 Figure 5. Distribution of the Bathonian facies associations on a theoretical downdip profile. This synthetic profile shows the succession of the facies across the platform, depending on bathymetry.

Figure 6 .

 6 Figure 6. A-Core picture of clays (facies F1b, lower offshore from the Vert-Le-Gr 4 well (VLG-4, 1525 m)). B-Core picture of

Figure 8

 8 Figure 8. A-Picture of the VSD1-1712,9m sample dominated by F3d facies (Rudstone to grainstone with bivalves and intraclasts) B-3D porosity representative volume elements (RVE) visualized using µCT scan at IFPEN on a mini-core sample from the VSD1 core (1712,9m sample in measured depth). Each individual pore (not connected to others) is depicted in a separate color. C-Thin section located in the VSD1 core, 20 cm beneath the sample, analyzed for 3D porosity (1712,9m sample). D-Graph of the porosity distribution in each of the three directions in the VSD1-1712,9m sample. Each color corresponds to a separate pore space.

Figure 9 .

 9 Figure9. A-Stratigraphic cross-section from SW to NE, between the Etrechy-1D well and the Berthuis-1 well. This cross-section is located in fig.1B. Each well is represented with depth, gamma ray values, and facies. Berthuis-1 is not in the model zone but is a reference well for the stratigraphic sequence in the Paris Basin[START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonatés: le jurassique moyen de l'Angleterre à la Méditerranée[END_REF]. B-Same stratigraphic cross-section as A, extracted from the 3D model in Petrel© (except Berthuis-1 well). This section represents the final facies model with well locations and stratigraphic surfaces (vertical exaggeration is set to 50).

Figure 10 .

 10 Figure 10. A-Stratigraphic cross-section from NW to SE, between Cachan Nord-1 well and BRIE-119 well. This cross-section is located in fig. 3B. Each well is represented with depth, gamma ray values and facies. B-Same stratigraphic cross-section as A, extracted from the 3D model in Petrel©. This section represents the final facies model with the location of wells and stratigraphic surfaces (vertical exaggeration is set to 50). TVD: True Vertical Depth.Below the Bt2 surface, the facies are mainly marl-limestone alternations (facies association

Figure 11 .

 11 Figure 11. Facies association proportion in well logs, in upscaled cells along each well and in the entire model (12.2 million cells).

Figure 12 .

 12 Figure 12. A-Facies map at the top of Bt5 maximum regressive surface with the thickness (black line) of the sequence MJXb deposits. B-Facies map at the top of Bt4 maximum regressive surface with the thickness of sequence MJXa deposits. C-View of the entire 3D model with cross-sections generated by Petrel© showing facies association architecture. Facies association is upscaled in each well. Only the bottom three layers (of 52 layers) are displayed.

Figure 13 .

 13 Figure 13. A-Mini-core porosity sorted by facies association. B-Mini-core permeability sorted by facies association. C-Minicore and NPHI log porosity upscaled data sorted by facies association. Note that values differ widely between A and C in the marl/limestone alternations and clay facies associations. This is a known phenomenon of uncalibrated Neutron Porosity (NPHI)

Figure 14 .

 14 Figure 14. A-Vertical (Pv) and horizontal (Ph) porosity distribution across each facies association. B-Porosity distribution in vertical and horizontal directions. C-Vertical (Kv) and horizontal (Kh) permeability distribution across each facies association (excluding clays, where no data could be found). D-Permeability distribution in vertical and horizontal directions.

Figure 15 .

 15 Figure 15. Porosity vs. permeability mini-core data and borehole (GCAH-2) well-logs (TCMR and KSDR) sorted by facies association (A for clays and marls, B for marl-limestone alternations, C for ooid/bioclast limestones and D for peloid to oncoidal limestones). Linear regression lines are plotted on the graphs. Regressive curves are computed only with mini-core porosity and permeability (not with NMR). Values of permeability equal to 0.1 mD on mini-core data are linked to physical detection limits of measuring devices. Φ are expressed in %.

Figure 16

 16 Figure 16. A-View of the porosity model from the southeast with vertical cross-sections showing facies variations. Stratigraphic surfaces are arrowed. B-View of the permeability model from the southeast with vertical cross-sections showing porosity variations.

Figure 17

 17 Figure 17. A-Variogram map of the porosity across the model showing the regional anisotropy of the property. B-Distribution of the raw porosity data in the ooid/bioclastic facies association. C-Example of the data analysis in Petrel© of a property (vertical semi-variance of the porosity in the lagoon) across the zone MJXb of the model.

Figure 18 .

 18 Figure 18. Porosity and permeability distribution in the different facies associations compared between the final model, upscaled cells in boreholes and original data in mini-cores given as a percentage of the total number of cells or measured data. The data distribution is displayed in darker colours for each facies, and the distribution of values in all cells of the model, independent of facies, is given in a lighter shade. Mini-core data are used as a reference dataset to match the property distribution in the model. The graphs are separated into 40 columns.

Figure 19

 19 Figure 19. A-Average porosity map of sequences MJVIIIb to MJXI. B-Average permeability map of sequences MJVIIIb to MJXI. C-Average porosity map of sequence MJXb. D-Average permeability map of sequence MJXb. E-Average porosity map of sequence MJXa. F-Average permeability map of sequence MJXa. Dashed lines delimit areas discussed in the main text.

Figure 20 .

 20 Figure 20. A-Temperature map at the top of the late Bathonian (Bt5) created from well data. B-Temperature gradient map at the top of the late Bathonian (Bt5) created from well data. C-Salinity map at the top of the late Bathonian (Bt5) created from well data D-Transmissivity map of reservoir between Sequences MJXI to MJVIIIb. E-Transmissivity map of Sequence MJXb. F-Transmissivity map of Sequence MJXa.

  mapped in the study zone (Fig.20D-F). The eastward deepening of the top of the reservoir from a vertical depth of 1460 m in Paris (temperature 58 °C) to 1660 m (temperature 82 °C) around Torcy is not sufficient to explain the 20 °C increase in temperature. It could only explain a drop of 4 °C to 8°C. Furthermore, the top of the reservoir is located at less than 1400 m depth near Arpajon (temperature 68 °C), 10°C hotter than under Paris (58°C at 1460 m),

(

  locally above 50 D.m in area 3, Fig 20D), providing trends and a better understanding of the variations at the regional scale. The focus on the sequences MJXa and MJXb highlights the best areas of geothermal energy in the southern part of the greater Paris area (Fig 20E-F).

  porosity and permeability values are controlled by(1) a high-energy environment, (2) low carbonate cement development, and (3) dissolution. A porous media micro-CT scan reveals marked pore space heterogeneity, as demonstrated for facies F3d(Fig 8). Facies F3d has a grainstone to rudstone texture and is composed of common intraclasts, gastropods, bivalves, brachiopods, and echinoderm fragments. Large pore spaces up to the millimetric scale and thin layers with millimetric bivalve and brachiopod shells probably induced strong horizontal permeability anisotropy (from 70 mD to 145 mD) and relatively high vertical permeability (135 mD). Heterogeneous millimeter-size blocky calcite crystals cementing the intergranular

FA3)

  are very risky(around Evry, Figs. 20 and 21). At the multi-doublet scale, our model allows us to view 3D heterogeneity in the grid. Eight doublets in the area around Evry provide a highresolution picture of facies, porosity, and permeability distribution. At this scale, the model presented in this study clearly shows that porous areas form patches of about 1500 m x 1000 m and thickness of 6 m on average in ooid/bioclast limestones, separated by lowporosity (<10%) units (Fig.21A-C). Stratigraphically, these porous patches are mainly located just below maximum regressive surfaces Bt5 and Bt4, in sequences MJXa and MJXb (Fig.21C).

Figure 21 .

 21 Figure 21. Close look at the zone of the Grigny Bathonian reservoir highlighted in Fig 3. A-West-east stratigraphic crosssection through the modeled area, crossing near cities of Arpajon, Bre tigny, Evry, Auvernaux and Maincy. The two main reservoir sequences MJXa and MJXb prograde westward. B-Facies distribution. C-Porosity distribution. D-Permeability distribution. The red rectangle highlights the reservoir thinning discussed in the main text.

  3. A high-resolution sequence stratigraphy scheme was constructed on the ammonite biozone scale, with the identification of 12 third-order depositional sequences (MJXIVb to MJVIIIa), that can easily be integrated into existing stratigraphic schemes on the basin scale and more generally on the scale of the western European basins. Sequences MJXa and MJXb (hodsoni and orbis Zone) form the two targeting sequences for geothermal development. 4. Facies association distribution defined in the 168 wells are upscaled to populate the grid using the Truncated Gaussian with Trends (TGwT) algorithm, applied sequence by sequence. The facies model presents facies heterogeneity in 53 layers to maintain a very good vertical resolution (average of 3.73 m). This 3D facies model locates the ooid/bioclast limestone wedge, prograding from east to west. 5. Porosity and permeability along wells are upscaled to the model layering by calculating the arithmetic average for the porosity and the geometric average for the permeability. Of the 12.2 million cells of the entire grid, 6384 cells (XYZ cell dimension of 100 m x 100 m x 3.73 m on average) have upscaled facies, 1686 cells have upscaled porosity and 469 cells have upscaled permeability. 6. The temperature of geothermal water ranges from 51 °C to 79.3 °C. A temperature gradient map shows heterogenous geothermal gradients, from high values (37 to 41 °C/km) in an area between Etrechy in the southwest and Croissy-Beaubourg in the northeast, to low gradients of about 27 °C/km, around Paris itself. The area presenting the highest temperature gradients also presents average good petrophysical reservoir qualities and transmissivities. 7. Geostatistics show the occurrence of permeable patches of 1600 m x 1100 m (XY dimensions) and a thickness of 4 m (Z dimension) elongated perpendicular to the depositional slope in the ooid/bioclast limestone. In the peloidal-oncoidal limestone, permeable patches consist of porous granular facies with a size of 850 m x 750 m x vertical thickness of 6 m (XYZ dimensions).

Table 1 :

 1 Facies table for the facies associations with sorted samples from cores. VR = Very rare < 5% R = Rare < 10%; C = Common 10 -20%;

	Depositional
	environments and corresponding facies and/or facies	Lithofacies	Non-bioclastic componen ts	Bioclastic components	Energy and Depositional environment	Texture (Dunham, 1962)	Sortin size g and grain	Porosity	Diagenesis	Reservoir potential
	association
					Very low
	Lower Offshore,	F1a -Marls	Silts, clay	-	energy, Lower	Mudstone	-	No effective porosity	-	Non-existent
	Palaeodepth >40				Offshore
	m (Facies
	association FA1 -				Very low
	Clays, marls)	F1b -Clays	Clay	-	energy, Lower	Mudstone	-	No effective porosity	-	None-existent
					Offshore
				Foraminifer
	Upper Offshore, Palaeodepth 20-40 m (Facies association FA2 -Alternating	F2a -Alternating marls and peloidal wackestone to packstone	Peloids (A), Quartz (C), Clay (C)	a (F), Bivalves (F), Echinoderm fragments (C), Crinoids (R), Gastropods	Low to moderate energy, Upper Offshore	Packstone, Wackeston e	Poorl d y sorte	No matrix macroporosit y in micritic y / very low microporosit	Micritization, pyrite early stylolithization, microsparitic cement,	Unusable
	marls and mud-			(R)
	dominated limestones)	F2b -Siliciclastic bioturbated packstone	Peloids (A), Quartz (C), Clay (C)	Foraminifer a (F), Bivalves (F), Echinoderm fragments	Low to moderate energy, Upper Offshore	Packstone	Poorl y sorte d	Sparse intraparticle macroporosit microporosit y,	Micritization, cement, early microsparitic	Unusable

F = Frequent 20 -50%; A = Abundant > 50% 4.2 Facies architecture and sequence stratigraphy

Table 2 :

 2 Anisotropy values used in the model, obtained from the Data Analysis. All measurements are in meters. (N.P. = Not Present in the zone, L.D. = Low Data, meaning no clear statistical distribution is identified).

  . The dataset of Hamm et al. (2015) is not publicly available, including geolocalization of wells, facies along wells, producing level and flow rate, for a reopening on GDM or other geomodeler such as Petrel©. Here, we propose a strong update of the geological data of the geothermal reservoir, providing an openaccess 3D digital model (XYZ dimension of the entire model of about 40 km x 50 km x 0.2 km).

Acknowledgements

This work is supported by the UPGEO "UPscaling and heat simulations for improving the efficiency of deep GEOthermal energy" program funded by the Agence Nationale de la Recherche [ANR-19-CE05-0032-01]. We would like to dedicate this work to Andro Mikelic, who left us in 2020, and who was a driving force behind the construction of the UPGEO project and who more generally developed many research activities in mathematics on upscaling and fluid flow in porous media. This work is also supported by the "Geothermal water circulation in

issued from the reports or well logs are present as hard data and can be observed in the corresponding well path. Considering the size of the modeled area, a variance map of the average porosity of all sequences flattened to the Bt5 surface is produced using 40 lags of 1 km (Fig. 17). Since this model is large, variogram maps cannot be used to determine the anisotropy accurately, thickness, length, and width of possible carbonate geobodies (porous layers of oolitic dunes, emerged islands, or less-cemented units; Brigaud et al., 2014). For the level of precision aimed at, variogram mapping is only used to detect the main direction of geobodies in a wide area (40 km in X, 40 km in Y). A main direction of N110° is detected (Fig. 16A). This is the main direction of the oolitic barriers around Paris at the end of the Bathonian, and this direction is, in general, perpendicular to the slope of the ramp [START_REF] Enay | Synthèse paléogéographique du Jurassique français[END_REF][START_REF] Mégnien | Synthèse géologique du bassin de Paris[END_REF][START_REF] Purser | Sédimentation et diagenèse des carbonates néritiques récents: les éléments de la sédimentation et de la diagenèse[END_REF][START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger, Bassin parisien[END_REF]. To better determine the anisotropy values, an individual variogram is calculated for each facies association using porosity and permeability distributions from well-logs (Fig. 17C, Table 2). The anisotropy of the different facies associations is determined first across the entire model and is then refined at the scale of the two targeted, high-quality reservoirs deposited during sequences MJXb and MJXa (Fig. 17C, Table 2). This step enhances the accuracy of the petrophysical data propagation across these two important sequences, especially for high reservoir-quality facies associations (1) ooid limestones and (2) lagoonal facies (Table 2). University Paris-Saclay benefits from the Schlumberger Software Donation Program. We are grateful to Schlumberger for the academic license for Petrel© software and for very valuable support in its utilization. We also thank Anaëlle Guillevic, Maë-Berinda Babadoudou, Maoli Compay (bachelor students), Vivien Folliot, Anne-Sophie Canivet, Myriam Mansour, Codjo Thomas Florent Essou (master students) for the digitalization of wells using Neuralog.

We also thank the four anonymous reviewers who helped with their comments to increase the clarity of the manuscript.

Data availability

Datasets of the 3D model are available as a .rar file which can be found at https://data.mendeley.com/datasets/kd6dyfwg5f/1 (doi:10.17632/kd6dyfwg5f.1), an opensource online data repository hosted at Mendeley Data. This file can be read only with Schlumberger software Petrel©, version 2019.2 and later. It contains an extracted grid of the model including zones, facies, porosity and permeability propagated in the volume, as well as temperature and salinity with their representative maps. Supplementary Table A summarizes the wells manager available in Petrel©. It contains the well names, positions, maximum depth, type of well, origin of the data and a summary of the available well logs. Supplementary Table B, an Excel sheet, resumes all the petrophysical data available from mini-cores and integrated in the model, as well as the temperatures and salinities.

Appendices