Analysis of Deformation Measurements for Reservoir Managemen

Analysis of Deformation Measurements for Reservoir Management — If reservoir deformation measurements can be analyzed to give consistent and coherent information on the volume changes and shear distortions taking place in the reservoir, data may be used for reservoir management and optimization of production and injection operations. Deformations may be measured at surface or at depth using a variety of technologies with different costs, ease of data collection, precision, areal coverage, and so on. The two most common techniques are the precision laser level survey, and the installation of geophysical tilt meters. Design of a suitable monitoring network for specific cases requires forward modeling using solutions that vary from spatial numerical integration of simple Green’s functions to a full nonhomogeneous three-dimensional finite element model. Oil & Gas Science and Technology – Rev. IFP, Vol. 57 (2002), No. 5


INTRODUCTION
Pressure and temperature changes in liquid-dominated reservoirs generate volumetric and shear distortions.Furthermore, injection or production of large amounts of solids will also generate a displacement field (Rothenburg et al., 1994) that can be sampled in various ways at the surface and at depth.A sufficient number of punctual measurements of adequate precision without systematic errors can, in principle, be analyzed to deconvolve the location and magnitude of the sources at depth.If the deformations at depth are transmitted to the surface purely through linear elastic deformations, in theory there is a unique solution.In practice, random error and a limited number of data points may make the solutions nonunique in the strict sense, but mathematically tractable and practically useful nonetheless.
If a successful deconvolution of reservoir deformation is achieved, it can help provide answers to many questions that arise in typical production management activities (Wallace and Pugh, 1992;de Rouffignac et al., 1995;Bruno and Bovberg, 1992;Houtenbos, 2000;Dusseault et al., 2001): -Where are produced fluids coming from and where are injected fluids going?-What are the geometric parameters of an induced hydraulic fracture?-Are there shear strain foci that exist and could develop into casing shearing difficulties?-Can the recompaction behavior observed in cyclic steam stimulation activities be quantitatively linked to deformation rates in the reservoir?-Are injected solids being contained within the target horizon?-Is exploitation activity leading to reactivation of largescale faults?Although local deformation sources may be plastic in nature (e.g.irreversible compaction or slip along a discrete plane), it is necessary to invoke elastic overburden behavior to achieve a solution.This is a corollary of d'Alembert's Principle: plastic strains at a point in a 3D setting must be translated into elastic strains a very small distance away from the point.The assumption that the overburden behaves elastically is always valid because overburden strains are small (ε ij < 10 -4 ), even if surface vertical movement, {∆z} s , is large (m).Much of the deformation that is transmitted to the surface is through direct deformation translation, rather than through strains that are related to stress changes.An elastic behavior assumption can be relaxed in the case of specific information that permits incorporation of shear bands or plastic volumetric deformation sites in the overburden (an active normal fault plane for example).This complex class of problems will not be addressed here.The article will also be restricted to simple cases for illustrative purposes: more complex problems can be tackled in practice.

PROCESS DEPTH AND DEFORMATION FACTORS
Because typical reservoir geometries involve reservoir widths (W) far larger than their thickness (H, Fig. 1), vertical movements at the surface are always larger than horizontal movements, usually by a factor of three or more.The reservoir width-to-depth ratio is a vital factor: if W/Z > is greater than 1, vertical deformations in the reservoir will be largely transmitted to the surface, regardless of depth, and they will be linearly proportional to the change in height of the reservoir, ∆H.
Where a volumetric strain (∆V) occurs in a volume where the lateral dimensions are small compared to depth (approximately < 1/20 th of Z, Fig. 2), a point deformation source is typically assumed.This "nucleus of strain" is defined by four parameters: X, Y, Z and ∆V.Green's functions exist that give the field of {∆d i } = {∆x, ∆y, ∆z} for the halfspace z + .In the simplest case, most commonly used in practice, only {∆z} is measured and inverted.The Green's function that relates {∆d i } to ∆V was developed by Mindlin and Cheng (referred to by Geertsma, 1973) for a linearly thermoelastic source (Dusseault et al., 1993).
The magnitude of induced vertical surface movement (∆z) in the case of a point ∆V source is ∆z ~ ƒ(1/Z 2 ).For a plane strain volumetric potential such as a rectangular source that is long with respect to depth, ∆z ~ ƒ(1/Z) except above the ends of the source.Where both reservoir dimensions exceed the depth, ∆z] s approaches ∆z] r , as with the Ekofisk (North Sea) reservoir with a W/Z ratio of ~1.8-2, and the surface deformation magnitude of ~9 m reflects ~85% of the 540 Rigorous deformation analysis falls into two categories: direct inversion and optimization of a forward model through error minimization.Three approaches are discussed: a direct inversion based on a nucleus-of-strain formulation, a multiparameter optimization of a single source function for hydraulic fracture analysis, and a displacement discontinuity forward optimization technique using a limited population of elements.Interpretation cannot be done in isolation: other data sources, including the project history, must be integrated to maximize the utility of the deformation analyses.As a final step, the data are used to help refine mathematical stress-flow reservoir models, which in turn become better predictors of deformation as well as oil production.measured reservoir compaction.In this case ∆z above the reservoir approaches the limit 1:1 (1/Z b , where b approaches 1).Thus, the deformation field magnitude depends strongly on depth and the geometric aspect ratios W/Z, L/Z, whereas it depends linearly on the source magnitude.There is a weak Poisson's ratio effect in cases where W > Z, but the Green's function contains no stiffness terms, therefore the surface distortion is independent of the stiffness characteristics of the overburden.As we will see, this is a simplification that can be overcome mathematically.
A shear distortion at depth gives rise to a general distortion field that is a function of the source magnitude and the geometric characteristics (Fig. 3).The decay functionalities are similar to those for the nucleus-of-strain, and are related to the source size, geometry and depth.For the inversion of data to yield the characteristics of the shear distortion at depth, it is mathematically convenient to represent a complex reality as a single displacement discontinuity (DD) shaped as a rectangle (Iwasaki and Sato, 1979).To fully define the simplest DD rectangle mathematically, ten parameters are needed (Fig. 4): -X,Y,Z, the global coordinates of a defined point on the DD rectangle; -L and W, the length and the width of the DD rectangle;   -Θ and Ψ, the dip and dip direction (dip azimuth) of the DD rectangle; and, -δx, δy, δz, the shear distortion magnitudes parallel to the L and W dimensions, and the volumetric strain normal to the DD plane.
If δx, δy, and δz are constants along the plane (a necessary simplifying assumption), a Green's function can be written that describes the general field deformations arising from the DD plane at depth (Okada, 1985).Even in the case of surface deformations only, these are relatively complex.
An hydraulic fracture source at depth can be treated in a similar manner (Okada, 1985;Davis, 1983).An elliptical penny-shaped fracture leads to a unitary Green's function, as does an ellipsoidal Sneddon-type crack or a Barenblatt circular horizontal crack.It is also possible to represent the fracture with a rectangular displacement discontinuity or some other discretized integral solution.The differences in the deformation fields arising from either assumption are of second-order magnitude; we carried out simulations using rectangular DD sources and elliptical crack sources of equal volume, and found that the differences in the shape of the surface deformation curve were so small that the deviation of superimposed plots could barely be detected visually (Rothenburg et al., 1994).

