Predicting hydraulic tensile fracture spacing in strata-bound systems

A model is presented which predicts the spacing of tensile-fractures due to ﬂ uid pressure increase in a multilayered sedimentary sequence comprising different typical sedimentary deposits such as mudstones, siltstones and sandstones. During normal burial and tectonic conditions, strata will undergo both extensional forces and an increase in ﬂ uid pressures. This model addresses the effects of the diffuse ﬂ uid pressure increase, and is useful for engineered applications such as the injection of ﬂ uid into a reservoir that may cause an increase of ﬂ uid pressure beneath a caprock, and for sedimentary sequences during normal digenetic processes of burial and fault activation. Analytical and numerical elastic stress strain solutions are compared to provide a robust normalised standard relationship for predicting the spacing of fractures. Key parameters are the local minimum horizontal stress, variability of the tensile strengths of the layers of a sedimentary sequence and the thickness of the beds. Permeability and storage are also shown to affect the fracture spacing. The model predicts many of the ﬁ eld observations made regarding strata-bound fracture systems, and should also prove useful in consideration of the impact of raised reservoir ﬂ uid pressures on caprock integrity. & 2013 The Authors. Published by Elsevier Ltd. All rights reserved.


Introduction
In the evaluation of the integrity of caprocks, and of analogue seals, the fracture spacing is of vital importance. In proposed CO2 storage sites, it is not the intact matrix of the caprock that causes concern for the retention of the injected CO2 rich fluids, or pure dense phase CO2. Rather, it is the presence of fractures at a series of scales which need to be quantified and analysed in terms of their connectivity and transport properties. During the characterisation of a reservoir for storage, the fluid pressure history and digenetic analysis of the caprock play an important role in understanding how the caprock will react to the presence of increased chemically aggressive fluid pressure loading beneath it. Indeed, the results of Rutqvist et al. [1] illustrate that hydraulic fracturing can be expected in the lower layers of a caprock after a relatively short period of time of fluid injection if the pressure is not controlled properly. It is generally accepted that hydro-fracturing will occur when the pore fluid pressure below the top seal equals or exceeds the minimum horizontal stress plus the tensile strength of the caprock [2].
Here we present a model that examines the impact of increased fluid pressure in multilayered sedimentary systems, the physical requirements for fluid driven fracturing of the strata in these layered systems, and the horizontal spacing between vertical fractures these systems could show. The model emphasises the importance of the local stress distributions on the formation of the fractures. It can be used to predict the likely fracture patterns of fluid driven (hydro-fracturing) in strata bound systems. A caveat to the model is that the presence of pre-existing fracture sets will influence the distribution of fluid pressure and impact on the predicted spacings. However, the model can be used for a firstorder assessment.

