Effective stress history and the potential for seismicity associated with hydraulic fracturing of shale reservoirs

The influence of effective stress history on reservoir stability and the response to hydraulic fracturing are simulated for two hypothetical shale reservoirs. It is numerically demonstrated that the effective stress history influences the present-day stability of faults and natural fractures. Any assumption of a homogeneous stress state is shown to be unrealistic in the majority of fractured shale reservoirs. The simulations demonstrate that, depending upon the effective stress history, shear displacements can be induced by stress changes associated with stimulation without a change of fluid pressure in the destabilized fractures. Elastoplastic shale and other fractured reservoirs are shown to have a memory of past geomechanical states. These numerical findings demonstrate the value of interpreting reservoir effective stress history when planning fluid injection. A previous case history is summarized to illustrate the overwhelming dependence of stress history analyses on a wide spectrum of earth science disciplines. It also illustrates the application of stress history analyses to activities other than hydraulic fracturing. The current inexperience of effective stress history description means that new approaches must be identified to optimize the necessary multidisciplinary investigations. Supplementary material: Details of the simulations are available at www.geolsoc.org.uk/SUP18743.

The stability of faults can be affected by a number of human activities. A previous study of seismicity induced by water impoundment identified stress history as a key influence. The focus here is on the influence of stress history on seismicity induced by fluid injection, in particular for purposes of hydraulic fracture stimulation in shales. The need to apply a wide range of geological disciplines in the analysis of effective stress history is emphasized.
The potential for seismicity induced by hydraulic fracturing can be analysed by comparing the balance of disturbing forces and resisting forces acting on known or assumed fault planes, incorporating the calculated pore pressure increase resulting from injection (Cornet et al. 2006). A homogeneous stress state is normally assumed in comparisons of the disturbing and resisting forces. As will be illustrated here, this approach has a number of disadvantages related to the difficulty, or impracticality, of accurately describing the stress state. Successful hydraulic fracture completions in shale reservoirs depend upon the presence of pre-existing natural fractures, either open or filled (King 2010). The stress state may exhibit substantial variation in such fractured, heterogeneous hydrocarbon reservoirs (Ketter et al. 2006;Harper 2011) and will be a function of the history of slip on faults and natural fractures. Determination of the reservoir fluid pressure, necessary for determination of effective stress, is notoriously challenging in such lowpermeability rocks. The presence and geometry of faults may be ill-defined unless the quality of the 3D seismic data is of an exceptional standard. The allocated or assumed properties governing the sliding resistance of a fault are typically uncertain or, at least, imprecise.
Segall (1989) analytically examined stress changes in a linear poroelastic medium and showed that stresses are generated where there are no changes in pore fluid content. A series of papers by Heffer and other researchers concerning oilfield waterflood correlations (e.g. Heffer & Dowokpor 1990;Heffer et al. 1994;Zhang et al. 2007) prompted Maillot et al. (1999) to develop a numerical model of the initiation, growth and radiation of shear failures induced by pore pressure changes in a poroelastic medium. Otherwise, simulation of stress transmission is either confined to the seismogenic source structure and greatly simplified using a Burridge-Knopoff model (Baisch et al. 2009) or omitted (e.g. Hummel & Shapiro 2010) with few exceptions (Pine & Batchelor 1984;Harper & Last 1989;Last & Harper 1990). Some processes that might occur are not normally taken into account, such as small stress changes induced in the reservoir at locations remote from the pressure disturbance or a time-dependence brought about by low rates of fluid flow into dilating fractures in low-permeability rocks such as shales. Stress history is ignored (except in the sense of stress changes local to faults and attributable to slip; e.g . Cundall 1990;Heffer et al. 1994).
For these and other reasons, a 'traffic light system' (Bommer et al. 2006) based on monitoring microseismicity has been proposed and will be implemented in the UK (DECC 2013) during hydraulic fracture operations. Threshold event magnitudes are predefined and operations conducted accordingly, including termination of the fracturing operation if a threshold event magnitude is reached. The traffic light approach normally assumes that instability of the microseismic source fracture depends upon change of the effective normal stress brought about by pore pressure diffusion away from the injector. Ceasing injection limits the maximum pore pressure to which the potential source structure is subjected.
Typically, the maximum event magnitude is experienced some hours after injection ceases (Charléty et al. 2007;Eisner et al. 2011), explicable in terms of pore pressure diffusion away from the injector well (Shapiro et al. 2003;Parotidis et al. 2004;Rozhko 2010). The elapsed time required for observation after the end of injection and prior to the next phase of injection would normally be selected on the basis of calculations of pore pressure diffusion.
There are examples of unexpectedly large events occurring after some surprising lengthy periods of time after injection has ceased Effective stress history and the potential for seismicity associated with hydraulic fracturing of shale reservoirs (Charléty et al. 2007;Bommer et al. 2006). McKinnon (2006) examined hard rock mining-induced seismicity and noted that while the majority of events are close to the excavation, events have been recorded hundreds of metres from the excavation. These tend to be of larger magnitude. They appear to be uncorrelated in time with mining activities and result from slip on pre-existing discontinuities such as faults. These considerations jointly suggest that the traffic light approach should be augmented by additional analyses.
By far the majority of hydraulic fracturing operations (Warpinski et al. 2012) and fluid waste injection wells (Frohlich 2012) appear to be aseismic in terms of felt events. Only a small minority (Nicol et al. 2011;Frohlich 2012) have been identified as probable causes of felt induced seismicity. This observation invites an explanation.
The first objective here was to reintroduce the hypothesis of Snow (1972), which appears to have been largely ignored. Snow (1972) was concerned with seismicity induced by water impoundment behind dams. He hypothesised that reservoirs may be aseismic if the underlying rocks have previously experienced a minimum effective normal stress lower than that induced by the impoundment of water. In the current study Snow's hypothesis was tested numerically using a computational mechanics code to examine the stability of two fractured shale reservoirs subjected to different hypothetical effective stress histories before stimulation. Each reservoir was subjected to simulated multi-stage hydraulic fracturing to investigate the effect on the distribution of induced slips and the influence of effective stress history. No attempt was made in either Snow's analysis or in the current analysis to distinguish between aseismic, stable sliding and (potentially seismic) stick-slip or to estimate the magnitudes of the seismic events which might be generated as a result of slip events.

