# Fluid Flow and Heat Transport Computation for Power-Law Scaling Poroperm Media

Authors: Peter Leary; Peter Malin; Rami Niemi
Journal: Geofluids
 (2017)

---

## Abstract

In applying Darcy’s law to fluid flow in geologic formations, it is generally assumed that flow variations average to an effectively constant formation flow property. This assumption is, however, fundamentally inaccurate for the ambient crust. Well-log, well-core, and well-flow empirics show that crustal flow spatial variations are systematically correlated from mm to km. Translating crustal flow spatial correlation empirics into numerical form for fluid flow/transport simulation requires computations to be performed on a single global mesh that supports long-range spatial correlation flow structures. Global meshes populated by spatially correlated stochastic poroperm distributions can be processed by 3D finite-element solvers. We model wellbore-logged Dm-scale temperature data due to heat advective flow into a well transecting small faults in a Hm-scale sandstone volume. Wellbore-centric thermal transport is described by Peclet numberPe ≡ a0φv0/D (a0 = wellbore radius, v0 = fluid velocity at a0, φ = mean crustal porosity, and D = rock-water thermal diffusivity). The modelling schema is (i) 3D global mesh for spatially correlated stochastic poropermeability; (ii) ambient percolation flow calibrated by well-core porosity-controlled permeability; (iii) advection via fault-like structures calibrated by well-log neutron porosity; (iv) flow Pe ~ 0.5 in ambient crust and Pe ~ 5 for fault-borne advection.


---

## Body

## 1. Introduction

Numerical models of fluid flow in crustal rock describe a process which cannot be readily observed on any but the smallest scales. The statistical nature of subsurface rock physical property fluctuations is thus a key aspect of crustal flow modelling. Flow and heat transport simulations computed using faulty rock property statistical distributions can lead to faulty conclusions about the unseen fluid flow processes. It is almost certainly the case that the statistical character of rock properties currently attributed to subsurface crustal flow is improperly influenced by the statistics of rock property distributions observed at the surface. In consequence, simulations of subsurface flow and heat transport are often, even typically, poor representations of crustal flow reality.The root observational snag began with the work of Darcy and Dupuit [1, 2]. Darcy and Dupuit jointly concluded that important fluid flow structures in the Paris Basin limestone formations occurred via the infrequent and narrow but laterally persistent clastic formations with flow properties that resembled those of unconsolidated sands used by municipal hydraulic engineers to filter groundwater [3]. In an era in which fluid mechanics was coming of age, fluid flow through unconsolidated sands was recognised through the work of Darcy and Dupuit to be markedly different from fluid flow though pipe-like channels regarded as characteristic of the Paris Basin limestones and chalks. Darcy’s law, for which loss of hydraulic head was directly proportional to flow velocity in unconsolidated filtration column sands and in clastic groundwater formations, was the formal outcome of the Darcy-Dupuit groundwater observations.It was not understood at the time and remained overlooked a century later (e.g., [4–8]) that unconsolidated sands with little or no internal structure are a faulty model for poroperm properties of subsurface clastic rock. An array of present-day well-log, well-core, and well-production evidence indicates that the process of sediment consolidation imprints on sedimentary crustal formations a set of scale-independent random spatial correlations that control formation permeability at all scales [9–16]. Of particular significance is the existence of through-going fluid flow paths at the largest spatial correlation scales. The highly attested spatially correlated nature of crustal rock flow properties violates the central limit theorem, the master tenet of geological thinking, and thus invalidates the standard statistical sampling approach to managing fluid flow at depth in the crust. The fundamental error concerning the statistical nature of crustal fluid flow property distributions remains generally unrecognized (e.g., [17, 18]).The following discussion undertakes to present a straightforward physically motivated computational framework in which proper stochastic procedures represent the spatially correlated distributions of unknown crustal flow properties in a manner which allows finite-element solvers to provide more realistic flow and heat transport simulations. We first review the mathematical and physical origins of the normal distribution used to represent stochastic processes in which significant degrees of spatial or temporal correlation are absent. We then introduce the vast array of well-specific empirical evidence that spatial correlations exist throughout the stochastic processes controlling fluid flow in crustal rock and indicate how the observed spatially correlated flow properties can be represented in numerical meshes amenable to finite-element simulation of fluid flow and heat transport of the crust. Last, we implement our finite-element modelling methodology for a Dm-scale fracture sequence embedded in a Hm-scale crustal volume of tight gas sands. Within this crustal volume a temperature survey in a newly drilled wellbore gave evidence of heat transport by fault-borne fluid flow. Our simulations duplicate the observed wellbore temperature distribution and assign Peclet numbers to the ambient and fault-borne advection flows that are consistent formal analytic solutions for wellbore-centric advective heat flow. We conclude that crustal flow simulation is best performed using a single global numerical mesh in which to embed a suitable scale range of spatially correlated crustal flow structures.


## 2. Spatial Correlation in Crustal Flow Property Distributions

### 2.1. Normal Distributions versus the Statistical Character of Spatially Correlated Crustal Porosity

Crustal reservoir flow modellers can be said to face a problem similar to that which Gauss formally solved in his 1809 treatise on celestial motion [19]. A collection of celestial observations xi,i=1⋯n, was taken to represent a “true” celestial mechanics value x0, measurements which were disturbed by a range of unknown physical processes that generated the observational spread xi. The problem was to determine what estimated value amidst the scattered data xi best approximates the actual physical value x0. In devising what has come to be called the “least squares” solution to the problem of inexact data, Gauss saw that the population distribution of “unbiased” random processes centered on an assumed value x_ formalised by de Moivre in 1738 for “fair” games of chance [20],(1)pxi∣x_=exp⁡-ηxi-x_2,could formally give a clear probability maximum if the unknown parameter x_ is the arithmetic mean of xi, x_=1/n∑xi. Maximizing the probability worked, however, if, but only if, the value of each event xi is independent of all other events. That is, if (1) is assumed to represent the probability pxi∣x_ of event xi given the unknown parameter value x_, then for independent event values xi, the probability P of observation set xi is the product of all event probabilities,(2)P=∏pxi∣x_=exp⁡-η∑xi-x_2.Evaluating (2) for the condition ∂P/∂x_ = 0 maximizes P if x_=1/n∑xi. Normalizing the total probability to unity, the factor η in (1)-(2) is related to the width σ of the normal distribution centered on mean value x_, η = 1/(2σ2).The statistical distribution of independent random events (1) is conveniently summarized by the familiar Gaussian normal distribution controlled by two parameters, the sample mean x_ and the sample standard deviation σ. It follows immediately that the sample values of the mean and standard deviation converge on the actual mean and standard deviation in direct proportion to the square root of the number of samples, √n. Gauss’s least squares solution procedure thus appears to be both compact and complete and as such pervades scientific and practical applications [21].In 1860 Maxwell gave forceful physical testament to the role of event independence by effectively using (1) to establish the statistical physics of ideal gases using molecular velocities as events xi. Application of (1) to the collective microscopic actions leading to macroscopic gas pressure and temperature assumes that each molecule velocity before and after collisions is independent of the velocities of all other molecules. In the words of Maxwell [22], “If experiments on gases are inconsistent with the hypothesis of these propositions, then our theory, though consistent with itself, is proved to be incapable of explaining the phenomena of gases.” In making the reasonable assumption that molecular motion is suitably independent before and after collisions, Maxwell accurately established the micromechanical nature of gas pressure and temperature uniformity throughout a volume of an ideal gas.Traditional crustal reservoir modelling proceeds in a like manner. As with Gauss and Maxwell, the mean and standard deviation of porosity sampleφi are assumed to reliably reflect the mean and standard deviation for the reservoir at large [18]. Following Maxwell’s analysis of an ideal gas, we can put the microscale independence assumption to the test. For a well-log sequence φi of porosity values measured at equal intervals ℓ ~ meter along a crustal wellbore, we ask if the physical events producing porosity sample φi are independent. To answer this test question, we note that a condition for physical property estimates φi to be independent is that the autocorrelation function of the data sequence φi(3)Aj∝∑φiφi+jis zero for all nonzero lags, j≠0. Assuming this property for (3), cosine transformation of an autocorrelation function [23] yields for Aj a constant power-spectrum for crustal sample sequence φi as a function of spatial frequency k,(4)Sk∝∑jAjcos⁡2πjk=A0cos⁡2π·0·k=A0~constant.We thus find that event independence applied to any well-log sequence φi of n values taken at uniform sample length l over length L=nl implies that the power-spectrum of the sequence is constant across the spatial frequency range 1/L<k<1/2l.It is straightforward to determine that well-log power-spectra for crustal porosity, and for crustal rock properties in general, consistently violate spectral condition (4) at all relevant scale lengths [9, 10]. Well-log spectra are observed worldwide to scale inversely with spatial frequency k over five decades of scale length,(5)Sk~1kβ,β~1,1/km<k<1/cm.Heeding Maxwell, we see from empirical observation (5) that crustal porosity sample sets φi are spatially correlated rather than spatial uncorrelated. It follows that, at least formally, sample well-log porosity sequences φi cannot be meaningfully interpreted as unbiased representations of the porosity distribution in the formation at large.


