Geothermal wells drilled into a supercritical geothermal resource, where temperature exceeds the critical temperature of water (~374 ºC) could generate an order of magnitude higher power output than the majority of conventional geothermal wells. However, the best strategy for supercritical resource utilization - direct production of single phase fluid, or injection of cold water to expand the potential resource, or a combination of both - has remained unclear. Due to the strong, non-linear temperature and pressure dependence of fluid properties at supercritical conditions as well as the strongly varying rock properties near a magmatic heat source, numerical simulation of such resource utilization has remained challenging as these characteristics preclude common simplifications typically used in geothermal reservoir modelling.

In this work we address these challenges and present reconnaissance 3D numerical simulations of a geothermal system evolution from the time of magma emplacement and through a subsequent formation of a supercritical geothermal resource. The simulations of two-phase liquid-vapor flow with boiling are performed using the 3D extension of the Control Volume Finite Element Method (CVFEM) [1] within the CSMP++ software library [2] together with an accurate equation of state for H2O-NaCl valid in a broad temperature (0 to 1000 ºC) and pressure (1 to 5000 bar) range [3]. The CVFEM scheme is locally mass-conservative and is able to capture strong gradients in fluid properties arising at the contact between a magmatic intrusion and the host rock. The magmatic heat source is explicitly represented as a region of an unstructured mesh and its permeability varies with temperature mimicking the brittle-ductile transition of a basaltic rock.

The potential supercritical resource is targeted by a well or a group of wells, direct production from and cold water injection into the geothermal reservoir is modelled. To do so we have implemented a Peaceman-like well model for the 3D CVFEM. We benchmarked our numerical well model against an analytical solution and tested it for mesh size convergence at supercritical conditions. Our reconnaissance simulations provide the first glimpse into the response of a supercritical resource to operation with wells.

[1] P. Weis et al., Geofluids 2014, 14, 3.
[2] S.K. Matthäi et al., Geological Society Special Publication 2007, 292.
[3] T. Driesner, C.A. Heinrich, Geochimica et Cosmochimica Acta 2007, 71.

The Ravn oilfield is located in licences 5/06 and 2/16 in the Central Graben of the Danish Offshore. It is a low permeable Upper Jurassic HP/HT oil reservoir which is developed using long horizontal wells with multiple hydraulic fractures. The Ravn oil is a light oil however with asphaltene in it. SARA analysis of the oil sample from one of the exploratory wells showed 3.3 wt% of asphaltene. The asphaltene inhibitor was injected downhole to mitigate potential asphaltene issues in the well tubing and surface facilities. The reservoir aquifer support was deemed very weak due to the low reservoir permeability and reservoir compartmentalization. Pressure maintenance, i.e., water injection, was considered not viable for the same reasons. When reservoir pressure drops below Asphaltene Onset Pressure (AOP), asphaltene will precipitate and likely cause formation plugging, especially when a reservoir has a low permeability sand with small pore throats. The lab AOP measurements were deemed suspicious due to the unrepresentative oil samples from the production well. The modelled AOP was estimated instead after the field was brought in production with the help of surface oil sampling data and production well pressure response.
In order to simulate the asphaltene formation damage which occurred in the reservoir, a characterized fluid was developed in the PVT modelling software PVTsim. The Equation of State (EOS) parameters were calibrated with the PVT experiments performed on the oil sample from one of the exploratory wells. The AOP and asphaltene solution with pressure variation were modelled with the characterized fluid model in PVTsim. Thereafter the characterized fluid model was exported to the dynamic reservoir modelling software tNavigator which utilized relevant asphaltene formation damage keywords together with the exported compositional model to simulate the asphaltene precipitation, deposition and reservoir permeability impairment. Single parameter analysis showed flocculation rate and permeability reduction versus asphaltene deposits saturation table were the most sensitive parameters in achieving a good history match.
So far, there has been no sufficient production data to pinpoint the location of the asphaltene formation damage. Two alternative scenarios were simulated: (1) asphaltene formation damage happened in both reservoir and hydraulic fractures; (2) asphaltene formation damage happened only in the reservoir. Dynamic models based on the two scenarios could both match the well production history quite well, however the production forecasts based on the two scenarios were different and indicated a range of uncertainty.

Surface complexation reactions are frequently used to describe the adsorption of charged ions onto mineral surfaces. However, they can be challenging to implement robustly into numerical software, especially when coupling geochemical equilibrium calculations to transport.

To compensate for the buildup of charge at a mineral surface, the composition of the electric diffuse layer next to the surface must balance the net charge of the surface. The mathematical relationship between surface charge and surface potential is obtained from the Poisson-Boltzmann equation. To calculate the composition of the diffusive layer explicitly, a time-consuming numerical integration must be performed inside each grid block at each Newton iteration. In addition to slowing down the simulations, the implementation in popular simulation software such as PHREEQC can sometimes fail to converge.

In this paper we show that for electrolytes containing neutral, monovalent, divalent ions and/or complexes with valence less than three (i.e., most electrolytes encountered in practice), we can derive analytical expressions for the electrostatic terms in the governing equations. This allows us to determine the composition of the diffusive layer analytically, without evaluating the integral numerically. To our knowledge, this formulation is novel. For the cases tested so far, the new method greatly increases numerical stability and speed of convergence during Newton-Raphson iterations. Even in cases where diffuse layer concentrations are not required directly, incorrect solutions may be found. The relationship between surface charge and surface potential is then given by the Grahame equation, which has two solutions, only one of which is physical: the solution where surface charge and surface potential have the same sign. Converging to the correct solution during ordinary Newton-Raphson iterations is not guaranteed but depends strongly on the choice of initial guess for the solution to the geochemical system. To avoid this complication, we use our new formulation to derive a version of the Grahame equation having only a single solution, the physically correct one.

Finally, we present a model that allows for the preferential accumulation of ions in the diffusive layer. The model, which is implemented mathematically by including ion exchange sites with variable exchange capacity, is consistent with recent data for ion transport in chalk.