Terminology
In common with Suckale (2010), the term 'induced' is used here in relation to seismicity without implying that the reservoir stress state is predominantly a result of injection. Quite the contrary, the processes discussed here are almost entirely a consequence of the release of reservoir strain energy. (Any seismicity associated with the release of 'virgin' reservoir strain energy has been classified as 'triggered' in some other reports.) Productive shale reservoirs are naturally fractured (King 2010). The term 'limit equilibrium' is used here to denote that the ratio of disturbing forces to resisting forces is unity for a given discontinuity (normally pre-existing weakness planes such as natural fractures, faults and bedding planes). This follows the terminology used in soil and rock mechanics literature. The phrase 'marginal equilibrium' is used to refer to bodies of rock (shale reservoirs here) that include regions that are in a state of limit equilibrium (normally because of the presence of one or more discontinuities that are in limit equilibrium). Some researchers use the term 'critical' in place of both limit equilibrium and marginal equilibrium. This can lead to confusion with usage of the term 'critical' to denote phase transitions and percolation thresholds, as noted by Grasso & Sornette (1998).

The effective stress threshold hypothesis
The impoundment of water following dam construction has commonly induced seismicity (e.g. Kumar et al. 2012). Forty years ago, in an attempt to explain why some such construction projects induced seismicity and others did not, Snow (1972) introduced the concept of a threshold effective stress state. Snow hypothesized that the key to whether seismicity was induced was whether the modified effective stress state, reduced by water impoundment, reached the (geologically) historical effective stress minimum for the region. He argued that if suitably oriented potential source faults had been subject to a more unstable condition (effective stress state) in the past, then they would remain stable and seismicity would not be induced despite the increased water pressure.
In the broadest of terms, effective stress changes occur predominantly because of changes of pore pressure, changes of normal stress or changes of shear stress. Typically, each of these will influence the others. Normal stress changes can result simply from geological displacements, either horizontal or vertical. Shear stress changes can be imposed by tilting or bending, such as during postglacial rebound. The changes resulting from pore pressure rise and an imposed shear stress, each promoting failure, are illustrated in Figure 1. Snow (1972) conducted a plane strain analysis of the effective stress changes in a fractured medium subject to loading by a body of impounded water. For normal faulting, wrench faulting and thrust faulting environments, Snow analysed water impoundment above a potentially unstable fault and a rise of the water table distant from the surface impoundment. Most importantly, he suggested that earthquakes can be triggered only if stress falls below some previous threshold, 'developed progressively or at some time in the past'. He noted that any area with historical earthquakes or even with recent surface expressions of faulting is probably in a state of stress close to failure, and that 'in seismically active areas, stresses are apparently so close to the strength limits that local load changes of only a few bars may trigger brittle failure'. Snow stated that this hypothesis 'can only be tested by collecting case histories consistently demonstrating its qualitative correctness'. An effective stress history for the Kremasta Lake area in Greece was constructed as an example case history.
The identification of seismically active areas may be critically restricted by the duration of the available record of historical seismicity. For example, the Bowland Basin in the UK is not recognized as a seismically active area yet induced seismicity has been associated with hydraulic fracturing (Eisner et al. 2011).
The hypothesis of Snow (1972) implies that rocks have a memory of past events. This is consistent with the behaviour of dynamical systems (Mees 1981), a description of geological processes suggested previously in the context of stress state evolution (Harper & Szymanski 1991;Ribiero 2002).
Snow's threshold might be compared with the Kaiser effect observed in uniaxial test specimens of metals, which is an irreversible occurrence of acoustic emission (AE) once the previous peak stress has been exceeded. The Kaiser effect is commonly thought to relate to a damage process (cracking of intact material). It was discussed in the context of rocks by Sondergeld & Estey (1981). Those researchers noted that the Kaiser effect placed no restrictions Fig. 1. Examples of stress change (Mohr circles) leading to failure and a state of limit equilibrium on a fracture (with a linear failure envelope). Left, pore pressure increase; right, shear stress increase. τ, shear stress; σ n , effective normal stress. on the nature of the applied stress field and they failed to find an effect in rocks that is exactly comparable with the Kaiser effect. A systematic migration of seismicity away from the injection well in geothermal reservoirs has been described as a Kaiser effect (e.g. Baisch et al. 2009). Because shearing on pre-existing fractures does not necessarily entail significant rock damage, this terminology may be inappropriate. More likely is a short-term reduction of cohesive strength that may not endure for time scales characteristic of many geological processes.
Many researchers (Hast 1967;McGarr & Gay 1978;Harper & Szymanski 1981) have suggested that stress state is controlled by the frictional strength of crustal rocks, either the intact rock or planes of weakness. Some widely used analyses of stress state (Zoback et al. 2003) imply that we can identify a stress state that is homogeneous throughout the region of a given fault and its surroundings. yet, as will be shown here by examples, in pervasively fractured reservoirs such as productive shale reservoirs the stress state is likely to be far from homogeneous. These analyses (Zoback et al. 2003) involve identifying bounds on the possible stress states by assuming that the stress state is limited by slip on faults. This implicitly assumes that the stress state around planes of weakness does not change with the history of slip (unless the physical properties governing fault strength change) and that limit equilibrium is continually restored by geological driving forces. Such analyses should be used to analyse the potential for induced seismicity only if the proposition of Snow (1972) is unfounded and if the heterogeneity of the stress state is minor. Otherwise, the concept of a uniform stress state controlled by rock strength and insensitive to geological history is a convenient simplification that cannot be assumed to be precise and universal. Clearly, Snow recognized this. Cornet et al. (2006) examined the stress state and microseismic response to fluid injection in a granite body. These researchers' findings are supportive of the suggestion of Snow (1972) in that they concluded that the stress state at this location is not controlled simply by the frictional strength of the main faults.
It appears, therefore, that there are limitations to any analysis based on a comparison of solely the current stress state (even with allowance for the natural lithological heterogeneity) and fault properties, total reliance on the traffic light system (as currently founded on pore fluid pressure diffusion in a reservoir with a simplified stress characterization) or identification of seismically active regions as means to define the potential for induced seismicity. Consequently, if the hypothesis of Snow (1972) of the critical role of any prior threshold failure state could be corroborated, it could lead to an enhanced ability to predict and analyse induced seismicity. It might help to explain why so many regions do not experience significant (felt) seismicity during hydraulic fracturing yet a few do. Implementation of Snow's approach would require the acquisition of geological information and its interpretation in terms of the geomechanical evolution of the subject reservoir (partly illustrated by Harper & Pritchard 2006. Fortunately, suitable geological information is often at least partly available. For example, a majority of oil and gas operators acquire much of the required information as part of normal exploration practice before well completions are designed. Snow (1972) stated that the threshold hypothesis 'can only be tested by collecting case histories consistently demonstrating its qualitative correctness'. Regrettably, rather than testing Snow's proposal, the focus of most subsequent analyses of induced seismicity appears to have been concerned with present-day geomechanical conditions, other than perhaps the work by Cornet et al. (2006) discussed above. Kumar et al. (2012) reported some evidence that may be relevant but is incompletely supportive without an analysis of effective stress history. Those researchers concluded that the factor most conducive to high seismicity in the Koyna-Warna region of India is that impounded water levels exceed previous impounded water levels.

