Spontaneous formation of fluid escape pipes from subsurface reservoirs

Ubiquitous observations of channelised fluid flow in the form of pipes or chimney-like features in sedimentary sequences provide strong evidence for significant transient permeability-generation in the subsurface. Understanding the mechanisms and dynamics for spontaneous flow localisation into fluid conductive chimneys is vital for natural fluid migration and anthropogenic fluid and gas operations, and in waste sequestration. Yet no model exists that can predict how, when, or where these conduits form. Here we propose a physical mechanism and show that pipes and chimneys can form spontaneously through hydro-mechanical coupling between fluid flow and solid deformation. By resolving both fluid flow and shear deformation of the matrix in three dimensions, we predict fluid flux and matrix stress distribution over time. The pipes constitute efficient fluid pathways with permeability enhancement exceeding three orders of magnitude. We find that in essentially impermeable shale (10−19 m2), vertical fluid migration rates in the high-permeability pipes or chimneys approach rates expected in permeable sandstones (10−15 m2). This previously unidentified fluid focusing mechanism bridges the gap between observations and established conceptual models for overcoming and destroying assumed impermeable barriers. This mechanism therefore has a profound impact on assessing the evolution of leakage pathways in natural gas emissions, for reliable risk assessment for long-term subsurface waste storage, or CO2 sequestration.

Fluid and melt expulsion from porous rocks occurs at many scales on Earth, from the mantle to the shallow subsurface [1][2][3][4][5][6][7][8][9][10] . In the shallow subsurface, distinct features are observed [11][12][13][14][15][16] of vertical chimneys associated with localised continuous or episodic buoyancy-driven fluid flow. These vertical fluid migration pathways are located in most of the studied sedimentary basins. As these regions are of economic interest, they are widely screened with geophysical methods, which provide excellent data that clearly exhibit fluid migration pathways in three dimensions 17 . These features are particularly recognisable by their specific signature on seismic cross-sections (Fig. 1a), while pockmarks and circular-shaped craters (Fig. 1b) are their related expression on the seafloor (Fig. 2a,b). These focussed flow pathways provide an important and efficient transport mechanism for fluid migration, yet the physical controls on their formation are not well understood and still debated [18][19][20] . We propose hydro-mechanical coupling as the dominant underlying physical process for the formation of high-permeability pathways, based on the following observations: (1) Vertical chimneys occur in various sedimentary basin lithologies 12,14,21 and appear unaffected by rock composition, refuting chemical reactions as the main formation mechanism. (2) The vertical chimneys develop through existing formations and are not influenced by sedimentary layering or structural features such as faults 13,21,22 . These observations also suggest that they mainly do not form by reactivation of pre-existing structures, but are the outcome of a self-sustaining dynamic process independent of the inherited geological setting. However, pipes or chimneys may depart from fault planes or root at a similar stratigraphic level. In addition, we do not exclude chemical reactions and hydro-chemo-mechanical coupling to contribute to the generation and propagation of pore fluids [23][24][25] . Since chimney formation cannot be inferred from static models nor geological history, there is a need to develop predictive models to resolve their spontaneous formation and to better constrain the parameters that govern their propagation.
We address this issue by considering a poro-visco-elastic model of rock deformation and flow to quantify the process of channelised fluid migration associated with the deformation of permeable rocks 26 . We use a continuum mechanics model to compute the porous matrix deformation using a non-linear bulk and shear rheology.
Resolving the mechanics is required to calculate the total pressure, strain rates and corresponding stresses, which implies evaluating the entire deviatoric stress tensor. We couple the mechanical solver with a Darcy flow solver and calculate the fluid fluxes using the non-linear Carman-Kozeny relation 27 for dynamic permeability. We streamline a supercomputing approach to resolve the coupled fluid-rock interactions in high resolution 28 . This hydro-mechanical model permits us to predict high-permeability chimney formation and propagation, and the stress distribution in the deforming poro-visco-elastic matrix in three dimensions (Fig. 3). The fully resolved hydro-mechanical coupling generates and propagates solitary waves 29-32 of porosity within a specific parameter range. A pressure sensitive viscous bulk rheology further triggers a significant flow focussing mechanism, 'decompaction weakening' 1 , at the top of the solitons. The usual spherically shaped solitons 32,33 turn into elongated chimney-shaped features 28 (Fig. 3). Solitary waves have been proposed as mechanism for enhanced fluid transport in deep crustal and mantle rocks, as well as for primary hydrocarbon and methane migration from source rocks into sedimentary basins 1,4,25,31,32,34,35 . These deep rock formations are generally considered to exhibit viscous or creep behaviour, in contrast to shallow sedimentary rocks, which are believed to deform in a more  elastic and brittle way. However, recent laboratory experiments on major sedimentary reservoir rock types indicate that time-dependent deformation must be considered 36,37 . The clay content of sediments strongly enhances creep and self-sealing capabilities while hindering the propagation of brittle fractures 37 . Thus, 'soft' clay-rich rocks are regarded as natural barriers 21,37,38 and are recommended as caprock for storage operations. However, although brittle processes are inhibited, a time-dependent creep rheology may lead to localised flow of the sedimentary material in response to applied stress, proportional to bulk viscosity values in the range of 10 16 -10 13 Pa.s 37,39 .

