Reservoir-scale transdimensional fracture network inversion

The Waiwera aquifer hosts a structurally complex geothermal groundwater system, where a localized thermal anomaly feeds the thermal reservoir. The temperature anomaly is formed by the mixing of waters from three different sources: fresh cold groundwater, cold seawater and warm geothermal water. The stratified reservoir rock has been tilted, folded, faulted, and fractured by tectonic movement, providing the pathways for the groundwater. Characterization of such systems is challenging, due to the resulting complex hydraulic and thermal conditions which cannot be represented by a continuous porous matrix.

By using discrete fracture network models (DFNs) the discrete aquifer features can be modelled, and the main geological structures can be identified. A major limitation of this modelling approach is that the results are strongly dependent on the parametrization of the chosen initial solution. Classic inversion techniques require to define the number of fractures before any interpretation is done.

In this research we apply the transdimensional DFN inversion methodology that overcome this limitation by keeping fracture numbers flexible and gives a good estimation on fracture locations. This stochastic inversion method uses the reversible-jump Markov chain Monte Carlo algorithm and was originally developed for tomographic experiments. In contrast to such applications, this study is limited to the use of steady-state borehole temperature profiles – with significantly less data. This is mitigated by using a strongly simplified DFN model of the reservoir, constructed according to available geological information.

We present a synthetic example to prove the viability of the concept, then use the algorithm on field observations for the first time. The fit of the reconstructed temperature fields cannot compete yet with complex three-dimensional continuum models, but indicate areas of the aquifer where fracturing plays a big role. This could not be resolved before with continuum modelling. It is for the first time that the transdimensional DFN inversion was used on field data and on borehole temperature logs as input.

1 Introduction

Aquifer systems in fractured rocks can be modelled by using either continuum or discrete models. Continuum models substitute the fractured rock with an equivalent porous media (Day-Lewis et al., 2000; Hao et al., 2008; Illman et al., 2009; Illman and Neuman, 2003; Sahimi, 2011; Vesselinov et al., 2001; Zha et al., 2015). These methods often fail to model the exact location of the fractures due to their strong averaging behavior and they are considered most suited to problems with high fracture density (Long et al., 1982).

Discrete fracture network models (DFNs) are more realistic representations of the fractured media, where fractures and the rock matrix are modelled separately (Dorn et al., 2013; Jang et al., 2008; Neuman, 2005; Niven and Deutsch, 2012). This separation allows to solve flow and transport problems effectively, limiting these expensive computations to the open fracture space. The use of DFN models for inversion purposes however is not as straightforward as classic equivalent porous medium models.

DFN models behave in a nonlinear, non-continuous way, as small changes in fracture locations could open new pathways for flow, or completely isolate large parts of the domain. Fractures are parametrized separately, leading to models with large numbers of parameters, leading to ill-posed inversion problems. And this number is usually an unknown, making any classic inversion approach impossible to use directly.

Substituting the complex DFN model with a simplified model has been proven a good solution to this problem. Simplified DFNs limited to a few main fractures were used to interpret hydraulic, tracer and temperature observations in several studies (Le Borgne et al., 2007; Le Goc et al., 2010; Klepikova et al., 2013, 2014). The simplifications made possible the identification and characterization of the preferential flow paths, but could not reconstruct complex fracture patterns.

The use of complete DFN models is typically limited for forward modelling purposes, or to the inversion of statistical DFN parameters instead of exact geometries. The unknown number of model parameters was a main obstacle, as standard inversion methods require a predefined parameter set to adjust.

Hestir et al. (1998) and Niven and Deutsch (2012) used activation-deactivation of fractures to keep the number of parameters fixed, while Dorn et al. (2013) stochastically generated fracture networks with a fixed number of fractures. Hui et al. (2019) and Pan et al. (2016) used the statistical parameters of the fracture network in combination with connectivity to model the dynamic behavior of fractured reservoirs. Yao et al. (2018) changed the parametrization of the fracture network into a set of continuous parameter fields by using Hough-transformation, making them suitable for inversion methods requiring Gaussian distributions.

