Structural Controls on Basin and CrustalScale Fluid Flow and Resulting Mineral Reactions
View this Special IssueResearch Article  Open Access
Tamara de Riese, Paul D. Bons, Enrique GomezRivas, Till Sachau, "Interaction between CrustalScale Darcy and Hydrofracture Fluid Transport: A Numerical Study", Geofluids, vol. 2020, Article ID 8891801, 14 pages, 2020. https://doi.org/10.1155/2020/8891801
Interaction between CrustalScale Darcy and Hydrofracture Fluid Transport: A Numerical Study
Abstract
Crustalscale 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 powerlaw size distribution. A Hurst factor of ~0.9 indicates that the system selforganizes. The abundant smallscale hydrofractures organize the formation of largescale 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 crustalscale 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 “crackseal 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 crackseal process thus implies strongly varying fluid fluxes, as well as stress states that must repeatedly reach the failure criterion to induce fracturing. Tensional crackseal 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 cracksealing 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 largescale connectivity of fracture systems, while dissolution may increase them [19]. Fracturing can also occur locally due to reactiondriven 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 wellstudied 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^{2} 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 530 km^{3}.
(a)
(b)
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^{2}) 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 onedimensional form in the vertical direction (positive downwards).
Permeability varies by >16 orders of magnitude in geological materials (from in crystalline rocks to in wellsorted 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 topographydriven flow, variations in fluid density due to temperature or salinity differences may result in convective flow (e.g., [1]). In topographydriven 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 orebearing 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^{2} to 10^{1} 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 socalled 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, longlived (10^{3}10^{6} years) hydrothermal systems exist [32, 63, 64], which indicate that processes such as hydrofracturing seem to be an important mechanism [65]. Largescale 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 permeabilitydepth 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 selforganization, which has been proposed for, e.g., the faultvalve behaviour in fractures [51, 52, 67], hydraulic fracturing [68], and in magmatichydrothermal 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]. Crustalscale permeability is therefore a dynamically selfadjusting 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 powerlaw 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 frequencysize distributions of rock fragments [32, 76], in the size distribution of fractures [77–79], the magnitudefrequency 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 particlelattice 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 selforganized 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 selforganized critical type of transport when transport in hydrofractures is activated. Miller and Nur [58] used a crustalscale 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 longrange 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 endmember 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^{1}) of the fluid. We can express this with:
is the reference volume (one m^{3} 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^{2}/s).
Following Fick’s second law, pressure gradients can dissipate in a diffusional fashion with an effective diffusivity (in m^{2}/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^{5} to 10^{1} m^{2}/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 lowpermeability rocks. Note that all used symbols and parameters are given in Table 1.

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, 2dimensional 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^{4} Pa/m with depth (assuming a fluid density of 1000 kg/m^{3} and a gravitational acceleration, , set at 10 m/s^{2}). 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.
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^{3}, 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 crackseal 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 healingrate 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 Darcianflow dominated (high ) behaviour.
2.3. Scaling
Scaling of the model values (subscript ) to realworld 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^{4} Pa/m. From this, we obtain for the scaling between and the realworld :
We assume a fixed fluid flux , which is typical for metamorphic fluid flux in the crust at a depth of 1015 km [39]. This influx enters one element along a window of m^{2} (assuming the model is 1 m thick in the third dimension) and is added to the fluid residing in the m^{3} 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 realworld diffusivity () and permeability ():
For the given and m, we obtain . Assuming , the chosen range of modelled values from 10^{7} 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.

2.4. Analysis
Frequency distributions of hydrofracture sizes have been evaluated. Simulation visualisation and calculation of the timeaveraged vertical pressure profile have been done with ImageJ [92]. The timeaveraged vertical pressure profile can be calculated via ImageStacks–Zproject and AnalyzePlot Profile. The effective diffusivity is the slope of the timeaveraged vertical pressure profile. Besides this, the rescaled range () statistics was calculated for the pressure fluctuations to quantify the degree of selfsimilarity. 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 selfsimilarity 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^{4}. 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.
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 fastflowing 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 lowpressure 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 faultvalve 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.
Whereas Figure 4 shows the mean pressure of the whole model area as a function of time, Figure 5 shows the timeaveraged 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 pressuredepth 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 shows that mean fluidpressure 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 timeaveraged 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 lowpressure 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 rescaledrange 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 selforganized 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)
3.2. Statistical Properties of Hydrofracture Size Distributions
Similar models and experiments [56, 72, 73] showed that flow avalanche magnitudes tend to develop powerlaw 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)
The resulting bestfit values for , the frequency of the smallest possible hydrofracture of one single element, and the powerlaw exponent () are listed in Table 3. The deviation of the powerlaw 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.