Results
We show the spontaneous development of high-permeability (>3 orders of magnitude over background values) chimneys from a fluid-enriched source region in three dimensions (Fig. 3). The buoyant pore-fluid triggers local decompaction of the porous medium and enables upward-migration within self-organised chimneys. This contrasts with Darcian flow models in non-deforming porous media that predict diffusive fluid flow and spreading of fluids. The focussed flow patterns have a tubular shape (Fig. 3a,b), which can only be resolved by 3-D models. Both the vertical and horizontal 2-D expression (Fig. 3c,d) of the contoured highly permeable regions (Fig. 3b) reproduce natural seismic pipes or chimneys ( Fig. 1) and pockmark observations (Fig. 2). The physical mechanism leading to chimney formation is a natural outcome of time-dependent creep deformation of the fluid-rich porous matrix interacting with the non-linear flow of the pore-fluid. The difference between total pressure and fluid pressure (Fig. 4c) affects the bulk viscosity distribution (Fig. 4b) in a non-linear way and leads to a significant permeability increase (Fig. 4a). The upward-migration of the highly permeable chimneys is sustained by active fluid drainage from the immediate surrounding regions, leading to localised compaction. The fluid flux vectors (white arrows) pointing inward into the pipes support this hypothesis (Fig. 4a,e). The resulting consolidated chimney rim is characterised by decreased permeability and increased viscosity values (Fig. 4a,b,e,f). This phenomenon may indicate that the consolidated rims correspond to the downward-bending horizontal reflectors diagnostic for chimney occurrence in seismic cross-sections ( Fig. 1). The consolidated rim confines pressure deviations ( Fig. 4c) and high fluid fluxes to within the chimney. Thus, pressure measurements outside the chimneys may show no significant perturbations. In contrast to pressure, the localised shear deformation may be detected outside of the chimneys, and the stress envelope runs slightly ahead of the propagating chimney. The second invariant of the deviatoric stress tensor (τ II ) quantifies the magnitude of shear deformation recorded by the porous matrix (Fig. 4d). Variations in stresses may be measurable with seismic methods and their increase in time may further trigger micro-seismic events. The transient fluid expulsion pulse induces irreversible alteration of the permeability and bulk viscosity distributions. These alterations may be responsible for the preservation of the dormant features recognisable in the field (Fig. 2). In contrast, variation in pressure and stresses are restricted to the transient deformation accompanying the fluid pulse.
Recording the fluid flux over time at a specific horizontal cross-section of the model confirms the ability of high-permeability chimneys to enhance vertical flow rates (Fig. 5). Using an initial permeability value of 10 −19 m 2 , representative of typical clay-rich sealing sequences 37 , we predict that the formation of high-permeability chimneys leads to maximal flow rates up to 0.1 m/yr (mean flow rates of about 0.01 m/yr) through a 1 m 2 horizontal cross-section of clay-rich shale (Fig. 5a). Thus, the preferential flow paths (chimneys) within clay-rich rocks enable flow rates of only one order of magnitude below expected diffusive Darcian flow through typical reservoir-quality sandstones with a permeability of 10 −15 m 2 (Fig. 5a). Prior to first chimney arrival, vertical fluid flux values are constant over time and reflect the four orders of magnitude discrepancy in background permeability values. At first chimney breakthrough (Fig. 5b), vertical fluid fluxes in the shale significantly increase and reach their maximal value over an extremely short time (Fig. 5c). Although the flow rates through the chimneys in the shale horizon decreased during the last year of the simulation (Fig. 5d), the highly conductive chimneys continuously enable high fluid fluxes at rates three orders of magnitude higher than the background.

