Crustal Reservoir Flow Simulation for Long-Range Spatially-Correlated Random Poroperm Media

Peter Leary ^1,2,†,*, Peter Malin ^2,†

1. GeoFlow Imaging, 43 High St, Auckland 1010, New Zealand

2. Advanced Seismic Instrumentation & Research, 317 E. Millbrook, Raleigh, NC 27609, USA

† These authors contributed equally to this work.

* Correspondence: Peter Leary

Academic Editor: Zhao Yang Dong

Special Issue: Modeling of Geothermal Systems

Received: November 30, 2020 | Accepted: March 15, 2021 | Published: March 19, 2021

Journal of Energy and Power Technology 2021, Volume 3, Issue 1, doi:10.21926/jept.2101013

Recommended citation: Leary P, Malin P. Crustal Reservoir Flow Simulation for Long-Range Spatially-Correlated Random Poroperm Media. Journal of Energy and Power Technology 2021; 3(1): 013; doi:10.21926/jept.2101013.

© 2021 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

The fundamental performance feature of geothermal energy extraction is production wellbore connectivity to spatially erratic convective cell flow structures. Wellbore power given by Q ~ ρC x T x V for fluids of volumetric heat capacity ρC scales with resource temperature T and wellbore volumetric fluid flow rate V = 2πr₀φv₀ℓ. For standard values of wellbore radius r₀ and open hole wellbore length ℓ, geothermal well performance depends directly on formation bulk inflow rate φv₀. Assuming adequate reservoir porosities φ(x,y,z), Darcy’s law gives wellbore power Q in terms of formation permeability structure at the wellbore, v₀ = k(x₀,y₀,z₀)/μ ∂_rP(r)|_r0, with ∂_rP(r)|_r0 = P₀/r₀/log(R/r₀) for suitable flow system radius R. Geothermal heat energy production rate Q is thus irreducibly controlled by permeability structures in the vicinity of the wellbore.

As we detail below, the empirical spatial variability of ambient crust permeability can be expressed by the mathematical form k(x,y,z) ~ exp(αφ(x,y,z)) within a radius R, √((x-x₀)² + (y-y₀) ² + (z-z₀)²) < R, surrounding a wellbore site at location (x₀, y₀,z₀). In this mathematical form, porosity φ(x,y,z) is normally distributed between φ_mn < φ(x,y,z) < φ_mx and is spatially correlated across cm-to-km length scales. The normally distributed porosity field generates an associated lognormally distributed permeability field for large empirical values of poro-connectivity parameter α. Small poroperm parameters α reduce spatial permeability spatial fluctuations to k(x,y,z) ~ 1 + αφ(x,y,z) that are consistent with thermal property fluctuations [1]. Such permeability structures are compatible with standard numerical flow simulation procedures developed for heat conduction in solids [2,3,4,5]. Empirical α values for the ambient crust are, however, of order 20 to 30 as set by the condition 3 < αφ < 5 for φ ~ 0.2 mean crustal reservoir porosities [6]. As such they are compatible with lognormal well-productivity data observed for ambient crustal geofluids worldwide. In consequence of empirical values 20 < α < 30, large-scale spatially-rapid variations in crustal permeability vastly exceed the magnitude and spatial complexity of rock thermal property distributions, thus destabilising standard crustal fluid flow modelling procedures.

To address numerical stability problems arising from large amplitude and rapidly varying ambient crustal permeability fluctuations, we here validate two non-traditional flow modelling procedures. These numerical procedures allow flow modelling to accurately connect geothermal production wellbore interval inflow velocity v₀ at wellbore sites (x₀,y₀,z₀) to long-range spatially erratic poro-connectivity permeability structures exemplified by the empirical expression k(x,y,z) ~ exp(αφ(x,y,z)). Accurate wellbore inflow modelling of complex convective geothermal flow systems can interact efficiently with large-scale flow imaging data acquired by surface seismic sensor arrays to guide production well drilling [7]. The driver for accurate wellbore geothermal system flow modelling lies in the lognormal distribution of well productivity throughout the ambient crust [6]. It is proving uneconomic for convective geothermal energy production to drill large numbers of unproductive wells before locating productive and sustainable flow structures for targeted drilling [7]. Using advancing observational surface seismic sensor array data processing to remotely identify productive crustal flow systems will significantly reduce both exploration and exploitation drilling costs. An emerging use for precision flow modelling in ambient crustal volumes lies in tracking permeability enhancement (EGS) by systematic injection of wellbore fluids into deep hot crustal volumes [8,9].

Nomenclature

Q ≡ Heat energy production rate in watts for a wellbore

ρC ≡ Volumetric heat capacity of wellbore fluid in joules per cubic meter of fluid

