Dynamic mesh optimisation for geothermal reservoir modelling

Abstract

Modelling geothermal reservoirs is challenging due to the large domain and wide range of length- and time-scales of interest. Attempting to represent all scales using a fixed computational mesh can be very computationally expensive. Application of dynamic mesh optimisation in other fields of computational fluid dynamics has revolutionised the accuracy and cost of numerical simulations. Here we present a new approach for modelling geothermal reservoirs based on unstructured meshes with dynamic mesh optimisation. The resolution of the mesh varies during a simulation, to minimize an error metric for solution fields of interest such as temperature and pressure. Efficient application of dynamic mesh optimisation in complex subsurface reservoirs requires a new approach to represent geologic heterogeneity and we use parametric spline surfaces to represent key geological features such as faults and lithology boundaries. The resulting 3D surface-based models are mesh free; a mesh is created only when required for numerical computations. Dynamic mesh optimisation preserves the surfaces and hence geologic heterogeneity. The governing equations are discretised using a double control volume finite element method that ensures heat and mass are conserved and provides robust solutions on distorted meshes. We apply the new method to a series of test cases that model sedimentary geothermal reservoirs. We demonstrate that dynamic mesh optimisation yields significant performance gains, reducing run times by up to 8 times whilst capturing flow and heat transport with the same accuracy as fixed meshes.

Introduction

Numerical modelling is an essential tool for the development of geothermal reservoirs. Modelling of these complex and data-poor environments is important to gain insight into reservoir processes as well as predict the response of the reservoir to production. Numerical modelling can be used to optimise well number, placement and geometry for maximum heat recovery, predict thermal breakthrough, and study the interaction of production with background hydraulic flow. Modelling the response to production is essential to assess the lifetime of a given reservoir and identify the best management strategy (Arthur et al., 2010; Blöcher et al., 2010; Herbert et al., 2013; Dussel et al., 2016; Willems and Nick, 2019).

Current approaches to geothermal reservoir simulation can be broadly separated into two categories based on the approach used to discretize space: those that use k-orthogonal structured grids or unstructured meshes, and those that use non k-orthogonal unstructured meshes. In the former, the governing equations are typically discretised using the finite difference method (FDM) or finite volume method (FVM); k-orthogonal grids are required to accurately calculate fluxes between grid elements using the two-point flux approximation (Eskilson and Claesson, 1988, Eskilson and Claesson, 1988; Pruess, 1990; Sliwa and Gonet, 2005; Diersch et al., 2011a, Diersch et al., 2011b; Nabi and Al-Khoury, 2012). Simulators which use the FDM/FVM include TOUGH2 (Pruess et al., 1999), ECLIPSE (Schlumberger, 2018), STARS (CMG, 2019) and TETRAD (ThinkGeoEnergy, 2019).

The use of k-orthogonal grids (typically Cartesian) or close to k-orthogonal grids (typically pillar grids, in which pillars extend from the top to the base of the model domain) often requires excessively high resolution to accurately capture geological heterogeneity (Wu and Parashkevov, 2009, Jackson et al., 2015). Features such as dipping faults, meandering channels or pinched-out layers introduce unrealistic stair-stepping effects and associated numerical artefacts in numerical solutions. The use of pillar grids also restricts flexibility of model resolution, because all layers in the model must have the same areal grid resolution. This can be very computationally expensive when trying to resolve a large range of length-scales within a large model domain, as is often the case when modelling geothermal reservoirs (O’Sullivan et al., 2001; Cho et al., 2015).

The use of non k-orthogonal meshes can therefore be beneficial, as it allows more accurate representation of geologic heterogeneity with fewer mesh elements (Geiger et al., 2004; Paluszny et al., 2007; Jackson et al., 2015; Zehner et al., 2015). Such meshes have been extensively used to model coupled geothermal processes in fractured porous media (Weis et al., 2012, Weis et al., 2014; Scott et al., 2015, Scott et al., 2016; Yapparova et al., 2014; Nissen et al., 2018; Patterson et al., 2018; Patterson and Driesner, 2019; Berre et al., 2019, Berre et al., 2020; Keilegavlen et al., 2020, Keilegavlen et al., 2021). The finite element method (FEM) is typically used to accurately calculate fluxes across the elements. Simulators based on this approach include FEFLOW (FEFLOW, 2019), PANDAS (Xing et al., 2019) CSMP (Nick and Matthai, 2011), FEHM (Zyvoloski et al., 1997) and OpenGeoSys (Kolditz et al., 2012). However, use of the FEM can also involve numerical challenges: solutions are typically more expensive per element than FDM/FVM and poor quality meshes can lead to convergence failures or non-physical solutions. This latter issue can be a significant problem when modelling geothermal reservoirs containing high-aspect ratio geological layers.