### 2.2. Spatial Correlations in Crustal Permeability

The impact on fluid flow of spatially correlated porosity distributions in crustal volumes emerges from considering well-core poroperm systematics. Well-core poroperm sequences worldwide [11, 15] show that changes in the logarithm of core permeability κi closely track spatial changes in core porosity φi,(6)δφi∝δlog⁡κi,i=1⋯n.Empirical relation (6) is heuristically plausible. If porosity is spatially correlated at all scales (i.e., porosity spatial fluctuations obey spectral scaling relation (5)), then connectivity between pores is plausibly spatially correlated as well. In mathematical terms, if porosity is expressed as a numerical density, φi∝Ni, the ability of pores to link together to produce larger scale permeability is plausibly proportional to the combinatorial factor Ni!~Ni(Ni-1)(Ni-2)(Ni-3)⋯1. Spatial fluctuations in pore density δNi then create changes in pore connectivity δlog⁡(Ni!), giving a simple physical interpretation to empirical relation (6) through Stirling’s formula [24], log⁡N!~N(log⁡N-1).The large-scale effect of spatial correlations in porosity and permeability is seen by integrating well-core-scale spatial fluctuation correlation relation (6) from the scale of sample interval ℓ ~ meter between individual cores samples to the Hm-scale of reservoir formations spanned by well-core sequences. Hm-scale crustal permeability spatial fluctuations can thus be expressed as(7)κi≈κ_exp⁡αφi-φ_,where α is a dimensionless integration constant to be determined empirically.Observed values of the integration constantα are typically 20–40 for crustal reservoir formations with porosity φ ~ 0.1–0.3 [14, 15]. For basement rock with lower porosity φ ~ 0.01, α ~ 300 [16]. As crustal porosity is often quasi-normally distributed over typically small ranges, for example, φ ~ 0.1–0.3, the observed values of α ~ 20–40 mean that (7) is effectively a lognormal distribution for permeability at the scales of crustal well-production. Lognormal well-production distributions consistent with (7) are observed worldwide for conventional oil/gas fields [25], unconventional oil/gas fields [26], geothermal fields [27], basement rock groundwater aquifers [28, 29], and fossil flow aqueous mineral deposits [30–32].The trio of essentially universal crustal property empirics—well-log spatial fluctuation spectral scaling (5), well-core spatial correlation (6) between porosity and the logarithm of permeability, and lognormal distributions of well-production (7)—implies that spatially correlated randomness at all scale lengths conditions fluid flow throughout the crust. Contrary to the commonly accepted prescription that a given geological formation is effectively uniform in flow properties [4–8], the statistical tactic of spatial averaging over poroperm properties cannot properly represent the significant degree of flow heterogeneity inherent in crustal rock. Because of power-law scaling (5), the deviation from a mean or medial permeability background grows with the scale of the flow system. To be relevant to crustal flow heterogeneity, numerical realisations of crustal flow properties must simulate the crust’s inherent spatial correlation and in particular allow for the largest stochastic spatial fluctuations to occur at the largest scales.


### 2.3. Numerical Representation of the Grain-Scale Physical Character of Crustal Permeability

Numerical realisation of the combined crustal poroperm empirics (5)–(7) has a simple grain-scale pictorial basis.While power-law spectral scaling (5) loosely resembles a fractal distribution [33], it is in fact a considerably more powerful physical statement than is implied by fractals. Fractal scaling in a physical system is a statement about population numbers within a system without implications for spatial organization within the system. In contrast, spectral scaling (5) is a physical statement about spatial organization of the crustal poroperm volume elements: in crustal rock everywhere, there is a tendency for crustal volume elements at any given scale to have poroperm properties similar to the poroperm properties of neighbouring crustal volume elements of the same scale [9].Spatial organization such as implied by (5) for crustal rock is observed in a wide range of physical systems [34]. Such systems can be described as undergoing a “critical state order-disorder” phase change. In the crust, (5) implies that an order-disorder phase change occurs throughout the brittle-fracture crustal section lying between the “ordered” ductile lower crust and the “disordered” disaggregated uppermost crust [9, 10]. In the ordered lower crust, porosity is largely absent as fluids tend to be absorbed in hydrated mineral complexes, with deformation occurring by plastic dislocation within an elastic continuum. In the disordered uppermost crust, porosity reaches a form of “critical value” associated with spatial dissociation of grains, with fluids moving freely through the medium [18]. The fluid flow properties of the order-disorder transition state in the brittle-fracture crust are radically distinct from those of lower crustal ductile state and uppermost-crustal disaggregation state.A prominent feature of critical state order-disorder phase transitions is that power-law scaling exponents are independent of the physical nature of microscale elements [34]. With this independence from microscopic physical properties, we can conceptualize rock in terms of a binary physical state that is easily represented numerically. These states are (i) intact grain-grain cement bond contacts forbidding passage of pore fluids and (ii) disrupted grain-grain cement bond contacts permitting passage of pore fluids.In this binary-state perspective, crustal rock critical state order-disorder phase transition properties are equivalent to a percolation lattice for which through-going fluid flow pathways are highly improbable for defect densities below a threshold value and highly probable for defect densities above the threshold value [35]. The effective physical parameter describing rock as a critical state binary population is the fraction N of defective grain-grain contacts relative to the population of intact grain-grain contacts. It is well known for percolation lattices that the transition between nonflowing and flowing lattice states occurs over a narrow range of defect densities centered on a “critical state” threshold defect density Ncrit [35].Regarding crustal rock as a critical state percolation lattice, rock in the ductile lower crustal represents the rock state with defect density below the percolation threshold, and rock in the cohesionless uppermost-crust represents the rock state with defect density above the percolation threshold. The intervening crustal domain of cohesive brittle-fracture rock is perpetually maintained at the critical state percolation densityNcrit by the opposing actions of two on-going crustal processes: (i) damage injection by tectonic finite strain at grain-grain contact that tends to drive the crust to the disordered state of uppermost-crustal disaggregation and (ii) healing/sealing by aqueous chemical disposition at grain-grain contacts that tends to drive the crust to the ordered state of lower crustal ductility.A steady-state percolation system balanced between gaining and losing grain-scale defects is consistent with studies of the fracture state of the crust. Laboratory experiments and numerical modelling identify strain levels of orderε ~ 0.003 and above as creating an elastic damage within the crustal fabric [36]. At the same time, field data assess ambient crustal strain levels in the seismically inactive Fennoscandian crust of order ε ~ 0.003 and above [37]. A binary population picture of crustal rock focused on cement-contacts is consistent with studies of the role cements play in crustal rock mechanics. Macroscopic elastic moduli are often directly proportional to the moduli of the bonding material [38]. Compression tests on a solid skeleton of glass beads demonstrate that small amounts of epoxy resin cement act to localize strain in the cements, thus preventing crushing of the glass beads [39]. Numerical simulation of composite materials shows that the effective elastic moduli are strongly affected by the presence of cement material in the range of 1% to 10% of the composite medium [40].In sum, a broad range of crustal studies is consistent with on-going finite strain damage in rock naturally concentrating in the grain-grain cement bond contact sites that can be associated with crustal poroperm empirics (5)–(7). Numerical representation of the fundamental spatial correlation properties of crustal rock is achieved though (i) allowing numbers between zero and one to represent the range of defect densities within elementary crustal volumes at the nodes of a computational mesh, (ii) generating power-law scaling spatial correlations between those nodal values across the mesh, and (iii) specifying permeability within the mesh by taking the exponent of the nodal effective porosity values scaled by the empirical parameter α appearing in (7). This representation of crustal properties is most effectively achieved for a single global mesh that allows the maximum representation of the crust’s long-range spatial correlation physics at all scales from mm to km. For the purpose of modelling crustal flow and heat transport, crustal poroperm structures controlling fluid flow can be embedded in a single global mesh.


