Seismic Imaging of Convective Geothermal Flow Systems to Increase Well Productivity

Peter Leary ^{1, 2, †, *}, Peter Malin ^{2, †}, Graeme Saunders ^{1, †}, Charles Sicking ^{3, †}

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

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

3. Ambient Reservoir Monitoring, Peninsula Place, Houston, TX 77459, USA

† These authors contributed equally to this work.

* Correspondence: Peter Leary

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): 012; 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.

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 high-resolution flow structure mapping can eliminate many of the 20-30 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 Gutenberg-Richter (G-R) frequency-size relation for microseismicity in the ambient crust. We first note that observed ambient crust microseismicity has a lognormal frequency distribution rather than the traditional power-law 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 well-log, well-core and well-productivity 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 microseismicity-permeability interactions across mm-km 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 hydrocarbon-bearing 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 flow-imaging 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], well-log, well-core and well-productivity 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 multi-scale nature of fluid flow in the ambient crust. The long-held 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 well-logs, well-cores, and well-production 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 multi-channel seismic data to detect low-level 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 low-cost remote seismic sensing can eliminate a great deal of uncertainty in drilling high-cost 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 ever-present 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 Fluid-Rock 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 mid-19^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 small-scale 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 quasi-uniform “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 small-scale and sparsely sampled well-log and well-core 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 well-log and well-core flow property sampling can be “upscaled” for reservoir-scale 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 well-known 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 N-point 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 power-spectrum of χ(x) is independent of spatial frequency k, |Χ(k)|² ~ 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 small-scale sampling of reservoir-scale 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. Well-log rock physical property sequences in geological settings worldwide systematically violate the |Χ(k)|² ~ const spectral condition required for Figure 1. Instead, well-log Fourier spectra everywhere obey the power-law scaling function |Χ(k)|² ~ 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].

Well-log spectral empiric |Χ(k)|² ~ 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 long-standing empirical Gutenberg-Richter (G-R) 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 G-R empirical constant b typically has a narrow range of values 0.75 < b < 1.5, the G-R 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 G-R relation applied to seismic activity in Indonesian geothermal fields [42]. The empirical G-R relation parameter b is evaluated as the slope of the straight-line power-law fit to the observed decline trend (red-squares) 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 Gutenberg-Richter 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 G-R 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 b-value parameter. The cumulative data format disguises the event number decline for which the standard G-R 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 small-magnitude event number decline in fitting the G-R 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 G-R power-law 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 G-R 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 power-law 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 low-magnitude and high-magnitude 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 magnitude-frequency 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 αφ ~ 3-4 observed for crustal porosity-permeability 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 G-R power-law 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 Gutenberg-Richter power-law scaling relation fails to represent a fundamental set of fluid-rock interaction physical processes described in [2,3,4] that occur at all scales in the ambient crust and generate the spectral empiric |Χ(k)|² ~ 1/k [7,8]. The observed 1/k power-law scaling empiric arises as crustal rock undergoes consolidation from high-porosity disordered uppermost crustal sediment granularity to low-porosity ordered lower-crust ductile metamorphic granularity. Crustal fluids play an essential role in consolidation by mediating the change in pore-connectivity, i.e., permeability, of the consolidating rock. Pore connectivity is a statistical phenomenon that conditions the transition from grain-scale disorder in the disaggregated uppermost crust to grain-scale order in the aggregated lower crust. Disordered granularity in the disaggregated uppermost crust is characterised by uncorrelated random fluctuations with flat spectral flat power-spectra |Χ(k)|² ~ 1/k⁰ ~ const. Ordered granularity in the aggregated lower crust is characterised by highly-correlated random fluctuations with steeply declining spectral power-spectra |Χ(k)|² ~ 1/k² arising from the macroscopic ductile flow interfaces of lower crust metamorphics. The order-disorder 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 2-Figure 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 well-core poroperm relation, δφ ~ δlog(κ), for which the proportionality constant α is widely attested to obey the constraint αφ ~ 3-4 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 αφ ~ 3-4 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 G-R interpretation of half of microseismicity event populations, the Figure 3 fit to Figure 2 microseismicity data is essentially a zero-free-parameter prediction for the range of event magnitudes derived from fundamental crustal rock-fluid 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 magnitude-frequency distribution. First, microseismicity is inherently coupled to fluid-rock properties of crustal rock rather than to purely elastic fracture-mechanical 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 large-scale 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₀

