NO20110444A1 - Stress and fracture modeling using superposition principles - Google Patents

Abstract

Voltage and crack modeling using the superposition principle is provided. A system simulates linearly independent far-field voltage models for a subsurface earth volume, and calculates voltage, strain and offset values based on the superposition of independent voltage tensors. Based on the predetermined values, the system generates real-time restoration of paleo-voltage values or voltage, strain, and offset parameters at any point in the subsurface volume while the user varies the far-field voltage values. The system recovers one or more tectonic events, or a voltage tensor represented by a relationship between principal values and associated orientation, using fault geometry, wellbore data (crack orientation and secondary fault plane data), GPS, InSAR, folded and discarded horizons, inclination gauges, and inclination gauges. The system uses various geological data from seismic interpretation, wellbore readings and field observation to provide a variety of results, such as predicted crack propagation based on perturbed voltage field.

Description

RELATED APPLICATIONS

[0001] This patent application takes priority from U.S. Pat. provisional patent application 61/317,412 to Maerten et al., filed March 25, 2010, and from U.S. patent application 13/052,327 to Maerten et al., filed on 21 March 2011, both of which are incorporated herein by reference.

BACKGROUND

[0002] A fault can be considered a finite, complex three-dimensional surface discontinuity in a volume of soil or rock. Fractures, such as joints, veins, veins, pressure-dissolution seams with stylolites, etc., can be intentionally widened to increase permeability in formations such as shale where optimizing the number, location and size of fractures in the formation increases the yield of resources such as shale gas .

[0003] Stress, within continuum mechanics, can be considered a measure of the internal forces acting within a volume. The tension can be defined as a measure of the average force per unit area on a surface within the volume where the internal forces act. The internal forces typically arise between the particles in the volume in response to external forces applied to the volume.

[0004] An understanding of the origin and development of faults and the tectonic history of faulted areas can be achieved by relating fault orientation, slip direction, geological and geodetic data to the state of stress in the lithosphere. In traditional inverse problems, the directions of the distant principal stresses and the relationship between their absolute values are delineated by analyzing field data on fault orientation and slip direction derived from features such as striations on exposed fault surfaces.

[0005] Traditional methods for inverting stresses, which use measured scour strips and/or fault jumps, are largely based on the assumptions that the stress field is uniform within the rock mass containing the faults (it is assumed that there is no perturbed stress field) and that the fault jump ( slip on fault) has the same direction and "sense" as the determined far-field stress on the fault plane. However, it has been shown that slip directions are affected by: anisotropic fault elasticity caused by irregular slip line geometry; anisotropy in the fault friction (surface corrugations); heterogeneity in the stiffness of the host rocks; and perturbation of the local stress field mainly as a result of mechanical interactions between adjacent faults. Mechanical interactions resulting from complex fault geometry in heterogeneous media should be taken into account in the inversion of the stress. When this is done, determining the parameters of this paleostress in the presence of multiple interacting faults requires running a number of simulations, and therefore extremely long computation time to fit to the observed data. The traditional parameter space must be searched for all possibilities and for each simulation the model and post processes must be recalculated.

[0006] Equations of motion are typically not used in traditional methods, and perturbations of the local stress field during fault jumps are not taken into account. The mechanical role played by the faults in tectonic deformation is not incorporated explicitly in such analyses. Nevertheless, a fairly complete mechanical treatment is used in traditional inversion of paleostresses. However, the results could be greatly improved if more types of data could be used to better define the inversion.

SUMMARY

[0007] Stress and crack modeling using

the superposition principle is provided. An example system simulates linearly independent far-field stress models for an underground

soil volume, and calculates stress, strain and/or displacement values based on superposition of independent stress tensors. Based on the pre-calculated values, the system can generate real-time results, such as recovery of paleostress values or stress, strain and displacement parameters for any point in the subsurface volume while the user varies the far-field stress value. The system can recover one or more tectonic events, or a stress tensor represented by the relationship between principal values and associated orientation, using only fault geometry, borehole data (comprehensive information on fracture orientation and secondary fault planes), GPS, InSAR, folded and faulted horizons, dip gauges, slip - and slip lines (slip and slickenlines) in faults. The system can use different types of geological data from seismic interpretation, borehole readings and field observations to provide a variety of different results, such as predicted crack propagation, based on a perturbed stress field.

[0008] This summary is not intended to provide a complete description of stress and crack modeling using the superposition principle or to provide an exhaustive list of details and elements. A detailed description with examples of designs follows.

BRIEF DESCRIPTION OF THE DRAWINGS

[0009] Figure 1 is a diagram of an example of a stress and crack modeling system.

[0010] Figure 2 is a block diagram of an example of a

computing environment to perform stress and crack modeling using the superposition principle.

[0051] Figure 3 is a block diagram of an example stress and crack modeling engine.

[0052] Figure 4 is a block diagram comparing paleostress restoration methods.

[0053] Figure 5 is a diagram of an example of a method applied to crack and conjugate fault planes using data sets without value information.

[0054] Figure 6 is a diagram comparing an example of stress inversion for a normal fault and a thrust fault.

[0011] Figure 7 is a diagram of cost functions for the voltage inversions in Figure 6.

[0012] Figure 8 is a diagram of an example model setup showing InSAR data points and a fault surface.

[0013] Figure 9 is a diagram comparing fringes from an original InSAR grid and an InSAR grid restored by an example stress and crack modeling method using the superposition principle.

[0014] Figure 10 is a diagram showing a view from above and a perspective view from the front of a plot of the sweep surface for an example of a method applied to the InSAR example in Figures 8-9, with sweep resolution indicated.

[0015] Figure 11 is a diagram showing results from an example of a method applied to a scenario with a flattened horizon.

[0016] Figure 12 is a flow diagram of an example of a method for stress and crack modeling using the superposition principle.

[0017] Figure 13 is a flowchart of an example of a method for stress and crack modeling using the superposition principle and a cost function.

DETAILED DESCRIPTION

Overview