## 2.1. Normal Distributions versus the Statistical Character of Spatially Correlated Crustal Porosity

Crustal reservoir flow modellers can be said to face a problem similar to that which Gauss formally solved in his 1809 treatise on celestial motion [19]. A collection of celestial observations xi,i=1⋯n, was taken to represent a “true” celestial mechanics value x0, measurements which were disturbed by a range of unknown physical processes that generated the observational spread xi. The problem was to determine what estimated value amidst the scattered data xi best approximates the actual physical value x0. In devising what has come to be called the “least squares” solution to the problem of inexact data, Gauss saw that the population distribution of “unbiased” random processes centered on an assumed value x_ formalised by de Moivre in 1738 for “fair” games of chance [20],(1)pxi∣x_=exp⁡-ηxi-x_2,could formally give a clear probability maximum if the unknown parameter x_ is the arithmetic mean of xi, x_=1/n∑xi. Maximizing the probability worked, however, if, but only if, the value of each event xi is independent of all other events. That is, if (1) is assumed to represent the probability pxi∣x_ of event xi given the unknown parameter value x_, then for independent event values xi, the probability P of observation set xi is the product of all event probabilities,(2)P=∏pxi∣x_=exp⁡-η∑xi-x_2.Evaluating (2) for the condition ∂P/∂x_ = 0 maximizes P if x_=1/n∑xi. Normalizing the total probability to unity, the factor η in (1)-(2) is related to the width σ of the normal distribution centered on mean value x_, η = 1/(2σ2).The statistical distribution of independent random events (1) is conveniently summarized by the familiar Gaussian normal distribution controlled by two parameters, the sample mean x_ and the sample standard deviation σ. It follows immediately that the sample values of the mean and standard deviation converge on the actual mean and standard deviation in direct proportion to the square root of the number of samples, √n. Gauss’s least squares solution procedure thus appears to be both compact and complete and as such pervades scientific and practical applications [21].In 1860 Maxwell gave forceful physical testament to the role of event independence by effectively using (1) to establish the statistical physics of ideal gases using molecular velocities as events xi. Application of (1) to the collective microscopic actions leading to macroscopic gas pressure and temperature assumes that each molecule velocity before and after collisions is independent of the velocities of all other molecules. In the words of Maxwell [22], “If experiments on gases are inconsistent with the hypothesis of these propositions, then our theory, though consistent with itself, is proved to be incapable of explaining the phenomena of gases.” In making the reasonable assumption that molecular motion is suitably independent before and after collisions, Maxwell accurately established the micromechanical nature of gas pressure and temperature uniformity throughout a volume of an ideal gas.Traditional crustal reservoir modelling proceeds in a like manner. As with Gauss and Maxwell, the mean and standard deviation of porosity sampleφi are assumed to reliably reflect the mean and standard deviation for the reservoir at large [18]. Following Maxwell’s analysis of an ideal gas, we can put the microscale independence assumption to the test. For a well-log sequence φi of porosity values measured at equal intervals ℓ ~ meter along a crustal wellbore, we ask if the physical events producing porosity sample φi are independent. To answer this test question, we note that a condition for physical property estimates φi to be independent is that the autocorrelation function of the data sequence φi(3)Aj∝∑φiφi+jis zero for all nonzero lags, j≠0. Assuming this property for (3), cosine transformation of an autocorrelation function [23] yields for Aj a constant power-spectrum for crustal sample sequence φi as a function of spatial frequency k,(4)Sk∝∑jAjcos⁡2πjk=A0cos⁡2π·0·k=A0~constant.We thus find that event independence applied to any well-log sequence φi of n values taken at uniform sample length l over length L=nl implies that the power-spectrum of the sequence is constant across the spatial frequency range 1/L<k<1/2l.It is straightforward to determine that well-log power-spectra for crustal porosity, and for crustal rock properties in general, consistently violate spectral condition (4) at all relevant scale lengths [9, 10]. Well-log spectra are observed worldwide to scale inversely with spatial frequency k over five decades of scale length,(5)Sk~1kβ,β~1,1/km<k<1/cm.Heeding Maxwell, we see from empirical observation (5) that crustal porosity sample sets φi are spatially correlated rather than spatial uncorrelated. It follows that, at least formally, sample well-log porosity sequences φi cannot be meaningfully interpreted as unbiased representations of the porosity distribution in the formation at large.


## 2.2. Spatial Correlations in Crustal Permeability

The impact on fluid flow of spatially correlated porosity distributions in crustal volumes emerges from considering well-core poroperm systematics. Well-core poroperm sequences worldwide [11, 15] show that changes in the logarithm of core permeability κi closely track spatial changes in core porosity φi,(6)δφi∝δlog⁡κi,i=1⋯n.Empirical relation (6) is heuristically plausible. If porosity is spatially correlated at all scales (i.e., porosity spatial fluctuations obey spectral scaling relation (5)), then connectivity between pores is plausibly spatially correlated as well. In mathematical terms, if porosity is expressed as a numerical density, φi∝Ni, the ability of pores to link together to produce larger scale permeability is plausibly proportional to the combinatorial factor Ni!~Ni(Ni-1)(Ni-2)(Ni-3)⋯1. Spatial fluctuations in pore density δNi then create changes in pore connectivity δlog⁡(Ni!), giving a simple physical interpretation to empirical relation (6) through Stirling’s formula [24], log⁡N!~N(log⁡N-1).The large-scale effect of spatial correlations in porosity and permeability is seen by integrating well-core-scale spatial fluctuation correlation relation (6) from the scale of sample interval ℓ ~ meter between individual cores samples to the Hm-scale of reservoir formations spanned by well-core sequences. Hm-scale crustal permeability spatial fluctuations can thus be expressed as(7)κi≈κ_exp⁡αφi-φ_,where α is a dimensionless integration constant to be determined empirically.Observed values of the integration constantα are typically 20–40 for crustal reservoir formations with porosity φ ~ 0.1–0.3 [14, 15]. For basement rock with lower porosity φ ~ 0.01, α ~ 300 [16]. As crustal porosity is often quasi-normally distributed over typically small ranges, for example, φ ~ 0.1–0.3, the observed values of α ~ 20–40 mean that (7) is effectively a lognormal distribution for permeability at the scales of crustal well-production. Lognormal well-production distributions consistent with (7) are observed worldwide for conventional oil/gas fields [25], unconventional oil/gas fields [26], geothermal fields [27], basement rock groundwater aquifers [28, 29], and fossil flow aqueous mineral deposits [30–32].The trio of essentially universal crustal property empirics—well-log spatial fluctuation spectral scaling (5), well-core spatial correlation (6) between porosity and the logarithm of permeability, and lognormal distributions of well-production (7)—implies that spatially correlated randomness at all scale lengths conditions fluid flow throughout the crust. Contrary to the commonly accepted prescription that a given geological formation is effectively uniform in flow properties [4–8], the statistical tactic of spatial averaging over poroperm properties cannot properly represent the significant degree of flow heterogeneity inherent in crustal rock. Because of power-law scaling (5), the deviation from a mean or medial permeability background grows with the scale of the flow system. To be relevant to crustal flow heterogeneity, numerical realisations of crustal flow properties must simulate the crust’s inherent spatial correlation and in particular allow for the largest stochastic spatial fluctuations to occur at the largest scales.


## 2.3. Numerical Representation of the Grain-Scale Physical Character of Crustal Permeability

