Hydroshearing poorly connected preexisting fractures in the presence of stress anisotropy as a random percolation process

Activation of natural fractures by ﬂuid injection is used to enhance the permeability of geological reservoirs. A fracture slips and activates when the ﬂuid pressure is larger than a critical value that depends on its orientation with respect to the in situ stress. Using large-cell Monte Carlo renormalization group, we show that activating poorly connected fractures belongs to the same universality class as random percolation despite the propensity of stress anisotropy to activate favorably oriented fractures. A crossover that does not change the universality class is identiﬁed as the size of the activated network exceeds the preexisting fractures’ correlation length. DOI


I. INTRODUCTION
Hydroshearing, the enhancement of rock's permeability by fluid injection to activate preexisting fractures, is used to extract natural gas from shale formations and geothermal energy from hot dry rocks [1][2][3][4].Preexisting fractures are fractal and a crucial requirement for a successful hydroshearing process is to create a network of convective transport paths that fill space effectively.Anisotropy in the overburden stresses favors activation of fractures oriented in a preferred direction, possibly leading to the creation of a low-dimensional network.In this work, we show using a large-cell Monte Carlo renormalization group (RG) approach [5] that the activation of preexisting fractures is a random percolation process in spite of stress anisotropy and long-range correlations due to the fractal nature of the preexisting fractures [6].Building on this result, we present constitutive laws that predict the effects of fracturing process conditions on the permeability evolution of stimulated rocks.
Preexisting fractures activate by shearing when the frictional force holding the fractures' surfaces in place is reduced [7,8].The induced mismatch, due to slippage, between the asperities of the two surfaces provides space for hydraulic flow and the fractures become activated [9].The fractures are typically modeled as planar objects and in two dimensions, Mohr-Coulomb's criterion yields a critical fluid pressure, P c , required to induce slippage [10,11]: Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
where σ 1 and σ 3 are the far-field maximum and minimum principal stresses, respectively, η is the friction coefficient, and θ is the angle of the fracture with respect to the 3 axis.Following Mohr-Coulomb's criterion, an inherently connected network of hydraulically conductive activated fractures can form when the fractures are slipped by the arrival of an injected fluid with a viscosity higher than that of the reservoir fluid.Due to the percolative nature of this mechanism, the connectivity of the formed cluster of activated fractures is predominantly governed by (1) the connectivity of the underlying network of preexisting fractures, i.e., its correlation length ξ 0 ; (2) the variability in the fractures' critical pressures which is proportional to σ 1 − σ 3 ; and (3) the injected fluid's viscous pressure drop over the fractures' length, l, controlled by the fluid's viscosity, the injection rate, and the fractures' aperture.
When the viscous pressure drop is small compared to the variability in the critical pressures, the fluid activates and flows through the least resistant accessible fractures.Each newly activated fracture is the one with the smallest value of P c , calculated from (1), that is connected to the network of activated fractures.Thus, the primary source of randomness is the interconnections between fractures associated with their random positions and orientations.
To provide a framework for understanding the hydraulic stimulation of rocks and the transport processes after stimulation, we will determine the universality class of the hydroshearing process.On first thought, one might expect anisotropy in the stress field to produce a quasi-one-dimensional network when percolating paths through the preferentially oriented fractures are available.Upon closer examination, one may realize that the fracturing process resembles site-bond percolation in which a certain number of sites on a lattice are filled with particles that are available for bonding.A sparse network of natural fractures would then correspond to the case where these particles are near their percolation threshold.Most site-bond percolation problems fall within the regular percolation universality class including those in which the probability of bond activation between well-connected sites is anisotropic and those with isotropic bond activation among FIG. 1.An example of a sparse network of preexisting fractures.The red fractures are the activated fractures that form a percolating network, while the blue fractures are the unactivated preexisting fractures.
sparse sites.One may then infer that hydroshearing is in the percolation universality class.However, this postulate is uncertain because no previous studies have considered the closest analogy to hydroshearing among site-bond percolation problems, i.e., one in which anisotropy persists and evolves as the number of sites approaches the percolation threshold.Furthermore, unlike all bond percolation problems, hydroshearing is a deterministic process in which all fractures with a given orientation are opened at the same injection pressure.
In this paper, we will use RG to show that, indeed, the hydroshearing process under the aforementioned assumptions belongs to the same universality class as random percolation.First, we will present the fracture network model used to represent the preexisting fractures.Then, we will describe the renomalization group method by which we calculate the universal exponents that determine the problem's universality class.Finally, the implications of the results for predicting the permeability of fractured reservoirs are discussed.

