Clément Lacroix
email: clement.lacroix@promes.cnrs.fr

Mahmoud Hacheche

Maxime Perier-Muzet
email: maxime.perier-muzet@univ-perp.fr

Driss Stitou
email: driss.stitou@promes.cnrs.fr

Solar-driven thermo-hydraulic process for reverse osmosis desalination

Keywords: Solar desalination, thermal energy, reverse osmosis, thermo-hydraulic cycle, dynamic modeling

Existing distillation-based desalination processes are highly thermal energy consuming. Reverse osmosis (RO) technique is more efficient than thermal-based processes but it remains a solution that still induces high operating and maintenance costs. In this paper, an innovative thermally powered RO-based desalination process is presented. This new RO thermo-hydraulic process enables the pressurization of the salty water beyond its osmotic pressure to allow the permeation water through a semi-permeable membrane, thanks to a piston or an elastic bladder that is set in motion in a reservoir by a working fluid following a thermodynamic engine cycle similar to an Organic Rankine Cycle. The evaporator is heated by low grade heat (70 to 80°C) such the one delivered by plate solar collectors, while the condenser is cooled by the concentrated salty water. In order to enable a continuous drinkable water production, this process needs to implement two reservoirs, alternatively connected either to a high pressure evaporator or to a low pressure condenser. Such installation, designed here for brackish water desalination (5 g/liter), should enable an average daily production of 300 liters of drinkable water per m² of solar collectors with a production cost below 4€/m 3 . That technology seems to be relevant for small scale (5 to 10 m 3 /day) the daily water needs of people living in remote areas, in accordance to the location and the solar resource. A modeling of the whole process, considering a quasi-steady state approach has been developed in order to study its dynamic behavior, optimize its design and maximize its performances. This paper presents the preliminary results relative to the performance of such solar-driven desalination process.

Introduction

Growth of the demand in seawater and brackish water desalination over the past years is motivated by the increasing water stress in more and more countries. Even if the major part of desalination systems is powered by fossil energies, desalination systems driven by renewable energies are now well known and available on the market. Solar desalination is particularly interesting considering the most water lacking areas have an important sun resource. Solar driven desalination can be led by two different ways: distillation or membrane separation. In solar distillation systems, the water is evaporated by the solar heat. The most current distillation based technologies are the multi-effect distillation (MED) and the multi-stage flash distillation (MSF). Those processes have a specific energy consumption (SEC) commonly ranging between 30 and 80 kWh of thermal energy per m -3 of produced drinking water [START_REF] Ng | Solar-assisted adsorption cycle for the production of cooling effect and potable water[END_REF]. Membranes-based desalination consists in separating the salt from water through a semi-permeable membrane. The most implemented membrane separation processes are electro-dialysis and reverse osmosis. In electro-dialysis processes, the salt is drained through an ion selective membrane by the way of two electrodes submitted to electrical potential. This process needs electrical energy to operate and leads to an energy consumption lower than 2 kWh elec .m -3 . This process is more appropriate for brackish water desalination with a salinity below 2 g/l [START_REF] Malek | Decentralized desalination of brackish water using an electrodialysis system directly powered by wind energy[END_REF].

Reverse osmosis (RO) enables to desalinate water with higher salinity through a semi-permeable membrane, thanks to the pressurization of the salty water beyond its osmotic pressure, which is proportional to the water salinity and temperature. Such a process needs mechanical energy to run, which is commonly provided by a high pressure electrical pump. The RO technology leads to specific energy consumptions ranging from 2 to 10 kWh mec .m -3 , and depending on the water salinity and temperature. This low energy consumption is obtained by recovering the mechanical energy of the pressurized concentrated water at the membrane output. Distillation-based systems are easier to set up but their energy needs are higher than membranebased technologies. Membranes-based desalination enables to obtain low energy consumption in comparison to distillation processes but they have some major operating issues, such as membrane deterioration or fouling, that lead to higher operating costs. Nevertheless, reverse osmosis is an interesting technology because of its low energy needs and its potential to be implemented at wide scale. Exploiting solar energy via photovoltaic panels in order to produce electricity that is later converted into mechanical energy, implies the use of expensive mechanical energy conversion devices and electricity storage devices and leads to an important loss of efficiency. Solar heat-driven RO-based systems obtained by coupling a di-thermal power cycle with a membrane separation unit could be an interesting solution to avoid this energy conversion loss and obtain better installation efficiency. Manolakos et al. [START_REF] Manolakos | On site experimental evaluation of a low-temperature solar organic Rankine cycle system for RO desalination[END_REF] studied and experimented such kind of desalination process by implementing an Organic Rankine Cycle (ORC), and obtained relatively good performances. With the same idea, a new likely process is investigated here in this paper. The new thermo-hydraulic process presented enables to use directly the expansion work of a fluid that follows a thermodynamic engine cycle, in order to pressurize the saline water in the RO unit. This solution enables to reduce importantly the energy conversion losses occurring during the mechanical work conversion and improves the global process efficiency. The process enables also to implement simple hydraulic and thermal components in order to be cost effective, easy to maintain, and to operate autonomously in off-grid remote areas by using low temperature solar heat provided by common flat plate solar collectors. In this paper, such a new solar heat driven RO desalination process is described. A modelling of the whole process is carried out in order to evaluate its performance. The results are discussed in order to assess its relevancy.

System description

The thermo-hydraulic desalination process studied in this work and described Fig. 1 results from the direct coupling of a solar-powered Rankine-like engine cycle for the work generation with a reverse osmosis unit for the desalination. A set of two hydraulic cylinders are implemented for recovering the hydraulic energy of the brine flowing out the RO unit and making circulate the saline water in the process and the thermodynamic fluid of the engine cycle . The engine cycle is mainly composed of an evaporator and a condenser. The evaporator, which is heated by the solar heat, generates vapors of the working fluid at high pressure. This vapor flows into a first reservoir that is filled with saline water and pressurizes it above its osmotic pressure. The condenser, which is cooled by the saline water to be desalinated, enables the condensation of the working fluid previously accumulated in a second low pressure reservoir. The pressurization of the water in the tank is realized thanks to a free mobile piston set in motion by the vapor of the working fluid (Butane) that is produced at high pressure (approximately 15 bars) by evaporation in the evaporator of the thermodynamic cycle thanks to solar heat delivered by simple flat plate collector at about 70 to 90°C. The hydraulic cylinders are set in movement thanks to the pressurized concentrated water flowing out of the RO module, recovering thus its hydraulic energy. A first hydraulic cylinder is used to pump and fill with water the low pressure reservoir and a second hydraulic cylinder is used to transfer the working fluid from the low pressure condenser to the high pressure evaporator. Each cylinder is composed of three chambers in order to make it behave in a symmetric way in both directions of movement, that is to say in order to ensure that the pumped and the pushed back volumes of water inside the chambers are equals. The pressurization is done in the middle chamber that contains the driving piston, which transmits the mechanical energy and set in movement the piston in the two lateral chambers in order to realize the pumping in a first step and backflow in a second step of the water to be desalinated, or the working fluid into the evaporator. Thanks to the first hydraulic cylinder, the feed water to be desalinated is firstly pumped at ambient temperature at a pressure just below the atmospheric pressure, then go through the condenser, acting thus as a cold source for the thermodynamic cycle and is temporarily accumulated in one of the lateral chamber of the hydraulic cylinder. Simultaneously, as the piston moves, the water to be desalinated, which was previously accumulated in the second lateral chamber of the hydraulic cylinder, is pushed out in order to fill the low pressure water tank connected to the condenser. At the same time, the water to be desalinated that is pressurized in the other tank by the high pressure vapor of the engine cycle, flows through the RO membrane. Clean water (permeate) is produced and the more concentrated water, which is still at high pressure, is used to drive the above mentioned hydraulic cylinders. Several distributors and valves allow the role inversion of the reservoirs and the brine distribution in the process, as well as the control and the cycling of the different operating steps of the process. The retentate flow is controlled at the cylinders output by a constant flow rate valve, and regulating thus moving velocity of the cylinders pistons. All process elements requires to be salt-resistant.