Numerical realisation of the combined crustal poroperm empirics (5)–(7) has a simple grain-scale pictorial basis.While power-law spectral scaling (5) loosely resembles a fractal distribution [33], it is in fact a considerably more powerful physical statement than is implied by fractals. Fractal scaling in a physical system is a statement about population numbers within a system without implications for spatial organization within the system. In contrast, spectral scaling (5) is a physical statement about spatial organization of the crustal poroperm volume elements: in crustal rock everywhere, there is a tendency for crustal volume elements at any given scale to have poroperm properties similar to the poroperm properties of neighbouring crustal volume elements of the same scale [9].Spatial organization such as implied by (5) for crustal rock is observed in a wide range of physical systems [34]. Such systems can be described as undergoing a “critical state order-disorder” phase change. In the crust, (5) implies that an order-disorder phase change occurs throughout the brittle-fracture crustal section lying between the “ordered” ductile lower crust and the “disordered” disaggregated uppermost crust [9, 10]. In the ordered lower crust, porosity is largely absent as fluids tend to be absorbed in hydrated mineral complexes, with deformation occurring by plastic dislocation within an elastic continuum. In the disordered uppermost crust, porosity reaches a form of “critical value” associated with spatial dissociation of grains, with fluids moving freely through the medium [18]. The fluid flow properties of the order-disorder transition state in the brittle-fracture crust are radically distinct from those of lower crustal ductile state and uppermost-crustal disaggregation state.A prominent feature of critical state order-disorder phase transitions is that power-law scaling exponents are independent of the physical nature of microscale elements [34]. With this independence from microscopic physical properties, we can conceptualize rock in terms of a binary physical state that is easily represented numerically. These states are (i) intact grain-grain cement bond contacts forbidding passage of pore fluids and (ii) disrupted grain-grain cement bond contacts permitting passage of pore fluids.In this binary-state perspective, crustal rock critical state order-disorder phase transition properties are equivalent to a percolation lattice for which through-going fluid flow pathways are highly improbable for defect densities below a threshold value and highly probable for defect densities above the threshold value [35]. The effective physical parameter describing rock as a critical state binary population is the fraction N of defective grain-grain contacts relative to the population of intact grain-grain contacts. It is well known for percolation lattices that the transition between nonflowing and flowing lattice states occurs over a narrow range of defect densities centered on a “critical state” threshold defect density Ncrit [35].Regarding crustal rock as a critical state percolation lattice, rock in the ductile lower crustal represents the rock state with defect density below the percolation threshold, and rock in the cohesionless uppermost-crust represents the rock state with defect density above the percolation threshold. The intervening crustal domain of cohesive brittle-fracture rock is perpetually maintained at the critical state percolation densityNcrit by the opposing actions of two on-going crustal processes: (i) damage injection by tectonic finite strain at grain-grain contact that tends to drive the crust to the disordered state of uppermost-crustal disaggregation and (ii) healing/sealing by aqueous chemical disposition at grain-grain contacts that tends to drive the crust to the ordered state of lower crustal ductility.A steady-state percolation system balanced between gaining and losing grain-scale defects is consistent with studies of the fracture state of the crust. Laboratory experiments and numerical modelling identify strain levels of orderε ~ 0.003 and above as creating an elastic damage within the crustal fabric [36]. At the same time, field data assess ambient crustal strain levels in the seismically inactive Fennoscandian crust of order ε ~ 0.003 and above [37]. A binary population picture of crustal rock focused on cement-contacts is consistent with studies of the role cements play in crustal rock mechanics. Macroscopic elastic moduli are often directly proportional to the moduli of the bonding material [38]. Compression tests on a solid skeleton of glass beads demonstrate that small amounts of epoxy resin cement act to localize strain in the cements, thus preventing crushing of the glass beads [39]. Numerical simulation of composite materials shows that the effective elastic moduli are strongly affected by the presence of cement material in the range of 1% to 10% of the composite medium [40].In sum, a broad range of crustal studies is consistent with on-going finite strain damage in rock naturally concentrating in the grain-grain cement bond contact sites that can be associated with crustal poroperm empirics (5)–(7). Numerical representation of the fundamental spatial correlation properties of crustal rock is achieved though (i) allowing numbers between zero and one to represent the range of defect densities within elementary crustal volumes at the nodes of a computational mesh, (ii) generating power-law scaling spatial correlations between those nodal values across the mesh, and (iii) specifying permeability within the mesh by taking the exponent of the nodal effective porosity values scaled by the empirical parameter α appearing in (7). This representation of crustal properties is most effectively achieved for a single global mesh that allows the maximum representation of the crust’s long-range spatial correlation physics at all scales from mm to km. For the purpose of modelling crustal flow and heat transport, crustal poroperm structures controlling fluid flow can be embedded in a single global mesh.


## 3. Numerical Modelling of Fracture-Borne Heat Transport in Crustal Rock

Crustal rock poroperm empirics (5)–(7) interpreted in terms of a generic percolation lattice of flow/no-flow sites are readily represented on a node-based numerical mesh. Finite-element solvers operating on such meshes can then compute Darcy pressure-gradient flow and the fluid transport of heat or solutes. As empirics (5)–(7) focus on long-range spatial correlations supporting critical state percolation backbone flow through a model crustal volume, numerical realisation of the crustal state is most efficiently achieved by making the support mesh encompass the entire model volume. Such global meshes effectively encode the scale-independent spatial correlation connectivity of grain-scale cement bond contact defects over the largest range of scales possible for a given mesh node count. Given that numerical mesh computation cannot match the 20-octave or greater physical process bandwidth of natural rock, feasible computational meshes represent degrees of coarse-graining over the sub-mm to mm scales of rock poropermeability. However, because of the observed five decades of well-log spectral power-law scaling (5) from cm to km, the necessary numerical coarse-graining from sub-mm to cm to dm and above does not significantly compromise numerical representations of large-scale crustal fluid flow processes, provided the spatial correlation structures are allowed to span the entire computational volume in a global mesh.


### 3.1. Illustration of Global Mesh Representation of Crustal Poroperm Structures

Wellbore-based experiments conducted in western Colorado Cretaceous sedimentary tight gas sands provide an opportunity to apply empirics (5)–(7) to modelling to Hm-scale crustal flow and heat transport. The 1995 MWX Multiwell Experiment project [41] drilled a quartet of wellbores to investigate hydrofracture stimulation of tight gas-bearing sands. The potential for active heat transport in the crustal study volume was observed in one of the wellbores as a sequence of thermal gradient spikes in a 200 m section of wellbore at 2.1–2.3 km depth. The thermal gradient spikes coincide with neutron porosity well-log fluctuation spikes in the same well. Interpreting the porosity spikes as evidence for fluid-conductive faults, it is logical that the spatially coincident thermal gradient spikes are due to advective heat entering the wellbore via fracture-borne fluids. The MWX field situation is pictured in Figure 1 as an ambient effective porosity crustal volume with a sequence of embedded fluid-conductive fracture horizons.Figure 1
Effective porosity distribution in 2 Hm-scale simulation of the MWX crustal volume transected by a vertical wellbore. The colour-bar codes the degree of porosity in the section. Cool colours denote the 3D ambient distribution of spatially correlated porosity in accord with spatial fluctuation power-law spectral scaling (5) for well-log porosity range 0.1 < φ < 0.3. Warm colours denote elevated effective porosity with elevated poroconnectivity parameter α controlling permeability via (7). These layers represent crustal faults that intersect the transecting wellbore at depths recording spikes in neutron porosity well-log.Regional erosion relief of order ~850 meters due to the nearby Colorado River makes it plausible that a regional hydraulic head can drive fracture-borne fluids into a newly drilled wellbore on a sufficiently large scale to generate an observable thermal signature. MWX project data for the specific 2 Hm-scale crustal volume intersected by the observation wellbore are shown in Figures2–6. The MWX support data include a well-log porosity sequence φi to fix spatial fluctuation spectral scaling (5), detailed well-log spatial thermal data to fix the system advected heat flow structure, and sizeable amounts of core-scale poroperm data to fix the poroperm-spatial correlation (6) and the poroconnectivity parameter α (7).Figure 2
MWX crustal empirical property (5): wellbore neutron porosity log (a) and porosity log power-spectrum with power-law scaling exponent β ~ 0.98 (b).
(a)
(b)Figure 3
(a) MWX empirical property (6): spatially correlated well-core porosity (blue) and logarithm of well-core permeability (red); poroperm data normalized to zero-mean/unit-variance. (b) MWX empirical property (7): well-core porosity against log(permeability) (blue) with best-fit slope for α ~ 24 (red).
(a)
(b)Figure 4
Neutron porosity wellbore log of MXW well between 1675 and 2520 meters.Figure 5
Thermal gradient wellbore log of MXW well between 1675 and 2520 meters.Figure 6
MWX wellbore temperature profile (blue) overlying wellbore porosity profile (red); the porosity data are exaggerated by a factor 2. Temperature maxima occur at the porosity maxima.In Figure1, a global mesh supports a stochastic numerical realisation of the MWX poroperm parameters and structures needed to simulate flow and/or heat transport of a crustal flow regime:(i)
Systematic spatial correlation of porosity at all scales controlled by power-law scaling exponentβ ~ 1 in (5)(ii)
Control of node-scale permeability (6) by node-scale porosity fluctuations in the range 0.1 < φ < 0.3(iii)
Control of multi-node-scale permeability by poroconnectivity parameter 20 <α < 60 in (7)(iv)
Control of sedimentary flow structures as needed by formation-based changes in porosity and poroconnectivity parameter.The temperature field recorded by a vertical wellbore through the Figure1 crustal volume can be computed by assuming that (i) regional topographic head drives fluids from the crustal volume periphery into the wellbore in proportion to the permeability of the crustal section at each depth along the wellbore and (ii) system heat withdrawal at the wellbore varies along the wellbore with the degree at which fluid flows into the wellbore from the fixed temperature at the volume periphery.With reference to Figure1, this model flow scenario is imposed by setting no-flow boundary conditions on the top and bottom of the crustal volume, setting constant pressure and temperature boundary conditions at the volume periphery, and fixing a constant withdrawal of advected heat from the central wellbore. By assuming a stratified flow approximation consistent with the well-log porosity spike sequence, it is possible to set a wellbore-centric flow Peclet number Pe≡a0φv0/D for radial flow through each section of the model [16], where v0 is the characteristic fluid flow velocity at the wellbore radius a0, φ is flow medium porosity, and D is the characteristic thermal diffusivity for the rock-water system. For low-permeability sections of the crustal volume, fluid flow will be limited to low Peclet numbers implying little advected heat, while for high permeability sections of the crustal volume, that is, the faults, heat flow will achieve higher Peclet numbers.