II. NETWORK MODEL DESCRIPTION
For simplicity, we examine the activation of a twodimensional network of preexisting fractures with equal lengths, l, and isotropic orientations as illustrated in Fig. 1.The connectivity of a network of equal-length fractures with isotropic or anisotropic orientations has been shown to belong to a random percolation universality class [12][13][14].Its correlation length, ξ 0 , is given by ξ 0 = k ξ (n − n c ) −ν n where n is the number density of preexisting fractures, n c is the percolation threshold number density, and k ξ is a proportionality constant.The value of n c , for isotropic preexisting fractures, has been estimated to be about 5.637 [15,16] and the value of ν n is 4/3.When 0 < n − n c n c , ξ 0 is much larger than l but is finite.Thus, the network looks fractal at length scales smaller than ξ 0 where long-range correlations determine fluid propagation, but it is homogeneous and stress anisotropy controls the activation process on scales larger than ξ 0 .This means the hydroshearing behavior is different at these different scales.
The lengths of natural fractures, typically, follow a powerlaw distribution, P(l ) ∼ l −a .The fractal dimension of the preexisting fractures depends on the spatial correlations between the fractures [17].Our model is applicable to geological formations such as the Gulf of Suez [18] with power-law exponents, a, greater than or equal to 3 whose connectivity is dominated by the smallest fractures.The fractal dimension of the random fracture model used in this paper is also believed to characterize systems of natural fractures such as The Geysers geothermal field [19].Regardless of the fractal dimension of the preexisting fractures, the main result that the stress anisotropy does not control the fracturing process of fractal networks should apply as long as the lengths of the fractures are of the same order of magnitude.
The procedure used to prepare the network of preexisting fractures is similar to that described in [15].To ensure the formation of a homogeneous network, the fractures' centers of mass are distributed homogeneously within a box of size (L + 1) × (L + 1).Free boundary conditions are used to define the system and fractures near the edge of a box of size L × L are truncated to prevent intersections of the fractures outside the domain and thus the formation of a percolating path outside the box.Under such conditions, the percolation behavior of the network of preexisting fractures is controlled by their dimensionless number density defined as where l i is the length of each fracture and N f is the total number of fractures.Given a specified number density of the fractures, fractures are added until the total length squared of the fractures after truncating the ones near the boundaries is equal to or larger than nL 2 .The maximum error in the number density scales as 1/L 2 .For systems larger than 100, the maximum deviation from the set number density was found to be less than 1%.The results presented in the paper correspond to system sizes that are larger than 100 × 100.
The critical pressures of the fractures are calculated using (1).By normalizing the critical pressures with their standard deviation, becomes where . σ is the normalized average of the fractures' critical pressures.The system is defined by two dimensionless parameters, η and σ .One can easily show that the average value of the critical pressures does not change the percolation behavior of the activated fractures.As the value of the friction coefficient, η, decreases, the anisotropy of the distribution of critical pressures increases.We found that varying the friction coefficient does not change the qualitative behavior of the activation process but slightly alters the value of the onset of percolation when the effects of stress anisotropy are pronounced, i.e., the natural fractures form a well-connected network.The presented results are based on the value η = 0.85, which is close to the experimental friction coefficients of most rocks.

