Mechanical origin of b-value changes during stimulation of deep geothermal reservoirs

Abstract

Background

The distribution of induced-earthquake magnitudes in deep geothermal reservoirs is a classical tool for monitoring reservoirs. It typically shows some important fluctuations through time and space. Despite being a very crude information (i.e., a scalar quantity) of very complex mechanical stress evolution, understanding these variations could still give us insights into the mechanics of the reservoir. Here, we analyze the output of a simple quasi-static physical model of a single fault and propose a new way to describe bursts that could be compared to seismic events.

Methods

Our model is an elastic fiber bundle model describing multiple ruptures along a single major fault of the reservoir. It consists of a set of fibers with various strengths, arranged in a planar L×L two-dimensional array, linking together two elastic half-spaces. It mimics a fault zone with various asperities. During load, the fibers break quasi-statically according to a stress threshold distribution. Contrary to classical fiber bundle model, here, when a fiber breaks, it redistributes the load on the surviving fibers through long range elastic interactions. Interestingly, the elasticity of the half-spaces which changes the range of the stress distribution, characterizes two distinct regimes. In a stiff regime, we find an effective ductile deformation regime at large scale since the damage distribution is very difuse. Conversely in softer systems, the distance between two consecutive breaking fibers gets smaller and failure of the interface is localized, exhibiting an effective brittle regime.

Results

We analyze two types of burst distributions: a classical one built from the statistics of the broken bonds during each failure step and a new one defines from a waiting time matrix of the fracture front propagation. The first one reproduces several known results for this type of model. The new one evidences the existence of effective creeping advances of the front with statistics that follow a Gutenberg-Richter distribution, in particular, in the ductile regime (stiff systems).

Conclusions

We proposed a new definition of bursts in a fault model, based on the local fracture front velocity. We find that b _v -values for the distribution of the velocity clusters are very consistent with Gutenberg-Richter distribution of induced seismicity. We link the b _v -value fluctuations in the reservoirs to the influence of the velocity threshold level that could be related to recording limitation. b _b -values obtained from the broken bond statistics are hardly comparable to seismic events because of the lack of space contiguousness of broken fibers during bursts. In the light of our results, we discuss the implications of b-value changes in geothermal reservoir in terms of fault asperities and normal stress evolution.

Background

Earthquakes recorded worldwide from various natural fault environments exhibit a specific magnitude distribution. This distribution is well represented by the Gutenberg-Richter empirical law which states that:

$$ \log_{10}N(M) = a -bM, $$

((1))

where N(M) is the number of events with a magnitude greater or equal to the magnitude M, and a and b are the parameters that are inverted from the earthquake magnitude dataset [Gutenberg and Richter (1944)]. The parameter a scales with the number of earthquakes in the investigated catalog while the parameter b is an exponent which quantifies the overall proportion between larger and smaller magnitude events in the earthquake catalog. The validity of Equation 1 has been observed over a wide range of scales. The scale-free nature of the Gutenberg-Richter relation suggests that laws governing seismic activity at large scales (i.e., fault scales) are the same as at small scales (i.e., laboratory scale) for rock samples [Scholz CH (2002)]. Although the empirical relation of Equation 1 is robust and tested against a large number of dataset of natural earthquakes and acoustic emissions in laboratory experiments, the exponent b exhibits fluctuations. While generally close to b=1.0, e.g., Schorlemmer et al. (2005) and Kwiatek et al. (2010), several variations of the b value are also noticed, e.g., Schorlemmer et al. (2005), Enescu and Ito (2001), and Amitrano (2003). These fluctuations attest for a possible physical control on the value of this exponent and its observed fluctuations.

Apart from natural systems and laboratory experiments, the Gutenberg-Richter distribution of earthquakes is also observed for induced seismicity. One example of the later is the case of the stimulation of deep geothermal reservoir. Indeed, for the development of deep geothermal reservoirs (Enhanced Geothermal Systems (EGS)) like at Soultz-sous Forêts, France [Dorbath et al. (2009)] or Basel, Switzerland [Bachmann et al. (2011)], hydraulic injections have been used to stimulate the existing fracture network both in the sediments and in the granite massif in order to enhance the permeability of the reservoir. For instance, during the hydraulic injection of the Soultz-sous-Forêts GPK2 well in 2000, a significant seismic activity was generated, mapping the development of the fault network [Cuenot et al. (2008)]. Up to 30,000 earthquakes were recorded with a magnitude ranging from −0.9 to 2.6. Earthquake magnitude distribution was found to follow the typical Gutenberg-Richter distribution. During the 2000 injection, the b-value was measured to be 1.29 for the entire period [Cuenot et al. (2008)]. However, when searching for temporal variations, the b-value was observed to evolve from 1.0 to 1.7. High injection rates led to higher b-value. Similar variations of the b-value during injection were observed in Basel during the 2006 experiment. During the injection, b-value was measured to be 1.58 but lowered to 1.15 after shut-in [Bachmann et al. (2011)]. These temporal and spatial fluctuations of the b-value suggest to use this information as a tracer of the evolution of the the reservoir (in the case of induced seismicity) or the fault system (in natural environment).

