Hadrien Thomas

Benjamin Brigaud

Thomas Blaise

Bertrand Saint-Bezar

Elodie Zordan

Hermann Zeyen

Simon Andrieu

Benoît Vincent

Hugo Chirol

Éric Portier

Eric Portier

Emmanuel Mouche

Contribution of drone photogrammetry to 3D outcrop modeling of facies, porosity, and permeability heterogeneities in carbonate reservoirs (Paris Basin, Middle Jurassic)

Keywords: Carbonates, Digital outcrop modeling, Drone, Photogrammetry, Facies modeling, Reservoir analogue, Pore type

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In carbonate rocks, the interplay between depositional environments and diagenetic processes can lead to significant variations in reservoir quality making it difficult to predict fluid pathways within and across geobodies in three dimensions (3D) (Ehrenberg et al., 2006a(Ehrenberg et al., , 2006b;;[START_REF] Borgomano | Stratigraphic well correlations for 3-D static modeling of carbonate reservoirs[END_REF][START_REF] Ronchi | Depositional setting and diagenetic processes and their impact on the reservoir quality in the late Visean-Bashkirian Kashagan carbonate platform (Pre-Caspian Basin, Kazakhstan)[END_REF]. Well-described field analogues can help understanding the geometry of subsurface reservoirs [START_REF] Shekhar | Modelling and simulation of a Jurassic carbonate ramp outcrop, Amellago, High Atlas Mountains, Morocco[END_REF]. For this purpose, limestone outcrops have been extensively studied, along with numerical modeling using different software packages (e.g. Ehrenberg et al., 2006aEhrenberg et al., , 2006b;;[START_REF] Rankey | Form, function and feedbacks in a tidally dominated ooid shoal, Bahamas[END_REF][START_REF] Harris | Delineating and quantifying depositional facies patterns in carbonate reservoirs: Insight from modern analogs[END_REF][START_REF] Palermo | Three-dimensional facies modeling of carbonate sand bodies: Outcrop analog study in an epicontinental basin (Triassic, southwest Germany)[END_REF][START_REF] Amour | Capturing and modelling metre-scale spatial facies heterogeneity in a Jurassic ramp setting (Central High Atlas, Morocco): Metre-scale facies modelling[END_REF][START_REF] Gasparrini | A multidisciplinary modeling approach to assess facies-dolomitizationporosity interdependence in a lower cretaceous platform (northern Spain)[END_REF][START_REF] Tomassetti | Modeling lateral facies heterogeneity of an upper Oligocene carbonate ramp (Salento, southern Italy)[END_REF]. [START_REF] Amour | Capturing and modelling metre-scale spatial facies heterogeneity in a Jurassic ramp setting (Central High Atlas, Morocco): Metre-scale facies modelling[END_REF][START_REF] Amour | Outcrop analog for an oolitic carbonate ramp reservoir: A scale-dependent geologic modeling approach based on stratigraphic hierarchy[END_REF] in 3D to better understand the meter-scale facies architecture of the deposits and their heterogeneity [START_REF] Lipinski | Reservoir analog model for oolite-microbialite sequences, Miocene terminal carbonate complex, Spain[END_REF][START_REF] Benson | Workflows for incorporating stratigraphic and diagenetic relationships into a reservoir-analogue model from outcrops of Miocene carbonates in SE Spain[END_REF][START_REF] San Miguel | Facies heterogeneity of a Kimmeridgian Carbonate Ramp (Jabaloyas, eastern Spain): a combined outcrop and 3D geomodelling analysis[END_REF][START_REF] Tomassetti | Modeling lateral facies heterogeneity of an upper Oligocene carbonate ramp (Salento, southern Italy)[END_REF][START_REF] Brandano | Facies heterogeneities and 3D porosity modelling in an Oligocène (Upper Chattian) carbonate ramp, Salento peninsula, southern Italy[END_REF]. The development of drone photogrammetry in recent years has eased the quantitative analysis of geological outcrops, even along cliff faces [START_REF] Schmitz | From 3D photogrammetric outcrop models to reservoir models: an integrated modelling workflow[END_REF][START_REF] Cabello | From conventional outcrop datasets and digital outcrop models to flow simulation in the Pont de Montanyana point-bar deposits (Ypresian, Southern Pyrenees)[END_REF][START_REF] Rubi | Processes and typology in Gilberttype delta bottomset deposits based on outcrop examples in the Corinth Rift[END_REF][START_REF] Nesbit | 3-D stratigraphic mapping using a digital outcrop model derived from UAV images and structure-from-motion photogrammetry[END_REF][START_REF] Siddiqui | High resolution facies architecture and digital outcrop modeling of the Sandakan 53 formation sandstone reservoir, Borneo: Implications for reservoir characterization and flow simulation[END_REF]. However, most of these digital meter-scale resolution captures were applied to sandstone outcrops (e.g., [START_REF] Hayes | Analyzing Dune Foreset Cyclicity In Outcrop With Photogrammetry[END_REF][START_REF] Cabello | From conventional outcrop datasets and digital outcrop models to flow simulation in the Pont de Montanyana point-bar deposits (Ypresian, Southern Pyrenees)[END_REF][START_REF] Rubi | Processes and typology in Gilberttype delta bottomset deposits based on outcrop examples in the Corinth Rift[END_REF]. While geological photogrammetric models exist, complete static geological models including distribution of properties like porosity and permeability in complex sedimentary systems have not yet been created in limestones using these technologies [START_REF] Madjid | Drones in carbonate geology: Opportunities and challenges, and application in diagenetic dolomite geobody mapping[END_REF][START_REF] Marques | Virtual and digital outcrops in the petroleum industry: A systematic review[END_REF].

As prototype of such a study, we used the outcropping clinoform formation of the Massangis quarry located in the Middle Jurassic carbonates of the Paris Basin (Fig. 1). In carbonate rocks, clinoforms are inclined depositional surfaces corresponding to paleo-bathymetric profiles, which are well expressed within progradational sedimentary systems and sometimes form prolific reservoir geobodies (e.g. [START_REF] Plink-Björklund | Shallow-water deltaic clinoforms and process regime[END_REF]. A clear understanding of the factors controlling the vertical and lateral heterogeneity of porosity and permeability within clinoforms is fundamental for improving conceptual reservoir models of these sedimentary features [START_REF] Bryn | Exploring the reservoir potential of Lower Cretaceous Clinoforms in the Fingerdjupet Subbasin, Norwegian Barents Sea[END_REF][START_REF] Steel | What is the topset of a shelf-margin prism?[END_REF]. Detailed static architectural models of such clinoforms are rare and still require further research.

The choice of this site is motivated by the fact that the Basin is well studied, mainly because of its importance as petroleum and geothermal reservoirs and for CO 2 storage [START_REF] Mougenot | Seismic imaging of a carbonate reservoir; the Dogger of the Villeperdue oil field, Paris Basin, France[END_REF]Brosse et al., 2010;[START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF][START_REF] Brigaud | Characterization and origin of permeability-porosity heterogeneity in shallow-marine carbonates: From core scale to 3D reservoir dimension (Middle Jurassic, Paris Basin, France)[END_REF][START_REF] Regnet | Influence of microporosity distribution on the mechanical behavior of oolithic carbonate rocks[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]. In particular, Massangis quarry in Burgundy has been studied as a reservoir analogue of the Oolite Blanche Formation [START_REF] Purser | Dedolomite Porosity and Reservoir Properties of Middle Jurassic Carbonates in the Paris Basin, France[END_REF][START_REF] Casteleyn | Interrelations of the petrophysical, sedimentological and microstructural properties of the Oolithe Blanche Formation (Bathonian, saline aquifer of the Paris Basin)[END_REF][START_REF] Casteleyn | An integrated study of the petrophysical properties of carbonate rocks from the "Oolithe Blanche" formation in the Paris Basin[END_REF][START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF][START_REF] Regnet | Influence of microporosity distribution on the mechanical behavior of oolithic carbonate rocks[END_REF], a major geothermal reservoir in the Paris Basin [START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France[END_REF]. Although the sedimentological, diagenetic, and petrophysical frameworks of this reservoir analogue are now well documented [START_REF] Regnet | Influence of microporosity distribution on the mechanical behavior of oolithic carbonate rocks[END_REF], no high-resolution 3D reservoir model is available until now for assessing the origin of reservoir heterogeneity. The main objective of this study is to integrate facies and petrographic descriptions, drone photogrammetry, and structural and geostatistical modeling to understand the facies heterogeneity of a carbonate outcrop dominated by clinoforms in the framework of a newly proposed depositional model. The resolution of photogrammetric models provides in sub-meter-scale (X, Y dimensions) to sub-decimetric (Z dimension) observations of facies trends and geobody geometry. We combine the photogrammetric observations with facies interpretation at different scales, from thin-section (micrometers) to core-scale, and to reservoirscale. These quantitative data and interpretations are then integrated into a reservoir model of a quarry of the Middle Jurassic carbonates of the Paris Basin using Petrel®.

Sedimentological, stratigraphic, and petrophysical frameworks

Massangis quarry is located in the south-eastern part of the Paris Basin, about 150 km south-east of the main geothermal reservoir currently in production around Paris (Fig. 1). The quarry modeled in this study is composed of two production sites: the Vaurion site, the main part of the quarry, and the Malpierre site, north of the D311 road (Fig. 1C). The model area, comprising both sites, is contained in a rectangle 900 m long and 500 m wide. The quarry contains limestones of Mid to Late Bathonian age. During the Bathonian, the Paris Basin was an epicontinental sea located in subtropical latitudes at about 30° to 25°N and open to the Atlantic, Tethys, and Northern Oceans [START_REF] Contini | Evolution paléogéographique de la France au Jurassique moyen[END_REF][START_REF] Thierry | Middle Callovian, map 9[END_REF].

Bathonian sedimentation across the Paris Basin consisted of four main types of facies association and corresponding depositional environments (Fig. 2): protected lagoon carbonates, high-energy carbonates (shoreface), low-energy carbonates (upper offshore), and lower offshore marls [START_REF] Purser | Tidal sediments and their evolution in the Bathonian carbonates of Burgundy, France[END_REF][START_REF] Purser | Plates-formes carbonatées-exemple du Jurassique moyen du bassin de Paris[END_REF]. The protected lagoon environments are found in the eastern part of the Paris Basin (Fig. 2). They consist of mudstones with oncoids and benthic foraminifera, wackestones/packstones with oncoids, and peloidal grainstones (Comblanchien Limestones, Fig. 2). This lagoon was adjacent to high-energy shallow carbonate sedimentation domains characterized by grainstone facies composed of ooids, peloids, bivalve, bryozoan, and gastropod debris (Oolite Blanche Formation, Fig. 2). Shallow carbonate sedimentation also occurred along the Armorican Massif in the western part of the basin (Armorican Platform, Fig. 2). An open marine domain characterized as lower offshore was located in the central part of the Paris Basin, consisting mainly of shales and marls (Fig. 2). A storm-dominated upper offshore domain occurred between the shallow carbonate platform and the deeper zone (Fig.

2).

Based on previous studies of the Massangis quarry, three main parts have been recognized [START_REF] Loreau | Géométrie et couleur des couches dans l'oolithe blanche de Massangis (Yonne)[END_REF][START_REF] Casteleyn | Interrelations of the petrophysical, sedimentological and microstructural properties of the Oolithe Blanche Formation (Bathonian, saline aquifer of the Paris Basin)[END_REF][START_REF] Casteleyn | An integrated study of the petrophysical properties of carbonate rocks from the "Oolithe Blanche" formation in the Paris Basin[END_REF][START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF]. The lower part is mainly composed of tidal deposits. The middle part is made up of limestones displaying large clinoforms up to 15 m-thick representing the main part of the rocks quarried commercially. These two parts of the quarry belong to the Oolithe Blanche Formation. The upper part consists of lagoonal deposits, named Comblanchien Limestones, capping the large prograding geobodies. [START_REF] Loreau | Géométrie et couleur des couches dans l'oolithe blanche de Massangis (Yonne)[END_REF] identified two beds including characteristic brachiopod species constraining stratigraphic framework: Burmirhynchia icaunensis at the top of the lower part and Burmirhynchia thierachensis at the base of the upper part. Burnirhynchia icaunensis is typical of the ammonite hodsoni Zone and Burnirynchia thierachensis of the orbis Zone (Fig. 2, [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]. The foraminifer Meyendorfina bathonica in the upper part is an additional marker, although of less importance, for the Late Bathonian age of the Comblanchien Limestones. The quarry deposits are thus dated to the Late Bathonian, hodsoni to orbis Zones (Fig. 2A), spanning from 166.1 Ma to 165.1 Ma.

Despite several detailed sedimentological studies, the genetic origin of the 15 m-thick prograding clinoform complex unit is still debated. [START_REF] Casteleyn | Interrelations of the petrophysical, sedimentological and microstructural properties of the Oolithe Blanche Formation (Bathonian, saline aquifer of the Paris Basin)[END_REF] and [START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF] [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](Jacquin et al., , 1998;;[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]Gaumet et al., 1996;[START_REF] Brigaud | Characterization and origin of permeability-porosity heterogeneity in shallow-marine carbonates: From core scale to 3D reservoir dimension (Middle Jurassic, Paris Basin, France)[END_REF]. Depositional sequences are composed of transgressive and regressive systems tracts bounded by sequence boundaries [START_REF] Catuneanu | Sequence stratigraphy: methodology and nomenclature[END_REF]. The maximum flooding surfaces (MFS) of each sequence, marking the boundary between deepening-upward to shallowing-upward trends, are well dated by the occurrence of ammonite fauna and, more often, specific brachiopod-rich beds 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] Garcia | The validity of two biostratigraphic approaches in sequence stratigraphic correlations: brachiopod zones and marker-beds in the Jurassic[END_REF]. The brachiopod association found in the lower part of Massangis quarry corresponds to the Concinna MFS (hodsoni Zone) and that of the upper part of the quarry marks the Digona MFS (orbis Zone, Fig. 2). The sequence boundaries Bt3 and Bt4 located in this time window are identified in the quarry on the basis of sequence stratigraphic analysis [START_REF] Vincent | Large-scale prograding oolitic geometries in the Middle Jurassic "Oolithe Blanche" Formation[END_REF]. All these surfaces are placed within the standard European Ammonite Zones (Fig. 2, [START_REF] Gradstein | The geologic time scale[END_REF].

9

Petrophysical data from outcrops of the Massangis quarry were acquired during the last decade.

Three types of porosity and permeability data from this quarry were compiled for the present study. [START_REF] Casteleyn | Interrelations of the petrophysical, sedimentological and microstructural properties of the Oolithe Blanche Formation (Bathonian, saline aquifer of the Paris Basin)[END_REF][START_REF] Casteleyn | An integrated study of the petrophysical properties of carbonate rocks from the "Oolithe Blanche" formation in the Paris Basin[END_REF][START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF] presented a total of 36 porosity and 10 permeability measurements. 18 total porosities come from plug samples taken from five large oriented blocks using water saturation within about 20 cubic decimeters each. On the same blocks, the 10 permeabilities were measured using N 2 as fluid. 9 macroporosities (corresponding to pores exceeding 10 µm) were measured on thin-sections using a threshold surface and 9 microporosities were deduced. Additional petrophysical data, acquired by [START_REF] Casteleyn | Interrelations of the petrophysical, sedimentological and microstructural properties of the Oolithe Blanche Formation (Bathonian, saline aquifer of the Paris Basin)[END_REF][START_REF] Casteleyn | An integrated study of the petrophysical properties of carbonate rocks from the "Oolithe Blanche" formation in the Paris Basin[END_REF] and [START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF] in two nearby quarries (Ravières and Bierry-les-Belles-Fontaines), were added: 85 porosities (50 total porosity, 17 macroporosity and 18 microporosity) and 33 permeabilities. Data from Ravières and Bierry-les-Belles-Fontaines are considered to be usable in this model as both of them have the same microfacies and depositional environments as shown by [START_REF] Casteleyn | Interrelations of the petrophysical, sedimentological and microstructural properties of the Oolithe Blanche Formation (Bathonian, saline aquifer of the Paris Basin)[END_REF][START_REF] Casteleyn | An integrated study of the petrophysical properties of carbonate rocks from the "Oolithe Blanche" formation in the Paris Basin[END_REF] and [START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF]. The Bierry-les-Belles-Fontaines quarry is 17 km distant from Massangis at 47°36'42.2" N 4°10'45.3"E. The Ravières quarry is 23 km distant from Massangis at 47°43'41.8"N 4°13'36.2"E. The lateral continuity of the studied platform allows transposing the petrophysical properties from these 2 other quarries to Massangis. A total of 164 petrophysical values (121 porosities and 43 permeabilities) from literature were thus added to populate the reservoir model.

Material and method

Facies analysis was conducted including a description of lithology, texture, bedding pattern, sedimentary structures, and faunal content (Fig. 1). Field observations included sample collection, facies description of all accessible cliffs, logging at seven separate locations in the quarry (sampling points A, B, C, F, G, H and J, Fig. 1C), and drone photogrammetry. Paleodepths were estimated for each facies by examining their sedimentary structures, textures, and allochems.

Macrofacies analysis was supplemented by detailed observation of pore types and diagenetic features on 35 blue resin-impregnated thin sections using the JMicrovision® software (JMicrovision®, [START_REF] Roduit | JMicroVision: Image analysis toolbox for measuring and quantifying components of high-definition images[END_REF]. Petrophysical data from the studies of [START_REF] Casteleyn | Interrelations of the petrophysical, sedimentological and microstructural properties of the Oolithe Blanche Formation (Bathonian, saline aquifer of the Paris Basin)[END_REF][START_REF] Casteleyn | An integrated study of the petrophysical properties of carbonate rocks from the "Oolithe Blanche" formation in the Paris Basin[END_REF] and [START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF] were assigned to specific facies based on the depositional environment described by the authors (total porosity, macroporosity, microporosity, and permeability). Measures that could not be associated with confidence were excluded from the dataset.

As pointed out by [START_REF] Catuneanu | Sequence stratigraphy: methodology and nomenclature[END_REF], surfaces selected as sequence boundaries vary from one sequence stratigraphy approach to another and the choice typically depends on which surfaces are best expressed in each situation. The surfaces described in this study underline a transition between different depositional environments and correspond to maximum regressive surfaces. The surfaces Bt4 and Bt3 (Fig. 2) were therefore used to characterize the evolution of facies in the geological model. Moreover, we did not differentiate between Highstand Systems Tract and Lowstand Systems Tract and they are grouped here within a single Regressive Systems Tract (Embry, 2009;[START_REF] Catuneanu | Sequence stratigraphy: methodology and nomenclature[END_REF].

The photogrammetric model was created from 7900 photos taken with a Mavic Pro 2 (Dji™) drone equipped with a 20 Mega Pixel resolution camera. Image resolution was about 2.3 cm/pixel. Planar overflights were necessary to cover most of the quarry. Vertical flights were used to cover the vertical quarry walls more accurately, allowing closer study of the sedimentary formations. The study area covers 0.4 km 2 (40.04 ha). All pictures were georeferenced with an on-board GPS, gyroscope, and accelerometer. The photogrammetric model was created with Pix4D® software. The software outputs are a point cloud file (.xyz) composed of 67 million points, and a file containing mesh triangles (.obj) at high resolution (752 020 nodes, 1 510 551 triangles). Because all points are represents the orthoprojection of the pictures in the dataset oriented and located on the triangles of the mesh (covered area of point with RGB color code: 0.94 cm² in Pix4D, size of the texture file computed: 1,01 GPixels, size of the texture file used in VRGS® : 268 MPixels, excluding some parts of the photogrammetric model of the quarry) allowing us to reconstruct a 3D Digital Outcrop Model Geological Studio (VRGS®), a software package developed by Hodgetts (2010) (Fig. 3) to interpret the sedimentary structures. The sedimentary and stratigraphic surfaces observed in the DOM were manually digitalized to create polylines exported in polygon (.ply) format (Fig. 3).

Petrel® software (Schlumberger) was used to build the geomodel. The georeferenced point cloud extracted from Pix4D® (.xyz) was imported in Petrel® as a pointset to recreate a DOM, adding the offset used by the software as a transition in Petrel® to precisely locate the point cloud of the quarry. These points were interpolated to generate a topographic surface using the "make surface" tool built-in Petrel®. The geocellular model is 917 m long, 380 m wide, and 49 m-thick. The top and the base were defined by the highest and lowest topographic points, which also correspond to the top and bottom of the sedimentary section. Polylines were imported from VRGS® software as guidelines for creating surfaces and clinoforms in Petrel®. To generate the surfaces, the polylines were converted to point clouds and replicated in a 2D window, perpendicular to the direction of the accretion of the clinoforms (Fig. 3, Supplementary Material Table S1). Surfaces were then modeled using the "Make surface" tool of Petrel® and using the point clouds generated previously as input.

This process has been repeated for the 16 original polylines transformed to surfaces. The modeled zones between the surfaces were set to "conformable", except for a bioturbated surface below the clinoform surfaces used as "base" and the surface on top of the clinoforms used as "erosional". In addition to 7 logs measured at the sampling sites in the quarry, 23 synthetic numerical logs were created in the photogrammetric model to position the sedimentological logs and samples. The cell size of the resulting 3D grid was set to 1 m x 1 m for the XY dimensions. The Z-size depends on the thickness of the sub-zones, averaging 0.88 m. The resulting 3D grid is composed of a total of 60 632 040 cells. The model framework is composed of three zones delimiting the main sequence boundaries. The "Truncated Gaussian with Trends" algorithm was used to accurately model facies on the cliffs of the quarry based on the upscaled facies logs created using the samples and observations on the photogrammetric model. Four petrophysical properties were created (total porosity, macroporosity, microporosity, and permeability). They were modeled individually for each facies using the "Gaussian Random Function Simulation" algorithm. Specific distributions were assigned to each facies, including minimum, maximum, mean and standard deviation. These values were determined using a facies by facies statistical analysis of the petrophysical data. The anisotropy of the properties could not be constrained and has therefore not been used. For closer control of the distribution of the properties related to the facies ("facies bias"), co-kriging relationships were established between total porosity and the other two porosity properties (macroporosity and microporosity). Permeability was co-kriged to total porosity as a trend and macroporosity as a collocated co-kriging property. Readers will need to bear these points in mind for a better understanding of the workflow in the following sections.

Sample observations were made on thin sections and macroporosity acquired by image analysis. Previous petrophysical data

were collected from the literature [START_REF] Casteleyn | Interrelations of the petrophysical, sedimentological and microstructural properties of the Oolithe Blanche Formation (Bathonian, saline aquifer of the Paris Basin)[END_REF][START_REF] Casteleyn | An integrated study of the petrophysical properties of carbonate rocks from the "Oolithe Blanche" formation in the Paris Basin[END_REF][START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF]. Images were processed using Pix4D® and the resulting Digital Outcrop Model interpreted using Virtual Reality Geological Studio (VRGS®). Facies, porosity, and permeability were modeled using Petrel® software. Note the stereonet (upper right) showing the dips and azimuths observed on the beddings of clinoforms (data in the Supplementary material Table S1).

Results

Facies

Ten facies were identified for the entire stratigraphic interval under study, defining three depositional environments (facies associations) ranging from upper offshore to lagoon (Fig. 4).

Observations and descriptions are summarized in Table 1. A synthetic sedimentological log of the quarry illustrates the stratigraphic succession and the depositional environments (Fig. 5). In this section, we describe each facies association in detail from the most distal to the most proximal facies. The same colour code for facies types has been used in subsequent figures..

Clinoforms (facies association FA1)

Description -Four facies are grouped in the "clinoforms" facies association based on dominant texture and allochems (Fig. 4): ( 1) marls (facies F1a), (2) bioturbated peloidal crinoidal packstones (facies F1b), (3) peloidal packstones with crinoids (facies F1c), and (4) ooids and peloidal grainstones with crinoids (facies F1d). Facies F1a is composed of marls, containing few bioclasts, principally crinoids and echinoderm fragments (Fig. 4, Table 1, Fig. 5, Fig. 6A). In the quarry, marls appear to be very rare and have been observed only once at the base of some clinoforms. Marls were deposited at the bottom of the clinoforms and subhorizontal planar laminations are observed (Fig. 6B). Facies F1b is a bioturbated peloidal crinoidal packstone (Fig. 4, Fig. 5, Fig. 6 C&E, Table 1). Some layers contain bivalves and crinoids. Intense bioturbation is also observed in this facies with skolithos and cruziana ichonofacies (Fig. 6E). It is located at the base of the clinoforms. Pores of 50 to 100 µm size form the macroporosity and have a geometric rhomb shape created by diagenetic dedolomitization (Fig. 6D).

This pore type may be characterized as moldic. Facies F1c is a peloidal packstone with crinoids (Fig. 4, Table 1, Fig. 6F). The main difference with facies F1b is the absence of bioturbations. It also differs from facies F1b in that it is found higher up in clinoform sets (Fig. 6G). Pores are 50-400 µm geometric rhombs and are formed by the same dedolomitization process as F1b. Facies F1d is an ooid and peloidal grainstone with crinoids (Apiocrinus sp., [START_REF] Loreau | Géométrie et couleur des couches dans l'oolithe blanche de Massangis (Yonne)[END_REF] and bivalves (Fig. 5, Fig. 6G&H, Table 1). Bioclast alignments underline the bedding. Blocky calcite and syntaxial cement fill the interparticle space, leaving only microporosity within the peloids and ooids.

Interpretation -Marls (facies F1a) and bioturbated peloidal crinoidal packstones (facies F1b) are interpreted as having been deposited in an upper offshore environment below the fair-weather wave base, as indicated by their muddy textures, bioturbations (skolithos and cruziana ichonofacies), and their basal position within clinobeds. As peloidal packstones (F1c) were located higher up in the clinobed sets, they probably formed in upper offshore to lower shoreface environments. Ooid and peloidal grainstones (F1d) formed in shoreface environments, as indicated by their high-energy granular texture and their position high up in the clinoforms.

Wave to tide-dominated ooid grainstones (facies association FA2)

Description -This facies association encompasses several grain-supported facies mainly composed of ooids and crinoids (Fig. 4, Fig. 5). Facies F2a is an ooid grainstone containing crinoids (Apiocrinus sp.), bivalves, brachiopods, and oncoids. Two perforated surfaces and one bioturbated surface (Thalassinoides) have been observed in this facies (Fig. 7A-D, Table 1). Sigmoid and planar crossbedding in marine carbonate sand medium dunes measuring several tens of centimeters is common (Fig. 5). Porosity is widespread due to dedolomitization preferentially affecting the sedimentary infill of perforations in hardgrounds. Facies F2b is an ooid grainstone with crinoids, gastropods, bryozoans, peloids, and herring-bone stratification (Fig. 7A,B,E, Table 1). Observed porosity corresponds mainly to microporosity within ooids (Fig 7F). No bioturbations or perforations are observed in facies F2b (Fig. 7B,E). Facies F2c is a bioclastic grainstone to rudstone, with concentric ooids, peloids, bivalves, echinoderms, gastropods, brachiopods and rare bryozoans, serpulids, and coral fragments (

), and bioturbations (sample MAS-4). Perforations are the only porosity in the samples; they may cause vugs, increasing facies porosity. E-Photograph of facies F2b with cross-stratification and herring-bones (arrows). F-Facies F2b ooid grainstone with bryozoans (Bryo.), echinoderms (Ech), and bivalves (Bi.), sample MAS-BAS-1. G-Photograph of facies F2c with perforations of the hardground and stylolites. H-Facies F2c bioclastic

grainstone to rudstone (sample MAS-1) with peloids (Pel.), bivalves (Biv.), echinoderms (Ech.), superficial ooids (sup. Ooid) and well-developed blocky calcite cement. Refer to figure 4 for colour code of facies types.

Lagoonal environments (facies association FA3)

Description -Facies association FA3 (Fig. 4) encompasses three facies: (1) peloidal grainstone with isopachous cements (facies F3a, Fig. 5, Fig. 8A,B), (2) ooid and oncoidal grainstone (facies F3b, Fig. 5, Fig. 8C, D), and (3) oncoidal packstone to floatstone (facies F3c, Fig. 5, Fig. 8C,E). Facies F3a is composed of ooids, bioclasts, and foraminifera Meyendorfina bathonica (Fig. 5, Fig. 8A,B, Table 1). It is located at the top of the quarry (Fig. 5, Fig. 8A,B, Table 1). Facies F3b (Fig. 5, Table 1) is rich in ooids and contains oncoids, peloids, and bioclasts (foraminifers e.g. Meyendorfina bathonica, bivalves, coral fragments). This facies also displays cross-stratification (Fig. 5) and microporosity is dominant. Facies F3c has a mud-supported texture (packstone to floatstone) dominated by centimeter-size oncoids with bioclasts (bivalves, foraminifers Meyendorfina bathonica, miliolids) (Fig. 5, Fig. 8C,E, Table 1).

Interpretation -The dominance of mud-supported textures and planar stratification both indicate low-energy conditions. The presence of peloids and very well-preserved gastropods is evidence of a sheltered lagoonal environment (Fig. 4). The scattered grain-supported layers are either local tidal deposits (channels, bars) or storm wash-over deposits.

Stratigraphic framework

Three depositional sequences are identified along the studied stratigraphic interval (Fig. 5).

Considering a time range of about 1 My for the interval under study based on [START_REF] Gradstein | The geologic time scale[END_REF], the average duration for each cycle is 0.3 My. Only the upper part of the regressive systems tract of the lowest sequence is observed (Sequence I, Fig. 5). It corresponds to the ooid grainstone facies with herring bones under tidal influence (facies F2b, Fig. 7B,E). This regressive systems tract is capped by the Bt3 maximum regressive surface (MRS), which is well expressed by a hardground (Fig. 5, Fig. 7G).

Sequence I is sub-horizontal and easy to follow on the outcrop. The transgressive systems tract of Sequence II is characterized by a 2 m-thick ooid grainstone deposit with cross-stratifications of several tens of centimeters in small to medium dunes. The top 50 cm is bioturbated.

The maximum flooding surface of Sequence II is located along a very well expressed and developed bioturbated surface (MFS Concinna, Fig. 5, Fig. 6E). The surface forms a very good marker at the quarry scale. The regressive systems tract of this Sequence II marks a shift from upper offshore

Facies modeling

To model the 3D facies heterogeneity of the carbonate deposits as accurately as possible, we adopted an iterative approach using stochastic methods such as the Petrel® "Truncated Gaussian with Trends", "Object Modeling", and "Assign Value" algorithms. The 3D geomodel encompasses the three sequences and four surfaces identified by sequence stratigraphy. The two maximum regressive surfaces, Bt3 and Bt4, and the two maximum flooding surfaces, MFS Concinna and MFS Digona, were digitized on VRGS®. Sixteen clinoform surfaces were also digitized on VRGS® as polylines. As some clinoforms could not be observed in some parts of the quarry, initial clinoforms were duplicated (29 times) with the same tangential dip (20°) and azimuthal (N70°, Fig. 3, Table S1) inclination to create a complete and continuous layout of the clinoforms from south to north, bringing the number of clinoform surfaces to 45. A total of 52 surfaces (45 clinoform surfaces and 7 structural modeling surfaces) were thus integrated into the model. Note that there is no difference between clinoform surfaces and structural modelling surfaces from the software point of view. We used these surfaces to distinguish three zones (Petrel zones) in the geomodel, representing the sequences, to populate the facies properties. The first zone (zone 1) groups Sequence I and the transgressive systems tract of Sequence II, both displaying a horizontal sedimentary pattern. The top of zone 1 corresponds to the maximum flooding surface MFS Concinna (Fig. 5). The second zone (zone 2) corresponds to the clinoforms, which belong to the regressive systems tract of Sequence II and are bounded by the MFS Concinna at the base and the maximum regressive surface Bt4 at the top (Fig. 5). Because of the particular geometry of the clinoforms, we modeled this unit separately to make the geometry of the cells fit the clinoform geobodies. All the clinoform surfaces (set as "conformable"), the Bt4 surface on top (set as "erosive") and the MFS Concinna (set as "base") were used as inputs in the "horizons" 29 window of the "structural modelling" panel to limit the toes of clinoforms in Petrel®. The third zone (zone 3) is composed of Sequence III, from the Bt4 maximum regressive surface to the top of the quarry. With the 52 surfaces created in the model, the three zones were divided into 51 sub-zones.

For Zones 1 to 3, the "simple grid" algorithm was used. Based on fieldwork and the 3D photogrammetric model on Pix4D, 30 measured (7) and synthetic ( 23) well logs were created to populate the model, providing matching points for the facies modelling (Fig. 9). The 23 synthetic logs were created according to the facies observations made in the photogrammetric model on VRGS® to map the vertical facies heterogeneity observed in the photogrammetric model. Once the facies were upscaled, we used the "Truncated Gaussian with Trends" algorithm implemented in Petrel® to closely control the distribution of the facies in the prograding trend in zone 2. It was also partially used in zone 1 and zone 3. Using this algorithm involves following a specific workflow for each zone: (1) choose the facies and check the facies order, (2) edit the transition of facies within the 51 sub-zones and (3) decide on the spatial extension of facies belts. For zone 1 (tidal to subtidal facies) and the lagoonal zone 3, other algorithms were marginally used coupled with statistical data analysis to obtain a satisfactory result for the facies distribution when compared to the photogrammetric model. The "Object Modeling (stochastic)" function was used to create geobodies (e.g. lenticular objects of tidally influenced facies F2c observed on top of the clinoforms in zone 2, Fig. 9). These lenticular geobodies cap the clinoform deposits and are rarely preserved but they have been seen in a few spots in the quarry (Fig. 9). The dimensions of the lenticular geobodies have been conservatively determined by measuring them on the photogrammetric model to reduce possible exaggerations on the semi-variogram range. They are created by wave influence at the top of the clinoforms due to a lower sea level. In addition, the "Assign Value" algorithm was used, where only one facies could be seen in a few sub-zones (F2b, F3c, Fig. 9).

Petrophysical characteristics and modeling

Four pore types, illustrated in Fig. 10, are observed in thin sections: (1) intergranular pores (Fig. 10D)

(2) vuggy, cavernous network linked to the dedolomitization of burrows within bioturbation surfaces or within perforations on hardgrounds (Fig. 10C), (3) moldic pores linked to the dissolution of dolomite rhombs (Fig. 10A), (4) and microporosity (Fig. 10C). We observe that bioclasts went through a micritization phase (ooids, peloids, micritized bioclasts, oncoids, foraminifers). Micritized grains are generally microporous throughout the log section. Micropores measure less than 10 µm and are well expressed by the greenish/bluish colour of resin in micritic grains (after the definition of [START_REF] Deville | Development of Water-Based Drilling Fluids Customized for Shale Reservoirs[END_REF]. Microporosity is homogeneously present within the model. Some microporosity is present even within the facies with the highest macroporosity. Four petrophysical properties are thus included in the model: (1) total porosity, (2) macroporosity,

(3) microporosity, and (4) permeability. They populate the grid according to their statistical distribution in each facies (Fig. 11).

In the model and accordingly to the petrophysical data, total porosity ranges from 11.3 to 20.6%, macroporosity from 0 to 15%, microporosity from 4.8 to 20.6%, and permeability from 0.1 to 5 mD (Table 2 and Supplementary Table S2). Macroporosity is variable in facies association FA1. No significant macroporosity is attributed to marls (facies F1a). The permeability data were not previously measured and a low range of total porosity (0.1 to 5%) and permeability (0.1 to 0.2 mD) is assigned from literature references (wackestone and mudstone zone of [START_REF] Nurmi | The land of giants[END_REF].

Measurements made on equivalent marl deposits in the Middle Jurassic of the Paris Basin [START_REF] Brigaud | Acoustic Properties of Ancient Shallow-Marine Carbonates: Effects of Depositional Environments and Diagenetic Processes (Middle Jurassic, Paris Basin, France)[END_REF]) support these values. In facies association FA1, the macroporosity ranges from 0.2 to 9% in F1b, from 8.5% to 15% in F1c, and from 0 to 3.3% in F1d. In this facies association, macropores are dominated by moldic pores forming euhedral rhombs or approximate rhombs. They are clearly related to the dissolution of dolomite crystals. It appears in thin-sections that facies F1c is the most porous and has been modeled this way in the static model. Facies F1b, which has also some of the highest macroporosity values of the dataset, is present at the base of the clinoforms. It thus appears that facies F1b and F1c are connected and both present good macroporosity, enlarging the area that could be qualified as fair reservoir, if we consider high permeability values. Total porosity in facies F1c ranges from 13.3 to 20.3%, demonstrating that microporosity contributes to the pore network. In facies association FA2, macropores range from 0.2 to 9.7% in F2a from 0 to 8.1% in F2b, and from 3.6% to 4.3% in F2c. Facies F2a displays high porosity in thin layers of a few centimeters on the hardgrounds, corresponding to the dedolomitization in perforations (mean of 2.6%). The dedolomitization of this facies is associated with bioturbation [START_REF] Purser | Dedolomite Porosity and Reservoir Properties of Middle Jurassic Carbonates in the Paris Basin, France[END_REF][START_REF] Loreau | Géométrie et couleur des couches dans l'oolithe blanche de Massangis (Yonne)[END_REF].

Dedolomitization of burrows generates vuggy pores along hardgrounds or bioturbated surfaces [START_REF] Purser | Dedolomite Porosity and Reservoir Properties of Middle Jurassic Carbonates in the Paris Basin, France[END_REF]. In lagoonal facies association FA3, macroporosity ranges from 0 to 2% in F3a, from 0.5% to 1.4% in F3b, and from 0% to 8% in F3c. Even if the different types of porosities vary across the geomodel, the permeability is quite low in general, ranging from 0.1 to 5 mD in most facies.

The petrophysical distribution "Gaussian Random Function Simulation" (GRFS) is used to populate the entire grid with properties, considering each facies independently. To populate all cells, even where no direct measurements exist, statistical values extracted from the data distribution (range, mean, and standard deviation) are defined and used to constrain the algorithm (Table 2). To add a degree of heterogeneity in the final petrophysical model it is possible to input anisotropy in the property body to create variations of said property in a constant facies when the GRFS alogorithm is applied. To illustrate the natural relationship between the facies and the petrophysical properties, all properties are upscaled to the model scale using the "arithmetic averaging" and "facies bias" options. The "total porosity" distribution is generated for each facies in all zones. The same method as for total porosity is used to upscale macroporosity and microporosity. However, to improve control of the property distributions, macroporosity and microporosity are co-kriged with the volume of the total porosity using the local varying mean correlation, which treats secondary data as a location-dependent mean. Finally, as with porosity properties, "permeability" is upscaled and populated in the grid using the same methodology. Permeability is co-kriged against total porosity as a trend, and macroporosity as a collocated co-kriging property which uses a correlation coefficient of the previously modeled secondary variable estimated internally in Petrel®. [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]from Nurmi, 1986).

**= ranges for F1c -peloidal packstone with crinoids, ranges were determined by projecting the macroporosity values measured on the samples of facies F1c onto the existing literature datasets [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]from Nurmi, 1986;[START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF].

Discussion

Depositional model

Three sedimentary systems developed successively during the Middle to Late Bathonian in Massangis quarry (Fig. 5 and Fig. 13). The first sedimentary system, developing during the deposition of the regressive of Sequence I and the transgressive systems tract of Sequence II, is dominated by intertidal to shoreface ooid environments, characterized by less than 10-15 m water depth considering classic depths for wave zonation [START_REF] Sahagian | Eustatic curve for the Middle Jurassic-Cretaceous based on Russian Platform and Siberian stratigraphy: Zonal resolution[END_REF][START_REF] Plint | Wave-and storm-dominated shoreline and shallow-marine systems[END_REF][START_REF] Andrieu | Linking early diagenesis and sedimentary facies to sequence stratigraphy on a prograding oolitic wedge: The Bathonian of western France (Aquitaine Basin)[END_REF]; Fig. 13A). The facies of this first sedimentary system are topped by the highly bioturbated Concinna maximum flooding surface, which is visible over the entire Burgundy platform [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] and as far as the western Paris Basin [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]. Clinoforms formed during the regressive systems tract of Sequence II (Late Bathonian) and correspond to sedimentary system 2 (Fig. 5 and Considering that the Bathonian formations had an overburden of 1000 to 1500 m [START_REF] Barbarand | Lower Paleogene denudation of Upper Cretaceous cover of the Morvan Massif and southeastern Paris Basin (France) revealed by AFT thermochronology and constrained by stratigraphy and paleosurfaces[END_REF][START_REF] Mangenot | An emerging thermochronometer for carbonate-bearing rocks: ∆47 /(U-Pb)[END_REF], the decompaction factor for grainstones is about 1.2 [START_REF] Goldhammer | Compaction and Decompaction Algorithms for Sedimentary Carbonates[END_REF]. The initial thickness of the clinoforms thus ranged between 18 and 20 m. Previous interpretations suggested that clinoforms correspond to a shoal complex forming the transition between a protected lagoon and a more open marine domain [START_REF] Casteleyn | Interrelations of the petrophysical, sedimentological and microstructural properties of the Oolithe Blanche Formation (Bathonian, saline aquifer of the Paris Basin)[END_REF][START_REF] Casteleyn | An integrated study of the petrophysical properties of carbonate rocks from the "Oolithe Blanche" formation in the Paris Basin[END_REF][START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF]. From this oolitic shoal, oolitic lobes prograded to a shallow tidally dominated lagoon. However, lagoon deposits and prograding oolitic lobes have never been observed or characterized in local studies simultaneously [START_REF] Casteleyn | Interrelations of the petrophysical, sedimentological and microstructural properties of the Oolithe Blanche Formation (Bathonian, saline aquifer of the Paris Basin)[END_REF][START_REF] Casteleyn | An integrated study of the petrophysical properties of carbonate rocks from the "Oolithe Blanche" formation in the Paris Basin[END_REF][START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF].

The thickness of the clinoforms in Massangis quarry does not match the expected thickness from prograding platform wedges (e.g. [START_REF] Pomar | Depositional model for a prograding oolitic wedge, Upper Jurassic, Iberian basin[END_REF]. The height of the clinoforms is consistently lower than for the Bathonian oolitic prograding platform wedge of Combe Brune quarry (Aquitaine basin, [START_REF] Andrieu | Linking early diagenesis and sedimentary facies to sequence stratigraphy on a prograding oolitic wedge: The Bathonian of western France (Aquitaine Basin)[END_REF]. Indeed, Combe Brune quarry presented clinoforms with an initial height of 35 to 60 m [START_REF] Andrieu | Linking early diagenesis and sedimentary facies to sequence stratigraphy on a prograding oolitic wedge: The Bathonian of western France (Aquitaine Basin)[END_REF]. A Middle Jurassic oolitic shoal system with a lagoon has been studied by [START_REF] Pierre | Spatial and temporal distribution of ooids along a Jurassic carbonate ramp: Amellago outcrop transect, High-Atlas, Morocco[END_REF] along the Amellago cliffs (High Atlas, Morocco). Clinoforms observed there prograde to the basinal domain, and the progradation of oolitic lobes toward the lagoon is very flat [START_REF] Pierre | Spatial and temporal distribution of ooids along a Jurassic carbonate ramp: Amellago outcrop transect, High-Atlas, Morocco[END_REF]. We propose here, following [START_REF] Vincent | Large-scale prograding oolitic geometries in the Middle Jurassic "Oolithe Blanche" Formation[END_REF], that the clinoforms correspond to sets of very large marine dunes of about 15-20 m height prograding N70° across the platform during a regressive systems tract (Fig. 13B). The paleobathymetry at the base of these very large dunes at the beginning of the regressive system tract could be 4 times the height of the very large dunes, i.e. 60-80 m [START_REF] Yalin | Mechanics of Sediment Transport[END_REF][START_REF] Dalrymple | Chapter 13 Estuarine Dunes and Bars[END_REF]. Such sedimentary structures are currently observed in the English Channel or on the Georges Bank on the USA/Canadian Atlantic Continental Shelf, for example [START_REF] Todd | International Year of Planet Earth 5. Applications of Seafloor Mapping on the Canadian Atlantic Continental Shelf[END_REF][START_REF] Ferret | Migration and internal architecture of marine dunes in the eastern English Channel over 14 and 56 year intervals: the influence of tides and decennial storms[END_REF]. Very large dunes on Georges Bank are 12 m heigh and they develop at a water depth of 50 m (in the troughs between very large dunes, [START_REF] Todd | International Year of Planet Earth 5. Applications of Seafloor Mapping on the Canadian Atlantic Continental Shelf[END_REF]. The third sedimentary system corresponds to the lagoonal environment of Sequence III (Fig. 13C). [START_REF] Enge | From outcrop to reservoir simulation model: Workflow and procedures[END_REF][START_REF] Buckley | Terrestrial laser scanning in geology: data acquisition, processing and accuracy considerations[END_REF][START_REF] Rotevatn | Overlapping faults and their effect on fluid flow in different reservoir types: A LIDAR-based outcrop modeling and flow simulation study[END_REF][START_REF] Marques | Virtual and digital outcrops in the petroleum industry: A systematic review[END_REF], drone photogrammetry allows to take pictures from unreachable places. Drones are easier to use and move around in remote places than the other technologies actually used in the field (e. g.

LiDAR). Quarries producing building stone and having vertical fronts cut at right angles are particularly suitable for 3D modeling. Kilometer scale drone acquisition is easy, relatively rapid (2 days of acquisition in this quarry) and the photogrammetric models are rather precise (mean ground sample of the photography dataset = 1.72 cm; note that only the most accurate pictures were used to get the value of the final model down to 0.94 cm). The complex geometries of clinoforms from hundred meters to kilometer scale within a carbonate platform can be easily observed, analysed, and interpreted after fieldwork with the photogrammetric model. The combination of direct outcrop observations and numerical models allows us to observe the complex geometries of the progradation within the stratigraphic framework in real 3D (Fig. 13). Coupled with a detailed and complete facies description, the photogrammetric model provides a continuous representation of architectures along the cliffs. Using VRGS®, it is then possible to manually draw surfaces on outcrops. All the datasets are easily imported into a geomodeler software package such as Petrel® (e.g., .ply, .obj., .las format files).

Continuous spot observations of facies in the photogrammetric model are used to create a reliable model using the "Truncated Gaussian With Trends" algorithm. Recent investigations have shown the benefits of the "Truncated Gaussian With Trends", "Truncated Gaussian Simulation", or "Sequential Indicator Simulation" algorithms when working with shallow carbonate ramp systems [START_REF] Brigaud | Characterization and origin of permeability-porosity heterogeneity in shallow-marine carbonates: From core scale to 3D reservoir dimension (Middle Jurassic, Paris Basin, France)[END_REF][START_REF] Tomassetti | Modeling lateral facies heterogeneity of an upper Oligocene carbonate ramp (Salento, southern Italy)[END_REF]. When continuous facies architecture observation and interpretation along a cliff is possible, "Truncated Gaussian with Trends" is suitable for creating a close-to-reality model because it is possible to control the geometry of facies transitions through the use of facies belts. This algorithm allows us to model reliable facies transitions and to easily reproduce heterogeneity when continuous images are available on the photogrammetric model. It helps to create a facies order with a configured facies geometry (facies belts) for each of the 51 zones of the model. This tool can be used to produce a highly deterministic facies model along the cliffs of the quarry where observations from the photogrammetric model are made to control the transitions, and the main directions of the clinoform geometry are used to extend the facies. Facies volumes can therefore be considered to be reliable at the quarry resolution XYZ 1 km x 0.4 km x 0.05 km, which is a valuable approach for reservoir exploration and exploitation (Fig. 13). Observations made on wellconstrained clinoforms from the photogrammetric model transcribe the geobodies to a higher level in the grid when compared to a simple "layer-cake" grid model without any implemented architecture, creating a built-in geobody geometry induced by the clinoform surfaces [START_REF] Cabello | From conventional outcrop datasets and digital outcrop models to flow simulation in the Pont de Montanyana point-bar deposits (Ypresian, Southern Pyrenees)[END_REF]. Coupling drone acquisition to an outcrop study permits facies observation on larger scales with a high degree of fidelity, and allows to obtain better visuals (orthophotographs), measure beddings (dip and azimuth) and use polylines to export isolated geobodies [START_REF] Hodgetts | Collection, processing, interpretation and modelling of digital outcrop data using VRGS: An integrated approach to outcrop modelling[END_REF][START_REF] Tomassetti | Modeling lateral facies heterogeneity of an upper Oligocene carbonate ramp (Salento, southern Italy)[END_REF][START_REF] Brandano | Facies heterogeneities and 3D porosity modelling in an Oligocène (Upper Chattian) carbonate ramp, Salento peninsula, southern Italy[END_REF][START_REF] Marques | Virtual and digital outcrops in the petroleum industry: A systematic review[END_REF].

Origin of the reservoir heterogeneities and possible reservoir analogues

The reservoir qualities of the studied carbonates are mainly controlled by two pore types: (1) microporosity and (2) secondary rhombic moulds. Blocky calcite cements and syntaxial overgrowths from crinoids reduced the intergranular porosity, which is rarely and poorly preserved in all the samples, even in grainstone textures. However, microporosity is developed in the micrite inside of grains (or due to micritization) and more rarely in the matrix, giving pore volumes of about 10 to 20%. In facies association FA2 and FA3, characterized by high microporosity, the porositypermeability relationship fits well with the relationship in microporous reservoirs, with permeability less than 5 mD [START_REF] Nurmi | The land of giants[END_REF][START_REF] De Periere | Morphometry of micrite particles in cretaceous microporous limestones of the Middle East: Influence on reservoir properties[END_REF] (Fig. 15). The rhombic pores are due to a dedolomitization process implying dissolution of a dolomite cement, precipitated as non-interlocking euhedral to subhedral rhombs between grains, and creating secondary pore space [START_REF] Makhloufi | Dolomitization of the Upper Jurassic carbonate rocks in the Geneva Basin, Switzerland and France[END_REF]. [START_REF] Mangenot | An emerging thermochronometer for carbonate-bearing rocks: ∆47 /(U-Pb)[END_REF] dated dolomite crystals using U/Pb geochronology in cores from the Middle Jurassic oolitic reservoir located 150 km north-west of Massangis and concluded that the dolomite probably precipitated during the Early Cretaceous. In the eastern Paris Basin at 150 km north east of Massangis, oolithic facies with dolomite associated to syn-burial precipitation were estimated as Early Cretaceous too [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]. Dolomite dissolution sometimes tends to mimic intergranular porosity. However, irregular grains with flat surfaces and small depressions indicate that dolomite crystals were dissolved (Fig. 10D). The existence of reservoirs principally characterized by moldic dedolomitized porosity implies a partial dolomitization of intergranular space within oolitic or peloidal deposits as a cement and the preferential dissolution of dolomite, preserving calcite grains [START_REF] Purser | Dedolomite Porosity and Reservoir Properties of Middle Jurassic Carbonates in the Paris Basin, France[END_REF]. Precipitation of dolomite cement and later dedolomitization are identified mainly in FA1, especially in facies F1b and F1c, and are also well-expressed at a metric and Oolite Blanche reservoir around Saint-Martin-de-Bossenay (PICOREF Sector, [START_REF] Delmas | Petrophysical Properties of the Middle Jurassic Carbonates in the PICOREF Sector (South Champagne, Paris Basin, France)[END_REF].

Conclusion

Massangis quarry allows us to observe a cross-sectional analogue of a prograding oolitic clinoform depositional system. Our work leads us to the following conclusions: (2) The use of drone photogrammetry makes it possible to constrain the facies architecture and heterogeneities even along inaccessible vertical faces of the quarry. Integration with field observations can greatly reduce the uncertainties of a model that arise from the structural framework or crude property modeling. The reduced uncertainty of the facies model reduces the uncertainty of the porosity and permeability properties in the final model. Another interesting aspect of the use of drones is the possibility of precisely studying outcrops from different perspectives after field work without the need of carrying heavy equipment. More outcrops will be observable with ease and at low cost.

(3) Using a combination of Pix4D®, VRGS®, and Petrel® software packages, it is possible to model complex geobodies in reservoirs with a highly detailed architecture serving as a basis for improved reservoir simulations. This enables us to better understand and model the complexity of a clinoformtype depositional environment. The use of surfaces generated in the grid helps in creating a realistic architecture model at cell level. This makes it possible to create a petrophysical model of clinoform geobodies in carbonates with a resolution never obtained before. By following the workflow exposed in this paper, it is possible to model rapidly any kind of complex-shaped geobody in a 3D space. The association of a fast acquisition method coupled with photogrammetric modeling allows extracting additional quantitative information at any time in the laboratory within a virtual outcrop, such as logs, facies observations and mapping, dips and azimuth measurements of sedimentary structures.

These measurements and quantitative information are some keypoints of the benefit of drone photogrammetry in this modeling workflow. This approach enables in a further step accurate hydrodynamic simulation in complex and heterogenous reservoirs.

(4) The relationship between the different porosities and permeabilities indicates that the Massangis quarry model is a good analogue for reservoirs exhibiting both microporosity and secondary porosity associated with dedolomitization. While Massangis quarry has often been studied as an analogue of the Oolite Blanche geothermal reservoir of the center of the Paris Basin, the porosity types and resulting permeabilities are different. However, some oil-reservoir units of the Middle Jurassic limestones of the Paris Basin exhibit similar moldic rhombohedral pore spaces. This model could thus serve as an analogue for a carbonate reservoir dominated by rhomb-molds.

  presented a model of facies heterogeneity combining digitalization by LiDAR technologies and geostatistical algorithms in Petrel® for a carbonate ramp analogue in Morocco. Other examples include the outcropping Miocene and Kimmeridgian carbonate ramps in southern Italy and Spain which have been modeled

Figure 1 .

 1 Figure 1. A-Geological map of the Paris Basin with the quarry location in the south-eastern part of the Basin (Burgundy, France). B-Detailed map of the quarry and its surroundings. C-Aerial view of Massangis quarry with its two sites, Vaurion on the left and Malpierre on the right, with the location of sampling sites and the number of samples taken from each site.

  interpreted the clinoform as a shoal system separating a protected lagoon from a more open marine environment. However, marls typical for a laterally protected lagoon and open marine facies are absent in local studies. More recently, Vincent et al. (2017) re-interpreted these clinoforms as giant (very large) migrating dunes on a large gently sloping bank flanking the Burgundy platform. Stratigraphic sequences and sedimentary geometries of the Bathonian deposits in the Paris Basin have been recognized since the pioneering work of Vail et al. (1987), and have been refined by various studies

Figure 2 .

 2 Figure 2. A-Chrono-stratigraphic framework of a west-east transect across the Paris Basin and location of Massangis quarry (profile A-B on the map below). B-Paleogeographic map of France during the Late Bathonian with the location of the quarry.

Figure 3 .

 3 Figure 3. Overview of the workflow for generating the 3D facies, porosity, and permeability model. First the pictures were captured by a drone equipped with 20 Mega Pixel camera and GPS during two days of field work at Massangis quarry.

Figure 4 .

 4 Figure 4. Distribution of the Bathonian facies on a theoretical downdip profile from the bottom of the clinoforms to the lagoonal environment. This synthetic profile shows the succession of the facies on the platform, depending on bathymetry.

Figure 5a .

 5a Figure 5a. Synthetic sedimentological log of Massangis quarry illustrating facies, paleoenvironment interpretation, and sequence stratigraphy.

Figure 6 .

 6 Figure 6. A-Photograph of a vertical face of excavation of the quarry showing facies F1a marls, F1b bioturbated peloidal crinoidal packstone, F1d ooid and peloidal grainstone with crinoids, and F2a ooid grainstone with bioturbations and hardgrounds. The bioturbated surface is arrowed (Thalassinoides). B-Detail showing facies F1a marls with some crinoids. C-Close-up of facies F1b bioturbated peloidal crinoidal packstone with a redox front. D-Facies F1b bioturbated peloidal (Pel.) crinoidal packstone with echinoderms (Ech.), bivalves (Bi.), and ooids under microscope (sample MAS2.01). Note the secondary macroporosity (moldic porosity) due to a dedolomitization phase. E-Blow-up of the bioturbated surface and lithoclasts between F1b and F2a made by Thalassinoides. F-Facies F1c peloidal (Pel.) packstone with crinoids, bivalves (Bi.), and echinoderms (Ech.) (sample MAS-10). Note the porosity of this facies is due to a more intense phase of dedolomitization (Dedol.). G-Example of facies F1d from the photogrammetric model of the quarry showing clinoform surfaces. H-Thin section example of facies F1d, peloidal (Pel.) and ooid grainstone with crinoids (sample MAS-15). Note the blocky calcite between the grains. Refer to figure 4 for colour code of facies types.

  Fig 5, Fig. 7G, 7H). Blocky calcite cements are well-developed in this facies and the main pore type is microporosity. 20 Interpretation -The widespread perforations, grain-supported textures, and bioclast fragments indicate high-energy environments. In facies F2a, sigmoid and planar cross-bedding in small to medium dunes indicate shallow high-energy wave-dominated conditions in a shoreface environment. Besides, herring-bone stratification in facies F2b indicates tidal influence in high-energy shoreface to foreshore environments.

Figure 7

 7 Figure 7. A-Photograph of the tidal facies at the base of the quarry in the Oolite Blanche Formation. B-Sketch from photo A showing sedimentological features and facies. C-Detailed photograph of facies F2c with bioclasts. D-Facies F2a ooid

Figure 8 .

 8 Figure 8. A-Facies F3a peloidal (Pel.) grainstone with echinoderms (Ech.) and isopachous cements (sample MAS-23). B-Facies F3a peloidal (Pel.) grainstone with foraminifers (Foram.), bivalves (Bi.) and isopachous cements (sample MAS-22). C-Photograph of the outcrop with facies F3b and F3c. D-Facies F3b ooid and oncoidal grainstone with Oncoids (Onco.), ooids, peloids (Pel.), bivalves (Bi.) and dedolomitization (Dedol.) (sample MAS-21). E-Facies F3c oncoidal (Onco.) rudstone with gastropod (Gast.), bivalves (Biv.), and echinoderm spines (Ech. Spine) presenting interparticle pores (sample MAS-20). Refer to figure 4 for colour code of facies types.

(

  facies F1a) to shoreface facies (facies F1d). It corresponds to a 15-20 m-thick succession of clinoforms composed of peloidal and ooid or bioclastic grainstone (facies association FA1, Fig.5). The main direction of progradation of the clinoforms is N70°. The bottom of the clinoforms is composed mainly of a highly bioturbated facies F2b, and locally of marls (facies F1a). Near the toe of the clinoforms, deposits are tangential to the base, and the dip of the sets is 7°. The dip angle steepens along the clinoform sets reaching a maximum of 25°. The top of the clinoforms exhibits a hardground corresponding to the maximum regressive surface Bt4.The transgressive systems tract of Sequence III is composed of 7 m-thick lagoonal facies (Facies F3b) topped by a bioturbated surface corresponding to the maximum flooding surface (MFS Digona, Fig.5, Fig.8C). The regressive systems tract of Sequence III is characterized by the development of lagoonal packstone to floatstone (facies F3c, Fig.10B, D, and Fig.11A).

Figure 9 .

 9 Figure 9. A-Panoramic view of the 3D model in VRGS®. B-Panorama of the quarry painted, showing the facies heterogeneity observed in the quarry. C-Relief-extruded part of the same part of the quarry showing the correlation between the field and photographic observations and the final facies model when controlled by the logs using synthetic logs. D-Close-up snapshot of the cells of the model with their dimensions. Note the clinoform polylines appearing out of the cliffs of the modeled quarry.

Figure 10

 10 Figure 10. A-Thin section of sample MAS-19 Illustrating some moldic macroporosity (dedolomitization-induced) and microporosity. B-Thin section showing porosity totally blocked by cements (MAS-10). C-Thin section showing truncated

Figure 11

 11 Figure 11 Total porosity and permeability distribution of the different facies in the quarry along the synthetic depositional profile. In the upper part of the figure, 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.

Figure 12 .

 12 Figure 12. Relief-extruded views of the different petrophysical properties in the geomodel (note the logs representing the facies). A -Total porosity, B -Macroporosity, C-Microporosity, D -Permeability. Refer to figure 4 for color code of facies types.

Fig. 13 )

 13 Fig.13). The bottomset of the clinoforms is marked by marls (facies F1a) or intensively bioturbated

Figure 13 .

 13 Figure 13. Depositional model detailed in three steps: (A) Sedimentary system 1 corresponding to Sequence I and the of Sequence II; (B) End of sedimentary system 2 corresponding to the Regressive Systems Tract of Sequence II; and (C) Sedimentary system 3 corresponding to Sequence III. Please refer to Fig. 5b for the legend of organic/inorganic elements.

Figure 14 .

 14 Figure 14. Complete model of the Massangis quarry with continuous properties.

  scale in facies F3b, in relation to the bioturbated surface below MFS Digona. The moldic porosity due to dedolomitization is also observed in much thinner layers in facies F2a due to dissolution in the perforations of the hardground characterizing Bt3. These diagenetic phases closely control the macroporosity in our model and produce a major control over reservoir heterogeneity in the Massangis (Fig.15;[START_REF] Delmas | Petrophysical Properties of the Middle Jurassic Carbonates in the PICOREF Sector (South Champagne, Paris Basin, France)[END_REF]. These variations in properties between the different stratigraphic analogues found in the center (Vert-Saint-Denis 1 well data) and on the edge of the Paris Basin (Massangis quarry, Saint-Martin-de-Bossenay by[START_REF] Delmas | Petrophysical Properties of the Middle Jurassic Carbonates in the PICOREF Sector (South Champagne, Paris Basin, France)[END_REF]) could be explained by the variation in the depositional environment, the diagenetic history of these zones and by different burial conditions.

Figure 15 -

 15 Figure 15 -A -Nurmi's abacus defining the relations between porosity and permeability in different carbonate textures (Middle East). B -Graphical representation of the samples of the quarry and data from the Vert-Saint-Denis VSD-1 borehole

( 1 )

 1 A new depositional model for the Late Bathonian clinoforms in Massangis quarry is proposed: the clinoforms correspond to sets of very large marine dunes of about 15-20 m in height prograding N70° across the platform. Examples of such sedimentary structures are currently observed in the English Channel or on the Georges Bank in Canada. This outcrop analogue for a world class carbonate geothermal reservoir in the Paris Basin could be used as a proxy for clinoform-based geobodies around the world.

Table 1 -

 1 Facies table of the samples from the quarry. R = Rare < 10%; C = Common 10 -20%; F = Frequent 20 -50%; A = Abundant > 50%

Table 2 -

 2 Data ranges used for the population of petrophysical properties * = ranges for F1a -marls, theoretically determined by the texture and the property behaviour of this type of rocks

	Properties

Figure 5b. Legend of the log of the quarry (Fig 5a)

Max: 2 430 µm

Acknowledgements

We thank Polycor for allowing us access to the quarry and for authorizing drone photogrammetry. This work is supported in part by the "Geothermal water circulation in heterogeneous aquifers" project funded by the "Paris-Saclay Initiative d'Excellence", a state aid program managed by the Agence Nationale de la Recherche (ANR-11-IDEX-0003). This work is also supported by the program UPGEO "UPscaling and heat simulations for improving the efficiency of deep GEOthermal energy" funded by the Agence Nationale de la Recherche (ANR-19-CE05-0032-01). 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 throughout its utilization. We also thank Virtual Reality Geological Studio (special thanks to David Hodgetts) and Pix4D for providing annual academic licenses for VRGS® and Pix4D®. We thank Sketchfab® for providing a free online repository. We sincerely thank the two anonymous reviewers who helped to greatly improve the manuscript and Laura Tomassetti as the guest editor of this special issue.

Data Availability

Datasets of the 3D model of Massangis quarry containing the Digital Outcrop Model (DOM), construction surfaces, facies and petrophysical data are included in the rar file, which can be found at https://data.mendeley.com/datasets/fcdg3jvjdp/1, an open-source online data repository hosted at Mendeley Data. This file can be read only with Schlumberger software Petrel®, version 48 2019.2 and later. The 3D photogrammetric model of Massangis quarry can be found on Sketchfab at https://skfb.ly/6RYGF with highlights.

limestone of this study. Dedolomitization and development of microporosity in micrite are probably related to flow of freshwater during telogenesis consecutive to the uplift of the limestones during the Cenozoic [START_REF] Barbarand | Lower Paleogene denudation of Upper Cretaceous cover of the Morvan Massif and southeastern Paris Basin (France) revealed by AFT thermochronology and constrained by stratigraphy and paleosurfaces[END_REF].

Porosity -permeability (phi-k) data from the Vert-Saint-Denis VSD-1 borehole (C.E.P., 1964, Fig. 1A), which is drilled into the geothermal reservoir of the Mid to Late Bathonian Oolithe Blanche Formation near Paris, are presented in a phi-k cross-plot together with the phi-k of Massangis quarry (Fig. 15;[START_REF]Vert-Saint-Denis 1 EMSCO GB 250[END_REF]. At 1600 m depth, close to the center of the Paris Basin, the Oolite Blanche formation displays different properties than at Massangis quarry, with permeabilities exceeding 10 mD (C.E.P., 1964). Porosity is dominated by intergranular pores remaining partially uncemented during burial and allowing high permeabilities from 100 mD to 1-3 Darcy [START_REF] Robelin | Differenciation des ciments carbonates par cathodoluminescence; cas des reservoirs geothermiques du Dogger d'Aulnay-sous-Bois (Seine-Saint-Denis, France)[END_REF][START_REF] Rojas | Caractérisation et modélisation du réservoir géothermique du Dogger[END_REF]. In contrast, the Oolite Blanche reservoir at the subsurface 70 km north of Massangis (Fig. 1A properties are rather poor because of pore throats smaller than 1 µm [START_REF] Delmas | Petrophysical Properties of the Middle Jurassic Carbonates in the PICOREF Sector (South Champagne, Paris Basin, France)[END_REF]. Some highly permeable reservoirs of the Paris Basin subsurface are characterized by dedolomitization processes resulting in porosities of 10-20% and permeabilities of 10 mD to 1 D [START_REF] Purser | Dedolomite Porosity and Reservoir Properties of Middle Jurassic Carbonates in the Paris Basin, France[END_REF]. Small oil fields of the Paris Basin (Coulommes, Chailly-en-Bière, and Saint-Martin-de-Bossenay) are locally dominated by secondary moldic rhombohedral pore spaces in 2-5 m-thick ooid units [START_REF] Purser | Dedolomite Porosity and Reservoir Properties of Middle Jurassic Carbonates in the Paris Basin, France[END_REF].

Moreover, the Late Jurassic geothermal reservoirs in the Geneva Basin (Switzerland) are controlled by dedolomitization either through calcite replacement or dissolution in granular facies [START_REF] Makhloufi | Dolomitization of the Upper Jurassic carbonate rocks in the Geneva Basin, Switzerland and France[END_REF]. Dedolomitization is also well-expressed in the Late Jurassic of the French and Swiss Jura Mountains [START_REF] Rameil | Early diagenetic dolomitization and dedolomitization of Late Jurassic and earliest Cretaceous platform carbonates: A case study from the Jura Mountains (NW Switzerland, E France)[END_REF]. Hence our model appears to be an interesting analogue for subsurface ooid reservoirs in which microporosity and moldic porosity is due to dedolomitization. The reservoir of Saint-Martin-de-Bossenay is described as such, and a correlation can be found between the characteristics measured in the samples of the study and the samples analysed in the quarry of

Supplementary material

Table S1. Dataset of the beddings orientations extracted from the VRGS ® software and providing the data to reconstruct the stereonet in Fig. 3. Table S2. Dataset of facies and petrophysical properties (total porosity, macroporosity, microporosity, permeability) of samples collected from bibliography (Castelyn et al., 2010,2011[START_REF] Makhloufi | Impact of sedimentology and diagenesis on the petrophysical properties of a tight oolitic carbonate reservoir. The case of the Oolithe Blanche Formation (Bathonian, Paris Basin, France)[END_REF].