### 3.2. MWX Crustal Poroperm Parameter Values

Figures2-3 display MWX well-log and well-core data used to control the Figure 1 MWX crustal volume poroperm distributions consistent with empirical poroperm properties (5)–(7). Figure 2 shows the wellbore neutron porosity log φi over wellbore depth interval 1200–2500 m (a) and the log-log plot of porosity fluctuation spectral power as a function of scale length in the spatial frequency range 1 cycle/km < k < 900 cycles/km (b). The observed power-law spectral scaling exponent of (5) is β ~ 1. Figure 3 validates the strong spatial correlation empirics (6) for well-core porosity φ and the logarithm of well-core permeability κ (a) and evaluates the poroconnectivity integration constant of (7), α ~ 24 (b).The observed MWX wellbore porosity and temperature gradient phenomenology is displayed in Figures4-5 for an 800 m crustal interval.Figure6 expands the view of the 200 m thermally active interval between 2.1 and 2.3 km and identifies the close spatial relationship between detrended wellbore temperature (blue) and well-log porosity (red). The coincidence of wellbore temperature maxima and wellbore porosity spikes is flagged by rectangles.The Figure6 spatial correspondence of temperature and porosity maxima implies a simple fluid flow model in which the ambient crustal volume of empirical poroperm properties fixed by Figures 2-3 data hosts a series of fluid flow horizons with poroperm properties consistent with observed weight and location of well-log porosity spikes. It is assumed that drilling the wellbore generated fluid flow paths within an undisturbed steady-state crustal volume and that the wellbore-disturbed fluid flow and heat advection paths correspond to the observed porosity spikes. For simplicity, we assume that the thermal mass of the rock is great enough and the fluid flow disturbance is mild enough that the wellbore fluid is in local thermal equilibrium with the flow/heat advection structure. The observed Figure 6 spatial distribution of wellbore temperature disturbance is then computed on the modelling assumption that heat is extracted from the crustal section at the wellbore and that the resulting temperature profile along the wellbore is due to varying degrees of fluid flow into the central wellbore in response to the steady withdrawal of crustal heat via the wellbore.


### 3.3. Finite-Element Modelling of Fracture-Borne Heat Transport in an Empirically Constrained Stochastic Poroperm Crust

The equation of thermal energy transport arises from combining Darcy’s law of fluid flow velocity proportional to the fluid pressure gradient,(8)v=κμ∇P,with conservation of matter,(9)1B×∂P∂t=∇·v,and conservation of thermal energy, (10)Crρr∂tT=∇·K∇T-CwρwTv.In (8), κ is spatially variable permeability [m2] and μ is constant fluid viscosity [Pa/s]; in (9) B is the essentially constant elastic bulk modulus of the poroperm medium [Pa]; in (10) K is the essentially constant thermal conductivity of rock [W/m·°C] and Crρr and Cwρw [J/m3·°C] are, respectively, the essentially constant volumetric heat capacity of rock and water.Computation of temperature fields (10) subject to thermal boundary conditions and fluid flow within the system requires specifying the system permeability structure κ(x,y,z) as illustrated by the Figure 1 stochastic porosity distribution on a single global mesh. For steady-state fluid flow, that is, setting ∇·v=0, the pressure field is determined by solving the pressure field for given fluid flow boundary conditions,(11)∇·κx,y,z∇Px,y,z=κx,y,z∇·∇Px,y,z+∇κx,y,z·∇Px,y,z=0.The pressure boundary conditions are uniform pressure at the sides of the crustal section, zero-flow across the upper and lower crustal layers, and advected heat loss at the central wellbore. The evolving temperature is then determined by solving for given thermal boundary and initial conditions,(12)∂tTx,y,z,t=∇·D∇Tx,y,z,t-ηTx,y,z,tvx,y,z.In parallel with the pressure boundary conditions, the thermal boundary conditions are uniform temperature at the sides, zero-flow at the upper and lower layers, and a fixed heat flow from the crustal section into the central wellbore. The essentially constant thermal properties for water and rock, rock thermal conductivity K, and volumetric heat capacities Cwρw and Crρr give for (12) constant rock thermal diffusivity D≡K/Crρr~3 W/m·°C/840 J/kg·°C/2400 km/m3 ~ 1.5·10−6 m2/s and constant heat capacity ratio η ≡ Cwρw/Crρr ~ 4280 J/kg·°C/1000 km/m3/(840 J/kg·°C/2400 km/m3) ~ 2.Figure7 shows the computed model thermal curve (red) in relation to the MWX (detrended) wellbore temperature profile from Figure 6 (blue). The model temperature distribution is fixed by siting the four principal porosity spikes at the locations given by the Figure 4 wellbore porosity log, with the porosity-controlled permeability parameter weight associated with each porosity spike consistent with Figure 4 data. The overall magnitude of the model temperature field is matched to the data by adjusting the (unknown) average heat extraction rate at the model wellbore.Figure 7
MWX wellbore temperature profile (red) overlying model temperature profile (blue) computed by finite-elements solution of (11)-(12) on global mesh of Figure 1 stochastic poroperm property distribution.Figure8 shows the fluid velocity field for a uniform pressure external boundary condition draining fluid into the central wellbore via the Figure 1 poroperm distribution. Blue arrows indicate the spatially variable flow in the two “strong” poroperm fault channels, while red arrows indicate spatially variable flow of the ambient crustal flow including flow in the two lesser poroperm fault channels. The time-evolving thermal field (12) is computed for the Figure 8 velocity field, yielding the central wellbore temperature distribution of Figure 7.Figure 8
MWX crustal section fluid flow velocity vector distribution from Figure1 poroperm structure. Blue arrows denote fluid flow in the major fault horizons of Figure 1. Red arrows denote fluid flow in the remaining crustal volume. It can be noted that Figure 1 spatial poroperm heterogeneity leads to azimuthally erratic flow paths within the faults and within the embedding ambient crustal volume. In the present horizontally stratified flow structure, heat transport modelled by spatially averaged flow will probably provide a reasonable wellbore temperature distribution. In the presence of more complex faulting, however, a spatially averaged poroperm medium would provide a poor approximation to the resultant erratic wellbore temperature distribution.Using the Figure8 fluid flow velocity field within the Figure 1 crustal section, the model extracted heat can be compared with the fluid flow velocity field for a mean model flow v0 at a given wellbore-centric radius a0. These model data combine with rock system thermal diffusivity D to estimate the wellbore-centric heat advection Peclet number Pe≡a0φv0/D. Identifying the model Peclet number allows in turn comparison of the wellbore-centric model heat transport system with analytic solutions to wellbore-centric advective flow [16].It can be noted in Figure8 that fluid flow in the poroperm medium is not azimuthally homogenous within the nodal planes of the computational mesh. There is no reason to suppose that such flow homogeneity exists in actual crustal rock, whether in the ambient rock mass itself, or in any faults embedded in the ambient rock mass. The standard modelling tactic is to ignore all such azimuthal flow heterogeneity by assuming spatial averaging over the heterogeneity eliminates any problems flow heterogeneity might cause. In many crustal flow cases, spatial averaging returns adequate modelling results. This is particular likely in low-velocity flow transport of heat because the thermal properties of crustal rock vary a little in space and with rock composition. In the stratified flow structure of Figure 1, spatially averaged flow bringing heat to a central wellbore provides an adequate transport picture. A different situation can, however, be easily imagined if the sequence of faults is taken to be horizontal planar features if Figure 1 were in fact a more spatially complex volume of flow paths. The central wellbore temperature field could in that instance be far more erratic along the wellbore axis and would effectively defeat attempts at modelling through spatially averaged model poroperm distributions.