Numerous models have been proposed to explain the variability of the b-value. One of the first model is a theoretical study of rock fracture developed by Scholz (1968). This model was based on a heterogeneous distribution of material strength in a sample under load. It was found that the b-value is a function of the applied stress, and that the frequency of fracture size was scale free. Similar results were obtained by Wyss (1973) suggesting that low b-values correlate with high stresses. Although these models bring a possible explanation of the observed b-value fluctuations, some care should be taken when interpreting stress from the b-value alone since it is a scalar quantity that cannot characterize completely the stress tensor. Many other models concerning the connection between earthquakes and micro-fractures have been developed through the years in the language of statistical physics, e.g., Holliday et al. (2008). In particular, approaches that connect fiber-bundle models to acoustic emissions during rock fracturing have been done [Turcotte et al. (2003); Pradhan et al. (2010)]. For an overlook of this particular field, we refer to the review by Rundle et al. (2003).

Several attempts have also been made in order to model the b-value and its fluctuations specifically for induced seismicity in geothermal system, e.g., Hakimhashemi et al. (2014). Baisch et al. (2010) considered a single fault to model the induced seismicity at Soultz-sous-Forêts. They indeed noticed that hypocenters are aligned along a preexisting planar fault structure which concentrates most of the recorded seismicity. Their quasi-static model includes sliders connected by springs that incorporate elastic stress redistribution. The effect of fluid is simply modeled as a change of pore pressure bringing each path of the model closer to failure in the Mohr-Coulomb sense. Based on this simple model, they reproduced most of the seismicity characteristics: migration of the seismicity, Kaiser effect and magnitude distribution. Variation of the b-value was also explained as the result of a dominant portion of smaller cracks close to the injection well, compared to the rock mass away from the well (far-field), e.g., Zang et al. (2014). Some more complicated models designed to incorporate more realistic fluid transfer were also proposed. In particular, Yoon et al. (2014) proposed a model of aggregated circular particles where both mode I and mode II ruptures are taken into account to compute fluid-induced seismicity in naturally fractured granitic EGS reservoirs. Their modeling was capable of reproducing most of the observed features of the seismicity including the frequency magnitude distribution. Similar results were also obtained by Gischig and Wiemer (2013) while using a completely different model. They incorporate a stochastic seed model generating synthetic seismic catalogs reproducing the Basel seismicity data. Other continuous models that include elastic long range interactions together with poro-elastic behavior have also been developed and show self-induced seismicity when shear rupture is dominant over the fluid migration process [Aochi et al. (2014)].

It then appears that many different models from very different physical mechanisms are able to reproduce the frequency magnitude distribution of the seismicity in EGS reservoirs. It therefore questions about how to interpret the change of b-values in the light of so different modeling. We decided here to study another aspect: the evolution of damage and its link with seismic events considering only a quasi-static model, where damages are represented on a single fault plane and where the fluid source is not explicitly modeled. Although our approach is oversimplified compared to the rich complexity of induced seismicity of geothermal system, it allows us to test the direct influence of several parameters such as the stiffness of the system. Our model is inspired by the model developed by Amitrano et al. (1999) that has shown very promising results in connection to observed seismicity. In this model, load is applied to an initially intact sample. When the stress exceeds the local strength, the sample is damaged. The effect of this damage is represented by a change in local Young’s modulus. Thus, increasing the stress in the vicinity of the damage zone increases the probability of damage to the neighboring area. This model has shown to have a brittle to ductile transition as a function of the internal friction angle of the Mohr-Coulomb failure criterion [Amitrano (2003)]. In the ductile regime (i.e., low internal frictional angle), the damage in the fracture process is highly distributed, while the failure in the brittle regime is highly localized for a high internal frictional angle. The transition also affects the b-value of the model. Completely ductile failure gives a b-value of 0.63, while for brittle fractures, b=1.2. A recently developed bottom-up model of fracture by Stormo et al. (2013) has shown to have similar diffuse to localized damage distribution transition. This model can address the competition between elastic stress redistribution and damage development related to geological heterogeneity. Here, we refine and expand the model proposed by Stormo et al. (2013) with the attempt to describe fracture development along a pre-existing weak interface. We observe in our model a similar effective ductile to brittle transition at large scale as observed by previous authors. We focus on the analysis of burst distributions. We notably explore the impact of the definition of the event in the results of the frequency magnitude distribution. We indeed introduced a new definition of a burst based on the local velocity of the crack front propagation. We show that this definition provides an interesting perspective for comparison with large scale fluctuations of the b-values.