Description of the operating cycle

The operating cycle of the process consists in several steps that are described in Fig. 1. An initialization phase enables to put the feed tank, which is already completely filled with feed water, at the high pressure P h by connecting the feed tank to evaporator. This high operating pressure level has to be preferably beyond two to three times the feed water osmotic pressure π to overcome transfer resistance throughout the RO membrane (see eq. ( 14)) The first step (Fig. 1a) consists in realizing the feed water desalination under ideally constant pressure. The high pressure feed water contained in the pressurized tank goes through the RO module and produces clean water. The still pressurized output brine is then used by the first hydraulic cylinder to pump the saline water and fill the second tank at the condenser pressure. The water filling of this tank discharges the vapors of the working fluid to the condenser where they condense at low pressure. The second step (Fig 1b) is also a desalination phase, but under a decreasing pressure. The high pressure tank is disconnected from the evaporator, stopping supplying the tank with high pressure vapor. Then the pressure tank decreases until the osmotic pressure, enabling a decreasing clean water production to carry on until this point. At this step end, the second tank is completely filled with feed water that will be later desalinated and all of the working fluid previously contained in it has been expelled into the condenser. As the pressure in the first tank becomes lower than the osmotic pressure (Fig. 1c), the remaining feed water cannot be desalted any longer. Nonetheless, this still pressurized water is used by a second hydraulic cylinder to transfer the liquid working fluid from the low pressure condenser to the high pressure evaporator. The vapor contained in this first tank depressurizes till the condenser pressure. At this point, a half-cycle has been achieved: the first tank is then fully filled with vapor and the second is full of feed water, both at low pressure, i.e. the working fluid saturation pressure Figures 1 represents the evolutions of the volumes of the vapor of the working fluid contained in the reservoirs V α and V β , V tank ,. They respectively correspond the feed tank volume at the end of phase 1 (at high pressure P h ), the feed tank volume at the end of phase 2 (at osmotic pressure ) and to the fully filled tank volume (at low pressure P l ). An isentropic expansion of the working fluid has been considered in order to calculate these volumes (1) :

      k l h V P V V P tan .   , (1)
The design of the hydraulic cylinders is currently led for optimizing the mechanical energy recovery. Their technology is under development in order to optimize the pumping depending on the energy recovered.

A first evaluation of the behaviour of this cyclic process is led here throughout a simple model. This modelling has been carried out by considering in first approximation a quasi-static operating mode, with the aim to first assess the process performances and its competitiveness by taking into account the variability of the solar resource and its cyclical operating.

Process modeling

The process behavior has been modeled by considering a nodal approach. All components of the system are assumed to be at a uniform temperature, pressure and concentration. The following mains assumptions have been considered:  Thermal equilibrium between the liquid and gaseous phases in the evaporator and the condenser  Frictions forces are negligible in the hydraulic cylinders and pipes  Condensation effects are negligible in the tanks and pipes

Solar collector coupled to the evaporator

The evaporator of the system is heated by a heat transfer fluid circulating in a solar collector with a given mass flow rate, corresponding approximately to the commonly used value of 0.0139 kg/s.m².

Fig. 2 : Solar panel and evaporator

The solar irradiation that is absorbed by the collector, of area A c , and transmitted to the evaporator by the heat transfer fluid can be written as follows considering the collector efficiency:

  in out p c c c T T C m IA       (2)
By taking into account an exchanger effectiveness value  of 0.7 for the evaporator (3), the temperature difference of the heat transfer fluid can then be expressed in the equation (4) as a function of the collector average temperature c T and the evaporator temperature T liq .

liq out c in c out c T T T T    _ _ _  (3)   liq c evap evap in c out c T T T T       2 2 _ _ (4)
The mean collector temperature c T is calculated iteratively by combining equations ( 2) and ( 4) and considering the commonly used expression for the collector efficiency c  , which is a function of the optical efficiency  0 , the mean fluid temperature and the solar irradiation,

                     I T T a I T T a IA T T C m amb t c amb t c c liq t c p 2 2 1 0 2 2     (5)
where a 1 and a 2 are respectively the linearly and quadratic thermal loss coefficients of the solar collector in W.m -2 .K -1 and W.m -2 K -2 .

T c_in T c_out I P ev T ev N liq P ev T ev N g c T g n  liq n  c 
Furthermore, the molar quantities of working fluid that are present in gaseous and liquid phase in the evaporator, can be obtained by considering energy and mass balances ( 6) and ( 7). An iterative procedure enables then to determine at each time step, the amount ev n  of liquid that is evaporated in the evaporator.

g g liq liq c ev h n h n dt dU       (6) Where       amb ev ev wall v ev ev g g ref ev pl l ev T T c m RT h N T T C N U       _ (7) ev liq liq n n dt dN     and g ev g n n dt dN     (8)-(9)
The molar flow rate g n  of the vapor that leaves the evaporator to the reservoir in order to pressurize the feed water is determined by the characteristics of the valve, i.e. its flow coefficient K v , and the pressure difference existing then between the evaporator and the reservoir

res ev v g P P K n   . 
A state equation for the working fluid is also considered in order to determine the pressure evolution in the evaporator and the pressurized reservoir, and then the pressure difference between these two components. The behavior of the working fluid vapor is in first approximation described by a virial equation of state (10):