T ≡ Centigrade temperature of geothermal fluid entering a production wellbore

V ≡ Velocity of wellbore fluid in cubic meters per second

r₀ ≡ Radius of open production wellbore in meters

v₀ ≡ Darcy fluid flow velocity at wellbore radius in meters per second

ℓ ≡ Axial length of open production wellbore in meters

(x₀,y₀,z₀) ≡ Central location of a production well open interval in meters

φ(x,y,z) ≡ Spatially variable porosity of geothermal flow system near wellbore

k(x,y,z) ≡ Spatially variable permeability of geothermal flow system near wellbore

μ ≡ Dynamic viscosity of wellbore fluid in Pascal-seconds

P(r) ≡ Radial pressure of convective fluid near the wellbore in Pascal

R ≡ Effective radius in meters of open wellbore connectivity to convective geothermal system

α ≡ Empirical proportionality coefficient relating changes of logarithm of well-core permeability to changes in well-core porosity, δlogκ ~ αδφ

φ ≡ Mean porosity of crustal volume of radius R hosting a geothermal flow structure

αφ ≡ Critical state product obeying empirical condition 3 < αφ < 5; for reservoir porosity range 0.1 < φ < 0.3, 15-20 < α < 25-30

β ≡ Exponent of well-log power-law scaling Fourier power-spectra, S(k) ∝ 1/k^β; β ~ 0 ≡ white noise; β ~ 1 ≡ pink noise; β ~ 2 ≡ red noise.

2. Materials & Methods-Computing Fluid Flow For Permeability k(x,y,z) ~ exp(αφ(x,y,z))

The fluid flow material properties of ambient crustal rock are set by a trio of widely attested empirical field relations (i)-(iii) derived respectively from crustal well-log, well-core, and well-productivity data. Closely associated with these crustal poroperm statistical and spatial distribution empirics are microseismic event empirics (iv)-(vi) observed in crystalline basement, geothermal fields, and active deformation zones. Building on the appropriate flow simulation methodology, we can accurately model crustal reservoir flow systems derived from microseismic emission data embedded in numerical realisations of ambient crustal permeability structures κ (x,y,z) ~ exp(αφ(x,y,z)) as mandated by the following ambient crustal flow structure empirics:

1. Well-log spatial sequences of rock properties establish the existence of long-range power-law-scaling spatial correlations as a fundamental feature of the ambient brittle-fractured crust. Well-log Fourier power-spectra, S(k) ∝ 1/k^β, β ~ 1, span five decades of spatial frequency, 1/cm < k < 1/km [10,11,12]. Of particular interest are neutron porosity logs from hydrocarbon reservoirs, geothermal formations, and basement rock at up to 6km depth. For contrast, the de facto standard assumption about crustal heterogeneity is lack of significant spatial correlation as denoted by ‘’white noise’’ (β = 0) power-spectra, S(k) ∝ 1/k⁰ ~ const. 
2. Well-core data from hydrocarbon reservoirs, geothermal formations, and deep basement rock show close spatial correlation systematics, δlogκ ~ αδφ, for spatially varying porosity φ and permeability κ [6]. On well-core evidence, rock permeability is a physical feature of fracture-connectivity between spatially-correlated porosity structures at all scales. The constant α in the well-core poroperm spatial correlation is given by the empiric condition αφ ~ 3-5 across two decades of mean porosity, 0.3% < φ < 30% spanning the basement-to-reservoir formation spectrum [13,14]. The well-core empiric for porosity and permeability follows from considering that connectivity between n pores scales as n-factorial, n! = n(n-1)(n-2)….1, so that changes in pore count n → n+δn lead to changes in the logarithm of permeability δlogn! ~ log((n+δn)!) - log(n!) ~ δn log(n), giving the observed well-log poroperm relation. The constant of proportionality α governs the degree to which pores are connected, with high degrees of pore connectivity leading to channelised flow structures.
3. Well-productivity lognormality follows directly from the integrated poroperm relation κ ~ κ₀ exp(αφ). The empirical condition αφ ~ 3-5 guarantees the lognormal well-productivity distributions observed everywhere for groundwater aquifers, conventional and unconventional hydrocarbon reservoirs, convective geothermal systems, and mineral deposition abundances in fossil flow systems [15,16].
4. Ambient crust microseismicity event magnitude distributions in the magnitude range -1 < m < 2 are lognormal rather than power-law/fractal, i.e., the standard Gutenberg-Richter relation [17,18]. The observed lognormal magnitude distribution is congruent with crustal permeability lognormality (iii), indicating that ambient crustal microseismic slip events are closely associated with crustal permeability structures.
5. Ambient crust permeability and microseismic event separations are spatially correlated according to two-point spatial correlation function Γ(r) ~ 1/r^1/2 [19,20]. The similar spatial correlation functions and scaling exponents indicate that microseismic slip events tend to occur on permeability structures rather than hypothetical fault structures.
6. Ambient crust basement rock microseismic event seismic waveform spectra are dominated by high frequencies, 200 Hz < f < 800 Hz, characteristic of seismic slip events on spatially complex crustal permeability structures rather than on hypothetical planar fault-like dislocation structures [9].