### 3.4. Wellbore-Centric Flow Heat Advection Peclet NumberPe

Setting up and executing the modelling task posed by applying (11)-(12) to the Figure 1 flow system geometry is simplified by working in a wellbore-centric radial flow approximation in which fluid flow is assumed to be effectively radial within and between the “fault” flow channels designated by the Figure 4 wellbore porosity data. Equation (11) flow in Figure 1 radial sections can be approximated, v(r)~r0v0/r. For radial component divergence operator ∇·A(r)=1/r∂r(rAr), (12) becomes(13)∂tTr,t=D∂r2Tr,t+1-ηPer∂rTr,t,where Pe ≡ r0φv0/D is a wellbore-centric radial flow Peclet number fixing the ratio of advected heat flow r0φv0 to thermal conduction heat flow parameter D. Steady-state advective flow T(r) for advective radial flow between inner radius r0 at temperature T0 and outer radius r1 at temperature T1 has the analytic expression [16](14)Tr=T0+T1-T0r/r0Pe-1r1/r0Pe-1steady-state.The analytic expression for (13) time-evolving temperature field T(r,t) is given by Carslaw & Jaeger [42]. For a line-source heat pulse of energy Q joules per unit length, the time-evolving temperature field is controlled by Peclet number as γ ≡ Pe/2,(15)Tr,t=Q4πKtΓγ+1r24Dtγexp⁡-r24Dt,with Γ(γ+1) the gamma function of argument γ+1 [24]. The solution for Figure 1 advective flow approximated by a line source of radius a at temperature T0 into an infinite medium initially at temperature zero is [42] (16)Tr,t=T0+2π·T0·raγ·∫dkkexp⁡-k2Dt·JγkrYγka-YγkrJγkaJγ2ka+Yγ2ka.Approximations to radial flow temperature distribution (16) computed for the fluid flow velocity field illustrated in Figure 8 give the model and simulation T(r,t) field approximations compared in Figures 9-10 for, respectively, low and high poroperm flow sections of the Figure 1 crustal section.Figure 9
MWX ambient crustal section fluid flow heat advection spatiotemporal profiles for three values of Peclet numberPe = 2γ = a0φv0/D. Blue curves are numerical model temperature fields; red curves are analytic temperature fields (16). Panels are titled by effective Peclet number Pe = 2γ determined by analytic wellbore-centric advective flow temperature field for a sequence of model evolution times. Away from the high porosity sections of the wellbore where fluid flows into the wellbore, the effective Peclet number Pe ~ 0.5 implies advective heat flow on the order of thermal conduction heat flow. The analytic curves are reasonable fits to the numerical curves because the model fluid inflow is more or less uniform along the wellbore axis as assumed for the analytic expression.Figure 10
MWX fault-controlled crustal section fluid flow heat advection spatiotemporal temperature profilesT(r,t) for three values of Peclet number Pe = 2γ = a0φv0/D. Blue curves are numerical model temperature fields; red curves are analytic temperature fields (16). At the high poroperm sections of the wellbore where fluid flows into the wellbore, the effective Peclet number Pe = 2γ ~ 5 implies that the effective fault-borne advective heat flow is of order 10 times the ambient crustal sections above and below the fault intersections. The numerical temperature field deviates from the analytic temperature field at larger radii because the fault-borne fluid radial inflow from the localized fault intersections along the wellbore expands axially along the wellbore where the analytic expression assumes no axial flow along the wellbore.In Figure9, the low poroperm flow section numerical model T(r,t) approximation (blue traces) can be fit to analytic expressions (red traces) for low values of wellbore-centric Peclet number Pe = 2γ = a0φv0/D ~ 0.5. In Figure 10, the high poroperm flow section numerical approximation (blue) is less well matched to the analytic expressions (red) because the numerical model flow is for a narrow axial range of fluid inflow at each fault intersection with the wellbore, while the analytic expression assumes axial symmetry. This disparity aside, the order of magnitude increase in Peclet number Pe = 2γ = a0φv0/D ~ 5 for the high poroperm sections of model flow in blue is in general agreement with the analytic expression in red.


## 3.1. Illustration of Global Mesh Representation of Crustal Poroperm Structures

Wellbore-based experiments conducted in western Colorado Cretaceous sedimentary tight gas sands provide an opportunity to apply empirics (5)–(7) to modelling to Hm-scale crustal flow and heat transport. The 1995 MWX Multiwell Experiment project [41] drilled a quartet of wellbores to investigate hydrofracture stimulation of tight gas-bearing sands. The potential for active heat transport in the crustal study volume was observed in one of the wellbores as a sequence of thermal gradient spikes in a 200 m section of wellbore at 2.1–2.3 km depth. The thermal gradient spikes coincide with neutron porosity well-log fluctuation spikes in the same well. Interpreting the porosity spikes as evidence for fluid-conductive faults, it is logical that the spatially coincident thermal gradient spikes are due to advective heat entering the wellbore via fracture-borne fluids. The MWX field situation is pictured in Figure 1 as an ambient effective porosity crustal volume with a sequence of embedded fluid-conductive fracture horizons.Figure 1
Effective porosity distribution in 2 Hm-scale simulation of the MWX crustal volume transected by a vertical wellbore. The colour-bar codes the degree of porosity in the section. Cool colours denote the 3D ambient distribution of spatially correlated porosity in accord with spatial fluctuation power-law spectral scaling (5) for well-log porosity range 0.1 < φ < 0.3. Warm colours denote elevated effective porosity with elevated poroconnectivity parameter α controlling permeability via (7). These layers represent crustal faults that intersect the transecting wellbore at depths recording spikes in neutron porosity well-log.Regional erosion relief of order ~850 meters due to the nearby Colorado River makes it plausible that a regional hydraulic head can drive fracture-borne fluids into a newly drilled wellbore on a sufficiently large scale to generate an observable thermal signature. MWX project data for the specific 2 Hm-scale crustal volume intersected by the observation wellbore are shown in Figures2–6. The MWX support data include a well-log porosity sequence φi to fix spatial fluctuation spectral scaling (5), detailed well-log spatial thermal data to fix the system advected heat flow structure, and sizeable amounts of core-scale poroperm data to fix the poroperm-spatial correlation (6) and the poroconnectivity parameter α (7).Figure 2
MWX crustal empirical property (5): wellbore neutron porosity log (a) and porosity log power-spectrum with power-law scaling exponent β ~ 0.98 (b).
(a)
(b)Figure 3
(a) MWX empirical property (6): spatially correlated well-core porosity (blue) and logarithm of well-core permeability (red); poroperm data normalized to zero-mean/unit-variance. (b) MWX empirical property (7): well-core porosity against log(permeability) (blue) with best-fit slope for α ~ 24 (red).
(a)
(b)Figure 4
Neutron porosity wellbore log of MXW well between 1675 and 2520 meters.Figure 5
Thermal gradient wellbore log of MXW well between 1675 and 2520 meters.Figure 6
MWX wellbore temperature profile (blue) overlying wellbore porosity profile (red); the porosity data are exaggerated by a factor 2. Temperature maxima occur at the porosity maxima.In Figure1, a global mesh supports a stochastic numerical realisation of the MWX poroperm parameters and structures needed to simulate flow and/or heat transport of a crustal flow regime:(i)
Systematic spatial correlation of porosity at all scales controlled by power-law scaling exponentβ ~ 1 in (5)(ii)
Control of node-scale permeability (6) by node-scale porosity fluctuations in the range 0.1 < φ < 0.3(iii)
Control of multi-node-scale permeability by poroconnectivity parameter 20 <α < 60 in (7)(iv)
Control of sedimentary flow structures as needed by formation-based changes in porosity and poroconnectivity parameter.The temperature field recorded by a vertical wellbore through the Figure1 crustal volume can be computed by assuming that (i) regional topographic head drives fluids from the crustal volume periphery into the wellbore in proportion to the permeability of the crustal section at each depth along the wellbore and (ii) system heat withdrawal at the wellbore varies along the wellbore with the degree at which fluid flows into the wellbore from the fixed temperature at the volume periphery.With reference to Figure1, this model flow scenario is imposed by setting no-flow boundary conditions on the top and bottom of the crustal volume, setting constant pressure and temperature boundary conditions at the volume periphery, and fixing a constant withdrawal of advected heat from the central wellbore. By assuming a stratified flow approximation consistent with the well-log porosity spike sequence, it is possible to set a wellbore-centric flow Peclet number Pe≡a0φv0/D for radial flow through each section of the model [16], where v0 is the characteristic fluid flow velocity at the wellbore radius a0, φ is flow medium porosity, and D is the characteristic thermal diffusivity for the rock-water system. For low-permeability sections of the crustal volume, fluid flow will be limited to low Peclet numbers implying little advected heat, while for high permeability sections of the crustal volume, that is, the faults, heat flow will achieve higher Peclet numbers.