g g g g V N T B RT N PV ) ( 1  (10)

Transfer Tanks

The volumes V f of saline water that is present in each transfer reservoirs are calculated thanks to a volumetric flow equation (11) knowing tank output mass rate of the feed water. The molar quantities of the gaseous working fluid in the reservoirs are also obtained throughout the molar balance (13). The polarity signs express the filling or emptying of reservoir with the saline water:

f s g res V V V   Re _ (11) f f f m dt dV     and g g n dt dN   
(12)-(13)

Reverse osmosis module

The permeate flowing through the osmotic membrane is determined by the mass transfer equation ( 14) established by Lonsdale et al. [START_REF] Lonsdale | Transport Properties of Cellulose Acetate Osmotic Membranes *[END_REF]:

  m m m w w p P A L m        (14)
where the trans-membrane pressure m P  is expressed as a function of the mass flow rate f m  , the pressure f P of the saline water, the pressure of the permeate that is supposed to be at atmospheric pressure, and the osmotic pressure  m that is function of the salt concentration and temperature of the feed water :

  f p f m RT C C    (15)

Condenser

The modelling of the condenser, which is cooled by the pumped feed water, is quite similar to the one developed for the evaporator. However, it is assumed that the liquid level in the condenser is kept at a minimal value, at about 1% of the total volume of the condenser, in order to fulfil the thermodynamic equilibrium condition (both states are required).

Fig. 4: Representation of the condenser

The excess liquid, which is produced by condensation of the vapour coming from the low pressure reservoirs, is transferred by gravity to a bottom receiver tank. The mass balance (15) is used to obtain the gas molar quantity present in the condenser.

liq g g n n dt dN     ( 15
)
An enthalpy balance ( 16) is also applied to the condenser to determine its condensing temperature T cond and pressure P cond .

cond liq liq g g h n h n dt dU       (16)
As the feed water flows through the condenser at a temperature T f close to the ambient one, in order to cool it, the condensation heat absorbed by the feed water can be expressed by considering the heat exchanger effectiveness cond  of the condenser :

    f cond cond f f in f out f f f cond T T cp m T T cp m       . . . . _ _   (17)
The flow rate of pumped the feed water is given by the displacement speed and the volume chamber of the pumping hydraulic cylinder.

Hydraulic cylinders

In order to keep an identical operating during the two half cycles, that is to say the same involved saline volumes, the hydraulic cylinders are designed with three chambers in which each of them moves a piston (Fig 4) The central motive piston is set in motion thanks to the pressurized concentrated saline water flowing out the RO module. The movement of the piston rod enables in one hand to pump the feed water and store it in a first chamber of the hydraulic cylinder. Simultaneously, it enables also to transfer to the second reservoir the feed water that was previously stored in the third chamber of the hydraulic cylinder and pressurize it at the condenser pressure. The displacement of the hydraulic cylinder is controlled by a control valve that keeps constant the mass flow rate of the exhaust water from the motive piston.

Results and discussions

This first modeling approach enables to simulate the process on a full day in order to analyze the dynamic behavior of this thermo-hydraulic solar driven desalination process. First simulations have been carried out with the physical parameters and simulation conditions summarized in table 1.

The thermodynamic properties of the working fluid (butane) and the saline water resulting from RefProp database (NIST database) have been considered [START_REF] Lemmon | NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties (REFPROP), Version 9.0[END_REF].

A preliminary results permit to obtain a first design of the process, i.e. solar collector area of 2 m² and a 112 m² RO membrane, corresponding the entry parameters of Table 1. 6 shows the pressure and volumes evolutions in the main components of the process operating under a solar irradiation of 800 W/m². These results show a highly dynamic system with half-cycles of about 125s, corresponding to a first step of 33s, a second step of 40s and third step of 52s. It appears in fig 6a that the evaporator enables efficiently the pressurization of the reservoir to nearly 15 bars, which is a sufficiently high pressure to ensure the desalination of the brackish water with a salt concentration of 5 g/l. It appears also that during the first step, which was previously considered to be taken place at a constant pressure, the pressure of the high pressure reservoir decreases showing thus that evaporator is not able to produce a sufficient amount of vapor to maintain the high pressure. The condenser pressure evolution shows that the mass flow rate of the pumped feed water is enough to cool it efficiently and maintain a stable condensing pressure. Permeate production rate over time is also computed to estimate the daily clear water production represented Fig. 7. This intermittent production results in the variation of the pressurized reservoir. When the reservoir pressure becomes lower than the osmotic pressure, then no production of clean water is possible. A full day simulation allowed estimating the daily production to 743 liters of clear water for the process that operated during 136 cycles in the climatic conditions of Perpignan (France) with only 2 m² of solar collectors. The specific thermal energy consumption is evaluated to about 6.7 kWh per m -3 of produced clear water.

Conclusions and perspectives

Solar-driven standalone reverse osmosis (RO) desalination may be a promising and cost effective technology with a high potential for facing water scarcity and meet the water demand at small scale level (5-10m 3 /day) of population living in remote and isolated area. In this objective, there are more and more innovative desalination processes that are developed in the world, especially solar driven systems that are particularly relevant for areas that benefit of an important solar resource and where water is lacking. Moreover, autonomy and robustness are key parameters for developing such kind of solar desalination systems foreseen to be installed in totally off-grid areas.

Directly powering an RO device by low temperature heat, through an integrated thermodynamic engine cycle without energy storage can leads to drastic reduction of installation, maintenance and operating costs in comparison to current implemented PV/RO solution. The solar driven thermohydraulic process described in this study meets these objectives and can be a real alternative for a remote community by its autonomy coupled to low cost hydraulic components. Simple technologybased components have been chosen to ensure low investments and maintenance cost and avoid salt corrosion problems, the only parts that may suffer from corrosion are distributors and valves. The first simulation results performances show that it is possible to produce approximately 750 liters per day from 5 g/l brackish water by such kind installation implementing only 2 m² of solar collectors. These results demonstrate that this innovation may be competitive in comparison to other renewable driven desalination systems. These numerical results correspond to a brackish water application. Deepest simulations have to be led and crossed with experimental results in order to confirm these promising results. A dynamic model of the RO membrane is currently developed and under experimental validation. Next simulations will use a more precise dynamic model and will be extended to seawater desalination.

Fig. 1 :

 1 Fig. 1 : Representation of the different steps during a half-cycle process: a) phase 1 :desalination at constant pressure, b) phase 2 :desalination at decreasing pressure, c)phase 3 : working fluid transfer

Fig. 3 :

 3 Fig. 3: schematic RO unit

  .

Fig. 5 :

 5 Fig. 5: Representation of the hydraulic cylinder

6 :

 6 Fig 6.b shows that at the end of phase 2, the discharge tank (tank 2) no longer contains gas, which means that it is full of feed water to be later desalinated. Phase 3 goes on until the feed tank (tank 1) is full of gas.

Fig. 7 :

 7 Fig. 7 : Evolution of the permeate production over time

  PROCEEDINGS OF ECOS 2018 -THE 31 ST INTERNATIONAL CONFERENCE ON EFFICIENCY, COST, OPTIMIZATION, SIMULATION AND ENVIRONMENTAL IMPACT OF ENERGY SYSTEMS JUNE 17-22, 2018, GUIMARÃES, PORTUGAL

		P l					Pressure
	Tank 1	P h	Tank 2		P h		Tank 1
					π
	evaporator				P l		Tank 2
					Phase n	0	1	2	3
							Volume
		RO unit	Permeate	V tank	Tank 2
					V β
	condenser		Brine	P h	V α	Tank 1
					Phase n	0	1	2	3
	P l		P h		Permeate production ṁ p
	Feed water
	15-20 C			Control	0
				valve	Phase n		0	1	2	3
		Brine
		P l					Pressure
	Tank 1	Tank 2	P h	Tank 1
					π
	evaporator				P l	Tank 2
					Phase n 0	1	2	3
							Volume
		RO unit		V tank	Tank 2
					V β
	condenser		Brine	V α	Tank 1
					Phase n	0	1	2	3
							Permeate production
	P l						ṁ p
				Control			0
				valve	Phase n	0	1	2	3
			Brine

  PROCEEDINGS OF ECOS 2018 -THE 31 ST INTERNATIONAL CONFERENCE ON EFFICIENCY, COST, OPTIMIZATION, SIMULATION AND ENVIRONMENTAL IMPACT OF ENERGY SYSTEMS JUNE 17-22, 2018, GUIMARÃES, PORTUGAL

Table 1 :

 1 Simulation parameters

	Description	Value
	Working fluid	Butane
	Feed water salinity	5 g.dm -3
	Feed water temperature	18°C
	Feed tanks volume	25 liters
	Membrane water permeability	2.7.10 -12 m 3 .s -1 .m -2 .Pa -1
	Membrane salt rejection	99%
	Water pumping pressure	0.9 bar
	Fig.

Volumetric permeate flow (m 3 .s -1 )

  PROCEEDINGS OF ECOS 2018 -THE 31 ST INTERNATIONAL CONFERENCE ON EFFICIENCY, COST, OPTIMIZATION, SIMULATION AND ENVIRONMENTAL IMPACT OF ENERGY SYSTEMS JUNE 17-22, 2018, GUIMARÃES, PORTUGAL

	1.6E-04 Nomenclature
	A C Cp h I L w L s M . m N n P R R s T	2.0E-05 4.0E-05 6.0E-05 8.0E-05 1.0E-04 1.2E-04 1.4E-04 concentration, kg/kg 300 area, m² mass specific heat, J/(kg K) 400 molar enthalpy, J/mol solar irradiance, W/m² water membrane permeability, m 3 /(m 2 s Pa) 500 Time (s) salt membrane permeability, m 3 /(m² s) mass, kg mass flow rate, kg/s molar quantity, mol molar flux, mol/s pressure, Pa ideal gas constant, J/(mol K) membrane salt rejection temperature, °C	600 Phase 2 Subscripts and superscripts 700 Phase 1 Phase 3 amb ambient conditions b brine c solar collector cond condenser ev evaporator f feed water g gas h high in input l low liq liquid m membrane out output p permeate
	T	collector mean temperature, °C	pump pumped
	U V	internal energy, J volume, m 3	s tank tank salt
	Greek symbols	v	vapour or control valve
	α	phase 1
	β	phase 2
	ε	heat exchanger effectiveness
	η	solar collector efficiency
	ϕ	incident heat flux, W
	π	osmotic pressure, Pa