Rates of wellbore heat energy extraction Q from convective system volumetric flow V =2πr₀φv₀ℓ thus depend explicitly on the details of production wellbore intersection flow v₀ given by Darcy’s law, v₀ = κ/μ ∂_rP(r)|_r0, applied to the open hole length ℓ within the complex crustal flow structures κ (x,y,z) ~ exp(αφ(x,y,z)) of radial dimension R. For accurate wellbore flow modelling, computational methods for convective geothermal flow systems must be able to robustly handle the high-amplitude multi-scale spatial variations arising from the crustal flow empiric κ (x,y,z) ~ exp(αφ(x,y,z)).

Our crustal fluid flow modelling concern readily emerges from including the crustal poroperm, empiric κ(x,y,z) ~ exp(αφ(x,y,z)) into Darcy’s law (1) and the conservation of fluid mass in steady state flow (2). Darcy’s law gives fluid flow velocity v in terms of the spatially variable formation permeability κ, fluid viscosity μ (here constant), and spatial gradient of fluid pressure P,

\[ \boldsymbol{v}(x,y,z)=\frac{\kappa(x,y,z)}{\mu}\mathbf{\nabla}P(x,y,z)\tag{1} \]

Steady-state reservoir flow models are mathematically derived from conservation of fluid mass constraint equation given by the spatial differential condition on fluid flow velocity requiring that zero net fluid moves into and out of any given crustal unit volume,

\[ \nabla\boldsymbol{v}(x,y,z)=0\tag{2} \]

Along with appropriate flow system boundary conditions, mass conservation condition (2) gives a closed mathematical system for numerical solution. For κ(x,y,z) = exp(αφ(x,y,z)) in (1)-(2), the crustal flow system poro-connectivity parameter α defines the flow modelling computational task as finding pressure fields P(x,y,z) that solve the equation

\[ \nabla^2P(x,y,z)+\alpha\mathbf{\nabla}\varphi(x,y,z)\cdot\mathbf{\nabla}P(x,y,z)=0\tag{3} \]

For small α and/or porosity spatial variations, the governing flow differential equation (3) is close to the mathematically and computationally well-known and highly tractable Laplace equation ∇²P(x,y,z) = 0 and its generalisation the Poisson equation with known source term ∇²P(x,y,z) = f(x,y,z) [2]. In consequence of the numerical stability offered by small values of α and/or small spatial variations in porosity, it quickly became routine to follow the lead of heat conduction numerical modelling in using the established ‘’strong’’ methodology of finite-difference numerical approximation to poroelastic media to model crustal flow [4,5].

Of historical note in ambient crust flow computations are the introduction and persistence of the representative elementary volumes (REVs) [21] derived from effective media approximations such as dual porosity/permeability [22,23,24]. These numerical constructs of crustal flow simulation are based on pre-emptive smoothing of poroperm spatial variability, i.e., the artificial reduction or elimination of troublesome porosity structure spatial gradients [25,26]. Such pre-emptive smoothing is typically justified on the assumption that spatial variations in crustal flow properties are spatially uncorrelated and therefore subject to smoothing through spatial averaging. Where substantial spatial flow heterogeneity is required in a reservoir model, the heterogeneity is taken to define a set of distinct (and suitably uniform) formations that are sutured together in an overall flow structure. As illustrated by the widely used Tough2 flow simulation code, ‘’strong’’ formulation method finite-difference numerical flow solvers based on pre-emptive generation of numerically tractable flow property distributions remain the mainstay of present reservoir fluid flow simulation [27,28,29,30].

Pre-emptive spatial averaging of reservoir flow property distributions for purposes of finite difference flow modelling numerical stability contradicts what we know about spatially correlated flow and wellbore fluid production. While spatially uncorrelated poroperm distributions have flat or constant spectra [12], the crustal empiric (i) shows that wellbore log spatial fluctuations, e.g., neutron porosity, have Fourier power-spectra that scale inversely with spatial frequency, S(k) ∝ 1/k^β, β ~ 1, over the five decades 1/cm < k < 1/km. Further, as given by crustal empirics (ii)-(iii), spatially correlated porosity generates wellbore production distributions that are lognormal rather than normal.

