Abstract
Estimation of spatially varying permeability within the two-phase porous media flow plays an important role in reservoir simulation. Usually, one needs to estimate a large number of permeability values from a limited number of observations, so the computational cost is very high even for a single field-model. This paper applies a nonlinear multigrid method to estimate the permeability field within the two-phase porous media flow. Numerical examples are provided to illustrate the feasibility and effectiveness of the proposed estimation method. In comparison with other existing methods, the most outstanding advantage of this method is the computational efficiency, computational accuracy, and antinoise ability. The proposed method has a potential applicability to a variety of parameter estimation problems.
1. Introduction
This paper investigates the estimation of a permeability field for the following diffusion equation within the two-phase porous media flow:with the initial-boundary condition
Equation (1) originates from the partial differential equation system, which can model two-phase immiscible flow of incompressible fluids in a porous medium with zero gravity effectswhere is the porosity, is the wetting-phase fluid saturation, is the fluid viscosity, is the fluid source term, is the capillary pressure, is the wetting-phase fluid pressure, and are the absolute and relative permeabilities, and subscripts mean —wetting fluid phase and —nonwetting fluid phase.
in (1) is a nonlinear function, which can model the main characteristics of the nonlinearity connected with some of the coefficient functions in (3). The nonlinearities in these two coupled nonlinear equations for and are related to both of the dependent variables. The influence of is through , and the influence of is through , so both and are considered as the independent variables in .
Determination of the permeability coefficient of the diffusion term has shown significant potential in reservoir simulation. It can help make important decisions about the management of petroleum reservoirs, for example, recovery method, well location, fluid production, and injection rates. However, the permeability is generally modelled as a piecewise constant function; namely, it is defined by a single value within each reservoir simulator grid cell. Thus one is interested in inferring as many parameters as there are grid cells in the simulator. This number is frequently more than for a single field-model. There are principally two difficulties about this.
Firstly, there is no enough sufficient information in measurable reservoir quantities (well pressures and flow rates) to estimate the permeability field with such a high resolution. The usual solutions are (1) to reduce the degrees of freedom in parameter space, or (2) to use prior knowledge about the permeability to aid in the estimation, or both. Such issues are not the theme of this paper; the interested reader is directed to referenced texts [1–4] for solution.
Secondly, the computational complexity of this problem is very large, even though the number of parameters is modest corresponding to the attainable resolution implicitly given by the measurements. It may take several hours to complete reservoir simulation for a single field-model, and it needs many such simulator runs for the estimation of a permeability field. In addition, the attainable permeability resolution is expected to increase in the near future, because of the appearance of the novel data acquisition techniques, such as time-lapse 3D seismic surveys. Therefore, there will be an urgent need for parameter estimation methodology able to handle a larger number of parameters within a reasonable time span. The aim of this paper is to develop a parameter identification method to overcome this difficulty.
There are many articles [5–10] focusing on the effective methods for the estimation of coefficient within the linear elliptic equation, which can describe single-phase porous media flow with constant fluid density and viscosityLess work [11–14] has been done on the estimation of coefficient within the linear parabolic equation, which can describe slightly compressible single-phase flow with constant compressibility, viscosity, and porosityIt is worth noting that the coefficient is just the permeability when (4) and (5) are used to describe one-phase flow processes in porous media. For the estimation of a permeability field in the practical application, the interested readers are invited to refer to the bibliography [15–20].
Multiscale techniques have been widely investigated to reduce computation for the forward problems, such as the elliptic problems [21], the parabolic homogenization problems [22], and the Navier-Stokes equations [23], and then were extended to solve the inverse problems [24]. Wavelet multiscale method is a specific form of multiscale techniques, which has recently been used to solve various parameter estimation problems. Successful applications of this method include the Bayesian tomography [25, 26], the Bayesian formulations of emission tomography [27], the thermal wave tomography [28], the diffuse optical tomography [29], the velocity estimation problems of a two-dimensional wave equation [30], the permeability estimation problems of a nonlinear convection-diffusion equation [31], and the parameter estimation problems of partial differential equations [32, 33]. It is shown in these papers that the performance of iterative parameter identification methods is much better with the help of the multiscale techniques.
Multigrid method, as another specific form of multiscale techniques, originally attracted interest because of its abilities to remove smooth error components and reduce the computational requirements of large numerical problems. This method works by recursively operating on the data at different grids, using the ideas of nested iterations and coarse grid correction. Multigrid method has been mainly used to solve partial differential equations, such as elliptic equation [34], modified phase field crystal equation [35], fractional diffusion equation [36], and Helmholtz equation [37]. Subsequently, researchers also applied it to a variety of imaging processing problems, for example, optical flow estimation [38], shape-from-X [39], and image matching [40].
More recently, multigrid method has been widely studied to solve parameter estimation problems [41]. In earlier work, Ranganath et al. [42] showed that the multigrid method could be applied to positron emission tomography. McCormick and Wade [43] applied a multigrid method to linearized electrical impedance tomography, and Spiliopoulos et al. [44] used a multigrid approach to estimate geometric anisotropy parameters from scattered spatial data that are obtained from environmental surveillance networks. Borcea [45] proposed a nonlinear multigrid approach for imaging the electrical conductivity and permittivity of a body, given partial, usually noisy knowledge of the Neumann-to-Dirichlet map at the boundary. Ye et al. [46] formulated the multigrid method directly in an optimization framework and used the method to solve Bayesian optical diffusion tomography. Nash [47] gave a multigrid method, which can solve a broad class of discretized optimization problems. Importantly, both the methods of Ye and Nash are based on the matching of cost functional derivatives at different grids.
In this paper, we propose a nonlinear multigrid method to estimate the permeability field within the two-phase porous media flow. This method is well suited to this permeability estimation problem for three reasons. First, the numerical examples indicate that the nonlinear multigrid method converges much faster than fixed grid methods. This is very important for the permeability estimation problem within the two-phase porous media flow because it is inherently three-dimensional. Second, the optimization at each grid is done in the space-domain where positivity constraints are easily enforced, so the positivity constraints can be implemented well by multigrid method. Generally, it is important that positivity can improve estimation quality when the problem is ill-posed. In our problem, the permeability parameters to be estimated have only positive values. Finally, multigrid method can effectively overcome disturbance of local minima, and permeability estimation problem within the two-phase porous media flow results in a nonconvex optimization problem; therefore, this robustness about local minima helps insure that a good solution can be reached.
A key innovation of our work is that the nonlinear multigrid method is directly formulated in an optimization framework by defining a sequence of optimization functionals at decreasing grids. In order to ensure the good convergence of the method, it is essential that the cost functionals at different grids be consistent. To achieve this, the coarse grid functionals are adapted by a recursive method, which guarantees that the exact fine grid solution is always a fixed point of the nonlinear multigrid method. This method greatly simplifies the application of multigrid to this permeability estimation problem and can be generalized to the other parameter estimation problems.
2. Nonlinear Multigrid Method
The forward problem is defined as ; that is, a density field is derived from solving (1). We take some density values , where are some fixed points in the spatial domain, as the observation data. The permeability estimation problem is formulated as ; that is, a permeability field is estimated from the limited number of observations. The parameter estimation problem is typically ill-posed, so we estimate a permeability field by minimizing the following cost functional:where is the regularization parameter. Then, we derive a specific expression of the nonlinear multigrid method for optimization of the cost functional in (6), by starting with the two-grid case and generalizing it for the V-cycle and full multigrid cases.
2.1. Two-Grid Method
Let denote the permeability field at the finest grid, and let be a coarser grid representation of with a grid sampling period of times the finest grid sampling period. Generally, can be computed from by the formulawhere is the restriction operator.
We make the hypothesis that a cost functional at grid is needed to be minimized, and there is an initial permeability field , which approximately minimizes this cost functional; namely,We aim at computing the permeability field at the next coarser grid and using it to improve or correct the fine grid permeability field by the following formula:where is the prolongation operator. To compute the coarse grid permeability field , a corresponding optimization problem at the coarse grid must be formulated. To achieve it, we first choose a coarse grid cost functional , which approximates well to . Without doubt, it is very important to choose this functional, which depends on the details of the considered problem. But for the moment only assume that is a reasonable approximation to the finer grid cost functional. Due to the possible discretization errors, we then solve an adjusted coarse gird optimization problemwhere is a row vector, which is introduced to correct the errors in the cost functional. In fact, the term is similar to the so-called residual term in the traditional multigrid, whose function is to adjust the errors between coarse and fine grids.
The question now is how should be chosen. It is ideal that the following equality holds for all values of :where the left side is the adjusted coarse grid cost functional, and the right side is the fine grid cost functional evaluated using the corrected result of (9). If (11) holds true, these two cost functionals will reach their minima at the same value of .
Generally speaking, there does not exist any row vector that makes (11) true, because the difference between the left and right sides of (11) is nonlinear. But a suitable choice of can match the derivatives of the two sides when . This condition yields the key expression for as follows:where , are, respectively, the gradient vectors of , . This condition is essential to assure that the optimum solution of (10) is a fixed point of this two-grid method and is illustrated graphically in Figure 1. It is worth noting that in (12) the restriction operator , which comes from the chain rule of differentiation, actually plays the role of a prolongation operator because it multiplies the gradient vector on the right. Most importantly, (12) holds for any choice of restriction and prolongation operators.
(a) The derivatives of the two sides in (11) are different
(b) The derivatives of the two sides in (11) are equal
2.2. Multigrid Method
Multigrid method can be implemented by recursively applying the two-grid method, and two recursions known as V-cycle and full multigrids are used in this section. Both these two multigrid methods move back and forth through coarse and fine grid in characteristic patterns as shown in Figure 2.
(a) V-cycle multigrid
(b) Full multigrid
The V-cycle multigrid method is a straightforward generalization of the two-grid method. In fact, the coarse gird optimization problem in the two-grid method is recursively solved by another two-grid method; then the V-cycle multigrid method can be got, as shown in Figure 2(a).
The full multigrid method, as shown in Figure 2(b), is on the basis of recursive calls of both the V-cycle and full multigrids. This structure causes this method to initially move to the coarsest grid. We prolongate the obtained coarsest grid permeability field to the second coarsest grid and use it as the initial permeability field for the second coarsest grid optimization problem, which is then solved by a V-cycle multigrid. Repeat this process until the final permeability field is got at the finest grid.
3. Simulation Results
In this section, we report simulation results on four numerical examples in two dimensions. We perform the simulations on a 2.6 GHz PC with 4 GB RAM in MATLAB 2012b environment under Windows 7. In our simulations, the permeability field is estimated both with the nonlinear multigrid method (NMGM) and with the fixed grid method (FGM, i.e., the Levenberg-Marquardt method [31]). For the estimation, we use the parameters as follows:
In the first two examples, the spatial and temporal step sizes , , and are, respectively, chosen to be , , and 0.0025 s, and the initial permeability field is chosen to be constant . The artificial Gaussian noise with level 1% is added to the measurement data. In the last two examples. The spatial and temporal step sizes , , and are, respectively, chosen to be , , and 0.005 s, and the initial permeability field is chosen to be constant 5. The level of Gaussian noise added to the measurement data is chosen to be 2%. For the nonlinear multigrid method, the regularization parameters in the coarsest and finest grids are, respectively, chosen as and . For the fixed grid method, the regularization parameter is chosen as . This is because the initial permeability field is not very good; a regularization parameter value smaller than will lead to the divergence of the fixed grid method.
Example 1. In the first example, the nonlinear function , and the exact permeability field to be estimated, , is piecewise constantFigure 3 displays the exact permeability field, the estimated permeability field using the nonlinear multigrid method, and the estimated permeability field using the fixed grid method. To validate the merits of the nonlinear multigrid method, the relative errors and CPU times in this example are listed in Table 1.
By observing Table 1, it is clear that the nonlinear multigrid method is much faster than the fixed grid method. Furthermore, by comparing Figures 3(b) and 3(c), it appears that the estimated permeability field using the nonlinear multigrid method is much better than the one using the fixed grid method. Thus, we conclude that the nonlinear multigrid method is efficient, accurate, and stable.
(a) Exact permeability field
(b) Estimated permeability field using NMGM
(c) Estimated permeability field using FGM
Example 2. In this example, the nonlinear function . In oil reservoirs the permeability field usually has greatly large jumps, so the exact permeability field is piecewise constantFigure 4 shows the exact permeability field, the estimated permeability field using the nonlinear multigrid method, and the estimated permeability field using the fixed grid method, and Table 1 also lists the relative errors and CPU times in this example.
It should be clear from this example that even if is highly nonlinear, the nonlinear multigrid method remains efficient, accurate, and stable. In fact, when the permeability field has greatly large jumps, the permeability field estimated with the fixed grid method is not satisfactory. This may be seen by comparing Figure 4(c), which presents an oversmooth permeability field, with Figure 4(b), which presents an acceptable permeability field. Therefore, we conclude that the multigrid procedure is effective and even necessary in oil reservoir simulation.
(a) Exact permeability field
(b) Estimated permeability field using NMGM
(c) Estimated permeability field using FGM
Example 3. In this example, the nonlinear function . The exact permeability field to be estimated, as shown in Figure 5(a), is a vertical stratified medium containing two interfaces, where the permeabilities from left to right are, respectively, 3.91, 5.03, and 6.18. The estimated permeability field using the nonlinear multigrid method is shown in Figure 5(b), and the relative errors and CPU times in this example are listed in Table 2.
(a) Exact permeability field
(b) Estimated permeability field using NMGM
Example 4. In this example, the nonlinear function . The exact permeability field to be estimated, as shown in Figure 6(a), is the model of two anomalous bodies in a homogeneous medium with a permeability of 6.27. The anomalous bodies have the permeabilities of 3.56 and 4.79, respectively. Figure 6(b) shows the estimated permeability field using the nonlinear multigrid method, and Table 2 also lists the relative errors and CPU times in this example.
It is evident from these two examples and the ones previously discussed that the merits of the nonlinear multigrid method, such as the computational efficiency, computational accuracy, and antinoise ability, still exist for the high-resolution and complicated permeability fields.
(a) Exact permeability field
(b) Estimated permeability field using NMGM
4. Discussion
Estimation of the permeability field within the porous media flow is not only a very important research topic itself but also has applications in a variety of domains, such as reservoir simulation and oil and gas exploration. Various methods had been presented to solve the problem of estimating the permeability field within the one-phase porous media flow. Within the growing body of work concerning multiphase porous media flow, few studies have focused on the issue of estimating the permeability field within the two-phase porous media flow. This paper has proposed a new method based on multigrid to estimate the permeability field for the nonlinear diffusion equation within the two-phase porous media flow. In this nonlinear multigrid method, the cost functionals at different grids are dynamically adjusted to be compatible, so that the exact fine grid solution can be always a fixed point of the nonlinear multigrid method. In comparison with the fixed grid method, the main advantages of this approach are that it can dramatically reduce the computational complexity and improve the estimation quality. Moreover, by taking into consideration observation noise, the robustness of the method has been demonstrated. This approach works entirely in an optimization framework; therefore, it provides a sensible way to estimate the permeability field within the two-phase porous media flow, which has a clear potential applicability to a wide variety of parameter estimation problems.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (11271102), the Science and Technology Funds for the Universities in Hebei Province of China (Z2015039), the Fundamental Research Funds for the Central Universities (N142303011), the PhD Foundation of Northeast University at Qinhuangdao (XNB2015002).