## 3.2. MWX Crustal Poroperm Parameter Values

Figures2-3 display MWX well-log and well-core data used to control the Figure 1 MWX crustal volume poroperm distributions consistent with empirical poroperm properties (5)–(7). Figure 2 shows the wellbore neutron porosity log φi over wellbore depth interval 1200–2500 m (a) and the log-log plot of porosity fluctuation spectral power as a function of scale length in the spatial frequency range 1 cycle/km < k < 900 cycles/km (b). The observed power-law spectral scaling exponent of (5) is β ~ 1. Figure 3 validates the strong spatial correlation empirics (6) for well-core porosity φ and the logarithm of well-core permeability κ (a) and evaluates the poroconnectivity integration constant of (7), α ~ 24 (b).The observed MWX wellbore porosity and temperature gradient phenomenology is displayed in Figures4-5 for an 800 m crustal interval.Figure6 expands the view of the 200 m thermally active interval between 2.1 and 2.3 km and identifies the close spatial relationship between detrended wellbore temperature (blue) and well-log porosity (red). The coincidence of wellbore temperature maxima and wellbore porosity spikes is flagged by rectangles.The Figure6 spatial correspondence of temperature and porosity maxima implies a simple fluid flow model in which the ambient crustal volume of empirical poroperm properties fixed by Figures 2-3 data hosts a series of fluid flow horizons with poroperm properties consistent with observed weight and location of well-log porosity spikes. It is assumed that drilling the wellbore generated fluid flow paths within an undisturbed steady-state crustal volume and that the wellbore-disturbed fluid flow and heat advection paths correspond to the observed porosity spikes. For simplicity, we assume that the thermal mass of the rock is great enough and the fluid flow disturbance is mild enough that the wellbore fluid is in local thermal equilibrium with the flow/heat advection structure. The observed Figure 6 spatial distribution of wellbore temperature disturbance is then computed on the modelling assumption that heat is extracted from the crustal section at the wellbore and that the resulting temperature profile along the wellbore is due to varying degrees of fluid flow into the central wellbore in response to the steady withdrawal of crustal heat via the wellbore.


## 3.3. Finite-Element Modelling of Fracture-Borne Heat Transport in an Empirically Constrained Stochastic Poroperm Crust

The equation of thermal energy transport arises from combining Darcy’s law of fluid flow velocity proportional to the fluid pressure gradient,(8)v=κμ∇P,with conservation of matter,(9)1B×∂P∂t=∇·v,and conservation of thermal energy, (10)Crρr∂tT=∇·K∇T-CwρwTv.In (8), κ is spatially variable permeability [m2] and μ is constant fluid viscosity [Pa/s]; in (9) B is the essentially constant elastic bulk modulus of the poroperm medium [Pa]; in (10) K is the essentially constant thermal conductivity of rock [W/m·°C] and Crρr and Cwρw [J/m3·°C] are, respectively, the essentially constant volumetric heat capacity of rock and water.Computation of temperature fields (10) subject to thermal boundary conditions and fluid flow within the system requires specifying the system permeability structure κ(x,y,z) as illustrated by the Figure 1 stochastic porosity distribution on a single global mesh. For steady-state fluid flow, that is, setting ∇·v=0, the pressure field is determined by solving the pressure field for given fluid flow boundary conditions,(11)∇·κx,y,z∇Px,y,z=κx,y,z∇·∇Px,y,z+∇κx,y,z·∇Px,y,z=0.The pressure boundary conditions are uniform pressure at the sides of the crustal section, zero-flow across the upper and lower crustal layers, and advected heat loss at the central wellbore. The evolving temperature is then determined by solving for given thermal boundary and initial conditions,(12)∂tTx,y,z,t=∇·D∇Tx,y,z,t-ηTx,y,z,tvx,y,z.In parallel with the pressure boundary conditions, the thermal boundary conditions are uniform temperature at the sides, zero-flow at the upper and lower layers, and a fixed heat flow from the crustal section into the central wellbore. The essentially constant thermal properties for water and rock, rock thermal conductivity K, and volumetric heat capacities Cwρw and Crρr give for (12) constant rock thermal diffusivity D≡K/Crρr~3 W/m·°C/840 J/kg·°C/2400 km/m3 ~ 1.5·10−6 m2/s and constant heat capacity ratio η ≡ Cwρw/Crρr ~ 4280 J/kg·°C/1000 km/m3/(840 J/kg·°C/2400 km/m3) ~ 2.Figure7 shows the computed model thermal curve (red) in relation to the MWX (detrended) wellbore temperature profile from Figure 6 (blue). The model temperature distribution is fixed by siting the four principal porosity spikes at the locations given by the Figure 4 wellbore porosity log, with the porosity-controlled permeability parameter weight associated with each porosity spike consistent with Figure 4 data. The overall magnitude of the model temperature field is matched to the data by adjusting the (unknown) average heat extraction rate at the model wellbore.Figure 7
MWX wellbore temperature profile (red) overlying model temperature profile (blue) computed by finite-elements solution of (11)-(12) on global mesh of Figure 1 stochastic poroperm property distribution.Figure8 shows the fluid velocity field for a uniform pressure external boundary condition draining fluid into the central wellbore via the Figure 1 poroperm distribution. Blue arrows indicate the spatially variable flow in the two “strong” poroperm fault channels, while red arrows indicate spatially variable flow of the ambient crustal flow including flow in the two lesser poroperm fault channels. The time-evolving thermal field (12) is computed for the Figure 8 velocity field, yielding the central wellbore temperature distribution of Figure 7.Figure 8
MWX crustal section fluid flow velocity vector distribution from Figure1 poroperm structure. Blue arrows denote fluid flow in the major fault horizons of Figure 1. Red arrows denote fluid flow in the remaining crustal volume. It can be noted that Figure 1 spatial poroperm heterogeneity leads to azimuthally erratic flow paths within the faults and within the embedding ambient crustal volume. In the present horizontally stratified flow structure, heat transport modelled by spatially averaged flow will probably provide a reasonable wellbore temperature distribution. In the presence of more complex faulting, however, a spatially averaged poroperm medium would provide a poor approximation to the resultant erratic wellbore temperature distribution.Using the Figure8 fluid flow velocity field within the Figure 1 crustal section, the model extracted heat can be compared with the fluid flow velocity field for a mean model flow v0 at a given wellbore-centric radius a0. These model data combine with rock system thermal diffusivity D to estimate the wellbore-centric heat advection Peclet number Pe≡a0φv0/D. Identifying the model Peclet number allows in turn comparison of the wellbore-centric model heat transport system with analytic solutions to wellbore-centric advective flow [16].It can be noted in Figure8 that fluid flow in the poroperm medium is not azimuthally homogenous within the nodal planes of the computational mesh. There is no reason to suppose that such flow homogeneity exists in actual crustal rock, whether in the ambient rock mass itself, or in any faults embedded in the ambient rock mass. The standard modelling tactic is to ignore all such azimuthal flow heterogeneity by assuming spatial averaging over the heterogeneity eliminates any problems flow heterogeneity might cause. In many crustal flow cases, spatial averaging returns adequate modelling results. This is particular likely in low-velocity flow transport of heat because the thermal properties of crustal rock vary a little in space and with rock composition. In the stratified flow structure of Figure 1, spatially averaged flow bringing heat to a central wellbore provides an adequate transport picture. A different situation can, however, be easily imagined if the sequence of faults is taken to be horizontal planar features if Figure 1 were in fact a more spatially complex volume of flow paths. The central wellbore temperature field could in that instance be far more erratic along the wellbore axis and would effectively defeat attempts at modelling through spatially averaged model poroperm distributions.