Testing the threshold hypothesis
The analyses here use computational mechanics to test the hypothesis in a different context: the fracture stimulation of shale reservoirs. Because hydraulic fractures are designed to create a mode 1 fracture by imposing a fluid pressure greater than the least principal stress, Snow's threshold will normally be attained in the local area of the treatment. Given that injection-induced seismicity occurs in only a small percentage of cases (Frohlich 2012) and that events associated with hydraulic fracturing of shale reservoirs typically are extremely small in magnitude, energy and size (Warpinski et al. 2012) it can be inferred that, if Snow's hypothesis is correct, a state of marginal equilibrium, or close proximity to such (a few bars), does not exist in most shale reservoirs exploited to date (assuming the presence of fractures potentially able to generate felt events).

Summary of the numerical tests
The analysis by Snow (1972) concerned changes of vertical load and pore pressure. In the context of fracture stimulation of hydrocarbon reservoirs, tests of the threshold hypothesis can be based upon the effect of pore pressure changes on the stability of natural fracture populations and a suitably oriented fault. Changes of pore pressure can result simply from burial and uplift processes during basin subsidence. In hydrocarbon reservoirs, pore pressure is commonly increased by oil and gas generation (see, for example, introductory discussion of Fan et al. 2012).
The effects of fracture stimulation treatments on the distribution of slips induced on natural fractures and a fault are compared for two hypothetical reservoirs. Different 'threshold effective stress states' are achieved simply by subjecting each reservoir to one or more transient phases of slightly elevated pore pressure prior to stimulation. The reservoir pressures were 35.29 MPa (equivalent to a gradient of 0.52 psi ft −1 ) in Reservoir A and 33.94 MPa (equivalent to a gradient of 0.5 psi ft −1 ) in Reservoir B. For each reservoir, multistage fracture stimulation was simulated both at a 'threshold effective stress' condition and in a more stable condition after evolution from two previous 'thresholds'. The two different 'thresholds' were achieved by transient increases of reservoir pressure by 0.68 MPa and by 1.36 MPa (equivalent to increments of gradient of 0.01 and 0.02 psi ft −1 ).
The stress distribution around mode 1 penny-shaped cracks in elastic media is well known (Pollard & Segall 1987;Healey et al. 2006) and the crack normal stress, which decays in inverse proportion to the square of distance from the crack centre, shows a characteristic bilobate pattern at each end. Hydraulic fracture stimulation treatments were simulated here simply by imposing an approximately steady-state (elliptical) region of pore pressure change around each injection point, which results in similar stress changes around the injection points (perforations). This aids the efficiency of the simulations and simplifies the interpretation of results by removing the inevitable uncertainty associated with the pressure distribution that develops around an injection location (Harper & Last 1990) in real situations. Ten sequential stages spaced by 100 m along a hypothetical lateral aligned with the azimuth of the minimum stress were simulated. Two-dimensional, 5 km × 5 km plan view models of the simulated shale reservoirs at a depth of 3000 m were constructed. First, the stress states were initialized in the Mohr-Coulomb continua with the required intact material properties and applied at the reservoir boundaries. Although a strong heterogeneity of elastic properties is typical of some if not most shale reservoirs, the option to include this was omitted to simplify interpretation of the simulation results. Once the reservoirs reached equilibrium the boundary conditions were changed to displacement-control before allowing the reservoirs to approach equilibrium. The FLAC (Itasca 2011) finite-difference formulation, used here, is fully coupled so that effective stress changes in response to changes of either fluid pressure or total stress. FLAC allows an elastic continuum to be populated with randomly distributed groups of zones with weakness planes of specified orientation and mechanical properties. After the 5 km × 5 km elastic continuum had reached equilibrium, three groups of vertical fractures were introduced to represent shear and extension fractures at different orientations occupying in total either 30% or 40% of the model zones. A single central fault inclined at 45° to the axes of applied maximum and minimum normal stress was added. A large number of cycles allowed the fractured medium to converge on an equilibrium state before either a transient increase of pore pressure or multistage stimulation. Effective stress computations were conducted without groundwater flow but with a change of total stress automatically implemented (assuming Biot's constant = 1) when a fluid pressure change is imposed by the user. In this mode, FLAC correctly represents the reservoir deformations induced by imposed pore pressure changes (such as during stimulation).
As noted, to allow the effective stress threshold hypothesis to be tested, each reservoir was subjected to two incremental, uniform increases in pore pressure each equivalent to 4.42 × 10 −4 MPa km −1 (0.01 psi ft −1 ). The reservoir was allowed to approach equilibrium at each stage and after the initial pore pressure had been restored. Therefore, for each shale reservoir base model, three states of identical mechanical properties and fracture distributions were available for fracture stimulation, one in a state of marginal equilibrium and two that had been exposed to lower effective stress 'thresholds' before restoring the initial reservoir pressure.

The two shale reservoirs
Two hypothetical shale reservoirs were constructed for the tests. The initial geomechanical properties and stress states of the reservoirs were very loosely based on published data but not intended to precisely represent any specific reservoir nor make any allowance for stress history. One reservoir (Reservoir A, representing a normal faulting environment) drew on some characteristics published for the Barnett Shale reservoir in Dunn and Mountrail counties, North Dakota, USA (Sorensen et al. 2010). ('Faulting environment' is used in the Andersonian sense.) The second hypothetical reservoir (Reservoir B) was chosen to represent a strikeslip faulting environment. This stress state is similar to that published for the UK Bowland Basin shale reservoir (Harper 2011). Sorensen et al. (2010) noted that most of the fractures in the Barnett Shale of Dunn and Mountrail counties were extension fractures relating to Laramide compression. The strike of the presentday maximum horizontal stress axis is inclined at c. 80° to that pertaining during the Laramide compression. A smaller number of mainly extension fractures are aligned with the present-day maximum horizontal stress azimuth; consequently, both sets are mainly subparallel to the present-day principal stress planes. Reservoir B was allocated a proportion of shear fractures inclined at a high angle to the principal planes to provide a contrasting reservoir condition. Figure 2 illustrates the distribution and orientation of fractures for the two reservoirs with a single, central fault.

Hydraulic fracture stimulation
Multistage stimulation in the central portion of each reservoir was simulated simply by sequential treatments at intervals of 100 m, along a hypothetical horizontal well aligned with the applied minimum stress. The magnitude of the maximum pressure at the injector location was chosen to be larger than the minimum applied total stress and smaller than the applied vertical total stress. In each case, the Bottom Hole Treating Pressure (BHTP) was chosen to be equivalent to a gradient of 18.1 MPa km −1 (0.8 psi ft −1 ) and held constant. The stress change around the hydraulic fracture is sensitive to the fracturing net pressure but less sensitive to the shape of the imposed elongate, elliptical pressure distribution. For purposes of this comparison, the injection pressure distribution was calculated on the basis of an ellipse with half-axes of 100 and 20 m. The mechanical response of the reservoir in the vicinity of the stimulation stage and elsewhere typically combined dominantly mode 1 (extension) failure near the injector point and slip on suitably oriented pre-existing natural fractures (ubiquitous joint zones).
A fault of 170 m length, inclined at 45° to the applied maximum stress, was centred along the lateral (most pronounced in Fig. 2b) and the sixth stage was immediately adjacent to the (right-hand) border

Reservoir stress state
The stress distribution throughout the reservoirs is heterogeneous as a consequence of slip on a proportion of the fractures. Figure 3 shows 2 km profiles of the maximum stress 1.5 km from the boundary of the two reservoirs along a line perpendicular to the applied maximum stress. The stress varies substantially along the line profile, varying typically 7 MPa and up to 15 MPa in Reservoir B. These variations represent approximately 10% and up to more than 20%, respectively, of the effective vertical stress at this depth. The maximum stress in Reservoir A varies less in magnitude and frequency (Fig. 3) largely because many more of the fractures in Reservoir A are aligned close to the axes of principal stress and do not slip. These magnitudes are coarsely comparable with the variation in the magnitude of the minimum stress of 12% of the overburden stress reported for a single lateral in the Barnett Shale by Ketter et al. (2006) and with the range of 9% reported for the same quantity in the Preese Hall well in the Bowland Shale (Harper 2011).
Although these normal stress magnitudes are indicative of the stress heterogeneity, the variation in shear stress provides a clearer demonstration. Figure 4 shows the shear stress acting parallel to the axes of the initial, applied maximum and minimum stress. It is clear that the principal stresses have rotated away from their initial orientations (otherwise the shear stress would be zero along the whole traverse). A strong fluctuation of shear stress has developed in both shale reservoirs.

Stability of the natural fractures and fault
Numerous increments of slip were observed during the large number of cycles required to approach equilibrium after the fractures were introduced in the model. Similarly, when the reservoir pressure was increased by what would be regarded as a minor amount in most hydrocarbon reservoirs (an increase in gradient of 0.23 MPa km −1 (0.01 psi ft −1 )), a large number of slip increments were triggered throughout much of the reservoir. When the reservoir fluid pressure was restored to the initial condition (by decreasing the pore pressure gradient by the same amount), all slip ceased. The reservoirs were no longer in a state of marginal equilibrium.
The increase in stability that is brought about by a transient rise of pore pressure can be illustrated by the changing stability of the fault placed in the centre of the models. A good approximation to the stability of the fault can be quantified based on the normal and shear stresses in each zone and assuming common, mean values of cohesion, friction angle and orientation for all of the zones forming the fault. The resisting forces and disturbing forces were calculated for profiles in the plane of the fault. Figure 5 compares the stability of the fault in the two reservoirs, before and after the transient rises of reservoir pressure. Figure 5 shows that the fault was close to limit equilibrium (resisting forces = disturbing forces) or, in geotechnical terms, close to a factor of safety of unity, prior to the transient rise of pore pressure. The stability of the fault was significantly increased by the pore pressure transients. The magnitudes of the transients are small relative to the values of overpressure commonly inferred for shale reservoirs.
The stability of the fault, intentionally aligned with the maximum shear stress, might also be calculated by ignoring the reservoir history and interaction with the reservoir natural fractures. The disturbing and resisting forces acting on a plane inclined at 45° to the applied maximum and minimum stresses can be calculated using the same Coulomb frictional description and rock properties. The results, which are constant and unaffected by any previous changes of effective stress, are shown in Figure 5 as grey dashed lines. Unsurprisingly (for a fault aligned with the maximum shear stress) the calculations suggest that the faults are unstable ( Fig. 5a  and b). This is incorrect because the stresses around the fault have evolved during the process of convergence on equilibrium, notably relieving the shear stress.
Each reservoir showed a distinctively different response to fracture stimulation before and after transient pore pressure increases. Figure 6 shows the responses expressed as the distributions of increments of slip on natural fractures, for an area of 1.21 km 2 , with and without prior pore pressure transients. The incremental pore  pressure changes corresponded to gradients of 0.22 MPa km −1 (0.01 psi ft −1 ) equivalent to 0.68 MPa (slightly less than 100 psi). The fractures are either isolated or clustered in small groups (Fig.  2) and the contours of slip reflect the localized nature of the isolated fractures or small clusters, appearing as little more than dots of differing densities at the scale of Figure 6.
For the cases in which many fractures were in limit equilibrium (no prior pore pressure transient), numerous slip increments were recorded remote from the stimulation locations ( Fig. 6a and d). In these cases (corresponding to the Snow (1972) threshold strength) stimulation-induced slip events extended to the boundaries of the model. Figure 7 shows an example of the marginal equilibrium condition plotted with contours of the change of stress difference induced by the multistage fracturing. Slip on the reservoir fractures tends to occur in the regions of increased stress difference (corresponding to increased maximum shear stress) resulting from the 10 fracture treatments.

Snow's threshold state
These results illustrate the process whereby an unstable fault or fracture becomes stable as a result of slip and the associated changes of stress around the fault. Therefore, if a reservoir effective stress state has been more adverse in the past, the discontinuities in the reservoir will be more stable than would otherwise be the case. The results presented in Figures 6 and 7 show that a major change in the response of a reservoir to disturbance occurs at the minimum effective stress experienced in the past and in this sense they are consistent with the hypothesis of Snow (1972). Snow's hypothesis applied to uniform changes of pore pressure (whereas these results correspond to a series of localized changes of the effective stress state).
It is concluded that these numerical results are qualitatively consistent with Snow's hypothesis of the significance of a 'threshold state' for a uniform rise of pore pressure. This threshold effective stress state is the reservoir state in which many fractures and faults are in limit equilibrium with the result that the reservoir is highly sensitive to disturbance. If this state is not reached during rise of the pore fluid pressure such as analysed for water impoundment by Snow (1972), reservoirs would be expected to remain aseismic.
For reservoirs in a state of marginal equilibrium, the results given here demonstrate that slip on fractures can occur remotely from injection wells in locations where pore pressure change does not occur. The triggering of these remote slips can only be a result of long-range stress transmission. A state of marginal equilibrium could arise because one fracture is in a state of limit equilibrium. It can also arise because many fractures are in limit equilibrium. As the effective normal stresses are reduced, for example by increasing pore pressure, increasing numbers of fractures reach the limit equilibrium state. The spectrum of behaviour ranges from stress changes that result in local slip on a single fracture or fault at a distance from the source of the stress change, to widespread slip. As more and more fractures in a reservoir become unstable (the density of the unstable fractures increases) a state that is typical of percolation thresholds (see, e.g. Kaye 1989) could arise. However, the mechanical response of the reservoir has not changed radically. It is simply that very small changes of stress can be propagated to large distances in a fractured, elastic medium and the induced slip distribution can reach a percolating state. This, of course, may be important in terms of the distance to which fluid (e.g. gas in a shale gas reservoir) can be transported. Because a state of marginal equilibrium embraces a wide range of states, with a corresponding range of instability from one fracture to many fractures, the only threshold is the state in which the first fracture reaches limit equilibrium. Snow (1972) applied the term threshold to a state of widespread instability. His analysis assumed fractures of a common orientation whereas a distribution of orientations and strengths is more realistic. These considerations imply that a single threshold is a less helpful concept than that of a marginal equilibrium state provided it is recognized that the distribution of those fractures that are in limit equilibrium (including lengths), and their frictional character, is relevant to induced seismicity potential.

Hydraulic fracture stimulation
In contrast to the uniform pore pressure changes considered by Snow (1972), fracture stimulation involves heterogeneous changes of pore pressure and of effective stress. A second major difference between hydraulic fracture stimulation and reservoir impoundment is that the former almost inevitably requires pressures that exceed the minimum stress but the latter does not. Fracture stimulation treatments are designed to exceed the minimum horizontal stress to create new mode 1 (extension) fractures. Because most fracture treatments exceed the  minimum total stress, slip may occur on any suitably oriented planes of weakness within the region of pore pressure change. However, slip is commonly limited to the regions of pore pressure change or close by (Fig. 6). This is qualitatively consistent with the conclusions of Shapiro et al. (2011) provided the response to stimulation exemplified by Reservoir B (Fig. 6) is applicable. A potentially far more serious risk of induced seismicity arises in marginally stable (percolating) reservoirs because long-range stress transfer can occur. Characteristically, a process of long-range interaction is set up as fractures in limit equilibrium, scattered throughout the reservoir, react to small stress changes and then in turn themselves becomes sources of stress change. Stability cannot be achieved instantaneously. The convergence towards a new state of marginal equilibrium is spatiotemporal. (This can be illustrated by movies obtained from explicit codes such as the FLAC code used here.) Below the Snow (1972) threshold, slip on fractures is not triggered remotely. This is demonstrated clearly for Reservoir B (Fig.  6d-f). It is less clear for Reservoir A (Fig. 6b and c, showing some shear displacements on fractures to a distance of about three fracture half-lengths from the injection point). This is almost certainly a response to the much higher net pressure applied in the treatments in Reservoir A, which resulted in much greater stress changes along and around the lateral (horizontal well section). For example, the increase of normal stress along the lateral, induced by the 10 stages of stimulation, was <10% of the applied minimum stress (acting parallel to the lateral) for Reservoir B. For Reservoir A the corresponding stress change was more than twice this amount.

Introduction
As indicated by Snow (1972), regions exhibiting pronounced historical seismicity are likely to be marginally stable. In the absence of seismic activity and/or obvious recent deformation, the investigations required to quantify effective stress history, or the trend thereof, depend upon a wide range of geological disciplines.

An example of reservoirs at or close to marginal equilibrium
The Cusiana oilfield and adjacent Cupiagua gas condensate field are located in the eastern foothills of the Colombian Andes mountains. This is a setting that displays a strong record of both historical and current seismicity. These fields are in the immediate vicinity of a major strike-slip fault structure defining the western margin of the South American plate, which extends from Ecuador to Trinidad and beyond (Audemard et al. 2005). A surface microseismic network was installed over the fields in 1992 and was in continuous operation at least until 2008 (Osorio et al. 2008). The network was originally set up to provide input to a geomechanical model of the fields for risk assessment and well and facilities siting. Surface microarrays were installed and the results were interpreted by beamforming. In 2008, Osorio et al. (2008) reported that an average of 1000 microseismic events per month was being recorded. This is a well-known seismically active region and abundant exposures are evident near the fields showing features associated with recent fault movements, including faulted and uplifted river gravels intermixed with underlying sediment, fault scarps and a sag pond. The driving forces are associated with east-northeasterly motion of the Nazca plate, subduction, uplift and development of gravitational potential associated with the Andes mountains immediately to the west. The planes of weakness appear to be mainly bedding planes and east-verging thrusts. Several of the latter are usually intersected by the wells.  reported many examples of casing ovalization with the short axis of the deformed casing closely aligned with the thrust direction. Buckling of a pipeline had been experienced elsewhere by another operator in the region.
In one case, a swarm of 33 events occurred just under 12 h after the onset of gas injection and these events extended to a distance of more than 10 km from the injector. Instances of induced seismicity have been associated with most forms of drilling and production operations, including drilling, injection and drawdown, yet many other examples of such activities were not detectably associated with induced seismicity despite the monitoring network. This variability has been partly illustrated by Osorio et al. (2008). One explanation may be that slip occurs aseismically in many cases, which would be consistent with the long-term, progressive nature of casing ovalization (Markley et al. 2002). Another possible explanation is that some proportion of the volume of the reservoir and overlying sediments has recently experienced stress relief, which is consistent with states of self-organized criticality.
The striking evidence of recent near-surface displacements and the history of seismicity clearly imply a state of marginal equilibrium at this location. This implication has been reinforced by the recent seismicity and the results of a programme of microarray deployment, global positioning system (GPS) displacement monitoring, fault movement dating, historical seismicity evaluation based on records of building damage and observations of casing damage.

A region that has experienced more adverse effective stress states in the past
In contrast to the high energy flux environment of the fields, the Lake Ontario region of upstate New york has displayed a very low level of historical seismic activity (Dames & Moore 1978b). During the construction of the Nine Mile Point Unit 2 nuclear generating plant, several faults that had displaced Holocene sediments were exposed in the bedrock. To assess the stability of the site, it was necessary to Fig. 7. Contours of slip increments on natural fractures and the fault occurring during 10 stages of stimulation of Reservoir A without a prior pore pressure transient. Plan view of 5 km × 5 km area. Contour density is identical to that used in Figure 6. The continuous lines that are also shown crossing the figure are contours that define the boundary between increased maximum shear stress and decreased maximum shear stress. Slip is largely confined to the regions of increased shear stress. implement a comprehensive investigation of the mechanics of the faulting process, Holocene geological evolution of the site and effective stress history. The programme is described in more detail than the Cusiana and Cupiagua setting to outline the requisite combination of geological specialities needed to evaluate effective stress history.
This summary is abstracted from the Dames & Moore reports prepared for the owner of the power plants (Dames & Moore 1978a,b,c). Where data and interpretations have been subsequently published, these are referenced.
The Nine Mile Point site is located on the southeastern shore of Lake Ontario. Tectonically, the region was described as very stable since the beginning of Palaeozoic time and generally free of major tectonic structures (Dames & Moore 1978a). The rocks of the region have been uplifted and eroded since the end of the Palaeozoic era and estimates of the amount of uplift since Cretaceous time range from 2400 to 3000 m (Dames & Moore 1978a). The beds dip gently southward. During the Pleistocene Epoch, several advances of continental ice sheets occurred in New york State. Each advance contributed to glacial-isostatic downwarping of the Earth's crust, and during each interglacial stage the crust rebounded.
The investigation involved three main disciplines: structural geology, geomorphology and rock mechanics. It relied on the principles and methods of many specialized fields including stratigraphy, geophysics, seismology, mineralogy, geochemistry, photogrammetry, hydrogeology, geochronology, palaeontology and geotechnical engineering.
The largest fault structures at the site strike WNW-ESE and have steep dips. Dames & Moore (1978a) identified three episodes of fault deformation. The fault showing the greatest displacements was initially developed as a vertical, left-lateral strike-slip fault with less than 1 m of stratigraphic offset. During this first phase of deformation, the fault was mineralized by calcite and sulphides, and analyses of fluid inclusions from the calcite indicated temperatures ranging from 170 to 220°C (corresponding to a depth of burial of 3-3.5 km). During the second phase of movement, the fault propagated downwards as a normal fault, experiencing a displacement of less than 1 m. Again the fault zone was mineralized with calcite and sulphides. Fluid inclusions in the calcite yielded homogenization temperatures ranging from 116 to 73°C, corresponding to a depth of burial of c. 2-3 km. The third episode of deformation was characterized by reverse slip, bedding plane slip and dilation confined to the top 60 m below the bedrock surface (Szymanski & Laird 1981). Field evidence clearly indicated two phases of movement during this third episode, the first of which preceded the deposition of thin overburden Pleistocene and Holocene sediments and the second of which post-dated their deposition. The Holocene episode naturally became the focus of the investigation and no further discussion of the preceding two phases of fault movement, which had occurred much earlier in the history of the basin, is given here. Figure 8 illustrates the nature of the deformation.
A detailed analysis of the mechanics of the third episode of deformation revealed that the buckling and associated dilation were dependent upon both layer-parallel effective normal stress and shear stress (Szymanski & Laird 1981). A detailed mapping programme was implemented to investigate the age of the sediments overlying the fault (Fig. 8). Measurements of ground water pressure revealed low pressures in the bedrock in the vicinity of the fault. An accompanying rock stress determination programme (Dames & Moore 1987c) had indicated that the maximum principal stress is subhorizontal and calculations based on the overcoring results (Dames & Moore 1978a) suggested that it is slightly inclined to the bedding (consistent with the offset of excavation blastholes at bedding planes and with the distribution of bedding slip around a water intake shaft (Harper & Szymanski 1983)). Recognizing that crustal warping associated with the Quaternary ice loads would result in an evolution of the beddingparallel shear stress, and hence possibly influence the deformation of the Quaternary sediments overlying the bedrock, a programme to examine the nature of the Holocene geomorphological changes and present-day differential vertical crustal movements was implemented.
The primarily geomorphological branch of the investigation (Dames & Moore 1978b) included careful mapping of the sediments overlying the fault and collection of samples of wood, peat and shell material for 14 C dating. The mapping revealed that the deformation of the overburden (Holocene) sediments was expressed as fluidized flow structures and folds and faults that deform the fluidization structures. Fluidization structures are often associated with vibratory ground motion but it was recognized that they could alternatively indicate one or more previous phases of transient, high upward-directed fluid pressure gradient. Samples were taken from the unconsolidated sediment for palynological analysis, to aid in correlations from trench to trench, and for grain-size analysis and X-ray diffraction, from the different coloured till units, from  (1978a). The Lake Iroquois and overlying sediments are of Holocene age. the fault filling and from the material filling the surrounding fractures and joints. Age dating bracketed the age of the (deformed) sediments immediately above the tills overlying the bedrock between 12600 ± 400 and 11750 ± 260 a BP. A field reconnaissance of glacial and postglacial lake shoreline features was undertaken to aid in determining the nature and amount of uplift resulting from glacial unloading.
A survey of the literature indicated that, during the Pleistocene, ice moved over the North American area in at least four major advances. The Great Lakes evolved from lakes that developed in the ice-eroded basins. At various times, drainage from these lake basins was restricted by the ice. At the time when ice blocked the northeasterly drainage of the Ontario Basin via the St. Lawrence River, these proglacial lakes drained westward.
An isobase map constructed for the main proglacial shoreline revealed strong tilting associated with the postglacial rebound. Geodetic relevelling data were reprocessed to obtain the distribution of present-day vertical crustal movements and conformed, mainly qualitatively, to the postglacial crustal movements deduced from the geomorphological programme.
Combining the postglacial uplift data with the late Wisconsin and Holocene history, Dames & Moore (1978b) constructed a time versus water level curve for the Lake Ontario Basin (reproduced by Harper & Szymanski 1983). It was estimated that the Nine Mile Point site had been submerged c. 93 m below the lake level before the ice receded to the north and opened progressively lower outlets that finally allowed drainage to the St. Lawrence River. It may be inferred that progressive unblocking of obstructions in the outlet to the NE resulted in the lake level dropping in a series of steps. The water level then rose progressively to its present elevation as the lake level rose and tilted in response to glacial isostatic rebound.
Analysis of the evolution of the two Holocene phases of deformation (Dames & Moore 1978a;Szymanski & Laird 1981), in the light of the multidisciplinary investigations, led to the following conclusions. The first phase of deformation was attributed to the combined effect of the high horizontal normal stresses resulting from basin uplift and the changes of bedding-parallel shear stress caused by glacial loading. The second phase of movement was attributed to the reduction in bedding-normal effective stress, probably amounting to a negative effective normal stress (fluid pressure greater than overburden stress) caused by the glacially induced incremental decrease of the water level in the Ontario Basin. The tilting associated with glacial-isostatic readjustment could be at present maintaining a limit state on bedding planes at shallow depth (Harper & Szymanski 1981, fig. 6) and contribute to or drive the phenomenon of rock squeeze well known in the region (Harper & Szymanski 1983). The subsequent tilting-induced rising water level in Lake Ontario implies that the vertical effective stress acting on the shallow beds at the site is progressively reducing (assuming a hydraulic connection to the lake). However, the magnitude of the present vertical (bedding-normal) effective stress is substantially higher than that experienced during ice retreat at times when abrupt removal of obstructions to drainage to the NE resulted in phases of rapid lowering of the lake level. For these reasons it could be concluded, in view of the geometry of the buckle structure, the tectonic environment and the decaying intensity of the glacial effects, that the dilatant fault structure could neither propagate further downwards nor generate vibratory ground motion (Dames & Moore 1978a) under present conditions. This brief summary of the Dames & Moore investigation at the Nine Mile Point site illustrates the dependence of analyses of effective stress history on a wide range of specialist geological disciplines. Not all of these have been mentioned. Mostly, the procedures required for combining these disciplines into efficient and reliable interpretations of stress history have yet to be developed.

Marginal equilibrium
Regions that are in a state of marginal equilibrium are unlikely to extend throughout the whole of a shale reservoir. Shale reservoirs may be structurally heterogeneous and much more extensive than conventional hydrocarbon traps. It can be speculated that reservoirs that have been subject to a complex synsedimentary and tectonic history (as, for example, the Bowland Shale in the UK), might comprise both marginally stable regions and stable regions. It can be speculated that such a 'domain structure' would change both with time (depending on geological processes) and with engineering activities such as closely spaced hydraulic fracturing.
Another consideration relevant to the extent of percolating regions (in a state of marginal equilibrium), or domains, is the reservoir rheology. The simulations reported here assume an elastoplastic behaviour with a reasonable allowance for damping incorporated in the mechanics. Processes that have not been simulated might change the state of the reservoir. For example, if the rate of any time-dependent dissipation of stress exceeds the rate of stress accumulation, a limit equilibrium state of fractures may be unsustainable. Investigations of this effect could also be conducted numerically. They would require representative constitutive descriptions of the rock behaviour and an evaluation of the historical strain rate. Snow (1972) emphasized the significance of his threshold effective stress hypothesis in relation to surface water impoundment. The simulations presented here are consistent with an increasing potential for fracturing-induced seismicity as reservoirs approach a state of marginal equilibrium. The further a reservoir is from a previous minimum effective stress state, the more stable are the fractures (including faults) likely to be.

Analyses of effective stress history
Hydraulic fractures are designed to reduce the least principal effective stress to less than that previously experienced by the reservoir. Localized slip on natural fractures is possible in the immediate vicinity of the hydraulic fracture, but may be prevented by the increased hydraulic fracture-normal stress, consistent with the slip distribution shown in Figure 6a-c (high fracture net pressure simulations). Hydraulic fracturing of shale reservoirs involves increases of shear stress at locations where the stress changes are influenced by the bilobate pattern of tensile stress and not dominated by increased crack-normal stress (Fig. 7). This increase of maximum shear stress gives rise to a potential for slip on natural fractures and faults. Given that the stability of faults is dependent upon the effective stress history (Fig. 5), so is the potential for induced seismicity.
Events of magnitude sufficient to be detected by regional seismic monitoring networks are only rarely induced (Nicol et al. 2011;Frohlich 2012). This implies that the large majority of reservoirs have experienced a previous geomechanical state that is less stable than the present one. By analysis of the effective stress history, somewhat similar to the Kremasta example of Snow (1972) and illustrated by the Nine Mile Point example, we can reasonably expect to distinguish between stable reservoirs and reservoirs that are potentially susceptible to induced seismicity. In all but the most obviously marginally stable regions, such as the eastern Cordillera of Colombia, this can be accomplished only by relying on geological data collection and interpretative techniques. Even the Colombian example required specialist dating techniques applied to sediments at fault scarps to quantify the history of slip.
A previous state of marginal equilibrium will not necessarily be of immediate interest when assessing the sensitivity of a reservoir to injection. Major changes of far-field stresses and displacements, such as the 80° rotation of the stress tensor reported by Sorensen et al. (2010) may well mean that, for example, fractures that previously influenced the stability of a reservoir become stable. Quantification of effective stress history is unlikely to be a trivial process but is extremely informative. Computational mechanics codes, combined with present-day descriptions of reservoir geomechanical state, are now a powerful tool for the interpretation of geological information in terms of reservoir stability.

Evaluation of the risk of induced seismicity
The highest risk of seismicity arises in reservoirs in the marginally stable state because the possibility of long-range triggering must be expected. This does not imply that displacements on faults will be large enough to induce events of a significant magnitude to cause damage. (Evaluation of the magnitudes of possible events is outside the scope of this work.) In the case of a reservoir that has experienced a previous phase during which the geomechanical conditions were far more favourable to slip on fractures, the potential for induced seismicity can theoretically be estimated by assuming that pore pressure diffusion controls instability, as demonstrated elsewhere (e.g. Parotidis et al. 2004;Rozhko 2010). However, determination of the stability of a fault prior to stimulation is not straightforward. The examples given here suggest that the common practice of replacing a representative description of the heterogeneous stress state by a single stress tensor deduced from few in situ determinations and neglecting the stress history is likely to result in imprecise estimates of fault stability (Fig. 5). The inadequacy of representing the heterogeneous stress state of a fractured medium by a single stress tensor is identical to McKinnon's observations in relation to mining-induced seismicity (McKinnon 2006). Such estimates may be inherently conservative if the fault has been subject to previously higher shear stress or higher pore pressure. The stability of a fault is substantially dependent on its geomechanical history so that a less conservative analysis would require quantification of the reservoir effective stress history.

Conclusions
Numerical computations of two hypothetical shale reservoirs subject to simulated fracture stimulation have yielded results that are consistent with the recognition by Snow (1972) of the significance of effective stress history. Snow's threshold effective stress is found to be a narrow range of (effective) stresses, rather than a single value. The range corresponds to the natural distribution of fracture orientations and strengths that are representative of real reservoirs. Elastoplastic reservoirs that are marginally stable are characterized by long-range responses to the small stress changes that occur around hydrofractures. The risk of seismicity induced by hydraulic fracturing of such marginally stable faulted reservoirs is potentially considerably higher than for reservoirs that are more stable.
The simulations given here, confirming Snow's proposal of the significance of effective stress history, are a demonstration that elastoplastic rock typically has a memory of its geomechanical evolution. This is consistent with the behaviour of ordinary dynamical systems.
Stress states in fractured reservoirs are liable to be highly heterogeneous. Analyses of the stability of faults based on comparisons of some assumed present-day uniform stress state, or stress state derived from the average of few stress determinations, with estimates of fault strength, are likely to be prone to significant uncertainty. Determination by measurements of the precise stress distribution around a fault is impracticable in an oilfield context with the instrumentation available today. The stability of a fault is a function of the history of slip and is influenced by previous effective stress states. In the absence of a knowledge of the effective stress history, the precision with which predictions of induced seismicity can be made on the basis of pore pressure diffusion is decreased.
When multiple hydraulic fractures are emplaced, typical of multistage stimulation, gross changes of stress state may arise and result in some shear displacements outside the region of pore pressure change. In the region of pore pressure change around an induced fracture, natural fractures in limit equilibrium prior to the treatment may not necessarily slip because of the increase of normal stress perpendicular to the induced mode 1 fractures.
Analysis of the effective stress history of a reservoir has the potential to be a powerful aid to assessments of the risk of induced seismicity associated with hydraulic fracturing. This conclusion applies also to other analyses of fault slip potential (or observed microseismicity) associated with diffusion of pore pressure away from an injector, such as waste fluid disposal or water impoundment. Because the present stability of a fault is a function of its effective stress history, there are many other types of engineering risk that can be assessed by such analyses. Both detailed and more generalized assessments of stress history can be greatly facilitated by employing available computational mechanics numerical tools. Detailed analyses require a wide range of specialist geological capabilities. Because the influence of stress memory on the response to injection (and other) operations has been largely ignored, procedures for application of the necessary range of disciplines are immature. A satisfactory level of experience and corresponding efficiency will require innovative combinations of existing techniques of data acquisition and interpretation.