Controls on fracture geometry
Several authors discuss joint formation mechanisms. Here we concentrate on opening mode fractures. The work of Price [3] discusses joint/fracture development wherever the effective tensile stress exceeds the tensile strength of the rock. Possible causes being a result of fluid overpressure, expansion of the rock mass under uplift and erosion, pull apart due to tension induced by a regional extension, salt diapirism and folding.
There are obviously several mechanisms that will lead to the formation of fractures. The dominant mechanism at any particular time, and the characteristics of the deposit (the packet of sediment and hard rock, including any existing fracturing) will influence the nature of the response of the deposit to the formation of fractures. Depending on the cause of fracturing, the fractures formed will exhibit different geometrical characteristics and be scale dependent.
Bonnet et al. [4] review several methods of scaling fracture systems, including the lognormal distributions, exponential distributions and gamma law distributions, and indicated a recent preference for the fractal approach. They point out that recent studies indicate lithological layering, from the scale of a single bed to the whole crust, is reflected in fracture system properties. This layering influences the scale range over which individual bed specific or fracture system specific scaling laws are valid. The above named distributions are mathematical fits of probability distributions, and to understand the cause of fracturing it is necessary to reference the mechanical constraints and drivers. In certain cases one model, with certain limiting factors, fits better than another, but there is no ubiquitous law to match the whole population of fractures.
In a typical geological deposit, several sets of fractures will be present. To understand the spacing of the fractures it is important to understand the mechanisms that have led to the development of the different fracture sets. The observation that lithological layering is reflected in the fracture systems suggests that a process operating at the scale of the lithological bed size is important in controlling the development of the fractures. Identifying the key processes behind fracturing as creating "separate fracture packets" or end members, will help in the analysis of the fracture spacing and the nature of the process leading to the fracturing.
Here we concentrate on strata bound fractures as opposed to fractures which cut across several formations, and hypothesise that hydraulic fracturing provides an important controlling mechanism for the development of strata bound fractures. Particularly the stress field developed during dynamic fluid fracturing significantly influences the location of the development of further fractures.
Bai et al. [5] summarised work from many authors to make the observation that "the fracture spacing in layered sedimentary rocks is roughly proportional to the thickness of the fractured layer, with a ratio of thickness's from less than 0.1 to greater than 10." They developed a finite element model describing fracture spacing as a result of a pull-apart model, and a transition of stress from one bed to another bed. From the results of this model, they subdivided the fracture spacing to bed thickness ratios into four categories, whereby they could explain two categories with their extensional model and the further two categories where the joint spacing was too tight to have been caused by the extensional mechanism explained. They concluded that the other sets of joint spacing ratios required flaws and fluid pressure to produce the spacing. They note that as the tensile stress increases between two existing fractures, eventually a fracture will be initiated. The location of this new fracture will be dependent on a result of a local heterogeneity, such as a pre-existing zone of weakness, or due to the increase in fluid pressure overcoming compressive strength. Bai et al. [6] note that experimental and field results indicate fracture spacing decreases approximately as the inverse of the applied strain in the direction perpendicular to fractures, by fractures forming between earlier formed fractures. Gross [7] used the term "sequential infilling" to describe this process. Bai et al. [5] developed the concept of a maximum fracture saturation distance, being related to the stress distribution caused by the presence of a fracture leading to an area of "stress shielding" and thereby setting a lower limit to possible fracture spacing. The stress shielding is caused by the compressive stress due to vertical shortening of the fractures and the horizontal constraint in the central area between two fractures.
Addressing multiple layer sequences, Schöpfer et al. [8] examined the role of the transfer of extensional stress between different layers, and focussed on the relationship between the tensile strength of an individual bed and the amount of stress that can be transmitted into that bed from adjoining beds as a function of the interface shear strength. The larger the tensile strength, the more tensional stress needs to be transmitted to cause fracturing which is satisfied either through wider fracture spacing or through a higher interface shear strength. Following [8], and the references therein, this is described as Price's model [3]. They show that different extensional models are applicable depending on the ratio of the tensile strength of the bed to the interfacial shear strength; however, the influence of fluid pressure is not addressed.
Boutt et al. [9] investigated both experimentally and numerically the formation of natural hydraulic fractures. By reducing the external fluid pressure in sandstone samples more rapidly than the internal pressure could equilibrate, they were able to generate hydraulic fractures considered to be a consequence of both the internal fluid pressure exceeding the confining pressure and tensile strength of the rock, and also to be a consequence of a strong pressure gradient existing within the sample. Numerically they were able to simulate this type of depressurisation of the sample and the density of the in filled fractures. The rate of pressure release within the samples is a function of the permeability and storage of the samples. They conclude that the processes they have observed are very important in the natural hydro fracture process found within the earth's crust.
Odling et al. [10] examined several high-quality data sets of fracture systems from four reservoirs and identified two end member types of fracturing, named as "strata bound" and "non strata bound". They suggest that in strata bound systems there is little mechanical coupling between the layers. The individual joints are confined to single layers, and there is a clear relationship between bed thickness and joint spacing. Such sequences are found in systems with strongly developed interbedded weak and strong layers, e.g. interbedded sandstones, limestones, mudstones and shales. They describe the system as having weak adhesion between the layers. Odling et al. [10] describe strata bound fracture systems as confined to single layers, the sizes are scale restricted and the spacing is regular. We note also from the observation of typical caprock analogues (unpublished in house) that fracturing may at times go slightly beyond the limits of the bed and into more plastic layers, and also that fractures extending only a partial distance in the fractured bed (half fractures) are also present.
The role of increased pore fluid pressure within the crust and the link to the development of natural hydraulic fracturing has long been accepted, (e.g. [11,12]), and there are several examples in the literature of natural fracture systems which are interpreted as being a consequence of hydraulic fracturing. The focus of this paper is on strata bound systems and the role of fluid overpressure, and we suggest that it plays a more significant role than previously acknowledged in the formation of strata bound systems.