III. RENORMALIZATION GROUP APPROACH
The universality class of the activation process will be identified by determining the value of the activated fractures' fractal dimension, D f , and the critical exponent ν.The exponent ν is related to the scaling of the correlation length, ξ , of the percolating network of activated fractures with the distance from the percolation threshold value.The percolation threshold value is defined as the fluid pressure required to activate a fraction F ∞ c of the percolating preexisting fractures and form an incipient infinite cluster (IIC) of the activated fractures.Since generating an infinite cluster using computer simulations is not viable, we use RG [20] to analyze the formation of percolating networks that span finite domains as described next.
RG is based on the premise that the IIC is self-similar and thus when renormalized with cells of a linear dimension, L, the correlation length of the renormalized IIC, ξ , is simply equal to that of the original network, ξ , scaled with L. Since ξ is a function of the probability to activate the renormalizing cell and the percolation threshold value does not change upon renormalization, one can compute the values of ν and F ∞ c once the dependence of cell activation on fracturing preexisting fractures contained in the cell is determined.A cell is considered activated when the preexisting fractures contained in the cell form a percolating path and the dependence of cell activation on fracturing of preexisting fractures is found through computer simulations.
To get a glimpse of the RG procedure, consider a homogeneous network of preexisting fractures.The correlation length of a network of activated fractures is given by ξ ∼ (F − F c ) −ν where F is the probability to activate a fracture and F c is the percolation threshold value.Replace a finite domain, of a linear size, L, with a supercell that has a probability to percolate, F .F is a function of the probability to activate the fractures, F , within the domain.The correlation length of the normalized network composed of supercells, ξ , is given by ξ ∼ (F − F c ) −ν .Since the connectivity of the network is reduced by L, one can relate ξ and ξ through ξ = Lξ .Since the value of F c and ν are the same for the original network and the normalized one, it can be easily shown that where * = dF dF at F = F c .Finally, the self-similarity of the network of activated fractures at all length scales is only valid when F = F = F c .Thus, one needs to find the relationship between F and F to calculate the value of F c .Using (4), one can calculate the value of the exponent ν after determining the value of F c .For example, in a triangular lattice of sites, one can replace every three sites that form a triangle with a supersite.The supersite percolates if all of the sites or two of them are activated.Thus, the closed form of 2 and the value of ν using (4) where L = √ 3, since 3 sites are used to form a supercell, is ν = 1.355 [20].
The procedure described above is only guaranteed to produce an exact result for F c and ν when L → ∞.Asymptotically, the value of ν calculated using a finite size of the supercell approaches the exact one.As the size of L is increased, it becomes increasingly difficult to find a closed form of F = f (F ).Thus, a Monte Carlo approach is used to sample the total probability to percolate the supercell.One way to sample the distribution is to perform several percolation simulations at different values of F to get the functional form of F [21].One should notice that F = f (F ) is the cumulative probability to percolate a supercell and thus * = dF dF is the probability distribution function of percolating the supercell evaluated at F c .Thus, one can perform percolation simulations using Newman and Ziff's algorithm [22] to calculate the distribution, , of the onset of percolation and evaluate * which reduces the number of simulations required to find Once the value of * is evaluated for different cell dimensions, the value of ν is given by [20] 1 while the value of Here, ν L and F c L are the cell-size-dependent values of the critical exponent, ν, and the percolation threshold, respectively.In fact, 1/ν L = ln( * ) ln(L) .In Newman and Ziff's algorithm, in the context of the activation process, preexisting fractures are successively activated starting from the ones with the lowest critical pressure until a percolating network is formed.The fraction of activated preexisting fractures, F c , which can be used as a percolation parameter, is recorded once a percolating network is formed.Using several realizations, one can find the distribution, , of the onset of percolation.By introducing an assumption about the functional form of the probability distribution, , of the onset of percolation, one only needs to store the mean, f c , and variance, σ f c , of the onset of percolation to evaluate ν [5].
As the system size increases, approaches a delta function and F c is the mode of the distribution.For a finite system, F c is expected to lie close to the mode of the distribution.This assumption has been found to produce a close value to the correct percolation threshold in other systems [5].Typically, a Gaussian distribution of f c is assumed and thus * scales with the inverse of the standard deviation of the distribution, 1/ √ σ f c .Since the probability to percolate is bounded between 0 and 1 and its distribution is skewed, we have used a β distribution as it produced a better fit of the distribution of the onset percolation as shown in Fig. 2. The β distribution is given by where the peak of the distribution and thus * can be related to f c and σ f c .In fact,