Hydrofractures that reach the surface are invariably large and their sizes do not follow the aforementioned powerlaw 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 lognormal 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 powerlaw distribution and spread fluid within the crust. These powerlaw distributions imply that rare but large events transport most of the fluid. Such distributions are often invoked to indicate a selforganized criticality (SOC) in the underlying dynamics of the system [72, 95]. The SOCtype 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 crackseal 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 greenschistfacies rocks. The lower part of the models with abundant hydrofractures may represent these veinrich 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 powerlaw distributions of hydrofractures not reaching the top of the model. These large events can be interpreted as “dragonkings” [99]. Such “dragon kings” are outliers which coexist with powerlaw distributions of event sizes but exceed the powerlaw extrapolations and are often associated with the occurrence of a bifurcation, catastrophes, or tipping points [99]. They also have been found in, for example, frequencyrupture area statistics of earthquakes [100]. Although not following a powerlaw distribution themselves, they are intimately related to the smaller hydrofractures that (self) organize the fluidpressure distribution and the initiation of large fluidescape 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^{2} Hidden Valley hydrothermal megabreccia (Figure 1(b)) formed over a period of ~150200 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 crustalscale 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 selforganizes. Abundant hydrofractures form at the base of the model and their sizefrequency distributions are power laws. The Hurst parameter of ~0.9, calculated from the mean pressure variations over time, supports the development of a selforganized critical state. The hydrofracture size distributions show “dragonking”like large hydrofractures that deviate from the powerlaw 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 crackseal 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 RYC2018026335I and research project PGC2018093903BC22).
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. View at: Publisher Site  Google Scholar
 S. A. Shapiro and C. Dinske, “Fluidinduced seismicity: pressure diffusion and hydraulic fracturing,” Geophysical Prospecting, vol. 57, no. 2, pp. 301–310, 2009. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 D. E. Huenges, Ed., Geothermal Energy Systems: Exploration, Development, and Utilization, WileyVCH 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. View at: Publisher Site  Google Scholar
 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 basementcover interfaces during basinrelated mineralization,” Economic Geology, vol. 101, no. 1, pp. 1–31, 2006. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Google Scholar
 C. Zhao, B. E. Hobbs, and A. Ord, “Convective and advective heat transfer in geological systems,” Geophysics & Geodesy, vol. 229, 2008. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar
 P. D. Bons, M. A. Elburg, and E. GomezRivas, “A review of the formation of tectonic veins and their microstructures,” Journal of Structural Geology, vol. 43, pp. 33–62, 2012. View at: Publisher Site  Google Scholar
 J. G. Ramsay, “The crack–seal mechanism of rock deformation,” Nature, vol. 284, no. 5752, pp. 135–139, 1980. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 A. Vass, D. Koehn, R. Toussaint, I. Ghani, and S. Piazolo, “The importance of fracturehealing on the deformation of fluidfilled layered systems,” Journal of Structural Geology, vol. 67, pp. 94–106, 2014. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 G. R. Watt, N. H. S. Oliver, and B. J. Griffin, “Evidence for reaction induced microfracturing in granulite facies migmatites,” Geology, vol. 28, no. 4, pp. 327–330, 2000. View at: Publisher Site  Google Scholar
 D. S. O’Hanley, “Serpentinites: records of tectonic and petrological history,” Oxford monographs on geology and geophysics, vol. 34, 1996. View at: Google Scholar
 K. Iyer, B. Jamtveit, J. Mathiesen, A. MaltheSørenssen, and J. Feder, “Reactionassisted hierarchical fracturing during serpentinization,” Earth and Planetary Science Letters, vol. 267, no. 34, pp. 503–516, 2008. View at: Publisher Site  Google Scholar
 B. Jamtveit, A. MaltheSørenssen, and O. Kostenko, “Reaction enhanced permeability during retrogressive metamorphism,” Earth and Planetary Science Letters, vol. 267, no. 34, pp. 620–627, 2008. View at: Publisher Site  Google Scholar
 R. C. Fletcher, H. L. Buss, and S. L. Brantley, “A spheroidal weathering model coupling porewater chemistry to soil thicknesses during steadystate denudation,” Earth and Planetary Science Letters, vol. 244, no. 12, pp. 444–457, 2006. View at: Publisher Site  Google Scholar
 B. Jamtveit, C. V. Putnis, and A. MaltheSørenssen, “Reaction induced fracturing during replacement processes,” Contributions to Mineralogy and Petrology, vol. 157, no. 1, pp. 127–133, 2009. View at: Publisher Site  Google Scholar
 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. View at: Google Scholar
 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. View at: Google Scholar
 N. H. Woodcock and K. Mort, “Classification of fault breccias and related fault rocks,” Geological Magazine, vol. 145, no. 3, pp. 435–440, 2008. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 M. Jébrak, “Hydrothermal breccias in veintype ore deposits: a review of mechanisms, morphology and size distribution,” Ore Geology Reviews, vol. 12, no. 3, pp. 111–134, 1997. View at: Publisher Site  Google Scholar
 A. Weisheit, P. D. Bons, and M. A. Elburg, “Longlived crustalscale fluid flow: the hydrothermal megabreccia of Hidden Valley, Mt. Painter Inlier, South Australia,” International Journal of Earth Sciences, vol. 102, no. 5, pp. 1219–1236, 2013. View at: Publisher Site  Google Scholar
 S. Staude, P. D. Bons, and G. Markl, “Hydrothermal vein formation by extensiondriven dewatering of the middle crust: an example from SW Germany,” Earth and Planetary Science Letters, vol. 286, no. 34, pp. 387–395, 2009. View at: Publisher Site  Google Scholar
 B. F. Walter, M. Burisch, and G. Markl, “Longterm 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. View at: Publisher Site  Google Scholar
 P. D. Bons, “The formation of large quartz veins by rapid ascent of fluids in mobile hydrofractures,” Tectonophysics, vol. 336, no. 14, pp. 1–17, 2001. View at: Publisher Site  Google Scholar
 P. D. Bons, T. Fusswinkel, E. GomezRivas, G. Markl, T. Wagner, and B. Walter, “Fluid mixing from below in unconformityrelated hydrothermal ore deposits,” Geology, vol. 42, no. 12, pp. 1035–1038, 2014. View at: Publisher Site  Google Scholar
 S. E. Ingebritsen and T. Gleeson, “Crustal permeability: introduction to the special issue,” Geofluids, vol. 15, no. 12, pp. 1–10, 2015. View at: Publisher Site  Google Scholar
 J. Bear, Dynamics of Fluids in Porous Media, Courier Corporation, New York, 2013.
 S. E. Ingebritsen and C. E. Manning, “Geological implications of a permeabilitydepth curve for the continental crust,” Geology, vol. 27, no. 12, pp. 1107–1110, 1999. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 S. E. Ingebritsen and C. E. Manning, “Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism,” Geofluids, vol. 10, no. 12, pp. 193–205, 2010. View at: Google Scholar
 M. K. Hubbert, “The theory of groundwater motion,” The Journal of Geology, vol. 48, 8, Part 1, pp. 785–944, 1940. View at: Publisher Site  Google Scholar
 S. F. Cox, “The application of failure mode diagrams for exploring the roles of fluid pressure and stress states in controlling styles of fracturecontrolled permeability enhancement in faults and shear zones,” Geofluids, vol. 10, no. 12, pp. 217–233, 2010. View at: Google Scholar
 N. H. S. Oliver, M. J. Rubenach, B. Fu et al., “Graniterelated overpressure and volatile release in the mid crust: fluidized breccias from the Cloncurry District, Australia,” Geofluids, vol. 6, no. 4, pp. 346–358, 2006. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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 ZnPb deposit, Yunnan, China,” Geoscience Frontiers, vol. 3, no. 1, pp. 73–84, 2012. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 A. Okamoto and N. Tsuchiya, “Velocity of vertical fluid ascent within veinforming fractures,” Geology, vol. 37, no. 6, pp. 563–566, 2009. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 R. H. Sibson, “Crustal stress, faulting and fluid flow,” Geological Society, London, Special Publications, vol. 78, no. 1, pp. 69–84, 1994. View at: Publisher Site  Google Scholar
 R. H. Sibson, “Fluid involvement in normal faulting,” Journal of Geodynamics, vol. 29, no. 35, pp. 469–499, 2000. View at: Publisher Site  Google Scholar
 R. H. Sibson, “Tectonic controls on maximum sustainable overpressure: fluid redistribution from stress transitions,” Journal of Geochemical Exploration, vol. 69, pp. 471–475, 2000. View at: Google Scholar
 A. Nur and J. Walder, “Hydraulic pulses in the Earth’s crust,” in International Geophysics, vol. 51, pp. 461–473, Academic Press, 1992. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 P. D. Bons and B. P. van Milligen, “New experiment to model selforganized critical transport and accumulation of melt and hydrocarbons from their source rocks,” Geology, vol. 29, no. 10, pp. 919–922, 2001. View at: Publisher Site  Google Scholar
 E. GomezRivas, 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. View at: Publisher Site  Google Scholar
 S. A. Miller and A. Nur, “Permeability as a toggle switch in fluidcontrolled crustal processes,” Earth and Planetary Science Letters, vol. 183, no. 12, pp. 133–146, 2000. View at: Publisher Site  Google Scholar
 J. Weertman, “Theory of waterfilled crevasses in glaciers applied to vertical magma transport beneath oceanic ridges,” Journal of Geophysical Research, vol. 76, no. 5, pp. 1171–1183, 1971. View at: Publisher Site  Google Scholar
 T. Dahm, “On the shape and velocity of fluidfilled fractures in the Earth,” Geophysical Journal International, vol. 142, no. 1, pp. 181–192, 2000. View at: Publisher Site  Google Scholar
 J. D. MartínMartín, E. GomezRivas, D. GómezGras 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. View at: Publisher Site  Google Scholar
 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 61–ECV 617, 2002. View at: Publisher Site  Google Scholar
 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. 78, pp. 766–771, 1997. View at: Publisher Site  Google Scholar
 A. Weisheit, P. D. Bons, M. Danisik, and M. A. Elburg, “Crustalscale 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. View at: Google Scholar
 S. Rojstaczer, S. Wolf, and R. Michel, “Permeability enhancement in the shallow crust as a cause of earthquakeinduced hydrological changes,” Nature, vol. 373, no. 6511, pp. 237–239, 1995. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 R. H. Sibson, F. Robert, and K. H. Poulsen, “Highangle reverse faults, fluidpressure cycling, and mesothermal goldquartz deposits,” Geology, vol. 16, no. 6, pp. 551–555, 1988. View at: Publisher Site  Google Scholar
 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. 12, pp. 321–337, 2015. View at: Publisher Site  Google Scholar
 P. Weis, “The dynamic interplay between saline fluid flow and rock permeability in magmatichydrothermal systems,” Geofluids, vol. 15, no. 12, pp. 350–371, 2015. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 P. Bak, C. Tang, and K. Wiesenfeld, “Selforganized criticality,” Physical Review A, vol. 38, no. 1, pp. 364–374, 1988. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 D. L. Turcotte, “Fractals and fragmentation,” Journal of Geophysical Research: Solid Earth, vol. 91, no. B2, pp. 1921–1926, 1986. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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 threedimensional hydraulic fractures,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 21, no. 4, pp. 229–240, 1997. View at: Publisher Site  Google Scholar
 P. C. Papanastasiou, “A coupled elastoplastic hydraulic fracturing model,” International Journal of Rock Mechanics and Mining Sciences, vol. 34, no. 34, pp. 240.e1–240.e15, 1997. View at: Google Scholar
 J. Vychytil and H. Horii, “Micromechanicsbased continuum model for hydraulic fracturing of jointed rock masses during HDR stimulation,” Mechanics of Materials, vol. 28, no. 14, pp. 123–135, 1998. View at: Publisher Site  Google Scholar
 F. Kraaijeveld, J. M. Huyghe, J. J. C. Remmers, and R. De Borst, “Twodimensional mode I crack propagation in saturated ionized porous media using partition of unity finite elements,” Journal of Applied Mechanics, vol. 80, no. 2, 2013. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 P. Grassl, C. Fahy, D. Gallipoli, and S. J. Wheeler, “On a 2D hydromechanical lattice approach for modelling hydraulic fracture,” Journal of the Mechanics and Physics of Solids, vol. 75, pp. 104–118, 2015. View at: Publisher Site  Google Scholar
 T. Sachau, P. D. Bons, and E. GomezRivas, “Transport efficiency and dynamics of hydraulic fracture networks,” Frontiers in Physics, vol. 3, p. 63, 2015. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar
 A. Liebscher, “Aqueous fluids at elevated pressure and temperature,” Geofluids, vol. 10, no. 12, pp. 3–19, 2010. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 H. E. Hurst, “Longterm storage capacity of reservoirs,” Transactions. American Society of Civil Engineers, vol. 116, pp. 770–799, 1951. View at: Google Scholar
 B. B. Mandelbrot and J. R. Wallis, “Noah, Joseph, and operational hydrology,” Water Resources Research, vol. 4, no. 5, pp. 909–918, 1968. View at: Publisher Site  Google Scholar
 D. L. Turcotte, “Selforganized criticality,” Reports on Progress in Physics, vol. 62, no. 10, pp. 1377–1429, 1999. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 C. M. Breeding and J. J. Ague, “Slabderived fluids and quartzvein formation in an accretionary prism, Otago Schist, New Zealand,” Geology, vol. 30, no. 6, pp. 499–502, 2002. View at: Publisher Site  Google Scholar
 Å. 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. View at: Google Scholar
 D. Sornette, “Dragonkings, black swans and the prediction of crises,” 2009, http://arxiv.org/abs/0907.4290. View at: Google Scholar
 M. K. Sachs, M. R. Yoder, D. L. Turcotte, J. B. Rundle, and B. D. Malamud, “Black swans, power laws, and dragonkings: earthquakes, volcanic eruptions, landslides, wildfires, floods, and SOC models,” The European Physical Journal Special Topics, vol. 205, no. 1, pp. 167–182, 2012. View at: Publisher Site  Google Scholar
 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. View at: Google Scholar
 K. Pfaff, R. L. Romer, and G. Markl, “UPb ages of ferberite, chalcedony, agate, 'Umica' and pitchblende: constraints on the mineralization history of the Schwarzwald ore district,” European Journal of Mineralogy, vol. 21, no. 4, pp. 817–836, 2009. View at: Publisher Site  Google Scholar
Copyright
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.