Discussion
Our results show that non-linear hydro-mechanical coupling provides a self-consistent mechanism for rapid and localised fluid expulsion even through a priori impermeable rocks. The chimneys are tubular features generated by transient fluid pulses that exhibit an increased permeability core and a compacted rim, which are preserved in the geological record. They display a characteristic size, spacing and propagation velocity. This stands in contrast to brittle fractures, which are planar features that occur instantaneously. Further, hydro-fractures require high fluid pressure to initiate and remain open, maintaining fluid conductivity. Hydro-fracturing results in a fluid pressure gradient that pushes fluid from the fracture interior into the host rock. Chimneys do not require an over-pressured source region and grow by collecting surrounding fluids into their under-pressured core. This self-sustained fluid collection process results in efficient fluid drainage (the white arrows in Fig. 4a,e), and the associated shear deformation generates a characteristic transient stress pattern.
Our results reproduce the natural observations of fluid migration conduits in many different locations (Figs 1 and 2) and geological settings 12,15,20,22,40 and may have significant implications in other fields of geosciences such as hydrothermal systems 41,42 or gas emissions in different tectonic environments 43 . Our findings suggest that the dominant mechanisms responsible for the spontaneous formation of fluid escape pipes in the subsurface are viscous creep of the porous matrix, decompaction weakening, and hydro-mechanical coupling. Numerical simulations predict that the activated creep leads to the rapid formation of chimneys in low-viscosity shales, expected to be natural flow barriers 13,21 . Thus, resolving the variations in flow patterns, viscous compaction and associated stresses is vital to evaluate storage potential and geological barrier integrity in oil and gas operations 13,44 and waste sequestration 38 . Our results also suggest chimneys and pockmarks to be the expression of rapid migration of natural gas such as methane out of hydrate-rich sediments 10, 15,45 . Understanding the migration mechanism is a  prerequisite for accurately constraining these natural greenhouse gas fluxes towards the atmosphere, with potential implications for long-term climate modification, the evolution of Earth, and for society.