Acknowledgments

The authors wish to acknowledge the LABEX Solstice laboratory for its financial support in this project. This work was supported by the Programme "Investissements d'avenir" (Investment for the Future) of the Agence Nationale de la Recherche (National Agency for Research) of the French State under award number ANR-10-LABX-22-01-SOLSTICE.


Nicolas Clausolles

Pauline Collon

Guillaume Caumon

Generating variable shapes of salt geobodies from seismic images and prior geological knowledge

Keywords:

We have developed an implicit method to automatically generate several possible models of salt top surfaces with varying geometries and topologies. This method can be conditioned to available data such as well markers and seismic picks. As seismic imaging of salt is prone to velocity uncertainty and Fresnel zone effects, the input of the method is a seismic image that is segmented into three regions: salt, sediments and uncertain. The uncertain region contains the salt boundary and all the further computations focus within this zone. A monotonic scalar field, ranging from zero on the edge of the certain salt body to one on the maximal possible salt boundary, defines the cumulated probability for a point to be outside the salt body. A random scalar field, also bounded between zero and one, is then used as a threshold for the first scalar field. The salt boundary is implicitly defined by the zero isovalue of the difference between the two fields, and can be further extracted using marching cubes. The random field parameters have a geometrical and topological impact on the simulated salt volumes. They can be adapted to reproduce specific geological features, but their inference remains difficult. The application of the method to the reconstruction of diapir boundaries from a set of partly interpreted sections shows that it is well suited to honor numerous constraining data while proposing different possible structural scenarios in the uninterpreted parts.

I