A. Universality class of the activation process
Since the sparse network of natural fractures is fractal when L ξ 0 while it is homogeneous when L ξ 0 , a crossover in the scaling of ν L with L is expected as L/ξ 0 ∼ O(1).To probe this crossover, using reasonable cell sizes, one can rescale the system size in (6) Here, ν is the calculated exponent related to the connectivity of the IIC in either the fractal or homogeneous regimes.ν n is the critical exponent related to the connectivity of the preexisting fractures that is equal to 4/3 [12].If ν = ν n , the ln(n − n c ) terms on the two sides of the equation cancel retrieving (6).The value of the constant C should not be a function of n when the fractal and homogeneous regimes are fully developed but it is different in the two regimes since the activation behavior is different in these regimes.Thus, a universal curve for all number densities that satisfy 0 < n − n c n c will be obtained when plotting ) for different number densities and a range of cell sizes.As expected, a crossover in the scaling of ν L occurs when L(n − n c ) −ν n = O (1).The slope of both the dash-dotted and solid lines is equal to 3/4 indicating that the value of ν is equal to that for percolation problems.The value of the proportionality constant, C, in both regimes is given by the intercept of the fit lines and is equal to 0.255 and −0.42 in the fractal and homogeneous regimes, respectively.FIG. 3.This plot is used to calculate the value of ν.The values of n range between 5.638 and 5.9 while 100 L 800.Each value of n is represented by a different color.The slope of the two lines is equal to 3/4 while the intercept of the dashed-dotted line is equal to 0.255 and the intercept of the solid line is equal to −0.42.
By estimating the volume of the activated fractures from core samples, one can employ the scaling of the cluster mass, M, to estimate the radius of the stimulated reservoir.The cluster mass is defined as the total length, scaled with λ, of the activated fractures.Since the cluster mass scales with the system size as Here, V inj is the injected volume and V f and λ are the average volume and length of the activated fractures.To determine the fractal dimension in the two regimes, the scaled total length, M, of the activated fractures is calculated for different domain sizes of the IIC. Figure 4 shows a plot of the cluster mass, M, scaled with the number density of preexisting fractures for different system sizes.Unlike the scaling of ν L , no transition is observed as the system size becomes of the same order of magnitude as the preexisting fractures correlation length.The slope of the solid line is 91/48 indicating that the fractal dimension in the two regimes is the same as the fractal dimension of a FIG. 4. This plot is used to calculate the value of D f .The values of n range between 5.638 and 5.9 while 100 L 600.The coloring follows Fig. 3.The slope of the solid line is 91/48 and the intercept is equal to 1.33.percolating network formed using random percolation.The intercept of the fit line, I d = 1.33, can be used to calculate the proportionality constant, k D .Since Fig. 4 is generated by rescaling the system size with the preexisting fractures correlation length, As discussed in Appendix, the value of k ξ is found to be equal to 0.43 and thus, k D = 5.35.If one fits the cluster mass as a function of the different number density of fractures without rescaling of the system size with ξ 0 , the average value of k D is found to be equal to 4.95 ± 0.71 in good agreement with the estimation using the calculated value of k ξ .
Results from Figs. 3 and 4 show that the effects of stress anisotropy are irrelevant details due to the large population of singly connected fractures (known as red links in percolation theory) the fluid has to activate in order to propagate through the rock.Moreover, the activation process is still in the same universality class as random percolation even when the activation process is controlled by topological long-range correlations, due to the fractal nature of the network of preexisting fractures.This result is consistent with the Harris criterion which states that the universality class of a process does not change upon dilution when dν > 2 where d is the Euclidean dimension of the network [23].It is worth mentioning that this criterion is not rigorously proven and there are reported systems that violate it [24].
For an optimized fracturing process, the stimulated region should be of the same size as the hydrocarbon-rich zone.Since the radius of the stimulated region is given by R/λ = 0.42(V inj /V f ) 48/91 , one can approximate the required amount of fluid to be injected prior to commencing the fracturing process.For instance, if the radius of the payzone is of the order of 1 km, the average length of the fractures is equal to 1 m, and the volume of the activated fractures is 5 × 10 −3 m 3 , the total volume to be injected to stimulate a region of the same size as the payzone is approximately equal to 2000 m 3 of a low-viscosity fluid.Interestingly, the stimulated radius scaling, when the viscous pressure drop is negligible, is universal for rocks with uniform preexisting fractures.That is, the required injected volume is only a function of the volume of the fractures and is independent of their number density and connectivity.