The Figure 1 well production frequency distributions are computed for a 2D reservoir section with standard normal distribution of reservoir porosity ranging over 0.1 < φ < 0.3 with mean φ ~ 0.2. The poro-connectivity parameter values 0 < α < 30 control the degree of flow connectivity and hence flow channelisation within the reservoir section. Low values of α correspond low degrees of flow channelisation, for which any given well has little chance of missing existing flow structures. The black fill histograms in the upper tier of Figure 1 show that for poro-connective parameter values 5 < α < 15, large percentages of wells produce at >70% of normalised unit peak flow. Low degrees of crustal flow channelisation given by the 5 < α < 15 poro-connectivity parameter range correspond to the standard hypothetical concept of spatially uncorrelated reservoir flow heterogeneity and flow modelling. In such flow media, wells can be drilled indiscriminately (e.g., on a regular drill grid as is commonplace in hydrocarbon fields) with statistically assured well production levels. Cyan fills show, however, that high α values penalise indiscriminate drilling.

Figure 1 Six well productivity distributions for six values of empirical poro-connectivity parameter α in crustal reservoir permeability distributions κ(x,y,z) ~ exp(αφ(x,y,z)) for porosities 0.1 < φ < 0.3 with mean φ ~ 0.2. Well production for each α is normalised to unit peak value. Black/cyan histogram fills denote well production respectively above/below 70% peak productivity. For 5 < α < 15 most or all wells flow > 70% peak production. For 20 < α < 30–values observed in crustal reservoirs, however, 50% to 97% of wells flow < 70% peak production. While the impact of reduced well flow is moderate for hydrocarbon fluid production–both cyan and black wells pay–the impact on geothermal reservoir fluid production is critical–only the small number of lower tier black wells flow strongly enough to pay.

Following ambient crustal empirics (i)-(iii), actual poro-connectivity parameters are in the range 20 < α < 30 for which reservoir fluid flow is increasingly channelised. With reservoir fluid flow restricted to ever higher degrees of poro-connectivity in ever narrower flow channels, the chance that wells fail to intersect existing flow channels rises substantially. Cyan fill portions of Figure 1 lower tier histograms show that for α ~ 20 half of wells flow at productivity less than 70% of peak, with higher values 25 < α < 30 corresponding to spatially correlated flow channelling for which only 3% to 10% of wells produce greater than 70% unit peak production.

While the Figure 1 well productivity progression has only moderate impact on hydrocarbon production pay (low well productivity means that hydrocarbon pay is delayed not cancelled), the impact of crustal flow empirics on geothermal well production is critical. Only the most productive geothermal wells can turn power-producing turbines. All other such wells have little or no economic worth.

Figure 1 numerical simulations of reservoir flow structures depend on accurate flow modelling capability appropriate to flow channelisation encountered for the poro-connectivity parameter range 20 < α < 30 indicated by empirics (ii)-(iii). Solving the flow constraint equation (3) by the ‘’strong’’ methods of finite difference expression evaluations on discrete numerical meshes faces numerical stability problems when the porosity gradient term ∇φ(x,y,z) is ill-behaved and/or the poroperm empirical parameter α is large. Modelling channelised geothermal flow thus requires the use of the more mathematically-oriented ‘’weak’’ solution methods developed to handle a wide class of constitutive functions such as κ(x,y,z) ~ exp(αφ(x,y,z)) for empirical values of poro-connectivity 3 < αφ < 5 (e.g., [31,32,33,34,35]).

As no mathematical expressions exist against which to test numerical flow simulations of κ(x,y,z) ~ exp(αφ(x,y,z)) media, we base our wellbore flow validation on testing simple cases of flowing fluid through 2D numerical permeability sections κ(x,y) ~ exp(αφ(x,y)) for Figure 1 well productivity simulation statistics . We first test if the spatial distribution of fluid flow velocity computed by two weak-formulation flow solver algorithms agree, then establish that the weak-formulation flow solutions retain the spatial correlation scaling nature of the host permeability distribution established by crustal flow and microseismicity empirics (i)-(vi).

It is logical to expect that fluid flowing in the physical equivalent of our numerical simulation permeability distribution κ(x,y,z) ~ exp(αφ(x,y,z)) tracks the physical permeability structure. Flow simulations that preserve the host medium flow structure can be construed to give an ‘’accurate’’ account of wellbore coupling to the flow medium. Similarly, numerical simulation flow structures that do not preserve these properties must be presumed to be ‘’inaccurate’’.

3. Results- Wellbore Interaction with κ(x,y) ~ exp(αφ(x,y)) Permeability Distributions