Parametric controls on hydraulic fracturing
There is a large body of literature particularly from the hydrocarbon industry examining the parametric controls on the development of hydraulic fractures in layered sequences. They deal particularly with a localised increase in fluid pressure due to fluid injection in a borehole, as opposed to a more regional increase in fluid pressure as would be the case in burial or a generic build up of pressure under a caprock. The key area of interest of this literature is the prediction of the length of the fractures generated and containment within different layers. There is some discussion as to the transfer of stress between different geological layers, key parameters being addressed include the contrast of the elastic modulus and Poisson's ratio between beds. Simonson [13] showed analytically that the effect of different elastic properties should have a very significant control on the vertical development of fractures. However, several authors have shown experimentally that this is not the case in nature. Smith et al. [14] discussed the impact of the modulus in layered sequences and showed that the contrasting modulus did not control growth of the fractures, but that the modulus was very significant in terms of determining the fracture aperture. Fung et al. [15] modelled a layered sequence using both an analytical solution and a numerical model to show that fracture development can be satisfactorily predicted without recourse to the elastic modulus even in cases where the modulus contrasts were a factor of 5. Amongst other authors Van Eekelen [16] came to the conclusion that the minimum stress profile was the key factor in controlling the growth of hydraulic fractures. This is related to the burial profile and Poisson's ratio of individual layers [17]. A number of times, shales have been shown to arrest hydraulic fractures, and it is suggested that this is due to the fact that mudrocks have a higher Poissons ratio than more indurated rocks. This leads locally to a higher minimum principal stress and therefore a higher fluid pressure necessary to cause fracturing. Conversely lower permeability leads to fluid pressures being maintained longer in the beds, and also a lower permeability will reduce fluid pressure drainage and therefore increase the overall fluid pressure in the system. Several authors (e.g. [18,19]) agree that in situ stress contrast is the dominant parameter controlling fracture growth. Warpinski et al. [19] suggest that material property interfaces are shown to have little effect after the examination of a number of large-scale hydraulic fracturing field tests and subsequent excavation.
The impact of material property interfaces is shown to be related to the amount of effective normal stress across that interface. This means that at deeper burial depths the system is more likely to be acting as a ridged body than at shallower depths, but still allows for slippage should enough fluid pressure be introduced into the interface zone. Zhang et al. [20] investigate hydraulic fracture propagation across frictional interfaces, and have developed a numerical model allowing the development of off-set fractures.
A synopsis of the above is that there appears to be a number of competing factors that influence the formation of fracture systems. For hydraulic fractures the key parameters which have been identified in order of importance are firstly the direction and magnitude of the minimum stress within the bed and secondly the tensile strength of the rock. Other important parameters include parameters relating to fluid flow, e.g. permeability and storage coefficient, and parameters relating to the elastic coefficients.

A simplified hydraulic fracturing model
Here we suggest that the extensional model provided by [5], amongst other authors, and developed by [8], presents one end member of possible mechanisms leading to bed thickness related fractures, and that a fluid pressure driven model can also provide another end member for these opening mode fractures with similar characteristics.
An advantage of a fluid pressure model is that it does not require the complex issues of extensional stress transfer from weaker beds to stronger beds in heterogeneous faulted systems, whilst still providing a method whereby fracture geometry consistent with field observations (e.g. [6,10] and references therein) is predicted.
Here we present a simplified hydraulic fracturing model based on the elastic interaction of fractures within the same bed. The model does not attempt to reproduce fracture tip stresses and predict the development of fracture pathways at a grain scale, rather at a larger scale it uses a simple pressure criteria to allow the development of a fracture within a bed as a function of the tensile strength of the material and the minimum principle stress. We assume that this fracture is fluid filled and do not at this stage try to represent transient flow conditions. Likewise at this stage the possible plastic and creep deformation to accommodate stress is not included, and is considered of secondary importance to the short-term stress distribution assumed in this model. This is because elastic response is not time dependent, plastic response is time dependent. The fractures are more or less confined to the strata beds as seen in strata bound systems; however, we assume that there is significant transfer of stress within the beds due to the pressurising of the fractures. In this way we integrate the main characteristics noted in the literature as controlling fracture development. To calculate the transfer of stress within the layered sedimentary system we apply two separate methods, an analytical approach and a numerical model and compare their results.
Although the model is a simplified representation of reality, it predicts and is validated by the main characteristics of strata bound fractures discussed in the extensive body of field observations in the literature cited above (e.g. [6,10] and references therein). The model also predicts, under the right loading conditions, the formation of synchronous orthogonal fracture systems. The formation of orthogonal fractures makes the model suggested here distinct from a pull-apart model. By normalising the relationships found in the model we are able to suggest a standard relationship to be expected from hydraulic influenced fractures dependent on the bed thickness and the tensile strength of the rock.
It is reasonable to assume that during the normal digenetic development of a sediment, flexing and burial, both pull apart and hydraulic fracturing models will be applicable, and reinforcing one another.