INFLUENCE FUNCTIONS FOR REAL CASES
For the layered elastic strata case (Jovanovich et al., 1974), an equivalent Green's function relating {∆z} to {∆V} is developed numerically using an axisymmetric finite element model and a unit volumetric strain at depth centered in the model (a ∆V "nucleus", Fig. 5).The elastic moduli for the FEM model are first estimated from seismic data and geophysical sonic and density logs, then the influence function can be refined in practice once actual deformation data become available.An identical approach is followed for flat-lying overburden strata that exhibit anisotropy with vertical rotational symmetry (E x , ν x = E y , ν y ≠ E z , ν z ).
In the more complex case of uniformly dipping strata, the finite element model must be three dimensional, but the stratigraphic model remains relatively straightforward.A limited number of mathematical nuclei are needed in the vertical plane-strain section normal to the strike direction.
In extremely complex, fully three-dimensional cases with no planes of symmetry, such as near a salt diapir or in a region that is intensely cut with normal faults, the stratigraphic model must be laboriously developed with seismic data and geophysical log data from available wells, and each volume element in the reservoir may have its own influence function.Nevertheless, influence functions coefficients are assumed to remain constant with time because they represent the overburden behavior, which is rarely changing with reservoir deformation (negligible strains, no pore pressure changes, etc.).Influence functions from finite elements.

THE "SNAPSHOT" APPROACH TO ANALYSIS
For reservoir processes involving slow changes in p and T, deformations occur slowly; it is not practical or mathematically feasible to attempt to track reservoir processes in "real time" in most cases.Furthermore, {∆d i } may be measured only episodically; e.g. a full {∆z} survey of a network of monuments can take up to two days, and is repeated perhaps every four weeks or every four months.Between two times t 2 and t 1 , the magnitude of {∆d i } must be large enough to warrant analysis: if the magnitude is too small, the random error of measurement can obscure genuine trends.It is also possible that some systematic errors could arise in a short time period, and if the nature and source of the error are obscure, it may not be possible for a reasonable correction to be made.In processes such as hydraulic fracture injection, fracture orientation changes can take place in a very short time, less than a minute, and if an analysis is carried out over too long a time interval, or from the beginning of the operation, the physical process changes will be missed, or at the least obscured.Finally, deformation measurements are often found to be particularly useful if they bracket an "event", such as a steam injection period (Bilak et al., 1991).
For these reasons, analysis is based on a time-differenced data set.Even in the case of "real time" hydraulic fracture mapping, data are chopped up into appropriate segments of time so that changes in behavior can be tracked.The actual deformation set being analyzed is the difference between two sets of measurements taken at a time interval that is Unit ∆V is applied centrally and its resulting surface deformation is an influence function used for inversion appropriate for the process: {∆z} t 2 -t 1 = {z(t 2 )} -{z(t 1 )}.If a continuous set of measurements or many frequent measurements are available, time series are plotted to try and identify changes in behavior: these become the beginning and the end of the "snapshot" (see example in Fig. 6).Alternatively, if an "event" is observed, such as the rupture of a casing (Dusseault et al., 2001), a set of deformation data can be collected as soon as feasible, differenced from the previous set, and analyzed to help determine the failure mode (e.g.shallow shear planes, steep shear movement along joints, thread popping from casing extension, bucking from compaction, etc.).
The segmentation process is generally an empirical process based on visual observations of time-series plots, on knowledge of the process, or on the available data set.

