A Peaceman-type well model for the 3D Control Volume Finite Element Method and numerical simulations of supercritical geothermal resource utilization

While supercritical geothermal resources receive increasing attention, it has so far remained unclear if and how they could be utilized. In order to provide a tool that can simulate both the natural transient evolution of a supercritical resource near a magma intrusion and its response to possible operation schemes, we augmented the CSMP++ simulation platform with a Peaceman-type well model. For the purpose of generic porous medium simulations of the supercritical resource’s response to direct production, injection or other operation schemes, only a single in-/outflow interval per well is considered and flow in the wellbore is not simulated. The model’s semi-analytical source/sink function accounts for the rapid change in pressure near the wells in the reservoir simulation. As the implementation is based on the 3D Control Volume Finite Element Method it allows the use of unstructured computational grids in order to be able to include geologically realistic geometries. We validate that the well model provides robust solutions of well pressures and rates. Such robust solutions are also obtained for supercritical conditions where water is highly compressible and, therefore, does not strictly fulfil the assumptions used in the derivation of the Peaceman model. To illustrate how a magma-related supercritical resource responds to well operation, we first simulate the evolution of a geothermal system around and above a 2 km deep, explicitly represented magma body. During the system’s hottest phase, ca. 2700 years after its initiation and during progressive inward cooling of the magma body, a well completion interval is activated at 2.1 to 2.2 km depth where temperatures then reach 420 to 490 ◦C. For the relatively low host rock permeability used in the model (10−15 m2), injection and production rates remain quite limited at several kg/s. Yet, the evolving patterns of temperature, enthalpy, density and flow rates provide some first insights into supercritical resource response to operation.


A B S T R A C T
While supercritical geothermal resources receive increasing attention, it has so far remained unclear if and how they could be utilized. In order to provide a tool that can simulate both the natural transient evolution of a supercritical resource near a magma intrusion and its response to possible operation schemes, we augmented the CSMP++ simulation platform with a Peaceman-type well model. For the purpose of generic porous medium simulations of the supercritical resource's response to direct production, injection or other operation schemes, only a single in-/outflow interval per well is considered and flow in the wellbore is not simulated. The model's semi-analytical source/sink function accounts for the rapid change in pressure near the wells in the reservoir simulation. As the implementation is based on the 3D Control Volume Finite Element Method it allows the use of unstructured computational grids in order to be able to include geologically realistic geometries. We validate that the well model provides robust solutions of well pressures and rates. Such robust solutions are also obtained for supercritical conditions where water is highly compressible and, therefore, does not strictly fulfil the assumptions used in the derivation of the Peaceman model. To illustrate how a magma-related supercritical resource responds to well operation, we first simulate the evolution of a geothermal system around and above a 2 km deep, explicitly represented magma body. During the system's hottest phase, ca. 2700 years after its initiation and during progressive inward cooling of the magma body, a well completion interval is activated at 2.1 to 2.2 km depth where temperatures then reach 420 to 490 • C. For the relatively low host rock permeability used in the model (10 −15 m 2 ), injection and production rates remain quite limited at several kg/s. Yet, the evolving patterns of temperature, enthalpy, density and flow rates provide some first insights into supercritical resource response to operation.