The class of poroelastic flow materials under consideration is illustrated in Figure 2 crustal flow sections showing porosity spatial fluctuations and Figure 3 showing the associated permeability spatial fluctuations in accord with the crustal flow empirics (i)-(iii). From left to right, the degree of spatial correlation of porosity and permeability fluctuations increases from uncorrelated ‘’white noise’’ (well-log spectral exponent β=0), to moderately correlated ‘’pink noise’’ (exponent β=1), to strongly correlated ‘’red noise’’ (exponent β=2). Applications of two weak-formulation flow simulation methodologies to the Figure 2 and Figure 3 poroperm sections are compared in Figure 4 and Figure 5. Figure 6, Figure 7, Figure 8 and Figure 9 validate that the two flow simulation methods generate velocity flow fields that preserve the spatial correlation properties of the empirical crustal permeability distributions κ(x,y) ~ exp(αφ(x,y)) for poro-connectivity parameter α = 25.

Figure 2 Numerical realisations of poroelastic crustal sections having (left to right) increasing degrees of power-law-scaling spatial correlation. (Left) S(k) ∝ 1/k^β, β ~ 0, uncorrelated ‘’white noise’’. (Middle) S(k) ∝ 1/k^β, β ~ 1, moderately correlated ‘’1/f’’ or ‘’pink noise’’). (Right) S(k) ∝ 1/k^β, β ~ 2, highly correlated Brownian or ‘’red noise’’). Numerical wellbores across each section return the designated power-law-scaling spectral exponent β.

Figure 3 Numerical realisations of crustal permeability sections corresponding to Figure 2 porosity sections. Permeability is given by κ(x,y) ~ exp(αφ(x,y)) for α = 25 and porosity 0.1 < φ < 0.3 with mean φ ~ 0.2. Permeability range shown at right is reduced overall by factors 4 and 2 for the left and middle permeability sections, respectively.

Figure 4 Channelised numerical flow amplitude distributions for Figure 3 middle permeability distribution with constant left/right pressure difference and zero-flow on upper/lower boundaries. (Left/right) Matlab-based finite-element [31,32,33,34,35] and MRST [36,37] flow solvers.

Figure 5 Channelised numerical flow amplitude spatial distributions for Figure 3 right permeability distribution with constant left/right pressure difference and zero-flow on upper/lower boundaries. (Left/right) Matlab-based finite- element [31,32,33,34,35] and MRST [36,37] flow solvers.

Figure 6 Sequence from upper tier left to lower tier right of 2D porosity sections according to degree of well-log power-law scaling correlation S_φ(k) ∝ 1/k^β with corresponding exponent β given above each panel. 2D permeability sections for flow computation are given by κ(x,y,z) ~ exp(αφ(x,y,z)).

Figure 7 Plot groups for sequence of 2D permeability sections corresponding to Figure 6 2D porosity sections. Loglog plots within each group show the power-law trend of two-point spatial correlation function fit to power-law exponent p. The upper plot group shows the spatial correlation of permeability input to the flow Matlab finite element flow modeller [31,32,33,34,35]. The modelled flow field reproduces the host permeability field to a high degree of fidelity. Virtually identical results proceed from the MRST flow solver shown in the lower plot group [36,37].

Figure 8 Selected wellbore-centric permeability distributions of 250m diameter sections extracted from the central Figure 6 km-scale permeability section with well-log fluctuation power spectral scaling exponent, S_φ(k) ∝ 1/k^β, β ~ 1.1. The top pair, illustrating the low permeability sections, and the bottom pair, illustrating high permeability sections, summarise 169 random well-locations from the Figure 6 permeability section. Sidebars show nominal permeability scales. Figure 10 gives the lognormal wellbore production distribution for the 169 sections.

Figure 9 Production wellbore-centric inflow filamentary log(velocity) distributions for Figure 8 permeability sections. Logarithm of nominal inflow range shown in side-bars. As seen in Figure 10, the distribution of wellbore-centric inflow amplitudes is lognormal.

Figure 10 Distribution of normalised wellbore inflow amplitudes for a km-square permeability section sampled by 169 randomly located 250m-radius flow sections as illustrated in Figure 8. The central section of Figure 6 shows the km-square source permeability section with well-log spectral power scaling, S_φ(k) ∝ 1/k^β, exponent β ~ 1.1. Red curve is the lognormal distribution fit to the histogram of wellbore flow amplitudes as illustrated in Figure 9.