Methods
Summary. We utilised a thermodynamically consistent set of two-phase equations coupling non-linear Darcian flow with a mechanical poro-visco-elastic solver 26 . We adopted a Carman-Kozeny relation 27 to capture the non-linear porosity-dependent permeability using a power-law exponent of 3. We utilised a viscous bulk and shear rheology for the matrix. The bulk viscosity η φ is inversely proportional to the porosity and reduced by a factor R in the regions of fluid overpressure to account for decompaction weakening 1 : The solid shear viscosity µ s is a non-linear function of the strain rate and tends to a reference value for negligible strain rates. We discretised the system of coupled hydro-mechanical partial differential equations using the Finite Difference Method on a regular Cartesian grid in 3-D. We reached an implicit solution of the stencil-based algorithm using an iterative approach in parallel on graphical processing units (GPUs). We utilised a high-resolution numerical grid of 500 × 500 × 1,000 grid points in x, y and z direction, respectively, required to accurately resolve chimney formation and propagation.
Hydro-mechanical coupled model formulation. We employed a system of two-phase equations 26 is the material derivative with respect to the solid. The momentum balance equations for the matrix and the pore-fluid were: where τ ij are the components of the stress deviator, δ ij is the Kronecker delta, g i are the components of the downward pointing gravity acceleration vector, s is the Darcy flux vector (the relative flux of the fluid relative to the solid), µ f is the pore-fluid viscosity and p, p f are the total and fluid pressures, respectively. The total porosity averaged density: s f includes constant solid and fluid densities ρ s and ρ f , respectively. The Carman-Kozeny relation 27 defines the porosity φ dependent permeability κ φ : 0 0 3 where κ 0 is the reference permeability and φ 0 the reference porosity. The strain rate tensor and non-linear viscous creep rheology are expressed as:  ij is the strain rate tensor, τ ij and τ II are the deviatoric stress tensor and the square root of its second invariant, respectively, n is the stress exponent and A is a pre-exponential constant that is equal to the inverse of the solid shear viscosity in the linear viscous case = n 1. The system is closed by a final constitutive equation accounting for viscous (de)compaction: k k s e where η φ is the bulk viscosity. The reference bulk compaction viscosity, η C , at reference porosity φ 0 and  p p f is: where C is the pore geometry dependent coefficient. At other porosity and fluid pressure values, the bulk viscosity η φ is inversely proportional to the porosity and drops with the increase of the effective pressure = − p p p e f to account for decompaction weakening, parametrised by a hyperbolic tangent function in the numerical implementation: where λ p is the sharpness of the transition zone between the decompacting and the compacting regime and R is a rheological constant 1 that quantifies the ratio of compaction η ( ) C over decompaction bulk viscosity. The effective solid shear viscosity µ s is a non-linear function of the strain rate and is implemented as: Simulation details and initial conditions. The hydro-mechanical simulation initial configuration consists of a rectangular box with a dimensionless extent of 30 × 30 × 60 in the x, y (horizontal) and z (depth) directions, respectively. The initial porosity follows an anisotropic Gaussian random-field distribution, with a standard deviation equal to 1, and correlation lengths of 5, 5, 1 in the x, y and z directions, respectively. Further, a high-porosity cylindrical ellipse is located at the first ¼ of the domain height and represents a fluid-rich source  region or reservoir. Permeability values within the source region were nine times higher compared to normalised background values of 1 (Fig. 6). The computational domain was subjected to the downward-pointing gravity field and affected by a background horizontal strike-slip shear-deformation of similar magnitude than buoyancy forces. The pore-fluid is twice as buoyant as the solid. The mechanical problem was solved using free-slip (no shear stress) boundary conditions on all sides of the box. For the fluid flow problem, we applied no flux boundary conditions on all vertical sides of the box, and fixed flux value at the bottom and top boundaries to satisfy the condition = p 0 e (no compaction or decompaction). A total number of 10,000 time steps was necessary to obtain the results. The dimensionless parameters used in the code (Table 1) allowed us to optimally converge the numerical simulation.
Solving strategy. We reached an implicit solution of the coupled set of equations using an iterative pseudo-transient relaxation approach. Acceleration of the residual convergence rates was achieved via a problem-specific damping strategy [46][47][48][49][50] . The numerical algorithm was written in C-CUDA and ran on Nvidia GPUs to efficiently process in parallel 0.25 billion grid points. We relied on Message Passing Interface (MPI) libraries to parallelise the application on distributed memory machines and supercomputers. We calculated the presented high-resolution numerical results in a five-day run on 128 GPUs (Nvidia GTX Titan X Maxwell) in parallel on the octopus supercomputer, in-house designed for such high-performance computations, hosted by the Swiss Geocomputing Centre, Institute of Earth Sciences, University of Lausanne, Switzerland.
Scaling relationships. We used three independent scales: c C c and their dependent combinations such the characteristic velocity δ τ = v / c c c to normalise all the variables, resulting in a dimensionless form of the governing equations. The characteristic length scale δ c is also referred to as the compaction length 51 . τ c is the characteristic time and p c is the characteristic pressure or stress (buoyancy force). Using the scaling relations (equations (13)(14)(15)), the dimensionless model results were scaled to dimensional values representative of reservoir-type rocks: shale, limestone and sandstone. The spacing and the width of the chimneys were mainly controlled by the compaction length δ c . The effective size and propagation speed of high-permeability chimneys ranged from centimetre-sized to metre-sized features in low-permeable shale to features in the hundreds of metres in permeable sandstones. The propagation speed varied from 30 centimetres per year in shale representative for caprock to 1 kilometre per year in permeable sandstones ( Table 2).