Methods

Our model is a variation of the fiber bundle model, i.e., a model of a single fracture described as an interfacial failure along a discrete set of connecting fibers. We assume that the fault network in the geothermal reservoir is dominated by a discrete set of pre-existing faults as it was shown by (Sausse et al. 2010) at Soultz-sous-Forêts. We even push further the hypothesis and reduce the fault network for seismicity modeling to a single main fault zone as previously proposed by Baisch et al. (2010). The response of the main fault is seen as a collection of weak fibers attached to the fault boundaries displaced apart by a distance D because of hydraulic pressure increase. Our variation of the fiber bundle model is to deal with the elastic fiber bundle model, first proposed by Batrouni et al. (2002b). It is an attempt to include non-local elastic mechanics to the fiber bundle models by replacing the rigid boundaries of the global load sharing configuration by deformable elastic half-spaces. Actually, this configuration of two facing elastic half-spaces is equivalent to a single elastic half-space with equivalent elastic moduli facing a rigid block as discussed in Batrouni et al. (2002b). Variability of the fracture aperture could have been introduced by considering a specific morphology of the rigid part [Hansen et al. (2000)], but this is not the case in the present study. However, we consider the presence of a static fluid in the fracture by dealing with normal effective stresses on the fibers but neglect other specific fluid effects like viscous pressure or dynamic effects [Jaeger et al. (2009)].

The elastic fiber-bundle model

In our model that follows Batrouni et al. (2002b), the fibers are arranged in a planar two-dimensional square array of L×L=N fibers distanced apart from their nearest neighbors by a distance a. Our square array is considered periodic in both in-plane directions. The fibers are attached to two infinite half-spaces, representing the clamps holding the fibers. One of these half-spaces is supposed to be elastic with a Young’s modulus E and Poisson’s ratio ν. The other one is assumed to be perfectly rigid. This configuration is actually equivalent to two facing elastic half-spaces as discussed in [Batrouni et al. (2002b)]. Figure 1 gives an illustration of the model setup.

Figure 1

Setup of the elastic fiber bundle model. On the left, a sketch of the fiber bundle. The upper half-space deflects in response to the forces applied by the fibers (red), while the infinitely rigid lower half-space is displaced downward of a distance D. To the right, a schematic top view of the entire fiber bundle where each red circle represents a fiber at the center of a mesh of size a×a.

Each fiber is modeled as a Hookean spring, and thus exercises an opposite force, f, on a square a ² of each facing half-spaces:

$$ f_{j}=-k(D-u_{j}). $$

((2))

Here, k is the spring constant of the fiber, D is the explicit imposed displacement to the rigid half-space with respect to the elastic half-space as seen in Figure 1, and u is the local displacement of the surface of the elastic half-space from equilibrium due to the force exercised by all the other fibers in the array. j is the index of the fiber. The spatially dependent displacement of the elastic medium owing to this force was calculated by 1929. He obtains the displacement of the fiber j as:

$$ u_{j}=\sum_{i\ne j}^{N} G_{i} f_{i}. $$

((3))

The Green’s function is:

$$ G_{i}=\frac{1-\nu^{2}}{\pi E a^{2}}\int^{+a/2}_{-a/2}\int^{+a/2}_{-a/2}\frac{dx'~dy'}{|\vec{r}_{0}-\vec{r}_{i}(x',y')|}, $$

((4))

$\vec {r_{0}}$ is the position vector along the interface (i.e., the fault) where the displacement u _j is computed, and $\vec {r_{i}}$ is the distance vector centered around the contributing fiber i. The x ^′ and y ^′ vectors have their origin at the ith fiber. If we measure r in units of the discretization a, we can find a single parameter controlling the impact of the Green’s function. We call this, the reduced Young’s modulus [Stormo et al. (2013)]

