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

Leary, Peter; Malin, Peter; Niemi, Rami

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

Geofluids

On this page

Abstract Introduction Conflicts of Interest Acknowledgments References Copyright Related Articles

Special Issue

Flow and Transport in Porous Media: A Multiscale Focus

View this Special Issue

Research Article | Open Access

Volume 2017 | Article ID 9687325 | https://doi.org/10.1155/2017/9687325

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

Peter Leary,¹Peter Malin,¹and Rami Niemi²

Academic Editor: Mohamed Fathy El-Amin

Received22 Feb 2017

Accepted20 Aug 2017

Published19 Oct 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 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.

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 Figures 2–6. The MWX support data include a well-log porosity sequence 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).

(a)

(b)

(a)

(b)

In Figure 1, 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 Figure 1 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 Figure 1, 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 for radial flow through each section of the model [16], where is the characteristic fluid flow velocity at the wellbore radius , is flow medium porosity, and 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

Figures 2-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 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 < < 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 Figures 4-5 for an 800 m crustal interval.

Figure 6 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 Figure 6 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, with conservation of matter,and conservation of thermal energy, In (8), κ is spatially variable permeability [m²] and μ is constant fluid viscosity [Pa/s]; in (9) is the essentially constant elastic bulk modulus of the poroperm medium [Pa]; in (10) is the essentially constant thermal conductivity of rock [W/m·°C] and and [J/m³·°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 as illustrated by the Figure 1 stochastic porosity distribution on a single global mesh. For steady-state fluid flow, that is, setting , the pressure field is determined by solving the pressure field for given fluid flow boundary conditions,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,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 , and volumetric heat capacities and give for (12) constant rock thermal diffusivity W/m·°C/840 J/kg·°C/2400 km/m³ ~ 1.5·10⁻⁶ m²/s and constant heat capacity ratio η ≡ ~ 4280 J/kg·°C/1000 km/m³/(840 J/kg·°C/2400 km/m³) ~ 2.

Figure 7 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 8 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 Figure 1 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 Figure 8 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 at a given wellbore-centric radius . These model data combine with rock system thermal diffusivity to estimate the wellbore-centric heat advection Peclet number . 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 Figure 8 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 Number

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, . For radial component divergence operator , (12) becomeswhere ≡ is a wellbore-centric radial flow Peclet number fixing the ratio of advected heat flow to thermal conduction heat flow parameter . Steady-state advective flow for advective radial flow between inner radius at temperature and outer radius at temperature has the analytic expression [16]The analytic expression for (13) time-evolving temperature field is given by Carslaw & Jaeger [42]. For a line-source heat pulse of energy joules per unit length, the time-evolving temperature field is controlled by Peclet number as γ ≡ ,with the gamma function of argument [24]. The solution for Figure 1 advective flow approximated by a line source of radius at temperature into an infinite medium initially at temperature zero is [42] Approximations to radial flow temperature distribution (16) computed for the fluid flow velocity field illustrated in Figure 8 give the model and simulation 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 number = 2γ = . Blue curves are numerical model temperature fields; red curves are analytic temperature fields (16). Panels are titled by effective Peclet number = 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 ~ 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 profiles for three values of Peclet number = 2γ = . 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 = 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 Figure 9, the low poroperm flow section numerical model approximation (blue traces) can be fit to analytic expressions (red traces) for low values of wellbore-centric Peclet number = 2γ = ~ 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 = 2γ = ~ 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 ≡ .

In MWX crustal rock of ambient mean porosity φ, fluid flow of Peclet number ~ 5 into a wellbore for an observed fault interval ℓ ~ 5 m produces wellbore fluid outtake = = m³/s ~ 0.1 L/s across each 5 m fracture interval. Each of the observed fracture intervals may be inferred to produce ~ 10 m³ of fluid per day. For an MWX 2 km wellbore of volume ~ 60 m³, 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 with ~ 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 for clastic reservoir rock to for basement rock, and poroconnectivity parameter for clastic reservoir rock is shifted to 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.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The research presented in this paper wholly was funded by St1 Deep Heat Ltd, Purotie 1, 00381 Helsinki, Finland.

References

H. Darcy, Les Fontaines publiques de la ville de Dijon, Victor Dalmont, Paris, France, 1856.
J. Dupuit, Etudes Theoriques et Pratiques sur le Mouvement des Eaux dans les Canaux Decouverts et a Travers les Terrains Permeables, Dunod, Paris, France, 2nd edition, 1863.
G. O. Brown, “Henry Darcy and the making of a law,” Water Resources Research, vol. 38, no. 7, pp. 1–12, 2002.
View at: Publisher Site | Google Scholar
M. Muskat, The Flow of Homogeneous Fluids through Porous Media, McGraw-Hill, 1937.
M. A. Biot, “General theory of three-dimensional consolidation,” Journal of Applied Physics, vol. 12, no. 2, pp. 155–164, 1941.
View at: Publisher Site | Google Scholar
D. R. Horner, “Pressure build-up in wells,” in Third World Petroleum Congress, pp. 503–523, The Hague, 1951.
View at: Google Scholar
M. K. Hubbert, “Darcy's law and the field equations of the flow of underground fluids,” International Association of Scientific Hydrology. Bulletin, vol. 2, no. 1, pp. 23–59, 1957.
View at: Publisher Site | Google Scholar
J. Bear, Dynamics of Fluids in Porous Media, American Elsevier, New York, NY, USA, 1972.
P. C. Leary, “Special section-assessment of schemes for earthquake prediction: Rock as a critical-point system and the inherent implausibility of reliable earthquake prediction,” Geophysical Journal International, vol. 131, no. 3, pp. 451–466, 1997.
View at: Publisher Site | Google Scholar
P. C. Leary, “Fractures and physical heterogeneity in crustal rock,” in Heterogeneity in the Crust and Upper Mantle: Nature, Scaling, and Seismic Properties, J. A. Goff and K . Holliger, Eds., pp. 155–186, Kluwer Academic/Plenum Publishers, New York, NY, USA, 2002.
View at: Google Scholar
P. C. Leary and F. Al-Kindy, “Power-law scaling of spatially correlated porosity and log(permeability) sequences from North-Central North Sea Brae oilfield well core,” Geophysical Journal International, vol. 148, no. 3, pp. 426–442, 2002.
View at: Publisher Site | Google Scholar
P. C. Leary and L. A. Walter, “Crosswell seismic applications to highly heterogeneous tight gas reservoirs,” First Break, vol. 26, no. 3, pp. 49–55, 2008.
View at: Google Scholar
Leary P. C., P. E. Malin, and J. A. Pogacnik, “Computational EGS – heat transport in 1/F-noise fractured media,” in Proceedings of the 37th Workshop on Geothermal Reservoir Engineering, Stanford University, 30 January–1 February, 2012.
View at: Google Scholar
P. Malin, P. Leary, E. Shalev et al., “Flow Lognormality and Spatial Correlation in Crustal Reservoirs: II – Where-to-Drill Guidance via Acoustic/Seismic Imaging,” in Proceedings of the World Geothermal Congress 2015 (WGC '15), Melbourne, Victoria, Australia, 19–24 April 2015.
View at: Google Scholar
J. Pogacnik, P. Leary, P. Malin, P. Geiser, R. Rugis, and B. Valles, “Flow Lognormality and Spatial Correlation in Crustal Reservoirs: III – Natural Permeability Enhancement via Biot Fluid-Rock Coupling at All Scales,” in Proceedings of the World Geothermal Congress 2015 (WGC '15), Melbourne, Victoria, Australia, 19–24 April 2015.
View at: Google Scholar
P. Leary, P. Malin, and R. Niemi, “Finite Element Modelling of Wellbore-Observed Fracture-Borne Heat Advection – Application to EGS Stimulation in Basement Rock,” in Proceedings of the 42nd Workshop on Geothermal Reservoir Engineering, Stanford University, 13–15 February 2017.
View at: Google Scholar
S. Ingebritsen, W. Sanford, and C. Neuzil, Groundwater in Geological Processes, Cambridge University Press, 2006.
G. Mavko, T. Mukerji, and J. Dvorkin, The Rock Physics Handbook, Cambridge University Press, 2009.
D. A. Sprott, “Gauss's contributions to statistics,” Historia Mathematica, vol. 5, no. 2, pp. 183–203, 1978.
View at: Publisher Site | Google Scholar | MathSciNet
H. M. Walker, “De Moivre on the law of normal probability,” in A Source Book in Mathematics, D. E. Smith, Ed., p. 78, Dover, 1985.
View at: Google Scholar
S. M. Stigler, “Gauss and the invention of least squares,” The Annals of Statistics, vol. 9, no. 3, pp. 465–474, 1981.
View at: Publisher Site | Google Scholar | MathSciNet
J. C. Maxwell, “Illustrations of the dynamical theory of gases,” in Kinetic Theory: Selected Readings in Physics, S. G. Brush, Ed., vol. 2, pp. 148–171, Pergamon Press, 1965.
View at: Google Scholar
R. Bracewell, The Fourier Transform and Its Applications, McGraw Hill, 2000.
M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, vol. 55 of National Bureau of Standards Applied Mathematics Series, 1972.
US Energy Information Administration, “Distribution and production of oil and gas wells by state,” 2011, https://www.eia.gov/naturalgas/archive/petrosystem/petrosysog.html.
View at: Google Scholar
IEA, “Special Report: Golden Rules for Golden Age of Gas, International Energy Agency,” 2012, http://www.iea.org/publications/freepublications/publication/WEO2012_GoldenRulesReport.pdf.
View at: Google Scholar
IFC, Success of Geothermal Wells: A Global Study, International Finance Corporation, Wash, USA, 2013.
G. Gustafson, “Strategies for groundwater prospecting in hard rocks: a probabilistic approach,” Norges geologiske undersøkelse Bulletin, vol. 439, pp. 21–25, 2002.
View at: Google Scholar
D. Banks, P. Gundersen, G. Gustafson, J. Mäkelä, and G. Morland, “Regional similarities in the distributions of well yield from crystalline rocks in Fennoscandia,” Norges geologiske undersøkelse Bulletin, vol. 450, pp. 33–47, 2010.
View at: Google Scholar
H. J. de Wijs, “Statistics of ore distribution. Part I: frequency distribution of assay values,” Journal of the Royal Netherlands Geological and Mining Society, vol. 13, pp. 365–375, 1951.
View at: Google Scholar
H. J. de Wijs, “Statistics of ore distribution Part II: theory of binomial distribution applied to sampling and engineering problems,” Journal of the Royal Netherlands Geological and Mining Society, vol. 15, pp. 124-125, 1953.
View at: Google Scholar
D. G. Krige, Lognormal-de Wijsian Geostatistics for Ore Evaluation, South African Institute of Mining and Metallurgy, Johannesburg, South African, 1981.
B. B. Mandelbrot, The Fractal Geometry of Nature, WH Freeman and Co., 1983.
J. J. Binney, N. J. Dowrick, A. J. Fisher, and M. E. J. Newman, The Theory of Critical Phenomena: An Introduction to the Renormalization Group, Oxford University Press, 1995.
D. Stauffer and A. Aharony, Introduction to Percolation Theory, CRC Press, Taylor & Francis, 1994.
Y. Hamiel, O. Katz, V. Lyakhovsky, Z. Reches, and Y. Fialko, “Stable and unstable damage evolution in rocks with implications to fracturing of granite,” Geophysical Journal International, vol. 167, no. 2, pp. 1005–1016, 2006.
View at: Publisher Site | Google Scholar
T. J. Craig, E. Calais, L. Fleitout, L. Bollinger, and O. Scotti, “Evidence for the release of long-term tectonic strain stored in continental interiors through intraplate earthquakes,” Geophysical Research Letters, vol. 43, no. 13, pp. 6826–6836, 2016.
View at: Publisher Site | Google Scholar
B. C. Trent and L. G. Margolin, “A numerical laboratory for granular solids,” Engineering Computations, vol. 9, no. 2, pp. 191–197, 1992.
View at: Publisher Site | Google Scholar
H. Yin and J. Dvorkin, “Strength of cemented grains,” Geophysical Research Letters, vol. 21, no. 10, pp. 903–906, 1994.
View at: Publisher Site | Google Scholar
J. Dvorkin, A. Nur, and H. Yin, “Effective properties of cemented granular materials,” Mechanics of Materials, vol. 18, no. 4, pp. 351–366, 1994.
View at: Publisher Site | Google Scholar
J. C. Lorenz, “MWX - The Multiwell Experiment in the Piceance Basin, Colorado: Reprise from 30 Years Ago,” AAPG Search and Discovery, Article ID 70129, 2012.
View at: Google Scholar
H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, Clarendon Press, Oxford, UK, 1959.
J. E. Warren and P. J. Root, “The Behavior of Naturally Fractured Reservoirs,” Society of Petroleum Engineers Journal, pp. 245–255, 1963.
View at: Publisher Site | Google Scholar
A. C. Gringarten, P. A. Witherspoon, and Y. Ohnishi, “Theory of heat extraction from fractured hot dry rock,” Journal Geophysical Research, vol. 80, pp. 1120–1124, 1975.
View at: Google Scholar
M. C. Cacas, E. Ledoux, G. de Marsily et al., “Modeling fracture flow with a stochastic discrete fracture network: Calibration and validation: 2. The transport model,” Water Resources Research, vol. 26, no. 3, pp. 491–500, 1990.
View at: Publisher Site | Google Scholar
J. W. Tester, B. J. Anderson, A. S. Batchelor et al., The Future of Geothermal Energy: Impact of Enhanced Geothermal Systems (EGS) on the United States in the 21st Century, Massachusetts Institute of Technology, 2006.
K. Pruess, C. Oldenburg, and G. Moridis, “TOUGH2 User's Guide, Version 2.1,” Tech. Rep. LBNL-43134, Lawrence Berkeley National Laboratory, 2012.
View at: Google Scholar
D. Sutter, D. B. Fox, B. J. Anderson, D. L. Koch, P. R. von Rohr, and J. W. Tester, “Sustainable heat farming of geothermal systems: a case study of heat extraction and thermal recovery in a model EGS fractured reservoir,” in Proceedings of the 36th Workshop on Geothermal Reservoir Engineering, Stanford University, 31 January–2 February 2011.
View at: Google Scholar
C. I. Voss and A. M. Provost, “Sutra - A model for saturated-unsaturated, variable-density ground-water flow with solute or energy transport,” Tech. Rep. 02-4231, USGS Water-Resources Investigations, 2010.
View at: Google Scholar
Matlab, Partial Differential Equation Toolbox User's Guide, The MathWorks Inc., Natick, Mass, USA, 2015.

Copyright

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.

PDF Download Citation

Download other formats

Order printed copies

Views

989

Downloads

1076

Citations