Hydraulic fracturing in strata-bound systems
Fluid pressure build up will occur naturally during the development of a multilayered sedimentary deposit as a result of burial and compression, a fluid charge from a deeper source or sudden settlement events such as ongoing tectonic activity. For burial to cause a sustained fluid overpressure, fluid migration in the layers needs to be restrained due to lower permeability layers. A stack of sediments will typically comprise sandstones, mudstones and siltstones. The model we present shows that the local minimum horizontal stress in the beds, the difference in the tensile strength of the beds, the difference in the permeability of the beds and the difference in the thicknesses of the beds are controlling factors behind the development of the strata bound fracture systems, the parameters thereof influencing the density of the fractures.
We start by postulating a simplified sedimentary sequence as a cut out from a typical multilayered sequence (Fig. 1). The sequence is saturated, and from base to top, there is a permeable sandstone or carbonate rock (maybe a reservoir rock), above this are two less permeable layers, whereby the tensile strength of the lower bed is less than the tensile strength of the upper bed, for instance the lower bed may be a mudstone, the upper bed a siltstone. The lower tensile stress bed has fractured normal to the minimum principal stress direction, s h . Possible mechanisms for this are discussed below, but for now what is important is consideration of the impact a fracture in this layer will have on the development of the local stress field in response to further hydrostatic pressure increase.
If higher pressure fluid is injected into the newly formed fractures, the geometry of the fractures will cause them to exert the hydrostatic pressure within them normal to the fracture walls.
The key behind the stress influence of the fractures on each other is the ability of the fracture wall to act as a load bearing surface in relation to the influx of extra fluid and increase in pressure within the fracture. The lower the permeability of the fracture walls (being in a low permeability deposit) the higher the load will be that is sustained and applied throughout the matrix, and the more localised the impact of the higher fluid pressure. In contrast in a higher permeability matrix the extra fluid pressure in the fracture will quickly be transferred to the matrix and be seen as a pore pressure increase.
The increase in fluid pressure in the fracture exerts a directional compressive stress on the bed normal to the fracture wall. This causes the development of a stress field represented by the Boussinesq bulbs of pressure [21] sketched in Fig. 1. As the fluid pressure continues to rise eventually a new fracture will propagate in the area of least horizontal stress development and depending on heterogeneities present, the most likely location being halfway between the existing fractures. This is because the "bulbs of pressure" dissipate the loading over an ever-increasing volume with distance from the extra (fluid pressure on fracture wall) loading.
At the location where a new fracture will be formed, the pore pressure is being driven from the shortest drainage path, either from the bed itself due to compaction or from the fluid source coming through the layer either above or below it. For the existing (vertical) fracture, the increase in fluid pressure is working against the fracture walls to increase the amount of horizontal stress in the bed. The extra compression due to the fracture loading works against the increased expansive force of the fluid pressure at the location where the new fracture is to be formed (and all other locations in the bed). The two forces are not equal however, due to the distribution of stress from the fracture wall, eventually the expansive force at the location of the new fracture overrides the minimum horizontal tectonic stress, the compressive force of the existing fractures and the tensile strength of the bed. This causes the formation of a new fracture normal to the minimum horizontal stress. As noted in the literature, e.g. [16], the key control on the location of the formation of the fracture is the local variation in the minimum principal stress.

Conditions causing hydro-fracturing
The initial fluid pressure u i (Pa) necessary to first open a fracture of the individual layers of tensile strength s t (Pa) in a confining stress field of s h (Pa) can be given as a first approximation as that is to cause a tensile fracture to develop both the confining stress and the tensile strength of the rock need to be overcome. If the fluid pressure exceeds this value then a tensile fracture must develop. The effective stress s e (Pa) is given by The detailed small-scale mechanism by which the fracture propagates through the concentration of stress at the fracture tip, and the exploitation of various heterogeneities and other weaknesses in the rock is outside the scope of this paper. What we consider is the fact that at a larger scale, a fracture will develop generally normal to the minimum principle stress once the tensile strength of the rock and the minimum principle stress has been exceeded.
In a draining medium the amount of effective stress is a measure of the amount of drainage occurred. In the case where u is increasing the effective stress becomes tensile, and failure occurs where it exceeds the absolute value js t j. The permeability of a bed influences the rate of the change in effective stress proportional to the drainage path length. In simplistic terms the fluid is either trying to get out of the bed (drain) or make space for itself.
Following Terzaghi [22], if the fluid pressure is caused by compaction of the fracturing bed then the higher fluid pressures are likely to be developed in the centre of the layer, as this has the longest drainage path to the higher permeability zones. Should the fluid be a charge, assumed from underneath, then we envisage the highest fluid pressure at the start of the fracturing located at the contact of the fracturing rock with the reservoir rock. We note that during compaction the fracturing layer could be above or below the reservoir layer.
Once the fracturing is initiated, it propagates normal to the least principal stress, s h . Fluid migrates into the fracture developing in the caprock until the fracture reaches the overlying layer of higher tensile strength and possibly different local s h . At this point the fluid in the fracture is at a higher pressure than the fluid in the surrounding matrix of the fracturing layer and at a higher pressure than the fluid in the confining layer. Drainage of this pressure occurs both through the overlying layer and into the fracturing layer depending on the relative permeabilities of these layers to each other, and the rate of recharge of the fracture fluid. Should there be a rapid rate of recharge then a higher pressure in the fracture can be expected and vice versa.
For the model, what is important is that the confining layer (Caprock facies II) retains a higher pressure than is required to fracture the fracturing layer (Caprock facies I). The confining layer does not fail until even more fluid pressure is applied.
Once a fracture has been developed in the fracturing layer, this fracture exerts the fluid pressure normal to the least principal stress. For the case where there is a fluid charge from beneath, vertically there is no differential stress seen in the fracturing layer as we assume the plan view extent to the layer is significantly more than the thickness of the layer. If the sequence is mechanically restrained vertically then the increase in the horizontal stress as a result of the increase in the fluid pressure s hΔu (Pa) is given by where υ is Poissons ratio, with Δu being the increase in the fluid pressure in the reservoir layer which we assume pushes the fracturing layer, but does not enter the fracturing layers matrix, i.e. the layer is compressed. If the strata are not mechanically restrained vertically, then some uplift will occur without a significant increase in the vertical stress. Including the effect of this uplift on u i , and allowing ϖ ¼ υð1 þ υÞ=1Àυ 2 we can write That is we assume that the confining stress s h comprises all rock mechanical contributions. To cause the rock to fracture we need to exceed s h by an increase in fluid pressure equivalent to s t . As per Eq. (3) this increase in fluid pressure will express itself in an increase in the confining pressure by ϖs t assuming vertical restraint, therefore the term u i needs to accommodate this. In the case where there is no vertical restraint the term ϖs t is not applied as we assume simple uplift. We introduce the term r as a value between 0 and 1 to indicate the degree to which vertical restraint is applied. For a completely restrained system r ¼ 1, for a totally free moving system r ¼ 0. For ease later we define the pressure above the minimum principal stress for fracture formation u f f to be At the point of fracturing we have higher fluid pressure in the fracture than in the fracturing layer. This will equilibrate with time, the rate being dependent on the permeability of the bed, and if we allow some mechanical pore deformation then it is also inversely proportional to the storage of the bed, the combined effect being the pressure diffusivity, e.g. [9,23]. At the edges of the area if there is little constraint on the layer, plastic and elastic strain accommodation of stress will occur. Outside of the boundary region this release is not available; therefore, there is an increase in the least principal stress experienced by the fracturing rock.
The increase in the least principal stress caused by the loading of the fracture walls by the fluid pressure is fundamental to the location of the development of the fractures. We compare both an analytical solution and a numerical solution to evaluate the distribution of stress under this fluid loading (bulbs in Fig. 1).