Figure 2 and Figure 3 poroperm sections are 1024x1024 node numerical grids representing porosity in the range 0.1 (cool tints) to 0.3 (warm tints) and the corresponding permeability given by κ(x,y) ~ exp(αφ(x,y)) for α = 25. The middle panel β = 1 ‘’pink noise’’ poroperm sections represent most actual ambient crustal reservoir sections. Flow model computation applied to these sections should show that Darcy flow velocities increase/decrease in areas of high/low poroperm properties, respectively. The Figure 4 and Figure 5 flow velocity distributions do not, however, reproduce the poroperm distributions because fluid flow measures poroperm connectivity rather than poroperm values. The degree of poroperm connectivity given by fluid flow solvers is controlled by the parameter α, thus giving the degree of flow channelisation. The flow continuity measure of spatially correlated poroperm properties given by the flow model solvers marks a fundamental departure from standard reservoir modelling that is based on uncorrelated poroperm properties. Standard modelling routinely assumes small values of α, while actual crustal flow structures have high degrees of channelisation that systematically go unrecognised in standard reservoir modelling. Accurate flow continuity modelling involving channelled flow is paramount for incorporating microseismicity emission imaging data into structural profiles of convective geothermal flow systems [7].

Figure 4 and Figure 5 show the Darcy velocity amplitude fields generated from Eq. (1) evaluated for the weak-formulation pressure field solutions of constraint Eq. (3) using two Matlab-based solver implementations [31,32,36,37]. These flow sections are applied respectively to the Figure 3 middle and right-hand spatially correlated permeability distributions. The flow model boundary conditions are constant differential pressure between the left/right vertical section boundaries for zero fluid flow condition across the upper/lower section boundaries.

Flow channelisation effects introduced by the spatially correlated permeability fields are the outstanding feature of Figure 4 and Figure 5, with the lesser/greater spatial correlation of Figure 3 middle/right creating the lesser/greater the channelisation of Figure 4 and Figure 5, respectively. Fluids entering at the left sides of the flow sections pass across to the right side of section via spatially erratic channelised poro-connectivity pathways. Flow velocity amplitudes denoted by the colour-bar vary by an order of magnitude. It is immediately evident from Figure 4 and Figure 5 that ideally positioned production wellbores will intersect high velocity flow structures (warm tints) and avoid low velocity flow structures (cool tints). The well productivity distributions associated with the high and low flow velocities of Figure 4 give the productivity frequency distributions appearing in Figure 1 for the α = 25 value of poro-connectivity parameter.

We note in Figure 4 and Figure 5 that the two flow solvers agree in their response to the gross scale distribution of permeability of Figure 3 sections. Where there are high permeabilities in the Figure 3 distributions, there are corresponding high flow velocities in Figure 4 and Figure 5, and equivalently for low permeabilities and flow velocities. The two solvers differ in the spatial discrimination of the flow response to the permeability distributions. The Matlab-based finite-element solver [31,32,33,34,35] maintains a higher degree of contrast between high and low permeability domains than does the Matlab-based MRST solver [36,37]. The difference in velocity amplitude resolution is due to how the two solvers distribute the κ(x,y,z) ~ exp(αφ(x,y,z)) permeability heterogeneity between the computational mesh cells. The Matlab finite element solver evaluates the trial solutions globally between all mesh cells, while the MRST solver evaluates the trial solutions locally between nearest neighbour mesh cells. In principle, the finite element solver of Figure 4 and Figure 5 (left) is more accurate than the MRST solver of Figure 4 and Figure 5 (right), but the difference is likely to be less significant than other sources of modelling and observational error.

Figure 6 and Figure 7 assess the degree to which the two “weak” flow structure solver methodologies preserve the spatial correlation nature of ambient crust poroperm structures. From the crustal flow observational empiric (v), the two-point spatial correlation properties of crustal section permeability have the form Γ(r) ~ 1/r^p, where the power-law exponent can in principle range from 0 < p < 1. The observed two-point spatial correlation power-law exponent p ~ ½ of (v) applies to volumetric distributions of permeability and microseismic events where the many-plane superposition of 2D distributions lowers the effective power-law exponent range. In the single-plane numerical simulation cases given in Figure 6 for 2D permeability sections, the exponent parameter range is 0 < p < 1.

The full 0 < p < 1 range of controlling numerical porosity spatial distributions in Figure 6 begins with the uncorrelated porosity S_φ(k) ∝ 1/k^β, β ~ 0, in the upper left and proceeds to the maximally correlated porosity S_φ(k) ∝ 1/k^β, β ~ 2 at lower right. The plot arrays of Figure 7 give the corresponding sequence of two-point 2D section spatial correlation power-law fits Γ(r) ~ 1/r^p for permeability event-pair offsets r for the two solver methods.