Figure 2 and Figure 3 give microseismicity-based evidence that permeability spatial variations κ(x,y,z) ~ exp(αφ(x,y,z)) pervade crustal volumes at all scales to influence the distribution of ambient fluid-related slip event magnitudes. Building on the microseismicity evidence for crustal fluid-rock 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₀φv₀ℓ is controlled by spatially variable crustal fluid flow v = φv₀ 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)|² ~ 1/k0 ~ const. Crustal flow simulation in Figure 4 illustrates how spatial flow connectivity structures controlled by the Fourier spectral power-law scaling function |Χ(k)|² ~ 1/kβ, β ~ 1.2 ± 0 affect crustal permeability distribution κ(x,y,z) ~ exp(αφ(x,y,z)). Well-log 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₀ (red, blue arrows) for crustal permeability sections generated by two random number realisations of spatially correlated random porosity fields with |Χ(k)|² ~ 1/k spectral power scaling [7,8]. While such velocity fields have the same population of high- and low-flow areas, the spatial organisation of the velocity fields is subject to essentially infinite variety. Drilling in the vicinity of high-flow clustering is more productive than drilling in the vicinity of low-flow clustering. No sparse small-scale sampling algorithm suffices to represent the illustrated large-scale flow variations significant to convective geothermal wellbore production. Instead, to locate the high-flow clusters for optimal drilling, crustal fluid flow variations v = φv₀ 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₀ 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 high-flow 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 high-connectivity 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 fracture-connectivity 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₀ amplitudes for four spatially- correlated random porosity fields. Red histogram bars denote the high-flow velocity structures v = φv₀ that would be logical targets for drilling high production Q ~ ρC x T x V wellbores for V = 2πr₀φv₀ℓ. Histogram axes are normalised to maximum value of computed crustal flow velocity v = φv₀. 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 well-log/well-core 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 well-productivity 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₀ 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₀φv₀ℓ structure by drill bit [48].

Three consequences of lognormal distributions of crustal permeability affect wellbore flow velocity heterogeneity V = 2πr₀φv₀ℓ controlled by Figure 4 crustal flow variations v = φv₀:

(i) Early drilling is inefficient as adequately resolved flow-structure information v = φv₀ is lacking;

(ii) The physical origin of lognormal distributions of permeability imply that spatial flow structure variations v = φv₀ 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₀φv₀ℓ controlled by crustal flow variations v = φv₀ 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₀ can be reliably and cost- effectively imaged by multi-channel 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 multi-channel 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 state-the-art stand-alone seismic sensor stations, TFI has demonstrated ability to acquire and process large amounts of seismic data to reliably detect and locate small-level 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 large-scale poro-connectivity structures (Figure 4) that emit low-level seismic energy in direct proportion to spatial and temporal flow-connectivity [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 down-going 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 time-adjusted 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², often generating terabytes of seismic data. These data are processed by modern computational power into detailed seismic-velocity stratigraphic sequences of order 10-25m 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 128-256, 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 time-adjusted seismic wave energy from sensor-source ray paths common to specific reflection point. The proper time-adjustments to create accurate seismic images is found by trail-and-error computation [16]. Building the reflection image fold requires extensive computation to determine the seismic wave travel-times assigned to each geological layer within a crustal section. Determining seismic wave travel-times 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 point-by-point 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 source-sensor 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 well-log 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 well-logs are point-structures driven by porosity variations in fracture-density and, more importantly, by variations in fracture connectivity. Flow-controlling formation permeability fluctuations along spatially erratic poro-connectivity structures within formations do not act coherently to form seismic wave reflection amplitudes. Flow-specific 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 point-by-point velocity distribution as generated for seismic reflection surveys. It is sufficient for TFI to have a comprehensive set of travel-times 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 cross-table TT(i,j).

Volcanic terrains, however, do not support reflection seismic imaging. In consequence, TFI in convective geothermal systems depends on the source-sensor geometry sketched in Figure 10, where travel-times 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 cross-table travel-time 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 low-level 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 flow-related seismic emission in any given time interval Δt. If the voxel is emitting seismic energy in the time interval Δt, the previously-determined cross-table of seismic wave travel-times TT(i,j) between the specified voxel location i and surface seismic sensor locations j=1…Nsns is used to add up the appropriately time-shifted seismic sensor interval data to build an Nsns ~ 1000 strong seismic emissions signal stack. The signal stack is built by summing sensor data for time-interval Δt with each sensor record adjusted according to the arrival times of a voxel-emitted energy given by the cross-table travel-time TT(i,j) from voxel i to sensor j. If during the given time interval Δt, the given voxel emits flow-related seismic energy, TFI processing coherently sums the seismic energy across all surface sensors to give an Nsns-fold 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Δt-fold 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 flow-related 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 travel-time 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 travel-time 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 travel-time 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 triple-stack crustal flow structure count Nfold = Nvox x NΔt x Nsns built up from Nvox spatial connected flow-active 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 fluid-flow-generated seismic emission signals, the Figure 11 25m image spatial resolution accurately locates the fracture specific fluid-active flow structures denoted by the red lines. The time-sequence of map-view flow-structures then shows the black patches of time-evolving seismic emission voxels following strands of the previously mapped fluid-active 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 red-line fracture connectivity structure. Red-line fracture connectivity structures inferred from TFI images prior to wellbore stimulation. Wellbore stimulation validated fracture-interpretation of TFI-mapped red-line 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 Travel-Time Cross-Table 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 travel-time data TT(i,j) exist for crustal volume hosting the convective flow system. As seen in Figure 12, however, travel-time 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 long-settled 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]. Well-log 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 travel-time cross-table TT(i,j) via seismic refraction as sketched in the Figure 13 version of the Figure 10 wellbore-seismic-source and surface seismic sensor array geometry.