$$ e \equiv \frac{a\cdot E}{L}. $$

((5))

In this way, we can re-scale the results of our calculations. The elastic response from the Green’s function is solely depending on the value of e, and thus changing the Young’s modulus, is equivalent to changing the discretization size a. Note that the system size is =a L. Accordingly, the reduced Young’s modulus can be written including the system size:

$$ e\equiv a^{2}E/{\cal{L}} $$

((6))

The model is quasi-static. Each fiber has a strength ${f_{j}^{t}}$ at which it breaks irreversibly if the applied force exceeds. The thresholds are following a random distribution. For the sake of this paper, the distribution is uniform. We impose the displacement D, calculate the forces on the fibers and find which fiber has the largest f/f ^t ratio. We then re-scale D so that this fiber is the only fiber with a force higher or equal to its threshold and break. We then redistribute the forces according to Equation 4 following an event-driven algorithm. This is done iteratively until all fibers of the system are broken. From this iterative process, it is possible to calculate the global stress Σ, as a function of the global displacement D of the system. The global stress is calculated for each iteration as the sum of force applied on the fibers divided by the area of the system:

$$ \Sigma=\frac{1}{a^{2} N}\sum_{i=1}^{N} f_{i} $$

((7))

We count time p as the ratio between broken fibers, n, and the total amount of fibers in the system N: p=n/N. A typical plot of Σ as a function of the damage is given in Figure 2. When Σ reaches its maximum, the system cannot support any additional load. We refer to the damage ratio at this point as the critical damage, p _c .

Figure 2

Typical load-displacement curve. Left: Typical Σ(n) force evolution with ‘time’ (i.e., number of broken fibers) in the elastic fiber bundle. Right: Zoom of force evolution indicated by the red area. Burst sizes δ are defined as the number of fibers broken before an increase of global load Σ is obtained. They are illustrated by horizontal red lines.

When a fiber breaks, a new equilibrium corresponding to a new u displacement field is computed, leading to a new load on the surviving fibers. Depending on the reduced Young’s modulus e, the redistribution of the load can be diffused or highly localized. The model contains two regimes of different behavior in this regard [Stormo et al. (2013)]. When the Young’s modulus is so high that the Green’s function is essentially turned off, the model is equivalent to the equal load sharing fiber-bundle model [Daniels (1945)]. We will call this regime the stiff regime. The system undergoes critical damage at $p_{c}^{\text {stiff}}$. As the system gets softer, a cross-over regime appears as the spatial influence described by the Green’s function begins to become relevant compared to the imposed displacement. At the end of this cross-over, we have the soft regime, where the Green’s function is completely turned on. This system undergoes critical damage at $p_{c}^{\text {soft}}<p_{c}^{\text {stiff}}$ [Stormo et al. (2013)].

There is no clear value of e separating the stable regimes from the cross-over, but we will use the values for stiff and soft systems that were used by Stormo et al. (2013), as they are safely outside the boundaries of the cross-over: e _stiff=32 and e _soft=7.63·10⁻⁶. The force felt by a fiber is proportional to u−D. For stiff systems, D≫u and the influence of u is negligible. For soft systems u∼D and the load distribution is highly dependant on the Green’s function. When we compare D and u for a system with no broken fibers, we get u=0.0162D for e=32 and u=0.9999D for e=7.63·10⁻⁶.

The ductile to brittle transition

There are several ways to define the difference between ductile and brittle damage. Just as Amitrano (2003), we will use the amount of inelastic, i.e., the non-linear deformation the system sustains before critical damage occurs as a measure of ductility. If there is instead a strongly linear load-displacement all the way until failure, we name it a brittle system. We then study the load-displacement curves from the stiff (e=32) and soft (e=7.63·10⁻⁶) systems shown in Figure 3. On top, we can see the load-displacement for stiff systems. It is clear that the load-displacement curve is highly non-linear until critical damage (indicated by the star). This corresponds to our definition of ductility. On the bottom, we see the same curve for a soft system. In contrast to the stiff system, there seems to be no sign of non-linear deformation before failure (indicated by the star). Thus, we have a strong indication of a ductile to brittle transition as we move from the stiff regime to the soft one.

Figure 3