[0018] Stress and crack modeling using the superposition principle is described here. Given various inputs, such as fault geometry and selectable or optional data sets or data targets, including fault jumps, dip or slip line directions, stress measurements, fracture data, secondary fault plane orientations, GPS (Global Positioning System) data, lnSAR (lnterferometric Synthetic Aperture Radar) data , geodetic data from surface inclinometers, laser ranging, etc, the example system can quickly generate or recover a variety of results. The systems and methods described herein apply the principle of superposition to fault surfaces of complex geometry in 3D (not just planar), and the faults are, by their nature, of finite dimension and not infinite or semi-infinite. The results are often rendered in real-time and may include, for example, real-time stress, strain and/or displacement parameters in response to a user request or an updated parameter; remote stress states for a variety of tectonic events; prediction of planned future fracturing; distinguishing pre-existing cracks from induced cracks; and other. The various inputs can be derived from wellbore data, seismic interpretation, field observation, etc.

[0019] The examples of systems and methods described below can be used for many different reservoir and

underground operations, extensive exploration and production operations for natural gas and other hydrocarbons, natural gas storage, hydraulic fracturing and matrix stimulation to increase reservoir production, water resource management, comprehensive development and environmental protection of aquifers and other water resources, capture and underground storage of carbon dioxide (C02) and other .

[0020] In an exemplary embodiment, a system uses a 3-dimensional (3D) boundary element method based on the superposition principle that applies to linear elasticity for heterogeneous, isotropic full- or half-space media. Based on pre-calculated values, the example system can determine a cost function to generate real-time results such as stress, strain, and displacement parameters for any point in a subsurface volume as the user varies the value of the far-field voltage. In one embodiment, the system uses only fault geometry and wellbore data, including e.g. crack orientation, information on secondary fault planes and/or in-situ stress measurement at hydraulic fractures, to restore one or more tectonic events or a stress tensor represented by a relationship between principal values and associated orientation. The system can use many different types of geological data from seismic interpretation, borehole readings and field observation to provide a variety of different results, such as predicted crack propagation based on a perturbed stress field.

[0021] Figure 1 shows an example of voltage and

crack modeling system 100. The system 100 can solve both common and more unusual geomechanical problems. The geometry of the faults is often known (and possibly measured fault jumps and imposed inequalities, such as normal force, tangential force etc. may be known). The user can typically have access to data from boreholes (fracture orientation, in-situ stress measurements, secondary fault plane), geodetic data (InSAR, GPS and inclinometer) as well as interpreted horizons. An example of stress and crack modeling engine 102 and/or equivalent

method examples can recover the distant stress state and tectonic regime for one or more relevant tectonic events, in addition to displacement discontinuity in faults, and for example estimate the displacement and perturbed strain and stress fields anywhere within the system.

[0022] Using the superposition principle, the example stress and crack modeling system 100 or engine 102 can perform each of three linearly independent simulations of stress tensor models in constant time regardless of the complexity of each underlying model. None of the models need to be recalculated. As indicated above, applications of the system 100 may include stress interpolation and crack modeling, reconstruction of one or more tectonic events, quality control of interpreted faults, real-time calculation of perturbed stress and displacement fields when the user performs parameter estimation, prediction of crack propagation, discrimination between pre-existing cracks and induced cracks, and a variety of other applications.

Example of environment

[0023] Figure 2 shows the system example 100 in Figure 1 within a data processing environment where stress and crack modeling using the superposition principle can be performed.

[0024] In the shown embodiment, a data processing device 200 comprises a component, such as the example of stress and crack modeling engine 102. The stress and crack modeling engine 102 is illustrated as software, but can be realized as hardware or as a combination of hardware and software instructions.

[0025] In the illustrated example, the data processing device 200 is communicably connected via sensing devices (and control devices) to a point of use, for example, an actual underground soil volume 202, a reservoir 204, a sedimentation basin, the seabed, etc. and associated wells 206 to produce a petroleum resource , for water resources management or for carbon storage services, etc.

[0026] The data processing device 200 can be a computer, a computer network or another device that has a processor 208, memory 210, data storage 212 and other associated hardware, such as a network interface 214 and a media station 216 for reading from and writing to a removable storage medium 218. The removable storage medium 218 can be, for example, a CD; DVD; flash memory etc. The removable storage medium 218 contains instructions which when executed by the computing device 200 cause the computing device 200 to perform one or more exemplary methods described herein. The removable storage medium 218 can thus contain instructions for realizing and running the example of stress and crack modeling engine 102.1 at least some parts of the stress and crack modeling engine 102 can be stored as instructions on a given realization of the removable storage medium 218, a removable device or in local data storage 212, for loading into the memory 210 for execution by the processor 208.

[0027] Although the illustrated example of stress and crack modeling engine 102 is shown as a program stored in memory 210, a stress and crack modeling engine 102 could be realized as hardware, such as an application specific integrated circuit (ASIC) or as a combination of hardware and software.

[0028] In this system example, the data processing device 200 receives incoming data 220, such as fault geometry and many other types of data, from multiple sources, such as wellbore measurements 222, field observations 224, and seismic interpretation 226. The data processing device 200 can receive many types of data sets 220 via the network interface 214, which can also receive data from the Internet 228, for example GPS data and InSAR data.

[0029] The data processing device 200 can calculate and generate modeling results, simulator results and control results, and a display controller 230 can output images of geological models and images and data from simulations to a display device 232. The images can be 2D or 3D simulation 234 of stress and crack results obtained using the superposition principle. The example stress and crack modeling engine 102 can also generate one or more graphical user interfaces for entering and displaying data.

[0030] The example of stress and crack modeling engine 102 can also generate or finally produce control signals to be used through control devices, e.g. drilling and exploration equipment or well control injectors and valves, for managing the reservoir 204, transport and delivery networks, surface facilities and others.

[0031] An example of system 100 may thus comprise a data processing device 200 and an interactive graphics display unit 232. The computer system as a whole may consist of simulators, models and the example of stress and crack modeling engine 102.

Example of engine

[0032] Figure 3 shows the example stress and crack modeling engine 102 in more detail than in Figure 1 and Figure 2. The illustrated embodiment is only one example embodiment for description purposes, to introduce details of and components of an engine that performs the example stress - and crack modeling using the superposition principle. The illustrated components are examples only. Different types or combinations of components than those shown may be used to perform the stress and crack modeling functions, and other or additional components may also be used. As indicated above, the example stress and crack modeling engine 102 can be realized in hardware or in combinations of hardware and software. Illustrated components are communicably connected to each other for communication as needed. The arrows are only shown to indicate process flow or data flow, as the components can communicate with each other as needed.