The effect of vertical fractures on the horizontal stress distribution
Integrating the Boussinesq [21] equation for a point load, e.g. [24], it is possible to obtain a number of elastic solutions for different geometrical conditions. For an analytical solution, referring to Fig. 2, we approximate the fracture which has been developed in the fracturing layer as an infinite strip foundation with a width of 2B (Fig. 3) exerting pressure normal to the strip, in the fracture case normal to the fracture wall. The standard solution and geometrical arrangement for the elastic solution of this stress field development for a semi-infinite layer are presented in Fig. 3 [25].
When applying the Boussinesq approach, the elastic modulus is taken as not having a significant impact on the stress distribution. As discussed in the Introduction, several authors show that contrasting moduli introduce minor quantitative differences in fracture spacing. In applying this strip solution we are interested only in the stress distribution, and assume full transfer of stress between the individual layers. Also we assume that there is a smooth frictionless contact between the fluid pressure in the fracture and the matrix, and that the side of the walls of the fracture is flexible and detached from one another.
The semi-infinite layer assumption suggests that the stress seen at the fracture will be seen in some way throughout the whole of the fracturing layer. As a first assumption this method is useful in understanding the distribution of stress in the fracturing layer, and the principle of stress superposition can be applied for subsequent fractures. However, to take into account that fractures in the fracturing layer will be developing parallel to each other and significantly influencing each other the closer they are together, it is necessary to select an analytical solution that encompasses this.
The increase in horizontal pressure as a result of the increase in vertical stress via Poisson's ratio is considered ubiquitous as the horizontal area of the bed is assumed significantly more than the thickness of the bed. Therefore, this stress increase will not be dissipated. However, to calculate the dissipation of the increase in horizontal stress caused by the loading at the fracture walls we apply a coupled processes finite element simulator, "OpenGeoSys" [26] and also use an analytical solution developed by Poulos [27] for a foundation underlain by an adhesive rigid base. The numerical model solved by OpenGeoSys is a standard elastic solution [28]. Although the approximation of the adhesive base in the analytical solution is incorrect, it is interesting to note the similarity between the numerical solution and the analytical solution. This indicates that the key processes have been addressed.
In foundation engineering numerous authors have used the symmetry of the elastic solutions to produce dimensionless relationships between the size of the bearing surfaces and various quantities located in the elastic half space such as directionally orientated stress or strain. Often given an initial loading of the bearing surface the quantity required is given by multiplying the initial loading value by an "Influence Factor" I. Poulos et al. [25] contains several examples of this approach, behind the influence factor is the spatial integration of the Boussinesq [21] equation for a point load, which then enables the ratio of the initial loading to  the value required at a point located in the elastic half space to be calculated. For analytical solutions, as the stresses and strains can be superimposed various standard solutions for rectangles, or strip loading, and points under corners or under the middle can be combined to provide ready access to required variables, e.g. [29].
Here we use the approach of using a numerical model that satisfies all the mechanical balance equations for the elastic solution, and an approximate analytical solution to calculate the influence factor of loading on the fractures, and the impact this will have on the minimum principal stress.
The analytical solution used is for the stress increase with depth (z) under the corner of a strip foundation upon a finite layer of thickness (h) underlain by a rigid base [27]. The rigid base is taken as the point of meeting of the influence of two fractures on each other, acting as a rigid base (see Fig. 2). Strictly this is incorrect as the stress fields will superimpose on each other and there will be stress and strain developed normal to the minimum principal stress. This is accounted for in the numerical model used.
We approximate the influence factor for the analytical solution I st shown in Fig. 4 as a polynomial function of B/h, and then estimate the increase in stress s Δu f as a result of the additional  fluid pressure loading ðΔuÞ on the fracture following [27] as: The numerical modelling approach is illustrated in Fig. 5. As discussed above the results can be presented in a dimensionless form. A high-density mesh (circa 32,000 elements) was used to represent the geometric space between two fractures. The boundaries of the mesh were far enough away from the interacting fractures so as to assume that the boundary conditions were negligible. The fractures were loaded as flexible sections with fluid pressure and the finite element scheme used to solve the elastic equations within the mesh. The edges of the mesh were set such that no movement in the y direction was possible, only the fracture faces were allowed to deform, and the fractures were not allowed to extend into the surrounding beds. The stress developed in the y direction in this model then represented the increase in horizontal stress. The stress at points P1 and P2 (Fig. 5a) was then used to evaluate the behaviour of the system and derive a function for the influence factor both in the middle of the bed and at the edge of the bed. For the centre of the bed substituting B/h with x the influence factor I for the numerical solution was found to be I ¼ À0:0313x 5 þ 0:3039x 4 À1:0394x 3 þ 1:3443x 2 þ 0:0393x þ 0:0015 for 0 ≤x ≤2 ð7Þ In this case and following we evaluated the influence factor such that Both the imperfect analytical models and the numerical model predict that the stress seen at the edge of the bed will be less than the stress increase seen in the centre of the bed, an important fact we will readdress later.
Using the influence factor we find that the most likely position for another fracture to develop will be halfway between two existing fractures, as per sequential infilling observed in the field. The influence existing fractures have upon the minimum principal stress between them can be calculated, and therefore the fluid pressure required to cause hydro-fracturing and the generation of a new fracture at this location can be evaluated.

