Krylov-Secant methods for estimating subsidence parameters
in hydrocarbon reservoirs

Adolfo Rodriguez
Center for Subsurface Modeling
The University of Texas, Austin TX 78712
adolfo@ices.utexas.edu
Hector Klie, Mary F. Wheeler

It is a well known fact that under primary fluid production conditions, pore pressures of hydrocarbon reservoirs and aquifers tend to decline. This decline may lead to severe deformations inside and around the reservoir. The deformation observed at the surface is known as subsidence and can produce negative environmental effects as well as damage to surface facilities and infrastructures. Positively speaking, subsidence observations can be combined with inverse algorithms in order to assess reservoir behavior and detect non-depleted regions that may be subject to further exploitation. The estimation of subsidence usually requires the performance of flow simulations coupled to mechanical deformation. These simulations are computationally intensive and, for this reason, they are seldom performed.

In this work we introduce a Krylov-secant inversion framework for estimating the deformation produced as a consequence of pore pressure reduction in the reservoir. The formulation is based on a solution of the equilibrium equation where perturbations due to pore pressure reduction and elastic modulus contrasts are introduced. The resulting equation for the strain is given in the form of the Lippmann-Schwinger integral, i.e.,

$\displaystyle \qquad \qquad \qquad \qquad {\bf {e}}\left( {\bf {x}} \right) = -... ... {\bf {\alpha }}\Delta p} \right\}d{\bf {x'}}, \qquad \qquad \qquad \qquad (1) $

where $ \Gamma \left( {{\bf {x - x'}}} \right)$ is the modified half space Green's function; $ \Delta\bf {C}$ is the elastic module contrast; $ \Delta p$ is the magnitude of the pore pressure drop; and $ \alpha$ is the Biot's poroelastic constant, assumed here to be a tensor. Both $ \Delta\bf {C}$ and $ \Delta p$ are zero outside the reservoir but can be functions of $ \bf {x'}$. The integration in (1) is performed over the reservoir domain.

A Born-type approximation is implemented, where the total field is assumed to be the incident field by analogy to electromagnetic theory.

Based on the discretization of (1) the resulting inverse problem can be stated as the minimization of the following mismatch functional:

$\displaystyle \qquad \qquad \qquad \qquad \qquad \qquad \qquad \phi =\frac{1}{2} \Vert\hat{d}-d\Vert _W^2, \qquad \qquad \qquad \qquad \qquad \qquad \qquad (2) $

where $ d$ and $ \hat{d}\in C^n$ are the predicted and observed data vectors respectively, and $ W$ is a diagonal weighting matrix, based on the inverse of the covariance of the measurements that compensates for the noise present in the data. In the problem addressed here, the observed data are the subsidence observations at some points, while the predicted data vector $ d$ is determined through solution of the forward model (1). The minimization problem (2) is large and ill-posed, especially in the event of high heterogeneities and few measurements. To perform the inversion, a Newton-Krylov procedure based on Krylov-secant updates is proposed.

The Krylov-secant framework entails a recycling or extrapolation of the Krylov information generated for the solution of the current Jacobian equation to perform a sequence of secant steps restricted to the Krylov basis. In other words, the Newton step is recursively composed with Broyden updates constrained to the reduced Krylov subspace. This is repeated until no further decrease of the nonlinear residual can be delivered, in which case a new nonlinear step yielding another Jacobian system is performed.

The proposed framework includes dynamic control of linear tolerances (i.e., forcing terms), preconditioning, and regularization to achieve both efficiency and robustness. Furthermore, this approach may optionally accommodate the latest deflation or augmented Krylov basis strategies for further efficiency. The framework has been previously applied for the solution of several nonlinear PDEs under Newton-Krylov implementations, but the present work explores further issues with respect to inexact Gauss-Newton methods based on Krylov iterative solutions such as LSQR.

Numerical experiments indicate that the current approach is a viable option for performing fast inversion implementations. Comparisons are made against traditional quasi-Newton and gradient-based implementations. It is concluded that the proposed approach has the potential for application to electromagnetic, radar, seismic and medical technologies.