Geofluids
Volume 2020 (2020), Article ID 8891801, 14 pages
https://doi.org/10.1155/2020/8891801

Research Article

Interaction between Crustal-Scale Darcy and Hydrofracture Fluid Transport: A Numerical Study

Tamara de Riese,¹ Paul D. Bons,1,2 Enrique Gomez-Rivas,³ and Till Sachau¹

¹Department of Geosciences, Eberhard Karls University Tübingen, Tübingen, Germany
²China University of Geosciences (Beijing), Xueyuan Road 29, Haidian district, 100083 Beijing, China
³Department of Mineralogy, Petrology and Applied Geology, University of Barcelona, Barcelona, Spain

Correspondence should be addressed to Tamara de Riese; tamara.de-riese@uni-tuebingen.de

Received 10 April 2020; Revised 24 June 2020; Accepted 15 October 2020; Published 5 November 2020

Academic Editor: Andrew H. Manning

Copyright © 2020 Tamara de Riese 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.

Crustal-scale fluid flow can be regarded as a bimodal transport mechanism. At low hydraulic head gradients, fluid flow through rock porosity is slow and can be described as diffusional. Structures such as hydraulic breccias and hydrothermal veins both form when fluid velocities and pressures are high, which can be achieved by localized fluid transport in space and time, via hydrofractures. Hydrofracture propagation and simultaneous fluid flow can be regarded as a “ballistic” transport mechanism, which is activated when transport by diffusion alone is insufficient to release the local fluid overpressure. The activation of a ballistic system locally reduces the driving force, through allowing the escape of fluid. We use a numerical model to investigate the properties of the two transport modes in general and the transition between them in particular. We developed a numerical model in order to study patterns that result from bimodal transport. When hydrofractures are activated due to low permeability relative to fluid flux, many hydrofractures form that do not extend through the whole system. These abundant hydrofractures follow a power-law size distribution. A Hurst factor of ~0.9 indicates that the system self-organizes. The abundant small-scale hydrofractures organize the formation of large-scale hydrofractures that ascend through the whole system and drain fluids in large bursts. As the relative contribution of porous flow increases, escaping fluid bursts become less frequent, but more regular in time and larger in volume. We propose that metamorphic rocks with abundant veins, such as in the Kodiak accretionary prism (Alaska) and Otago schists (New Zealand), represent regions with abundant hydrofractures near the fluid source, while hydrothermal breccias are formed by the large fluid bursts that can ascend the crust to shallower levels.

1. Introduction

Fluid flow through rocks and sediments plays a crucial role in many geological and geomechanical processes. They are responsible for a variety of geological phenomena, for example, the formation of veins and hydraulic breccias. Fluids carry dissolved components and heat that advect together with the fluid and are a critical control on many geological processes (e.g., [1]). Understanding the fundamental controls on fluid flow is of primary importance for many applications, e.g., for hydrocarbon migration into reservoirs and during production [2–4], geothermal energy extraction [5, 6], and hydrothermal ore formation and alteration [7]. Classic crustal-scale flow models generally assume continuous flow models that take place by slow fluid percolation through pores, typically driven by topography [8–10], or density differences normally associated with thermal instabilities [11, 12]. However, fluid flow can be dynamic, which, for example, is of major practical relevance as seismicity may be triggered by sudden changes in fluid flow due to fluid injection during fracking or geothermal energy production [2, 13].