Dynamic system
As fluid pressure is building up at the base of the layer undergoing hydro-fracturing, so also the fluid pressure is increasing in the fractures within this layer which have already formed. This leads to a dynamic system illustrated in Fig. 6.
Illustrating this, let us assume that we require an overpressure of 1 MPa (we call this value u f s ), to cause the first tensile fracturing uninfluenced by surrounding fractures (illustrated in Fig. 6). This pressure we call u f s (where s (set) ¼1 being the first set of fractures to form) is the fluid pressure required to cause the first set of fractures, initiated due to the heterogeneities in the rock at random weak locations in addition to the in situ minimum principal stress s h . If we assume that the minimum principal stress is 10 MPa, from (4), assuming that we allow uplift, i.e. no vertical restraint so that ϖ ¼ 0, we know that u i ¼ 11 MPa. Therefore, For the first set u f f ¼ u f s . Should the pressure now remain at 11 MPa, we can predict the extra pressure now required to cause the next set of fracturing u f ðsþ1Þ , being 1 MPa (u f s ), plus the extra horizontal stress across the location where the next fracture set is to be formed which is evaluated as Iu f s (arrow 1 in Fig. 6), plus s h . Let us speculate that this value is 1.1 MPa ðu f s þ Iu f s Þ plus s h . Should the fluid pressure now increase to 1.1 MPa (arrow 2 in Fig. 6) plus s h , we need to take account of the fact that the 0.1 MPa increase in pressure in the fractures will also have a further compressive effect across the new fracture location, arrow 3 in Fig. 6. This can also be evaluated as in Eq. (8), however, we note that we now have a dynamic system, where the local minimum principal stress across the new fracture location is increasing with the fluid pressure increase in the existing fractures. The pressure u f ðsþ1Þ in addition to s h required for a new fracture f ðs þ 1Þ to form can be expressed as a power series which is a convergent geometrical series as long as I ≤1, expressed as The sum of a convergent geometrical series is given by The actual real pressure p for the formation of the fracture set s +1 is given by Once the pressure has been reached for the next set of fractures to infill, u f ðsþ1Þ the pressure for the following set of infilling fractures can be evaluated as s¼ 2.
To include the effect of the vertical stress developed if the strata sequence is mechanically restrained vertically, the horizontally induced component of the vertical stress is calculated as in Eq. (3). We set the first tensile stress fracture of the layer at a defined value as per Eq. (1), and assume that s h in the layer at this moment contains all the resolved stress components. Further increases in the fluid pressure in the underlying rock now act also to uplift and further compress the fracturing rock, leading to an increase in the horizontal compression. This is included in Eq. (10) as follows: Again using the sum of a convergent geometrical series this can be expressed as For graphical representation this equation can be normalised against u f f , Eq. (5), and values of I calculated as a function of B=h, such that We select a distance such that I≈0, e.g. 1024 Â bed thickness, and assume at this distance that for the first open fracture set u f s ¼ u f f . The actual value of I is calculated and the equation solved for the value u f sþ1 ð Þ , this being the pressure additional to s h required to cause a fracture halfway between two existing fractures to open, i.e. fracture set s+1. This procedure is now repeated whereby the calculated value of u f sþ1 ð Þ in the previous iteration is now the value used for u f s . Repeating this procedure for the interval for which the fitted function of I is valid allows the construction of Fig. 7.