Recent advances in seismic acquisition methods and processing algorithms enable the imaging of complex geological settings [START_REF] Jones | Seismic imaging in and around salt bodies[END_REF]. The increasing size and quality of the new seismic data sets have motivated numerous works in seismic interpretation, and especially on the development of automatic interpretation methods. In this global context, a specific attention has been paid to salt tectonics contexts. Salt structures are known, however, to often present challenges for seismic imaging. These challenges stem from three reasons: (1) specific physical properties of evaporites (high velocity and strong acoustic impedance contrasts with surrounding sediments [e.g., [START_REF] Fossen | Structural geology[END_REF]), (2) internal heterogeneity of salt bodies [e.g., Jackson and Lewis, 2012], and (3) complex geometries and topologies that develop when salt and surrounding sediments undergo deformation [e.g., [START_REF] Hudec | Terra infirma: Understanding salt tectonics[END_REF]. This leads to a variety of difficulties for seismic imaging, such as poorly illuminated zones, complex travel paths or anisotropy representation [START_REF] Jones | Seismic imaging in and around salt bodies[END_REF], which raises many challenges and uncertainties during structural interpretation [START_REF] Dellinger | The Garden Banks model experience[END_REF].

Salt interpretation is usually done in two complementary (but not mandatory) steps: the interpreter first tries to highlight the salt geobodies (or their boundaries) in the seismic data, and then extracts the boundaries by creating three dimensional surfaces. The most common way to highlight salt is to compute seismic attributes that are sensitive to the image texture. Some authors [e.g., [START_REF] Gao | Volume texture extraction for 3D seismic visualization and interpretation[END_REF][START_REF] Berthelot | Texture attributes for detection of salt[END_REF][START_REF] Di | Multi-attribute k-means clustering for saltboundary delineation from three-dimensional seismic data[END_REF] use pixel-based statistical descriptions such as the Gray-Level Co-occurrence Matrix [START_REF] Haralick | Textural features for image classification[END_REF] or the Gradient of Texture [START_REF] Shafiq | A texture-based interpretation workflow with application to delineating salt domes[END_REF]. Other authors focus on character-izing the local feature orientations. This characterization can be achieved by computing the image gradient using either the gradient covariance matrix [START_REF] Berthelot | Salt diapirs without well defined boundaries -a feasibility study of semi-automatic detection[END_REF] or the structure tensor [Wu, 2016], whose eigen decomposition provides a measure of the local image anisotropy. Similarly, the reflection dip provides information on the local strata orientation [e.g., Lomask et al., 2007, Halpert and[START_REF] Halpert | Salt body segmentation with dip and frequency attributes[END_REF]. It can be computed using numerous methods, such as plane-wave destruction filters [START_REF] Claerbout | Earth soundings analysis: Processing versus inversion[END_REF] or the ones previously mentioned for the gradient computation. A third alternative to pixel-based statistical descriptions and local orientations characterization consists in characterizing the salt boundaries by solving an edge detection problem [e.g., Aqrawi et al., 2011, Amin and[START_REF] Amin | A new approach for salt dome detection using a 3D multidirectional edge detector[END_REF]. As one single attribute may not be sufficient to highlight the salt boundary in the whole seismic cube, many authors now propose to combine a set of attributes [e.g., [START_REF] Halpert | Salt body segmentation with dip and frequency attributes[END_REF][START_REF] Berthelot | Texture attributes for detection of salt[END_REF][START_REF] Di | Multi-attribute k-means clustering for saltboundary delineation from three-dimensional seismic data[END_REF].

Once the choice of the input attribute(s) has been made, the interpretation of the salt boundary can be seen as a general image segmentation problem. Image segmentation consists in partitioning an input image into different subsets of pixels according to some similarity criterion. In other words, it aims at assigning a label to each image pixel, pixels sharing the same label sharing also some characteristics. This problem has been thoroughly studied in the fields of computer vision and image processing, and numerous methods have been adapted to extract salt boundaries. [START_REF] Lomask | Application of image segmentation to tracking 3D salt boundaries[END_REF], [START_REF] Halpert | Salt body segmentation with dip and frequency attributes[END_REF] use normalized cuts image segmentation [START_REF] Shi | Normalized cuts and image segmentation[END_REF] which is a graph partitioning-based approach. Other authors use methods based on partial differential equations (PDEs) to fit a surface to the salt boundary. [START_REF] Haukås | Automated salt body extraction from seismic data using the level set method[END_REF][START_REF] Haukås | Automated salt interpretation, Part II : Smooth surface wrapping of volume attribute[END_REF] use methods inspired from active-contour models [START_REF] Kass | Snakes: Active contour models[END_REF] and level sets [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] All these methods aim at automatically determining the "best" possible predicate for the salt boundary. But as noted by [START_REF] Long | A comparative study of texture attributes for characterizing subsurface structures in seismic volumes. Interpretation[END_REF], the interpreter should keep in mind that such automatic labeling process can only lead to a first approximation of the salt boundary and not to the exact surface. Generating a single deterministic salt geometry can lead to biases in velocity modeling, hydrocarbon trap determination, and reservoir gross rock volume assessment. The topic of seismic interpretation uncertainties has been raised for a long time in classic stratigraphic contexts [e.g., [START_REF] Abrahamsen | Bayesian kriging for seismic depth conversion of a multi-layer reservoir[END_REF], Thore et al., 2002[START_REF] Bond | What do you think this is? "Conceptual uncertainty" in geoscience interpretation[END_REF], however only a few papers broach this topic and propose solutions to handle it in the presence of salt. [START_REF] Rojo | The use of seismic attributes to enhance imaging of salt structures in the Barents Sea[END_REF] discuss the impact of the attribute selection and tuning on seismic noise reduction along salt boundaries. Most authors also propose to integrate manually picked information to guide the interpretation, either directly [e.g., Wu, 2016] or during the training phase for methods relying on classification algorithms [e.g., [START_REF] Berthelot | Texture attributes for detection of salt[END_REF][START_REF] Di | Multi-attribute k-means clustering for saltboundary delineation from three-dimensional seismic data[END_REF], Waldeland et al., 2018[START_REF] Shi | SaltSeg: Automatic 3D salt segmentation using a deep convolutional neural network[END_REF]. Finally, [START_REF] Haukås | Automated salt interpretation, Part II : Smooth surface wrapping of volume attribute[END_REF] propose an interesting approach to avoid overinterpreting seismic images: they extract surface patches only in well-defined boundary zones, and then smoothly connect the patches in the poorly resolved image parts. Such approach remains, however, deterministic as one single interpretation is produced in the end.

To tackle the problem of interpretation bias reduction, we propose a general framework for characterizing the uncertainties underlying a seismic image. In practice, it allows for the exploration of different structural scenarios that can be admissible considering the imaging uncertainties. This framework is designed in two parts: the first one consists in identifying the image parts in which the interpretation is ambiguous, and the second one focuses on the generation of salt body boundaries with both variable geometries and topologies. The first part is, therefore, specific to seismic interpretation uncertainty characterization. It has already been investigated by several authors, and numerous works on automatic interpretation can be adapted to address this question (see above references). As it is not the main topic of this paper, we confine the presentation to the strategy we retained to model this uncertainty. The second part of the framework, which represents the main contribution of this paper, is more general and can be applied to the modeling of any type of geobodies (salt geobodies, ore geobodies, igneous intrusions, etc.). It addresses the problem of generating complex geobody boundaries while allowing for the integration of data (e.g., seismic picks and well markers) and prior geological knowledge (e.g., type and connectivity of the structures to interpret).

The paper is organized as follows. We first introduce the workflow on a 2D seismic section to provide the reader with a broad overview of the method and the different elements that take part in the modeling process. We then present various simulation results obtained from 3D synthetic data to illustrate the potential of the method, the impact of the different simulation parameters, and data conditioning. We finally discuss the topics of uncertainty characterization on seismic images, integration of prior geological knowledge and parameter inference from available data.

M

The most common approach for sampling structural uncertainties in geomodeling consists of two steps [e.g., [START_REF] Abrahamsen | Bayesian kriging for seismic depth conversion of a multi-layer reservoir[END_REF][START_REF] Lecour | Modelling of stochastic faults and fault networks in a structural uncertainty study[END_REF]. First, an explicit reference structural model is built. This reference model can be seen as the "best" possible model considering the limited amount of available data. Then, the different elements of the model (usually horizons and faults) are perturbed around their reference position to sample different possible geometries. This approach presents, however, some limitations. As noted by [START_REF] Caumon | Elements for Stochastic Structural Perturbation of Stratigraphic Models[END_REF], perturbing an explicit reference model only allows for exploring geometrical uncertainties and not topological ones. Moreover, perturbing complex surfaces (such as salt geobody boundaries) introduces the risk of generating self-intersecting surfaces [Wellmann and [START_REF] Irakarama | Finite Difference Implicit Modeling of Geological Structures[END_REF].

To overcome these limitations, we adopt an implicit modeling strategy. Implicit methods are now well-known in geomodeling and have been widely used in a variety of applications. The seismic grid directly serves as modeling support. The geobody boundaries are represented as a level set of a scalar field [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF]. It prevents classic topological pitfalls, such as self-intersections and invalid contacts, encountered when perturbing explicit models [Wellmann and [START_REF] Irakarama | Finite Difference Implicit Modeling of Geological Structures[END_REF]. It also offers a frame to incorporate conceptual knowledge, and especially locally varying anisotropy, into the modeling [e.g., [START_REF] Martin | Iterative refinement of implicit boundary models for improved geological feature reproduction[END_REF]. This last point is addressed later in the paper.

Characterizing the uncertainties

The input of our workflow is a seismic image (Fig. 1.a). The first step consists in identifying the uncertainties in this image. The idea is to discriminate the image parts in which we consider the interpretation as sure from those which require more focus. In practice, we propose to classify each pixel of the image into one of the following three categories (Fig. 1.b): "sediments" (in red), "salt" (in blue), and "uncertain" (uncolored). This type of partitioning has also been used by [START_REF] Haukås | Automated salt interpretation, Part II : Smooth surface wrapping of volume attribute[END_REF] to avoid over-interpreting the image parts in which the seismic signal is too ambiguous. In our application, the pixels classified as "uncertain" form a volume (referred to as "uncertainty envelope") between the "sediments" and "salt" regions that contains the salt boundary.

In the synthetic examples we present after, the uncertainty envelope is manually defined as we do not have associated seismic images. More details about existing methods that can help to automatically generate the uncertainty envelope are provided in the discussion.

Generating salt boundary interpretations

The next step is to generate salt boundary interpretations in the uncertainty envelope. As stated previously, an interpretation is implicitly defined by a level-set of a scalar field. The generation of this scalar field is inspired by the Object-Distance Simulation (ODSIM) method [START_REF] Henrion | ODSIM: An Object-Distance Simulation Method for Conditioning Complex Natural Structures[END_REF][START_REF] Rongier | Simulation of 3D karst conduits with an object-distance based method integrating geological knowledge[END_REF]]. An interpretation is defined by the combination of two scalar fields: a pseudo-distance field D (Fig. 1.c) and a spatially correlated random field ϕ (Fig. 1.d). The main difference with the ODSIM method is that we perturb the probability of being in the salt body rather than the actual distance to the salt. This approach offers a way to efficiently explore the solution space, and allows for a control of the geometry and topology of the generated boundaries while honoring local salt boundary observations.

The pseudo-distance field D (Fig. 1.c) is a monotonically varying scalar field that ranges from 0 at the contact with the "salt" region and the base of sediments, to 1 at the contact with the external "sediments" region. It represents a relative normalized distance of each sample to salt. In other words, it can be seen as the probability for each sample to be located outside the salt boundary (i.e., to be sediments). This scalar field is obtained by: (1) setting to 0 the value of the samples in contact with the "salt" region and the base of sediments, (2) setting to 1 the value of the samples in contact with the "sediments" region, and (3) interpolating the intermediate values in the remainder of the "uncertain" region. In this paper, we use the finite difference-based interpolation method proposed by [START_REF] Irakarama | Finite Difference Implicit Modeling of Geological Structures[END_REF] to generate this scalar field.

The second scalar field ϕ (Fig. 1.d) is a spatially correlated random field which is applied to the pseudo-distance field D. The values of ϕ must lie in the same range as the values of D (between 0 and 1). Two parameters control the random field generation: a distribution model and a spatial correlation model. The distribution model describes the repartition of the values that will be drawn (and must therefore be bounded by 0 and 1). The spatial correlation model describes the local variability of the random field. The choice of the random field parameters is mainly expert-driven. In this paper, we use a sequential Gaussian simulation (SGS) to generate this random field ϕ. The spatial correlation model is thus represented by a variogram model.

Once the two scalar fields D and ϕ have been generated, we define the salt boundary S B (Fig. 1.e) as the location of the 0-level surface of the scalar field defined as the difference between D and ϕ:

S B : {(t, x, y) | D(t, x, y) -ϕ(t, x, y) = 0}.
As we want to sample structural uncertainties, we need to generate many realizations. In our case, all the realizations share the same pseudo-distance field D. In order to simulate different salt boundaries, we just need to generate one random field per realization. Figure 1.f presents three different salt boundaries generated using the same simulation parameters (i.e., distribution and variogram models).

S

The above workflow is based on several elements: the uncertainty envelope, the pseudo-distance field D and the random field ϕ. Each of these elements requires the definition of some parameters, which raises several questions. How do parameters impact the realizations? Can we introduce conceptual knowledge into the modeling to reproduce specific geological features? Can we infer parameters from data and prior knowledge? And overall, can we appraise the different interpretations we simulate? The remainder of the paper aims at providing elements to answer these questions. In this section, we use an example of 3D synthetic salt diapir to illustrate the impact of the simulation parameters on the geometry and topology of the salt boundaries. Using the same example, we then present two examples illustrating how conceptual knowledge can be introduced by an appropriate selection of parameters to reproduce specific geological features. Finally, we explain how local salt observations or interpretations can be integrated in the process.

Generating variable salt shapes

Figure 2 illustrates the generation of different salt boundary geometries. In this example, we consider that the segmentation of the image resulted in a clearly defined salt diapir (below the blue surface in Fig. 2.a). The uncertainties are thus limited to the position of the boundary in the envelope. The distance field D is computed in the uncertainty envelope (Fig. 2.a). The SGS parameters used to generate the random field ϕ are: (1) a triangular distribution model ranging from 0 to 1, with a mode at 0.5, and (2) an isotropic Gaussian variogram model with a range of 800 meters. These parameters are used as "default" parameters in our applications as they allow for the generation of relatively smooth boundaries (the global model size is about 3000 meters) that are centered within the uncertainty envelope. It is important to note that the generation of the random field ϕ is intrinsically prone to statistical biases, due to the relatively large variogram ranges as compared to the uncertainty envelope dimensions. Several strategies can therefore be considered to generate it. The one we selected consists in simulating the random field in the whole seismic image, and keep only the sub-sample corresponding to the uncertainty envelope. Three realizations are presented in Figs.

2.b, 2.c, and 2.d. Depending on the random field values, the simulated boundaries are generated everywhere in the uncertainty envelope and can be either close to salt (blue surface parts) or sediments (red surface parts).

Impact of the simulation parameters

The impact of the parameters used for the simulation of the random field ϕ is relatively intuitive and can be described in terms of oscillations of the simulated boundary around a reference position. As a reminder, two parameters control the generation of ϕ: a distribution model and a variogram model.

The mean of the distribution model determines the reference position of the boundary in the envelope. If the mean is close to zero, the boundary will be statistically close to the "salt" region and, conversely, if the mean is close to one, the boundary will be statistically close to the "sediments" region. It controls therefore the total volume of salt that is simulated. For a given spatial correlation model, the standard deviation of the distribution determines the amplitude of the boundary oscillations as compared to its reference position: the higher the standard deviation, the larger the amplitude of the oscillations. The principal ranges of the variogram model control the frequency of the oscillations around the reference position, and thus the smoothness of the boundary: the larger the ranges, the smoother the boundary. Very short variogram ranges may also trigger the apparition of isolated salt blobs. Figure 3 summarizes the impact of the parameters on the boundary geometry.

It is important to note that the uncertainty envelope, the spatial correlation model and the distribution model all affect the boundary smoothness: the envelope shape controls the global shape of the distance field D (and thus the reference boundary position and its initial shape), the spatial correlation model range determines the global wavelength of the boundary oscillations within this envelope, and the distribution model standard deviation determines how much the boundary can move away from its reference position.

Generating more complex salt shapes

In geology, it is common to have some prior knowledge about the structures being interpreted. Such knowledge can be inferred from many sources of information, such as available data, global geological context, known field analogs, etc. Therefore, an interpreter may want to tune the simulation parameters in order to reproduce specific features. We present here two examples to illustrate how to select the simulation parameters according to a given target salt boundary shape.

Halokinetic sequences are sedimentary objects that are often associated with salt diapirs. These thin packages of strata are folded upward during the diapir growth, which results in nonconformable relationships between the different packages of strata [Giles and Rowan, 2012], and typical cusp-shaped patterns for the salt boundary. Salt Christmas trees are diapirs with very specific shapes. They are characterized by allochthonous salt "wings" located

Generating variable salt topologies

For the moment, we have only considered the generation of variable boundary geometries. Depending on the configuration of the regions interpreted during the image segmentation phase, it is also possible to generate salt boundaries with variable topologies. In particular, if the "salt" region is made of one basal part and a top teardrop-like part, then both connected and disconnected salt boundaries can be generated. This is illustrated in Fig. 6.a: the "salt" region contains a basal layer part, and a salt blob has been isolated above. The intermediate part has, however, been classified as uncertain, which allows for the interpretation of either a diapir stem or a weld. To characterize the topological changes between realizations, the approach of Thiele et al.

[2016] could be use. They propose a framework for characterizing the topological relationships between the different volume entities of a geomodel. In particular, they distinguish between structural and lithological topologies which describe the relations between respec-tively the different (structural) volumes of conformable stratigraphic units, and the different lithologies (independently of their spatial location). This classification is not perfectly suited for salt tectonics, as diapirs usually have a limited extent at the model scale (they do not cross the whole model) and halokinetic sequences define complex stratigraphic relations (they can be in the same time conformable in the basin center and unconformable in the diapir vicinity). We consider, however, that it provides an interesting way of describing the different topologies that our method can generate. Our method, in itself, can generate variable structural topologies (depending on whether the simulated bodies are connected or not). It does not generate variable stratigraphic topologies, as we do not model the sediments. Modeling halokinetic sequences would be, however, an interesting extension of the method and would allow for the generation of variable stratigraphic topologies. Indeed, as the salt boundary geometry varies from one realization to another, modeling different disconnected salt bodies implies modeling different welds, and thus modifying the lithological contacts through the weld.

Conditioning to available data

The method proposed in this paper has the advantage of resting on a simple principle. This makes it easy to adjust to more specific purposes. Especially, it can be constrained to honor various types of geological data. The data conditioning is ensured by constraining the random field generation.

The easiest type of data to handle is punctual information on the salt boundary location (e.g., well markers or deterministic seismic picks indicating the boundary). As the salt boundary is defined by the set of points for which the value of the pseudodistance field D is equal to the value of the random field ϕ, such observations can directly be treated as conditioning data for the SGS (by setting at data locations the value of the random field ϕ to the value of the normalized distance field D). Figure 7 illustrates this conditioning. Two well markers (Fig. 7.a) are added to the example presented in Fig. 6. The simulated salt boundaries (Fig. 7.b, 7.c and 7.d) exhibit variable geometries, but honor the two markers.

This type of conditioning does not restrict the number of data points that can be taken into account. Our method can, therefore, be applied in cases where some image parts have been interpreted, to reconstruct the boundary in the missing image parts. Figure 8 presents an example that illustrates this application. We consider that a seismic interpreter has made a fast interpretation of a salt diapir. The interpreter picked the salt boundary every twenty lines, in the regions where it is clearly defined (Fig. 8.a). The top of the diapir and the top of the salt layer away from the diapir, which are generally well imaged, are therefore already interpreted. The diapir flanks are, however, poorly resolved and remain to interpret. Figures 8.b, 8.c, and 8.d present three automatically generated interpretations of the possible diapir boundary, which all honor the manually picked samples but exhibit very variable geometries and topologies.

We can also imagine integrating lithological observations (e.g., punctual well log information like "salt" or "sediments"). We do not illustrate this point here, but Henrion et al. [2010] showed that it only requires the use of a more sophisticated sampling method than direct sampling. They use Gibbs sampling [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF] in their application. The main idea is that the simulated value at conditioning data location should be in accordance with the observation: if we observed salt, the simulated value should satisfy the condition Dϕ < 0, and conversely if we are in sediments (then Dϕ > 0). The Gibbs sampling ensures that the simulated data still verify the spatial correlation model.

D

In this section, we develop some aspects that have been mentioned above. In particular, (1) we discuss the problem of the definition of the uncertainty envelope and the pseudo-distance field D, and their impact on the topology of the simulated boundaries; (2) we summarize the different elements about integration of prior knowledge and geological concepts in the modeling; and (3) we develop some points about parameter inference.

Workflow applications and anisotropy modeling

A common question that arises when modeling geological structures is the question of anisotropy modeling. Concerning salt modeling, there are two types of anisotropy: the anisotropy related to the large-scale structures (e.g., salt domes, walls, sheets, etc.) and the anisotropy related to the local features of the salt boundary (e.g., smooth, cusp-shaped, with wings, etc.). It is important to note that, depending on the application, we are not interested in integrating the same level of details into the models. As an example, in geomodeling we are interested in building models that integrate as much information as possible (such as precise halokinetic sequence contacts) in order to make predictions on potential hydrocarbon traps, reservoir heterogeneity, well planning, etc. In seismic imaging, using velocity models with such small scale details may, however, lead to more distorted images than using smoother velocity models, as they introduce details that drastically complicate the travel paths. We propose therefore to integrate separately structurerelated anisotropy and boundary-related anisotropy into the model. Structure-related anisotropy is handled during the definition of the uncertainty envelope, and boundary anisotropy is controlled by the variogram model.

Generation of the uncertainty envelope and influence on the simulated boundaries

As previously discussed, the topology of the uncertainty envelope impacts the topology of the simulated salt bodies. As a reminder, if the initial image segmentation phase produced several independent "salt" regions, then the simulation may produce one or several distinct geobodies. Otherwise, the simulation will produce one single geobody (provided that the variogram principal ranges are large enough to prevent the apparition of isolated blobs). The shape of the uncertainty envelope also has a direct impact on the geometry of the simulated boundaries, as it controls the shape of the pseudo-distance field D. It introduces a trend in the interpretations as it defines the "reference" boundary position together with the distribution model mean. We estimate, however, that this trend mainly determines the type of structures that will be interpreted (e.g., salt diapirs, salt walls or allochthonous salt sheets). When interpreting seismic images, most of the uncertainties relate to the exact location of the boundaries or the potential salt body connectivities, rather than the type of salt body to interpret. We consider, therefore, that the uncertainty envelope should be as smooth as possible so most of the geometrical properties of the simulated boundary are determined by the simulation parameters rather than the pseudo-distance field D itself.

As it determines the reference position of the boundary, the definition of the uncertainty envelope is a crucial step. Several strategies can be considered to achieve this step. A first solution consists in coupling an automatic interpretation method proposed in the literature (see the introduction review) with a confidence criterion. The most straightforward method could be to use the salt indicator function proposed by Wu [2016] and to define the three different categories by thresholding it. As an example, [START_REF] Haukås | Automated salt interpretation, Part II : Smooth surface wrapping of volume attribute[END_REF] use this type of approach. If enough data points are available (either lithological or boundary observations), it is also possible to consider building the uncertainty envelope directly from these data, using for example the method proposed by [START_REF] Martin | Iterative refinement of implicit boundary models for improved geological feature reproduction[END_REF]. It would provide a good way to integrate locally varying structural anisotropy. Finally, we can also consider applying velocity uncertainty to this envelope [Thore et al., 2002], or re-migrating it [Osypov et al., 2013, Fomel and[START_REF] Fomel | Structural uncertainty of time-migrated seismic images[END_REF] to account for imaging parameters uncertainty. To go further, it could be interesting to consider methods from the surface reconstruction domain in computer graphics to provide control on the topology of the uncertainty envelope [e.g., [START_REF] Huang | Topology-controlled Reconstruction of Multi-labelled Domains from Cross-sections[END_REF].

Despite the above, it may not be easy to generate a smooth uncertainty envelope from automatic image segmentation techniques. To handle envelopes that would be too irregular, we suggest to use a second envelope that encompasses the first one (thus larger, but smoother). This second envelope is used to generate a smooth distance field D, and the points of the second envelope that do not belong to the first one (i.e., those that were first classified as "sediments" and then as "uncertain") can be used as lithological constraints (using a Gibbs sampler) to ensure that the simulated boundaries remain within the initial uncertainty envelope.

Another difficulty a user might face in some configurations consists in generating a pseudo-distance field D allowing for the interpretation of variable salt topologies. As an example, in Fig. 1 we had to constrain the values of the field D within the diapir stem in order to generate both connected and disconnected salt bodies. Indeed, as we impose its value to 1 on each side of the stem (at the contact with the "sediments" region), the interpolation populates the stem with values extremely close to 1, which makes the probability of simulating a connected salt diapir almost null. To constrain the values within the stem, we manually picked several samples within the stem and imposed the distance field values to be equal to 0.5 at the sample location. It is possible to automate this conditioning, for example by extracting the medial axis of the envelope and imposing its value.

To go further, it is possible to impose the topology of the simulated bodies. A potential strategy is to analyze the topology of the pseudo-distance field D to detect the presence of critical points. These points are characterized by a vanishing gradient of D, and identify the locations (and the associated values of D) where topological changes will occur when extracting a level set (i.e., when moving from a single to disconnected bodies). The values of ϕ can then be adapted to impose a given topology: if ϕ < D at a critical point, disconnected salt bodies are simulated, and conversely if ϕ > D a single salt body is simulated.

Integrating prior knowledge and geological concepts

Knowledge can be integrated either in the uncertainty envelope definition or in simulation process, depending on the type of information. The uncertainty envelope constrains the global shape of the structures. Its definition offers a way to integrate information about regional context, such as the regional structure orientations, the type of structures and their connectivity. The distribution model defines the average volume of salt that is generated, and how much the boundary can move away from its reference position in the uncertainty envelope. The variogram model offers a control on the local shape of the boundary. Depending on the application, a user might choose a long range isotropic variogram (to generate smooth boundaries, Fig. 2), or an anisotropic one to reproduce specific patterns (Fig. 4). Locally varying anisotropy can be introduced (and especially varying orientations) to reproduce, for example, depositional or salt spreading processes (Fig. 5). The choice of the variogram model also determines the influence of punctual data (e.g., well markers ans seismic picks) in the simulation. Choosing a Gaussian model (which is characterized by a large correlation for short distances) ensures that the interpretation in the neighborhood of the data points is consistent with the primary information. The selection of the variogram ranges then offers a way to control the area of influence of the data points.

To go further, it would be interesting to integrate geological concepts in the modeling. The most direct example is maybe the modeling of salt welds. When the simulation produces disconnected salt bodies, we know from a geological point of view that an allochthonous salt body cannot "float" among younger sediments. It has to have been connected to a former salt source through a feeder, whose existence is indicated by the presence of a weld. Modeling welds should be of interest in the oil and gas industry. They play indeed a major role in basin compartmentalization and in the deformation processes that occur at basin scale, and they can severely affect seismic imaging [START_REF] Jackson | Salt tectonics: Principles and practice[END_REF].

Simulation parameter inference

The definition of the simulation parameters is a key step of the workflow. It should take into account both the geometrical and topological specificities of the salt boundaries to interpret, when such characteristics can usually not be observed on seismic images. Therefore, strategies to calibrate the model parameters need to be developed depending on the type of salt geobodies to be simulated. [START_REF] Henrion | ODSIM: An Object-Distance Simulation Method for Conditioning Complex Natural Structures[END_REF] discuss some aspects of such model parameter calibration using prior geological knowledge. Among the available sources of information, well data and manual seismic picks can be used to obtain information about the actual salt distribution, geometry and topology. Analogs such as lab experiments [e.g., Vendeville and Jackson, 1992, [START_REF] Dooley | Inflation and deflation of deeply buried salt stocks during lateral shortening[END_REF], outcrop descriptions [e.g., [START_REF] Ringenbach | Salt tectonics in the Sivas Basin, Turkey: outstanding seismic analogues from outcrops[END_REF] or interpreted seismic data sets from similar geological contexts usually represent another source of information, as they can generally be used together with field data to derive some set of model parameters [START_REF] Guardiano | Multivariate geostatistics: beyond bivariate moments[END_REF]. The use of analogs is, however, not directly possible in our method.

This can be explained by considering the mathematical aspects of our problem. Let D per t be the scalar field defined as D per t = Dϕ. From a geostatistical point of view, we have three random fields (D , D per t , and ϕ) related by a linear combination. Considering that D and ϕ are independent and intrinsically second order stationary, we can relate their spatial correlation models: γ D = γ D per t + γ ϕ . This means that, knowing the variogram models γ D and γ D per t of D and D per t , we can deduce a variogram model γ ϕ for the simulated field ϕ. As we know the pseudo-distance field D, we can compute an experimental variogram from it and define an analytical variogram model γ D . The field D per t , whose variogram model represents the true variability of the salt boundary, is unknown yet (as the field ϕ is not known). We might consider using some analog to deduce an analytical model γ D per t . For the analog spatial correlation model (and thus γ D per t ) being consistent with γ D , we shall define an uncertainty envelope on the analog. Yet, this is not possible on lab and outcrop analogs. On analog seismic data sets, we can generate an uncertainty envelope. It is, however, not straightforward to estimate the analog boundary variability relatively to the envelope, as we only know the zero-level of the analog perturbed distance field, and not this whole field. Moreover, the fields D and D per t are non-stationary and exhibit zonal anisotropy, with varying local principal orientations. In practice, it is, therefore, complex to define variogram models that actually describe these fields. Finally, the spatial distortion effect induced by the varying width of the uncertainty envelope adds another degree of complexity to the variogram model definition.

C

The proposed workflow allows for the generation of stochastic interpretations of salt geobody boundaries reflecting seismic imaging uncertainties. It can simulate both connected and disconnected salt bodies while honoring preexisting geological observations, such as manual seismic picks, well markers and lithological observations. The parameters of the method can be calibrated to reproduce specific geological features; however, quantitative parameter inference remains difficult. To address the impact of geometrical and topological salt interpretation uncertainties, further work should focus on appraising the various interpretations we generate with our method to inject them back into the seismic processing chain.

Moreover, the workflow is generic and can be used in other types of applications. In particular, it is well suited for the generation of any type of geological body boundaries (salt bodies, ore bodies, igneous intrusions, etc.), even in the presence of nu-merous constraining data. As an example, it could be of great use for estimating ore deposit extent in mining applications, where numerous data (and especially drillholes) are available.

A

This work was performed in the frame of the RING project (http://ring.georessources.univ-lorraine.fr) at Université de Lorraine. We would like to thank for their support the industrial and academic sponsors of the RING-GOCAD Consortium managed by ASGA. We also acknowledge Paradigm-Emerson for the SKUA-Gocad software and development kit. The softwares corresponding to this paper are available as the Goscope plugin of SKUA-Gocad and the stand-alone SIGMA library. We would finally like to thank the associate editor B. Wallet, and E. de Kemp, X. Wu and two anonymous reviewers for their work and comments that greatly improved the paper.

Figure 1

 1 Figure 1 Simulating stochastic interpretations of salt boundaries. The input seismic data [a, image from Jackson and Lewis, 2012] are segmented into three categories (b): sediments (red), salt (blue) and uncertain (uncolored). A normalized "distance" field D (c) and a normalized random field ϕ (d) are then generated in the uncertain region. The salt boundary S B (e, dashed curve) is obtained by extracting the zero isovalue of the difference between D and ϕ (e). Different salt boundaries (f) can be obtained by generating different random fields (the three realizations share the same distribution and variogram models).

  Figure 4.a illustrates an example of such an interpreted boundary from Giles and Rowan [2012]. To reproduce these features, we propose to use an anisotropic variogram model, with horizontal ranges higher than the vertical range (Fig. 4.b). The shorter vertical range reduces the vertical correlation and allows the boundary position to vary more rapidly within the envelop, reproducing the cusp-shaped patterns.

Figure 2

 2 Figure 2 Simulating variable salt top boundary geometries. (a) Input data: the "salt" region (in blue), the "sediments" region (in red) and the "uncertain" region (painted with the pseudo-distance field D). The dashed line represents the contact between the external sediments and the uncertainty envelope on the sections. (b), (c) and (d) Three different simulated salt boundaries painted with their respective random fields ϕ. ϕ follows a triangular distribution, and the variogram model γ is isotropic.

Figure 3

 3 Figure 3 Impact of simulation parameters on the boundary geometry. Each sinusoid represents an idealized boundary within the uncertainty envelope. (a) The distribution model affects the mean position of the boundary. (b) The variogram ranges affect the smoothness of the boundary. The variogram's principal ranges linearly increase from the left to the right of the uncertainty envelope. Input data from Fig. 2.a.

Figure 4

 4 Figure 4 Reproducing halokinetic sequences contact. (a) Example of the typical cusp-shaped patterns observed at the contact between halokinetic sequences (in red) and a salt diapir (in blue). Image and interpretation from Giles and Rowan [2012]. (b) An example of salt boundary generated using a transversely isotropic variogram model with short vertical range. Input data from Fig. 2.a.

  Figure 5.d presents a realization example of such a Christmas tree salt diapir.

  Figures 6.b, 6.c, and 6.d present three examples of simulated salt boundaries, generated using the same random fields as in Figs. 2.b, 2.c, and 2.d. Depending on the random field values, we generate either connected (Fig. 6.d) or disconnected (Figs. 6.b and 6.c) salt boundaries, which would need to be connected by a weld.

Figure 5

 5 Figure 5 Reproducing Christmas tree-like salt diapirs. Example of a salt Christmas tree diapir wing on uninterpreted (a) and interpreted (b) seismic sections. Image and interpretation from Mohr et al. [2007]. (c) Input data from Fig. 2.a, with the relative geologic time function used to align the variogram axes displayed on the vertical sections. (d) An example of salt boundary generated using a uniform distribution model and a transversely isotropic variogram model, with a short "vertical" range and larger "horizontal" ranges aligned with the relative time function.

Figure 6

 6 Figure 6 Simulating variable salt top boundary topologies. (a) Input data: the "salt" region (in blue), the "sediments" region (in red) and the "uncertain" region (painted with its pseudo-distance field D). The dashed line represents the contact between the external sediments and the uncertainty envelope on the sections. (b), (c) and (d) Three different simulated salt boundaries, painted with their respective random fields ϕ.

Figure 7

 7 Figure 7 Conditioning the simulations to well data. (a) Input data from Fig. 6.a, with a synthetic well and two markers (white diamonds) indicating the salt boundary position. (b), (c) and (d) Three different simulated salt boundaries (painted with their respective random fields ϕ) obtained using the same simulation parameters. The boundaries still have variable geometries, but honor the two well markers.

Figure 8

 8 Figure 8 Sampling variable diapir geometries and topologies while honoring incomplete manual interpretation. (a) Manual seismic picks of the salt boundary (painted with the pseudo-distance field D), used as constraints with the input data from Fig. 6. (b), (c) and (d) Three simulated salt boundaries (painted with their respective random fields ϕ) obtained using the same simulation parameters. The boundaries honor all the seismic picks, but exhibit variable geometries and topologies in the unconstrained part of the diapir.

  ], Wu [2016] solves a screened Poisson equation to compute an implicit salt indicator function and Wu et al. [2018] extract an optimal path by solving the eikonal equation. Finally, several authors have also proposed to use various machine learning algorithms [e.g., Berthelot et al., 2013, Di et al., 2018, Waldeland et al., 2018, Shi et al., 2019].

S. T. Thiele, M. W. Jessell, M. Lindsay, V. Ogarko, J. F. Wellmann, and E. Pakyuz