Somogyvári et al. (2017) implemented a transdimensional DFN inversion approach, that resolved the problem of parameter numbers by keeping them flexible throughout the inversion. The method was capable to explore different DFN geometries and find the preferential transport trajectories. The method has been used to interpret synthetic tracer tomography experiments, and to characterize cross-borehole aquifer conditions.

The presented study uses this transdimensional DFN inversion algorithm to characterize the geothermal reservoir of Waiwera in New Zealand. After introducing the field site, we show how the methodology needed to be modified for applicability in a large-scale non-tomographic setting using temperature borehole data. We demonstrate the applicability on a synthetic example for proof-of-concept, then present the results of the interpretation of the field data, to infer the fault geometry of the Waiwera geothermal reservoir. This is the first application of this methodology on field observations.

2 Field site

Waiwera is a small east coastal township north of Auckland on the northern island of New Zealand. A low temperature (max. 50 ^∘C) geothermal reservoir is located under the city, making the area a popular recreational resort. At this unique location, artesian flow from bores occurred until the end of the 1960s and thermal springs were emanating next to oceanic beaches until the mid of the 1970s. As a result of overproduction the warm water has to be pumped from the reservoir since then (Kühn and Stöfen, 2005).

A detailed description of the local geological and hydraulic setting can be found in (ARWB, 1980).

In the aquifer the upflowing geothermal water is mixing with fresh groundwater and seawater. These complex flow conditions are forming a steady state temperature anomaly at the site, which is shown in Fig. 1b. The site is well studied, as well based on borehole information (Kühn and Altmannsberger, 2016).

Figure 1 Geological and thermal conditions of the Waiwera geothermal site. (a) Structural elements of the aquifer, (b) interpolated temperature distribution.

The aquifer body is formed of the Waitemata sandstone formation, which is a strongly stratified sandstone with a hydraulic conductivity of $\mathrm{2.3}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ (average of entire setting determined by a pumping test). The spacing of the stratification is unknown (Kühn et al., 2016). The stratification is not horizontal, because the formation has been folded, fractured and faulted (Kühn and Schöne, 2017). It is also intersected by two high inclination fault sets (Fig. 1a). One of these faults provide the pathway for the geothermal water at the aquifer bottom. From hydraulic tests and well drill activities the hydraulic conductivity of faults and fractures is estimated (in the range of 10⁻⁴–10⁻⁶ m s⁻¹) and used in the inversion (ARWB, 1987).

Previous studies modelled the aquifer using continuum approaches, combining the hydraulic properties of the fractures, the stratification and the aquifer matrix (Kühn et al., 2016; Kühn and Altmannsberger, 2016; Kühn and Stöfen, 2005). These models however do not include the discrete features and could not describe the thermo-hydraulic behavior of the aquifer. The aim of this study is to interpret the borehole temperature profiles with DFN models, to infer the geometry of the faults.

Figure 2 The transdimensional DFN inversion algorithm. Every iteration of the inversion consists of two phases. In the model update phase, the DFN is modified randomly, via fracture deletion, addition or horizontal translation of a fault. In the model evaluation phase, the updated model is evaluated by simulating the temperature field in the aquifer.

3 Methodology

3.1 Inversion

In this study we use the transdimensional DFN inversion introduced by (Somogyvári et al., 2017), based on the reversible-jump Markov chain Monte Carlo (rjMCMC) algorithm (Green, 1995). This iterative technique uses subsequent random geometry perturbations to generate and test possible model configurations (Fig. 2). Every iteration consists of two phases: an update and an evaluation phase.

In the update phase, the current DFN realization (θ_n) is altered by fracture addition, deletion or movement. The type of the update and its properties are chosen randomly, with a predefined probability ($q\left({\mathit{\theta }}_n\to {\mathit{\theta }}^{\ast }\right)$). Fractures could only be inserted to specific discrete points along the existing DFN. This is an important condition to preserve the reversibility of the chain (see Somogyvári et al., 2017).

In the evaluation phase the temperature distribution in the aquifer is modelled on the updated DFN realization (θ^∗) using the forward model (details in Sect. 3.2). The model evaluation is done using the Metropolis-Hastings-Green acceptance criterion (Green, 1995):