Results and discussion
The evaluation of the impact of bed thickness and fracture spacing can be normalised against u f f and the ratio of bed thickness to distance to the next fracture to provide a standard relationship, (16). This is given in Fig. 7a for the numerical solution and Fig. 7b for the approximate analytical solution.
The curves presented in Fig. 7a can be used in a graphical fashion to determine pressure and spacing relationships. For a more quantitative approach these curves have been fitted such that for the y¼f(x) relationship And for the inverse x ¼f(y) relationship where x is the fracture spacing/bed thickness. Parameters and validity range are given in Table 1. The accuracy of these curves is presented in Fig. 8.
As an example, let us say we want to evaluate the fluid pressure required to create strata bound fractures at a spacing of 20 m in a bed with a thickness of 2 m. Let us postulate that the tensile strength of the bed is 5 MPa, and that there is an overlying bed with a higher tensile strength, and low enough permeability to cause the necessary pressure build up in the bed we are looking at.
The bed thickness to fracture spacing ratio is 10, therefore from Fig. 7a or Eq. (17) the fluid pressure to tensile strength ratio is circa 1.03. This means that the fluid pressure required is then the tensile strength of the bed 5 MPa multiplied by 1.03, giving 5.15 MPa plus the horizontal stress.
Comparing Fig. 7a and b shows that although the two solutions do not provide identical results, the main features of the behaviour have been captured. The numerical solution more closely satisfies the boundary conditions and the distribution of stress within the model and is taken as to be the preferred curve. It is apparent that the distribution of stress within the models is the key feature for the development of the fracturing location.  It is interesting also to note that the elastic modulus is not included in the development of the stress field due to loading in analytical solutions [21,27]. This is reflected in the numerical model. The impact of Poisson's ratio on the stress field during the uplift caused by the fluid pressure and local changes in the minimal principal stress is estimated as described above in Eq. (16) and presented in Fig. 9.
Given the limited divergence due to possible differences in Poissons ratio and also the similarity between the mathematically correct numerical model and the approximate analytical solution used, the relationship described seems to be fairly robust. The processes it illustrates can be applied to understand a number of phenomena.
The difference between the fluid pressure required to fracture the centre of the bed, and that required to fracture the edge of the bed, is minimal until the spacing is reduced to about four times the bed thickness. This suggests that fractures that do not fully transect the bed will develop in the later stages of fracturing at higher pressures.
As fluid pressure increases, we note that the horizontal stresses should be of a similar size, as may be expected during early burial, that it is possible to fracture the systems orthogonally. Differences in s h and s H would be reflected in the spacing of the fracture sets. Fig. 10 demonstrates the development of fracturing with a fracture spacing to bed ratio of down to circa 1.5. If we postulate the tensile strength of the bed fracturing is 1 MPa, then the layer causing the pressure build up in this case has a tensile strength of at least 2.5 MPa, illustrated in Fig. 10 at 3.6 MPa. As the fluid pressure increases the degree of heterogeneity in the fracturing layer will determine initially the location of the first sets of fractures (moving on Fig. 10 from right to left along the bottom). However, as soon as the difference caused by the heterogeneity is less than the stress superposition of the fracture systems, the general heterogeneity will play less of a determining role in the location of the fractures. Obviously the presence of already existing fractures and other significant planes of weakness may dominate the location of all the fracturing. From Fig. 7a, it can be seen that the influence of stress interference becomes more significant under a strata bound fracture spacing of circa 20 bed thicknesses.
Following again as the pressure increases so the next four sets of sequential infilling fractures arise until the fluid pressure for the next smallest set of fractures exceeds the tensile strength of the confining layer.
The spacing and therefore the number of fractures within a bed will be a function of the relative tensile strength of the bed in comparison to the other beds within the stratigraphic sequence. As soon as the fluid pressure has been able to rupture the confining layer, hydraulic fracturing in that bed will stop until a higher pressure can be retained.
This means for more general observations there will be no hard and fast rule for the spacing of fractures as a function of the tensile strength, rather in a sequence the lower tensile strength deposits will be more densely fractured, and in sequences with high amounts of tensile strength variability this will be reflected in the increased variability of the strata bound fracture spacing.
If for a given fluid pressure (p g ) an estimate of the expected fracture spacing for a bed is to be evaluated, first the horizontal stress confining stress component is removed Secondly the fluid pressure u f f is calculated for the bed (5) and the pressure u g is normalised (divided by) u f f . Third using Eq. (16) and the relationship illustrated in Fig. 7a, the bed thickness to fracture spacing ratio is determined by setting u f ðsþ1Þ =u f f ¼ u g =u f f and reading of the graph, or using Eq. (18). Fourth this is converted to a real number by multiplying by the bed thickness.
As discussed previously the key behind the stress influence of the fractures on each other is the ability of the fracture wall to act as a load bearing surface in relation to the influx of extra fluid into the fracture and the matrix when fluid is present within it to act under hydrostatic stress. If, due to the high permeability of the bed, the fracture wall is not able to act as a load bearing surface the process will be arrested, and therefore there will be less control on the location of new fractures. This suggests that there should be more variability seen in the spacing of strata bound fractures within higher permeability deposits than within lower permeability deposits once normalised against the bed thickness and the tensile strength of the rock.
Additionally the amount of fluid pressure generated in the sedimentary profile will also be a function of the rate of the source supply, the permeability of the individual beds and their thickness. Allowing normal Darcy flow, the amount of flow is a linear function of both the pressure gradient and the permeability. Therefore, if a source is defining how much flow there is to be through a system, this will define the pressure gradient in the system as a function of the contrasting permeabilities of the beds to each other. The pressure gradient across a bed is a linear function of its permeability and an inverse function of its thickness. Therefore, the thicker a bed and the lower its permeability the higher the pressure will be necessary to sustain a constant flow rate. In a source term driven system, the source term is forcing fluid through a sequence and if the rate of the source term increases then the pressure gradient has to increase to  accommodate this. This increase in pressure could be enough to trigger the hydraulic fracturing described above.