The load Σ as a function of displacement D . Top: Load displacement for stiff, ductile systems (e=32). Bottom: Load displacement for soft, brittle systems (e=7.63·10⁻⁶). The point of critical damage is indicated by a star.

In order to see how the brittle-ductile transition affects the microscopic damage process in the two regimes, Figure 4 shows a map of damage at $p_{c}^{\text {soft}}$ and $p_{c}^{\text {stiff}}$ and the distribution of 〈Δ r ²〉^1/2. The red and the blue vertical lines in the 〈Δ r ²〉^1/2 map indicates $p_{c}^{\text {soft}}$ and $p_{c}^{\text {stiff}}$, respectively. The yellow line on top of the 〈Δ r ²〉^1/2 map shows the ensemble average |〈Δ r ²〉^1/2|. Figure 5 shows the same for a soft system.

Figure 4

The damage process of stiff ( e =32) systems. Top: Damage maps at $p_{c}^{\text {soft}}$ and $p_{c}^{\text {stiff}}$. The yellow areas indicate intact fibers, while the black ones indicate damaged fibers. Bottom: Distribution of 〈Δ r ²〉^1/2 as a function of fiber-damage n for L=128 (reprinted from Stormo et al. (2013) Ⓒ2012 American Physical Society). Darker colors indicate higher density. The yellow line shows |〈Δ r ²〉^1/2|.

Figure 5

The damage process of soft ( e =7 . 63·10 ⁻⁶ ) systems. Top: Damage maps at $p_{c}^{\text {soft}}$ and $p_{c}^{\text {stiff}}$. The yellow areas indicate intact fibers, while the black ones indicate damaged fibers. Bottom: Distribution of 〈Δ r ²〉^1/2 as a function of fiber-damage n for L=128 (reprinted from Stormo et al. (2013) Ⓒ2012 American Physical Society). Darker colors indicate higher density. The yellow line shows |〈Δ r ²〉^1/2|.

In the stiff systems (e=32), the localization of damage is independent from iteration to iteration. The average square distance between two consecutive broken fibers is $|\langle \Delta r^{2}\rangle ^{1/2}|=L/\sqrt {6}$ [Stormo et al. (2013)]. We can see from Figure 4 that the average is slightly above 50, consistent with $L/\sqrt {128}=52.26$. In soft systems (e _soft=7.63·10⁻⁶), the influence from the Green’s function ensures that the fibers closer to already broken fibers have a larger chance for failure, and thus |〈Δ r ²〉^1/2| dramatically decreases as seen in Figure 4.

The diffuse distribution of damage in the stiff systems and the highly localized damage of the soft systems led to describe the two regimes as ductile and brittle, respectively, see Anderson (2005).

Results and discussion

Burst from constant force steps

Damage burst definition

We now introduce a first definition of burst that is typically employed in quasi-static models in order to compare quasi-static time p to laboratory time. It is assumed that during a simulation, if a fiber breaks and the resulting load redistribution induces more fibers exceeding their force threshold f ^t, then the latter fibers are also damaged at the same time. Following the Σ(n)-curve in Figure 2, we can see local maxima where the subsequent fibers fail at a lower load. In laboratory time, we assume the number of fibers:

$$ \delta=\Delta n, $$

((8))

that break before Σ is increased, fail in an avalanche or burst.

Damage burst distribution, b _b

In order to compare our model to the one presented by Amitrano (2003), we want to find the distribution of temporal damage bursts. We can then extract a burst distribution exponent b _b -value from our system using the event-size distribution. We recall from Equation 8 that we define the size of an event δ as the number of fibers that are broken for each increase of the load Σ. In Figure 2, we see a zoomed version of Σ(p), and we can identify δ. Burst data ends at p _c , as it is the point of maximum Σ. Following Amitrano (2003), we assume that the size of events are comparable to the magnitude of acoustic emissions, and we thus get the δ- b _b connection:

$$ \log_{10}[\!N(\delta)]\propto -b \log_{10}[\!\delta], $$

((9))

where N(δ) is the number of bursts of size δ or higher.

We run the elastic fiber-bundle model for both the stiff (e=32) regime with diffuse ductile failure, and the soft (e=7.63·10⁻⁶) regime with localized brittle failure. This is done for both regimes at L= [ 32,64,128,256].

Figure 6 shows the burst distribution for stiff systems.

Figure 6

Burst distributions for stiff systems. The cumulative distribution of burst sizes for stiff (e=32) systems. The curves for different system sizes have been collapsed into a single curve. The best collapse has been obtained for α=1.6 and γ=1.2.

