Crustal Reservoir Flow Simulation for LongRange SpatiallyCorrelated Random Poroperm Media
Peter Leary ^{1,2,†,*}, Peter Malin ^{2,†}

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

Advanced Seismic Instrumentation & Research, 317 E. Millbrook, Raleigh, NC 27609, USA
† These authors contributed equally to this work.
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 LongRange SpatiallyCorrelated 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.
Abstract
In interest of accurately modelling wellbore production in convective geothermal systems, we draw attention to two approaches to numerical simulation of crustal fluid flow. The historical ‘’strong’’ approach simulates the fluid mass conservation condition on a localised finite difference mesh approximation to poroelastic continua. The more flexible ‘’weak’’ approach globally enforces fluid mass conservation by in effect integrating over all numerical mesh cells. Strong formulation crustal flow simulation emerged from numerical simulations of mathematically tractable conductive heat flow in media typified by normallydistributed/spatiallyuncorrelated thermal property distributions. Strong formulation numerical stability problems quickly intrude, however, when confronted with the lognormallydistributed/spatiallycorrelated crustal permeability distributions that pervade the ambient crust. Weak flow conservation solver formulations accurately handle abrupt wideranging crustal permeability spatial fluctuations arising from these crustal permeability structures. We use two Matlabbased weak formulation flow simulation schemes to model wellborecentric production fluid connectivity within largescale spatially erratic convective geothermal flow structures described by the empirical poroperm relation k(x,y,z) ~ exp(αφ(x,y,z)) acting over cmtokm scale lengths. Through weak formulation reservoir flow modelling of convective geothermal systems, the advancing ability to support largescale crustal flow structure models with multichannel surface seismic data acquisition/processing can guide geothermal field exploration and developmental drilling via feedback loops between crustal flow imaging data and accurate wellbore production modelling data.
Keywords
Convective geothermal flow; reservoir flow modelling; crustal permeability; finite element flow modelling
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_{0}φv_{0}ℓ. For standard values of wellbore radius r_{0} and open hole wellbore length ℓ, geothermal well performance depends directly on formation bulk inflow rate φv_{0}. Assuming adequate reservoir porosities φ(x,y,z), Darcy’s law gives wellbore power Q in terms of formation permeability structure at the wellbore, v_{0} = k(x_{0},y_{0},z_{0})/μ ∂_{r}P(r)_{r0}, with ∂_{r}P(r)_{r0} = P_{0}/r_{0}/log(R/r_{0}) 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, √((xx_{0})^{2} + (yy_{0})^{ 2} + (zz_{0})^{2}) < R, surrounding a wellbore site at location (x_{0},_{ }y_{0},z_{0}). In this mathematical form, porosity φ(x,y,z) is normally distributed between φ_{mn} < φ(x,y,z) < φ_{mx} and is spatially correlated across cmtokm length scales. The normally distributed porosity field generates an associated lognormally distributed permeability field for large empirical values of poroconnectivity 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 wellproductivity data observed for ambient crustal geofluids worldwide. In consequence of empirical values 20 < α < 30, largescale spatiallyrapid 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 nontraditional flow modelling procedures. These numerical procedures allow flow modelling to accurately connect geothermal production wellbore interval inflow velocity v_{0} at wellbore sites (x_{0},y_{0},z_{0}) to longrange spatially erratic poroconnectivity 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 largescale 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_{0 }≡ Radius of open production wellbore in meters
v_{0 }≡ Darcy fluid flow velocity at wellbore radius in meters per second
ℓ ≡ Axial length of open production wellbore in meters
(x_{0},y_{0},z_{0}) ≡ 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 Pascalseconds
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 wellcore permeability to changes in wellcore 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, 1520 < α < 2530
β ≡ Exponent of welllog powerlaw scaling Fourier powerspectra, S(k) ∝ 1/k^{β}; β ~ 0 ≡ white noise; β ~ 1 ≡ pink noise; β ~ 2 ≡ red noise.
2. Materials & MethodsComputing 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 welllog, wellcore, and wellproductivity 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:
 Welllog spatial sequences of rock properties establish the existence of longrange powerlawscaling spatial correlations as a fundamental feature of the ambient brittlefractured crust. Welllog Fourier powerspectra, 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) powerspectra, S(k) ∝ 1/k^{0} ~ const.
 Wellcore data from hydrocarbon reservoirs, geothermal formations, and deep basement rock show close spatial correlation systematics, δlogκ ~ αδφ, for spatially varying porosity φ and permeability κ [6]. On wellcore evidence, rock permeability is a physical feature of fractureconnectivity between spatiallycorrelated porosity structures at all scales. The constant α in the wellcore poroperm spatial correlation is given by the empiric condition αφ ~ 35 across two decades of mean porosity, 0.3% < φ < 30% spanning the basementtoreservoir formation spectrum [13,14]. The wellcore empiric for porosity and permeability follows from considering that connectivity between n pores scales as nfactorial, n! = n(n1)(n2)….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 welllog 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.
 Wellproductivity lognormality follows directly from the integrated poroperm relation κ ~ κ_{0} exp(αφ). The empirical condition αφ ~ 35 guarantees the lognormal wellproductivity distributions observed everywhere for groundwater aquifers, conventional and unconventional hydrocarbon reservoirs, convective geothermal systems, and mineral deposition abundances in fossil flow systems [15,16].
 Ambient crust microseismicity event magnitude distributions in the magnitude range 1 < m < 2 are lognormal rather than powerlaw/fractal, i.e., the standard GutenbergRichter 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.
 Ambient crust permeability and microseismic event separations are spatially correlated according to twopoint 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.
 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 faultlike dislocation structures [9].
