Seismic Imaging of Convective Geothermal Flow Systems to Increase Well Productivity
Peter Leary ^{1, 2, †, *}, Peter Malin ^{2, †}, Graeme Saunders ^{1, †} and Charles Sicking ^{3, †}

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

Advanced Seismic Instrumentation & Research, 1311 Waterside, Dallas, TX 75218, USA

Ambient Reservoir Monitoring, Peninsula Place, Houston, TX 77459, USA
† These authors contributed equally to this work.
Academic Editor: Catalin Teodoríu
Special Issue: Advancement of Geothermal Technology for Sustainable Energy Production
Received: December 03, 2019  Accepted: July 26, 2020  Published: July 29, 2020
Journal of Energy and Power Technology 2020, Volume 2, Issue 3, doi:10.21926/jept.2003012
Recommended citation: Leary P, Malin P, Saunders G and Sicking C. Seismic Imaging of Convective Geothermal Flow Systems to Increase Well Productivity. Journal of Energy and Power Technology 2020;2(3):28; doi:10.21926/jept.2003012.
© 2020 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
A core feature of convective geothermal resource production is wellbore energy flow Q ~ ρC x T x V. E.g., for wellbore fluid of volume heat capacity ρC ~ 4.3MJ/m^{3}∙^{o}C, temperature T ~ 230^{o}C, and volumetric flow V ~ 50L/s, wellbore heat energy production is Q~ 50MW_{th} ~ 5MW_{e}. Wellbore fluid flow V =2πr_{0}φv_{0}ℓ for open wellbore length ℓ is given in turn by the spatially variable product crustal porosity times crustal fluid velocity v ≡ φv_{0} at the wellbore radius r_{0}. For a geothermal wellbore to be productive (nominal Q ~ 5MW_{e}), locally variable bulk inflow rates v = φv_{0} across crustal volumes of dimension ℓ must be adequate to sustain high wellbore flows (nominal V ~ 50L/s). Wideranging crustal well productivity statistics show that few crustal wells flow at these rates. This is not surprising as local bulk flow ~ 10^{2} m/s needed for production wellbores is decades greater than ambient bulk fluid flow ~ 10^{8}10^{7}m/s characteristic of natural convective geothermal systems. Such rare high flow locales must be found. While existing crustal surveys generally fix resource temperatures T with nominal 500m spatial resolution, as yet no survey technique provides data to locate productiongrade flow rates at nominal ℓ ~ 50m spatial resolution. In consequence, many costly unproductive wells are drilled before sites of productive local geothermal fluid bulk flow v = φv_{0} ~ 10^{2} m/s are located. A seismic survey methodology sensitive to crustal flow v = φv_{0} at ℓ ~25m resolution has evolved from multichannel seismic reflection technology applied to the production of shale hydrocarbons. The fieldproven seismic flowimaging methodology can be adapted from sedimentary terrains to volcanic terrains through appropriate seismic refraction means for generating effective seismic velocity models for convective geothermal flow volumes. Numerical simulation for characteristic ambient crustal heterogeneous flow distributions v = φv_{0} at ℓ ~50m resolution shows that traveltime data for seismic source energy refracted from deep geothermal wellbores to surface seismic sensor arrays can replicate the multichannel seismic flowimaging capability demonstrated for shale formations. Multichannel seismic reservoir flow monitor detection and mapping of flowinduced seismic emissions in shale formations can be adapted to achieve similar mapping of convective geothermal flow system structures with v = φv_{0} ~ 10^{2}m/s at nominal ℓ ~ 50m resolution. Such seismic emission flow structure maps can greatly reduce the uncertainty and hence the cost of drilling effective production wells at brownfield sites and assessment/development wells at greenfield sites
Keywords
Geothermal energy; crustal fractures; crustal fluid flow; microseismicity; reservoir seismic monitoring; geocriticality; multichannel seismic data processing; tomographic fracture imaging; reduced drilling costs
1. Introduction
To bring greater certainty to drilling convective geothermal systems, we adopt a new perspective on microseismicity through which to map the inherent and often strong crustal permeability spatial heterogeneity that creates drilling uncertainty. The new perspective can be expressed as “Meqs ~ Permeability”, with the sense that in the ambient crust the spatial/magnitude distributions of microseismicity closely follow the spatial/magnitude distributions of permeability. An observational procedure using surface seismic array data to locate microseismicity associated with permeability can identify high permeability flow sites with sufficient spatial resolution to target drilling in convective flow systems. Creating drilling certainty by highresolution flow structure mapping can eliminate many of the 2030 low productivity “exploration” wells now typically required for geothermal resource development. In addition to reducing direct drilling costs, mastering well production uncertainty allows more rapid and reliable costing of appropriate power generation facilities at brownfield sites, and lowers the cost threshold for developing smaller resource plays at greenfield sites.
Our new microseismicity perspective unfolds through reinterpreting the traditional GutenbergRichter (GR) frequencysize relation for microseismicity in the ambient crust. We first note that observed ambient crust microseismicity has a lognormal frequency distribution rather than the traditional powerlaw or fractal frequency distribution [1]. We then note that the newly observed microseismicity lognormal spatial and size distributions are statistically congruent with the lognormal spatial and size distributions of crustal permeability derived from welllog, wellcore and wellproductivity data across a wide range of geological settings [2,3,4]. The observed statistical congruence of ambient crust microseismicity with crustal permeability spans five to six orders of microseismicity scale magnitude m_{eq} ~ 3 to m_{eq} ~ +3. Where once narrow ranges of individually identifiable microseismic events have been interpreted to imply isolated faults as fluid conduits [5], our new perspective focuses the evidence for systematic microseismicitypermeability interactions across mmkm range of permeability scale sizes in the ambient crust [6,7,8].
Systematic interaction between crustal microseismicity and crustal permeability has been recently observed in hydrocarbonbearing shale formations [9,10,11,12,13,14]. Surface seismic array data acquired during stimulation of shale formations can be routinely processed into detailed maps of crustal flow structures at 25m spatial resolution. Proceeding from our perspective on ambient crustal microseismicity, we seek here to adapt the sedimentary terrain flowimaging methodology to volcanic terrains hosting convective geothermal flow systems [15]. While the seismic wave transmission properties of sedimentary and volcanic terrains vary significantly [16,17], welllog, wellcore and wellproductivity data from the two terrains show similar relations between microseismicity and permeability [1,2,3,4,6]. Our stated technology adaptation task thus focuses overcoming the difference in seismic wave transmission in orderly sedimentary sequences versus disorderly volcanic intrusion/extrusion sequences, with the goal of effectively mapping convective geothermal flow system heterogeneity.
In our discussion, §2 seeks first to be precise about the multiscale nature of fluid flow in the ambient crust. The longheld engineering view based on ignoring heterogeneity by averaging over crustal property fluctuations proves to be deeply incompatible with the empirics of ambient crustal fluid properties recorded by welllogs, wellcores, and wellproduction data. In particular, the engineering view cannot explain the observed statistics of ambient crust microseismicity. By contrast, our empirical view makes direct quantitative connection between observed crustal microseismicity and observed crustal permeability at all scales. With fully incorporated ambient crustal flow empirics in our discussion, it straightforward for §2 to quantify how detecting microseismicity associated with crustal flow at all scales can generate drilling cost savings for production wells at brownfield sites and exploration wells at greenfield sites. In §3 we outline the systematics of acquiring/stacking multichannel seismic data to detect lowlevel reflection imaging signals in sedimentary sections, then show via numerical computation how to achieve the parallel signal acquisition and processing systematics for volcanic terrains. With affirmative §3 results in hand, §4 cites standard drilling costs to quantify how lowcost remote seismic sensing can eliminate a great deal of uncertainty in drilling highcost geothermal wells in order to mobilise investor funds for winning clean baseload energy from convective geothermal sites worldwide.
2. Materials and Methods
Crustal rock with its everpresent fluid component is very far from a traditional engineering material. Engineering materials are typically valued for strong degrees of spatial uniformity and predictability, neither of which are conspicuous features of crustal flow systems [1,2,3,4,6]. Expressing crustal flow heterogeneity in terms of engineering materials obscures our understanding of the physical origin and nature of the heterogeneity and hinders development of operational means for managing crustal flow heterogeneity. Our new perspective on crustal microseismicity abandons the engineering sense that large isolated microseismic events dominate crustal flow properties, and instead takes microseismicity emissions to be relevant to crustal flow on all scales. Large seismic arrays open to the full range of microseismicity data offer geothermal operators the means of effectively surveying the inherent crustal flow heterogeneity of their resource.
2.1 Microseismicity as Signifier of the FluidRock Interaction Empirics of Crustal Permeability
The conventional view of crustal fluid flow treats rock as an elastic continuum punctuated by small scale porosity that passively hosts the crustal fluid [18,19]. Darcy’s mid19^{th}century groundwater flow empirical relation viewed crustal permeability as analogous to the permeability of unconsolidated sands used to filter municipal water supplies [20,21]. The unconsolidated sands analogy for crustal permeability was maintained during the early 20^{th}century groundwater survey support of expanding western US agriculture [22,23], and passed on to 1930s groundwater and hydrocarbon reservoir flow modelling based on the mathematics of heat conduction in thermally uniform solids [24,25,26,27].
For the centenary of Darcy’s law, Hubbert produced the Figure 1 sketch in making a case for ignoring smallscale groundwater flow complexities [28,29]. The historical perception of crustal permeability is here made explicit: any complication of groundwater flow at small scales (left portion of sketch) in a geological unit inevitably averages out to a quasiuniform “effective” flow property characteristic of that unit; only at the formation scale (right portion of sketch) do flow properties undergo significant changes that cannot be averaged over. Hubbert’s Figure 1 sketch was later formalised into the familiar concept of the representative elementary volume (REV) as an aid to numerical modelling of reservoir flow [30,31].
The underlying but largely unacknowledged reason for Figure 1 and the REV is the desire that smallscale and sparsely sampled welllog and wellcore data be accepted as properly, or at least “effectively”, representative of the groundwater aquifer or hydrocarbon reservoir formation flow properties at the reservoir scale. In effective medium theory, small scale welllog and wellcore flow property sampling can be “upscaled” for reservoirscale flow modelling [32]. Despite considerable early cautionary evidence for the “flaw of averages” [33,34,35,36], Figure 1 spatial averaging paradigm remains the default position for dealing with crustal flow heterogeneity.
The Figure 1 sketch implicitly assumes that spatial averages are effectively bounded at all significant scales. Such bounded behaviour for random spatial fluctuation sequences occurs only if the spatial fluctuations are uncorrelated in the appropriate scale range (the central limit theorem [6]). Uncorrelated random fluctuations have wellknown statistical and spectral conditions. For a random fluctuation sequence of physical values χ such as porosity at a sequence of physical locations x, χ(x), to be spatially uncorrelated (and therefore suitably bounded), the Npoint autocorrelation function A(y) ≡ 1/N ∑_{i}_{=1…N} χ(xi)χ(xi + y) must obey the twin conditions that A(y=0) = 1 and A(y≠0) = 0 [6]. The associated spectral condition is that the Fourier powerspectrum of χ(x) is independent of spatial frequency k, Χ(k)^{2} ~ const in the appropriate scale range kmin < k < kmax.
Figure 1 Sketch of Hubbert spatial averaging process for characterising crustal properties related to fluid flow in geological formations [28]. At granular scales (left) and formation scales (right), spatial fluctuations shown along the vertical axis can be arbitrarily large over a given range of crustal volumes ∆V shown along horizonal axis. At reservoir scales (center volume range), spatial fluctuations average out to an equivalent flow medium with stationary mean value (dashed line). Accordingly, it is tacitly assumed that sparse smallscale sampling of reservoirscale variations provides a statistically robust measure of formation flow properties.
With the physical and statistical conditions lying behind Figure 1 and the REV made explicit, we can see unambiguously that no such condition is observed for the actual crust. Welllog rock physical property sequences in geological settings worldwide systematically violate the Χ(k)^{2} ~ const spectral condition required for Figure 1. Instead, welllog Fourier spectra everywhere obey the powerlaw scaling function Χ(k)^{2} ~ 1/kβ, where β ~ 1.2 ± 0.1, over five decades, 1/km k < 1/cm, of scale range relevant to all aspects of geological formation and reservoir flow [2,3,7,8].
Welllog spectral empiric Χ(k)^{2} ~ 1/k shows that the Figure 1 hypothesis that crustal permeability as equivalent to an unconsolidated sand is fundamentally wrong. Whether for sedimentary, igneous or metamorphic rock, spatial averages over crustal flow property fluctuations are not bounded as hypothesised in Figure 1. The spatial correlation features of crustal rock flow properties cannot be adequately represented by equivalent or effective media approximations such as the REV or elastic continua that depend on hosting crustal fluids in spatially uncorrelated voids [6,7,8].
The Figure 1 spatial averaging hypothesis has effectively conditioned interpretations of crustal microseismicity. In a Figure 1 elastic continuum view of crustal porosity and permeability, rock failure and fracturing occur largely without reference to hosted fluids [18]. The Griffith microfracture concept of metal failure applied to rock fracture processes effectively limited the role of crustal fluids to pressure relief of confining stresses acting on rock interfaces [18,37]. In an elastic continuum framework, crustal microseismicity becomes a structureless agglomeration of frictional slips on spatially uncorrelated fracture surfaces of random size and orientation [18,38]. The longstanding empirical GutenbergRichter (GR) relation, N(m) ~ 10(b x m), became regarded as a standard organisational principle for crustal seismicity by giving the statistical expectation number N(m) of seismic slip events larger than successive event magnitudes m [1,39]. As the GR empirical constant b typically has a narrow range of values 0.75 < b < 1.5, the GR relation came to be seen as a fractal, N(L) ~ 1/LD, relating event numbers to event dimension L with fractal dimension D [1,40,41].
Figure 2 illustrates the standard GR relation applied to seismic activity in Indonesian geothermal fields [42]. The empirical GR relation parameter b is evaluated as the slope of the straightline powerlaw fit to the observed decline trend (redsquares) for high magnitude cumulative event numbers log(N(m)) at successive event magnitudes m. Of greater interest to us, however, is the raw numbers of events given in Figure 2 by the open triangles. Not only do raw event numbers steadily decline for successively larger event magnitudes, but raw event numbers steadily decline for successively smaller event magnitudes. The standard GutenbergRichter relation has no physical explanation for the event number decline at small magnitudes (beyond assuming observational deficiencies [1]). Moreover, while the cumulative data format given by the red squares manifestly has no physical meaning, the format somewhat masks the messy raw event count that conspicuously fails to match the standard GR relation at low event magnitudes.
Figure 2 Microseismicity event numbers log10(N) as a function of event magnitudes m observed at Indonesia geothermal fields [42]. Triangles denote recorded event numbers for each magnitude interval. Squares denote cumulative recorded event numbers N>m, the number of events having magnitudes ≥ m. The cumulative data format has no physical significance but produces a more stable curve with which to evaluate the purely statistical bvalue parameter. The cumulative data format disguises the event number decline for which the standard GR relation had only a loose speculative explanation in terms of large numbers of supposed events too small to be detected [1]. M_{C} in the plot indicates the magnitude at which the seismic network is “complete”, i.e., the magnitude at which seismic networks are supposed to begin to lose event count due to observational shortcomings.
The smallmagnitude event number decline in fitting the GR relation to microseismicity is historically ascribed to failure of seismic networks to record the very large number of supposed events that occur according to the standard GR powerlaw trend at low magnitudes [1,42,43,44,45,46]. In Figure 2, the quantity labelled M_{C} is the event magnitude at which it is estimated that the event catalogue is “complete”. By hypothesis, all microseismic events of magnitude m > M_{C} that occurred are recorded, while for events of magnitude m < Mc, progressively fewer events than expected for the standard GR are detected.
Figure 3 presents our alternative physically based analysis of the crustal microseismicity event size distribution [1,2,3,4]. Instead of the Figure 2 powerlaw fit restricted to the high magnitude limb of the microseismicity number distribution, Figure 3 shows that a class of lognormal distributions track the entire range of observed microseismicity event number distributions, duly accounting for event numbers for both lowmagnitude and highmagnitude events. While microseismic network observational deficiencies likely reduce the number of low magnitude events detected [1], the fidelity of the Figure 3 lognormal curve fits over the entire magnitude range indicates that incidental observational defects are not the primary cause of microseismicity event number count declines seen in Figure 2 triangle data (cf. [1] for further details).
Figure 3 Figure 2 Indonesia geothermal field microseismicity magnitudefrequency distribution (red circles) matched to a sequence of three crustal permeability distributions (blue traces). The crustal permeability frequency distribution curves are given by the permeability function κ ~ exp(αφ) evaluated for empirical parameter values α = 16:3:22 for the fixed normal porosity φ distribution shown in the lower right plot. The best fit distribution curve parameter, α ~ 19 with φ ~ 0.2, corresponds to the empirical condition αφ ~ 34 observed for crustal porositypermeability data worldwide [2,3,4]. The decline of microseismicity event numbers for low magnitudes occurs for reasons of crustal physics rather than seismic network observational deficiency [1]. Confining the standard GR powerlaw relation to high magnitude events poorly represents the underlying physics of ambient crust microseismicity [1].
The Figure 3 lognormal fit to the Figure 2 distribution of microseismicity event sizes is not an incidental result of statistical curve fitting. Rather, from Figure 2 and Figure 3, we can infer that the standard GutenbergRichter powerlaw scaling relation fails to represent a fundamental set of fluidrock interaction physical processes described in [2,3,4] that occur at all scales in the ambient crust and generate the spectral empiric Χ(k)^{2} ~ 1/k [7,8]. The observed 1/k powerlaw scaling empiric arises as crustal rock undergoes consolidation from highporosity disordered uppermost crustal sediment granularity to lowporosity ordered lowercrust ductile metamorphic granularity. Crustal fluids play an essential role in consolidation by mediating the change in poreconnectivity, i.e., permeability, of the consolidating rock. Pore connectivity is a statistical phenomenon that conditions the transition from grainscale disorder in the disaggregated uppermost crust to grainscale order in the aggregated lower crust. Disordered granularity in the disaggregated uppermost crust is characterised by uncorrelated random fluctuations with flat spectral flat powerspectra Χ(k)^{2} ~ 1/k^{0} ~ const. Ordered granularity in the aggregated lower crust is characterised by highlycorrelated random fluctuations with steeply declining spectral powerspectra Χ(k)^{2} ~ 1/k^{2} arising from the macroscopic ductile flow interfaces of lower crust metamorphics. The orderdisorder transition in crustal granularity observed in the seismic crust is an example of second order thermodynamic phase changes such as the ferromagnetic state, water at its critical point, and critical opalescence in binary fluids [3].
We can interpret the lognormal distribution of microseismicity in Figure 2Figure 3 in terms of the crustal permeability empiric κ(x,y,z) ~ exp(αφ(x,y,z)) [2,3,4,6] providing a physical basis for microseismicity event number distributions. The crustal permeability empiric derives from the highly attested wellcore poroperm relation, δφ ~ δlog(κ), for which the proportionality constant α is widely attested to obey the constraint αφ ~ 34 for mean formation porosity φ [3]. Figure 3 gives the best fit to Figure 2 microseismicity magnitude distribution for α ~ 19; for a porosity distribution with φ ~ 0.2, α ~ 19 yields αφ ~ 4. The Figure 3 crustal poroperm property αφ ~ 34 for geothermal domain microseismicity thus accords with crustal microseismicity observed in crustal domains from basement rock with mean porosities φ < 1% to reservoir rock with high mean porosities φ > 25% [3]. In contrast to the standard GR interpretation of half of microseismicity event populations, the Figure 3 fit to Figure 2 microseismicity data is essentially a zerofreeparameter prediction for the range of event magnitudes derived from fundamental crustal rockfluid interaction physics occurring at all scales [2,3,4,6,7,8].
Two important crustal factors emerge from the Figure 3 interpretation of the Figure 2 geothermal domain microseismicity magnitudefrequency distribution. First, microseismicity is inherently coupled to fluidrock properties of crustal rock rather than to purely elastic fracturemechanical properties as has been long and widely supposed [18,19]. Second, lognormal permeability/microseismicity distributions κ ~ exp(αφ) are intrinsically associated with the spatial connectivity property of crustal flow systems in general and geothermal flow systems in particular [3]. In light of these crustal empirical factors, we can posit the need to directly observe the largescale spatially erratic continuity flow structures as a routine feature of managing geothermal resource production. §2.2 expands on this feature of ambient crustal flow distributions.
2.2 Quantifying Crustal Flow Velocity Heterogeneity v = φv_{0}
Figure 2 and Figure 3 give microseismicitybased evidence that permeability spatial variations κ(x,y,z) ~ exp(αφ(x,y,z)) pervade crustal volumes at all scales to influence the distribution of ambient fluidrelated slip event magnitudes. Building on the microseismicity evidence for crustal fluidrock interaction, we can estimate the impact of crustal fluid flow velocity heterogeneity on geothermal reservoir wellbore energy production Q ~ ρC x T x V, where wellbore fluid flow V = 2πr_{0}φv_{0}ℓ is controlled by spatially variable crustal fluid flow v = φv_{0} in the vicinity of the wellbore.
The Figure 2 and Figure 3 lognormal microseismicity event magnitude distributions invalidate the Figure 1 assumption that crustal spatial flow properties exhibit the spatial fluctuation spectral condition Χ(k)^{2} ~ 1/k0 ~ const. Crustal flow simulation in Figure 4 illustrates how spatial flow connectivity structures controlled by the Fourier spectral powerlaw scaling function Χ(k)^{2} ~ 1/kβ, β ~ 1.2 ± 0 affect crustal permeability distribution κ(x,y,z) ~ exp(αφ(x,y,z)). Welllog data attest that this crustal spatial correlation property holds across the five decades of scale range 1/km < k < 1/cm relevant to geothermal resource production [7,8]. The overlaid red and blue synthetic 2D flow distributions in Figure 4 are controlled by the spatially correlated random porosity fluctuations that are inherent material features of the ambient crust.
Figure 4 Overlay of two fluid velocity fields v = φv_{0} (red, blue arrows) for crustal permeability sections generated by two random number realisations of spatially correlated random porosity fields with Χ(k)^{2} ~ 1/k spectral power scaling [7,8]. While such velocity fields have the same population of high and lowflow areas, the spatial organisation of the velocity fields is subject to essentially infinite variety. Drilling in the vicinity of highflow clustering is more productive than drilling in the vicinity of lowflow clustering. No sparse smallscale sampling algorithm suffices to represent the illustrated largescale flow variations significant to convective geothermal wellbore production. Instead, to locate the highflow clusters for optimal drilling, crustal fluid flow variations v = φv_{0} need to be mapped at the scale of interest.
It is visually evident from the Figure 4 velocity field V that some arbitrarily located wells will be good producers (i.e., have large Q), while many will be poor producers (i.e., have low Q). Figure 5 histograms show the result of statistical evaluation of Figure 4 flow field v = φv_{0} realisations. The population of velocity amplitudes is lognormally distributed, as expected from the underlying crustal permeability distribution κ(x,y,z) ~ exp(αφ(x,y,z)).
The highflow velocity structures flagged as red histogram bars in Figure 5 may be statistically correlated with high magnitude microseismicity [5], but it is not clear how strongly microseismicity is correlated with the highconnectivity crustal flow structures that constitute a convective geothermal flow system. No claims for wellbore drilling targeting are made for high magnitude microseismic event locations [5]. Further, standard microseismicity events are typically too poorly located to function as a reliable map of crustal fractureconnectivity fluid flow structures that reservoir operators can target for production well drilling [5].
Figure 5 Normalised lognormal distributions of crustal fluid velocity field v = φv_{0} amplitudes for four spatially correlated random porosity fields. Red histogram bars denote the highflow velocity structures v = φv_{0} that would be logical targets for drilling high production Q ~ ρC x T x V wellbores for V = 2πr_{0}φv_{0}ℓ. Histogram axes are normalised to maximum value of computed crustal flow velocity v = φv_{0}. There may be a statistical correlation between the high flow events and high magnitude microseismic events, but no such correlation has been established, and microseismicity event location is too poor to provide useful drilling targets [5].
Figure 6 compares the distribution of Figure 4 fluid flow velocity structure amplitudes (left) with the distribution of convective geothermal production well outputs Q ~ ρC x T x V (right) generated by producing fields worldwide [47]. Combining Figures 2, 3, 4, 5, 6, we see that the physical and statistical features of crustal flow velocity heterogeneity V are quantitatively consistent across a range of observations (microseismicity and welllog/wellcore data), scales (cm to km), and locations (geothermal fields, basement rock, aquifers, hydrocarbon reservoirs).
Figure 6 Lognormal distribution of Figure 4 synthetic fluid flow velocity amplitudes at left plotted in format of convective geothermal field lognormal wellproductivity distribution at right [47]. Velocity distribution horizonal axis (left) normalised to observed megawatts of electrical energy production MW_{e} on horizontal axis (right). Figure 4 crustal flow velocity field amplitude distributions are in quantitative agreement with the observed distribution of convective geothermal well productivity worldwide.
Figure 7 The average wellbore output Q ~ ρC x T x V MW as a function of number of wells drilled in Indonesia geothermal fields [48].
Our discussion of wellbore flow Q ~ ρC x T x V culminates with the Figure 7 exhibit of the operational consequences of the standard treatment of crustal fluid flow heterogeneity in convective geothermal resources. In absence of reliable forms of convective geothermal flow structure survey technology to locate the high value flow locations v = φv_{0} illustrated in Figure 4, historical wellbore performance Q ~ ρC x T x V exploration is limited to the slow and costly method of locating high production well flow V = 2πr_{0}φv_{0}ℓ structure by drill bit [48].
Three consequences of lognormal distributions of crustal permeability affect wellbore flow velocity heterogeneity V = 2πr_{0}φv_{0}ℓ controlled by Figure 4 crustal flow variations v = φv_{0}:
(i) Early drilling is inefficient as adequately resolved flowstructure information v = φv_{0} is lacking;
(ii) The physical origin of lognormal distributions of permeability imply that spatial flow structure variations v = φv_{0} inherently raise the issue of resource survey spatial resolution relevant to drilling targets;
(iii) Unresolved flow structure variations force average well production for all wells shown in Figure 7 to be nearly 50% of the average well production for “successful” wells; i.e., with inadequate knowledge of flow spatial variation more than 50% of geothermal drilling expenditure is a “sunk cost”.
If follows that in principle surface sensor surveys that can accurately identify the wellbore flow structure spatial variations V = 2πr_{0}φv_{0}ℓ controlled by crustal flow variations v = φv_{0} illustrated in Figure 4 can give geothermal operators prospects for:
(i) Systematically reducing sunk costs;
(ii) Reducing many early stage development drilling failures;
(iii) Promoting small field development through lowered exploratory well drilling costs.
3. Results
Seismic emissions from fluid flow crustal flow velocity structures v = φv_{0} can be reliably and cost effectively imaged by multichannel seismic data processing developed for monitoring the stimulation and production of shale formation hydrocarbons at spatial resolutions of order 25m [9,10,11,12,13,14,15]. Tomographic Fracture Imaging (TFI) builds on the multichannel seismic reflectivity data acquisition and processing technology perfected in the search for and production of hydrocarbon hosting stratigraphic traps in layered geological basin sedimentary sequences [16]. By using large arrays of a statetheart standalone seismic sensor stations, TFI has demonstrated ability to acquire and process large amounts of seismic data to reliably detect and locate smalllevel seismic emissions associated with spatially and temporally persistent crustal flow structures. It should be explicitly noted that TFI flow imaging precisely targets the inherent crustal largescale poroconnectivity structures (Figure 4) that emit lowlevel seismic energy in direct proportion to spatial and temporal flowconnectivity [6,7,8,9,10,11,12,13,14].
Application of TFI from sedimentary settings to convective geothermal flow structure settings requires adaptation of reflection seismics technology. In this section, we first sketch existing field proven TFI acquisition/processing achieved via seismic velocity models gained from reflection seismic imaging in ordered sedimentary crustal sections, then describe the data acquisition/processing adaptation needed to deploy TFI methodology in disorderly convective geothermal terrains that do not admit of seismic reflection velocity structures.
3.1 The principles of TFI Seismic Array Crustal Velocity Structure Imaging in Sedimentary Terrains
Sedimentary basin sequence geological layers typically extend laterally for kms and are demarcated by small but laterally consistent changes in seismic wave speed across the layer interfaces. The changes in seismic wave speed cause downgoing seismic waves to reflect wave energy back towards the crustal surface [16]. The seismic reflection sequence is sketched in Figure 8.
Figure 8 Seismic wave reflection action at formation interfaces in layered sedimentary sequences. The widths of the wave trajectories scale with wave amplitude. Layer seismic velocities tend to increase with depth due to compaction (V4 > V3), but seismic wave energy also reflects when layers decrease in velocity (V3 < V2).
Seismic reflection surveys map the reflectivity sequence of sedimentary layers using thousands of independent sensors deployed along the surface of a crustal survey volume. As outlined in Figure 9, each sensor returns a sequence of reflection signals generated by seismic waves created by surface sources of known positions at known times. By systematically determining the seismic velocity of each layer, weak reflection signals received by the many sensors are timeadjusted and added together to create a “signal stack” for the reflector sequence below the sensors. Reflection signal stack counts are typically of order 50 to 100 per sensor location. Such surveys are conducted across sensor arrays of 50 to 100 sensors per km^{2}, often generating terabytes of seismic data. These data are processed by modern computational power into detailed seismicvelocity stratigraphic sequences of order 1025m vertical and horizon resolution [16]. Drilling the surveyed geological sections routinely validates the seismic velocity stratigraphic imaging process.
The central seismic reflection survey parameter is signal stack count or “fold”. In the Figure 9 seismic reflection survey sketch, survey fold is the number of seismic traces with information about the reflectivity strength of a given patch of crust at a given depth. The larger the number of sensors, larger the fold, and the better the estimate of the reflectivity strength and the better the seismic image. A typical reflection seismic trace fold is half the number of sensor stations times the station interval per source interval, Nfold = Nsns x Dsns / 2 x Dsrc. I.e., for equal source and sensor spacing, Dsns = Dsrc, with a source point at each sensor offset, Nsrc = Nsns = 6, the survey fold is 6/2 = 3 for reflector points below the source and sensor. Shot points too far from a sensor to be likely to register useful reflector point information do not contribute to the sensor fold. Reflection seismic sensor fold is thus computed for “active” rather than total sensors and sources. Sources that are active for a given sensor typically number in the range of 128256, giving an effective seismic reflectivity image fold of 64 to 128.
Figure 9 How seismic reflection survey sensors and sources build reflector signal count or “fold”. Sensors record reflected wave energy according to ray paths specific to a given source point. Signal fold is built by summing properly timeadjusted seismic wave energy from sensorsource ray paths common to specific reflection point. The proper timeadjustments to create accurate seismic images is found by trailanderror computation [16]. Building the reflection image fold requires extensive computation to determine the seismic wave traveltimes assigned to each geological layer within a crustal section. Determining seismic wave traveltimes in reflection seismics is possible because the location and times of all seismic sources are known.
From Figure 9, each of Nsns sensor traces contributing to image fold for each crustal patch at any given depth is associated with active source points of known location and time. Knowledge of each active source and sensor location and time for each seismic trace in a signal stack gives sufficient information for seismic data processing to construct a model of the seismic velocity structure at each point of the survey volume [16]. A reflection seismic image is thus effectively a pointbypoint seismic velocity model of the crustal geological formation sequence.
TFI proceeds according to the signal stacking principle of Figure 8 and Figure 9, with signal stack fold depending directly on summing data over large numbers of surface sensors Nsns. TFI signal fold building depends, however, on microseismic energy sources of unknown locations and times and must follow a different data acquisition and processing path from the seismic reflection imaging process sketched in Figure 8 and Figure 9. More particularly, TFI requires knowing in advance the seismic velocity structure of the survey volume. Processing TFI data for sedimentary terrains deploys the same degree of computing power as reflection seismics, but computation is directed towards finding unknown seismic emission sources in a known velocity field rather than finding unknown seismic velocities in a known sourcesensor field.
In more physical terms, TFI acquisition and processing differs from reflection seismic acquisition and processing because the target crustal fluid flow structures are treated as a series of point seismic sources which are different from the laterally extensive layers of sedimentary sequences reflecting seismic energy from (known) seismic sources. Sedimentary formation welllog sonic seismic velocities typically vary by no more than 6% from a formation mean velocity, and long wavelength seismic reflection amplitudes are spatially averaged over most of this velocity variation. Reflection seismic image stacks comprise coherent sums of seismic energy along the length of sedimentary interfaces.
For crustal fluid flow structures, by contrast, the characteristic 6% variations in seismic velocity measured by sonic welllogs are pointstructures driven by porosity variations in fracturedensity and, more importantly, by variations in fracture connectivity. Flowcontrolling formation permeability fluctuations along spatially erratic poroconnectivity structures within formations do not act coherently to form seismic wave reflection amplitudes. Flowspecific spatial variations make no contribution to Figure 8 and Figure 9 seismic reflected energy sequences. Rather, as learned from the §2 discussion of crustal microseismicity, TFI empirics establish that fluid flow in crustal permeability structures are sources of seismic energy emission by which TFI seismic flow images are generated [9,10,11,12,13,14]. The principles of TFI data acquisition and processing are thus to identify and locate the sources of spatially and temporally persistent seismic energy emission generated by fluid flowing along erratic percolation pathways within a crustal volume [Figure 4].
Because TFI signal acquisition lacks both source times and source locations from which to build a velocity model as in Figure 8 and Figure 9, TFI acquisition requires that an effective velocity model be generated in advance of image data processing. The TFI effective velocity model need not, however, be pointbypoint velocity distribution as generated for seismic reflection surveys. It is sufficient for TFI to have a comprehensive set of traveltimes TT(i,j) for seismic wave paths through the survey volume from a subsurface source point i to a surface sensor j. TFI in shale formations can rely on reflection seismic velocity models for the crosstable TT(i,j).
Volcanic terrains, however, do not support reflection seismic imaging. In consequence, TFI in convective geothermal systems depends on the sourcesensor geometry sketched in Figure 10, where traveltimes TT(i,j) connect the i^{th} source point ( the red box) in the survey volume to the j^{th} sensor at the surface (blue dots). Acquisition and utilisation of crosstable traveltime sets TT(i,j) in volcanic terrains is described in §3.2.
Provided with an effective velocity model of the crustal survey volume, TFI processing proceeds as indicated by the Figure 10 survey geometry. Let us consider a specific seismic energy source volume element or “voxel” represented by the red box within the Figure 10 survey volume. The key TFI interest is knowing if that particular voxel is emitting a temporal succession of lowlevel seismic energy due to fluid flow through the test voxel in flow communication with adjacent voxels. Any particular voxel 1 < i < Nvox may or may not be a source of flowrelated seismic emission in any given time interval Δt. If the voxel is emitting seismic energy in the time interval Δt, the previouslydetermined crosstable of seismic wave traveltimes TT(i,j) between the specified voxel location i and surface seismic sensor locations j=1…Nsns is used to add up the appropriately timeshifted seismic sensor interval data to build an Nsns ~ 1000 strong seismic emissions signal stack. The signal stack is built by summing sensor data for timeinterval Δt with each sensor record adjusted according to the arrival times of a voxelemitted energy given by the crosstable traveltime TT(i,j) from voxel i to sensor j. If during the given time interval Δt, the given voxel emits flowrelated seismic energy, TFI processing coherently sums the seismic energy across all surface sensors to give an Nsnsfold signal stack for each voxel for the given time interval. The same stacking procedure is then performed for each voxel over sequences NΔt of time intervals to give an Nsns x NΔtfold fold signal for each voxel.
Figure 10 The seismic sensor element of the TFI signal stacking fold is the number Nsns of surface sensors denoted by the grid of blue dots. For TFI, the surface sensors record flowrelated source activity emitted over short time intervals Δt by a source “voxel” (volume element) represented by the red box. TFI data stacks for any given source voxel are a double sum over a sequence of NΔt time intervals and Nsns sensors. The seismic wave traveltime TT(i,j) from any given source voxel 1 ≤ i ≤ Nvox to any given surface sensor 1 ≤ j ≤ Nsns is computed from a known velocity model or from traveltime data acquired for known seismic source activity at the wellbore depth indicated by the red asterisk. TFI performed over shale formation stimulation crustal volumes establishes crustal flow structure images connecting Nvox seismic emitting voxels acquired by surface seismic reflection geometry illustrated in F Figure 8 and Figure 9. For volcanic terrains, the necessary observational traveltime data TT(i,j) are acquired as discussed in §3.2.
Figure 11 illustrates the final TFI stacking step covers the entire set of voxels. The Figure 11 image produced by TFI data acquisition and processing is essentially that of a triplestack crustal flow structure count Nfold = Nvox x NΔt x Nsns built up from Nvox spatial connected flowactive voxels each defined by NΔt x Nsns seismic emission samples. We note in Figure 11 that 25m spatial resolution of standard seismic stratigraphy can distinguish the most significant static geological structural features of a crustal volume. By acquiring fluidflowgenerated seismic emission signals, the Figure 11 25m image spatial resolution accurately locates the fracture specific fluidactive flow structures denoted by the red lines. The timesequence of mapview flowstructures then shows the black patches of timeevolving seismic emission voxels following strands of the previously mapped fluidactive fracture structures.
Figure 11 Time sequence of shale formation fluid flow seismic activation (black smudges) generated by frack stimulation of wellbore (black line) at intersection of the wellbore and TFI redline fracture connectivity structure. Redline fracture connectivity structures inferred from TFI images prior to wellbore stimulation. Wellbore stimulation validated fractureinterpretation of TFImapped redline fracture connectivity structures. Time sequence intervals in minutes; wellbore ~ 1km long; spatial resolution of image elements is 25m. The unpublished TFI data were recorded in the Barnett shale of Texas.
3.2 Estimating the TFI TravelTime CrossTable TT(i,j) in Volcanic Terrains
TFI imaging of convective geothermal flow systems proceeds in the same manner as TFI in shale formations illustrated in Figure 11 provided that the necessary seismic traveltime data TT(i,j) exist for crustal volume hosting the convective flow system. As seen in Figure 12, however, traveltime data TT(i,j) cannot be acquired through the reflection seismic imaging of sedimentary sequences that embed shale formations. In contrast to seismic waves traveling in longsettled orderly basin sedimentary layer sequences sketched in Figure 8 and Figure 9 and illustrated in the upper panel of Figure 12, seismic waves traveling in volcanic and magmatic terrains illustrated in the lower panel are disrupted at all scales [17]. Welllog physical properties within sedimentary layers are disrupted at ~ 6% deviation from the mean, but when these spatial variations are averaged over seismic wavelengths, orderly differences between individual formation layers generally exceed disorderly differences within formations. No such lateral order is imposed over the internal disorder for intrusive magmatic and extrusive volcanic terrains.
Figure 12 Comparison of sedimentary layering (above) and absence of layering in magmatic intrusive (below). Views are to NW and SE, respectively, across the Rio Grande River plain in the San Juan volcanic fields near Creede, Colorado, USA.
In accordance with the Figure 12 crustal rock order/disorder contrast, TFI applied to convective geothermal terrains does not attempt to process surface reflection seismic array data in terms of orderly seismic wave fronts implicit in Figure 8 and Figure 9. Instead, a geothermal terrain TFI survey needs to estimate the seismic traveltime crosstable TT(i,j) via seismic refraction as sketched in the Figure 13 version of the Figure 10 wellboreseismicsource and surface seismic sensor array geometry.
Figure 13 shows the refraction seismic elements needed for estimating TT(i,j) crosstable data for conducting TFI surveys in volcanic terrains. A 30MJ charge of propellant is ignited in an exploration or production wellbore to generate refraction seismic wave energy (green blob) that passes through the crustal survey volume (black dots) to the surface sensor array (blue dots). The observed refraction data, denominated TT(0,j) for sensors 1 ≤ j ≤ N_{sns}, are then processed to estimate the crosstable TT(i,j) by which to compute TFI stacks for each of the voxels 1 ≤ i ≤ N_{vox}. Figure 14, Figure 15 and Figure 16 illustrate the numerical simulations that establish the TT(i,j) estimation procedure specific to volcanic terrains [17].
Figure 13 Elaboration of Figure 10 detailing the sourcesensor geometry for seismic traveltime data TT(0,j) acquisition from which to estimate the voxelsensor traveltime crosstable TT_{est}(i,j), 1 ≤ i ≤ Nvox and 1 ≤ j ≤ Nsens within the crustal volume outlined by black dots [17]. The green dot represents a seismic wave sourced from 30MJ of propellant energy generated in the wellbore (red line) to be recorded at the surface seismic sensor array (blue dots).
Figure 14 (left) illustrates a 100 x 100 x 100 node generic crustal numerical seismic velocity cube characteristic of the ambient crust [6,7,8]. The generic crustal velocity variations derive from two geological processes:
(i) Spatial fluctuations with 68% RMS deviation due to variations of porosity distributed according to the powerlaw scaling spectral distribution Χ(k)^{2} ~ 1/k, 1/km < k < 1/Dm, observed in welllogs worldwide [6,7,8];
(ii) Nominal vertical and lateral velocity gradients of amplitude [gx,gy,gz].
Taking Figure 14 velocity cube to be the “true” velocity structure defining the fieldobserved sourcetosurface traveltimes TT_{true}(0,j), we proceed to search for generic trial velocity cubes that provide estimated sourcetosurface traveltime sets TT_{est}(0,j) that approximate the TT_{true}(0,j); i.e., we look for trial velocity distributions such that TT_{est}(0,j) ~ TT_{true}(0,j). Figure 14 (right) displays an example of the residual distribution TT_{est}(i,j)  TT_{true}(i,j) from such searches. For sourcetosensor travel times of order 1 second, the typical RMS residual difference TT_{est}(i,j)  TT_{true}(i,j) is ~ 3msec with ~ 3msec standard deviation. 3msec temporal fluctuations are equivalent to 15m spatial fluctuations.
Figure 14 (Left) A 100 x 100 x100 node distribution of seismic velocities comprising (i) spatial fluctuations with powerlawscaling power that scales inversely with spatial frequency k, X(k)^{2} ~ 1/k, and (ii) vertical and horizontal gradients; the colourbar codes velocities in m/s; the velocity cube defines the cross/table TT_{true}(i,j). (Right) The distribution of traveltime residuals TT_{est}(i,j)  TT_{true}(i,j) for the lefthand velocity distribution; the colourbar codes the residuals in msec; the RMS residual is ~ 3msec with 3msec standard deviation.
Figure 15 illustrates the extent to which it is possible to find trial seismic velocity cubes for which the estimated seismic traveltimes j = 1. N_{sns} for a nominal voxel source i TT_{est}(i,j) in red closely overlie the “true” traveltimes TT_{true}(i,j) in blue.
Figure 15 Overlay of Nsns surface sensor traveltimes from a nominal source voxel i computed from “true” velocity field (blue) and from an estimated velocity field (red). As the red traveltimes closely approximate the blue traveltimes, TFIprocessed Nsnsfold stacks of voxel seismic emissions will be accurate representations of actual seismic emissions recorded by the surface sensor array for the i^{th} voxel.
Figure 16 gives the result of using suitable estimation velocity cubes to locate nine instances of arbitrary source voxels denoted by the red asterisk within the “true” velocity cube of Figure 14 (left). For each asterisk source location of traveltime crosstable TT_{true}(i,j), six black circles denote a sequence of estimated positions for the source voxel derived from an estimated crosstable TT_{est}(i,j). The bestfit black circles correspond to the actual source location with an effective accuracy of 1.3 nodes. The location procedure can be repeated for multiple estimated travelcross tables TT_{est}(i,j) to strengthen the statistical certainty of the source location. The Figure 16 caption gives details.
Figure 16 Display of sourcevoxel search algorithm performance. For the Figure 14 (left) “true” velocity field with sourcesensor traveltime crosstable TT_{true}(i,j), nine random sourcevoxel locations are selected, i_{m}, m = 1….9. For each selected sourcevoxel, the location i_{m} is identified by a red asterisk. Six values of estimated source positions 1 ≤ i_{est} ≤ 6 are plotted as black circles. The six i_{est} locations are those that most nearly minimise the traveltimeresidual function ∑_{j=1…Nsns} TT_{true}(i_{m},j)  TT_{est}(i_{est},j). For each test voxel location i_{m}, the estimated source positions i_{est} with the smallest residuals successfully locate the test voxel. Estimated locations with higher residuals show the potential for source voxel location error. Performing the location search for several estimated traveltravel crosstables TT_{est}(i,j) reliably reduces the statistical sourcevoxel location error. The statistical accuracy of estimated sourcevoxel locations is ~1.3 nodes.
4. Discussion
The TFI timelapse sequence of flowinduced seismic emission imaging shown in Figure 11 illustrates that Nvox sets of Nsns traveltime tables repeated for NΔt timeintervals of TFI listening data can produce 1525m accurate crustal flow activity maps in sedimentary settings for which seismic reflection imaging provides accurate seismic velocity models. Figure 14, Figure 15 and Figure 16 illustrate the means by which seismic refraction data acquired in the Figure 13 sourcesensor geometry can provide adequate velocity models by which to perform TFI on convective geothermal systems at spatial resolution ~ 50m. Creating flowstructure images that resolve flow heterogeneity v = φv_{0} illustrated in Figure 4 at ℓ ~ 50m spatial resolution can directly give convective geothermal resource operators the information needed to efficiently drill production wellbores defined by thermal power output Q ~ ρC x T x V, where wellbore volumetric flow V =2πr_{0}φv_{0}ℓ is expressed in terms related to heterogeneous flow images giving v = φv_{0} at scale lengths ℓ ~ 50m.
The operational implications for presentday geothermal energy production are summarised in the two parts of Figure 17 [49,50]. At left is 50 years of steady rise in global energy consumption; at right is evidence of slowing rather than advancing geothermal power production growth. The slowing growth in the contribution of geothermal power to the global energy economy is arguably due to the high cost and risk of identifying a suitable convective geothermal resource and estimating its size and productivity. While there have been improvements in surface exploration technology that can in principle develop ample geothermal prospects in a range of markets, Figure 17 makes it evidence that to date these improvements have done little to mitigate the capital cost and risk of exploration drilling codified in the wellbore production expression Q ~ ρC x T x V. As we have discussed, current exploration techniques do not accurately predict where to drill. Neither do these surveys provide reliable estimates of the resource fluid flow capacity as measured by the ability to drive an electricitygenerating turbine. Without improved assessment of geothermal flow system structure denoted by volumetric flow broken down into terms relevant to flow heterogeneity, i.e., the simple and straightforward expressions Q ~ ρC x T x V and V =2πr_{0}φv_{0}ℓ connecting wellbore performance to crustal flow uncertainty, drilling production wellbores remains a significant investment risk. The risk is greatest at the initial stages of a project before knowing whether the geothermal resource has enough capacity to recover the drilling costs, e.g., Fig 0.1 [51]. Early stage test drilling can account for 15 percent of the overall capital cost at a point when the risk of project failure is still high. These upfront risks have a knockon effect throughout the project from both technical and financing perspectives.
Figure 17 (Left) Steady growth of global energy consumption for last 50 years [49]. (Right) Slowing growth of geothermal power production (red) and stationary growth in capacity (blue) over last 20 years (diamonds = United States; squares = Indonesia/Philippines) [50].
To put numbers to the “cost of exploration”, a representative capital outlay for initial drilling of a prospect is USD25M for two to four wells [51]. Further drilling costs of order USD50100M are predicated on (1) the initial drilling outcome being one successful test production well and one successful test reinjection well, and (2) a 75% success rate of subsequent drilling returning average production well output of Q ~ 7MWe. Figure 6 (right) indicates that the global norm for “successful” geothermal production wells at 34MWe is one half of this expectation. From inspection of the Figure 7 early stage drilling success rates, both the expected outcome of the initial USD25M test drilling and the success profile of the followon USD50100M development drilling are far from reassuring to an investor.
Given that there is a clear global demand for clean baseload electric power together with an equally clear availability of geothermal resources to supply that power, there is every reason to support the World Bank conclusion that the problem facing geothermal resource development lies with persistent uncertainty in drilling wellbores that access highflow zones in convective geothermal flow system [47]:
“There is no reasonable basis for forecasting the probability of success in the Exploration Phase of a project (first 5 wells). … This low average success rate across exploration wells serves to confirm the high risks of initial drilling, as well as driving high return expectations for earlystage investors. It is clear, therefore, that every effort should be made to reduce risk in the Exploration Phase.”
The cost profile of performing the adapted TFI flowstructure surveys on convective geothermal systems is a small percentage of the potential drilling cost saving. Established TFI data acquisition uses standard standalone reflection seismic sensor modules available at nominal bulk rental rates, and TFI data processing is performed at commercial rates set by reflection seismic data processing [9,10,11,12,13,14]. Performance of wellbore seismic refraction seismic sourcing for TFI surveys on convective geothermal terrains is based on established wellbore deflagration stimulation technology [15]. Processing wellbore seismic refraction traveltime data to acquire the TFI sourcevoxeltosensor traveltime crosstable is performed via Matlabimplemented standard Fast Matching Method wavefront propagation computation. Matlab performed all computations discussed in this paper. Neither field data acquisition technology/cost, nor data processing technology/cost, should be seen as impeding implementation and development of a TFI capability for progressing convective geothermal brownfield or greenfield operations.
We summarise the cost implications our discussion in terms of specific numbers from a specific class of convective geothermal resources, that of the island of Sumatra [52]. Developing Indonesia’s largely untapped ~1150MWe convective geothermal power production potential in Sumatra requires drilling ~ 160 production wellbores of mean productivity at 7MWe mean wellbore productivity. Judging from Figure 6 and Figure 7 present rates of drilling success, we might easily double the estimated number of wells required in realising the Sumatra promise without improved drilling success. At present and at ~ USD5M per production well, accessing Sumatra’s geothermal power potential requires an investment of order USD1.5B.
From the example of Sumatra alone, if TFI provides the technical means to reduce the risk of geothermal production well drilling by 25% to 50%, the benefits from raw project drilling cost saving and from downstream project planning and enabling are measured in the hundreds of millions of dollars. In the world beyond the Sumatra prime resource, there are a large number of geothermal sites of lesser potential where TFI drilling derisking can be proportionally more productive and enabling. Reducing convective geothermal drilling uncertainty through understanding and deploying TFI flow structure mapping technology can save many millions of dollars to the benefit of millions of people.
5. Conclusions
Crustal reservoir flowstructure imaging technology using multichannel seismic surface array data acquisition and processing is proven for sedimentary shale formations. We have shown how this shale formation flow imaging technology can be adapted to volcanic terrains hosting convective geothermal flow systems. The difference between sedimentary and volcanic settings is essentially only that sedimentary formations acquire the necessary seismic velocity structure information via seismic reflection data while volcanic terrains will require use of seismic refraction data. Simulated seismic refraction data acquired via subsurface propellant seismic sources recorded by surface seismic sensor arrays provide sufficient seismic travel time map detail to locate flowrelated seismic emission sites with 50m spatial resolution within a geothermal reservoir survey volume. Connecting the imaged seismic emission sites into convective flow structures with 50m spatial resolution gives highly accurate targets for exploration and production well drilling. As current development of convective geothermal systems depends on costly drilling with spatial targeting control of perhaps 500 meters, lowcost remote crustal permeability images giving productive/sustainable convection flow structures with 50m accuracy can greatly improve the return on drilling investment for both brownfield and greenfield convective geothermal resources.
Acknowledgments
The authors acknowledge with pleasure Tom Fleure, Ambient Reservoir Monitoring, Houston TX, and Ontowiryo Alamsyah, Bambang Pramono and Wahyuddin Diningrat, Star Energy, Jakarta ID, for extensive conversations on the content, procedures, and substance of this paper. Essential aspects of ambient crustal microseismicity magnitude distributions were established via cooperative assess to the st1 Deep Heat EGS stimulation project managed by Tero Saarno.
Author Contributions
PL, PM, GS and CS are jointly responsible for the material, interpretations, results and opinions presented in this paper. PL prepared the text; PM partners on field procedures; GS partners on capital acquisition; CS partners on data processing.
Competing Interests
The authors contribute this paper as an act of scientific advancement without commercial interest.
References
 Leary PC, Malin PE, Saarno T. A physical basis for the gutenbergrichter fractal scaling. Proceedings of the 45th Workshop on Geothermal Reservoir Engineering. 2020 February 1012; Stanford, California, US. Stanford: Stanford University.
 Leary P, Malin P, Saarno T, Heikkinen P, Diningrat W. Coupling crustal seismicity to crustal permeability  Powerlaw spatial correlation for EGSinduced and hydrothermal seismicity. Proceeding of the 44th Workshop on Geothermal Reservoir Engineering. 2019 February 1113; Stanford, California, US. Stanford: Stanford University.
 Leary P, Malin P, Saarno T, Kukkonen I. αφ~ αφcrit  Basement rock EGS as extension of reservoir rock flow processes. Proceeding of the 43th Workshop on Geothermal Reservoir Engineering. 2018 February 1214; Stanford, California, US. 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. DOI:10.3390/en10121979. [CrossRef]
 Maxwell S, Deere J. An introduction to this special section: Microseismic. Leading Edge. 2010; 29: 277. DOI:10.1190/1.3353722. [CrossRef]
 Leary P, Malin P, Niemi R. Fluid flow and heat transport computation for powerlaw scaling poroperm media. Geofluids. 2017; 2017. DOI:10.1155/2017/9687325. [CrossRef]
 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. Heterogeneity in the Crust and Upper Mantle: Nature, Scaling, and Seismic Properties. Boston: Springer US; 2003. pp. 155186. [CrossRef]
 Geiser P, Lacazette A, Vermilye J. Beyond ‘dots in a box’: An empirical view of reservoir permeability with tomographic fracture imaging. First Break. 2012; 30: 6369.
 Geiser P, Leary P. Tomographic Fracture Imaging (TFI): Direct 5D mapping of transmissive fracture/fault zones using Seismic Emission Tomography (SET). Proceeding of 39th Workshop on Geothermal Reservoir Engineering. 2014 February 2426; Stanford, California, US. Stanford: Stanford University.
 Sicking C, Vermilye J, Lacazette A, Yaner A, Klaus A, Bjerke L. Case study comparing microearthquakes, fracture volumes, and seismic attributes. SPE/AAPG/SEG Unconventional Resources Technology Conference. Denver, Colorado: Unconventional Resources Technology Conference; 2014. DOI:10.15530/urtec20141934623. [CrossRef]
 Ross J, Parrott K, Vermilye J, Klaus A. Tomographic fracture imaging: Examples of induced fracture and reservoirscale observations during wellbore stimulations, Niobrara and Bakken plays, USA. Leading Edge. 2017; 36: 437444. DOI:10.1190/tle36050437.1. [CrossRef]
 Sicking C, Vermilye J, Yaner A. Forecasting reservoir performance by mapping seismic emissions. Interpretation. 2017; 5: T451T459. DOI:10.1190/INT20150198.1. [CrossRef]
 Malin PE, Leary PC, Cathles LM, Barton CC. Observational and critical state physics descriptions of longrange flow structures. Geosciences. 2020; 10: 50. DOI:10.3390/geosciences10020050. [CrossRef]
 Leary P, Malin P, Saunders G, Fleure T, Sicking C, Pullammanappallil S. Flowimaging of convective geothermal systems  obtaining seismic velocity models needed for production well targeting. Proceedings of the 45th Workshop on Geothermal Reservoir Engineering. 2020 February 1012; Stanford, California, US. Stanford: Stanford University.
 Waters KH. Reflection seismology: A tool for energy resource exploration. New York: WileyInterscience Publication; 1981.
 Bannister S, Melhuish A. Seismic scattering and reverberation, Kaingaroa plateau, Taupo Volcanic Zone, New Zealand. N Z J Geol Geophys. 1997; 40: 375381. [CrossRef]
 Jaeger JC, Cook NG, Zimmerman R. Fundamentals of rock mechanics. 4th Ed. Oxford: Blackwell Publishing; 2009.
 Ingebritsen SE, Sanford WE, Neuzil CE. Groundwater in geologic processes. 2nd Ed. Cambridge: Cambridge University Press; 2006.
 Ritzi Jr RW, Bobeck P. Comprehensive principles of quantitative hydrogeology established by Darcy (1856) and Dupuit (1857). Water Resour Res. 2008; 44. DOI:10.1029/2008WR007002. [CrossRef]
 Brown GO. Subsurface HydrologyHenry Darcy and the making of a law. Water Resour Res. 2002; 38: 11. DOI:10.1029/2001WR000727 [CrossRef]
 Slichter CS. Theoretical investigation of the motion of groundwaters. The 19th Annu Rep. US Geophy Survey. 1899: 304319.
 Meinzer OE. The occurrence of groundwater in the United States, with a discussion of principles. Chicago: University of Chicago; 1923.
 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 AGU. 1935; 16: 519524. [CrossRef]
 Muskat M, Botset HG. Flow of gas through porous materials. Physics. 1931; 1: 2747. [CrossRef]
 Muskat M. The flow of fluids through porous media. J Appl Phys. 1937; 8: 274282. DOI:10.1063/1.1710292. [CrossRef]
 Muskat M. The flow of homogeneous fluids through porous media: Analogies with other physical problem. New York: McGrawHill; 1937. [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. Hydrol Sci J. 1957; 2: 2359. [CrossRef]
 Bear J. Dynamics of fluids in porous media. New York: American Elsevier Publishing Company, Inc; 2013.
 Bachmat Y, Bear J. On the concept and size of a representative elementary volume (Rev). Advances in Transport Phenomena in Porous Media. Dordrecht: Springer; 1987. DOI:10.1007/9789400936256_1 [CrossRef]
 Nordahl K, Ringrose P. Identifying the representative elementary volume for permeability in heterolithic deposits using numerical rock models. Math Geosci. 2008; 40: 753771. DOI:10.1007/s1100400891824 [CrossRef]
 Slichter C. Field measurements of the rate of movement of underground waters. US Gelogical Survey WaterSupply and Irrigation Paper. 1905; 140.
 Law, J. A statistical approach to the interstitial heterogeneity of sand reservoirs. Trans AIME. 1944; 155: 202222. [CrossRef]
 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]
 King Hubbert M, Rubey WW. Role of fluid pressure in mechanics of overthrust faulting: I. Mechanics of fluidfilled porous solids and its application to overthrust faulting. Geol Soc Am Bull. 1959; 70: 115166. [CrossRef]
 Hoek E. Brittle fracture of rock. Rock mechanics in engineering practice. Hoboken: John Wiley & Sons: 1968.
 Shi Y, Bolt BA. The standard error of the magnitudefrequency b value. Bull Seismol Soc Am. 1982; 72: 16771687. [CrossRef]
 Aki K. A probabilistic synthesis of precursory phenomena. Earthquake prediction: An international review. American Geophysical Union; 1981; 4: 566574. [CrossRef]
 Hirata T. A correlation between the b value and the fractal dimension of earthquakes. J Geophys Res Solid Earth. 1989; 94: 75077514. [CrossRef]
 Zaky DA, Nugraha AD, Sule R, Jousset P. Spatial analysis of earthquake frequencymagnitude distribution at geothermal region south of bandung. Proceeding of 9th International Workshop on Statistical Seismology; 2015 June 1418; Potsdam.
 Gibbs JF, Healy JH, Raleigh CB, Coakley J. Seismicity in the Rangely, Colorado, area: 19621970. Bull Seismol Soc Am. 1973; 63: 15571570.
 Utsu T. Representation and analysis of the earthquake size distribution: A historical review and some new approaches. Seismicity Patterns, their Statistical Significance and Physical Meaning. Basel: Birkhäuser; 1999. pp. 509535. [CrossRef]
 Wiemer S, Wyss M. Minimum magnitude of completeness in earthquake catalogs: Examples from Alaska, the western United States, and Japan. Bull Seismol Soc Am. 2000; 90: 859869. [CrossRef]
 Woessner J, Wiemer S. Assessing the quality of earthquake catalogues: Estimating the magnitude of completeness and its uncertainty. Bull Seismol Soc Am. 2005; 95: 684698. DOI:10.1785/0120040007 [CrossRef]
 IFC, Success of geothermal wells: A global study, international finance corporation. 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.
 Sanyal SK, Morrow JW, Jayawardena MS, Berrah N, Li SF. Geothermal resource risk in Indonesia–A statistical inquiry. Proceeding of 36th Workshop on Geothermal Reservoir Engineering. 2011 January 31  February 2; Stanford, California, US. Stanford: Stanford University.
 BP Statistical Review of World Energy. London: BP P.L.C; 2016. Available from: https://www.bp.com/content/dam/bp/businesssites/en/global/corporate/pdfs/energyeconomics/statisticalreview/bpstatsreview2019fullreport.pdf.
 Bertani R. Geothermal power generation in the world 2010–2014 update report. Geothermics. 2016; 60: 3143. [CrossRef]
 Gehringer M, Loksha V. Geothermal handbook: Planning and financing power generation. Energy Sector Management Assistance Program (ESMAP). Available from: https://www.esmap.org/sites/esmap.org/files/DocumentLibrary/FINAL_Geothermal%20Handbook_TR00212_Reduced.pdf
 Asian Development Bank and The World Bank. Unlocking Indonesia’s Geothermal Potential. © Asian Development Bank and The World Bank. 2015. Avaliable from: https://openaccess.adb.org; https://openknowledge.worldbank.org.