Advanced well model for superhot and saline geothermal reservoirs

.


Introduction
Magma-driven hydrothermal systems located in volcanic areas can give rise to the formation of superhot/supercritical water resources (Heřmanská et al., 2019;Jolie et al., 2021).These resources are shallow enough to be targeted using current well drilling technology and have a high energy production potential (Axelsson et al., 2014;Friðleifsson et al., 2014;Reinsch et al., 2017).Drilling into such resources remains a risky endeavor given the uncertainties surrounding the prospective well performance.Understanding how superhot resources would respond to well operation requires a combination of reservoir and well modeling.While exhaustive literature is available on simulation of conventional, high-enthalpy resources (O'Sullivan et al., 2001;McDermott et al., 2006;Driesner and Geiger, 2007;Gunnarsson et al., 2011), studies that explicitly represent a magma body as heat source are scarce and have mostly been limited to 2D (Hayba and Ingebritsen, 1997;Scott et al., 2015Scott et al., , 2016Scott et al., , 2017)).Only recently have studies appeared that combine reservoir simulation with an explicitly represented magma body and a simplified handling of wells as local source/sinks (Magnusdottir and Finsterle, 2015;Yapparova et al., 2022a).
Performing coupled well reservoir modeling is a fast and cheap means to address key questions emerging from the prospective to the production stages of a geothermal energy production project.It is an invaluable tool to assess the producibility of the geothermal resource and predict the long-term performance of a potential well.Well simulators have existed for many decades, progressively gaining more complex features.Early models (Chierici et al., 1981;Gould, 1974;Miller, 1979;Ortiz-Ramirez, 1983) were already quite advanced, solving coupled mass, energy, and momentum conservation equations but also had various limitations, such as: restriction to steady state solutions, a single feedzone or assuming a hydrostatic pressure gradient in the well.
Modern well models use implicit, transient, and fully coupled formulations, often solving conservation equations in a well discretized into well segments (i.e., control volumes), see for example: Holmes et al. (1998), Livescu et al. (2010) or Pan and Oldenburg (2014).These so-called multi-segments well models can usually deal with multiple feedzones, variable well inclination and cross-sectional area, and even well branching.Latest advances in well modeling are the inclusion of non-condensable gases like CO 2 (Battistelli et al., 2020), the coupling with geomechanical simulators (Rutqvist et al., 2018) the simulation of well scaling (Khasani et al., 2021;Zolfagharroshan and Khamehchi, 2020) and the simulation of airlifting procedures (Akbar et al., 2016).A review of current well models and their differences in implementations can be found in Tonkin et al. (2021).
Several limitations of these existing models motivated us to develop a new multi-segment well model.Firstly, current models often use simplified or limited equations of states to evaluate the fluids thermodynamic properties.Our model uses thermodynamic tables valid for pressures ranging from 0.1 to 500 MPa, temperatures from 5 • C to 1000 • C and salinity from 0 to 100 wt% NaCl (Driesner, 2007;Driesner and Heinrich, 2007).This allows our model to accurately simulate superhot systems as well as potential halite (i.e., solid NaCl) precipitation.For the latter, we could not find an equivalent feature in well models currently presented in the literature.Salinities in superhot reservoir can be significant, this is for instance the case in Reykjanes (Iceland) where the IDDP-2 well was drilled ( Ármannsson, 2016), this suggest that halite scaling could be an issue.
Secondly, we could not find well models incorporating a moving airwater interface in the well, at the top of the water column.At various stages of the well operating life an air phase can be present in the upper part of the well, for example during cold water injection or when the air pressurization technique is used to start the well.For this reason, our model includes an air phase and incorporates a simplified handling of the displacement of the air-water interface.Neglecting air density, this displacement essentially gives rise to a moving pressure boundary condition within the well.Air pressurization and airlifting are two different techniques that can be used to start a well, with air pressurization technique the air is present above the water column and is pressurized before being quickly depressurized, airlifting, as simulated by Akbar et al. (2016), consist in the injection of air within the water column to produce a low density mixture.
Lastly, the well model is developed for a control volume finite element (CVFE) scheme, which allows for unstructured meshes and provides great flexibility for the geometrical representation of complex reservoir permeability structures and well completions.We use a recent 3D upgrade of the 2D CVFE scheme of Weis et al. (2014) within the CSMP++ simulation platform.Another important advantage of this reservoir model is that we can explicitly represent magmatic intrusions as well as the potential degassing of magmatic volatiles (Lamy-Chappuis et al., 2020a;Scott et al., 2015).Magmatic intrusions act as localized and transient heat sources for the geothermal reservoir.Their explicit inclusion allows us to reproduce the natural state of the reservoir more accurately in magma driven hydrothermal systems.
In this paper we focus on the technical introduction of this new well model and on several generic applications demonstrating the insights that can be generated using its original capabilities.

Conservation equations
The well governing equations cover mass, energy, momentum and salt mass conservation.They are expressed in a single dimension, in the direction of the well.
The primary variables associated with the above equations are the fluid velocity, bulk specific enthalpy, fluid pressure and bulk salt mass fraction (NaCl salt)."Fluid" refers to either a single-phase fluid (liquid or vapor or supercritical) or to a mixture of vapor and liquid, "bulk" refers to the mixture of fluid and solid NaCl (halite).All fluid phases possibly contain dissolved NaCl which can precipitate in the form of halite.
We use a homogeneous flow model; therefore, in the case of multiphase flow, the vapor velocity is equal to the liquid velocity and is referred to below as "fluid velocity".Although upgrading the homogenous flow model to a drift flux model (Zuber and Findlay, 1965) is possible it seems that experimental calibrations of the drift flux model for the fluids and conditions relevant to our simulations are currently lacking and we have, therefore, not pursued this option.
The momentum conservation with the homogeneous flow model is written: where p is the fluid pressure (we do not consider capillarity and, therefore, all fluid components share the same pressure), ∂x is the space increment is in the direction of the well, ρ f is the fluid density, g is the gravity acceleration, θ is the well inclination angle, f is Darcy's friction factor, r eff is the effective radius of the well and u the fluid velocity.The effective well radius reflects a reduction of the flowing section of the well due to halite precipitation and r eff is taken equal to r w (1 − S h ) where S h is the halite saturation and r w the initial well radius.We neglect the air density and use ∂p ∂x = 0 for the air phase.The second term in Eq. ( 1) is the frictional pressure loss and the last two terms are the inertial pressure losses.Following Zigrang and Sylvester (1982), Darcy's friction factor is written as: where ϵ is the roughness height, and Re is the Reynolds number.The Reynolds number is written: where μ β and S β are the dynamic viscosity and saturation of phase β.The term in the denominator is an approximation of the fluid mixture viscosity obtained via the product of the viscosities to the power of the saturations.
For simplicity, below we show the conservation equations pertaining to water saturated segments, the treatment of air saturated or partially air saturated segments is explained in Section 3.
Fluid velocity is derived from fluid mass conservation: We use the following convention: u is positive for downflow and negative for upflow.The fluid source Q f can stem from water injection from the top of the well and/or exchange with the formation.The calculation of this source term is described in Section 2.2.
The bulk energy (fluid + halite) conservation is written: where E contains the internal and kinetic energy of the fluid and halite phases.The transport of energy is only done by the fluid as halite is assumed to precipitate on the pipe's wall and remain immobile.All components of the fluid are assumed to share the same temperature.The first term on the right-hand side is the advection of internal and kinetic energy done by the fluid, the second term is the change in the fluid's gravitational potential energy.The term Q h is a source term stemming from water injection from the top of the well and/or exchange with the formation.
The term Q loss is the conductive heat exchange between the well and the formation, neglecting the contribution of the air phase.The heat loss model currently implemented is similar to the one used in the T2WELL simulator (Pan and Oldenburg, 2014): where κ is the formation thermal conductivity, T well is a well temperature, T ∞ is the ambient formation temperature and α is the formation's thermal dispersivity.The term in the denominator is known as the Ramey heat loss function (Ramey, 1962).We neglect the temperature gradient between the outside casing and the center of the well.
Finally, the salt mass conservation equation: where C NaCl is the total salt content in kg.m − 3 (adding the contributions of the dissolved salt and halite) while salt transport is only done by the fluid.X f is the salt mass fraction of the fluid.Q NaCl is a source term stemming from the exchange with the formation.

Inflow performance relationship
The source terms Q f , Q h , Q NaCl are calculated using an inflow performance relationship based on the Peaceman formula (Peaceman, 1983), for instance for the mass source and for well inflow: For well outflow: where the subscript "res" stands for reservoir, W is the well index, P d is the pressure drawdown, λ β is the relative mobility of the fluid phase β and λ tot is the total relative mobility of the reservoir fluids (sum of the relative mobilities, i.e. λ tot = kr,v ).For well outflow, we use the reservoir total relative mobility rather than the individual mobilities, we also do not use mobilities that would be evaluated in the well.Doing so would lead to erroneous transfer rates.For example: if the well is saturated with vapor and the reservoir with liquid, the outflow of vapor should be limited by the mobility of the liquid phase present in the reservoir.
For the flow in the reservoir, we include a liquid residual saturation (S l,res ) below which the transport contribution from the liquid phase is null; in this study we use S l,res = 0.3 and a simple linear relationship between the phase's relative permeability and their saturations, such that: where S * β = Sβ 1− S h , we can then rewrite λ tot in term of saturations: The well index is calculated via: , where L c is the completion length, k, the formation permeability, S the skin factor (we used S = 0) and r e the equivalent radius.Chen et al. (2003) proposed an expression for the equivalent radius in the case of a 2D CVFE mesh comprised of triangular elements: the equivalent radius is , where A is the area of the control volume associated with the well completion.This expression can be used in a 3D CVFE mesh provided that, for consistency, we now surround the well with a sleeve of prism elements.We noticed that subsequently subdividing each prism into 3 tetrahedra elements would provide essentially identical results.
Inherent approximations associated with this formulation are the assumptions of incompressible, steady state and radial flow near a well.Another possible issue with this model is that with insufficient mesh resolution we would not capture the rapidly changing fluid properties in the vicinity of the well, especially at supercritical conditions.We address these possible limitations in the next section.
Finally, the source term Q h and Q NaCl are readily computed by multiplying the mass transfer rate Q f with the upwinded specific enthalpy and salt mass fraction, respectively (i.e., for these two properties we always use reservoir values for inflow and well values for outflow).

Solution method
The placement of the different variables on the well and reservoir CVFE mesh is shown in Fig. 1.
Well pressure is solved at the nodes; it is then linearly interpolated to the segment's center for the purpose of evaluating the thermodynamic properties of the fluid contained in the segment.The reservoir temperature is located on the mesh nodes, it is then linearly interpolated to the segment centers in order to obtain the reservoir temperature used for the calculation of Q loss .Completions are centered on a mesh node, their Here, 8 tetrahedral elements belonging to the reservoir mesh are represented.Nodes (black dots) are connected by the edges of the elements.Connecting the barycenters of the elements edges, faces and volume, we construct sector volumes.We represented one sector volume of the front facing and darkened tetrahedron, the one associated with the central node.The assemblage of all sector volumes around a given node forms a control volume.Well segments are placed between neighboring nodes, one such segment, labelled "Segment N" is shown, the entire well is composed of connected segments.We compute wellreservoir exchange terms on the bottom (deepest) node of the relevant segments, where both reservoir and well pressures coexist.
B. Lamy-Chappuis et al. lengths can take any value up to the sum of half the lengths of the attached well segments.
The partial differential equations governing the flow in the well are discretized using a time implicit and 1D upwind multi-segment discretization.First, the residuals of the momentum, mass, energy and salt mass conservation equations are evaluated on each segment.In the following equations, N is the segment index, K the time level (values taken at the previous time level are shown within square brackets), L the segment's length and A eff = A(1 − S h ) is the segment's effective cross area that can be modified by halite precipitation.For water saturated segments, the following residual equations can be found.
Momentum residual for segment N, R M,N : Continuity residual, R C,N , the equations below are for upflow, downflow, convergent flow and divergent flow in segment N.
Energy residual for segment N, R E,N , only shown for the upflow case: The first three terms represent the accumulation of fluid internal energy, halite energy and fluid kinetic energy, the next four terms represent the transport in and out of the segment of enthalpy and kinetic energy, the next term is gravitational potential energy gain or loss depending on velocity direction.
Salt mass residual for segment N, R S,N , only shown for the upflow case: The coupled system of equations is then solved iteratively using the well-known Newton-Raphson method.The procedure is as follow: After initialization Eq. ( 18)), Eqs. ( 19) and ((20) are solved in succession until a stopping criterion is satisfied.
Here, x is the solution vector, x 0 is an initial solution guess (we use the solution from the previous time level), δx is the increment of the solution, J is the Jacobian matrix and R is the residual vector.We assume convergence when all relative residuals R ′ satisfy R ′ < 10 − 4 .
For each segment the relative residual for momentum, fluid mass, fluid + halite energy and salt mass are, respectively: R M,N /p N , The general form of Eq. ( 19) is: With for instance: The location of the air-water interface is calculated at each iteration of the Newton-Raphson loop, it is determined by calculating the total volume of fluid in the well.Knowing the total mass of fluid present in the well at each iteration we can subtract the mass present in each segment (i.e., ρ f LA eff ) starting from the bottom of the well until a segment is only partially filled or we reached the top of the well.The displacement of the interface modifies the well solution (e.g., if modifies the height of the water column and therefore the location of the top pressure boundary condition).Consequently, the displacement also modifies the density used during the mass distribution step.The new location of the interface converges with the rest of the well solution during the Newton-Raphson loop.
When the interface moves from a segment to another, the salt mass fraction in this latter segment is taken equal to the one of the segment of origin.The specific enthalpy is taken equal to the one of the segment of origin and is corrected for the gravitational potential energy change caused by the displacement (the correction is ±g▵h where g is the acceleration of gravity and ▵h is the displacement).
For fully or partially air saturated segments we neglect the weight of air and set the pressure to the wellhead boundary condition pressure.The enthalpy and salt mass fractions of air saturated segments are set to zero.
The well-reservoir coupling is done in a partially implicit manner.From the initial well and reservoir pressures the well-reservoir exchange terms are calculated (Q f , Q h , Q NaCl and Q loss for the well and their opposite for the reservoir), then, updated flow rates and pressures in the well and reservoir are computed iteratively until convergence.We assume convergence when the absolute change in reservoir pressure does not exceed 0.02 bars from one iteration to another.
A flow chart of the global solution method is shown in Fig. 2. The reservoir transport computations are conducted using a control volume finite element scheme developed on the CSMP++ platform (Weis et al., 2014), extended to 3D for this work.The reservoir domain is currently solved using an implicit pressure explicit saturation (IMPES) procedure; with phase-by-phase upwinding.

Verification of the inflow performance relationship
We first tested our implementation of the inflow performance relationship for simple conditions (with constant temperature and fluid properties) by comparing our model to the commercial software Eclipse on a quarter five spots test case.The test case includes a producer well at a pressure of 200 bars and an injector running at 300 bars.The temperature is 140 • C, the constant fluid properties are: 2.10 − 4 Pa.s viscosity, 940 kg/m 3 density and 10 − 9 Pa − 1 compressibility.The reservoir has a thickness of 10 m, a permeability of 1 Darcy (≈ 10 − 12 m 2 ) and a porosity of 30%.The wells are located in the opposite corners of a 1 × km square.The pressure profile between the wells and the well flow rate obtained with our coupled well-reservoir model in CSMP++ matches the one obtained with Eclipse.The well flow rate is 3746.5 ± 0.5 tons per day with Eclipse (with mesh resolutions of 10 and 100 m) and 3745.5 ± 0.5 tons per day with our model (with mesh resolutions of and 100 m).The pressure results are shown in Fig. 3.
To assess the performance in a more realistic case with non-constant fluid properties we ran similar simulations at various mesh resolutions as this should uncover possible influences of rapidly changing fluid properties due to drawdown near the well at supercritical conditions.For this test we chose a temperature of 400 • C, an injector pressure of 320 bars and producer pressure of 230 bar (giving an average pressure of 280 bars).The permeability of the host rock is 10 − 14 m 2 and the porosity 5%.Under those temperature and pressure conditions, water is in supercritical state and is highly compressible (in the order of 10 − 7 Pa − 1 ), the density varies between the injector and producer from 420 to 220 kg.m − 3 .Even though different mesh resolutions will lead to different pressures being computed in the vicinity of the well we found that this difference was not large enough to significantly affect the fluid density and well flow rates.Four mesh resolutions were used, 10, 25, 50 and 100 m, the well flow rates calculated were all 25.6 ± 0.1 tons per day.The results from these simulations are compiled in Fig. S1 in the supplementary materials.Another study confirming the applicability of the Peaceman equation at supercritical conditions can be found in Yapparova et al. (2022b).

Verification of the well model against IDDP-2 well data
We tested the well model by replicating the cold-water injection test that was carried on the IDDP-2 well in January 2017 and during which temperature and pressure data was recorded (Friðleifsson et al., 2017;Saether, 2020).Details about the procedure to replicate both the state of the formation near the IDDP-2 well and the cold-water injection results is published in Lamy-Chappuis et al. (2020b).Fig. 4 shows the good match obtained between the well data and the simulated data.This simulation also demonstrates the need to include an air phase in the well model since the well is not fully saturated with water during the cold-water injection test: the water column in the well only rises to a depth of approximately 800 m.

Model setup
To illustrate the capabilities of the well model, we first conducted a series of simulations in which the superhot resource is formed due to an anomalously high heat flux.While such resources have not yet been drilled, this type of scenario allows focusing on the effect of the Fig. 2. Solution method for the well and reservoir.To improve stability and convergence within the Newton-Raphson loop, we limit the variations in fluid density: the update of the fluid density from one iteration to another is underrelaxed by a factor a = 0.1 (for instance: X n+1 = aX n+1 + (1 − a)X n means that the variable X is under-relaxed by a factor a between the iteration n and + 1).The well model can function in target rate mode whereby the well solution is converged repeatedly while corrections of the wellhead pressure boundary condition are applied to obtain the desired mass flow rate.formation's permeability, temperature and salinity on well starting and short-term performance without the additional complexities resulting from including a magmatic heat source.Simulations with a magmatic heat source are presented afterwards.The two different setups are shown in Fig. 5.
For all simulations, the well is assumed to have a flowing cross section diameter of 10 cm and a constant roughness height of 0.1 mm.We do not account for the friction caused by the holes of the slotted liner or the change in roughness due to halite precipitation.The well is divided into segments of about 100 m length.

High heat flux cases
For the high heat flux with no magmatic intrusion simulations, the model dimensions are 12 × 12 × 8.5 km with the top boundary being placed at 500 m depth, where the formation's water table is located.The choice of the water table location was made to mimic the conditions encountered at the IDDP-2 well.A single well is centered in this volume and extends vertically down to a depth of 5 km.Well-formation exchange occurs in the last 250 m of the well.
We apply top boundary conditions of 50 • C and 1 bar at 500 m depth  This translates to simple linear geothermal gradients ranging from 50 to 100 • C/km, respectively.We also look at the sensitivity to salinity by prescribing either no salinity (pure water) or 3 wt% NaCl, homogenously distributed in the formation at initialization.In places where the fluid phases present in the formation cannot contain 3 wt% of dissolved NaCl under thermodynamic equilibrium, it will be present in the form of solid halite.The formation permeability is either 10 − 14 or 10 − 15 m 2 and decreases log-linearly to 10 − 22 m 2 in the 500-600 • C temperature window to account for the effect of the rock's brittle-ductile transition, resembling the properties of basaltic host rock (Hayba and Ingebritsen, 1997;Scott et al., 2015).A limited set of simulations were conducted with a formation's permeability of 10 − 16 m 2 , only to confirm that this permeability would be insufficient to obtain economically viable well performance.The host rock is otherwise homogeneous with: 5% porosity, 860 J/K heat capacity, 2 W/(m.K) thermal conductivity, 2500 kg/m 3 density and 10 − 10 Pa − 1 compressibility.

Shallow magmatic intrusion cases
For the simulations with a shallow magmatic intrusion the dimensions of the model are (10.8× 10.8 × 7.5 km).The initially emplaced magmatic intrusion is 3 × 3 × 2 km in size, and its roof is located at 3 km depth.In this set of simulations, we try to evaluate the long-term prospects of a geothermal well drilled directly above a magmatic intrusion, the well extends to a depth of 3000 m and the well completion is between 2300 and 3000 m.Specifically, we want to know how a more localized and transient supercritical resource (typically just above a magmatic intrusion) would behave in the long term.
The magmatic intrusion is assumed to be dry and has an initial temperature of 900 • C. We let the system evolve for 15,000 years and perform system saves (snapshots of the simulation from which a secondary simulation can be started) every 1000 years.For these simulations we let the water table be at the surface (the top of the model) and choose the following boundary conditions: 15 • C and 1 bar at the top and a BHF of 150 mW.m − 2 at the bottom (corresponding to "cold" cases presented before).We assume that the reservoir recharge is sufficient to maintain the water table at the surface.For simplicity, we did not include salt in these simulations.The permeability of the formation (k) follows the following law: logk = − 14 − αlog 10 (D + 1), where D is the depth in kilometers and α is set to 1 in simulation A and 2 in simulation B and C. For simulation C we also add a 500 m thick lower permeability caprock at the surface (k = 10 − 15 m 2 ).The additional permeability dependence on temperature is kept as before.

Geothermal well starting and short-term performance for the high heat flux cases
First, we evaluate the feasibility of starting the well using the air pressurization technique.We apply air pressure at the wellhead and let the well equilibrate with the reservoir for the first 5000 min of the simulation (about 3.5 days).The well starting is then achieved by reducing the wellhead air pressure by 1 bar/minute until the water column potentially reaches the surface.As the water column rises in the well, hotter and less dense fluid fills the well, reducing the pressure gradient in the well with respect to the formation and in turn generating a pressure difference between the well and the reservoir.In order for well starting to succeed, sufficient pressure drawdown needs to be generated to overcome the 500 m final difference in water column height between the well and the formation.One limiting factor is that the well will tend to thermally re-equilibrate with the surroundings, therefore cancelling the density difference between the well and the formation at any given depth, and as a result, cancelling the pressure drawdown.A requirement is therefore that the flow rate generated is large enough to suppress conductive heat losses.Another limiting factor for the pressure drawdown is the formation's fluid density gradient, a smaller density gradient meaning that the warmer fluid entering the well needs to travel further for a sufficient density difference to occur.For this reason, there is a minimum initial air pressure required to start the well (i.e., the initial water column needs to be sufficiently deep).To some degree, the rate of air pressure change can also affect the process, although, we did not find it to be a limiting factor.
Fig. 6 shows the minimum air pressure required to start the well depending on the reservoir temperature gradient and for a permeability of 10 − 15 m 2 (the pressure requirement for a permeability of 10 − 14 m 2 would be marginally lower).It shows that super-hot and deep wells, tapping into low density supercritical resources should be particularly easy to start, requiring an initial air pressure of about 20 bars.A relatively colder formation, although already displaying a large temperature gradient of 50 C/km would be very hard to start (requiring air pressure of 150 bars).Fig. 7 shows temperature, pressure and density profiles with various bottom heat flux configurations.It shows the clear distinction in density profile, with a rapid density drop occurring at 4000 m for warmer formations (this corresponds to the onset of supercritical region).The figure also shows more realistic profiles resembling the reservoir state near the IDDP2 well.With this latter configuration, it is apparent that because of the density transition, a well drilled at depths deeper than 3500 m should be easy to start.We then conducted short term (one year long) production simulations, focusing on the 150 and 200 mW.m − 2 BHF cases.To understand the well behavior at various production rates we gradually increase the production rate in steps of 10 kg/s by decreasing the wellhead pressure, until the maximum production rate is achieved (at a wellhead pressure of 1 bar).
During well starting, we target a small production rate of 10 kg/s to mitigate the well blowup.The blowup of the well can occur when the well's fluid density suddenly drops throughout the well leading to a cataclysmic transient outflow.The blowup usually occurs after well starting because at this point a drop in density in the well can no longer be accommodated by a rise in the water column.The well then enters a Fig. 6.Minimum initial air pressure required to start wells depending on the formation geothermal gradient.For each geothermal gradient we run simulations and increase the initial air pressure by steps of 10 bars until the well starting is successful.self-reinforcing loop of pressure drop, density drop and increased inflow rate.
After this transient event, the well inflow rate remains permanently high due to newly established low fluid density and low bottom hole pressure, obtaining a production rate of 10 kg/s can then require very high wellhead pressures (more than 200 bars in the hot and high permeability scenario, see Fig. 8).
In Fig. 8, we compare the effect of the formation permeability on the well performance in the pure water and 200 mW.m − 2 case.In the 10 − 14 and 10 − 15 m 2 permeability cases the fluid flow through the well is almost isenthalpic (gravitational potential energy loss correspond to approximately 50,000 J/kg and heat loss by conduction is of similar magnitude but gradually becomes insignificant at higher flow rates), the produced fluid specific enthalpy is about 3 MJ/kg at the highest mass flow rates.With a formation permeability of 10 − 16 m 2 the well is displaying a meager production rate of less than 2 kg/s.With such a low mass flow rate, the conductive heat loss to the surroundings is more significant, the maximum produced fluid specific enthalpy is lower than  in the other case, at 2.8 MJ/kg.It is interesting to note that the decrease in wellhead pressure required to increase the mass flow rate is accompanied by a large drop in temperature, as the larger pressure drop experienced by the fluid as it rises through the well leads to more expansion cooling.We also observe that, largely due to the drop in density of the produced fluid, an increasingly larger pressure drop is required every time the production is ramped up by 10 kg/s.
Our conclusion is that wells drilled in superhot formations where the permeability is less than 10 − 16 m 2 would not be viable due to the low production rate, this conclusion is probably generally valid.Depending on the rock type and the brittle ductile transition temperature, such low permeabilities could be found at temperatures as low as 350 • C. For our particular configuration, a permeability of 10 − 15 m 2 could be just sufficient, with a thermal power output of 40 MW and a permeability of 10 − 14 m 2 would much more desirable with a thermal energy produced of more than 160 MW.
In the simulation where a thermal energy production of 160 MW and a mass flow rate of 56 kg/s have been reached, the wellhead pressure had to be lowered to one bar.More realistically this well would operate at the previous 50 kg/s step, providing an output of 150 MW for a wellhead pressure of 85 bars, this would limit scaling and provide sufficient pressure for the efficient functioning of the turbines.It is also possible that depending on the configuration of the power plant, an intermediate mass flow rate might offer the best compromise between the enthalpy and the temperature of the produced fluid.
Secondly, we investigated the effect of the formation temperature with simulations where the basal heat flux was 150 mW.m − 2 instead of 200 mW.m − 2 .The results are compiled in Fig. 9.With a colder formation, similar patterns emerge but the fluid produced has a lower temperature and higher density.An interesting consequence is that although the specific enthalpy at the wellhead is almost halved in the colder formation case, the higher mass flow rates achieved are such that the thermal power of the well is increased to 170 MW in the 10 − 14 m 2 permeability case and to 70 MW in the 10 − 15 m 2 permeability case.
Thirdly, we investigate the effect of adding NaCl salt to the produced fluid.Fig. 10 we show the production history results from a 10 − 14 m − 2 permeability and 200 mW.m − 2 basal heat flux simulation where we set the initial formation salt content to 3 wt% NaCl (a similar salinity to the one found a IDDP2).At the depth of production, the salt in mostly present in halite form in the pore space and the produced supercritical fluid has a salinity of 0.005 wt% NaCl.The outcome of adding salt is that the well rapidly becomes totally clogged with halite at the point of production.In the results presented here we arbitrarily set 70% as the maximum halite saturation above which the well-formation exchange will be halted.The true halite saturation value at which this would occur does not matter, as in fact, if the formation is in the vapor + halite field, any pressure drawdown between the formation and the well will be sufficient to provoke a steady and continuous halite buildup, eventually rendering the well inoperable.In this simulation, halite buildup is almost entirely occurring in the well, if we were to use a high resolution mesh near the well it would likely also appear in the formation.Well profiles of several variables of interests, in time intervals of 10 days, are shown in Fig. S2 in the supplementary materials.
Fig. 11 compares the pure water and 3 wt% NaCl case in the colder formation case (150 mW.m − 2 ).A remarkable result is that the addition of salt in the cold formation case does not significantly affect the well performance and in particular does not lead to well blockage, this is because a liquid phase remains in the well during production and carries the salt to the surface.In light of this, it seems safer to target lower temperature regions in saline formations.

Long-term well performance in the shallow magmatic intrusion cases
We ran the different formation permeability scenarios (A, B, and C), with the cooling of a shallow intrusion, until comparable states of the magmatic-hydrothermal system would be reached before the start of the well simulations.For simulation A we choose a starting time for well operations of 2000 years after the emplacement of the magmatic intrusion, for simulation B we choose 4000 years and for C, 5000 years.Fig. 12 shows the initial permeability and temperature on vertical slices of the model for all three cases as well as the temperature distribution at the start of the well simulations.wells are then operated with a constant wellhead pressure of 10 bars for 200 years.Fig. 13 shows the production history in all simulations.
Although the simulations are not directly comparable, they display similar results.There is a steady and steep decline in the temperature of   During the 200 years of production, the total thermal energy produced by the well goes up by 50% as the increasing mass flow rate overcompensate the drop in specific enthalpy.There are decades long oscillations of the well output that correspond to oscillations of the fluid density near the shallowest part of the well completion, this is illustrated in a video that can be found in the supplementary materials.Fig. 14 shows the states of the hydrothermal systems at the end of the well simulations compared to the case where no well operations were performed.It highlights the fact that well production for 200 years only had a very localized effect on the temperature of the hydrothermal system, there is no significant effect of the migration of cold meteoritic water towards the well.

Conclusions
Using a new multi-segment well model we were able to simulate well starting and production from deep and superhot geothermal wells.We found that for formation temperatures of 500 • C near the well  completion, a permeability of 10 − 14 m 2 would be favorable and 10 − 16 m 2 would not be viable.We did not assess the impact of possible local permeability enhancement in such tight formations.In this work we used large completion lengths.The results would be comparable if we were to divide the completion length by a given factor and simultaneously multiply the permeability by the same factor.For instance, a 50 m long completion with a permeability of 5.10 − 14 or a 250 m completion with a permeability of 10 − 14 would lead to similar well inflow performance.Such narrow and permeable completions are the most realistic targets of hot and deep wells; at great depths they would correspond to fractured or faulted zones in an otherwise impermeable formation.
While running well simulations, we did not consider powerplant efficiency or lifetime questions.For instance, fluid chemistry is expected to be a limiting factor for well performance and that a minimum well pressure threshold would be required to avoid the formation of hydrochloric acid due to recondensation of chlorine bearing steam as well as to avoid potential silica and halite scaling issues.We found that the addition of salt in the hottest cases (500 • C) would lead to well blockage due to halite scaling at the point of production.Keeping a well head pressure of 230 bars, therefore limiting production to about 10 kg/s would only delay the moment when the well becomes inoperable.Halite scaling is however not an issue for colder formations (400 • C).
Here we only investigated the water-NaCl system and it is possible that the results could be different if we were to use a more complex fluid chemistry.It is also unclear if the erosion and transport of halite crystals could significantly alter the results, but we do not expect halite crystals to be easily transported in the vapor phase present in the corresponding simulations.
We found that as the fluid ascends the well it decompresses and cools which means that higher mass flow rates obtained at lower wellhead pressures can be associated with a significant drop of the wellhead temperatures.We also found that lower temperature formations can yield a higher thermal energy output since the gain in mass flow rate allowed by higher fluid densities can overcome the loss in specific enthalpy.In short, depending on the power plant design, different optimums of flow rate, temperature and enthalpy would need to be found.Binary cycle power plants using heat exchangers could benefit from the temperature being the highest possible while this might not be the case for dry or flash steam power plants.
We found that wells drilled in hot formations should be easy to start with air pressurization technique even if the formation's water table is located at 500 m depth.If a single producer well were to be drilled, we expect the performance in term of specific enthalpy and temperature to decrease gradually but remain high for several decades.In future work we would like to evaluate how stopping production for a period of time would allow the recharge of the resource.
This study remains generic and more relevant results would stem from the simulation of more specific and realistic cases.Most notably it would be useful to include the whole well history, including cold water injection if relevant.Equally useful would be to represent more accurately the formation permeability structure, including the presence of faults and fractured zones.More complex improvement of the model could be the inclusion of non-condensable gas phases and the upgrade of the homogeneous flow model to a drift flux model.

Declaration of Competing Interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Benoit Lamy-Chappuis reports financial support was provided by Equinor Energy AS. Benoit Lamy-Chappuis reports financial support was provided by European Union's Horizon 2020, GEOPRO.

Fig. 1 .
Fig. 1.Placement of well segment and well variables on 3D host rock mesh.Here, 8 tetrahedral elements belonging to the reservoir mesh are represented.Nodes (black dots) are connected by the edges of the elements.Connecting the barycenters of the elements edges, faces and volume, we construct sector volumes.We represented one sector volume of the front facing and darkened tetrahedron, the one associated with the central node.The assemblage of all sector volumes around a given node forms a control volume.Well segments are placed between neighboring nodes, one such segment, labelled "Segment N" is shown, the entire well is composed of connected segments.We compute wellreservoir exchange terms on the bottom (deepest) node of the relevant segments, where both reservoir and well pressures coexist.

Fig. 3 .
Fig. 3. (a) Plot of fluid pressure on a line connecting the producer and injector wells.The current CSMP++ model results (lines) at resolutions of 100 and 25 m agree well with the high resolution (10 m) Eclipse results (triangles).(b) Fluid pressure contours in the reservoir together with the locations and pressures of the two wells.

Fig. 4 .
Fig. 4. Plot of simulated and well data for a 40l/s injection of 10 • C water in the IDDP-2 well.We assumed that three feedzones at present, at depths of 2300, 3400 and 4400 m.

Fig. 5 .
Fig. 5. General model setup for the high heat flux (a) and shallow intrusion (b) cases.The black and white line in the center represent the well, the white part represents the well completion.

Fig. 7 .
Fig. 7. Depth profiles of temperature, pressure and density with various linear geothermal gradients and in a more realistic case mimicking the highly convective hydrothermal system near IDDP2.

Fig. 10 .
Fig. 10.Production history in the hot and permeable formation case with 3 wt% NaCl salt added to the fluid in the formation.

Fig. 11 .
Fig. 11.Comparison of the well production history in the cold formation case (400 • C, 150 mW.m − 2 ) with or without salt and for different permeability of the formation.

Fig. 12 .
Fig. 12. (a)-(c): permeability structures in simulations A, B and C. (d)-(f): corresponding temperature distributions at initialization and at the start of the well simulations.The well is represented by a straight line in the center, with the part indicating the completion (from 2300 to 3000 m depth) in lighter color.

Fig. 13 .
Fig. 13.200 years well production history in simulations A, B and C.

Fig. 14 .
Fig. 14.Left to right: simulations A, B and C. Top: temperature distribution 200 years after the start of the simulation in the absence of a well, bottom: temperature distribution after 200 years of uninterrupted production.Illustrating the localized effect of the supercritical resource depletion.