B. Transport properties of fractured rocks
Having determined the universality class of the activation process, let us relate the transport properties of the fractured rock during the activation process to the injection fluid pressure, P inj , and the rock's statistical properties.Such relations can help simplify the design of the fracturing process and predict the permeability after stimulation.Moreover, relating the effective permeability of the stimulated rock to the pressure of the injected fluid, as argued in [25], can allow for continuum modeling of the stimulation process when the viscous pressure drop becomes important at large length scales.
When the viscous pressure drop is negligible, the permeability, K, and porosity, φS, of the hydraulically stimulated rock, should scale as K ∼ [(P ∞ cc − P inj )/(σ 1 − σ 3 )] and φS ∼ [(P ∞ cc − P inj )/(σ 1 − σ 3 )] β , respectively.Here, P ∞ cc is the threshold fluid pressure required to form an infinite percolating network.It corresponds to the inverse of the cumulative FIG. 5.This plot shows the scaling of the saturation, S, with the system size at the onset of percolation.When L ξ 0 , S is independent of L and it scales as S ∼ L −β/ν when ξ 0 L. The slope of the solid line is equal to −β/ν = −5/48.The coloring follows Fig. 3. distribution of the critical pressures, F , evaluated at F ∞ c .The form of F for randomly oriented fractures is given by where ), where θ min is the orientation of the fracture with the lowest critical pressure.The preferred orientation for fracture slippage is given by θ min = ±atan(η + η 2 + 1) + mπ (for any integer m) (9) is obtained by integrating (3) for θ from −π/2 to π/2 normalized by π .φ ∼ ξ −β/ν 0 is the porosity of the medium if all the percolating preexisting fractures are activated, while S is the fraction of preexisting fractures that belong to the percolating network of activated fractures.The values of the critical exponents , β are universal and in two dimensions equal to 1.3 and 5/36, respectively [20,21].However, the value of P ∞ cc depends on the connectivity of the network and the stress anisotropy.Hence, estimating the value of P ∞ cc using field data prior to stimulation helps to predict how the permeability will evolve as the rock is stimulated under negligible viscous pressure drop.
As a sample verification of the derived relations, we probed the scaling of the saturation of activated fractures, S, with the cell size.At the onset of percolation, the saturation should scale as S ∼ (L/ξ 0 ) −β/ν in the homogeneous regime while it becomes a constant in the fractal regime.As shown in Fig. 5, the crossover between the two behaviors occurs as L/ξ 0 ∼ O (1).Derivation of these scalings stems from the fact that the correlation length of the activated fractures, ξ , scales as ξ ∼ [(P ∞ cc − P inj )/(σ 1 − σ 3 )] −ν .At the onset of percolation, ξ is of the same order as the system size and when n = n c , ξ ∼ ξ 0 since the preexisting fractures form a quasi-one-dimensional network [26].Now, let us find the dependence of the percolating pressure P ∞ cc on the natural fractures' number density.Depending on the ratio of the correlation length of the preexisting fractures to the system size, the finite-size scaling of F c is different.
where k H is a proportionality constant.This is the typical finite-size scaling used for percolation problems on perfect lattices.Since the activation problem belongs to the same universality class as random percolation, the percolation behavior on a homogeneous network of preexisting fractures is expected to be similar to that on a perfect lattice.
The percolating network of preexisting fractures when n → n c has the same fractal dimension as a percolating network formed at the percolation threshold value on a lattice.For an infinite percolating network, there exists only one path connecting two points within the network and thus F ∞ c should be equal to 1 in order to percolate the network [26].Thus, the scaling of F c with the system size is expected to follow when ξ 0 L. Similar to the scaling of ν L with the system size where the value of the proportionality constant, C, is different in the fractal and homogeneous regimes, the values of k f and k H are not necessarily equal.
If one calculates the size-dependent percolation threshold, F c , for an arbitrary number density such that n − n c n c , the scaling of F c with the system size is expected to follow (11) for relatively small systems and as one increases the system size a transition occurs and F c starts to scale as (10).Assuming that the transition in the scaling of F c from (11) to (10) occurs sharply when L = O(ξ 0 ), one can replace L with ξ 0 = k ξ (n − n c ) −ν n in ( 11) and (10) and equate F c in the two equations to obtain Since both the mean and the mode of the onset of percolation follow the same finite-size scaling and approach the same value as L → ∞, one can choose either value to calculate F ∞ c [5].In our simulation, we only store the average, f c , and variance, σ f c , of the onset of percolation.To calculate the mode, we had to assume a functional form of the distribution of the onset of percolation.To avoid unnecessary assumptions, we base our calculations of F ∞ c on the size-dependent values of f c .
To realize the transition in the scaling of F c for a given number density such that n − n c n c , one needs to simulate very large system sizes.Since that is not a viable option using current computational power, we showed the transition by means of calculating F c for different number densities and a range of system sizes.By substituting (12) into (10), one can show that Therefore, plotting (1 − F c )L 1/ν n vs L 1/ν n (n − n c ) will yield a universal curve as shown in Fig. 6.Since the scaling of F c in the fractal regime is given by (11) where the value of k f is independent of the value of n, a horizontal line given FIG. 6.This plot shows how F c L scales with L in the two regimes.The slope of the dashed-dotted line is given by (k f − k H )k −1/νn ξ and is found to be equal to 0.161.The coloring follows Fig. 3.
by y = k f should be obtained when (n 1, the finite-size scaling should follow ( 13) and the transition between the two scalings occurs when (n − n c )L 1/ν n ∼ O(1).In Fig. 6, the solid line represents (11) while the dashed-dotted line represents (13).
The proportionality constants in (12) are universal for all rocks when the length of the preexisting fractures are of the same order.Thus, one only needs to measure the number density of preexisting fractures from core samples to determine To verify (12), the predicted values of F ∞ c are compared with those calculated from computer simulations in Fig. 7.The circles are for the case when anisotropy effects are taken into account.The squares represent the case where the critical pressures are random such that the activation process is analogous to random percolation.By fitting the values of F c L , shown in Fig. 6, using the finite-size scaling arguments, the values of the proportionality constants are estimated to be k f ≈ 0.57, FIG. 7.This plot shows how the percolation threshold for infinite systems depends on the number density of natural fractures.The squares represent the case where P c 's are random while the circles are for the case where P c 's are given by (1).FIG. 8.This plot shows how the minimum injection pressure required to percolate the rock depends on the rock's properties.The units of the principal stress field is MPa.The dashed-dotted lines are calculated from (12), while the circles are calculated from the simulated data of F ∞ c for the case where the anisotropic effects are taken into account.The asterisk symbols represent the case where the critical pressures are random and have the same range as in the case where P c 's are calculated from (3).
k H ≈ 0.19, and k ξ ≈ 0.43.The methods used to estimate the constants are discussed in Appendix.
The curve created by the circles in Fig. 7 represents the transition from conditions yielding a finite cluster of activated fractures to those producing an infinite percolating network.As shown in the figure, (12) holds when n − n c n c , the regime where the network is sparse, i.e., ξ 0 l.Deviations from (12) for larger number densities of the preexisting fractures are due to the invalidity of the scaling of ξ 0 far from threshold.Since the fracturing process belongs to the same universality class as random percolation, an infinitely large network of activated fractures is expected to be isotropic [14].That is, the percolation threshold value is independent of the percolation direction even though the fluid propagates through the preferably oriented fractures.As indicated by the red circles in Fig. 7, it is easier for the injected fluid to percolate through the fractures in the presence of an anisotropic stress field when compared to the case where the stress field is highly heterogeneous, i.e., the critical pressures are random.When n − n c n c , the percolation threshold value is not affected by the stress anisotropy due to the large population of singly connected fractures the fluid has to activate in order to propagate through the rock.
Figure 7 is analogous to the Griffith phase diagram developed for the dilute Ising model to predict the ferromagnetic Curie temperature given the concentration of impurities [27,28].If the stimulation process is conducted at a constant pressure, there exists a minimum injection pressure, P ∞ cc , needed to ensure the continuous propagation of the injected fluid.Estimating this value prior to stimulation is essential in order to ensure the success of the process.Figure 8 shows example values of P ∞ cc calculated using (9) for different stress conditions and number densities required to form a percolating network of activated fractures.The circles use the simulation data shown in Fig. 7 while the dashed-dotted lines were calculated using (12).As shown in the figure, the injection pressure required to form a percolating network decreases with increasing connectivity of the preexisting fractures.For well-connected fractures, the injected fluid can easily find paths with low resistance to propagate through the rock.
Now, let us discuss how the evolution of the injection pressure can be predicted using easily measurable quantities.If the cumulative distribution function of the critical pressures is linearized, it can be shown that ( The value of k P depends on the connectivity of the preexisting fractures and the orientation distribution of the preexisting fractures.For randomly oriented sparse fractures, it is equal to the derivative of (9) evaluated at p c = P ∞ cc /δ pc .By replacing the system size in ( 11) and ( 10) with the scaling of the stimulated reservoir radius, it can be easily shown that the injection pressure evolves as P inj = σ 1 − k f k P δ pc (V f k D /V inj ) 1/(νD f ) at initial times when R ξ 0 and changes with the total amount injected as P inj = P ∞ cc − k H k P δ pc (V f k D /V inj ) 1/(νD f ) as the cluster radius exceeds the value of ξ 0 .Because the injected fluid propagates away from the injection well, the injection pressure required to drive the injected fluid increases as the total injected volume increases.Moreover, higher stress anisotropy allows lower injection pressures because the injected fluid percolates through lowresistance fractures.Interestingly, one can probe the connectivity of the preexisting fractures through analysis of the transient injection pressure profile.If one plots the injection pressure as a function of the total injected amount, a critical injected volume, V c inj , at which the two scalings intersect can be identified.The correlation length of the preexisting fractures is of the same order of magnitude as (V c inj /V f k D ) 1/D f .By solving for V c inj from the derived scaling of the injection pressure, one can show that where the numerical prefactor is given by (k f − k H ) −D f ν /k D .Thus, one can estimate the value of P ∞ cc if statistical information about the rock is not known a priori.In fact, one can use this relation to estimate the number density of the preexisting fractures if the orientation distribution is known.
Finally, from the scaling of the permeability of a stimulated reservoir, derived earlier in this section, its dependence on the injected volume is given by K where w is the average fracture aperture.Similarly, the porosity of the stimulated reservoir, φS, is given by φS = k i k s (V f k D /V inj ) β/(νD f ) .The proportionality constant k i is equal to k f at initial times when R ξ 0 and it is equal to k H when the homogeneous regime, R ξ 0 , is developed.The values of k k and k s are independent of the number density of preexsting fractures but are different in the fractal and homogeneous regimes.At scales smaller than ξ 0 , the preexisting fractures form a fractal network and the permeability and porosity of stimulated reservoirs are independent of the number density of fractures.This means that the performance of the fracturing process for fractal preexisting fractures is only a function of the total injected volume of the fracturing fluid and the fractal dimension of the network.

V. CONCLUSION
In conclusion, we showed that the sparsity typical of networks of activated fractures overcomes stress anisotropy and leads to the formation of a fractal network of activated fractures with D f = 91/48.The insight that the fracturing process is in the universality class of random percolation provides a means to relate the performance of the process to the statistical properties of the rock.The dependence of the percolation threshold value on the number density of the preexisting fractures was derived.We showed how this relation can be used to estimate the minimum injection pressure required to successfully stimulate a rock.In general, rocks with higher connectivity of the preexisting fractures require less energy to stimulate.A simple criterion to estimate the required injected volume of fracturing fluid required to fracture a hydrocarbonrich region of known size was provided.We showed how the transport properties of the stimulated rocks evolve with variations in the volume of injected fluid and the connectivity of the preexisting fractures.For fractal preexisting fractures, the permeability and porosity of the stimulated rock are functions of the total volume of fluid injected.If the viscous pressure drop becomes important at large length scales, the derived scaling of the permeability and porosity, which depend on the fluid pressure, can be applied locally at length scales above which the viscous pressure drop becomes important.Consideration of the effects of the broad length distribution often observed in natural fractures [29] on the fractal dimension of the network fractures is an important direction for future research.

APPENDIX: ESTIMATION OF THE PROPORTIONALITY CONSTANTS
In this Appendix, we explain how we estimated the value of the proportionality constants, k ξ , k f , and k H .
The values of k f and k H are calculated by fitting the finitesize scaling of F c , given by ( 11) and (10), in the fractal and homogeneous regimes, respectively.Using an approximate value of k −1/ν ξ , the correlation length of the preexisting fractures was estimated to be equal to ξ 0 = [0.4(n− n c )] −ν .Data that satisfies the condition that ξ 0 is more than ten times L and L is larger than 110 were used to fit (11).Figure 9 shows the values of F c for n − n c n c used to estimate the value of k f .Figure 10 shows the calculated values of k f for the different densities.Since the constant k f is independent of n, an average value was used, i.e., k f = 0.571.
Similarly, k H was calculated for different densities where the data collected must satisfy ξ 0 10L and L 200 and an average value was used.Figure 11 shows the values of F c used FIG.9. Values of F c used to calculate the value of k f .These data were chosen such that (11) is satisfied, i.e., ξ 0 L. ξ 0 ≈ 3000 for n = 5.643 was a lower boundary for ξ 0 .The system sizes are in the range 110 L 800.Since k f is independent of n and F ∞ c = 1 in the fractal regime, all the data almost collapse.
to calculate k H and Fig. 12 shows the obtained values k H for different number densities.
The value of k ξ was calculated using two different methods.In the first method, the data in the region where L(n − n c ) ν n 1 in Fig. 6 were fitted with (13) to extract the value of k −1/ν n ξ from the slope of the dashed-dotted line that is given by (k f − k H )k −1/ν n ξ .By now, we know the values of k f and k H . Thus, the value of k ξ was found to be equal to 0.4232.Figure 13 shows the data used for this purpose along with the fitting line.
Since the values of F ∞ homogeneous regime, the comparison between the calculated values of F ∞ c and the predicted values using (12) can only prove the linear dependence but not the form of the proportionality constant.Thus, one needs to use a different set of data in order to calculate k ξ .We constructed an empirical function to fit all the data in Fig. 6 including the transition regime.The function was constructed such that it can fit the data in the transition regime and has the correct behavior in the fractal and homogeneous regimes.In the fractal regime, i.e., (n − n c )L 1/ν n 1, finite-size scaling shows that (1 − F c )L 1/ν n = k f while in the homogeneous regime, i.e., (n − n c )L 1/ν n 1, finite-size scaling gives ( where k f = 0.571 and k H = 0.1925.Therefore, k ξ = 0.4232.One should note that the intercept is equal to k H which was found to be equal to 0.193.two limits: (A1) By taking the limit of the above function as L → ∞, one can recover the dependence of F ∞ c on the number density n, i.e., (12).The values of k f and k H were obtained from previous calculations, i.e., k f = 0.571 and k H = 0.1925 and k ξ was used as a fitting parameter.The empirical function with k ξ = 0.406 was found to give a good fit for the simulated data in Fig. 14.This value is in good agreement with that, k ξ = 0.4232, found earlier.
FIG. 14.The solid line given by (A1) when fitting to determine k ξ .
FIG. 2. A comparison between normal and β distributions to fitthe simulated onset of percolation.In this example, L = 20 and n = 5.639.

c
FIG.10.Extracted values for k f from data in the regime where ξ 0 L. Using these values, k f is estimated to be k f = 0.571 ± 0.0002.

1 −
FIG.12.Extracted values for k H from data in the regime where l ξ 0 L where l is the length of the preexisting fractures.Using these values, k H is estimated to be k H = 0.1925 ± 0.01.