[0033] The example stress and crack modeling engine 102 illustrated in Figure 3 comprises a buffer for data sets 302 or at least access to the data sets 302, an initialization engine 304, stress model simulators 306 or at least access to the stress model simulators 306, an optimization parameter selector 308, a cost determination engine 310 and a buffer or output for results 312. These components are shown for description purposes. Other components, or other arrangements of the components, may enable different implementations of the stress and crack modeling engine 102. The functionality of the example stress and crack modeling engine 102 will be described in the following.

The operation of the system and engine example

[0034] Figure 4 shows different methods for restoring paleotension. Figure 4(a) shows a traditional technique that only uses fault jump distributions as data input, but gives unstable results. The technique in figure 4(a) is not capable of the fault jump distributions with alternative types of data input, such as GPS data and others. Figure 4(b) shows a traditional (2D) Monte Carlo method, but without optimization using the superposition principle. Figure 4(c) shows an example of a method as described here, including techniques that apply the superposition principle and dramatically reduce the model's complexity. The example method shown in Figure 4(c) can be realized by the example stress and crack modeling engine 102. For example, in one embodiment, the initialization engine 304 generates, through the stress model simulators 306, three pre-computed models of the far-field stress for a subsurface volume 202. For each of the three models, the initialization engine 304 precomputes, for example, displacement, strain, and/or stress values. The optimization parameter selector 308 iteratively scales the displacement, strain and/or stress values for each superimposed model to minimize the cost of the costing engine 310. The values optimized in real time in this manner are used to generate specific results 312.

[0035] In one embodiment, the three linearly independent far-field voltage parameters are: (i) north orientation, and (ii) & (iii) the two principal values. These far-field stress parameters are modeled and simulated to generate a set of variables for each of the three simulated models. Four variables can be used: fault jump, the displacement field at any data point or observation point, a strain tensor at each observation point, and the tectonic stress. The optimization parameter selector 308 selects an alpha variable for each simulation, i.e., a set of "alpha variables" for the three simulated stress models to serve as variable optimization parameters for iterative convergence toward values of these variables to minimize one or more cost functions, which will be described below. In one embodiment, the optimization parameter selector 308 randomly selects optimization parameters to begin converging the scaled strain, stress, and/or displacement parameters at the lowest possible cost. When the scaled (nearly optimized) parameters are considered to have the lowest possible cost, the scaled strain, stress and/or displacement parameters can be used to predict an outcome, such as a new tectonic stress.

[0036] Because the example method in Figure 4(c) uses pre-calculated values superimposed from its respective simulations, the example method or stress and crack modeling engine 102 can produce results quickly, even in real time. As noted above, the stress and fracture modeling engine 102 can quickly recover multiple tectonic events underlying current conditions in the subsurface volume 202 , or more quickly distinguish induced cracks from pre-existing cracks than traditional methods, or enable real-time parameter estimation while the user varies a stress parameter, or quickly predict cracking etc.

[0037] Although traditional paleostress inversion often uses a fully mechanical scenario, improving stress and

the crack modeling engine 102 traditional methods by applying several types of data. Data sets 302 to be used are generally of two types: those that only provide orientation information (such as cracks, secondary fault planes with internal friction angle, and scour stripes on faults, etc.) and those that provide information on sizes (such as

fault jump, GPS data, InSAR data etc). Traditionally, palaeostress inversion has been calculated using displacement measurements on the fault plane.

[0038] The traditional technique shown in Figure 4(a) performs an operation that goes over two steps: (i) determine the initial remote stress tensor aR for those fault elements that do not have relative displacement data and solve for the unknown relative displacements (bji Figure 4( a)); and (ii) use the calculated and known relative displacements to solve for aR (Figure 4(a)). An iterative solver is traditionally used which iterates between steps (i) and (ii) until convergence is achieved.

[0039] The traditional technique shown in Figure 4(b) is based on a Monte Carlo algorithm. However, this traditional technique turns out to be inapplicable since it requires a long computation time with a complexity of 0( n2 + p), where n and p are respectively the number of triangular elements modeling the faults and the number of data points. For a given simulation, a random far-field stress aR is chosen, and the associated displacement discontinuity b on faults is calculated. Then, as a post process in data points and depending on the type of measurements, cost functions are calculated using either the displacement, strain or stress field. In this traditional scenario, for hundreds of thousands of simulations, the best cost (close to zero) is kept as the solution.

[0040] The example method shown in the diagram of Figure 4(c), on the other hand, which can be performed by the stress and crack modeling engine 102, extends inversion to a variety of types of data and is a much faster modeling engine 102. As an example, a resulting fast and reliable voltage inversion described below. The different types of data can be weighted and combined. The stress and fracture modeling engine 102 can quickly recover the tectonic event(s) and displacement discontinuity in faults using different datasets and sources, and then obtain an estimate of the displacement and the perturbed strain and stress field anywhere within the medium, using data available from seismic interpretation, boreholes and field observations. Application of the superposition principle allows a user to perform parameter estimation in a very fast manner.

[0041] A numerical method for carrying out the examples of methods will be described in the following. Next, a reduced distance tensor used for simulation is described, and then the superposition principle is described. An estimate of the complexity is also described.

[0042] In one embodiment, a formulation used by the stress and crack modeling engine 102 can be performed using IBEM3D, a successor to POLY3D (POLY3D is described by F. Maerten, P. G. Resor, D. D. Pollard and L. Maerten, Inverting for slip on three-dimensional fault surfaces using angular dislocations, Bulletin of the Seismological Society of America, 95:1654-1665, 2005, and by A. L. Thomas, Poly3D: a three-dimensional, polygonal element, displacement discontinuity boundary element computer program with applications to fractures, faults, and cavities in the earth's crust, Master's thesis, Stanford University, 1995). IBEM3D is a boundary element code based on the analytical solution of an angular displacement in a homogeneous or inhomogeneous elastic whole or half space. An iterative solver is used for speed reasons and for parallelization on multi-core architectures.