The slope of 1.475 for the linear trend for δ/L ^γ≤1 is consistent with the theoretical 1.5 burst exponent for the equal load sharing fiber bundle model (see Pradhan et al. (2010), Hemmer and Hansen (1992) and Pradhan et al. (2005)). This gives us a size-independent stiff b _b -value of 1.475 for L=[ 32,64,128,256]. We also observe a good collapse of the data for the different size systems using the collapsing method of Girard et al. (2010) (Figure 6) and obtained similar values of γ=1.2 and α=1.6.

In contrast from the stiff systems, the soft systems have size-dependent b _b -values (see Figure 7). The slope of the burst distribution gives b _b (L)= [ 2.02,2.23,2.60,3.00] for L=[ 32,64,128,256], with no clear sign of convergence. As changing L also influences the elastic behavior of the system from Equation 5, we also tested the sole influence of e on the slope of the burst distribution, keeping L constant and taking L=64 (see Figure 8). Results obtained in these simulations are less precise that those obtained for e=32 and e=7.63·10⁻⁶, but they highlight that the slope is controlled by the value of e. At high e the slope reaches the value 1.475 independently of e as observed previously while the b _b -value goes to higher values as e decreases. It confirms that e≈10⁻² is at the transition between both regimes.

Figure 7

Burst distributions for soft systems. The cumulative distribution of burst sizes for soft (e=7.63·10⁻⁶) systems. The black line is a guide to the eye for the given slope of the distribution. The b _b -values extracted are b _b (L)= [ 2.02,2.23,2.60,3.00].

Figure 8

b _b -values of damage burst distribution. b _b -values obtained from the slope of the cumulative distribution of damage burst sizes for various e moduli. The system size is set to L=64. Each point is an average over 400 simulations.

Burst from thresholding the local fracture front velocity

Velocity burst definition

Another approach to characterize the burst activity of the model for a comparison with earthquake seismicity is to follow the analysis proposed by Grob et al. (2009). Here, we extend this approach from a single front monitoring, initially proposed in Måløy et al. (2006) to a full interface with multiple active fronts at the same time. In order to first calculate the velocity of the fracture front, we use the concept of waiting time matrix, W, introduced by Måløy et al. (2006). The waiting time matrix, W, is a matrix with dimension L×L, storing the number of times a fiber is part of a fracture front. Initially, every element of W is set to 0. At every time step, the positions of the fracture fronts are identified, and the corresponding elements of W are increased by one.

We define any surviving fiber that neighbors a damaged one as a part of an evolving fracture front.

When the fracture is complete, we create the velocity matrix by taking the inverse of the elements of W times the discretization length:

$$ V(x,y) = \frac{a}{W(x,y)}. $$

((10))

where V(x,y) is the velocity at which the fracture front passes trough element (x,y).

Finally, in order to identify velocity bursts, we introduce the threshold level C as in Måløy et al. (2006). We first calculate the average velocity of the fracture front, 〈V〉. We then go through the velocity matrix and create the velocity cluster map by setting every element of V(x,y) above C·〈V〉 to 1 and the rest to 0. We then go trough this cluster map and identify the cluster size distribution using a Hoshen-Kopelman algorithm [Hoshen et al. (1978)]. By doing this, we can identify the cluster size distribution p(s,C), where s is the size of any given cluster. Example of velocity map and velocity cluster map are given in Figure 9 for stiff systems and in Figure 10 for soft systems. We see that clustering emerges for soft systems with the highest velocities.

Figure 9

Velocity map for stiff systems. L=128, left: Map of log10[ V(x,y)] for a stiff system. The average velocity is 〈V〉=0.000147. Right: The corresponding C=2 cluster-map. We identify the velocity cluster by setting all elements above C·〈V〉 to one (yellow) and below to zero (black).

Figure 10

Velocity map soft systems. L=128, left: Map of log10[ V(x,y)] for a soft system. The average velocity is 〈V〉=0.000187. Right: The corresponding C=2 cluster-map. We identify the velocity cluster by setting all elements above C·〈V〉 to one (yellow) and below to zero (black).

Velocity burst distribution, b _v

We calculate the cluster distribution p(s,C) for both regimes. This distribution is given in Figure 11 for stiff systems and in Figure 12 for soft systems.

Figure 11