Introduction
Magmatic intrusions drive fluid convection in geothermal systems located in volcanically active areas (Stimac et al., 2015;Jolie et al., 2021). Close to the magmatic heat source, fluid temperature and enthalpy may exceed the critical values for pure water, resulting in the development of a supercritical geothermal resource (Scott et al., 2015;Heřmanská et al., 2019). Supercritical geothermal resources have been identified as potentially interesting future targets for geothermal exploration due to their high production potential Reinsch et al., 2017). While 2D numerical modelling has been used to study the influence of host rock permeability, magma emplacement depth, and fluid salinity on the development of supercritical geothermal resources (Scott et al., 2015(Scott et al., , 2016(Scott et al., , 2017 and references therein), the production potential of supercritical geothermal resources above this well suggested a power generation potential of 34 MWe (Axelsson et al., 2014).
Prior modelling studies of production from a supercritical reservoir have used idealized initial and boundary conditions (Yano and Ishido, 1998;Croucher and O'Sullivan, 2008;Magnusdottir and Finsterle, 2015). No prior study has investigated the production potential of a supercritical geothermal resource with an explicitly represented magmatic intrusion. Due to the challenges associated with drilling into high temperatures and producing fluid with high potential for silica scaling (Ingason et al., 2014;Kruszewski and Wittig, 2018), supercritical resources have also been identified as targets for injection rather than production wells, in order to stimulate deep permeability and provide pressure support to the overlying conventional geothermal system and possibly accessing a larger volume of hot rock, extending the lifetime of a geothermal operation.
Due to the fact that the discretization of grid blocks in a numerical reservoir simulation is usually much coarser than the diameter of wells, representing production or injection wells in numerical reservoir simulation requires the integration of a well model (Williamson and Chappelear, 1981). In order to reproduce the steep pressure gradients that develop around wells and to calculate accurate production/injection rates, these models introduce an additional ''well pressure'' variable for each block containing a well in addition to the reservoir fluid pressure. The default well model for reservoir simulation continues to be based on the approach of Peaceman (1978) and is widely used in industry (TOUGH2, Pruess et al., 1999) Wang et al., 2020).
The Complex Systems Modelling Platform (CSMP++; Matthäi et al., 2007) incorporates a numerical scheme based on the Control Volume Finite Element method (CVFEM; Forsyth, 1991;Swaminathan and Voller, 1992). The approach is tailored to ensure stable solutions in the presence of strong gradients in fluid properties and large fluid source terms (Weis et al., 2014), as might arise in scenarios with production/injection wells. In addition, CSMP++ allows for embedded lower dimensional elements, where fractures can be modelled as 2D planes in 3D reservoirs (Patterson et al., 2018) and wells can be represented as 1D lines in 2 and 3D models (Yapparova et al., 2014;Peters et al., 2018). While substantial grid refinement would need to be applied near a well in order to reproduce the steep pressure gradients in the near-vicinity of the well, such grid refinement would result in very small time steps, which in the context of the IMPES-like (implicit pressure, explicit saturation) CVFEM scheme would lead to extremely long simulation times. Moreover, as discussed in Salimzadeh et al. (2019), a well model is necessary even in the presence of significant mesh refinement around a well.
To this end, numerical simulations of supercritical resources have used continuum approaches, representing the resource as a porous medium. Permeability is typically considered as a temperaturedependent property, reflecting its likely response from brittle to ductile rock behaviour (Hayba and Ingebritsen, 1997;Scott et al., 2016) and discrete feed zones such as fractures or faults are usually not represented in the model.
In this paper we propose the Peaceman-type well model implementation for the Control Volume Finite Element Method, following previous work by Fung et al. (1992), Chen et al. (2006) and Chen (2007) on the 2D (Peaceman, 1978) model extensions for the CVFEM. This paper also presents first 3D simulation results obtained using the CVFEM scheme from Weis et al. (2014), which was previously used only in 2D (Weis et al., 2014;Weis, 2015;Scott et al., 2015Scott et al., , 2016Scott et al., , 2017Lamy-Chappuis et al., 2020). In line with the continuum representation of the resource and in order to reduce computational cost, we assume that the well model will be applied to a single interval of in-/outflow such that in-well flow can be ignored. After benchmarking the numerically computed pressure near a well against an analytical solution, we perform a reconnaissance study of generic production/injection scenarios for deep geothermal wells tapping supercritical energy resources.

Governing equations
The governing equations of multi-phase mass and energy conservation are solved using a continuum porous media approach with a pressure-enthalpy-based formulation in a Control Volume Finite Element Method numerical scheme. The implementation of this approach in 2D within CSMP++ has been described in Weis et al. (2014), and is only briefly discussed here. In the 3D extension of the CVFEM pressure unknowns are stored on the nodes of a tetrahedral mesh. In our well model implementation, the wells are represented by a connected set of 1D line elements embedded into the 3D mesh. The flow inside the well is not considered. Each well node represents a well completion in hydraulic connection with a thin reservoir layer perpendicular to the associated well segment. The flow in each such layer is assumed to be radial towards the well (production) or from the well (injection). The point source/sink associated with the well node represents the well completion rate (total rate of the well being the sum of all nodal rates).
Here, we review the governing equations for radial flow into a well, the extension of the Peaceman's model for the 3D CVFEM, and solution of the coupled system in the CSMP++ framework. A brief comparison of our well model and the deliverability model from TOUGH2 (Coats, 1977;Pruess et al., 1999) is presented in Appendix.

Analytical formula for the radial flow near a well
The analytical solution for the radial flow of a single phase fluid near a well in a thin infinitely large reservoir is the basis of the well model derivation (Muskat, 1946;Peaceman, 1978;Chen et al., 2006). This solution is obtained under the assumptions that the reservoir is homogeneous and isotropic (hence permeability is not a tensor but a scalar), and the fluid density and viscosity are constant (Chen et al., 2006). Fig. 1 shows a cylindrical region near a well ( is the well radius), ℎ is the reservoir thickness, is the well mass flow rate.
In order to derive the analytical solution for pressure we need to switch to the polar coordinate system ( , ). The mass conservation equation for steady-state flow without gravity in the Cartesian coordinate system (Chen et al., 2006): where is the volumetric velocity and is the Dirac delta function representing a well as a point source placed at the origin, is transformed to: in the polar coordinate system, where ( , ) = ( )( , ), and an additional condition, that the well mass rate is constant and equal to the mass flux through a cylindrical surface of height ℎ and radius , for any > 0: Darcy's law in polar coordinates in case of the radial flow into a well is: where is the fluid pressure, and hence the mass flow rate can be written as: Separating the variables, and integrating by parts from to , where 0 < < , and from = ( ) to ( ) we obtain: Eq. (8) is the analytical flow model near a well and the well radius is chosen as a reference point for convenience, is subsequently referred to as the ''well pressure'' and -the ''reservoir pressure'' or ''fluid pressure''. As Eq. (8) implies, the steepest pressure gradients develop in the near-vicinity of the well.

Extension of the Peaceman well model for CVFEM in 3D
To avoid the need for extreme mesh refinement near the well, we adapt the well model first introduced by Peaceman (1978) for the Finite Differences discretization (Fig. 2a,b) that includes a concept of an effective radius being a distance at which the numerically computed reservoir pressure value in a block containing a well ( ) is equal to the analytical steady state pressure.
In the Finite Differences discretization, the pressure within a cell is constant, whereas in the Finite Element or Control Volume Finite Element (CVFE) discretization (with linear elements) the pressure varies linearly across the cell (Fig. 2c,d). In case of our CVFEM discretization, the reservoir pressure at the well node is equal to the analytical steady state pressure at the distance of the effective radius.
We closely follow the extension of Peaceman's well model for the Finite Element method by Chen et al. (2006), but additionally extend it to 3D. We begin with a finite element solution of the steady state pressure equation in a 3D domain , where is the well rate evenly distributed across the well nodes. We consider ℎ -a finite element partition (mesh) of using tetrahedra { } and the associated space of piece-wise linear functions ℎ . Unknown pressure is approximated on using the piece-wise linear basis functions : The well is represented as a subset of the mesh containing several nodes (''well nodes'') connected by the edges of the tetrahedra (''well segments''). As we do not model in-well flow this approach should be limited to a single in-/outflow interval per well. In this interval, however, the well can be represented by several segments with the assumption that over the limited extent of the interval in-well pressure follows a near-hydrostatic distribution. For a single well node 0, the point source/sink 0 is the associated well completion rate, 0 is the unknown pressure at this node, and 0 is the basis function ( 0 ∈ ℎ ) with support 0 , where 0 ⊂ (Fig. 3).
Eq. (9) can be written in a weak form on 0 (Chen et al., 2006): where is the number of neighbouring nodes of the well node 0 (connected by an edge) and is the transmissibility between the two nodes, where is the number of tetrahedral elements that contain the edge between those nodes, | | is the volume of a tetrahedron. The main difference between our derivation for 3D and the derivation for 2D from Chen et al. (2006) is that the layer thickness ℎ does not appear in Eqs. (10)-(13), instead the pressure equation is discretized in a volumetric domain 0 .
We assume that the pressure in all nodes sharing an edge with the well node can be obtained exactly from the analytical solution (Eq. (8)): where is the edge length between the two nodes. The layer thickness ℎ is the sum of the halves of the well segment lengths attached to the well node. As the analytical solution was derived under the assumptions of the radial flow in a thin layer, naturally a restrictive condition is imposed on the mesh: lateral element edges have to be bigger than the layer thickness ( ≫ ℎ, ∀ ). Substituting Eq. (14) into Eq. (12) and rearranging for 0 : where is the analogue of the equivalent radius from (Peaceman, 1978): Eq. (17) relates the ''well pressure'' and ''reservoir pressure'' of the well node 0. The skin factor takes into account either the formation damage due to drilling (positive skin) or flow enhancement due to hydraulic fracturing or acid treatment of the well (negative skin), and can be incorporated into Eq. (17), rearranged for the well node rate: The correction for gravity is not needed in Eq. (19) as the ''well pressure'' and ''reservoir pressure'' always share the exact same location on the well node, unlike in the Integrated Finite Differences discretization where the actual location of the well completion does not always coincide with the center of a cell.
For a well comprising nodes the total well rate is: A. Yapparova et al.

Fig. 2.
The concept of the Peaceman's well model in 2D: for the finite differences discretization (a) a reservoir block containing a well and its neighbouring block, the reservoir pressure and the well pressure associated with a block containing a well, is the reservoir pressure of the neighbouring block; (b) analytical pressure distribution (dashed line) and numerically computed reservoir pressure (solid lines) in the blocks from (a), the equivalent radius is a distance at which is equal to the analytical pressure; for the finite elements discretization (c) both reservoir and well pressures are associated with a well node and is a reservoir pressure of a neighbouring node, that is connected by an edge with a well node; (d) analytical pressure distribution (dashed line) and numerically computed reservoir pressure (solid lines) for a well node and its neighbouring node from (c), the equivalent radius is a distance at which is equal to the analytical pressure.

Fig. 3.
Tetrahedral mesh around the well node 0, 0 is the pressure at this node, is the pressure at a neighbouring node, connected by an edge ; 0 contains all tetrahedra that include node 0; ℎ is the sum of half lengths of the edges connecting the well node 0 to two neighbouring well nodes.
where is the well node rate extended for two-phase flow under the assumption of identical phase pressures: where , , with = , for liquid and vapour, are the phase relative permeabilities, densities and viscosities, respectively, at well node ; permeability of the Control Volume centred at the node is calculated as a harmonic average of the permeabilities of all elements contributing to this volume; ℎ is the sum of the half lengths of the well segments attached to the node; the effective radius is calculated for each node separately.

Full system of CVFEM equations for multi-phase flow in porous media
The Control Volume Finite Element Method (CVFEM) implemented within the CSMP++ software library allows modelling of the multiphase flow of liquid water and vapour in porous media. Phase velocities are obtained using an extended two-phase form of Darcy's law: where is the gravitational acceleration vector. A linear relative permeability model with a liquid residual saturation of 0.3 and vapour residual saturation of zero is adopted (Hayba and Ingebritsen, 1997).
The governing equations for fluid (Eq. (23)), salt (Eq. (24)) and heat transport (Eq. (25)) are listed below, and the solution procedure is described in detail in (Weis et al., 2014): The conserved properties are: the total fluid mass per unit pore volume, where is the phase volumetric saturation, subscript ℎ stands for halite, is the porosity; the total mass of salt per unit volume, where is the salt mass fraction in the th phase; and the total enthalpy of a system, split into the fluid-and rock-dominated parts: where = ℎ + ℎ is the total fluid enthalpy and ℎ , = , are specific enthalpies of liquid and vapour phases, is the rock density, -the rock isobaric heat capacity, -the temperature, and -the thermal conductivity; H 2 O+NaCl , NaCl and ℎ are the source/sink terms for fluid, salt and heat, respectively. The addition of the wells is taken into account by the well rate term (Eq. (20)) entering Eqs. (23) and (29).
The transient pressure Eq. (29) is derived from Eq. (23): where is the total system compressibility, which we approximate as = + (1 − ) , with and being the fluid and rock compressibilities, respectively; and the term − | | | accounts for the fluid density change at constant pressure due to temperature and composition variations.
The system of equations is completed by the Driesner (2007) and Driesner and Heinrich (2007) EoS for the H 2 O − NaCl fluid system, which is valid in the temperature and pressure range from 0 to 1000 • C and 0 to 5000 bar.

Well control and solution of the coupled system
Wells are commonly simulated under the assumption of a given rate of fluid production/injection or a given bottom hole pressure (bhp) (Lie, 2019). In the case of the fixed rate control, injection/production rate directly enters the pressure Eq. (29) in terms of the mass source/sink . For the numerical solution of Eq. (29) the overall well rate is evenly distributed across all well nodes. The reservoir pressure is obtained using the CVFEM. The well pressure is calculated in a post-processing step (rearranging Eq. (21)): on a node basis for all well nodes. The hydrostatic pressure difference across the well nodes arises automatically as in the reservoir pressure solution the gravity terms are taken into account. For the bottom hole pressure control, the well pressure is known for the lowermost well node ( 1 = ℎ ) and the well pressures at the nodes above are calculated according to the hydrostatic pressure distribution, with density taken average between the two nodes: where the well nodes are numbered consequently from the bottom of the well towards the well-head, is the fluid density at the well node , is the length of the line segment between the nodes and −1. The unknown reservoir pressure additionally appears in the source term (Eq. (20)) of the pressure Eq. (29) for all well nodes. Explicit treatment of the well source term, where fluid properties and reservoir pressure are taken from the previous time step produces oscillations in the pressure solution. Therefore, a semi-implicit treatment of the source term was employed, where fluid properties are taken from the previous time step and the reservoir pressure at the new time step. After the reservoir pressure is obtained numerically, the well rate can be first calculated using Eq. (21) for each node, and then summed up for the total injection/production rate of a well (Eq. (20)); individual phase rates can be obtained trivially by using phase mobilities.
In the numerical simulations presented below, both reservoir and injected fluid were assumed to be pure water. However, in general the CVFEM is capable of modelling the multi-phase flow of H 2 O − NaCl brine including halite precipitation (Weis et al., 2014), and no additional modifications are necessary to perform such simulations with wells. In this case, the salinity of the produced fluid (including the distribution of NaCl between liquid water and vapour) will be calculated in the CVFEM and for injection wells, salinity of the incoming fluid will need to be specified as a boundary condition.

Quarter of the five-spot benchmark
The assumptions considered above for the development of the line well model are quite strict -constant rock and fluid properties. In reality, geothermal reservoirs are usually heterogeneous and have a complex permeability structure and varying fluid properties. Analytical solutions describing the transient fluid flow in such systems typically do not exist. This produces a challenge for the well model validation, as an idealized model setup has to be used.
One popular benchmarking case from the petroleum reservoir modelling tradition is the so-called quarter of the five-spot model. Oil & gas reservoirs are often produced with wells drilled in a regular wellpattern (e.g five-spot, seven-spot, nine-spot). In the five-spot pattern one production well is surrounded by four injection wells. Due to its repetitiveness and symmetry along both vertical axes, and the minimum number of wells in a single pattern, the five-spot pattern is often chosen for benchmarking, as only a quarter of the well-pattern with one injection and one production well needs to be computed.
In the following, we first perform a quarter of the five-spot model benchmarking study with idealized (constant) fluid properties to validate our technical implementation. This is followed by a mesh convergence study with real fluid properties near 400 • C at pressures where strong fluid property variations occur. We demonstrate that mesh convergence is achieved and that the well model should therefore deliver robust result even under supercritical conditions.

Model setup
The quarter of the five-spot model has been constructed in CSMP++ for the purpose of benchmarking the line well model (Fig. 4). The model comprises a 1000 m square domain with a thickness of 30 m. Two wells of 0.1 m radius are located in the opposite corners, an injection well with the bottom hole pressure (bhp) set to 300 bar and a production well with 200 bar bhp. The temperature follows a small vertical linear gradient between 140 and 140.75 • C across the model. The porous medium has a constant porosity = 0.3, permeability = 10 −12 m 2 , and rock compressibility = 5 ⋅ 10 −5 1/bar. The resulting flow is at single phase liquid conditions. The model domain is divided into three 10 m thick horizontal layers and the wells consist of three 10 m long line segments in the opposite model corners. Each layer is meshed using tetrahedral elements with a side length of 100, 50 and 25 m, producing three progressively finer meshes for comparison.

Explicit solution
The explicit solution by Muskat and Wyckoff (Muskat, 1946;Lie, 2019) for the pressure difference between the injection and production wells in the repetitive five-spot well placement is given by: where is the pressure difference between the injection and production wells, is the distance between the wells. Eq. (32) allows to obtain the value of the injection/production rate for the benchmark setup. Note that Eq. (32) assumes constant fluid and rock properties throughout the model.
We have used the density and viscosity values at the node in the middle of the model (15 m depth, = 251.38 bar, = 140.375 • C) for the   (32) we calculate the injection/production rate of = 163.2 kg∕s for each well. As we have only a quarter of the repetitive five-spot well pattern, in a numerical simulation each well has only one quarter of the analytical ratê= 40.8 kg∕s. To take this into account, a new factor, -the angle multiplier, has been added to Eqs. (17) and (18): = exp where = 0.25 for a well located in the model corner, and̂= . Analogously, = 0.5 can be used for a well located at the model boundary.

Numerical results
Using the calculated well production/injection rate obtained from Eq. (32), we can use Eq. (8) to retrieve the analytical solution for the reservoir pressure near the well. If constant pressure conditions are applied directly to the well nodes without using the well model, the numerical solution is mesh-dependent and fails to reproduce the rapid pressure drop in the near-vicinity of the well (Fig. 5). This highlights the need for a well model.
We conducted a total of six benchmark simulations with the well model, applying either constant bottom hole pressure or a constant rate for both wells on three different meshes. In either case, the numerical solution for reservoir and well pressures agrees well with the analytical solution (Fig. 6), reproducing the logarithmic drop in pressure in the near-vicinity of the well. Numerically-computed pressure curves do not show a strong dependence on the mesh resolution, particularly for the injection well (Fig. 6a,b). However, small discrepancies between the numerical and analytical solutions would be expected as the analytical solution formula assumes constant fluid properties, which is not the case in our simulations. Fluid properties vary more strongly around the production well than the injection well, as injection temperature is fixed at 140 • C. As a result the simulations with a fixed rate control (Fig. 6b,d) show computed well pressure values very close to the designed 300 and 200 bars (Table 2), and the calculated well rates for the simulations with a fixed bhp differ only slightly from the designed 40.8 kg/s These numerically-obtained results are consistent with the analytical solution regardless of whether a fixed bhp or fixed rate approach is used. Table 3 presents the effective radius values computed for well nodes at different depths for the three different meshes. The average value for the reservoir pressure at the production well is ∼226 bar, this pressure value for the analytical pressure solution is achieved at a distance of ∼10 m, as can be seen from Fig. 5a. This correspondence between the computed reservoir pressure and the analytical solution at the distance from the well given by the effective radius shows that the numerical calculation of the effective radius is reasonable.
In Section 2.2 we discussed a restrictive condition on the mesh ≫ ℎ, that arises from the assumption of a radial flow in a thin and horizontally extended reservoir necessary for the derivation of the analytical solution for fluid pressure near a single well. However, our numerical results show that for lateral mesh size to well segment length ratios ∶ ℎ of 100 ∶ 10, 50 ∶ 10 and even 25 ∶ 10 the well model is stable. At the same time we do not impose any regularity conditions on the mesh surrounding the well (as for example is described in Salinas et al. (2021), where they enclose the well in a triangular elements sleeve).

Mesh convergence study under supercritical conditions
We have also performed a mesh convergence study at higher temperature (400 • C) and pressure conditions (initial pressure 280 bar, injection well bhp = 330 bar, production well bhp = 230 bar), where the analytical solution is not applicable (Fig. 7). The fluid pressure in the reservoir stays almost constant and equal to the initial pressure value everywhere except the regions near both wells; the finer the mesh, the better the rapid pressure change near the well can be resolved with the CVFEM.
The well model shows accurate results even at supercritical conditions, the mesh dependence of the well model itself is minor. The computed reservoir pressure values at the well are almost the same for all three mesh resolutions (Table 4); the well rates are close to each other for 100 m and 50 m meshes and higher for the 25 m mesh (Table 4), supposedly the lateral size to segment length ratio is not sufficiently big in case of the 25 m mesh.
The good performance of the well model for supercritical conditions with significantly varying fluid properties may seem surprising given that the Peaceman formula was originally derived for constant fluid properties. We interpret the result of the mesh convergence study such that the spatial and time discretization of our approach approximates this assumption well enough for the individual computational cells. This is in line with the fact that Peaceman-type well models have successfully been used for a long time in simulators such as TOUGH2 to simulate flow of twophase liquid+steam mixtures around wells although such mixtures can have very high, saturation-and temperature-dependent compressibility (Grant and Sorey, 1979). is the average finite element size in horizontal direction, vertical mesh size (layer thickness) ℎ is 10 m in all cases.

Utilization scenarios for geothermal wells tapping supercritical resources
In this section we apply the well model in CVFEM to modelling of injection into and production from a supercritical geothermal resource. These represent a first reconnaissance application to illustrate well/supercritical resource interaction, while a comprehensive survey of parameters space is beyond the scope of this contribution.

Model setup
Our model setup was partially inspired by the conditions encountered in the IDDP-1 well in the Krafla geothermal system in Iceland (Axelsson et al., 2014;Scott et al., 2015). The 3D modelling domain is 5 km deep and 10 km 2 in vertical and horizontal dimensions, respectively (Fig. 8). An initial conductive heat flux of 75 mW/m 2 results in an initial geothermal gradient of 37.5 • C km −1 . Initial pressure distribution is hydrostatic and all boundaries are open to flow, with pressure and temperature fixed at ambient conditions (1 bar, 10 • C) at the top boundary.
A sill-shape magmatic intrusion, 1 km thick and of a horizontal diameter of 3 km, is instantaneously emplaced into the host rock with a top depth of 2 km depth. The intrusion is initially at 900 • C, but cools via conductive heat transfer across the brittle-ductile transition (BDT) to the convecting hydrothermal system in the surrounding host rock. We simulate the full evolution of the geothermal system from the incipient stage soon after the initial magma emplacement, passing through the main stage, when boiling connects from the surface to depth, through the waning stage, when boiling remains at the surface but the intrusion is cooled at depth (Scott et al., 2016).
The permeability of the host rock is initially set to 10 −15 m 2 , as this allows relatively large supercritical resources to develop (Scott et al., 2015) but, on the other hand, limits injection and production rates. The effect of more conductive feed zones has not been included in A. Yapparova et al. Fig. 6. Benchmarking results: numerically computed reservoir (lines) and well pressures (isolated markers) for three different meshes compared to the analytical solution; is the average finite element size in horizontal direction, vertical mesh size (layer thickness) ℎ is 10 m in all cases. these reconnaissance models. To mimic the reduction in differential stress at high temperatures and the closure of brittle fractures across the BDT, permeability is temperature-dependent (Hayba and Ingebritsen, 1997;Scott et al., 2016). As is appropriate for basaltic host rocks, the onset of the BDT occurs at 550 • C (Violay et al., 2012), above which temperature the logarithm of the permeability decreases linearly to a minimum of 10 −22 m 2 at 700 • C.

Natural geothermal system evolution
The explicit representation of the magmatic heat source leads to a transient natural evolution of the geothermal system without a steady state (Scott et al., 2016). Immediately after emplacement, the magma chamber begins cooling by conductive heat transfer to the surrounding host rock, heating the groundwater and initiating hydrothermal convection (Fig. 9a). At early times (1.5 ka after emplacement, Fig. 9a), only groundwater in the near-vicinity of the intrusion has been heated by the intrusion. The thermal structure of the developing geothermal system consists of a central down-flow region of cold meteoric water from the surface, surrounded by a uniform pipe-like upflow partially at two-phase boiling conditions. As the radius of the central downflow gradually decreases, the upflow pipe merges into a single upflow zone and boiling reaches the surface after 2.1 ka (Fig. 9b).
The geothermal system reaches its main stage (Figs. 9c and 10), by 2.7 ka after the emplacement of the intrusion the magma chamber has substantially cooled down, especially on the flanks (only the middle 1 km remains ductile). The upper ∼300 m of the original magma chamber have become permeable (Fig. 10d,h) and a potentially exploitable supercritical resource ∼300 m in radius develops here. In this region,  temperature is above 400 • C (Fig. 10e), fluid enthalpy is higher than 2.8 MJ/kg (Fig. 10f) and the fluid is at single phase vapour-like conditions (Fig. 10g).
In the following subsections, a vertical well completion interval (5 well nodes with 25 m spacing) is modelled in the middle of the superhot area between 2.1 and 2.2 km, and injection and production scenarios are simulated.
While this superhot single-phase area is sustained between 2.5 and 3 ka, by 3.2 ka the geothermal system enters its waning stage (Fig. 9d), as the intrusion has completely cooled down and a boiling extending from the surface to 1 km depth is underlain by the single phase liquid.

Injection scenario
At the well completion interval between 2.1 and 2.2 km depth, 2.7 ka into the geothermal system evolution (Fig. 11a,b), temperature is in the range of 420 to 490 • C and the reservoir pressure is ∼19 MPa, similar to the conditions encountered at IDDP-1 at Krafla (Axelsson et al., 2014;Scott et al., 2015). Direct production of fluid at such high temperatures is challenging to the production equipment and well casing integrity. Cold water injection could expand the geothermal resource, cool down the area near the well to temperatures more favourable to the equipment by mixing of injected and supercritical fluid, and increase fluid pressures in the overlying conventional geothermal resource.
To investigate the general response of the supercritical zone to injection we have modelled injection of 80 • C water into the supercritical zone at a rate of 3.2 kg/s. Fig. 11 shows the resulting distribution of temperature (c,e) and liquid saturation (d,f) around the injection well after 100 years of injection. A zone of single-phase liquid ∼180 m in radius forms around the injection well with a sharp boundary to the hot single phase supercritical fluid. This suggests limited mixing of the injected fluid and supercritical fluid, and contraction rather than A. Yapparova et al. expansion of the supercritical resource. Fluid pressure near the well builds up in response to injection (Fig. 12); our well model allows to capture the sharp pressure increase near the injection well, where the ''well pressure'' is ∼30 MPa and the ''reservoir pressure'' is ∼23 MPa.
Cold water injection increases the rate of cooling of the magma chamber (Fig. 13), and the depth of the brittle-ductile transition retreats ∼30 m further into the intrusion compared to what would be expected if the system evolved naturally. Given the small rates used in the model, the injected fluid was not able to significantly impact the reservoir pressure in the shallow part of the geothermal system.

Production scenario
We simulate production from the supercritical zone 2.9 ka after the magma emplacement, when the temperature along the vertical line from 2.1 to 2.2 km depth through the centre of the potential supercritical resource is 396-446 • C (Fig. 14a,c). One hundred years of fluid production at a constant bhp of 100 bar leads to a decrease in the temperature at the bottom of the well from 446 to 353 • C and the whole region surrounding the well falls into the two-phase liquid-vapour conditions (Fig. 14b,d). In comparison, during the natural system evolution the temperature at the bottom of the well within 100 years decreases from 446 to 419 • C only.
Remarkably, fluid production affects the liquid saturation in the whole upflow boiling zone, especially in the ∼500 m above the well, where the liquid saturation is increased from ∼0.7 to ∼0.9, compared to the liquid saturation of ∼0.8 after 100 years of natural system evolution.
Fluid production at a constant bottom hole pressure shows a monotonically increasing well mass rate (from 2.5 to 5.4 kg/s; Fig. 15c), due to the increase in production fluid density with decreasing temperature. Lowering the bottom hole pressure from 100 to 60 bar naturally results in a higher well rate (from 3.5 to 7.8 kg/s) at a cost of a more rapid decline in the temperature of the produced fluid (Fig. 15a). Two-phase conditions are encountered at the lower bhp well ∼20 years prior to the higher bhp well, which is reflected in the drop of enthalpy of the produced fluid (Fig. 15b).
The well thermal energy rate ( ℎ = ⋅ ℎ; calculated on the nodeper-node basis and summed up for total) is also increasing with time ( Fig. 15d), with growth in the mass rate overpowering the decrease in enthalpy of the produced fluid. Oscillations are attributed to the finite well discretization, with each drop in the energy rate associated with a single well node first encountering two-phase vapour-liquid conditions.

Discussion
In our simulations simple well operation protocols have been applied, an injection well was modelled under a fixed rate condition and a production well under a fixed bottom hole pressure control. Fixed pressure control appears to be more intuitive for fluid production, than the fixed mass rate, as the density of produced fluid changes in time. While operating a large geothermal field more sophisticated conditions could be applied for well control, e.g. constant thermal output, well head pressure control, pressure/rate changing with time dependent on the well output, etc. However, the fixed rate and fixed bhp used in this model can be considered end member cases of well operation.

Comparison of fixed rate and bottom hole pressure controls
We have performed two simulations complementary to the simulations described in sections above. For the injection well simulation with a fixed rate of 3.2 kg/s we calculated the corresponding bottom hole pressure at initial conditions and started a simulation with a fixed bhp of 300 bar (Fig. 16d). For the production well simulation with a fixed bhp of 100 bar we obtained the corresponding well rate and started a simulation with a fixed rate of 2.5 kg/s (Fig. 16b). In this way, two pairs of simulations -for a single injection and for a single production well start from the same initial conditions but have different well controlfixed rate or fixed bhp.  The temperature of the injected fluid is fixed at 80 • C and its density is very close to 1 kg/l, therefore the volumetric well rate is approximately equal to the mass rate and the conditions near the injection well are close to the constant fluid properties assumptions of the Peaceman's well model. In a simulation with a fixed injection rate the bottom hole pressure is increasing with time (Fig. 16a) and in a simulation with a fixed bottom hole pressure the injection well rate is decreasing with time (Fig. 16d); this behaviour agrees with an intuitive understanding of the process.
In contrast, the fluid properties near a production well experience strong variation. For the production well simulation with a fixed bottom hole pressure the mass rate increase with time is counterintuitive (Fig. 16b); in a conventional oil & gas, groundwater or low-T geothermal reservoir one would expect the well rate decreasing with time due to the pressure drawdown. However, in our case although the volumetric well rate is decreasing with time as expected (Fig. 16c), the mass rate more than doubles within 100 years, due to the increasing density of the produced fluid. The bottom hole pressure increases by 40 bar in 100 years for a fixed rate production simulation (Fig. 16e), as A. Yapparova et al. Fig. 11. Temperature (a) and liquid saturation (b) prior to injection, and after 100 years of cold water injection with rate = 3.2 kg/s in the vertical slice along the well (c,d) and the horizontal slice (e,f) through the well middle (2150 m depth). Well completion location is shown in white. a progressively smaller pressure difference (between the well pressure and the reservoir pressure) is needed to extract less and less volume of fluid. Even though the mass rate curves diverge significantly for the two simulations, the volumetric rate evolution curves follow the same trend. The density growth with time dominates over changes in the mass rate ( = ∕ ).
The energy rate ( ℎ = ⋅ℎ) is dominated by the growing mass rate for the simulation with the fixed bhp and is increasing despite decreasing enthalpy of the produced fluid (Fig. 16f). For a constant mass rate the energy rate is decreasing together with enthalpy. These two curves suggest that a production strategy of keeping a constant energy output of the well -slightly increasing the bottom hole pressure to keep the mass rate increasing just enough to overcome the effect of decreasing enthalpy of the produced fluid -possibly lies in between the two well control options-fixed bhp or constant mass rate. Dumkwu et al. (2012) provide an extensive review on the well models used in numerical reservoir simulation, discussing both the conventional Peaceman-type well models and more advanced wellbore models. They give an overview of the Peaceman-type well model extensions that allow to account the effects of (1) complex well geometry (horizontal, deviated, multilateral wells, wells located not in the center of the grid block, multiple wells in one grid block) and (2) permeability structure near a well (heterogeneous permeability, partial penetration of wells, skin effects).

Limitations of the Peaceman-type well model
In the context of the Control Volume Finite Element Method (CVFEM) within CSMP++ the simulations are performed on unstructured tetrahedral grids and the wells are explicitly represented as lower dimensional features of the mesh. This allows for high geometrical flexibility of well trajectories: vertical, horizontal, deviated, multilateral wells can be easily discretized as sets of well nodes connected by 1D line segments, embedded into the tetrahedral mesh. No modifications to the Peaceman-type well model presented in this paper are necessary to represent the flow near non-vertical wells. Regarding the second type of the well model limitations discussed in Dumkwu et al. (2012), heterogeneous rock permeability and skin-effects can be taken into account in our well model, but the simulations presented in this paper have focused on the homogeneous permeability near the wells.
Assumptions of the constant fluid properties made in the Peaceman well model derivation are violated even at low temperature conditions. For example, in our 140 • C benchmark study the fluid viscosity and density variations were negligible (0.6% and 0.3%, respectively), but the numerical results still deviated slightly from the analytical solution. Nevertheless, the Peaceman-type well models have been used successfully in numerical reservoir simulation over that past 50+ years, essentially due to the lack of an alternative model.
We have conducted an analysis of fluid properties variation at high temperature conditions in the near-vicinity of injection and production  wells in the simulation results presented above (Sections 4.3 and 4.4, 100 bar case). In case of an injection well the density of the injected fluid is constant, whereas the density of produced fluid varies with time. In the beginning of injection the sharp contrast in temperature between the colder injected fluid and the original reservoir fluid results in high density variations near the well (Fig. 17, left, results after 1 year). However, our results show that if the injection is sustained for a longer time, the reservoir fluid in the immediate vicinity of the well is progressively cooled down and the density differences near the well are reduced (Fig. 17, left). After 10 years of cold water injection, relative density variation within the 30 m radius of the well does not exceed 5% (Table 5). Fluid viscosity variation within the 30 m distance from the injection well also level with time, but remain high even after 10 years of cold water injection, due to viscosity dependence on pressure. In case of a production well, the density variation in the 30 m vicinity of the well remains close to 15% during the first 10 years of production; the viscosity of the produced supercritical fluid is infinitesimally low and has a limited variation near the well (Fig. 17, right). As discussed in Dumkwu et al. (2012) the conventional well models are ''used to provide the source and sink functions for numerical simulators'' and are ''not comprehensive enough to adequately represent the physics of the processes taking place in these wells''. The Peacemantype well model or a simple ''deliverability model'' from TOUGH2 do not model flow in the wellbore and do not take into account pressure losses in the wellbore due to friction, heat exchange between the wellbore and the reservoir, slip between the phases etc., unlike the more sophisticated wellbore and coupled wellbore/reservoir models (Livescu et   2020). Nevertheless, we believe that the simulations using a simple Peaceman-type well model together with an explicit magmatic heat source representation and a transient geothermal system evolution can provide insight into the possible response of a supercritical resource to direct production and cold fluid injection. Furthermore, the Peacemantype well model will serve as a basis for the development of a coupled wellbore/reservoir simulator.
Although we have benchmarked our well model against the analytical solution, performed a mesh convergence study and estimated the magnitude of fluid properties variation near the well, the only reliable way to assess its accuracy would be to compare our modelling results with the real production data from a supercritical well. Another possibility would be to compare simulations that use the Peaceman model to simulations on meshes highly-refined around a line representing a well (without the well model), including extreme cases such as cold injection into a supercritical reservoir. However, we consider such a computationally expensive study beyond the scope of the present work.

Outlook
In the simulations presented in this paper temperature-dependent permeability has been used in the region of the original magma chamber. The initial host rock permeability has been assumed homogeneous and depth-independent. Natural geothermal systems often comprise layers of contrasting permeability, cross-cut by faults and fractures, and the wells typically have an open wellbore or a slotted liner completion below ∼1 km and are hydraulically connected to the reservoir via several discreet feed zones. In contrast, in our simulations both wells have a single and rather short production interval of constant permeability. This could be a reason for rather low well rates that are observed in our simulations.
New deep geothermal wells tapping supercritical resources to be drilled in Icelandic geothermal fields would be surrounded by many existing conventional wells with a long history of production and reinjection. The pressure distribution and fluid flow near such wells will  be therefore very different from our idealistic cases of isolated wells. Future studies will include different well configurations.

Conclusions
We have implemented a semi-analytical Peaceman-type well model in 3D Control Volume Finite Element method (CVFEM) and benchmarked it against the analytical solution. Mesh convergence studies at supercritical conditions have exhibited numerically stable behaviour of the well model. This paper presented first reconnaissance utilization scenarios' simulations for geothermal wells tapping the potential supercritical resources. Our well model coupled to the CVFEM scheme allowed us to model injection of cold water into a region with initial temperatures exceeding 400 • C and supercritical fluid production from that area.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. and the Peaceman-type well model implemented in the Control Volume Finite Element Method (CVFEM) within CSMP++ (this paper): 1. The main idea of the well (deliverability) model is based on an analytical solution for pressure near a single well. Assuming single phase steady state radial flow towards a well (with a radius ) in a thin horizontal layer (of thickness ℎ) with constant fluid (density , viscosity ) and rock properties (permeability ), the well rate is proportional to the difference between the ''well pressure '' and ''reservoir pressure'' (eq. 12 in Coats (1977); eq. 10 in Peaceman (1978); eq. 28 in Pruess et al. (1999);Eq. (19), this paper): ( − ).
The major difference between the models is in the formula for the effective (equivalent) radius ; in Coats (1977) and in (Pruess et al., 1999): where is the grid block area, = for an areal Cartesian grid; (Peaceman, 1978) provides a different formula for Cartesian grids: The difference in two formulas stems from the fact that Coats (1977) starts his formula derivation by assuming that the grid block pressure is equal to the average grid block pressure calculated in the simulator, but Peaceman (1978) argues strongly that it is not the case. Peaceman (1978) proves both analytically and numerically that the grid block pressure is not the average pressure in the grid block, but the actual flowing pressure at the distance of the effective (equivalent) radius; explicitly showing that Coats (1977) formula for the effective radius is not accurate (eq. 12 and fig. 3 in Peaceman (1978)). Let us compare the two formulas for a square Cartesian grid.
In (Coats, 1977): In the CVFEM of CSMP++, 3D unstructured tetrahedral grids are used, and the Peaceman's formula for the effective radius on the Cartesian structured grids is not applicable. Chen et al. (2006) extended the Peaceman well model for the CVFEM in 2D and suggested a new effective radius formula, which we have modified for the 3D CVFEM (Eq. (18), this paper). 2. In Coats (1977), Pruess et al. (1999) and CSMP++ individual phase rates at two-phase boiling conditions are obtained using phase mobilities (eq. 14 in Coats, 1977;eq. 28 in Pruess et al., 1999;Eq. (21), this paper): where subscript = , is either liquid or vapour. 3. For wells spanning several model layers (nodes) = 1, and operated with a prescribed bottom hole pressure , the well pressures in the layers (nodes) above the lowermost layer are calculated in a recursive manner both in TOUGH2 (eq. 32 in Pruess et al., 1999) and CSMP++ (Eq. (31), this paper): 4. Wells operated under a fixed rate condition are treated similarly in TOUGH2 and CSMP++, the total well rate is divided between the well cells (nodes). The way the well is operated under a fixed rate condition (target rate) in Coats (1977) is more sophisticated than what is implemented in TOUGH2 and CSMP++. In Coats (1977) the target rate is not distributed between the layers, but one value is prescribed for the whole well. The individual layer rates become additional unknowns. A similar approach is implemented by Coats (1977) for wells operated under a fixed bottom hole pressure. The well pressures in the individual layers are not assigned using a simple hydrostatic pressure gradient like in TOUGH2 and CSMP++ but become additional unknowns. As a result, the numerical solution described in Coats (1977) has a high computational cost, and maybe this was one of the reasons why a simplified version was implemented in TOUGH2.
The brief analysis presented above allows us to conclude that the well deliverability model from Coats (1977) is very similar to the Peaceman (1978) well model, except in the formula for the effective (equivalent) radius calculation.