(See for example F. Maerten, L. Maerten and M. Cooke, Solving 3d boundary element problems using constrained iterative approach, Computational Geosciences, 2009). However, inequalities cannot be used as they are non-linear and the superposition principle therefore does not apply. In the chosen code, faults are represented by triangulated surfaces with discontinuous displacement. The advantage is that three-dimensional fault surfaces provide a better approximation of curvilinear surfaces and curved boundary lines ("tip-lines") without introducing overlaps or gaps.

[0043] Mixed boundary conditions can be prescribed, and when stress boundary conditions are specified, the initialization engine 304 solves for unknown Burgers components. After the system is solved, it is possible to calculate, anywhere within the full or half space, displacement, strain or stress at observation points as a post process. Specifically, the stress field at any point of observation is given by the perturbed stress field resulting from fault jump plus the contribution from the remote stress. Consequently, it is not enough to simply determine the perturbed stress field as a result of fault jump. Furthermore, the estimate of fault jump from seismic interpretation is only given along the dip direction. Nothing is known along the strike direction, and a fully mechanical scenario is required to recover the unknown components of the slip vector as it will affect the perturbed stress field. Changing the imposed far-field voltage (orientation and or relative absolute values) will change the slip distribution and consequently the perturbed voltage field. In general, a code such as IBEM3D is suitable for calculating full displacement vectors for faults, and has been heavily optimized using a Jf matrix method. The unknown that remains for modeling purposes is the estimation of the far-field voltage, which must be given as boundary conditions.

[0044] In an example of a method that can be performed by the stress and crack modeling engine 102, a model consisting of several fault surfaces is subjected to a constant far-field stress tensor aR defined in the global coordinate system of equation (1):

Assuming a subhorizontal far-field voltage (however, the technology here is not limited to this case), equation (1) simplifies to equation (2):

Since the addition of a hydrostatic stress does not change aR, the far-field stress tensor aR can be written as in equation (3):

[0055] Accordingly, a definition of a far-field voltage is obtained with three unknowns, namely {alv a22, a12}.

[0056] The far-field stress tensor, as defined in equation (3), can be calculated using only two parameters instead of the three {an, å22, a12}. By applying a spectral decomposition of the reduced aR, equation (4) can be obtained:

where, as in equation (5): the far-field stress tensor aR is the matrix of principal values, and in equation (6):

is the rotation matrix about the global z axis (since a subhorizontal stress tensor is assumed).

[0057] By writing, in equation (7): equation (4) can be transformed into equation (8):

[0058] If one omits the scaling parameter a± as a result of Property 1 discussed below (when ax = 8 \ Property 1), aR can be expressed as a function of two parameters, 6 and k, as in equation (9):

[0059] These two parameters are naturally delimited by the equations (10): assumed that uniaxial remote stress applies for k > 10. For k = 1, a hydrostatic stress tensor is obtained, which has no effect on the model. Furthermore, using a lithostatic far-field stress tensor (which is therefore a function of depth z) will not invalidate the technique shown, and equation (9) can be written as equation (11): which depends linearly on z. The simplified tensor definition given by equation (9) is used in the coming sections to determine (Ø,fc), or equivalently {on,^,^}, according to measurements.

[0060] It is worth noting that even when 2-dimensional parameter space is used for the Monte Carlo simulation using ( 0, k ), three components are still used for the far-field voltage, specified by the parameters ( alt a2, a3 ). The transformations are given by expression (12): where, in equations (13):

and { a±, a2, «3)©r9'tt v©d equation (20) further below.

The principle of superposition

[0061] The example stress and crack modeling engine 102 uses the superposition principle, a well-known principle in linear elasticity theory, to recover the displacement, strain and stress at any observation point P using the pre-calculated specific values from linearly independent simulations. The superposition principle states that a given value / can be determined by a linear combination of specific solutions.