$$\begin{array}{}\text{(1)}& \mathit{\alpha }\left({\mathit{\theta }}^{\ast }\mathrm{|}{\mathit{\theta }}_n\right)=min\left[\mathrm{1},{\displaystyle \frac{p\left({\mathit{\theta }}^{\ast }\right)}{p\left({\mathit{\theta }}_n\right)}}{\displaystyle \frac{L\left(\mathit{\xi }\mathrm{|}{\mathit{\theta }}^{\ast }\right)}{L\left(\mathit{\xi }\mathrm{|}{\mathit{\theta }}_n\right)}}{\displaystyle \frac{q\left({\mathit{\theta }}^{\ast }\to {\mathit{\theta }}_n\right)}{q\left({\mathit{\theta }}_n\to {\mathit{\theta }}^{\ast }\right)}}\left|J\right|\right].\end{array}$$

The likelihood-ratio ($L\left(\mathit{\xi }\mathrm{|}{\mathit{\theta }}^{\ast }\right)/L\left(\mathit{\xi }\mathrm{|}{\mathit{\theta }}_n\right)$) is calculated from the improvement of the RMS error of the simulated observations (ξ|θ). We consider the observations (ξ) normally distributed, hence the L functions have a gaussian shape. The proposal ratio ($q\left({\mathit{\theta }}^{\ast }\to {\mathit{\theta }}_n\right)/q\left({\mathit{\theta }}_n\to {\mathit{\theta }}^{\ast }\right)$) is calculated as the ratio of the probabilities of the reverse and the forward update. Because these updates are discrete, the calculation of these probabilities are straightforward. For example, the probability of a fracture deletion is one over the total number of fractures. A detailed description on the update probabilities can be found in (Somogyvári et al., 2017). The so-called Jacobian (|J|) is a consequence of the transdimensional updates, and it is required to maintain the stability of the rjMCMC. With discrete model updates however, its value is always one (Denison et al., 2002).

If a proposal gets accepted the new iteration starts with it, if it gets rejected the chain continues with the previous realization. The final product of rjMCMC inversion is not one calibrated model, but a set of realizations often referred to as ensemble. The ensemble is suitable for further statistical investigations and for uncertainty analysis (Somogyvári et al., 2017). To eliminate the effect of the arbitrary chosen initial model, the first part of the chain is discarded and not used in the ensemble (Somogyvári and Reich, 2019).

3.2 Forward model

To simulate the steady state temperature distribution in the DFN model, the tracer transport simulator developed by Jalali (2013) was used. This efficient algorithm could simulate temperature distributions in two dimensional DFN models within a few seconds, making it suitable for MCMC applications.

The forward model takes the DFN geometry and hydraulic and temperature boundary conditions as an input and returns the steady state temperature distribution in the fracture network. The simulator is transient, both the pressure and temperature simulations start from an initial field with perturbations by the boundary conditions, then the simulation is run until a steady pressure and temperature field is obtained.

First the DFN model is discretized to 1-D segments and nodes. Pressure distribution is simulated by using mass conservation, with an implicit finite difference method. Flow in the fractures are calculated by Darcy's law from the pressure gradient. For thermal transport, the advection-dispersion equation is solved considering a steady-state flow field, with an implicit upwind finite difference solver (Jalali, 2013; Ringel et al., 2019).

In this simple model no temperature transport is simulated within the impermeable rock matrix, only the temperature loss from the fracture space towards the rock matrix is simulated. Hence, cross-fracture heat transport (between non-connected near fractures) is not possible. The simulation is limited to 2-D and no complex hydraulic and thermohaline boundary conditions are considered.

The temperature field inside the non-fractured aquifer matrix is calculated via linear interpolation with triangulation. The interpolated field is then intersected with the boreholes, and the extracted 1-D temperature profiles were used as the simulated observations. In theory it would be possible to compare the interpolated 2-D temperature field, with the interpolated field from the site (see Fig. 1b), but we wanted to demonstrate the robustness of the method with only using temperature logs from a limited number of boreholes. The complete simulation could be done faster than a second on a laptop, thus this simple DFN transport model is efficient and fast to be used in an MCMC framework.