Distribution of velocity clusters for stiff systems. The velocity cluster distribution p(s,C) for L=128, C= [ 1,2,3,5,10]. The black 5/3 slope is added for comparison.

Figure 12

Distribution of velocity clusters for soft systems. The velocity cluster distribution p(s,C) for L=128, C= [ 1,2,3,5,10]. The black 5/3 slope is added for comparison.

In order to compare our model to seismic activity, we will follow the analysis of Grob et al. (2009). A link is made between the area of the velocity clusters and the seismic moment M ₀ following Kanamori and Anderson (1975):

$$ M_{0}=\mu s d. $$

((11))

where μ is the shear modulus and d is the slip offset. Since our model does not undergo slip but only mode I opening, we assume d to be constant. Accordingly, the probability density for seismic moment follows the power law:

$$ p(M_{0})\propto M_{0}^{-1-B}, $$

((12))

where B has been found to be stable around 2/3 for a large number of seismic regions [Scholz CH (2002)].

In both Figures 11 and 12, we show a 5/3=2/3+1 slope for comparison. We can see that for small C-values (pink color), the cluster distribution is consistent with B=2/3 for both regimes in particular for small events. From Scholz CH (2002), we get the connection to the Gutenberg-Richter distribution, as b=3/2B. We thus obtain for the front advance distribution, b _v =1 for both ductile and brittle systems. However, a clear finite size effect is visible for soft systems and large events. As C becomes larger, the difference between both stiff and soft regimes increases. In particular, we obtain an artificially large slope (large b-values) for stiff systems (Figure 11) for large C-values. This is not the case for soft systems for which the C-value only changes the pre-factor (i.e., a-value) but not the general behavior of the distribution. This observation suggests that the change in b-value can be related to the a posteriori clipping level and not to a physical parameter of the model. In this case, b-value changes are more related to recording instrumentation rather than different physical processes.

When we look at the velocity clusters, we can see that there is a possible b _v ≈1 value for both ductile and brittle systems. However, for the system sizes currently available, the large clusters get rare as we increase C. This happens faster for the stiff systems, as the velocity distribution is narrower (Figure 9) than for the soft systems (Figure 10). With this result, we extend to a new configuration (multiple fronts), the results obtained for a single directional front propagation reported by both Måløy et al. (2006), Grob et al. (2009), and Tallakstad et al. (2011).

The C-level is the minimum local crack velocity that we have in a velocity cluster. As C increase, the average local velocity of the clusters is larger. If we compare to shear rupture, an increase of C is analogous to a rising of the average slip velocity. In that sense, it will select only events with large velocities or amplitudes if measured with a velocimeter. This assumes a direct link between rupture velocity along the fault and wave magnitudes at the seismometer. Another discrepancy with our model and the seismic activity is the fact that our model is quasi-static and not dynamic (i.e., no dynamical rupture propagation). Thus, it is not straightforward to compare with dynamic rupture and dynamic effects. In order to get a firmer connection to the seismic activity, the model has to be further developed.

Our model displays the same type of change in b _b -value during ductile to brittle transition, i.e., a transition from diffuse to localized damage in the system, as described by Amitrano (2003). In the brittle regime, we get b _b -values that are twice the b-values observed during the hydraulic injections in the EGS experiments (average of 1.58 during injections, 1.15 after shut-in at Basel, Bachmann et al. (2011)) (average of 1.29, ranging from 1.0 to 1.7 in the GPK2 site in Soultz-sous-Forêts, [Cuenot et al. (2008)]).

In the analysis of Cuenot et al. (2008), they found that when they increased the injection rate from 40 to 50 ls ⁻¹ they initially saw a drop in the b-value. Eventually it recovered, before it climbed to a new peak. They speculate that this seemed to come from the fact that the weak points in the fault were already damaged, and that the strong points had to break for the water to seep trough to the rest of the network. Breaking the strong points would release more energy giving an initial drop in the b-value. When the fluids finally were able to seep through to the rest of the network, a lot more weak points were available. Thus, a significant increase in b was to be expected.

In our model, we see a similar increase in the b _b -value when we increase the number of fibers in the system (L). In the fiber bundle model, the weakest bonds have the highest risk of failure. As the system becomes larger, the average distance between the weak points increases due to the uncorrelated nature of the threshold distribution. When we compare small systems with fewer fibers, the build-up of stress concentration will be much denser. Thus, the chance of a failure triggering an avalanche nearby is higher. This explains why our b _b -value increases with system size. This could be linked to the Cuenot’s argument for the b-value fluctuation.

