Читать книгу Geophysical Monitoring for Geologic Carbon Storage - Группа авторов - Страница 28

When a volume of fluid is introduced under pressure at some depth within the Earth, it changes the state of rock‐fluid system around the injection site. The nature of the change is determined by the temperature, composition, pressure, and flow rate of the fluid and the initial conditions within the host formation. Generally, the injected fluid will displace the in‐situ fluid, leading to pressure and, consequently, volume changes that depend on the compressibility of the system as a whole. These volume changes will lead to strain within the Earth that will be transmitted outward from the injection site, ultimately reaching the Earth's surface. Given sufficient strain at depth, the resulting surface displacement can be large enough to produce a significant signal, observable by modern geodetic instruments such as a SAR satellite. Such signals can produce valuable information concerning the source of the deformation.

Given an observable pattern of surface deformation, one can attempt to infer properties of the source generating the deformation. That is, one can invert the observed data to estimate parameters describing the source. One of the most important aspects of an inversion for volume change is how the source model is parameterized. In order to effectively represent the source, one must know its basic geometrical properties and boundary conditions. For example, if the injected fluid is confined to a porous layer, then it is important to include the impermeable boundaries of the layer. Or, if the injected fluid induces the opening of a tensile fracture, then the strike and dip of the feature will need to be specified, perhaps by a data‐fitting procedure. For a tensile fracture, the volume change is completely determined by the aperture change over the fracture. Thus, the nature of the elemental source will be different from the volume expansion of a grid block. In fact, the expansion of a grid block can be modeled by three orthogonal tensile fractures, one along each axis. Other, more complicated source combinations are possible, such as slip induced along a tensile fracture. A correct formulation of the problem requires knowledge of the geology as well as of the general stress conditions at depth. Furthermore, one must be willing to modify the formulation in light of new information. For example, at In Salah, a double‐lobed pattern of surface deformation indicated a tensile feature at depth, giving rise to a modified source model.

Under favorable circumstances, surface deformation may be used to image the migration of the injected and displaced fluid at depth and the geological features that are controlling it. Needless to say, making such inferences will require knowledge of the properties of the fluids involved and of the overburden. However, the nature of the injected fluid is usually well known and the structure and properties of the overburden are also constrained by seismic data and well logs. One does not typically inject fluid blindly into the Earth. For the sake of our discussion, we will assume that the Earth behaves elastically at a sufficient distance from the source and for the time intervals in question, a few months to a few years. Furthermore, we will assume that the elastic properties can be reasonably estimated from the available data. Other rheological models, such as viscoelasticity, are certainly possible and are not a barrier to the approach that we shall lay out here. However, for the conditions associated with the geological storage of greenhouse gases, an elastic overburden is a reasonable model.

Geodetic methods do not really allow us to image fluid flow directly. In particular, displacements at the surface and within the overburden are responding to the effective volume changes induced by fluid injection and withdrawal. Imagine a box around that part of the reservoir influenced by fluid injection or withdrawal. If we cut out the reservoir and simply applied the same displacements as those caused by the fluid flow to the boundaries of the box, the resulting deformation in the overburden, and on the surface, would be the same. So, in reality, one can solve only for such effective displacements or volume changes. Other assumptions, such as constitutive laws for the reservoir, are required if we are to interpret the volume changes in terms of fluid pressure changes, thermal expansion, or nonelastic processes. Therefore, we will formulate our inversion in terms of reservoir fractional volume change in order to minimize the number of assumptions that must be invoked.

In an elastic Earth, the displacement at the surface, and in the overburden, is linearly related to the volume changes within the source region (Aki & Richards, 1980; Vasco et al., 1988). Thus, one can write the calculated displacements as an integral over with source volume V as

(2.7)

where Δ v(y) is the fractional volume change at source location y. The quantity G _i (x, y) is a Green's function representing the i‐th displacement component at location x that results from a point volume increase at y. The Green's function encapsulates the physics of the propagation of elastic deformation from the source to the observation point, obtained by solving the governing equation with a point source. For a simple medium with sufficient symmetry, such as a homogeneous half‐space, it is possible to produce an analytic Green's function. For a more general medium, one must resort to numerical approaches in order to compute G _i (x, y). Even a layered medium requires a numerical approach to compute the Green's function (Wang et al., 2006).

Equation (2.7) constitutes the forward problem whereby the injected volume is specified and the displacement at the observation points is calculated. For computational purposes, the source volume is usually decomposed into a discrete sum of N elementary volumes, such as rectangular grid blocks, and the total volume has the basis function representation

