Geofluids

Volume 2017 (2017), Article ID 9687325, 12 pages

https://doi.org/10.1155/2017/9687325

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

^{1}Advanced Seismic Instrumentation and Research, 1311 Waterside, Dallas, TX 75218-4475, USA^{2}St1 Deep Heat Ltd, Purotie 1, 00381 Helsinki, Finland

Correspondence should be addressed to Peter Leary

Received 22 February 2017; Accepted 20 August 2017; Published 19 October 2017

Academic Editor: Mohamed Fathy El-Amin

Copyright © 2017 Peter Leary et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 number ≡ ( = wellbore radius, = fluid velocity at , = mean crustal porosity, and = 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 ~ 0.5 in ambient crust and ~ 5 for fault-borne advection.

#### 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 , was taken to represent a “true” celestial mechanics value , measurements which were disturbed by a range of unknown physical processes that generated the observational spread . The problem was to determine what estimated value amidst the scattered data best approximates the actual physical value . 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 formalised by de Moivre in 1738 for “fair” games of chance [20],could formally give a clear probability maximum if the unknown parameter is the arithmetic mean of , . Maximizing the probability worked, however, if, but only if, the value of each event is independent of all other events. That is, if (1) is assumed to represent the probability of event given the unknown parameter value , then for independent event values , the probability of observation set is the product of all event probabilities,Evaluating (2) for the condition = 0 maximizes if . Normalizing the total probability to unity, the factor *η* in (1)-(2) is related to the width *σ* of the normal distribution centered on mean value , *η* = .

The statistical distribution of independent random events (1) is conveniently summarized by the familiar Gaussian normal distribution controlled by two parameters, the sample mean 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, . 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 . 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 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 of porosity values measured at equal intervals *ℓ* ~ meter along a crustal wellbore, we ask if the physical events producing porosity sample are independent. To answer this test question, we note that a condition for physical property estimates to be independent is that the autocorrelation function of the data sequence is zero for all nonzero lags, . Assuming this property for (3), cosine transformation of an autocorrelation function [23] yields for a constant power-spectrum for crustal sample sequence as a function of spatial frequency ,We thus find that event independence applied to any well-log sequence of values taken at uniform sample length over length implies that the power-spectrum of the sequence is constant across the spatial frequency range .

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 over five decades of scale length,Heeding Maxwell, we see from empirical observation (5) that crustal porosity sample sets are spatially correlated rather than spatial uncorrelated. It follows that, at least formally, sample well-log porosity sequences 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 closely track spatial changes in core porosity ,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, , the ability of pores to link together to produce larger scale permeability is plausibly proportional to the combinatorial factor . Spatial fluctuations in pore density then create changes in pore connectivity , giving a simple physical interpretation to empirical relation (6) through Stirling’s formula [24], .

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 aswhere *α* 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 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 [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 density 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.