#### 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^{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 5-30 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 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^{-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 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^{3}-10^{6} 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^{-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 low-permeability 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, 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^{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 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^{4} 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^{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 real-world 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 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^{-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 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.

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 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)**

##### 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)**

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.

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^{2} 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)*