## 3.4. Wellbore-Centric Flow Heat Advection Peclet NumberPe

Setting up and executing the modelling task posed by applying (11)-(12) to the Figure 1 flow system geometry is simplified by working in a wellbore-centric radial flow approximation in which fluid flow is assumed to be effectively radial within and between the “fault” flow channels designated by the Figure 4 wellbore porosity data. Equation (11) flow in Figure 1 radial sections can be approximated, v(r)~r0v0/r. For radial component divergence operator ∇·A(r)=1/r∂r(rAr), (12) becomes(13)∂tTr,t=D∂r2Tr,t+1-ηPer∂rTr,t,where Pe ≡ r0φv0/D is a wellbore-centric radial flow Peclet number fixing the ratio of advected heat flow r0φv0 to thermal conduction heat flow parameter D. Steady-state advective flow T(r) for advective radial flow between inner radius r0 at temperature T0 and outer radius r1 at temperature T1 has the analytic expression [16](14)Tr=T0+T1-T0r/r0Pe-1r1/r0Pe-1steady-state.The analytic expression for (13) time-evolving temperature field T(r,t) is given by Carslaw & Jaeger [42]. For a line-source heat pulse of energy Q joules per unit length, the time-evolving temperature field is controlled by Peclet number as γ ≡ Pe/2,(15)Tr,t=Q4πKtΓγ+1r24Dtγexp⁡-r24Dt,with Γ(γ+1) the gamma function of argument γ+1 [24]. The solution for Figure 1 advective flow approximated by a line source of radius a at temperature T0 into an infinite medium initially at temperature zero is [42] (16)Tr,t=T0+2π·T0·raγ·∫dkkexp⁡-k2Dt·JγkrYγka-YγkrJγkaJγ2ka+Yγ2ka.Approximations to radial flow temperature distribution (16) computed for the fluid flow velocity field illustrated in Figure 8 give the model and simulation T(r,t) field approximations compared in Figures 9-10 for, respectively, low and high poroperm flow sections of the Figure 1 crustal section.Figure 9
MWX ambient crustal section fluid flow heat advection spatiotemporal profiles for three values of Peclet numberPe = 2γ = a0φv0/D. Blue curves are numerical model temperature fields; red curves are analytic temperature fields (16). Panels are titled by effective Peclet number Pe = 2γ determined by analytic wellbore-centric advective flow temperature field for a sequence of model evolution times. Away from the high porosity sections of the wellbore where fluid flows into the wellbore, the effective Peclet number Pe ~ 0.5 implies advective heat flow on the order of thermal conduction heat flow. The analytic curves are reasonable fits to the numerical curves because the model fluid inflow is more or less uniform along the wellbore axis as assumed for the analytic expression.Figure 10
MWX fault-controlled crustal section fluid flow heat advection spatiotemporal temperature profilesT(r,t) for three values of Peclet number Pe = 2γ = a0φv0/D. Blue curves are numerical model temperature fields; red curves are analytic temperature fields (16). At the high poroperm sections of the wellbore where fluid flows into the wellbore, the effective Peclet number Pe = 2γ ~ 5 implies that the effective fault-borne advective heat flow is of order 10 times the ambient crustal sections above and below the fault intersections. The numerical temperature field deviates from the analytic temperature field at larger radii because the fault-borne fluid radial inflow from the localized fault intersections along the wellbore expands axially along the wellbore where the analytic expression assumes no axial flow along the wellbore.In Figure9, the low poroperm flow section numerical model T(r,t) approximation (blue traces) can be fit to analytic expressions (red traces) for low values of wellbore-centric Peclet number Pe = 2γ = a0φv0/D ~ 0.5. In Figure 10, the high poroperm flow section numerical approximation (blue) is less well matched to the analytic expressions (red) because the numerical model flow is for a narrow axial range of fluid inflow at each fault intersection with the wellbore, while the analytic expression assumes axial symmetry. This disparity aside, the order of magnitude increase in Peclet number Pe = 2γ = a0φv0/D ~ 5 for the high poroperm sections of model flow in blue is in general agreement with the analytic expression in red.


## 4. Discussion/Summary

We have interpreted a series of coincident wellbore porosity and temperature spikes in a Hm-scale block of tight gas sandstone as instances of thermal advection in localized fracture-based fluid flow systems intersected by the wellbore. Our modelling of the observed wellbore temperature events demonstrates how finite-element solvers using a single global mesh can use local well-log details to simulate 3D Darcy flow and heat transport in the complex heterogeneous poroperm flow structures that pervade crustal rock.A global mesh can represent both large- and small-scale stochastic spatial fluctuations in ambient crustal flow structure given by empirics (5)–(7). MWX tight gas sand well-log and well-core data validate these stochastic representations via three physical parameters:(i)
Well-log fluctuation power-law spectral scaling exponent,β ~ 1, observed across the cm-km range of scale lengths(ii)
Well-core and well-log porosity 0.1 <φ < 0.3 typical of clastic sedimentary rock, and 0.3 < φ < 0.5, observed for wellbore fracture-site intervals(iii)
Well-core poroconnectivity parameterα ~ 24 typical of ambient sedimentary rock, and α ~ 60 elevated for model localized fault structures in clastic formations.Our approach to flow and transport simulation in spatially correlated crustal poroperm media contrasts with conventional spatial averaging approaches such as effective-media (e.g., [4–8]) and/or fluid flow primarily conducted by (assumed) laminar flow planar fractures (e.g., [43–48]). Conceptually and numerically, the emphasis lies with flow heterogeneity due to long-range spatial correlations rather than flow involving numerical meshes representing flow between more or less uniform poroperm blocks. As illustrated in Figure 1, a single global 3D mesh allows a complete range of flow-physical properties by assigning numerical values to mesh nodes representing flow structures embedded long-range spatially correlated poroperm distributions. Embedded poroperm structures with higher porosities and higher poroconnectivity represent fluid flow and heat transport via macroscopic crustal fault structures. The global mesh numerical mesh solver used here is Sutra [49]. Comparable results for steady-state 3D flow and transport are obtained with the Matlab Partial Differential Equation Toolbox [50].Wellbore-specific crustal flow complexity naturally extends from the MWX crustal fluid flow to flow in low-porosity, low-permeability basement rock. Well-log, well-core, and well-production data validate crustal flow empirics (5)–(7) for basement rock porosities φ ~ .01 and associated poroconnectivity parameters α ~ 300–700 [16]. Application of data-driven modelling control can thus give a perspective on Enhanced/Engineered Geothermal System (EGS) flow stimulation in terms of wellbore-centric Peclet number Pe ≡ a0φv0/D.In MWX crustal rock of ambient mean porosityφ, fluid flow of Peclet number Pe ~ 5 into a wellbore for an observed fault interval ℓ ~ 5 m produces wellbore fluid outtake V = 2πr0φv0l = 2πPelD~10-4 m3/s ~ 0.1 L/s across each 5 m fracture interval. Each of the observed fracture intervals may be inferred to produce ~ 10 m3 of fluid per day. For an MWX 2 km wellbore of volume ~ 60 m3, the rate of fluid flow from one of the major observed MWX fractures displaces the wellbore fluid column in one week. At such low wellbore flow rates, the observable thermal effects of the wellbore disturbance of the ambient steady-state flow/thermal regimes are local to the wellbore-centric fault fluid flow.Using wellbore temperatures to estimate the Peclet number of a natural heat advection system helps to calibrate EGS wellbore flow stimulation. Extrapolating the natural fracture-borne fluid flow associated withPe ~ 5 from a 5 m interval to 1 km wellbore intervals, the total fluid production by the wellbore is 2 liters per second. An EGS stimulation producing 2 liters/second is roughly equivalent to achieving crustal stimulations of 200 MWX crustal fractures of 5 m thickness and 20-meter radius for each km of EGS production wellbore.For EGS stimulation of basement rock, the physical properties to be varied are porosityφ and poroconnectivity parameter α. Porosity is shifted from φ~10% for clastic reservoir rock to φ~1% for basement rock, and poroconnectivity parameter α~24 for clastic reservoir rock is shifted to α~300 for basement rock [16]. Fundamental crustal flow empirics (5)–(7) suggest that crustal flow stimulation naturally proceeds via increased poroconnectivity parameter α rather than increased porosity φ. Noting that increased flow via increased α requires far less work against confining stresses than does increased flow via increased φ, as characteristic of past EGS efforts [43–48], presents a clear implication for future EGS stimulation efforts.



---

*Source: 9687325-2017-10-19.xml*