The upper set of Figure 7 plots shows that the two-point section power-law spatial correlation distributions for the originating permeability field of Figure 6 are matched in detail by permeability fields provided by the Matlab finite element solver through reverse expressing Darcy’s law as κ(x,y) = |v(x,y)|/|∇P(x,y)|. As the input permeability field is essentially duplicated by the reverse permeability field, we are assured that the Matlab finite element solver is numerically robust for κ(x,y,z) ~ exp(αφ(x,y,z)) crustal permeability fields for reasonable values of poro-connectivity parameter α. Virtually identical results are provided by the MRST solver shown in the lower set of Figure 7 plots.

Two-point correlation power-law exponents p ~ 0 occur for spatially uncorrelated permeability derived from spatially uncorrelated porosity illustrated in Figure 6 top tier. There is no preferred pair-wise offset range for spatially uncorrelated distributions. As the controlling porosity spatial correlation increases towards the β ~ 1 porosity spatial correlation condition characteristic of crustal flow systems in the Figure 6 central tier, the pair-wise spatial correlation becomes optimally organised for fluid flow by favouring nearest-neighbour ranges before falling steeply for increasing pair-range with exponents p ~ 1. Intermediate values of two-point correlation exponent p ~ ½ occur for the high-correlation spectral distributions S_φ(k) ∝ 1/k^β, β ~ 2 given in the Figure 6 lower tier.

Review of the Figure 7 sequence of two-point spatial correlation distributions for the Figure 6 sequence of porosity sections shows that well-log fluctuation power spectral scaling exponent β ~ 1 observed in the ambient crust is the most efficient value for generating areas of high and low porosity given by warm/cool colours, respectively. As a general consequence of the reservoir flow structure empirical poroperm relation κ(x,y,z) ~ exp(αφ(x,y,z)), wellbores drilled into high porosity regions are systematically more permeable over wider crustal volumes and are therefore more productive than wellbores drilled into low porosity regions. This phenomenon is a mechanical picture of well productivity lognormal distribution of crustal flow empiric (iii) as seen in Figure 10.

From the ambient crust poroperm empirics (i)-(vi), we see how entrenched is the key drilling statistic observed world-wide in the form of lognormal well-productivity statistics for all crustal fluids. What we do not, however, at present know is how to proactively deal with the entrenched and inherent statistical uncertainty caused by flow channelling arising from universal spatially correlated poroperm randomness at all scales. As drilling outcome predictions based on the law of averages are poor in spatially correlated poroperm media, we can directly see from Figure 4 and Figure 5 channelised flow the need to acquire and quantify information about the large-scale reservoir flow structure if we are to reduce the overall cost of production well drilling in convective geothermal reservoirs [7].

Building on Figure 4 and Figure 5 channelised flow computations, we now return to the source of Figure 1 statistics arising from large-scale flow channelling structures in crustal reservoirs. Figure 8 and Figure 9 show 2D permeability section flow simulation data that visually quantify the drilling outcome problem posed by spatially correlated poroperm media. In composing Figure 8, we start with a km-scale poroperm distribution such as illustrated in Figure 6 for porosity spatial correlation well-log spectral scaling exponents in range 1.1 < β < 1.3 as given in flow empiric (i). Within the Figure 6 km-scale permeability section, we can compute wellbore fluid inflow for a sequence of production wellbore locations that are taken to be at the centre of 250m diameter drainage zones. Figure 8 compares two permeability distributions that are the least productive (above) with two of the most productive 250m diameter flow sections (below). Figure 9 illustrates the wellbore-centric filamentary distributions of log(flow) amplitudes for the Figure 8 permeability sections. The accompanying side-bar nominal amplitude scales show the flow sections range over a factor of 6 in permeability and factor 20 ~ e³ in wellbore flow. Figure 10 shows that the distribution of wellbore flows of the type illustrated in Figs 8-9 is lognormal, as observed for crustal reservoir flow systems worldwide as given in flow empiric (iii). In reference to Figure 1, Figure 10 shows that for a given permeability section (or volume) only the right-most handful of geothermal wells are productive in terms of spinning turbines. The need to acquire flow structure information in advance of drilling convective geothermal production wells is manifest.

4. Discussion - Wellbore Flow Evaluation of 3D Convective Geothermal Flow Structures

While Figures 1-10 are based on flow model computations for 2D sections embodying the spatial correlation properties of crustal flow systems, applying these spatial-correlation outcomes to actual reservoirs requires 3D computation. Via the essential universality of ambient crustal fluid flow empirics (i)-(vi), we can expand from 2D generic reservoir section flow simulations of Figs 1-10 intended for geothermal flow systems to 3D applications using recent field-scale fluid flow data acquired in crystalline basement rock [14] and hydrocarbon-bearing shale formations [7].