3.3 Conceptual model

The methodology introduced in (Somogyvári et al., 2017) needed to be modified to work on reservoir scale scenarios. The main modification is a simplification to mitigate the limited amount of observations available. Three fracture sets are defined: two for the major fault directions and one for the stratification. In this modified setting, fractures are completely intersecting the domain. This simplification is done based on the geological model of the aquifer (Fig. 1a). Hence, fracture lengths are not a free parameter of the inversion as they are calculated from the fracture position within the domain. In addition, fracture centers are aligned to the vertical or horizontal centerline of the domain (depending on the fracture set) – basically reducing the free parameters to one per fracture.

Some further modifications were also done. Fracture addition and deletion is only used for the low-inclination fracture set (stratification), and fracture movement only for the high inclination sets (the exact location of these faults is not known).

The hydraulic and thermal fracture properties were chosen after available borehole information. The two fault sets are simulated with an aperture of 0.5 mm. The stratification aperture is set to 0.8 mm. These values are chosen after (Kühn et al., 2016) using the cubic law to convert transmissivities to apertures. The values were further tuned by preliminary modelling trying to obtain similar extent of the temperature anomaly as the observed one. These tests also showed that the anisotropy has a stronger effect than the exact aperture values. The hydraulic conductivities of the fractures are calculated from the apertures using the cubic law. Note that both flow and heat transport is limited to 2-D in the model, and no 3-D effects are considered.

Figure 3 Boundary conditions of the forward modelling.

Boundary conditions are defined on the top and bottom borders (Fig. 3). In the DFN forward models, they appear when a fracture intersects with these borders. Note that the used pressure values are not absolute and were selected after the hydraulic gradient. For the source of the geothermal water, a 100 m long source region with increased pressure and temperature is defined in the bottom center of the model. The observed maximum temperature of 50 ^∘C was assigned. Pressure of the source region was arbitrarily selected after the parameter tests with inspecting the extent of the temperature anomaly (Fig. 3).

Note that the used methodology is limited to invert the geometry, and does not estimate the physical fracture properties. Introducing the hydraulic parameters of the fractures as free parameters into the inversion would result in an increase in model freedom, which without involving additional data would lead to different equivalent solutions (Somogyvári et al., 2017).

Figure 4 Reconstruction of the synthetic DFN example: (a) Synthetic DFN model, (b) Reconstructed fracture probability map.

4 Results

4.1 Synthetic modelling

The method was first tested on a synthetic dataset. The chosen synthetic DFN model is shown in Fig. 4a. The inversion here is limited to the fractures representing the stratification. The high inclination faults are considered known, and not estimated by the algorithm.

The transdimensional DFN inversion was implemented in python and was run on a standard office laptop with an Intel i5-7267U CPU. The chains were run for 1000 iterations which took 5 min on average. Because of the stochastic nature of the used inversion technique, we ran several simulations with the same initial parameters to see the stability of the results. All these simulations led to reconstructed geometries with the same general features, showing the robustness of the algorithm.

The result of the synthetic example is shown in Fig. 4b in the form of a fracture probability map. This plot visualizes the realizations of the rjMCMC chain by rasterizing the DFN models first, then counting the pixels with fractures along the ensemble. This is equivalent of taking the mean of the MCMC chain for visualization (Bodin and Sambridge, 2009). The fracture probability map shows perfect fits with the two fractures in the bottom part of the profile. The three fractures in the top are not captured exactly, but the extension of this fractured zone is well resolved. The absence of fractures in the center part is captured perfectly. This result is a proof of concept that the methodology is applicable on the proposed geological problem.

4.2 Field modelling

For field application, the three virtual boreholes defined as 1-D slices of the profile shown in Fig. 1b were chosen. One borehole is taken from the center of the anomaly, the other two from the two sides. These virtual boreholes were selected instead of real wells for practical reasons, as all temperature profiles from the site were already compiled and integrated to a 3-D temperature dataset. This removes any local effects of the actual boreholes and also provide room for sensitivity analysis of different borehole locations numbers along the profile, which is to be done in the future.