Geothermal reservoirs host a range of thermal-hydrological-geomechanical-chemical (THMC) processes and numerical simulators must solve a system of highly-coupled, non-linear partial differential equations which is very computationally expensive. To reduce cost, small domains can be modelled or rock properties can be upscaled to a coarser mesh. However, these approaches can both introduce new problems. Large domains may be necessary to avoid boundary effects and to include the whole convective system hosting a particular reservoir, while upscaling can result in a loss of geological fidelity (Cho et al., 2015; Burnell et al., 2015; O’Sullivan et al., 2001; Renard and De Marsily, 1997).

In other areas of computational fluid dynamics, high fidelity solutions have been obtained at lower computational cost by use of dynamic mesh optimisation (DMO), in which the resolution and geometry of the mesh varies during a simulation to minimize an error metric for one or more solution fields of interest such as pressure or velocity (Alauzet and Loseille, 2016; Xie et al., 2016; Hu et al., 2018). DMO varies the mesh resolution such that higher resolution is used in parts of the domain where the solution is complex, and lower resolution is used elsewhere. However, DMO has not been applied in geothermal reservoir simulation, and it is not trivial to adapt existing DMO methods because of the need to maintain a description of the underyling geologic heterogeneity. This heterogeneity causes spatial variability in rock properties such as permeability and porosity which can have a significant impact on flow.

Modelling the spatial variability of rock properties on a fine mesh can allow geologic heterogeneity to be properly captured, but is computationally expensive on a fixed mesh and does not allow efficient application of DMO because, each time the simulation mesh changes, the rock properties must be up-, cross- or down-scaled onto the new mesh. This scaling is trivial for volumetric properties such as porosity, but non-trivial for the tensor properties of fluid permeability and thermal conductivity, for which no straightforward scaling methods are available (Renard and De Marsily, 1997; Rühaak et al., 2015).

Here, instead of the classical grid-based approach to geological reservoir modelling, a surface-based approach is used (Pyrcz et al., 2005; Caumon et al., 2009; Jackson et al., 2015; Melnikova et al., 2016; Jacquemyn et al., 2019). In surface-based geologic modelling (SBGM), all geological heterogeneities of interest are represented by surfaces: these surfaces may capture faults, stratigraphic surfaces, boundaries between different facies or lithologies, boundaries between different diagenetic regions, and any other geologic controls on the spatial distribution of petrophysical properties. The surfaces define rock volumes which we term ‘geologic domains’; these domains have constant petrophysical properties, or simple, mathematically defined trends such as upwards or downwards increases in permeability (reflecting, for example, upwards or downwards coarsening of grain size). In this sense, surface-based modelling is analogous to grid-based modelling: any discretization requires properties to be homogeneous at some scale. Here, we identify the domains based on geology, rather than using an arbitrary mesh. It has been shown that capturing correlated variability using such geologically-defined domains is key for capturing flow (Osman et al., 2020).

The aim of this paper is to demonstrate and test, for the first time, the use of dynamic mesh optimisation in geothermal reservoir simulation, using surface-based modelling to capture geologic heterogeneity and the Double-Control-Volume-Finite-Element-Method (DCVFEM; (Salinas et al., 2017a) to solve the equations governing heat flow on the unstructured, dynamically adapting mesh. The combination of surface-based geologic modelling and DMO allows the unstructured mesh to adapt over time to solution fields of interest such as pressure and temperature (e.g. Fig. 1). Numerical simulations capturing heterogeneity and flow over many length-scales can be performed at lower computational cost and without the need for the user to manually adjust the mesh resolution.

We begin by validating the approach against an analytical solution for advection-diffusion to demonstrate results are accurate and convergent for both fixed and adaptive meshes. We then apply the approach to a number of realistic test cases, including a well doublet in a homogeneous reservoir, production from a sedimentary fluvial reservoir containing sinuous, channelized sandbodies interbedded with mudstones, and production from a complex faulted reservoir with several fault blocks separated by intersecting conjugate faults. These heterogeneous test cases represent common reservoir scenarios that can be difficult to represent in conventional geothermal reservoir modelling workflows. We focus here on thermal-hydrological (TH) processes in the low enthalpy geothermal systems that are used to source and store heat in many different locations around the world (Lopez et al., 2010; Gluyas et al., 2018; Willems and Nick, 2019). Extensions of the approach to include geomechanics and chemical reaction (MC) are discussed.