Conclusions
We identify the hydro-fracturing process as one possible mechanism for tensile fracture development and present a model for investigating the characteristics of tensile fracturing driven by fluid pressure increase in multilayered sedimentary systems. The model allows the derivation of a standard normalised rule applicable to all strata bound systems. We suggest that both extension and fluid pressure-fracturing can operate during a normal digenetic burial process. We also conclude that the fluid pressure model presented here will better explain certain features seen such as orthogonal fracturing. This model will also be applicable during larger scale engineered fluid injection into reservoirs under caprocks.
The model predicts that strata bound fracture systems will follow a standard curve during hydraulic fracturing which can be used to determine the pressures of fracturing as a function of the spacing of the factures and the tensile strengths of the beds. The key feature of the model is the interaction of fractures with one another, the buildup in fluid pressure and the tensile strength of the individual layer. The stress superposition caused by the fluid pressure loading in the fractures coupled with the heterogeneities in the beds defines the development of the spacing of the fractures. The model predicts that the following factors will be related to the spacing of the strata bound fractures, given below and summarised in Fig. 11.
First the bed thickness: this is directly related to the fracture spacing, as the thickness of the bed acts as the length of the load bearing surface which defines how far the effects of the fluid pressure increase at the fracture walls is transmitted into the bed fracturing.
Second the permeability: the higher the permeability of a bed, the more varied the possible fracture spacing, as with increasing permeability the fracture walls will act less efficiently as load bearing surfaces, spreading the fluid load distribution within the matrix and thereby increasing possible fracture spacing. Also the variability in the contrasting degree of permeability within the system will be related to the fluid pressure profile in the system. The larger the range in permeability the more there is the possibility of a low permeability layer causing fluid pressure to increase beneath it.
Third the bed thickness variability: the variability in the contrasting degree of bed thicknesses within the system, particularly the lower permeability beds will control the fluid pressure profile in the system. Thick low permeability beds will allow the buildup of higher fluid pressure them it.
Fourth the rate of source: the larger the rate of fluid charge, the less variability there will be in fracture spacing because the fluid pressure does not have as long to dissipate into the matrix before fracturing occurs. The larger the rate of the source the more efficiently the fracture walls will act as load bearing surfaces.
Fifth the tensile strength of the beds: the variability in the contrasting degree of tensile strength within the system (high tensile strength of beds will reduce fracture density, a high tensile strength bed of low permeability will cause the beds under it to hydro fracture) Sixth the size of the local principal horizontal stresses will define whether parallel or orthogonal fracturing will occur, and their relative spacing.
Additionally the model predicts that there is a minimum fracture spacing, however combined with an extensional regime fracture spacing can be reduced. Finally the model also provides an explanation for fractures which extend only partially through a bed, and suggests that they are formed at later stages and higher fluid pressures. Also in agreement with other work the model suggests that the elastic moduli (Young's modulus and Poisson's ratio) have little impact on the fracture spacing.