Figure 13 shows the refraction seismic elements needed for estimating TT(i,j) cross-table 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 cross-table 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 source-sensor geometry for seismic travel-time data TT(0,j) acquisition from which to estimate the voxel-sensor travel-time cross-table 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 6-8% RMS deviation due to variations of porosity distributed according to the power-law scaling spectral distribution |Χ(k)|² ~ 1/k, 1/km < k < 1/Dm, observed in well-logs 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 field-observed source-to-surface travel-times TT_true(0,j), we proceed to search for generic trial velocity cubes that provide estimated source-to-surface travel-time 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 source-to-sensor 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 power-law-scaling power that scales inversely with spatial frequency k, |X(k)|² ~ 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 travel-time residuals TT_est(i,j) - TT_true(i,j) for the left-hand 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 travel-times j = 1. N_sns for a nominal voxel source i TT_est(i,j) in red closely overlie the “true” travel-times TT_true(i,j) in blue.

Figure 15 Overlay of Nsns surface sensor travel-times from a nominal source voxel i computed from “true” velocity field (blue) and from an estimated velocity field (red). As the red travel-times closely approximate the blue travel-times, TFI-processed Nsns-fold 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 travel-time cross-table TT_true(i,j), six black circles denote a sequence of estimated positions for the source voxel derived from an estimated cross-table TT_est(i,j). The best-fit 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 travel-cross tables TT_est(i,j) to strengthen the statistical certainty of the source location. The Figure 16 caption gives details.

Figure 16 Display of source-voxel search algorithm performance. For the Figure 14 (left) “true” velocity field with source-sensor travel-time cross-table TT_true(i,j), nine random source-voxel locations are selected, i_m, m = 1….9. For each selected source-voxel, 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 travel-time-residual 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 travel-travel cross-tables TT_est(i,j) reliably reduces the statistical source-voxel location error. The statistical accuracy of estimated source-voxel locations is ~1.3 nodes.

4. Discussion

The TFI time-lapse sequence of flow-induced seismic emission imaging shown in Figure 11 illustrates that Nvox sets of Nsns travel-time tables repeated for NΔt time-intervals of TFI listening data can produce 15-25m 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 source-sensor geometry can provide adequate velocity models by which to perform TFI on convective geothermal systems at spatial resolution ~ 50m. Creating flow-structure images that resolve flow heterogeneity v = φv₀ 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₀φv₀ℓ is expressed in terms related to heterogeneous flow images giving v = φv₀ at scale lengths ℓ ~ 50m.

The operational implications for present-day 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 electricity-generating 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₀φv₀ℓ 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 knock-on 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 USD50-100M 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 3-4MWe 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 follow-on USD50-100M 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 high-flow 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 early-stage 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 flow-structure surveys on convective geothermal systems is a small percentage of the potential drilling cost saving. Established TFI data acquisition uses standard stand-alone 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 travel-time data to acquire the TFI source-voxel-to-sensor travel-time cross-table is performed via Matlab-implemented standard Fast Matching Method wave-front 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 flow-structure imaging technology using multi-channel 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 flow-related 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, low-cost 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.