In Figure 11 we illustrate the Matlab-implemented finite-element flow solver applied to 3D wellbore-centric flow modelling of heat advective flow transport at 1.5km depth in basement rock [14]. To match the steady-state temperature profile observed in the wellbore, inflow to the well was conjectured to occur at the observed temperature anomaly depths. The inflow structure was modelled as a planar zone of enhanced permeability as might be attributed to a standard ambient fracture/fault structure. Flow boundary conditions were adjusted until the model wellbore temperature profile matched the observed wellbore flow profile. The 3D finite-element flow modelling procedure provided good model matches for multiple wellbore temperature profiles, indicating that 3D fluid flow modelling in spatially complex media is feasible.

Figure 11 Quarter-section of 3D version of 2D wellbore-centric inflow modelling of crustal flow systems κ(x,y) ~ exp(αφ(x,y)) illustrated in Figures 1-10. The pictured volumetric flow distribution is dominated by a horizontal planar flow structure which transports heat to the wellbore fluid [14].

In a second 3D field-scale flow image, Figure 12 displays surface seismic emission acquisition/processing image data obtained during shale formation stimulation and production of a 1km-scale crustal volume [7]. The three-panel time sequence of flow stimulation of the shale formation shows tracking of seismic emission energy associated with fluid flow (white/black smudges). Fluid from a lateral wellbore (black line) enters the formation via the imaged fracture system (red lines) at the point of wellbore intersection (left panel), then spreads (mid/right panels) through the fracture-connectivity pathways to flow-stimulate seismic emissions elsewhere in the fracture structure of the shale formation. Crustal fluid flow empirics (i)-(vi) show that such seismic emission flow imaging data acquired for convective geothermal flow systems can provide flow-model structures within which to embed precision flow images V to give production well output Q ~ ρC x T x V. Production well drilling guided by combining enhanced-flow seismic emission data and precision wellbore-centric flow modelling generated from wellbore flow data can greatly increase the efficiency of convective geothermal system drilling.

Figure 12 Time sequence of fracture-system enhanced flow structure from hydrocarbon-bearing shale formation as depicted by surface seismic sensor array data acquisition/processing [7]. Red lines represent fracture-flow system. Black/white smudges represent seismic emissions generated by stimulation fluids injected into the crust from a horizontal well (black line). Seismic emissions progress as shown in the minute-by-minute sequence following the time of initial fluid injection at left. The fact that seismic emissions are confined to the red-line domain connected to the injection wellbore indicate that fluid is detected by seismic emissions associated with fluid flow along the connected red-line fracture system. This “channelled” field-scale flow image is equivalent to the Figure 4 and Figure 5 model flow images.

Figure 11 and Figure 12 stand as guide to using the Matlab-based finite element flow solver (or equivalent) to model flow into a convective geothermal production wellbore drilled into a km-scale crustal volume. A km-scale reservoir model volume with 256 nodes to a side can represent an actual geothermal field flow system at 4m spatial resolution through seismic emission data acquired by a multi-channel surface sensor array as outlined in [7]. The reservoir model comprises a generic 3D crustal empiric permeability number field κ(x,y,z) ~ exp(αφ(x,y,z)) that is embedded with a 3D enhanced flow structure image derived from field seismic emission data.

5. Conclusions

In convective geothermal reservoirs, channelised flow arising from highly attested spatially correlated crustal poroperm properties leads to lognormal distributions of well productivity. While hydrocarbon production can extract pay from both low and high productivity wells, the persistent high risk of low productivity non-pay geothermal wells greatly hinders investment in convective geothermal power projects [7]. Geothermal drilling risk can be addressed by the growing ability of multi-channel seismic emission data acquisition/processing to locate high-flow geothermal drilling targets and avoid low-flow areas.

Crustal flow channelisation empirics at scales from cm to km is given in terms of the spatial correlation properties of crustal permeability distribution κ(x,y,z) ~ exp(αφ(x,y,z)) for large values of the poro-connectivity parameter α. Flow modelling for geothermal production wells must handle the spatially-erratic high-amplitude spatial fluctuations arising from the κ(x,y,z) ~ exp(αφ(x,y,z)) empiric. Two ‘’weak formulation’’ flow solver algorithms implemented in Matlab demonstrate that these numerical methods give stable accurate reservoir flow modelling capability for the realistic range of poro-connectivity parameter α. We project that large scale seismic emission flow imaging data embedded in the crustal empiric permeability structure κ(x,y,z) ~ exp(αφ(x,y,z)) will allow precision wellbore flow modelling to refine reservoir flow structures and quantitatively model production wellbore flow for sustainable and progressive resource operation of convective flow systems.

Author Contributions

PCL performed the computations; PCL and PEM jointly prepared the manuscript.

Competing Interests

The authors have declared that no competing interests exist.