MEASUREMENT PRECISION, ERRORS AND INTERPOLATION
Precision is a characteristic of the nature of the measurement, and has to do with the resolvability of the source and the accuracy of the inversion in a relative sense, as well as an absolute sense.Detectability and resolvability (Fig. 7) are issues that must be studied before a specific technique to measure deformations is chosen.The expected volume changes are used in forward models to predict the deformations at the sites where they are desired.Attempting to achieve a degree of resolution at depth that is not warranted by the precision of the data, or by the "white noise" or random errors in the data, can lead to unnecessary expense and time delays.For example, the size of the reservoir blocks in the discretized reservoir at depth should not be too small, even if a great deal of data are available (as from a SAR (synthetic aperture radar) interferometric analysis, see Section 5.1, because small blocks in principle cannot be resolved mathematically with confidence, and too many blocks leads to underdetermined problems.

Precision
The measurement precision required for reliable mathematical inversion is a function of the signal magnitude.In most cases, a practical target is to have a precision of at least ± 5% of the maximum signal magnitude.This may not be easy to achieve in deep cases where ∆V is small.For example, in a 650 m deep steam flood case in Saskatchewan, only ~12 mm maximum deformation between "snapshots" was observed (average of ~8 mm), and the precision of the survey network was on the order of ± 0.7 mm.Clearly, the same basic survey precision under conditions of 200 mm vertical deformation would give an excellent precision.On the other hand, in cases of massive subsidence, say 200 to 400 mm between "snapshots", a precision of ± 10 mm would suffice for reasonably high-resolution inversion.

Errors
There are four types of errors that arise during measurement.
One is because of the fundamental precision limit of the measurement system; the second is random background error (white noise) arising from small perturbations and various random factors; the third type of measurement error is associated with a large perturbation of a measurement point, and the fourth type of error is the systematic error that is imposed on the entire array by some other process that may or may not have been identified.Time segmentation of real data.
The first two types of error can be easily handled; errors that are truly random (Gaussian) can be accommodated in analysis, provided that the measurement method is sufficiently precise.We carried out many inversions with theoretical data that were randomly perturbed, and found that random errors of 10% in {∆z} barely degraded the inversion process in a nucleus-of-strain method, and 25% random error still gave useful inversions, albeit somewhat degraded.
The third type of error may, for example, arise because a truck ran over a monument site, heavy equipment is parked right beside the site causing additional local deflection, or a poor installation results in a sudden settlement or increase in elevation during the time interval being analyzed.These errors are identified through quality control criteria during data processing, and human judgment is generally invoked in deciding whether to reject any specific piece of data.Outliers are individually examined and accepted or discarded depending on the past history of the project, the spatial location of the data site, and through comparison to the mean and the distribution of the deformation magnitudes.Highly improbable data must be removed.
Systematic errors, such as shallow but widespread ∆z from soil shrinkage in the dry season, are difficult to identify, isolate and remove.If there is the possibility of having two independent measurement systems based on different principles, one may be used to "check" the other, and help identify systematic error sources.In any case, any deformation monitoring output must be carefully and repeatedly evaluated for sources of systematic error.

Interpolations
Interpolation and extrapolation are used to generate intermediate "pseudo-data points" that are valuable in inversion analysis, and also to provide information as to the asymptotic nature of the decay in deformation with distance from the source.It is absolutely imperative that the interpolation scheme be rigorously faithful to the spectral characteristics of the deformation data, otherwise information is being lost (and new, incorrect information generated).No other option is possible if accelerated mathematical techniques such as FFT are to be used (Fig. 8).For the decay of the deformation curve asymptotically at considerable distances from the deformation source, a general decay function can be fitted to the system, thereby giving a means of generating pseudo-data points far from the centre of the system (e.g. at least at 1.5⋅Z from the edge of the deforming zone, Fig. 1), where the information content is extremely sparse.This extrapolation using known decay behavior helps achieve a better inversion that gives more precise total volume changes at depth, and also helps in the numerical inversion using Fourier transform methods.

General Methods
It is not practical to review in detail all methods of obtaining quantitative spatio-temporal deformation data; some methods and their characteristics will be briefly reviewed.Reference to Geophysical Journals (Journal of Geophysical Research, Geophysical Research Letters, Earth Science, GEOS, Seismological Research Letters, etc.) Johnson (1991) and Barends et al. (1995), various GPS and SAR conferences, etc.) is recommended.
Some methods, such as sea-floor depth gauges or pressure sensors (used in Ekofisk, Valhall and other cases of offshore subsiding reservoirs), underwater acoustic depth ranging, side-scanning sonar, or other approaches give useful data, but rarely are these data vectors precise enough for accurate inversion because the small second-order variations in the deformation field have been smoothed out by the standard error of measurement, or by the statistical error allocation method used to "smooth" or interpolate the data.
Earth satellite techniques involve two basic methods: direct measurements of differential displacement, and indirect methods using interferometric methods.Differential GPS (differential elevation) uses a network of fixed sites and literally millions of measurements over a short time period to give a deformation relative to a GPS site that is assumed to be spatially "fixed" (Dixon et al., 1997;Cacon', 1998).Precisions of sub-centimeter scale in horizontal positioning, and ~20 mm in ∆z are feasible, depending on the spacing and density of the differential GPS sites.
Interferometric SAR images can be made through using two sets of measurements taken at different times, several months or even years apart.The data are analyzed so that all other sources of displacement are accounted for, then the differenced sets give an interference pattern (in-phase and out-of-phase waves) spaced according to the wavelength of the radiation (Fielding, 1998), and deformations over the time interval can be deconvolved (Van der Kooij, 1997).Clearly, short wave-length probing gives a better chance at higher precision, but attenuation is also a larger problem.Interferometric measurements can, in principle, give centimeter scale resolution for {∆z} under ideal conditions; in practice, this technique can give a precision in the range of 8-80 mm, depending on the moisture content of the atmosphere and the ground reflectivity at the various times at which images are obtained.Two-frequency lasers with a very small frequency offset, operated from a "fixed" benchmark to measure the distance to a number of small fixed reflector targets, can give {∆x, ∆y} measurements of high precision, within ± 0.5 mm generally over long baselines.Because the magnitude of {∆x, ∆y} is far less than {∆z}, this method is seldom used except in geophysical deformation surveys, where lateral distortion across fault traces is being monitored.In principle, laser ranging could be used for {∆z} if a ranging site with good vertical offset from the measurement point array were available; in practice, this is never used.
Low-level aerial photography using multiple (up; to 10fold) overlapped images can be used to give precisions as good as 5-10 mm for {∆z}.Permanent visible targets (such as a 20 cm diameter white dot), protected from frost heave or other systematic error, can be installed an any number of sites across the field of interest, and relative elevations to a distant stable benchmark can be determined for any subset of highly resolved dots.
Various types of bubble or pendulum geotechnical inclinometers may be installed, giving precisions of 10 -6 to perhaps 2 × 10 -7 radians tilt change over a time interval (∆θ/∆t).
Downhole deformation can be measured through various borehole methods such as gamma-ray logging of radioactive bullets emplaced before the casing is cemented in place, casing collar logging using magnetic induction to estimate the axial ∆l of a casing joint (usually 10 m per joint), strain gages bonded on the outside of casing and directly wired to surface, piano wire extensometers to the surface, and so on.The precision of these methods depends on many factors, but values of ± 5 to 15 mm over the measurement baseline (∆L between points) are reasonable.Radioactive bullet logging was developed in the 1950's by Shell Oil in the massively subsiding reservoirs of Maracaibo lake, Venezuela.Casing collar logs were developed in the 1950's and 1960's in Long Beach and Wilmington, California, to measure compaction in production wells.Wire extensometers were used in measuring subsidence from deep groundwater withdrawal in the San Joachim Valley, California, as early as the 1960's, and were more recently used in Alberta by Imperial Oil Limited to measure cyclic heave and settlement atr depth in a 480 m deep large cyclic steam stimulation project.
Gravimeters (Grannell et al., 1979;Lambert et al., 1986), either placed in shallow protected sites or used as a precision geophysical logging system, can give deformation data, provided the density of the underlying strata is sensibly constant, which may not be the case when gas saturations are increasing through depletion, for example.

Laser Level Surveys
Precision laser level surveys are commonly used to determine the relative elevation of an array of monuments, perhaps in combination with GPS (Bitelli et al., 2000).Temperaturestable monuments are established by anchoring a steel bar at a depth of 6 m (Fig. 9) and filling the plastic protective casing with vermiculite.Multiple back-sightings and statistical error allocation are applied to the data set after removing spurious data (Verhoef et al., 1997) or systematic errors.Approximately ±0.7 mm precision can be routinely achieved over a baseline distance of 100 m.At least one "stable" benchmark for reference should be established > 2.5⋅Z from the edge of the area at depth that is of interest (Fig. 10).Sprectrum match for interpolated data.
Figure 9 Design of monuments for lase levelling surveys.
Figure 10 Array design for laser levelling surveys.
Because there are regions that contain more information than others, the spatial distribution of sites should reflect the density of information, and this is true for all of the methods mentioned that give punctual information.For example, in areas where d 2 z/dl 2  is large, there is much information, and the spacing of survey monuments should approach 0.20⋅Z; on the flanks, at a lateral distance from the ∆V zone being monitored, spacing drops off rapidly.No monuments are used farther than ~1.8⋅Z from the edge of the zone, and the spacing of the monuments can be optimized for usefulness using the forward model and trial pseudo data, minimizing some choice of error function.

Tilt Meters
The most powerful and precise deformation measuring approach is installation of precision tilt meters (Wright et al., 1998;Wright, 1998;Wright et al., 2001).These electronic bubble devices developed for geophysical purposes can achieve accuracies of 10 -9 radians in tilt change, equivalent to a 2-3 mm vertical displacement of one end of a rigid, 1 km long bar.This is the equivalent of about two orders of magnitude more precision than the best conventional inclinometer.Because of this great sensitivity and electronic output, "real-time" measurements are feasible, as sufficiently accurate data for inversion can be collected over a reasonably short time interval (providing the deformation changes are sufficiently large).
Tilt meters are particularly useful to monitor hydraulic fracture geometric parameters during liquid or solid-slurry injection.Because read-out is continuous, a set of data may be obtained at any time, and many differenced data sets generated to track short-term and long-term evolution of the injection process.For long-term deformation monitoring (many repeated measurements over a period of months or years), tilt meters may give problems because of zero-base drift that is negligible in the short term, but can accumulate in the long term to give a systematic error that swamps the long-term signal.Hence, if long-term analysis is needed, tilt meters can be installed along with a sparse levelling network.
Tilt meters are expensive, particularly so for a downhole tilt meter array placed with a wireline truck unit in a dedicated borehole.Thus, they typically will be used only for short-or moderate-term deformation monitoring.For hydraulic fracture or slurry injection monitoring, an array may consist of 12 to 24 tilt meters.

INVERSION OF NUCLEUS-OF-STRAIN METHODS
The coefficient matrix linking a strain nucleus array {∆V j } to the surface array measurements {∆z i } (Fig. 11) can be developed by writing down all the linear summation equations for ∆z i in the standard fashion for N nuclei (Geertsma, 1973): Repeating for all ∆z M leads to a standard integral (matrix) statement: However, this is typically an ill-conditioned problem where a small change in the input vector data can lead to a relatively large change in the specific ∆z values (although the total ∆V may match perfectly!).Some form of regularization, which is akin to conditioned removal of high-frequency noise, is needed.We use a form called Arsenin-Tikhonov regularization, where a Laplacian finite difference operator L is applied to orthogonal spatial data points equally spaced at the surface, and weighted with a coefficient α, in the following manner: Note that a higher-order symmetrical operator could be used, such as a 4 th -order biharmonic operator, or any other suitable spatial operator that does not introduce any nonphysical concept.In general, the coefficient matrix must be over-determined, and M, the number of surface measurement points, is substantially larger than N, the number of blocks in the discretized reservoir.
The choice of α is critical.If α is too large, an oversmoothed solution is given and the sum of squared errors is large; if α is too small, there is an excess of high-frequency variation in the solution (Fig. 12), but the sum of squared

Zone of interest at reservoir depth
Empirical choice of weighting coefficient.errors is smaller.We have found that the optimum value of α for an individual case may be based on an approach shown in Figure 13; choosing α at the point of sharpest curvature of such a plot almost invariably gives the most useful results.
Inversion of the weighted (smoothed) coefficient matrix is difficult if the deformation vector is scattered randomly in space.To accelerate the solution, orthogonal directions are divided into equally spaced distances, and the spectrallyfaithful interpolation used to generate pseudo-data points (Fig. 8).Furthermore, equally spaced distant pseudo-data points are generated using decay function extrapolations.It is now possible to convert the problem for flat-lying reservoirs with uniform overburden properties to two orthogonal fast Fourier transform solutions, which collapses the time for large problems from hours to less than a minute.This is possible because the problem for a flat-lying reservoir is really a 2D problem, and the FFT reduces the dimensionality by one, so that the solution becomes simply the addition and subtraction of scalar quantities, an extremely rapid process, even for very large arrays.
The nucleus-of-strain inversion approach is suitable for general determination of sources of volumetric deformation in a reservoir; it is not suitable for the deconvolution of data arising from a single concentrated planar or linear source, such as a hydraulic fracture plane or a heated linear source.Also, the nucleus-of-strain is not suitable for use to determine shear band sources, as there are no shear distortion terms in the basic Green's function.Therefore, in analyzing a subsiding reservoir, for example, it could give the spatial distribution of the volume changes with reasonable accuracy, but could not identify any slip distortion associated with reactivation of a high-angle fault, or bedding plane slip associated with compaction (Fig. 14).
In the general reservoir engineering case, volume changes arise over a distributed area; there is no strong singular deformation source, but a widely spread-out group of sources.In this case, the reservoir zone is discretized into equal-sided blocks, and it is assumed that the volume change in each reservoir block is a block-centered strain nucleus.Providing that the blocks are less than 1/20 th of the depth to the reservoir, the error generated by this discretization can be ignored.
Perhaps the strongest use of this method is to determine in a general way where injected fluids are going to, and where produced fluids are coming from (Bruno, 1997).However, this approach cannot distinguish the specific reason for the volume change; for example, if large amounts of solids in an aqueous slurry are being injected, volume changes may arise from the following processes (Fig. 15): -Emplacement of a discrete, solids-filled body, perhaps in a roughly planar shape such as expected for a conventional hydraulic fracture (Dusseault et al., 1997); Sum of errors squared Linking {∆V} to {∆z} arrays, nucleus-of-strain inversion.
Figure 12 Under-and over-smoothed inversion solution cases.
-The generation of shear dilation in the granular reservoir around the fracture, such dilation volume being filled with liquid from the injection process; -Volume changes associated with effective stress changes from the changing pressures and the injected volumes; -Volume changes associated with any temperature changes, such as cooling of warm rock by prolonged injection of cold fluids.If injected fluid simply is "taken away" rapidly by an extremely permeable zone without causing any change of pressure (∆σ′ ~ 0), there will be no significant volume change that could be detected by any monitoring method.This has turned out to be important in several cases of steam (and hot water) injection where the volume changes determined by analysis turned out to be a fraction (20-25%) of the injected volumes, indicating massive hot water losses through flow into a "thief zone".Note that this was identified well before it was verified by analysis of production data.

INVERSION OF HYDRAULIC FRACTURE SOURCE FUNCTIONS
It turns out to be difficult to find a direct matrix solution to invert for the case of an array of tilt meter measurements and a specific source function.A forward optimization technique is therefore used.Each tilt data point is a vector, with dip (tilt magnitude) and dip-direction obtained between two times t 2 -t 1 .The source function is a closed-form solution-e.g. a vertical or horizontal Sneddon or Barenblatt crack (Davis, 1983), circular or elliptical-or a semi-analytical solution-e.g. a rectangular displacement discontinuity (Iwasaki and Sato, 1979).The tilt data may be analyzed in terms of a single source function, or it may be decomposed into two or more source functions.
The procedure for the single source function (Φ) is described first.It is a classical forward optimization approach, and all of the techniques in such solution methods can be brought to bear as required (accelerations to convergence, etc.).
-Φ is stipulated to contain the injection depth (X, Y, Z) I , thereby constraining possible solutions to fractures intersecting the injection point.-A reasonable guess based on experience is made for the parameters that define Φ, such as ellipticity, inclination, mean aperture, etc. -A forward model is executed and the theoretical tilt at each site is computed.-The error is calculated; we use the sum of squares of the differences between real and predicted normalized tilt magnitude and direction, weighted differently.-A better estimate is made, and the process is iterated until a minimum solution is found.It is possible to minimize the predicted and observed error for any specific parameter.For example, assume that a DD solution was used for the forward problem, and that the operator wants a solution where the best possible estimate of the mean aperture is desired.A solution that minimizes the error of estimation for the aperture is then carried out by changing the weighting of the error sources in an algorithm that deals with specific parameters.In other words, the best DD solution for the aperture estimate will be sought.This process can be repeated to obtain solutions that minimize error for other parameters.(However, one may not assume that the best global solution is the linear combination of all of the individual best solutions because of interactions).
A decomposition model may also be used.Assume that only two simple solutions are available in the software: a vertical circular Sneddon crack, and a horizontal circular Barenblatt crack (Fig. 16).The goal is to seek the combination   of vertical and horizontal circular cracks that best represent the data set in linear combination.The procedure is approximately the following (several variations are possible): -The most closely fitting vertical circular Sneddon source is determined using minimization of a least-squares functional.-The contribution of this source to the tilt data is removed mathematically from the actual data, leaving a residual.-The residual is analyzed to see if it is more strongly composed of horizontal or vertical components.-The most closely fitting horizontal circular Barenblatt source is determined statistically in the same manner.-The procedure is repeated, giving a set of vertical and horizontal components that, in combination, can account for the observed deformation patterns.There exist many other ways of analysing a set of real data, and the goal of analysis may dictate which method is used.For example, if it is believed that vertical and horizontal sources dominate over a specific time interval, rather than generation of an inclined source, the decomposition method may be preferred, and a global minimization solution used.If it is relatively certain that shallow rising fractures are being generated, a single source function with a dip angle may be used.Furthermore, the single source solution may also incorporate slip along the inclined fracture plane; something that we know happens during the generation of shallow hydraulic fractures that tend to rise because the principal stresses are horizontal and vertical, but the parting plane is inclined.
Hydraulic fracture source functions of any type are not suitable for use in cases where the volume changes and shear distortions are distributed over a wide area and taking place at different times in different places.A more general solution is needed, one that is mathematically tractable, yet does not involve the laborious calculations associated with, for example, a forward finite element solution.

GENERALIZED SOLUTION USING DD SOURCES
It is not, nor will it likely ever be mathematically tractable to perform a generalized three-dimensional inversion that accounts for distributed volume changes plus localized planes of shear displacement with full solution freedom (no constraints) and full complexity.Some form of simplification is needed, and the freedom of the solution space must be constrained in manners that make physical sense.As mentioned in Section 1, it is possible to use a series of rectangular DD sources to help analyze a set of data, and appropriate simplifications and constraints can be naturally included.

The General Approach
A set of data may involve vertical and horizontal displacements at the surface or at depth, from various sources.Each additional high quality data point with only random error acts as a positive constraint on the solution space.
Assume that all the displacements measured arise from a series of DD elements, with ten unknown parameters for each element (Fig. 4).If M data points of sufficient precision are available, a solution can be pursued for N DD plane sources, where N < M/10 so that the problem is overdetermined.The solution is pursued using a forward optimization approach, as mentioned in the previous section.However, there now appear a number of additional issues that must be dealt with.Several of these issues will be discussed here, but it is not possible to present the approaches used for all sources of problems in analysis.
In the forward optimization, it is typical to have approximately 10-20 initial DD planes "seeded" in the zone where the deformations are expected to take place, depending on the number of data points and the degree of constraints on the solution.As before, a least squares error functional is minimized, and the error functional can contain terms of different statistical weight, to give more credence to certain types of information or certain locations.The choice of these weights is not arbitrary, but it is based on empirical factors and experience, and is thus more in the realm of "tradecraft" than in the realm of rigorous mathematics.
During the mathematical solution, there will be encountered a number of local minima that do not represent the best global solution.The procedure of forward optimization uses steepest descent gradient methods to go in the direction of the expected optimal solution, and includes algorithms to test whether a minimum is local, or likely to be global.To get out of local minima, the amoeba method, discussed in the optimization literature, is used.

Penalty Functions and Physical Constraints
In a generalized solution, the coordinates of the DD planes generally tend to rise out of the zone, to place themselves close under the measurement sites because there is simply too much freedom allowed in the analysis.To cope with this tendency, penalty functions are stipulated.For example, for a region of volume change, it is highly unlikely that the ∆V will take place in the shales above the reservoir, and certainly not in zones below the reservoir, therefore depth penalty functions are included.These functions give a strong penalty to any solution that tends to "wander" out of the dictated solution space.In other words, this is a method of introducing additional solution constraints that are based on our physical understanding of the problem.
The implementation of penalty functions based on depth is relatively straightforward, but it is a bit more difficult to decide on what are reasonable constraints on other parameters, and this depends on the nature of the process.For example, suppose the process is simply one of high-pressure water injection and fluids production from a series of wells.Not only can we introduce depth constraints on the DD sources, but it is possible to fix the center of the DD plane to the individual well locations, so as to link ∆V to individual wells.It may also reasonable to assume that if the injection and production pressures are not near fracture pressures, shear slip will be negligible, and in this case the shearing components could be set to zero (although this would almost destroy the goal of the DD solution).More reasonably, the dip angle of the planes could be set to zero if we are convinced that there will be no inclined slip planes, and δx, δy can be allowed to "float".
The point of this discussion is to emphasize that the complexity of the processes at depth is so large that as many additional constraints as possible must be introduced.These are not selected randomly: they are carefully thought out for each problem, based on the nature of the physical processes taking place at depth.

Human Intervention in the Solution
During the development of the numerical program that is used to seek a solution in the form of a number of DD sources, it became apparent that certain aspects could be accelerated by judicious human intervention in the optimization procedure.For example, if a DD source starts to elongate excessively or perhaps develop excessive shear distortion, and these factors have not been included explicitly in the solution penalty functions, it is highly beneficial on the solution time and the nature of the solution to partially reseed the solution space after some of the DD sources have been eliminated, or perhaps to eliminate "bad" sources altogether.Although it is possible to incorporate certain criteria in a general coding, it is also possible to examine the solution at several intermediate stages to see if the DD planes are well behaved and are likely to lead to a reasonable general solution, given the physics of the process being monitored.This brings in aspects of judgement and experience.

Example of a DD Solution
Figures 17 and 18 give an example of a solution for a 650 m deep case in Saskatchewan where a sequential cyclic steam stimulation process was being tried.In this pilot project, two rows were always being steamed, one or two rows were "soaking" and the remaining rows were being produced.The process of injection-soak-production moved from the bottom of the well rows to the top in steps, and then started again at the bottom.
Precision levelling surveys over 189 points were taken at two different times (4.5 weeks apart) and analyzed using both the nucleus of strain and the distributed DD methods.A few of the findings are listed here: -In the top of the well pattern, most of the steam was being lost out-of-pattern (verified during the production cycle).-Injection is almost invariably associated with shear movement (low σ′), production planes have no shear displacement whatsoever.-Injection planes are low angle (~12-15°) thrust planes.
-When the reservoir compaction ceased, oil flow from the producing wells virtually ceased (only water was then produced).-The production volumes could be linked (albeit not at 1:1) with the monitored volume changes.-The spatial distribution of the deformation solutions at depth correlated in an excellent manner with the production-injection activity.
In other cases, we have been successful in mapping the progressive growth of large volume emplacement of solids through high pressure injection, reservoir relaxation through pressure dissipation, and even the location of high concentrations of shear displacement that eventually led to the loss of several wells through shearing.Apparently, deformation monitoring and inversion can be used to understand the physics of processes, as well as for monitoring processes to try and optimize them.

FURTHER ISSUES
The classes of deformation problems addressed here are relatively simple, compared to what might occur in practice.For example, in complex stacked reservoirs, one may envision simultaneous injection and production into several zones, with volume changes and shearing taking place at different depths, combined with temperature and compaction effects.At some stage of complexity, it may become Figure 18 ∆V disposition from nuclei-of-strain solution, compared to the DD solution planes.impossible to deconvolve a data set in a meaningful way, no matter how many constraints can be introduced.However, for the great majority of real cases, including all those we have encountered to date, it is possible to make some meaningful inferences from the inverted solutions.It may be necessary to make some simplifications, but these are usually necessary in any case to cope with a limited data set, as in the case of tilt meters or levelling survey points.
The goal of a fully automated inversion program that can be used without human intervention seems, at first glance, worthy of pursuit.However, as appealing as this may be, it cannot account for a number of important aspects: -Experience and judgment are difficult to factor into a set of criteria or rules that can be exercised in the software environment.-In any of the methods, the optimum weighting criteria (e.g. for α in Section 6) may change as the influenced zone grows larger or changes shape, and this is also affected by judgment.-Criteria for a "bad" solution may differ from case to case, and these criteria may only become apparent after a number of analyses.-The physical processes giving rise to the distortions and deformations can change with time (for example, in a steam flood, T effects are dominant only after a long time).-Constraints may evolve with time in the process, and a rigid program of fixed constraints can result in bad solutions.
We therefore believe that the best approach is to combine the power of a number of solutions over different time intervals with the capacity of the human operator to integrate and draw subtle inferences.It is important to "re-test" the solutions and "re-visit" all the assumptions and constraints periodically to assure that conditions have not changed radically, something that could lead to a drift of the solutions farther and farther away from the "real case".
The nature of the interpretation of a solution is another issue that can be contentious.Both the nature of what constitutes an acceptable solution and the nature of the interpretation depends in part on what other supportive information is available, and may also depend substantially on the "view" of the interpreter.Because deformation and strain are geomechanical issues, we believe that the whole issue of deformation monitoring and analysis must reside in the hands of a rock mechanics engineer who understands processes such as bedding plane shear, plasticity, thermal expansion, compaction, and so on.This engineer must be reasonably sophisticated in reservoir analysis, but should be working in conjunction with reservoir engineers when interpretation and use of solutions takes place.
In the petroleum industry, there are almost always other sources of data that can help the analysis as well as the interpretation.For example, individual well data and the volume balance of liquids produced and injected are always known, and this provides additional constraints on the solution space.For example, if the net volume change for a pure injection phase in a project over ∆t is +10 000 m 3 , this limit can be introduced as an upper boundary constraint.In other words, if the solution to the measured deformation field exceeds this number, it is unlikely to be correct, and error sources must be sought.We advocate the use of deformation data in conjunction with other monitoring data such as produced and injected fluid volumes and locations, geophysical logging information, microseismic emissions, and so on.Limiting the interpretation to only the deformation data source is counterproductive: the analysis must be carried out in conjunction with other independent data sources to achieve maximum benefits.
The quantitative solution that is provided by analysis of a deformation problem is a statistically averaged solution in all cases.If a reservoir block is analyzed and shows a ∆V of 1000 m 3 , the specific vertical location and areal uniformity of the ∆V within the reservoir is uncertain; the choice of a certain reservoir block size constrains the answer to be an average within that block.Because reservoirs are usually thin with respect to the depth, these issues cannot be resolved without specific data collected by other means in the reservoir.This is both the power and the weakness of a remote method, compared, for example, with punctual measurements of p and T (or even ∆L) along a cased wellbore.The latter data are available only locally, and it is exceedingly difficult to extrapolate them far into the reservoir between wells.The spatially averaged data from deformation monitoring can be a great aid to this process, and the two data types used in conjunction (well data and remote spatially averaged solutions) are synergetic.
Finally, we note that the solutions that can be provided using the methods described herein provide powerful constraints to flow models (Carnec and Fabriol, 1999), and in particular, coupled geomechanics-flow models.The coupling between processes in geomechanics must be achieved through knowledge of the compliance of the system.For example, in the case of thermal expansion, a change in temperature leads to a change in volume, and this is predicted through the thermal compliance of the rock.Similarly, the effective stress compressibility (in the general sense) allows linking pressure changes to volume changes.Therefore, if the volume changes are measured, and if pressures are known or can be inferred, knowledge on the mechanical behavior of the reservoir can be extracted through back analysis, using finite difference or finite element techniques.Thus, these models can become far more powerful tools for reservoir engineers if they can be rapidly "calibrated" in situ.We believe that the only realistic way of achieving consistent and coherent calibration is through measuring and interpreting the deformation field.

CONCLUSIONS
Deformation analysis is largely a "solved" problem in the computational and mathematical sense.It can never be a fully solved issue because of the uncertainties involved and the statistical issues that arise in real cases (precision issues, limited data points, error sources, spatial and temporal resolvability of the processes, etc.What we find surprising, in fact, is that deformation data, which are now relatively economical to collect and analyze, seem to be largely ignored by the petroleum industry as a means of managing reservoirs.Particularly in cases of thermal recovery, large-scale subsidence, waste solids injection, and massive injectionproduction operations with attendant pressure changes, deformation analysis would seem to be an easy way to go towards better reservoir management.However, given the complexity of the class of processes that give rise to deformations, the usefulness of deformation analysis increases substantially with time as more data become available, as models become "calibrated", and as the physics of the processes are deconvolved.

Figure 4 A
Figure 4A displacement discontinuity rectangle.

Figure 15 Deformations
Figure 15Deformations can arise from several processes.

Figure 16 Two
Figure 16Two components of a general fracture.
Figure 17Production cycle and disposition of the displacement discontinuity solution planes.