[0062] In the stress and crack modeling engine 102, restoring a far-field stress involves restoring the three parameters (all a2>a3). The number of linearly independent solutions used is therefore three. In other words, in equation (14): where ( alta2, «3)© are real numbers and <x(<0>(for i = ltil 3) are three linearly independent remote stress tensors. If F is chosen as the Greens function for strain , stress, or displacement, the resulting values of e, a, and u in P can be expressed as a combination of three specific solutions, as shown below. Thus, the strain, stress, and displacement field for a tectonic load is a linear combination of the three specific solutions, and is given by equation (15):

[0063] Correspondingly, the use of (of a2, a3) enables the restoration of the displacement discontinuities in the faults, as in equation (16): and any far-field voltage is also given as a combination of the three parameters, as in equation (17):

Complexity estimate

[0064] Changing aR normally requires recalculating the entire model to determine the associated unknown displacement discontinuities. Then, at any observation point P, the voltage is determined as a superposition of the far-field voltage aR and the perturbed voltage field resulting from sliding elements.

[0065] For a model formed by n discontinuous triangular elements, calculation of the state of stress at point P first requires solving for the unknown displacement discontinuities on triangular elements (for which the complexity is ø(n<2>), and then carrying out approximately 350/7 multiplications using the standard method. If one applies the superposition principle, on the other hand, it is not necessary to recalculate the unknown displacement discontinuities on triangular elements, and only 18 multiplications must be performed by the example stress and crack modeling engine 102. The complexity is independent of the number of triangular elements in the model, and is constant in time.

[0066] Some direct applications of the method examples will now be described, such as evaluation of deformation and the perturbed voltage field in real time while a user varies a parameter for the far field voltage. Estimation of paleostress using different data sets 302 is also described in more detail below, and also an example of a method for restoring several tectonic phases as well as a description of how the method example can be used for quality control during fault interpretation.

Real-time calculation

[0067] Before we describe the method for inverting paleostress, another method is first described for real-time calculation of displacement discontinuity in faults and the displacement, strain and stress fields in observation points while the orientation and/or absolute value of the remote stress is varied.

[0068] If the tectonic stress aR is given and three independent solutions are known, there exists a unique triplet (alf a2, a3) that fulfills equation (17), and equations (15) and (16) can be used.

[0069] In matrix form, equation (17) is written in the form shown in equation (18): or, in compact form, as in equation (19):

[0070] Since the three relevant solutions a(0) are linearly independent, the system is invertible, which gives equation (20):

[0071] In equation (20), A' 1 is precalculated by the initialization engine 304. Given a user-selected remote stress, aR, the stress and crack modeling engine 102 recovers the three parameters { a^, a2,( x.- s), and then the fault jump becomes and the displacement, strain and stress field calculated in real time using equations (16) and (15) respectively. To do this, the three current solutions of the displacement, strain and stress are stored at initialization at each observation point, together with the displacement discontinuity in the faults. In one embodiment, the example stress and crack modeling engine 102 allows the user to vary the orientation and absolute value of oR and interactively display associated deformation and perturbed stress field.

Paleotension inversion using datasets

[0072] As indicated above, the main unknowns in forward modeling for estimating the distribution of fault jumps, and consequently the associated perturbed stress field, are the orientation and the relative absolute values of the far-field stress aR.

[0073] If field measurements are known in some given observation points (e.g. displacement, strain and/or stress, crack orientation, secondary fault planes formed near major faults, etc), it is possible to restore the triplet (alt a2, a3) and therefore also to restore the tectonic stress aR and the associated tectonic regime (see Appendix A, below). The next chapter describes the determination method and the cost functions used to minimize the cost of different types of data sets 302.

Determination method

[0074] Using a Monte Carlo method makes it possible to find the parameters (alf a2, a3) which minimize the cost functions given three independent far-field voltages (see equation 15). Although (a1,a2,a3) implies a 3-dimensional parameter space, this space can however be reduced to two dimensions (namely the two parameters 6 and k), as the transformation is given by equation (20) and (Ø,<fc> ) —*■ (<r00, a01, a1±) —*■ ( alta2, a3), where, in equations (21):

(see also Figure 4(c) and Algorithm 1 for a detailed description). Consequently, the search method (the search for optimal parameters) is accelerated by reducing the parameter space by one dimension.

[0075] A simple sampling method can be performed by considering a two-dimensional rectangular domain where the axes correspond to 6 and k. This 2D domain is sampled randomly with np points, and the associated cost function (which will be defined in the coming sections) becomes used to determine the point of least cost. A refinement is then determined around the selected point and the procedure is repeated with a smaller domain. Algorithm (1) shows a simplified version of the example procedure, for which there is no refinement. The example of sampling method shown here can be greatly optimized with different methods.

Data set

[0076] The special feature of this method lies in the fact that many different types of data set 302 can be used to refine the inversion. Two groups of data will be presented in the following sections: the first includes only orientation information and the second includes information about displacement and/or stress value.

Without value information from the data sets

[0077] For open cracks (eg joints, veins, folds), the orientation of the surface normal to the crack plane indicates the direction of least compressive stress in 3D (<r3). Similarly, the surface normals of pressure dissolution seams and stylolites indicate the direction of greatest compressive stress ( a±). Using measurements of the orientation of fractures, pressure-dissolution seams, and stylolites enables the stress and fracture modeling engine 102 to recover the tectonic regime that generated these features.

[0078] At any observation point P, the local perturbed stress field can be easily determined from a numerical perspective by using three linearly independent simulations. Figure 5 shows the crack plane and conjugate fault plane. Figure 5(a) shows the orientation of

<r3i relative to an open crack (joints, veins, times) given by its normal vector n in 3D. Figure 5(b) shows the same as Figure 5(a), except for the orientation of <r3i relative to a join given by its projected normal vector n (eg, trace on outcrop). Figure 5(c) and Figure 5(d) show that

same as Figure 5(a) and Figure 5(b), except shown for a stylolite. Figure 5(e) shows the orientation of a2 °9 ffi1 relative to conjugate fault planes given by one of the normal n in 3D and the internal friction angle 6. The aim is to determine the best fit of the far-field stress aR, and therefore parameters ( ax, a2, a3) , given some orientations of open crack planes for which the normals coincide with the directions of minimum compressive stress a(\ P, or equivalently for which the crack's plane contains the greatest compressive stress (crj, as in Figure 5(a) and Figure 5(b).

[0079] By varying (alt a2, a3) one can quickly calculate the state of tension in any observation point P using the three pre-calculated models. The cost function to be minimized is given in equation (22):

there"." is the dot product, n is the surface normal of a crack plane and m is the number of observation points.

Minimization of a function of the three parameters is expressed by equation (23):

[0080] Correspondingly, for pressure dissolution seams and stylolites, the cost function is defined as in equation (22) using minimum compression stress <r3 as in equation (24) (see figure 5(c) and figure 5(d)):

Use of secondary fault planes

[0081] The orientation of secondary fault planes that develop in the vicinity of larger active faults can be estimated using a Coulomb break criterion, defined by equation (25):

where 6 is the angle of the fracture planes in relation to the maximum compressive principal stress a1 and \l is the internal

the coefficient of friction. Two conjugate fault planes intersect along <r2, and the fault orientation is affected only by the orientation of the principal stresses and the friction value.

[0082] The cost function is therefore defined by equation (26):

where a1 is the direction of greatest compression stress and <r2 is the direction of middle principal stress. The first term on the right-hand side of equation (26) maintains orthogonality between the calculated a2°9 normal to the fault plane, while the second term ensures that the angle between the calculated a± and the fault plane is close to 0 (see figure 5(e)).

Example 1: Normal and thrust faulting

[0083] Figure 6 shows a constructed example that uses an inclined plane fault as one of two randomly chosen conjugate fault planes. The angle of the two conjugate fault planes is shown for a normal fault (figure 6(a)) and a thrust fault (figure 6(b)). The dip azimuth and dip angle of each conjugate fault plane is used to perform the inversion and the internal friction angle is 0 = 30. The largest active fault is represented by the rectangular dip plane 602.

[0084] Initially, the model is delimited by a far-field voltage in some observation points 604, where the two conjugate planes are calculated using an internal friction angle of 30 degrees. Then, for each observation point 604, one of the conjugate fault planes is randomly selected and used as input for the stress inversion.

[0085] Figure 7 shows the cost function for the constructed example from Figure 6. Figure 7(a) shows the cost function for the normal fault, and Figure 7(b) shows the cost function for the thrust fault. In both cases, the recovered regional stress tensor, the fault step, and the predicted conjugate fault planes provide good agreement with the initial constructed model.

Use of scour strips on faults

[0086] For scour strips on faults, the cost function is defined as in equation (27):

where dg and å™ respectively represent the normalized slip vector from a simulation and the measured slip vector.

Data set containing value information

[0087] The value of displacements can be used to determine not only the stress orientation, but also the value of the remote stress tensor, instead of just the principal stress ratio.

[0088] The procedure for doing this is similar to that described above. Given equations (15) and (16), however, it is clear that there exists a parameter 8 for which the calculated displacement discontinuity in faults and the displacement, strain and stress fields in observation points scale linearly with the applied far-field stress. In other words, as in equation (28):

[0089] This results in the following property:

[0090] Property 1: Scaling the far-field stress by 8 e 9t scales the displacement discontinuity in faults and the displacement, strain, and stress fields in observation points by 8.

[0091] Based on this property, all measurements in data points are globally normalized before calculation and the scaling factor is specified (the simulations are also normalized, but the scaling factor is insignificant). After the system is solved, the recovered far-field stress, displacement, and stress fields are scaled back by a factor of 8^-

Use of GPS data

[0092] For a GPS data set, the cost function is defined in equation (29):

where up"- is the globally normalized measured power delivery change at point P from the horizon and uP c is the globally normalized calculated power delivery change for a given set of parameters { of a2, a3). The first term on the right-hand side in equation (29) represents a minimization of the angle between the two displacement vectors, while the second term represents a minimization of the difference in norm.

Use of InSAR data

[0093] When an InSAR dataset is used, there are two possibilities. Either the global displacement vectors of the targets are calculated using the displacement u along the direction of the satellite line of sight s, in which case equation (30) is applied:

and the same procedure used for the GPS data set (above) is used with the calculated uPc, or alternatively the calculated displacement vectors are calculated along the satellite line of sight, in which case equation (31) is applied: where"." is the dot product. The cost function is therefore given by equation (32):

Example

[0094] Figures 8, 9 and 10 show a constructed example using an InSAR data set. Figure 8 shows a model layout showing the InSAR data points 802 and the fault surface 804. Figure 9 shows a comparison of the edges from the original InSAR grid 902 and the restored InSAR grid 904. Figure 10 shows a plot of the sweep surface as a function of 6 (x -axis) and k (y-axis). On the left is a section seen from above 1002 of the plot and on the right is a perspective sketch seen from the front 1004 of the plot. The minimized cost solution 1006 in each figure is marked with a small white circle (1006).

[0095] To apply an InSAR data set, a forward model is run using one fault plane 804 and one observation grid (Figure 8) at the surface of the half-space (see the surface containing the data points 802). A satellite direction is selected, and for each observation point 802 the displacement along the satellite line of sight is calculated. Next, an example of the stress and crack modeling method described here is applied using the second type of InSAR cost function given in equation (32). Figure 9 compares the original interferogram 902 (left) with the recovered interferogram 904 (right). Figure 10 shows how complex the cost surface can be, even for a simple, constructed model. In one case, the cost surface was sampled with 500,000 data points 802 (number of simulations), and took 18 seconds on an average laptop with a 2GHz processor and with 2GB of RAM running on Linux Ubuntu version 8.10, 32 bit.

Use of a flattened horizon

[0096] When applying the mean plane of a given seismic horizon (flattened horizon), the stress and fracture modeling engine 102 first calculates the change in elevation for each point forming the horizon. The GPS cost function is then applied, for which only the wz component is given, which gives equation (33):

If pre- or post-folding of the area is observed, the mean plane can no longer be used as a proxy. Therefore, a smooth and continuous conformation surface must be constructed that must remove the fault deformations while preserving the folds. The same procedure as for the mean plane is then used to estimate the paleostress. In some cases, and before the continuous fitting surface is defined, it is sometimes necessary to filter the input horizon from noise that includes high frequencies, such as corrugations and bumps while preserving natural deformations.

Example

[0097] Figure 11 shows results from applying an example stress and crack modeling method to a constructed example using a flattened horizon. Figure 11 (a) shows a model layout showing the horizon 1102 and the fault surface 1104. Figure 11 (b) shows a comparison of the original dip displacement 1106 (left) and the recovered dip displacement 1108 (right). Figure 11(c) shows a comparison of the original stroke displacement 1110 (left) and the recovered stroke displacement 1112 (right). Figure 11 (d) shows original vertical displacement 1114 (left) from the flattened horizon (left) and the recovered vertical displacement 1116 (right) from the flattened horizon.

[0098] In the example shown in Figure 11, a complex designed fault initially bounded by a far-field voltage will consequently shift to accommodate the far-field voltage. At each point in an observation plane that intersects the fault, the stress and fracture modeling engine 102 calculates the resulting displacement vector and deforms the computational grid accordingly. Then, inversion is carried out using the fault geometry. After flattening the deformed computational grid, the change in elevation level of each point is used to delineate the inversion and to restore the previously applied far-field stress and fault jump and displacement field. A comparison of the original and inverted dip displacement (figure 11 (b)) and the strike displacement (figure 11(c)) shows that they coincide well (same scale). A good coincidence is also observed for the displacement field in the observation grid (figure 11 (d)).

Use of fall displacement information

[0099] When drop displacement data is used, the cost function is defined as in equation (34):

where bf is the measured drop displacement value for a triangular element e and be c is the calculated drop drop value.

Use of all available information

[00100] The example stress and fracture modeling engine 102 may combine the cost functions described above to better delineate stress inversion using all available data (e.g., fault and fracture plane orientation data, GPS data, InSAR data, flattened horizon data, dip displacement measurements from seismic reflection, scour stripes on faults, etc). Furthermore, data can be weighted differently, and each individual data unit can also support a weight for each coordinate.

Several tectonic events

[00101] In the case of several tectonic events, it is possible to restore the largest ones, e.g. those for which the tectonic regime and/or orientation and/or absolute value are significantly different. Algorithm (2), below, shows a way to determine different events from fracture orientation (joints, stylolites, conjugate fault planes) measured along boreholes.

[00102] After carrying out a first simulation, a cost is associated with each observation point that shows the reliability of the restored tectonic stress in relation to the data associated with this observation point. Cost zero means good reliability, while cost one means poor reliability. See Figure 7 for an example of a cost plot. By only selecting data points that lie below a given threshold value and running another simulation with these points, it is possible to extract a more precise value for the paleostress. Then, the remaining data points above the threshold are used to run another simulation with this paleostress state to reconstruct another tectonic event. If the graph of the new cost shows differences, the example of the procedure above is repeated until satisfactory results are obtained. During the determination of the tectonic phases, the observation points are classified according to their respective tectonic event. However, the chronology of the tectonic phases remains undetermined.

Quality control of seismic interpretation

[00103] It is useful to have a quality control (QC) method for interpreted fault geometries from seismic interpretation. The basic idea is to use the fracture orientations from boreholes to recover the far-field stress and displacement discontinuities in active faults. The calculated displacement field is then used to deform the initially flattened horizons. The geometry of the resulting deformed horizons can be compared with the interpreted ones. If clear deviations are identified (e.g. interpreted uplift and calculated subsidence), the fault interpretation is possibly incorrect. For example, an interpreted fault could fall in the wrong direction. Note that an unfolded horizon can be approximated with its midplane, as described above in connection with flattened horizons.

Conclusion and perspectives

[00104] The example stress and crack modeling engine 102 uses superposition properties that always apply in linear elasticity for real-time calculation of the perturbed stress and displacement field around a complex faulted area, and the displacement discontinuity in faults. Furthermore, the formulation realized by the example stress and fracture modeling engine 102 enables rapid paleostress inversion using multiple types of data, such as fracture orientation, secondary fault planes, GPS, InSAR, fault jump, and fault slip lines. In one embodiment, using only fracture orientation and/or secondary fault planes from wellbore, the stress and fracture modeling engine 102 recovers one or more tectonic events, the recovered stress tensor being given by the orientation of and relationship between the principal values. The example stress and fracture modeling engine 102 and associated methods can be applied across a wide range of applications, including stress interpolation in a complex faulted reservoir, fracture prediction, quality control of interpreted faults, real-time calculation of perturbed stress and displacement fields simultaneously with interactive parameter estimation, fracture prediction, discernment between induced cracks and pre-existing cracks etc.

[00105] In a variation, another important application of the stress and fracture modeling engine 102 and associated methods is evaluation of the perturbed stress field (and therefore the tectonic event(s)) for "shale gas" extraction. Since shale has low matrix permeability, gas production in commercial quantities requires fractures that provide permeability. This is typically done by hydraulic fracturing to form extensive artificial cracks around wellbore, and therefore requires a good understanding of how cracks will expand as a result of the perturbed stress field.

Examples of procedures

[00106] Figure 12 shows an example of a method 1200 for stress and crack modeling using the superposition principle. In the flowchart, the operations are summarized in individual steps. The example method 1200 may be performed by hardware or combinations of hardware and software, such as by the example stress and crack modeling engine.

[00107] In step 1202, linearly independent stress models for a subsurface volume are simulated.

[00108] In step 1204, stress, strain and/or

displacement parameters for the subsurface volume calculated, based on a superposition of the linearly independent stress models.

[00109] In step 1206, an attribute of the subsurface volume is iteratively predicted based on the pre-calculated stress, strain and/or displacement values.

[00110] Figure 13 shows an example of method 1300 for stress and crack modeling using the superposition principle. IN

the flowchart is the operations summarized in individual steps. The example method 1300 may be performed by hardware or combinations of hardware and software, for example by the example stress and crack modeling engine.

[00111] In step 1302, fault geometry for a subsurface soil volume is received.

[00112] In step 1304, at least one data set for the subsurface volume is also received.

[00113] In step 1306, three linearly independent far-field stress tensor models are simulated in constant time to generate strain, stress and/or displacement values.

[00114] In step 1308, a superposition of the three linearly independent far-field stress tensor models is calculated to produce pre-calculated strain, stress and/or displacement values.

[00115] In step 1310, a post-processing part of the method is initiated, which can calculate different real-time results based on the superposition principle.

[00116] In step 1312, optimization parameters for each of the linearly independent far-field stress tensor models are selected.

[00117] In step 1314, the pre-calculated stress, strain and/or displacement values are scaled with the optimization parameters.

[00118] In step 1316, a cost associated with the scaled pre-calculated stress, strain and/or displacement values is determined. If the cost is not satisfactory, the method returns to step 1312 to select new optimization parameters. If the cost is satisfactory, the method continues to step 1318.

[00119] In step 1318, the scaled strain, stress and/or displacement values are applied to the subsurface volume, e.g. in connection with a query about the underground volume or in response to an updated parameter regarding the underground volume.

[00120] In step 1320, a query or updated parameter regarding the subsurface volume is received, which initializes or initiates generation of the post process results in the real time results portion (1310) of the method 1300.

C onclusion

[00121] Although examples of systems and methods are described with wording that is specific to structural features and/or method steps, it must be understood that the object defined in the appended claims is not necessarily limited to the specific features or steps described. On the contrary, the specific features and steps are shown as examples of how the systems, procedures and structures for which protection is required can be realised.

APPENDIX A

For a given far-field stress tensor, the tectonic regime can be easily determined using the principal values at and the vectors $t (arranged in decreasing order). Since one axis is vertical, the tectonic regime is given by equation (35):

Claims (20)

1. Executable program contained on a computer-readable storage medium, wherein the executable program contains instructions which, when executed by a computer, perform a process comprising: simulating three linearly independent far-field stress models for a subsurface soil volume, calculating stress, strain or displacement values for data points in said underground soil volume, based on a superposition principle applied to the three linearly independent far-field stress models, and predict a stress attribute at said underground soil volume in real time, based on the pre-calculated stress, strain and displacement values. 2. Executable program contained on a computer-readable storage medium according to claim 1, wherein the three linearly independent far-field stress models are based on different data sets, where each data set may be weighted differently and each data set comprises one of fault geometry data, crack orientation data, stylolite orientation data, secondary fault plane data , fault jump data, slip line data, GPS (Global Positioning System) data, lnSAR (lnterferometric Synthetic Aperture Radar) data, laser ranging data, inclinometer data, displacement data for a geologic fault, or stress value data for the geologic fault. 3. Executable program contained on a computer-readable storage medium according to claim 1, wherein the predicted stress attribute comprises one of a stress inversion, a stress field, a far-field stress value, a stress interpolation in a complex faulted reservoir, a perturbed stress field, a stress ratio and associated orientation, one or several tectonic events, a displacement discontinuity in a fault, a fault jump, an estimated displacement, a perturbed strain, a distribution of fault jumps, quality control of interpreted faults, crack prediction, prediction of crack propagation as a result of a perturbed stress field, real-time calculation of perturbed stress and displacement fields at the same time as performing interactive parameter estimation, or distinguishing an induced crack from a pre-existing crack. 4. A computer-executable method, comprising: simulating linearly independent far-field stress tensors in constant time based on multiple types of geologic data for a subsurface volume, calculating precalculated strain, stress, or displacement values based on a superposition principle applied to the simulations, and iteratively minimizing a cost function to determine optimization parameters to predict a stress attribute of the subsurface volume in real time based on the pre-calculated strain, stress and displacement values. 5. Computer-executable method according to claim 4, where the predicted optimization parameters are used during one of: adaptation of a locally perturbed stress field to a far-field stress value in the underground volume in real time, or calculation of a perturbed stress, a displacement field and a displacement in faults around a complex faulted area in real time. 6. Computer-executable method according to claim 4, where simulating linearly independent far-field stress tensors in constant time makes it possible to apply the superposition principle to the simulations regardless of the individual complexity of each model associated with a given simulation. 7. Computer-executable method according to claim 4, where the several types of geological data are weighted differently and where each data unit is able to support a weight for each coordinate. 8. Computer-executable method according to claim 4, wherein the several types of geological data are obtained from several sources, comprising one of seismic interpretation data, wellbore data or field observation data, and wherein the several types of geological data comprise one of fault geometry data, fracture orientation data, stylolite orientation data, secondary fault plane data, fault jump data, slip line data, GPS (Global Positioning System) data, lnSAR (lnterferometric Synthetic Aperture Radar) data, laser ranging data, inclinometer data, displacement data for a geological fault or stress value data for the geological fault. 9. Computer-executable method according to claim 4, wherein the predicted stress attribute further comprises one of a stress inversion, a stress field, a far-field stress value, a stress interpolation in a complex faulted reservoir, a perturbed stress field, a stress ratio and associated orientation, one or more tectonic events, a displacement discontinuity in a fault, a fault jump, an estimated displacement, a perturbed strain, a distribution of fault jump, quality control of interpreted faults, crack prediction, prediction of crack propagation as a result of a perturbed stress field, real-time calculation of perturbed stress and displacement fields simultaneously with performing interactive parameter estimation, or distinguishing an induced crack from a pre-existing crack. 10. Computer-executable method according to claim 4, further comprising using a Monte Carlo method to iteratively minimize the cost function. 11. Computer-executable method according to claim 4, further comprising: using fracture orientation data from wellbore data to restore a far-field stress and displacement continuity in active faults, deforming initially flattened horizons with a displacement field from the far-field stress and displacement continuity, comparing the geometry of the deformed horizons with a interpretation of the horizons, and detecting deviations between the deformed horizons and the interpreted horizons as a quality control. 12. Computer-executable method according to claim 4, where the several types of geological data only comprise crack orientation data and/or secondary fault plane data from wellholes, and the predicted stress attribute comprises one or more restored tectonic events, where a restored stress tensor is given by the orientation of and the relationship between principal values. 13. Computer-readable storage medium physically containing instructions that when executed by a computer perform a process comprising: receiving fault geometry data for one or more faults in a subsurface volume, receiving at least one data set containing data points of collected data comprising fault orientation data, displacement data and/ or stress value data in the subsurface volume, simulate three linearly independent far-field stress tensor models in constant time to generate three sets of strain, stress, and/or displacement value components for each of the data points, calculate strain, stress, and/or displacement values for each data point based on a superposition principle applied to the three stress tensor models, iteratively determine optimization parameters to fit the precalculated strain, stress, and displacement values to a far-field stress value for each data point in the subsurface volume, and apply the optimization parameters on the pre-calculated strain- , the stress and/or displacement values to predict stress and displacement parameters in the subsurface volume in real time. 14. The computer-readable storage medium of claim 13, further comprising instructions for receiving a user-selected tectonic stress for the subsurface volume, and calculating the resulting stress and displacement parameters in the subsurface volume in real time. 15. Computer-readable storage medium according to claim 13, further comprising instructions for representing faults in the underground volume as a model comprising triangulated surfaces with discontinuous displacements. 16. Computer-readable storage medium according to claim 13, wherein iteratively determining the optimization parameters further comprises: selecting a set of optimization parameters, scaling the strain, stress and displacement values with the optimization parameters, determining a cost associated with the scaled strain, stress and displacement values, where, when the cost is not approximately minimized during an iteration, selecting a next set of optimization parameters, rescaling the strain, stress, and displacement values according to the next set of optimization parameters, and again determining the cost in a subsequent iteration, and when the cost is approximately minimized , apply the scaled strain, stress, and displacement values to predict a stress result in the subsurface volume. 17. Computer readable storage medium according to claim 13, further comprising instructions for predicting a tectonic stress in real time based on the pre-calculated strain, stress and displacement values. 18. Computer-readable storage medium according to claim 13, wherein at least one data set is weighted and wherein at least one data set comprises one of fault geometry data, crack orientation data, stylolite orientation data, secondary fault plane data, fault jump data, slip line data, GPS (Global Positioning System) data, lnSAR (lnterferometric Synthetic Aperture Radar) data, laser ranging data, inclinometer data, displacement data for a geological fault or stress value data for the geological fault. 19. Computer-readable storage medium according to claim 13, wherein each of the three linearly independent far-field stress tensor models is based on independent data sources, where each data source is obtained from one of seismic interpretation data, well-hole data or field observation data. 20. Computer-readable storage medium according to claim 13, wherein the predicted stress and displacement parameters comprise one of a stress inversion, a stress field, a far-field stress value, a stress interpolation in a complex faulted reservoir, a perturbed stress field, a stress ratio and associated orientation, one or more tectonic events, a displacement discontinuity in a fault, a fault jump, an estimated displacement, a perturbed strain, a distribution of fault jump, quality control of interpreted faults, crack prediction, prediction of crack propagation as a result of a perturbed stress field, real-time calculation of perturbed stress and displacement fields simultaneously with performing interactive parameter estimation, or distinguishing an induced crack from a pre-existing crack.