To improve the estimation, background temperature trends were removed from the virtual boreholes (using the temperature trends from the two sides of the model). With the removal of these trends, we have also mitigated the effect of freshwater-seawater mixing, which is not simulated by the forward model.

Simulations were run for 1000 iterations, and similarly to the synthetic case the first half of the chain was discarded. The simulations were completed within 15 min, a little longer then the synthetic case. The results are presented in Fig. 5.

Figure 5 (a) Fracture probability map of the Waiwera geothermal aquifer, (b) mean temperature field of the reconstructions (c) observed borehole temperature profiles (black) and reconstructed borehole profiles (gray).

The reconstructed fracture probability map has the following three distinct features. Strong stratification is present close to the bedrock with high probability. Fractures in this area are close to the source and provide pathways for the heat to spread out horizontally. A non-stratified aquifer body present in the center of the domain. It could not only mean that this part is non-fractured, this could be also explained with the existence of fractures that do not play a role of forming the anomaly. This is supported by the similarity of additional synthetic examples where this feature was usually present. In contrast, there is high probability of stratification in the top. The role of these fractures is to close the anomaly, by allowing mixing with the colder waters. The inversion also placed faults close to the location of two boreholes for similar reasons.

Figure 5b shows the mean of the reconstructed temperature fields of the ensemble. Compared to the observed temperature anomaly, the reconstruction shows similar properties. The two anomalies have a comparable extent horizontally and vertically. The skewedness of the anomaly is reconstructed. Originally this was explained by the mixing of different waters, but our results show that it could be explained by the structure of the aquifer.

Figure 5c shows the fit quality of the borehole temperature profiles. Continuous models from the same site provide better fits for most of the temperature observations (Kühn and Stöfen, 2005) but cannot match some observations which are supposed to be based on the fractured characteristic of the reservoir rock. The main limitations here are the 2-D simulation and the lack of complex geometry options, which restricts the shape of the reconstructed anomaly. This could be addressed by generating and testing more complex fracture patterns, but can only be executed by involving additional data into the inversion.

5 Conclusions

Temperature observations are usually considered as complementary information for thecharacterization of aquifer structure (Bravo, 2002; Klepikova et al., 2014; Leaf et al., 2012). In this paper we have presented a methodology that treats steady-state borehole temperature observations as the main driver of the interpretation process. We have shown that a handful of these profiles could provide enough information to resolve the main structural elements of the investigated aquifer using a DFN based transdimensional inversion.

We used the inversion algorithm introduced by Somogyvári et al. (2017) which provides a data driven approach by keeping the number of the modelled fractures flexible throughout the inversion. The DFN model is highly conceptualized and suitable for solving the ill-posed inversion problem, given the limited amount of observations.

Using synthetic models, we have proven that the transdimensional DFN inversion, which was originally developed for tomographic observations, could be used in this significantly different scenario.

We also present the first field application of the methodology, for the characterization of the well-studied Waiwera geothermal reservoir. Transdimensional DFN inversion indicated possible locations of strongly stratified aquifer parts, and areas which are less likely to be fractured. These findings are to be validated by future explorational campaigns whose results will be also used to infer DFN models with more complexity.

Compared to other interpretations from the site using three dimensional porous medium models, the shape of the reconstructed anomaly, and the fit of the temperature profiles are less accurate. The simplified DFN model does not consider flow and transport in the rock matrix, it does not take into account the complex hydraulic and thermohaline conditions and it is also limited to two dimensions. Still despite these limitations it indicated discrete aquifer features which the continuum modelling was not capable to capture. Thus, for a complete aquifer characterization campaign, we recommend combining the two approaches, by enhancing the three-dimensional continuous models by the discrete features obtained from DFN modelling.

We did not explore the stochastic behavior of the methodology, the resulted ensemble could be the subject of further statistical analysis and uncertainty quantification (Somogyvári et al., 2017; Somogyvári and Reich, 2019). This aspect will be investigated in a future study. The presented interpretations are to be used to select aquifer sections for more detailed exploration, to characterize the hydraulic properties of the fractures and identify the locations of the faults.