Rates of wellbore heat energy extraction Q from convective system volumetric flow V =2πr_{0}φv_{0}ℓ thus depend explicitly on the details of production wellbore intersection flow v_{0} given by Darcy’s law, v_{0} = κ/μ ∂_{r}P(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 highamplitude multiscale 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} \]
Steadystate 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 poroconnectivity 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 wellknown and highly tractable Laplace equation ∇^{2}P(x,y,z) = 0 and its generalisation the Poisson equation with known source term ∇^{2}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 finitedifference 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 preemptive smoothing of poroperm spatial variability, i.e., the artificial reduction or elimination of troublesome porosity structure spatial gradients [25,26]. Such preemptive 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 finitedifference numerical flow solvers based on preemptive generation of numerically tractable flow property distributions remain the mainstay of present reservoir fluid flow simulation [27,28,29,30].
Preemptive 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 powerspectra 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 poroconnectivity 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 poroconnective 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 poroconnectivity 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 poroconnectivity 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 poroconnectivity 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 poroconnectivity 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 powerproducing 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 poroconnectivity 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 illbehaved and/or the poroperm empirical parameter α is large. Modelling channelised geothermal flow thus requires the use of the more mathematicallyoriented ‘’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 poroconnectivity 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 weakformulation flow solver algorithms agree, then establish that the weakformulation 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’’ (welllog spectral exponent β=0), to moderately correlated ‘’pink noise’’ (exponent β=1), to strongly correlated ‘’red noise’’ (exponent β=2). Applications of two weakformulation 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 poroconnectivity parameter α = 25.
Figure 2 Numerical realisations of poroelastic crustal sections having (left to right) increasing degrees of powerlawscaling 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 powerlawscaling 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 zeroflow on upper/lower boundaries. (Left/right) Matlabbased finiteelement [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 zeroflow on upper/lower boundaries. (Left/right) Matlabbased 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 welllog powerlaw 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 powerlaw trend of twopoint spatial correlation function fit to powerlaw 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 wellborecentric permeability distributions of 250m diameter sections extracted from the central Figure 6 kmscale permeability section with welllog 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 welllocations 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 wellborecentric inflow filamentary log(velocity) distributions for Figure 8 permeability sections. Logarithm of nominal inflow range shown in sidebars. As seen in Figure 10, the distribution of wellborecentric inflow amplitudes is lognormal.
Figure 10 Distribution of normalised wellbore inflow amplitudes for a kmsquare permeability section sampled by 169 randomly located 250mradius flow sections as illustrated in Figure 8. The central section of Figure 6 shows the kmsquare source permeability section with welllog 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 weakformulation pressure field solutions of constraint Eq. (3) using two Matlabbased solver implementations [31,32,36,37]. These flow sections are applied respectively to the Figure 3 middle and righthand 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 poroconnectivity pathways. Flow velocity amplitudes denoted by the colourbar 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 poroconnectivity 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 Matlabbased finiteelement solver [31,32,33,34,35] maintains a higher degree of contrast between high and low permeability domains than does the Matlabbased 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 twopoint spatial correlation properties of crustal section permeability have the form Γ(r) ~ 1/r^{p}, where the powerlaw exponent can in principle range from 0 < p < 1. The observed twopoint spatial correlation powerlaw exponent p ~ ½ of (v) applies to volumetric distributions of permeability and microseismic events where the manyplane superposition of 2D distributions lowers the effective powerlaw exponent range. In the singleplane 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 twopoint 2D section spatial correlation powerlaw fits Γ(r) ~ 1/r^{p} for permeability eventpair offsets r for the two solver methods.
The upper set of Figure 7 plots shows that the twopoint section powerlaw 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 poroconnectivity parameter α. Virtually identical results are provided by the MRST solver shown in the lower set of Figure 7 plots.
Twopoint correlation powerlaw exponents p ~ 0 occur for spatially uncorrelated permeability derived from spatially uncorrelated porosity illustrated in Figure 6 top tier. There is no preferred pairwise 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 pairwise spatial correlation becomes optimally organised for fluid flow by favouring nearestneighbour ranges before falling steeply for increasing pairrange with exponents p ~ 1. Intermediate values of twopoint correlation exponent p ~ ½ occur for the highcorrelation spectral distributions S_{φ}(k) ∝ 1/k^{β}, β ~ 2 given in the Figure 6 lower tier.
Review of the Figure 7 sequence of twopoint spatial correlation distributions for the Figure 6 sequence of porosity sections shows that welllog 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 worldwide in the form of lognormal wellproductivity 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 largescale 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 largescale 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 kmscale poroperm distribution such as illustrated in Figure 6 for porosity spatial correlation welllog spectral scaling exponents in range 1.1 < β < 1.3 as given in flow empiric (i). Within the Figure 6 kmscale 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 wellborecentric filamentary distributions of log(flow) amplitudes for the Figure 8 permeability sections. The accompanying sidebar nominal amplitude scales show the flow sections range over a factor of 6 in permeability and factor 20 ~ e^{3} in wellbore flow. Figure 10 shows that the distribution of wellbore flows of the type illustrated in Figs 89 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 rightmost 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 110 are based on flow model computations for 2D sections embodying the spatial correlation properties of crustal flow systems, applying these spatialcorrelation 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 110 intended for geothermal flow systems to 3D applications using recent fieldscale fluid flow data acquired in crystalline basement rock [14] and hydrocarbonbearing shale formations [7].
In Figure 11 we illustrate the Matlabimplemented finiteelement flow solver applied to 3D wellborecentric flow modelling of heat advective flow transport at 1.5km depth in basement rock [14]. To match the steadystate 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 finiteelement 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 Quartersection of 3D version of 2D wellborecentric inflow modelling of crustal flow systems κ(x,y) ~ exp(αφ(x,y)) illustrated in Figures 110. 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 fieldscale flow image, Figure 12 displays surface seismic emission acquisition/processing image data obtained during shale formation stimulation and production of a 1kmscale crustal volume [7]. The threepanel 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 fractureconnectivity pathways to flowstimulate 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 flowmodel 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 enhancedflow seismic emission data and precision wellborecentric flow modelling generated from wellbore flow data can greatly increase the efficiency of convective geothermal system drilling.
Figure 12 Time sequence of fracturesystem enhanced flow structure from hydrocarbonbearing shale formation as depicted by surface seismic sensor array data acquisition/processing [7]. Red lines represent fractureflow 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 minutebyminute sequence following the time of initial fluid injection at left. The fact that seismic emissions are confined to the redline domain connected to the injection wellbore indicate that fluid is detected by seismic emissions associated with fluid flow along the connected redline fracture system. This “channelled” fieldscale 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 Matlabbased finite element flow solver (or equivalent) to model flow into a convective geothermal production wellbore drilled into a kmscale crustal volume. A kmscale 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 multichannel 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 nonpay geothermal wells greatly hinders investment in convective geothermal power projects [7]. Geothermal drilling risk can be addressed by the growing ability of multichannel seismic emission data acquisition/processing to locate highflow geothermal drilling targets and avoid lowflow 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 poroconnectivity parameter α. Flow modelling for geothermal production wells must handle the spatiallyerratic highamplitude 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 poroconnectivity 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.
References
 Clauser C, Huenges E. Thermal conductivity of rocks and minerals. Rock physics and phase relations: A handbook of physical constants. Washington: American Geophysical Union; 1995.
 Carslaw HS, Jaeger JC. Conduction of heat in solids. Oxford: Clarendon Press; 1959. p. 510.
 Theis CV. The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage. Eos Trans Am Geophys Union. 1959; 16: 519524. [CrossRef]
 Muskat M. The flow of fluids through porous media. J Appl Phys. 1937; 8: 274282. [CrossRef]
 Kruger P, Ramey HJ. Workshop on geothermal reservoir engineering. Stanford: Stanford University Workshop; 1975.
 Leary P, Malin P, Saarno T, Kukkonen I. αφ ~ αφ_{crit}Basement rock EGS as extension of reservoir rock flow processes. Proceedings of 43th Workshop on Geothermal Reservoir Engineering; 2018 February 1214; Stanford, CA, USA. Stanford: Stanford University.
 Leary P, Malin P, Saunders G, Sicking C. Seismic imaging of convective geothermal flow systems to increase well productivity. J Energ Power Technol. 2020; 2: 012. [CrossRef]
 Leary P, Malin P, Saarno T, Heikkinen P. St1 deep heat egs projectcomputation of wellboretowellbore permeability stimulation. Proceedings of World Geothermal Congress; 2021 May 2126; Reykjavik, Iceland. Bonn: International Geothermal Association.
 Leary P, Malin P. MEQ ~ permeabilitymodelling high frequency emissions from stimulationinduced seismic activity in the ambient crust. Proceedings of World Geothermal Congress; 2021 May 2126; Reykjavik, Iceland. Bonn: International Geothermal Association.
 Leary PC. Rock as a criticalpoint system and the inherent implausibility of reliable earthquake prediction. Geophys J Int. 1997; 131: 451466. [CrossRef]
 Leary PC. Fractures and physical heterogeneity in crustal rock, in heterogeneity in the crust and upper mantle: Nature, scaling, and seismic properties. New York: Kluwer Academic/Plenum Publishers; 2002. p. 155p. 186. [CrossRef]
 Leary P, Malin P, Niemi R. Fluid flow & heat transport computation for powerlaw scaling poroperm media. Geofluids. 2017; 2017: 9687325. [CrossRef]
 Leary P, Malin P, Pogacnik J, Rugis J, Valles B, Geiser P. Lognormality, dk~ kdφ, EGS, and all that. Proceedings of 39th Stanford Geothermal Workshop, Stanford University; 2014 February 2426; Stanford, CA, USA. Stanford: Stanford University.
 Leary P, Malin P, Saarno T, Kukkonen I. Prospects for assessing enhanced geothermal system (EGS) basement rock flow stimulation by wellbore temperature data. Energies. 2017; 10: 1979. [CrossRef]
 Success of geothermal wells: A global study, international finance corporation [Internet]. Washington: International Finance Corporation; 2013. Available from: https://www.ifc.org/wps/wcm/connect/topics_ext_content/ifc_external_corporate_site/sustainabilityatifc/publications/publications_gpn_geothermalwells.
 Malin P, Leary P, Shalev E, Rugis J, Valles B, Boese C, et al. Flow lognormality and spatial correlation in crustal reservoirs: IIWheretodrill guidance via acoustic/seismic imaging. Proceedings of the World Geothermal Congress; 2015 April 1924; Melbourne, Australia. Bonn: International Geothermal Association.
 Leary P, Malin P, Saarno T. A physical basis for the GutenbergRichter fractal scaling. Proceedings of 45th Workshop on Geothermal Reservoir Engineering; 2020 February 1012; Stanford, CA, USA. Stanford: Stanford University.
 Zaky DA, Nugraha AD, Sule R, Jousset P. Spatial analysis of earthquake frequencymagnitude distribution at geothermal region in the south of Bandung, West Java, Indonesia. Potsdam: 9th International Workshop on Statistical Seismology; 2015.
 Leary P, Malin P, Saarno T, Heikkinen P, Diningrat W. Coupling crustal seismicity to crustal permeability  Powerlaw spatial correlation for EGSinduced and hydrothermal seismicity. Proceedings of the 44th Workshop on Geothermal Reservoir Engineering; 2019 February 1113; Stanford, CA, USA. Stanford: Stanford University.
 Leary P, Malin P. Correlation function Γ_{meq}(r) ~ 1/r^{1/2} coupling of microseismicity to permeability  The basis for fluid flow seismic image targeting for geothermal production wells. Reykjavik: Proceedings World Geothermal Congress; 2020.
 Bear J. Dynamics of fluids in porous media, New York: American Elsevier; 1972.
 Warren JE, Price HS. Flow in heterogeneous porous media. Soc Pet Eng J. 1961; 1: 153169. [CrossRef]
 Warren JE, Skiba FF. Macroscopic dispersion. Soc Pet Eng J. 1964; 4: 215230. [CrossRef]
 Warren JE, Root PJ. The behavior of naturally fractured reservoirs. Soc Pet Eng J. 1963; 3: 245255. [CrossRef]
 Hubbert MK. Motion of ground water. Trans N Y Acad Sci. 1941; 3: 3955. [CrossRef]
 Hubbert MK. Darcy’s law and the field equations of the flow of underground fluids. Int Assoc Sci Hydrol. Bull. 1957; 2: 2359. [CrossRef]
 Pruess K, Oldenburg C, Moridis G. TOUGH2 User’s Guide, Version 2.1 [Internet]. Berkeley: Lawrence Berkeley National Laboratory; 1999. Avilable from: https://www.osti.gov/biblio/751729tough2userguideversion. [CrossRef]
 Pruess K. The TOUGH codesA Family of simulation tools for multiphase flow and transport processes in permeable media. Vadose Zone J. 2004; 3: 738746. [CrossRef]
 Peaceman DW. Interpretation of wellblock pressures in numerical reservoir simulation. Soc Pet Eng J. 1978; 18: 183194. [CrossRef]
 Thomas K, Dixon TN, Pierson RG. Fractured reservoir simulation. Soc Pet Eng J. 1983; 23: 4254. [CrossRef]
 Partial differential equation toolbox user's guide [Internet]. Natick: The MathWorks Inc. Available from: http://www.math.chalmers.se/~mohammad/teaching/PDE1_TM/pde1_lab/lab2/toolboxguide.pdf.
 Hughes TJ. The finite element method: Linear static and dynamic finite element analysis. New York: Dover Publications; 2020.
 Leary P, Malin P, Niemi R. Finiteelement modelling of wellboreobserved fractureborne heat advection  application to EGS stimulation in basement rock. Proceedings of 41st Workshop Geothermal Reservoir Engineering; 2016 February 2224; Stanford, CA, USA. Stanford: Stanford University.
 Brehme M, Regenspurg S, Leary P, Bulut F, Milsch H, Petrauskas S, et al. Injectiontriggered occlusion of flow pathways in geothermal operations. Geofluids. 2018; 2018: 4694829. [CrossRef]
 Brehme M, Leary P, Milsch H, Petrauskas S, Valickas R, Kamah Y, et al. Natural and altered physical flow structures in the earth´s crust with applications for geothermal energy. Proceedings of 43rd Workshop on Geothermal Reservoir Engineering; 2018 February 1214; Stanford, CA, USA. Stanford: Stanford University. [CrossRef]
 Aarnes J, Gimse T, Lie KA. An introduction to the numerics of flow in porous media using Matlab. In Geometric Modelling, Numerical Simulation, and Optimization. Berlin: Springer; 2006.
 Lie KA. User guide for the MATLAB Reservoir Simulation Toolbox (MRST). Cambridge: Cambridge University Press; 2019.