Veins and breccias provide abundant evidence in the geological record, that natural fluid flow is not always steady state. Veins are dilatant structures, typically fractures, in which minerals precipitated from fluids [14]. Veins may feature microstructures that reveal that they formed by the “crack-seal mechanism” [15], where a fracture repeatedly opens and is subsequently sealed again due to mineral growth. The process of cracking and sealing can be repeated thousands of times in a single vein. Fractures can be highly conductive structures when open and connected but become barriers to flow when filled with minerals [14, 16, 17]. The crack-seal process thus implies strongly varying fluid fluxes, as well as stress states that must repeatedly reach the failure criterion to induce fracturing. Tensional crack-seal veins are common down to midcrustal levels. Tensional (dilatant) fractures at midcrustal levels can only form when elevated fluid pressure provides a significant contribution to reaching the failure criterion, making such fractures “hydrofractures” [18]. Repeated crack-sealing must thus, at least for a significant part, be due to fluctuating fluid pressures and not only to variations in tectonic stresses. In the diagenetic range (50 to 200°C), fracture development is affected by chemical reactions ([19], on the one hand through chemically assisted fracture growth [20] and on the other hand by dissolution/precipitation reactions, which effect the openness of a fracture. Cements can reduce the local permeability and the large-scale connectivity of fracture systems, while dissolution may increase them [19]. Fracturing can also occur locally due to reaction-driven expansion and has been described, e.g., in granulite facies migmatites [21], during serpentinization [22–24], and during spheroidal weathering [25]. Phase transformations as the replacement of leucite to analcime produce an increase in volume which generates local stresses in the material that can cause fractures, enabling new fluid pathways [26, 27]. However, stress changes due to mineral reactions are slow compared to those due to fluid pressure changes.

Breccias are rock masses composed of broken rock fragments or clasts. Breccias can form by deposition of transported clasts (namely sedimentary breccias [28, 29]) or due to tectonic diminution on faults [30], or as of interest here, hydraulic fracturing by overpressured fluids (namely hydraulic breccias [31, 32]). The Black Forest ore district in SW Germany is a well-studied region where fluid overpressure and dynamic fluid flow has been invoked to explain the formation of mineralised veins and hydrothermal breccias (e.g., [33, 34]). These veins are typically situated close to or at the basement/cover unconformity [33, 35]. The different stages of breccia formation that can be observed in breccias in the Teufelsgrund mine (Münstertal, SW Germany; Figure 1(a)) indicate hydraulic fracturing as the main formation mechanism [31]. This can be inferred from the arrangement of clasts, where some clasts lay next to others in their original cracked position, whereas other clasts are separated from each other. The elevated fluid temperature at which these breccias formed furthermore indicates rapid fluid transport from deeper crustal levels [36]. The Hidden Valley breccia (Figure 1(b)) in South Australia is a huge 10 km² hydrothermal breccia with clasts sizes from tens of microns to hundreds of meters, with mixing of clasts that were originally >km apart vertically [32]. According to Weisheit et al. [32], the formation of this breccia spanned about 200 Myrs with flow velocities up to the order of 1 m/s. As the authors point out, such high fluid flow velocities could not have been maintained for prolonged periods of time, as the total fluid flux would then by far exceed the available fluid budget of about 5-30 km³.

(a)

(b)

Figure 1: Examples of hydrothermal breccias. (a) Teufelsgrund mine, Black Forest, Germany (47°5037N, 7°4914E). Breccia clasts consist of local but altered host rock. The matrix mostly consists of fluorite and barite. (b) Hidden Valley, Mt. Painter inlier, Northern Flinders Ranges, South Australia (30°0525S, 139°3028E). Clasts are a mix of various, heavily altered lithologies, such as metasediments (layered clast on right), basalts (dark clast bottom left), and basement gneisses.

The permeability of rocks is a primary and critical geological parameter because migrating fluids in the Earth’s crust play a major role in mass and heat transfer and for crustal rheology. Permeability has been regarded as a dynamic parameter, as it changes in response to stress, fluid production, and geochemical reactions [37]. In Darcian flow through a porous medium permeability ( in m²) and fluid viscosity (in Pa·s) are the rate parameters that determine fluid fluxes at a given driving fluid pressure gradient ( in Pa/m). The flux () of a fluid through a permeable porous medium is described by Darcy’s law [38], here given in one-dimensional form in the vertical direction (positive downwards).

Permeability varies by >16 orders of magnitude in geological materials (from in crystalline rocks to in well-sorted gravels) [39]. Analyses of geothermal and metamorphic data indicate that the permeability of a tectonically active continental crust decreases logarithmically with depth [40, 41].

The driving pressure is the difference between the equilibrium fluid pressure and the actual fluid pressure [42]. Where fluid pressure is close to hydrostatic variations in surface elevation may result in topography-driven flow, variations in fluid density due to temperature or salinity differences may result in convective flow (e.g., [1]). In topography-driven and convective flow systems, regions exist with downwards fluid flow, for which must be negative. However, this is difficult to reconcile with the high, supralithostatic fluid pressures required to fracture rocks to form veins and breccias [14, 43]. When fluid pressure is close to lithostatic, is positive. Fluids that formed the Hidden Valley Breccia, or the abundant ore-bearing breccias in the German Black Forest, are therefore more likely to have been produced from below these breccias by upward flow. Possible sources are igneous intrusions [44], fluid release due to decompression [33, 36], or mineral dehydration reactions [32, 45].

Transport of clasts over longer distances, as well as structures, indicates that fluid flow rates may be as high as m/s [32, 44, 46]. Fluid ascent rates of 10⁻² to 10⁻¹ m/s have been estimated based on analyses of breccia fragments [47] and veins [48]. Such velocities are difficult to achieve over long time periods and indicate that fluid flow is not continuous but occurs in short bursts [32]. It is well known that high fluid pressures that result in opening of fractures are not constant, but rather intermittent [47, 49–57]. Fluid pressure builds up to exceed the tensional strength of the rock and causes failure, which causes the local permeability to suddenly increase, after which pressure decreases and flow can occur until the fracture permeability is sealed off again. This fluid flow mechanism is not permanent and only of short duration and localised in both time and space. Fluid flow is therefore a highly dynamical process and is able to efficiently transport large amounts of fluid via hydrofractures.

Intermittent flow is predicted to occur when the matrix permeability of a rock is insufficient to accommodate fluid flow, which leads to an increase in fluid pressure and opening of fractures [56, 58]. Fluid pulses released by this process may propagate as batches along so-called mobile hydrofractures [35], which propagate at the upper tip while simultaneously closing at their bottom end when exceeding a critical length [59]. These mobile hydrofractures propagate together with their contained fluid and can potentially reach very high velocities [32, 35, 59, 60] and transport deep fluids and heat up to very shallow crustal levels (e.g., [61]). During propagation, a hydrofracture can interact and merge with other hydrofractures [62], which can lead to mixing of different fluids [36]. Propagation pathways and arrest of mobile hydrofractures are controlled by the local fluid pressure, local material properties, and by the externally applied stress field [60].

Despite the tendency for crustal permeability to decrease with time, long-lived (10³-10⁶ years) hydrothermal systems exist [32, 63, 64], which indicate that processes such as hydrofracturing seem to be an important mechanism [65]. Large-scale crustal permeability adjusts to accommodate rates of internal and external forcing, and in the deeper crust, internal fluid production induced by metamorphism, magmatism, and mantle degassing are responsible [65, 66]. The logarithmic permeability-depth curve therefore reflects a dynamic competition between permeability creation due to fluid sourcing and rock failure, and permeability destruction due to compaction, mineral precipitation, hydrothermal alteration, and retrograde metamorphism [41]. This dynamic interaction can lead to intermittent behaviour and self-organization, which has been proposed for, e.g., the fault-valve behaviour in fractures [51, 52, 67], hydraulic fracturing [68], and in magmatic-hydrothermal ore deposits [69]. Diagenetic assemblages close to normal faults [70] could prove that normal faults can function as valves that intermittently transport fluid as a consequence of enhanced permeability after failure [51]. Permeability has been described as a toggle switch that can take on very high (effectively infinite) and very low (effectively zero) values [58]. Crustal-scale permeability is therefore a dynamically self-adjusting property.

Many transport systems show intermittent behaviour in experiments [71] and in numerical models [56, 58, 72–75]. Scale invariance is a typical feature of intermittent dynamical systems, leading to power-law distributions. Weisheit et al. [32] analyzed the clast size distribution of the Hidden Valley breccia and figured out that a single process was responsible to produce clast sizes over six orders of magnitude. Scale invariance is often found in nature, as in the frequency-size distributions of rock fragments [32, 76], in the size distribution of fractures [77–79], the magnitude-frequency distribution of earthquakes [80], and more.

A large number of studies modelled hydrofracture formation and dynamic fluid flow, using, e.g., finite element [75, 81–84] or particle-lattice models [17, 85–88]. Although hydrofracturing is a complex and highly dynamical phenomenon, the associated transport dynamics can be modelled with a very simple setup using a cellular automaton as introduced by Bak et al. [72], and already adapted to fluid flow by Miller and Nur [58] and Bons and van Milligen [56]. The sandpile model was the first model of a dynamical system displaying self-organized criticality [72]. The idea of this model is that one randomly adds particles to a grid. The pile grows and the slope increases, until it reaches a critical value. Additional particles then lead to avalanches that maintain the critical slope. The special feature of this model is that it evolves into a critical state with no intrinsic time or length scale, without detailed specifications of the initial conditions [72]. Bons and van Milligen [56] used a similar setup to simulate the production, accumulation, and transport of fluids within rocks and observed a self-organized critical type of transport when transport in hydrofractures is activated. Miller and Nur [58] used a crustal-scale cellular automaton model to show that increasing fluid pressure induces local hydrofractures, with hydrofracture size distributions following power laws when the system is at a critical state. Sánchez et al. [73] added a diffusive component to a standard sandpile, identifying that both transport mechanisms interact with each other and that diffusivity can erase the memory of the system. The transport characteristics can be distinguished by identifying the effects on long-range correlations, and by analysing differences in avalanche transport characteristics when diffusive transport is activated [73].

Steady Darcy flow and intermittent fracture flow both seem to be important mechanisms, and these two end-member fluid flow modes have been investigated extensively. There is abundant evidence that such transport systems are neither purely diffusive nor purely intermittent, but surprisingly, little work has been published on the transition between the two regimes (e.g., [71]), especially regarding fluid flow. In this contribution, we assume that fluid transport through the crust can be described as a bimodal transport system. We use a model that is analogous to the sandpile model of Bak et al. [72] in that falling sand grains are replaced by increments in fluid pressure and opening of, and sudden, transient discharge through hydrofractures represent the avalanches. As in Sánchez et al. [73], we add Darcian flow through pores as a second, much slower transport mode. To achieve this, we treat Darcian fluid flow (Eq. (1)) as a diffusive process in which the driving pressure is dissipated by diffusion. We particularly investigate the transition between the very slow (Darcian transport through pores) and very fast (intermittent fracture flow) transport mode.

2. Methods

2.1. Diffusional Pressure Dissipation

We assume a rigid matrix model in which porosity () is constant. This assumption is permissible as fluids are generally at least an order of magnitude more compressible than rocks. Pressure and amount of fluid in the pore space are related through the compressibility (in Pa⁻¹) of the fluid. We can express this with:

is the reference volume (one m³ times the porosity) of fluid at (hydrostatic pressure) and an additional volume of fluid that is added to the pore space. Fluid flux () can now be related to a change in pressure in time (e.g., [89]):

Note that this equation is similar to Fick’s second law for diffusion, with the pressure diffusion coefficient (in m²/s).

Following Fick’s second law, pressure gradients can dissipate in a diffusional fashion with an effective diffusivity (in m²/s). The compressibility of fluids varies with pressure and temperature [90], but for simplicity is here taken as a constant . Porosities in solid rocks also vary from <1 to ~30%. The viscosity of aqueous fluids varies as well, from ca. with increasing temperature [91]. Taking these ranges into account, one obtains a range of about 10⁻⁵ to 10⁻¹ m²/s for in rocks (but can be much larger in unconsolidated coarse clastic sediments). Pressure diffusion is therefore normally much faster than chemical diffusion, except in very low-permeability rocks. Note that all used symbols and parameters are given in Table 1.

Table 1: Symbols.

Symbol
Description
SI unit

Flux

Permeability
m²

Fluid viscosity
Pa·s

Compressibility
Pa⁻¹

Pressure diffusion coefficient
m²/s

Porosity

Solid pressure
Pa

Effective lithostatic pressure

Pa

Volume
m³

Solid density
kg/m³

Fluid density
kg/m³

Gravitational acceleration
m/s²

Time
s

Depth
m

Model variable: Average P over cluster elements

Model variable: Element size

(subscript)
Subscript indicates model parameters

Hurst exponent

2.2. Model Setup

In order to simulate bimodal fluid transport by the two mechanisms, Darcian flow and transport through fractures, we use a square, 2-dimensional orthogonal grid of elements of linear size (Figure 2). We only consider the case that fluid enters the model at the base and is transported to the top. The model is laterally wrapping, i.e., fluid leaving the model on one side enters on the other side. Instead of tracking fluid flow itself, we record the pressure for each element, where is defined as the difference between the actual fluid pressure () and the hydrostatic fluid pressure (): . is set at zero at the top row of elements and increases with 10⁴ Pa/m with depth (assuming a fluid density of 1000 kg/m³ and a gravitational acceleration, , set at 10 m/s²). Fluid flow is implicitly modelled by tracking the evolution of using the pressure diffusion equation (3). Fractures initiate and propagate when a failure criterion is reached. It is assumed that flow through a fracture is fast enough so that pressures within it are equalised effectively instantaneously.

Figure 2: Flow diagram of the numerical simulations. The model consists of a rectangular grid of elements in which the local fluid pressure in excess of the hydrostatic fluid pressure is recorded. Influx of fluid from below is simulated by an incrementally increasing the pressure in the bottom row of the model. The main loop consists of adding fluid and letting the pressures diffuse. Adding pressures is carried out one randomly selected element at a time after which possible hydrofracture initiation is considered. If this occurs, a subloop is initiated for the further propagation of the hydrofracture.

The basic loop in the simulation is as follows:

(1)One element in the bottom row is randomly selected and the pressure in that element is increased by to simulate the influx of fluid from below. Each element in the bottom row is selected once every time step in a random order, which introduces a certain level of random noise in the simulations.(2)After the pressure increase, the possible initiation of a hydrofracture is assessed for every element in the model. A fracture is initiated when the pressure exceeds the effective lithostatic pressure , with the rock density set at 2800 kg/m³, and set to zero at the surface, i.e., in the top row of elements. The hydrofracture is initiated by giving the element the label “broken.”(3)If a fracture is initiated, a propagation loop is started. As only one hydrofracture can exist at any one time in the simulation, all elements with the label “broken” form one connected cluster. First step of the fracture propagation subloop is to equalise the pressure to in all elements within the cluster. is simply the average of the individual pressures of all cluster elements. In the second step, the element on the edge of the cluster with the highest () is selected. If the failure criterion () in that element is reached, one randomly selected direct neighbour element that has not yet failed is added to the cluster. This subloop is repeated until either

(1)none of the elements in the cluster reaches the failure criterion, or(2)the cluster reaches the surface. In the latter case, the pressure in all elements within the cluster is set to zero, which implies that fluid pressure is reduced to hydrostatic and any excess fluid is released at the surface. Once fracture propagation is finished, all elements within the cluster are reset to “unbroken” (implying instantaneous closing or healing of the fracture) and the time and size (number of broken elements) of the cluster is recorded.

(4)Finally, once pressures are increased in all elements in the bottom row (step 2) and any resulting hydrofractures are dealt with (step 3), Darcian flow is simulated using an explicit, forward finite difference scheme. Pressures in the top row are set to zero.

In the current model, fractures completely close or heal after one calculation step. Healing is thus effectively instantaneous relative to the diffusional flow process. Microstructures of crack-seal veins do show that cracks can and do commonly seal faster than fluid pressure builds up. We therefore chose to model the healing as effectively instantaneous, which obviates the need to add a healing-rate parameter. It should, however, be borne in mind that this is an end member and that reducing the healing rate changes the patterns of fractures [17].

Equation (3) shows that the effective pressure diffusion coefficient, , is a function of porosity, fluid viscosity, compressibility, and permeability. is varied in the simulations, which implies a variation in permeability, as the other variables are kept constant. This allows us to explore the transition from hydrofracture (low ) to Darcian-flow dominated (high ) behaviour.

2.3. Scaling

Scaling of the model values (subscript ) to real-world values is determined by the size of individual elements and the pressure increase () at the base due to fluid influx. A unit model pressure gradient in () equals the effective lithostatic gradient ·10⁴ Pa/m. From this, we obtain for the scaling between and the real-world :

We assume a fixed fluid flux , which is typical for metamorphic fluid flux in the crust at a depth of 10-15 km [39]. This influx enters one element along a window of  m² (assuming the model is 1 m thick in the third dimension) and is added to the fluid residing in the  m³ pore space of the element. This results in a pressure increase of:

Assuming a porosity and a compressibility of , we obtain  Pa/s. We can now determine the time step :

In the model, the pressure at the base is raised by every time step. For an element size of  m, the time step is thus (almost 6 days). We can now finally determine the scaling of the model pressure diffusion coefficient () to the real-world diffusivity () and permeability ():

For the given and  m, we obtain . Assuming , the chosen range of modelled -values from 10⁻⁷ to 0.005 corresponds to a permeability range of to . The minimum permeability for which no failure would occur at the given lithostatic gradient and fluid flux is . The simulations therefore almost reach the permeability would never be activated. All used simulation settings are given in Table 2.

Table 2: Simulation parameters and movie names (in the supplementary material (available here)). Values of are used for all figures and are the equivalents in nature.

Permeability
Timesteps []
Movie name

0
0
0

1_D0.avi

0
0
0

10_D0.avi

0
0
0

100_D0.avi

10⁻⁷

10⁻²³

100_D2e-9.avi

10⁻⁴

10⁻²⁰

100_D2e-6.avi

10⁻³

10⁻¹⁹

100_D2e-5.avi

100_D6e-5.avi

10⁻⁴

100_D1e-4.avi

2.4. Analysis

Frequency distributions of hydrofracture sizes have been evaluated. Simulation visualisation and calculation of the time-averaged vertical pressure profile have been done with ImageJ [92]. The time-averaged vertical pressure profile can be calculated via Image-Stacks–Z-project and Analyze-Plot Profile. The effective diffusivity is the slope of the time-averaged vertical pressure profile. Besides this, the rescaled range () statistics was calculated for the pressure fluctuations to quantify the degree of self-similarity. To calculate the rescaled range () of a time series, one calculates , which is the difference between the maximum and minimum accumulated departure from the mean within time span , divided by the standard deviation () [93, 94]. Estimates of are calculated for subsets of time intervals. When plotting these values () versus time lags, the logarithmic slope of this plot gives the Hurst parameter [93]:

Where is a constant and is the time lag. A Hurst parameter of indicates random noise, and a value of indicates a fully correlated time series. It has been shown that a Hurst exponent of to indicates a strong self-similarity of the system [56, 73].

3. Results and Discussion

3.1. Simulation Results

Figure 3(a) shows snapshots of the modelled pressure distribution for ranging from 0 to 10⁻⁴. All the fluid is transported by hydrofractures when the diffusion coefficient is . The top region of the model has a zero overpressure (black in Figure 3(a)). High overpressures only occur near the base (bright colours in the figure). Since matrix flow is zero, fluid overpressures can only propagate upwards in hydrofractures and elements with and must have been reached by hydrofractures that did not reach the surface. Hydrofractures that extend more than about 1/3 of the model tend to reach the surface and drain all their fluid.

Figure 3: Snapshots of simulation results, in the form of maps of the pressure field. Movies of the model can be found in the supplementary material (available here). (a) Results for different diffusion coefficients. (b) Five steps of the simulation with , starting at a time where one hydrofracture reached the top of the model. Note how the initially sharp boundaries become increasingly fuzzy, as fluid flows (i.e., pressure diffuses) into the recently drained regions.

At , corresponding to the other end of the spectrum, the pressure field decreases approximately linearly, signifying that almost all fluid pressure is transported towards the top of the model by diffusional matrix flow. In between these two end members, matrix flow and hydrofracture flow interact with each other, giving rise to a complex behaviour. In general, overpressured fluid gets accumulated and is transported upwards in sudden and fast-flowing bursts within propagating hydrofractures. Hydrofractures that do not reach the top bring fluid up to create regions with nearly lithostatic pressure. Those that do make it to the surface drain all their contained fluid, resulting in low-pressure zones that emanate from the base. At low , the boundaries between the different regions are sharp, while they become fuzzy with increasing . This is illustrated in Figure 3(b), which shows five stages (for ), starting from the point in time where one hydrofracture just crossed the whole model.

The average pressure in the system (Figure 4) fluctuates strongly in the case of pure hydrofracture transport (), even though the input flux at the base of the model is constant. The mean pressure remains low, about a quarter of the maximum mean pressure (for pure diffusion resulting in a linear decrease from bottom to top). With an increasing diffusion coefficient, the pressure fluctuations become more periodical and hydrofracture events become less frequent. With fluctuations become highly periodical. Especially with , one sees an ever repeating regular cycle of pressure buildup, followed by a strong drop in mean pressure as a single hydrofracture drains a large area of the model. The average pressure plot for shows strong similarities with the pressure plot of the fault-valve process by Sibson [51], where a postseismic valving discharge occurs due to enhanced permeability postfailure when the normal fault zone is locally overpressured. In both situations, the pressure over time is highly periodical.

Figure 4: Mean pressure in the whole model as a function of time. No and little diffusion produces an intermittent and irregular behaviour. Higher diffusion coefficients result in more periodical pressure fluctuations.

Whereas Figure 4 shows the mean pressure of the whole model area as a function of time, Figure 5 shows the time-averaged pressure as a function of depth. The pressure decreases from a maximum at the base towards zero at the top. At low diffusion rates (), the average pressure decreases steeply, and at first approximately linearly, to about zero halfway up the model. This is visible in Figure 3 where the top half of the models is black at low . The average pressure-depth curves become smoother and less steep with increasing , trending towards a linear decrease in the case of pure diffusional transport. These results reveal the interesting effect that adding the diffusional transport mode to hydrofracturing does not lead to better draining of the system, but the opposite: on average, more fluid resides inside the model, and fluid pressures are higher at the top than in the absence of diffusion. The two transport modes are not simply additive. The additional diffusion inhibits the formation of hydrofractures, as it reduces the increase of pressure peaks that would lead to hydrofracturing.

(a)

(b)

Figure 5: (a) Time-integrated vertical pressure profile. Red line indicates theoretical graph for pure diffusion. (b) Effective diffusion coefficient with depth.

Figure 5 shows that mean fluid-pressure gradients vary with depth, even though parameters such as failure criterion and pressure diffusion (and hence permeability) are constant throughout the model. If such gradients would be observed in nature, they could be interpreted as the result of variations in permeability [39]. Considering that the flux is constant, one could determine an apparent by dividing the pressure gradient by the flux (based on Eqs. (1) and (3)). At high , the time-averaged pressure decreases steadily from the base of the model to the top, resulting in an apparent diffusion coefficient that decreases moderately with depth (Figure 5(b)). At low , the pressure decreases almost linearly from the base towards about the middle of the model. For this region, one would obtain a constant and high apparent . The upper half of the model has a very low-pressure gradient, which would lead to a very high apparent (Figure 5(b)). This exercise highlights the fact that the activity of hydrofractures affects the pressure gradient and may result in erroneous apparent pressure diffusivity and, hence, apparent permeability. This apparent permeability should not be employed in models that assume Darcian flow only.

The nature of the pressure fluctuations is illustrated with the rescaled-range analysis (; [93, 94]) as is shown in Figure 6(a). The slope of the analysis gives a Hurst exponent of for (Figure 6(b)). For higher diffusion coefficients, the slope remains relatively constant up to , decreases from to , and increases to almost for highest (Figure 6(b)). A Hurst exponent around to indicates some kind of self-organized criticality [56, 73]. The initially decreasing Hurst exponent reveals that the memory of the system stored in local inhomogeneities is erased due to the continuous smoothing driven by diffusion [73]. As soon as diffusion is the dominant transport mechanism, the Hurst exponent increases to , which accounts for a fully correlated signal.

(a)

(b)

Figure 6: Self-similarity analysis of pressure fluctuations with rescaled range analysis. (a) Rescaled range () versus time lag in model time steps. The slope of the curve indicates the Hurst exponent (, plotted in (b)). Dotted lines show and . (b) Hurst exponents versus for shown data-series, including those shown in (a). For , where indicates self-organized criticality [56, 73]. With a higher diffusion coefficient, the Hurst exponent first decreases, until it increases up to almost , indicating a fully correlated time series.

3.2. Statistical Properties of Hydrofracture Size Distributions

Similar models and experiments [56, 72, 73] showed that flow avalanche magnitudes tend to develop power-law frequency distributions. Here, we use the number of elements in a hydrofracture to define its size. The frequency of hydrofractures systematically decreases with their size down to sizes of about 2% of the model size (a few hundred elements in the standard model) (Figure 7). The frequency () distributions of hydrofracture sizes (expressed in number A of elements that form a hydrofracture) that do not reach the surface follow power laws, resulting in a straight line in a log() versus log () plot (Figure 7):

(a)

(b)

Figure 7: Absolute frequency of hydrofracture size distributions (a) for different element sizes () and . Frequencies of hydrofractures that do not reach the surface follow a power law, whereas the distribution for hydrofractures that reach the surface is approximately log-normal. (b) Same as in (a), but for different diffusion coefficients. An element size of was used for all simulations. Power-law fits for hydrofractures not reaching the surface for both (a) and (b) are given in Table 3.

The resulting best-fit values for , the frequency of the smallest possible hydrofracture of one single element, and the power-law exponent () are listed in Table 3. The deviation of the power-law distribution at large sizes is due to finite size effects (limited system size). When comparing the frequency distributions of hydrofracture batch sizes for different element sizes (Figure 7(a)), we observe that the probability to reach a certain batch size becomes higher in a larger system, but the slope does not change with increasing the system size. This simply reflects that there are fewer fractures in a small area than in a large one. The relative frequencies (slope in Figure 7(a)) of fracture sizes should, however, not change, as is the case.

Table 3: Data and power-law fits of data shown in Figure 7.

Volume of fractures reaching the surface (%)
Volume of fractures not reaching the surface (%)
No. of events reaching the surface
No. of events reaching the surface (%)
No. of events not reaching the surface
No. of events not reaching the surface (%)
Power-law fit for fractures not reaching the surface

14.81
85.19
10
0.52
1,906
99.48

2.78
97.22
34
0.19
18,197
99.81

1.87
98.13
401
0.16
253,050
99.84

1.87
98.13
393
0.15
251,940
99.84

1.88
98.12
459
0.17
265,630
99.83

2.93
97.06
447
0.25
179,091
99.75

85.14
14.86
159
18.73
690
81.27

97.99
2.00
15
75
5
25

Hydrofractures that reach the surface are invariably large and their sizes do not follow the aforementioned power-law trends. Fractures of about 10% of the model area are the most frequent. As the hydrofractures that reach the surface are always large, their number relative to the ones that do not reach the surface is small. The fraction of hydrofractures that reach the surface equates the change that one hydrofracture reaches the surface. This fraction is approximately constant at about 0.15% at low (Table 3). The number of hydrofractures per time step also remains approximately constant. This means that at , diffusion essentially plays no role and hydrofractures drain all fluid from the system.

With increasing , the hydrofracture frequency decreases (Figure 7(b)), in particular the frequency of hydrofractures not reaching the surface (Table 3). At the same time, the fraction and size of fractures that reach the surface increases (Figure 7(b) and Table 3). The frequency distribution of escaping hydrofractures reaching the top of the model seems to exhibit log-normal behaviour but is difficult to evaluate due to limited data points.

Hydrofractures that do not reach the top of the model (Figure 7) have a power-law distribution and spread fluid within the crust. These power-law distributions imply that rare but large events transport most of the fluid. Such distributions are often invoked to indicate a self-organized criticality (SOC) in the underlying dynamics of the system [72, 95]. The SOC-type behaviour arises from the existence of a critical slope, which triggers the removal of excess fluid when it is locally exceeded. The activation of a “ballistic” system locally reduces the drive, as the fluid escapes in hydrofractures. The system then (locally) switches back to diffusional transport. The criterion for activation of ballistic dissipation is itself scale independent; pressure has no scale. The activation of the ballistic mode thus happens on all scales at the same time.

The many hydrofractures that do not reach the surface are restricted to the bottom of the model, just above the level where fluid is produced. The expression in the geological record would be abundant crack-seal veins. Abundant veins are common in metamorphic rocks in accretionary prisms, such as the Kodiak accretionary complex, Alaska, USA, [96] or the Otago Schists, New Zealand [97, 98]. In both cases, fluids are sourced by dehydration reactions below, and veins are most abundant (up to 30% of the rock volume) in greenschist-facies rocks. The lower part of the models with abundant hydrofractures may represent these vein-rich zones.

Hydrofractures that do reach the top of the model (Figure 7) are those that transport all the fluid out of the model. This means that a relatively small number of large events actually drain the fluid. Their size distribution does not follow the power-law distributions of hydrofractures not reaching the top of the model. These large events can be interpreted as “dragon-kings” [99]. Such “dragon kings” are outliers which coexist with power-law distributions of event sizes but exceed the power-law extrapolations and are often associated with the occurrence of a bifurcation, catastrophes, or tipping points [99]. They also have been found in, for example, frequency-rupture area statistics of earthquakes [100]. Although not following a power-law distribution themselves, they are intimately related to the smaller hydrofractures that (self-) organize the fluid-pressure distribution and the initiation of large fluid-escape events [101]. In geological reality, these large hydrofractures would represent large fluid expulsion events that would rapidly transport large volumes of fluids from deep to shallow crustal levels where they may form hydrothermal breccias as shown in Figure 1 that are at relatively shallow crustal levels compared to their fluid source. Although fluid pressures would be expected to be close to hydrostatic here, as in the top region of the models, the ascending hydrofractures contain overpressured fluids, which are capable of inducing fracture propagation and even hydraulic brecciation. If fluid ascent is fast enough, as would be expected during hydrofracture propagation [35], fluids will not be able to cool to ambient temperatures during ascent. The resulting fluids thus arrive to the arrest site overpressured and hot (i.e., hydrothermal) and may carry high concentration of dissolved elements, making them prone to deposit ores and cause rock alterations.

The simulations shown here may give some indication on the character of escaping fluid batches. In case of high fluid fluxes (high fluid production rate and/or short duration of the fluid production), the system is expected to be close to the pure hydrofracturing end member. Pressure diffusion would be too low to significantly modify pressure buildup. Fluid escape events would be frequent, but escaping volumes relatively small. The maximum duration of the formation of Jurassic hydrothermal ore deposits in the Black Forest (Figure 1(a)) is ~55 Myrs [33, 102]. However, the actual duration may have been much shorter and related to pulses in basin formation. The Black Forest ore province is characterised by very many, but mostly small and presently uneconomic deposits, consistent with many, but relatively small fluid escape events.

The 10 km² Hidden Valley hydrothermal megabreccia (Figure 1(b)) formed over a period of ~150-200 Myrs, much longer than the Black Forest ore deposits. Lower fluid fluxes would allow more pressure diffusion. The system would be positioned in the transition regime between hydrofracture and diffusion dominated behaviour. The resulting fluid escape events are now expected to be fewer and more regularly spaced in time (Figure 4), but also relatively large. Very large fluid escape events may explain the extreme brecciation, as well as fluidisation and mixing of breccia clasts [32].

4. Conclusions

(i)A numerical model is presented to explore the effect of the relative contributions of Darcian porous flow and flow through hydrofractures on crustal-scale fluid flow. This is achieved through varying the fluid pressure diffusivity, a function of permeability, while keeping the hydrofracture initiation and fluid flux constant(ii)When hydrofracture transport dominates, the system self-organizes. Abundant hydrofractures form at the base of the model and their size-frequency distributions are power laws. The Hurst parameter of ~0.9, calculated from the mean pressure variations over time, supports the development of a self-organized critical state. The hydrofracture size distributions show “dragon-king”-like large hydrofractures that deviate from the power-law distribution. These large hydrofractures actually drain the fluid from the system(iii)With the increasing contribution of Darcian porous flow, pressure fluctuations become larger in magnitude and more regular and cyclical. The transitional regime to Darcian flow is thus characterised by fewer, but larger fluid expulsion events(iv)The observed fluid transport behaviour may explain the abundance of crack-seal veins in metamorphic rocks in, for example, accretionary complexes, as well as the development of hydrothermal hydraulic breccia deposits at shallower crustal levels

Data Availability

All the models and data in this article are based on the equations published here and parameters described in the main text and Tables 1 and 2.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

EGR acknowledges funding by the Spanish Ministry of Science, Innovation and Universities (“Ramón y Cajal” fellowship RYC2018-026335-I and research project PGC2018-093903-B-C22).

Supplementary Materials

Movies of pressure field in the model box over time, for the simulations indicated in Table 2. (Supplementary Materials)

References

G. Chi and C. Xue, “An overview of hydrodynamic studies of mineralization,” Geoscience Frontiers, vol. 2, no. 3, pp. 423–438, 2011.

S. A. Shapiro and C. Dinske, “Fluid-induced seismicity: pressure diffusion and hydraulic fracturing,” Geophysical Prospecting, vol. 57, no. 2, pp. 301–310, 2009.

J. Garland, J. Neilson, S. E. Laubach, and K. J. Whidden, “Advances in carbonate exploration and reservoir analysis,” Geological Society, London, Special Publications, vol. 370, no. 1, pp. 1–15, 2012.

S. M. Agar and G. J. Hampson, “Fundamental controls on flow in carbonates: an introduction,” Petroleum Geoscience, vol. 20, no. 1, pp. 3–5, 2014.

S. A. Shapiro, C. Dinske, C. Langenbruch, and F. Wenzel, “Seismogenic index and magnitude probability of earthquakes induced during reservoir fluid stimulations,” The Leading Edge, vol. 29, no. 3, pp. 304–309, 2010.

D. E. Huenges, Ed., Geothermal Energy Systems: Exploration, Development, and Utilization, Wiley-VCH Verlag, Berlin, 2010.

W. S. Fyfe, N. J. Price, and A. B. Thompson, Fluids in the earth’s Crust, Blsevier, Amsterdam, 1978.

J. Oliver, “Fluids expelled tectonically from orogenic belts: their role in hydrocarbon migration and other geologic phenomena,” Geology, vol. 14, no. 2, pp. 99–102, 1986.

N. H. S. Oliver, J. G. McLellan, B. E. Hobbs, J. S. Cleverley, A. Ord, and L. Feltrin, “100th Anniversary Special Paper: numerical models of extensional deformation, heat transfer, and fluid flow across basement-cover interfaces during basin-related mineralization,” Economic Geology, vol. 101, no. 1, pp. 1–31, 2006.

M. Person, A. Mulch, C. Teyssier, and Y. Gao, “Isotope transport and exchange within metamorphic core complexes,” American Journal of Science, vol. 307, no. 3, pp. 555–589, 2007.

S. K. Matthäi, C. A. Heinrich, and T. Driesner, “Is the Mount Isa copper deposit the product of forced brine convection in the footwall of a major reverse fault?” Geology, vol. 32, pp. 357–360, 2005.

C. Zhao, B. E. Hobbs, and A. Ord, “Convective and advective heat transfer in geological systems,” Geophysics & Geodesy, vol. 229, 2008.

J. L. Rubinstein and A. B. Mahani, “Myths and facts on wastewater injection, hydraulic fracturing, enhanced oil recovery, and induced seismicity,” Seismological Research Letters, vol. 86, no. 4, pp. 1060–1067, 2015.

P. D. Bons, M. A. Elburg, and E. Gomez-Rivas, “A review of the formation of tectonic veins and their microstructures,” Journal of Structural Geology, vol. 43, pp. 33–62, 2012.

J. G. Ramsay, “The crack–seal mechanism of rock deformation,” Nature, vol. 284, no. 5752, pp. 135–139, 1980.

N. H. S. Oliver and P. D. Bons, “Mechanisms of fluid flow and fluid–rock interaction in fossil metamorphic hydrothermal systems inferred from vein–wallrock patterns, geometry and microstructure,” Geofluids, vol. 1, no. 2, pp. 137–162, 2001.

A. Vass, D. Koehn, R. Toussaint, I. Ghani, and S. Piazolo, “The importance of fracture-healing on the deformation of fluid-filled layered systems,” Journal of Structural Geology, vol. 67, pp. 94–106, 2014.

M. A. Etheridge, “Differential stress magnitudes during regional deformation and metamorphism: upper bound imposed by tensile fracturing,” Geology, vol. 11, no. 4, pp. 231–234, 1983.

S. E. Laubach, R. H. Lander, L. J. Criscenti et al., “The role of chemistry in fracture pattern development and opportunities to advance interpretations of geological materials,” Reviews of Geophysics, vol. 57, no. 3, pp. 1065–1111, 2019.

B. K. Atkinson, “Subcritical crack propagation in rocks: theory, experimental results and applications,” Journal of Structural Geology, vol. 4, no. 1, pp. 41–56, 1982.

G. R. Watt, N. H. S. Oliver, and B. J. Griffin, “Evidence for reaction induced micro-fracturing in granulite facies migmatites,” Geology, vol. 28, no. 4, pp. 327–330, 2000.

D. S. O’Hanley, “Serpentinites: records of tectonic and petrological history,” Oxford monographs on geology and geophysics, vol. 34, 1996.

K. Iyer, B. Jamtveit, J. Mathiesen, A. Malthe-Sørenssen, and J. Feder, “Reaction-assisted hierarchical fracturing during serpentinization,” Earth and Planetary Science Letters, vol. 267, no. 3-4, pp. 503–516, 2008.

B. Jamtveit, A. Malthe-Sørenssen, and O. Kostenko, “Reaction enhanced permeability during retrogressive metamorphism,” Earth and Planetary Science Letters, vol. 267, no. 3-4, pp. 620–627, 2008.

R. C. Fletcher, H. L. Buss, and S. L. Brantley, “A spheroidal weathering model coupling porewater chemistry to soil thicknesses during steady-state denudation,” Earth and Planetary Science Letters, vol. 244, no. 1-2, pp. 444–457, 2006.

B. Jamtveit, C. V. Putnis, and A. Malthe-Sørenssen, “Reaction induced fracturing during replacement processes,” Contributions to Mineralogy and Petrology, vol. 157, no. 1, pp. 127–133, 2009.

S. Centrella, A. Putnis, P. Lanari, and H. Austrheim, “Textural and chemical evolution of pyroxene during hydration and deformation: a consequence of retrograde metamorphism,” Lithos, vol. 296, pp. 245–264, 2018.

D. W. Morrow, “Descriptive field classification of sedimentary and diagenetic breccia fabrics in carbonate rocks,” Bulletin of Canadian Petroleum Geology, vol. 30, no. 3, pp. 227–229, 1982.

N. H. Woodcock and K. Mort, “Classification of fault breccias and related fault rocks,” Geological Magazine, vol. 145, no. 3, pp. 435–440, 2008.

N. Keulen, R. Heilbronner, H. Stünitz, A. M. Boullier, and H. Ito, “Grain size distributions of fault rocks: a comparison between experimentally and naturally deformed granitoids,” Journal of Structural Geology, vol. 29, no. 8, pp. 1282–1300, 2007.

M. Jébrak, “Hydrothermal breccias in vein-type ore deposits: a review of mechanisms, morphology and size distribution,” Ore Geology Reviews, vol. 12, no. 3, pp. 111–134, 1997.

A. Weisheit, P. D. Bons, and M. A. Elburg, “Long-lived crustal-scale fluid flow: the hydrothermal mega-breccia of Hidden Valley, Mt. Painter Inlier, South Australia,” International Journal of Earth Sciences, vol. 102, no. 5, pp. 1219–1236, 2013.

S. Staude, P. D. Bons, and G. Markl, “Hydrothermal vein formation by extension-driven dewatering of the middle crust: an example from SW Germany,” Earth and Planetary Science Letters, vol. 286, no. 3-4, pp. 387–395, 2009.

B. F. Walter, M. Burisch, and G. Markl, “Long-term chemical evolution and modification of continental basement brines - a field study from the Schwarzwald, SW Germany,” Geofluids, vol. 16, no. 3, pp. 604–623, 2016.

P. D. Bons, “The formation of large quartz veins by rapid ascent of fluids in mobile hydrofractures,” Tectonophysics, vol. 336, no. 1-4, pp. 1–17, 2001.

P. D. Bons, T. Fusswinkel, E. Gomez-Rivas, G. Markl, T. Wagner, and B. Walter, “Fluid mixing from below in unconformity-related hydrothermal ore deposits,” Geology, vol. 42, no. 12, pp. 1035–1038, 2014.

S. E. Ingebritsen and T. Gleeson, “Crustal permeability: introduction to the special issue,” Geofluids, vol. 15, no. 1-2, pp. 1–10, 2015.

J. Bear, Dynamics of Fluids in Porous Media, Courier Corporation, New York, 2013.

S. E. Ingebritsen and C. E. Manning, “Geological implications of a permeability-depth curve for the continental crust,” Geology, vol. 27, no. 12, pp. 1107–1110, 1999.

C. E. Manning and S. E. Ingebritsen, “Permeability of the continental crust: implications of geothermal data and metamorphic systems,” Reviews of Geophysics, vol. 37, no. 1, pp. 127–150, 1999.

S. E. Ingebritsen and C. E. Manning, “Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism,” Geofluids, vol. 10, no. 1-2, pp. 193–205, 2010.

M. K. Hubbert, “The theory of ground-water motion,” The Journal of Geology, vol. 48, 8, Part 1, pp. 785–944, 1940.

S. F. Cox, “The application of failure mode diagrams for exploring the roles of fluid pressure and stress states in controlling styles of fracture-controlled permeability enhancement in faults and shear zones,” Geofluids, vol. 10, no. 1-2, pp. 217–233, 2010.

N. H. S. Oliver, M. J. Rubenach, B. Fu et al., “Granite-related overpressure and volatile release in the mid crust: fluidized breccias from the Cloncurry District, Australia,” Geofluids, vol. 6, no. 4, pp. 346–358, 2006.

J. Vry, R. Powell, K. M. Golden, and K. Petersen, “The role of exhumation in metamorphic dehydration and fluid production,” Nature Geoscience, vol. 3, no. 1, pp. 31–35, 2010.

G. Chi, C. Xue, H. Qing, W. Xue, J. Zhang, and Y. Sun, “Hydrodynamic analysis of clastic injection and hydraulic fracturing structures in the Jinding Zn-Pb deposit, Yunnan, China,” Geoscience Frontiers, vol. 3, no. 1, pp. 73–84, 2012.

P. Eichhubl and J. R. Boles, “Rates of fluid flow in fault systems; evidence for episodic rapid fluid flow in the Miocene Monterey Formation, coastal California,” American Journal of Science, vol. 300, no. 7, pp. 571–600, 2000.

A. Okamoto and N. Tsuchiya, “Velocity of vertical fluid ascent within vein-forming fractures,” Geology, vol. 37, no. 6, pp. 563–566, 2009.

R. H. Sibson, J. M. M. Moore, and A. H. Rankin, “Seismic pumping—a hydrothermal fluid transport mechanism,” Journal of the Geological Society, vol. 131, no. 6, pp. 653–659, 1975.

R. H. Sibson, “Crustal stress, faulting and fluid flow,” Geological Society, London, Special Publications, vol. 78, no. 1, pp. 69–84, 1994.

R. H. Sibson, “Fluid involvement in normal faulting,” Journal of Geodynamics, vol. 29, no. 3-5, pp. 469–499, 2000.

R. H. Sibson, “Tectonic controls on maximum sustainable overpressure: fluid redistribution from stress transitions,” Journal of Geochemical Exploration, vol. 69, pp. 471–475, 2000.

A. Nur and J. Walder, “Hydraulic pulses in the Earth’s crust,” in International Geophysics, vol. 51, pp. 461–473, Academic Press, 1992.

A. M. Boullier, B. Charoy, and P. J. Pollard, “Fluctuation in porosity and fluid pressure during hydrothermal events: textural evidence in the Emuford District, Australia,” Journal of Structural Geology, vol. 16, no. 10, pp. 1417–1429, 1994.

S. F. Cox, “Faulting processes at high fluid pressures: an example of fault valve behavior from the Wattle Gully Fault, Victoria, Australia,” Journal of Geophysical Research: Solid Earth, vol. 100, no. B7, pp. 12841–12859, 1995.

P. D. Bons and B. P. van Milligen, “New experiment to model self-organized critical transport and accumulation of melt and hydrocarbons from their source rocks,” Geology, vol. 29, no. 10, pp. 919–922, 2001.

E. Gomez-Rivas, P. D. Bons, D. Koehn et al., “The Jabal Akhdar Dome in the Oman Mountains: evolution of a dynamic fracture system,” American Journal of Science, vol. 314, no. 7, pp. 1104–1139, 2014.

S. A. Miller and A. Nur, “Permeability as a toggle switch in fluid-controlled crustal processes,” Earth and Planetary Science Letters, vol. 183, no. 1-2, pp. 133–146, 2000.

J. Weertman, “Theory of water-filled crevasses in glaciers applied to vertical magma transport beneath oceanic ridges,” Journal of Geophysical Research, vol. 76, no. 5, pp. 1171–1183, 1971.

T. Dahm, “On the shape and velocity of fluid-filled fractures in the Earth,” Geophysical Journal International, vol. 142, no. 1, pp. 181–192, 2000.

J. D. Martín-Martín, E. Gomez-Rivas, D. Gómez-Gras et al., “Activation of stylolites as conduits for overpressured fluid flow in dolomitized platform carbonates,” Geological Society, London, Special Publications, vol. 459, no. 1, pp. 157–176, 2018.

G. Ito and S. J. Martel, “Focusing of magma in the upper mantle through dike interaction,” Journal of Geophysical Research: Solid Earth, vol. 107, no. B10, pp. ECV 6-1–ECV 6-17, 2002.

L. M. Cathles, H. J. Erendi, and T. Barrie, “How long can a hydrothermal system be sustained by a single intrusive event?” Economic Geology, vol. 92, no. 7-8, pp. 766–771, 1997.

A. Weisheit, P. D. Bons, M. Danisik, and M. A. Elburg, “Crustal-scale folding: Palaeozoic deformation of the Mt. Painter Inlier, South Australia,” in Deformation Structures and Processes within the Continental Crust, vol. 394, pp. 53–77, Geological Society, London, Special Publications, 2013.

S. Rojstaczer, S. Wolf, and R. Michel, “Permeability enhancement in the shallow crust as a cause of earthquake-induced hydrological changes,” Nature, vol. 373, no. 6511, pp. 237–239, 1995.

S. A. Rojstaczer, S. E. Ingebritsen, and D. O. Hayba, “Permeability of continental crust influenced by internal and external forcing,” Geofluids, vol. 8, no. 2, pp. 128–139, 2008.

R. H. Sibson, F. Robert, and K. H. Poulsen, “High-angle reverse faults, fluid-pressure cycling, and mesothermal gold-quartz deposits,” Geology, vol. 16, no. 6, pp. 551–555, 1988.

G. Preisig, E. Eberhardt, V. Gischig et al., “Development of connected permeability in massive crystalline rocks through hydraulic fracture propagation and shearing accompanying fluid injection,” Geofluids, vol. 15, no. 1-2, pp. 321–337, 2015.

P. Weis, “The dynamic interplay between saline fluid flow and rock permeability in magmatic-hydrothermal systems,” Geofluids, vol. 15, no. 1-2, pp. 350–371, 2015.

S. D. Burley, J. T. Mullis, and A. Matter, “Timing diagenesis in the Tartan Reservoir (UK North Sea): constraints from combined cathodoluminescence microscopy and fluid inclusion studies,” Marine and Petroleum Geology, vol. 6, no. 2, pp. 98–120, 1989.

J. J. Roering, J. W. Kirchner, L. S. Sklar, and W. E. Dietrich, “Hillslope evolution by nonlinear creep and landsliding: an experimental study,” Geology, vol. 29, no. 2, pp. 143–146, 2001.

P. Bak, C. Tang, and K. Wiesenfeld, “Self-organized criticality,” Physical Review A, vol. 38, no. 1, pp. 364–374, 1988.

R. Sánchez, D. E. Newman, and B. A. Carreras, “Mixed SOC diffusive dynamics as a paradigm for transport in fusion devices,” Nuclear Fusion, vol. 41, no. 3, pp. 247–256, 2001.

E. Milanese, O. Yılmaz, J. F. Molinari, and B. Schrefler, “Avalanches in dry and saturated disordered media at fracture,” Physical Review E, vol. 93, no. 4, article 043002, 2016.

T. D. Cao, F. Hussain, and B. A. Schrefler, “Porous media fracturing dynamics: stepwise crack advancement and fluid pressure oscillations,” Journal of the Mechanics and Physics of Solids, vol. 111, pp. 113–133, 2018.

D. L. Turcotte, “Fractals and fragmentation,” Journal of Geophysical Research: Solid Earth, vol. 91, no. B2, pp. 1921–1926, 1986.

F. P. Agterberg, Q. Cheng, A. Brown, and D. Good, “Multifractal modeling of fractures in the Lac du Bonnet batholith, Manitoba,” Computers & Geosciences, vol. 22, no. 5, pp. 497–507, 1996.

B. Berkowitz and A. Hadad, “Fractal and multifractal measures of natural and synthetic fracture networks,” Journal of Geophysical Research: Solid Earth, vol. 102, no. B6, pp. 12205–12218, 1997.

E. Bonnet, O. Bour, N. E. Odling et al., “Scaling of fracture systems in geological media,” Reviews of Geophysics, vol. 39, no. 3, pp. 347–383, 2001.

B. Gutenberg and C. F. Richter, Seismicity of the earth, p. 310, 1954.

S. H. Advani, T. S. Lee, R. H. Dean, C. K. Pak, and J. M. Avasthi, “Consequences of fluid lag in three-dimensional hydraulic fractures,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 21, no. 4, pp. 229–240, 1997.

P. C. Papanastasiou, “A coupled elastoplastic hydraulic fracturing model,” International Journal of Rock Mechanics and Mining Sciences, vol. 34, no. 3-4, pp. 240.e1–240.e15, 1997.

J. Vychytil and H. Horii, “Micromechanics-based continuum model for hydraulic fracturing of jointed rock masses during HDR stimulation,” Mechanics of Materials, vol. 28, no. 1-4, pp. 123–135, 1998.

F. Kraaijeveld, J. M. Huyghe, J. J. C. Remmers, and R. De Borst, “Two-dimensional mode I crack propagation in saturated ionized porous media using partition of unity finite elements,” Journal of Applied Mechanics, vol. 80, no. 2, 2013.

D. Koehn, J. Arnold, and C. W. Passchier, “Fracture and vein patterns as indicators of deformation history: a numerical study,” Geological Society, London, Special Publications, vol. 243, no. 1, pp. 11–24, 2005.

I. Ghani, D. Koehn, R. Toussaint, and C. W. Passchier, “Dynamic development of hydrofracture,” Pure and Applied Geophysics, vol. 170, no. 11, pp. 1685–1703, 2013.

P. Grassl, C. Fahy, D. Gallipoli, and S. J. Wheeler, “On a 2D hydro-mechanical lattice approach for modelling hydraulic fracture,” Journal of the Mechanics and Physics of Solids, vol. 75, pp. 104–118, 2015.

T. Sachau, P. D. Bons, and E. Gomez-Rivas, “Transport efficiency and dynamics of hydraulic fracture networks,” Frontiers in Physics, vol. 3, p. 63, 2015.

G. I. Barenblatt, V. M. Entov, and V. M. Ryzhik, Theory of Fluid Flows through Natural Rocks, 1989.

R. A. Fine and F. J. Millero, “Compressibility of water as a function of temperature and pressure,” The Journal of Chemical Physics, vol. 59, no. 10, pp. 5529–5536, 1973.

A. Liebscher, “Aqueous fluids at elevated pressure and temperature,” Geofluids, vol. 10, no. 1-2, pp. 3–19, 2010.

C. A. Schneider, W. S. Rasband, and K. W. Eliceiri, “NIH Image to ImageJ: 25 years of image analysis,” Nature Methods, vol. 9, no. 7, pp. 671–675, 2012.

H. E. Hurst, “Long-term storage capacity of reservoirs,” Transactions. American Society of Civil Engineers, vol. 116, pp. 770–799, 1951.

B. B. Mandelbrot and J. R. Wallis, “Noah, Joseph, and operational hydrology,” Water Resources Research, vol. 4, no. 5, pp. 909–918, 1968.

D. L. Turcotte, “Self-organized criticality,” Reports on Progress in Physics, vol. 62, no. 10, pp. 1377–1429, 1999.

D. M. Fisher and S. L. Brantley, “Models of quartz overgrowth and vein formation: deformation and episodic fluid flow in an ancient subduction zone,” Journal of Geophysical Research, vol. 97, no. B13, p. 20043, 1992.

C. M. Breeding and J. J. Ague, “Slab-derived fluids and quartz-vein formation in an accretionary prism, Otago Schist, New Zealand,” Geology, vol. 30, no. 6, pp. 499–502, 2002.

Å. Fagereng and C. Harris, “Interplay between fluid flow and fault–fracture mesh generation within underthrust sediments: geochemical evidence from the Chrystalls Beach Complex, New Zealand,” Tectonophysics, vol. 612, pp. 147–157, 2014.

D. Sornette, “Dragon-kings, black swans and the prediction of crises,” 2009, http://arxiv.org/abs/0907.4290.

M. K. Sachs, M. R. Yoder, D. L. Turcotte, J. B. Rundle, and B. D. Malamud, “Black swans, power laws, and dragon-kings: earthquakes, volcanic eruptions, landslides, wildfires, floods, and SOC models,” The European Physical Journal Special Topics, vol. 205, no. 1, pp. 167–182, 2012.

P. D. Bons, J. K. Becker, M. A. Elburg, and K. Urtson, “Granite formation: stepwise accumulation or connected networks?” Earth and Environmental Science Transactions of the Royal Society of Edinburgh, vol. 100, pp. 105–115, 2010.

K. Pfaff, R. L. Romer, and G. Markl, “U-Pb ages of ferberite, chalcedony, agate, 'U-mica' and pitchblende: constraints on the mineralization history of the Schwarzwald ore district,” European Journal of Mineralogy, vol. 21, no. 4, pp. 817–836, 2009.