(2.8)

Note that the volumes may be taken to be quite small and sources are often decomposed into one or more point‐like sources. Substituting the volume expansion (2.8) into the integral 2.7, and using the linearity of integration, results in the expression

(2.9)

It is assumed that the fractional volume change given by Δ υ _i is constant over the elemental volume V _j . The function G _ij (x) represents the integral of the Green's function G _i (x, y) over V _j ,

(2.10)

For InSAR observations, the data will consist of range‐change values, given by the projection of the displacement vector onto the satellite look vector l = (l ₁, l ₂, l ₃). Thus, the range change will be given by r = l ⋅ u and the projection applied to equation 2.9 gives

(2.11)

where the range change kernel M _j (x) is

(2.12)

If we write the collective range‐change data for all of the pixels as a vector r, then the corresponding set of linear constraints, each in the form of equation 2.11, may be written as a matrix equation

(2.13)

In the inverse problem, observational data are used to estimate the distribution of volume change at depth. The most common approach is a least‐squares formulation in which one seeks a model v that minimizes the sum of the square of the residuals

(2.14)

While it is possible to try to estimate the volume change for a full three‐dimensional source model by minimizing S, the solution will most likely be nonunique or poorly determined. That is, there will be many possible solutions and solving the least‐squares constraint equations will not lead to a stable answer. One can stabilize the solution by adding penalty terms to the misfit function S, representing aspects of the model that are considered undesirable, a procedure known as regularization. For example, the magnitude of deviations from an initial or prior model,

(2.15)

is often included as a penalty term. A term penalizing model roughness is also commonly used to regularize the inverse problem. Note that if one takes v ₀ equal to zero, the norm penalty will have the undesired effect of biasing the solution to have the largest changes at the shallowest depth. That is, in order to minimize the magnitude of the solution, most significant anomalies are placed close to the surface where they will have the greatest effect on the observations.

Given fluid injection and production data, it is often useful to constrain the changes in particular formations to try to honor the fluid volume information in addition to fitting the geodetic data. For example, given the net injected volumes, ν _i , in a specific set of grid blocks, B _I , we can add the penalty term

(2.16)

to the data misfit function S, where f is a scaling factor to account for the fact that the injection of a cubic meter of water leads to a change that is a fraction of the injected volume. If no such information is available, or if one is uncertain about how to scale the injected fluid volume to fractional volume changes at depth, it is still possible to penalize changes that are far from the well. That is, we can bias the solution to try to put the most significant volume changes near the injection site, where the pressure changes should be the largest. For example, if a well is located at x _w , then the penalty function takes the form

(2.17)

where Δ(x _w , i) is a function measuring the distance between the well and the center of the i‐th grid block.

One of the advantages of geodetic observations is the relatively dense temporal sampling in comparison with other geophysical methods. For example, SAR data may be gathered at weekly to monthly intervals, leading to a time series of range change for each scatterer. Hence, it is possible to image the time evolution of the volume changes and to try and relate these changes to fluid movement and hydrological properties such as permeability (Vasco, 2004; Vasco et al., 2008; Rucci et al., 2010). The sequence of volume changes within each grid block of the reservoir model can be used to estimate an arrival or onset time, that moment at which the volume of a grid block is changing most rapidly. As discussed in Vasco et al. (2008), the onset time can be related to the travel time of the pressure front initiated by the start of injection. Specifically, for an elastic medium and a sharp injection profile resembling a step function, the onset of the peak rate of volume change, T _peak , is related to the phase σ of the propagating pressure front according to (Vasco et al., 2000). For propagation governed by the diffusion equation, the phase is given by the solution of an eikonal equation,

(2.18)

where γ (x) is the inverse of the hydraulic diffusivity (Vasco et al., 2000; Vasco & Datta‐Gupta. 2016, p. 138). The nonlinear, first‐order eikonal equation is equivalent to the system of ordinary differential equations,

(2.19

and

(2.20)

where x(s) is the flow path in the porous medium and s is the distance along the path. These expressions may be used to find the hydraulic diffusivity, and consequently the effective permeability, within the reservoir. As was shown by Rucci et al. (2010), permeability estimates based upon diffusive travel times are primarily sensitive to the kinematics of the pressure propagation and not sensitive to the coupling between the magnitudes of reservoir pressure and volume change.