The model does come with its share of limitations though. First of all, we do not get the exact same b _b -values as seen in the hydraulic injections. And even though we can see a change of b _b -value, the transition described by Cuenot et al. (2008) cannot at this stage be addressed if the change comes from the opening of the fracture network, as our model is dealing with a single fault. Connecting different samples and activating them as the fluid spreads through the system might be a way of addressing this issue.

Secondly, as illustrated in Figure 1, the load imposed on our fiber bundle is external. This is in contrast to the internal pressure experienced by the faults during the hydraulic injections. Rules for applying the development of hydro-pressure in the weak plane between the elastic half-spaces should be developed to account for this. This should not change the burst distribution in the stiff regime as it is controlled by the disorder, but might have an effect on the soft systems.

A third limitation of our model in connection to the EGS experiments, is the fact that the model is describing a pure mode I loading. This is in contrast to the typical slip faulting in earthquakes. Fortunately, the energy released for mode I, II, and III fracture only differ by a prefactor as described by Schmittbuhl et al. (2003). Thus, the transfer of load and the shape of the damage burst distribution should be applicable to those cases as well.

In our model, the length a (see Equation 5) refers to the distance between two neighboring fibers. In this sense, the quantity, 1/a can be understood as the number of fiber per unit length (or as a fiber linear density). Then, for a fixed system size, /a=L represents the space density of fibers which can be considered as the density of contact points between the facing surfaces of the fracture. In a geological context these contact points are referred as asperities along the fault plane. Changing the asperity density L along the fault plane (i.e., the heterogeneity of the interface), implies changing the important parameter of the model: e=(a·E)/L (from Equation 5). When e is small, the system is considered as soft, localization of the plastic deformation (i.e., fracture) is observed and b-values are untypical. When e is large, the system is stiff, deformation is diffused and b-values are more typical. It is of interest to see the Young E might have a inverse effect on e compared to the asperity density L. Elastically stiff systems (i.e., large E) are similar to a small density of asperities L in particular if the asperity size gets close to the system size: a→. Stiff fiber-bundle models can then be viewed as systems with large scale heterogeneities along the interface (i.e., large a and small L) for which the damage zone is notably extended.

The main difficulty to relate our model to natural fault system is to assess the degree of heterogeneity of the interface at depth event in the context of geothermal reservoir. If we were able to resolve such properly, we could then make a direct comparison between the b-values encountered in various interfaces with the results of our model. Although we have limited knowledge of the fault heterogeneity, we propose two observations that is linked to this feature and would help to draw some link between our model and natural reservoir. First, it has been shown that the contact area along a rough fault, is typically related linearly, at first order, to the normal stress along the fault [Batrouni et al. (2002a), Hyun et al. (2004)]. Accordingly, the asperity size a will tend to the system size as normal stress is increased. The system is then expected to become effectively stiffer (i.e., increase of e) as the normal stress increases along the fault. If so the b-value is expected to decrease the normal stress increases (see Figure 8). Second, one could use the out of plane distribution of seismicity as a function of the fault thickness. As shown for faults in California, this distribution can be related to the stress variations owing to the fault roughness and thus can be interpreted as a quantification of the heterogeneity of the interface [Powers and Jordan (2010)]. In this sense, it would then be directly related to the L parameter introduced in our system. Provided that earthquake locations are of sufficient accuracy to estimate precisely such a distribution, it would then make a link between the parameters of our model and observable quantities.

Conclusions

We have found that the elastic fiber bundle model shows promising features to understand the mechanical origin of the b-value fluctuations during stimulation of deep geothermal reservoirs. We proposed two types of events in our model that can be compared to seismic events: a classical measure, the damage burst, i.e., the number of broken fibers during an elementary rise of the loading stress, and a new measure, the velocity burst, i.e., the size of a front velocity cluster during its propagation. We show that both event definitions lead to different event distributions. The damage event distribution has a b _b value that is influenced by the normalized e modulus: changing at small e (soft systems) but constant at large e (stiff systems). The second definition based on the front velocity produces distributions with a b _v value that is evolving differently. We have shown the feasibility of connecting the velocity cluster distribution to the distribution of seismic moment M ₀ of earthquakes. The main change of the b _v -value is coming from the clipping level and not from the transition from stiff to soft (or ductile to brittle). It is then related to the observation sensitivity. It suggests that b-value changes could be related to measurement evolution.