A wavelet multiscale method for the inverse problem of a nonlinear convection–diffusion equation

Accepted Manuscript A wavelet multiscale method for the inverse problem of a nonlinear convection–diffusion equation Tao Liu

PII: DOI: Reference:

S0377-0427(17)30405-3 http://dx.doi.org/10.1016/j.cam.2017.08.016 CAM 11268

To appear in:

Journal of Computational and Applied Mathematics

Received date : 2 July 2015 Revised date : 17 January 2016 Please cite this article as: T. Liu, A wavelet multiscale method for the inverse problem of a nonlinear convection–diffusion equation, Journal of Computational and Applied Mathematics (2017), http://dx.doi.org/10.1016/j.cam.2017.08.016 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights (for review)

Highlights 1. Permeability identification of a nonlinear convection-diffusion equation is studied. 2. The finite difference scheme for nonlinear convection-diffusion equation is given. 3. A wavelet multiscale method is proposed for parameter identification. 4. This method is one of multiscale inversion method. 5. Numerical simulations testify the effectiveness of the wavelet multiscale method.

Manuscript Click here to view linked References

A wavelet multiscale method for the inverse problem of a nonlinear convection-diffusion equation Tao Liu∗ School of Mathematics and Statistics, Northeast University at Qinhuangdao, Qinhuangdao 066004, PR China

Abstract This paper is concerned with the problem of identifying the diffusion parameters in a nonlinear convection-diffusion equation, which arises as the saturation equation in the fractional flow formulation of the two-phase porous media flow equations. The forward problem is discretized using finite-difference methods and the inverse problem is formulated as a minimization problem with regularization terms. In order to overcome disturbance of local minimum, a wavelet multiscale method is applied to solve this parameter identification inverse problem. This method works by decomposing the inverse problem into multiple scales with wavelet transform so that the original inverse problem is reformulated to a set of sub-inverse problems relying on scale variables, and successively solving these sub-inverse problems according to the size of scale from the smallest to the largest. The stable and fast regularized Gauss-Newton method is applied to each scale. Numerical simulations show that the proposed algorithm is widely convergent, computationally efficient, and has the anti-noise and de-noising abilities. Keywords: wavelet multiscale method, inversion, Tikhonov regularization, nonlinear convection-diffusion equation, porous media flow, permeability 2000 MSC: 65F22, 65T60, 76R50, 76S05, 86A22

∗

Corresponding author Email address: [email protected] (Tao Liu)

Preprint submitted to Journal of Computational and Applied MathematicsJanuary 17, 2016

1. Introduction The immiscible displacement of oil by water in porous media with zero gravity effects can be described by the following set of partial differential equations (see [1]): ∇ · V = f1 (x, t), (1.1) V = −q(x)λ(x, u)(∇p − ρ(u)∇h),

(1.2)

φ(x)ut + ∇ · (f (u)V + fg (u)q(x)∇h) − ∇ · (q(x)N(u)∇u) = f2 (x, t). (1.3) Here, f1 and f2 are the injection and production wells, V is the total Darcy velocity, q is the permeability, λ is the total mobility of the phases, u is the saturation of the wetting phase, p is the global pressure, ρ is the density of the wetting phase, h is the height, φ is the porosity, and N is a nonlinear diffusion function. Eq. (1.3) is the fractional flow formulation of the mass balance equation for water, and f is a nonlinear fractional flow function, which is typically S-shaped. Furthermore, fg (u) = (ρ(u) − ρ0 (u))f (u)λ0(x, u), where ρ0 and λ0 , respectively, are the density and phase mobility of the non-wetting phase. Oil reservoir simulation based on the identification problem of permeability q(x) in Eq. (1.3) is a powerful tool to help make important decisions regarding the management of petroleum reservoirs, including selection of the type of recovery method, fluid production and injection rates, and well locations. Generally, this permeability identification problem has computational difficulties because the direct model is described by the solution of the nonlinear convection-diffusion equation, which is computationally demanding to solve. Furthermore, this problem is highly nonlinear, so fast iterative methods like Gauss-Newton or quasi-Newton require good initial estimates. As stated above, there is a pressing need for parameter identification methods able to solve these difficulties. This paper is concentrated on the further development of such parameter identification methods. A large amount of literature (e.g., [2-7]) is devoted to the development of appropriate methods for the recovery of q(x) within the linear elliptic equation −∇ · (q(x)∇u) = s(x), (1.4) 2

while much less work has been carried out on the parameter estimation in the linear parabolic equation ut − ∇ · (q(x)∇u) = s(x, t).

(1.5)

As for the identification of q(x) in Eq. (1.5), the interested readers are invited to refer to, e.g., [8-13]. Eqs. (1.4) and (1.5) can be considered as the descriptions of one-phase flow processes in porous media where q(x) is the permeability. For the permeability identification problems in the practical application, we refer the interested reader to, e.g., [14-17]. This paper investigates the identification of permeability coefficient q(x) in the following nonlinear convection-diffusion equation, which is closely relevant to two-phase porous media flow: ut +

∂ ∂ f (u) + g(u) − ∇ · (q(x)N(u)∇u) = s(x, t), ∂x ∂y

(x, t) ∈ Ω × (0, T ), (1.6)

with the initial condition u(x, 0) = u0 (x),

x ∈ Ω,

(1.7)

and the boundary condition u(x, t) = 0,

(x, t) ∈ ∂Ω × (0, T ),

(1.8)

where f and g are S-shaped flux functions of Buckley-Leverett type in the x− and y−direction respectively. The nonlinearity in the diffusion term, N, is a positive function. s(x, t) is a given source term, which is piecewise smooth. For simplicity, Ω is assumed to be the unit square, and the flux term is written in the following form ∇ · (f, g) =

∂ ∂ f (u) + g(u). ∂x ∂y

This identification process is based on the observed data k, that is u(xd , t) = k(xd , t),

d = 1, 2, . . . , D,

t ∈ (0, T ).

(1.9)

These observed data may contain noise. Eq. (1.6) closely resembles Eq. (1.3), except for the time derivative and convection terms. If porosity φ(x) = 1, the time derivative terms in Eqs. 3

(1.6) and (1.3) will be equal. The difference in the convection terms between Eqs. (1.6) and (1.3) is only that the convection term in Eq. (1.6) does not have permeability dependence and varying coefficient. Multiscale method was first studied to solve the forward problem, for example, the two-phase immiscible flow simulations in heterogeneous porous media [18,19], and then was extended to nonlinear problems [20,21]. Multiscale method has recently emerged in the field of inversion. This inversion strategy can speed up convergence, enhance stability of inversion, and avoid impact of local minimum. The main essence of multiscale inversion is to decompose the original problem into different scales or frequency bands, and to carry out inversion beginning at the smallest scale, because the objective function at a smaller scale has less local minimum and stronger convexity, so as to have a good opportunity of searching out the global minimum. The global minimum of the smallest scale is such a good approximation of the global minimum of the second smallest scale that it can be regarded as the initial estimate for the second optimization problem. Successively finding upward in this way can successfully search out the global minimum of the original problem. For the one-dimensional inversion problem, scale decomposition was shown to be very effective [22]. As a specific form of multiscale methods, multigrid method is used to reduce the computational requirements of large-scale numerical problems [23]. This method works by recursively moving between different scales, thereby propagating information between small and large scales. The multigrid method has been mainly used for solving partial differential equations [24], and more recently it has been applied to a variety of imaging problems such as optical diffusion tomography [25] and electrical impedance tomography [26]. Later, Oh et al. [27] constructed a general framework for nonlinear multigrid inversion, which can be employed for the solution of a large variety of inverse problems. The conclusion drawn from these works is that iterative inversion methods perform much better when there exists a scale decomposition. Wavelet has had many important applications in areas such as data compression, signal processing, and image analysis, however, its development in applied mathematics is fairly recent. The last few years have witnessed the intense activity and interest in the application of wavelet theory and its associated multiresolution analysis [28]. It is very interesting to use the wavelet multiresolution method to solve inverse problems. Recently, there exists a significant amount of research in the literature on the wavelet multiresolution 4

method for the parameter identification inverse problem. Successful applications of this method include the parameter identification inverse problem of an elliptical equation [29,30], the parameter identification inverse problem of a two-dimensional acoustic wave equation [31,32], the parameter identification inverse problem of a nonlinear diffusion equation [33], the parameter identification inverse problem of Maxwell equations [34], the parameter identification inverse problem of fluid-saturated porous media [35,36], the parameter identification inverse problem of the electrical capacitance tomography model [37] and the parameter identification inverse problem of the diffuse optical tomography model [38]. All these works showed the effectiveness of wavelet multiscale method on the parameter identification inverse problem. The aim of this paper is to develop a wavelet multiscale method for solving the permeability identification of a nonlinear convection-diffusion equation. The approach is based on the hierarchical approximation of wavelet, and combined with the regularized Gauss-Newton method, which is taken as the basic inversion method applied to each scale. This wavelet multiscale method has been proved to be widely convergent, computationally efficient, and have the anti-noise and de-noising abilities. This paper is arranged as follows: in Section 2, we give a brief introduction of wavelet and multiresolution analysis. In Section 3, we describe the discretized model and regularized Gauss-Newton method. In Section 4, we present the wavelet multiscale method in detail. In Section 5, we give the simulation examples. The conclusions are discussed in the final Section. 2. Wavelet and multiresolution analysis Definition 2.1. A function ψ ∈ L2 (R) is called a wavelet, if it satisfies the admissibility condition Z b 2 < ∞, Cψ = 2π dξ|ξ|−1|ψ(ξ)| (2.1) b where ψ(ξ) is the Fourier transform of ψ(x).

Definition 2.2. For any function f ∈ L2 (R), its continuous wavelet transform is defined by Z 1 x−b Wf (a, b) = dxf (x) · |a|− 2 ψ( ). (2.2) a 5

Its inverse transform is f (x) =

Cψ−1

Z

+∞ −∞

Z

+∞

−∞

dadb x−b Wf (a, b)ψ( ). 2 a a

(2.3)

Definition 2.3. A nest of closed subspaces {Vj , j ∈ Z} of L2 (R) is called a multiresolution analysis (MRA) for L2 (R), if it satisfies (1) Vj ⊂ Vj+1, ∀j ∈ Z; T Vj = {0}; (2) j∈Z S (3) ( Vj ) = L2 (R); j∈Z

(4) f (x) ∈ Vj ⇔ f (2x) ∈ Vj+1 , ∀j ∈ Z; (5) There exists a function φ ∈ V0 such that {φ(x − k), k ∈ Z} is an orthonormal basis for V0 . The dilates and translates of the scaling function φ comprise the basis function φj,k (x) in subspace Vj , and the dilates and translates of the wavelet function ψ comprise the basis function ψj,k (x) in subspace Wj , which is the orthogonal complement of Vj in Vj+1 . From Definition 2.3, it is easy to see VJ = VJ−1 ⊕ WJ−1 = VJ−2 ⊕ WJ−2 ⊕ WJ−1 = V0 ⊕ W0 ⊕ · · · ⊕ WJ−1 . (2.4) A given function can be expressed in terms of scaling and wavelet functions f=

X k

cJ,k φJ,k (x) =

X

c0,k φ0,k (x) +

J−1 X X j=0

k

dj,k ψj,k (x).

(2.5)

k

When the problem is decomposed into scale j from scale j +1, because the space Vj is only composed of the lower frequency part of the space Vj+1, the inversed solution in the space Vj should be an approximation to the original one, which does not contain the higher frequency part. Indeed, it is merely a local trend near the parameter of the model. Nevertheless, because of the higher frequency part contained in the subspace Wj of Vj+1, the information of higher frequency (which can’t be obtained at the smaller scale in the space Vj ) can be obtained through the inversion in the space Vj+1. Similarly, the information of further higher frequency, which can’t be obtained in the space Vj+1, can be obtained through the inversion in the space Vj+2 . Hence the resolution can be gradually enhanced. 6

3. Discretized model and regularized Gauss-Newton method 3.1. Finite difference discretization By using the classic finite difference scheme, Eqs. (1.6)-(1.8) become n−1 uni,j − ui,j n−1 n−1 n + ∇ · (f (ui,j ), g(ui,j )) − ∇ · (qi,j Ni,j ∇uni,j ) = sni,j , ∆t n = 1, 2, . . . , M; i = 1, 2, . . . , n1 − 1; j = 1, 2, . . . , n2 − 1, u0i,j = u0 (i∆x, j∆y), i = 0, 1, . . . , n1 ; j = 0, 1, . . . , n2 ,

(3.1)

un0,j = unn1 ,j = uni,0 = uni,n2 = 0, n = 1, 2, . . . , M;

i = 0, 1, . . . , n1 ;

j = 0, 1, . . . , n2 ,

n where uni,j = u(i∆x, j∆y, n∆t), qi,j = q(i∆x, j∆y), Ni,j = N(uni,j ), sni,j = s(i∆x, j∆y, n∆t), n1 = 1/∆x, n2 = 1/∆y, M = T /∆t. ∆x and ∆y are the step sizes of the rectangular grid in the x− and y−direction, respectively, n−1 n−1 n and ∆t is the time step size. ∇ · (qi,j Ni,j ∇uni,j ) and ∇ · (f (ui,j ), g(ui,j )) are respectively the discretizations of the diffusion term and the convection term, and to obtain their expressions, the following definitions are given. Let uni,j − uni−1,j uni+1,j − uni,j x n x n D− ui,j = , D+ ui,j = ∆x ∆x denote the discrete derivatives in the x−direction for the backward and forward difference approximations, respectively. For the the discrete derivatives in the y−direction, a corresponding notation is used. The mean values concerning with the discretized permeability q and nonlinear function N(u) are defined as

1 1 qxi,j+ 1 = (qi+ 1 ,j+ 1 + qi− 1 ,j+ 1 ), q yi+ 1 ,j = (qi+ 1 ,j+ 1 + qi+ 1 ,j− 1 ), 2 2 2 2 2 2 2 2 2 2 2 2 1 n 1 n x y n n + Ni,j ), (N )ni,j+ 1 = (Ni,j+1 + Ni,j ), (N )ni+ 1 ,j = (Ni+1,j 2 2 2 2 where qi+ 1 ,j+ 1 = q((i + 21 )∆x, (j + 12 )∆y). Then, the discretization of the 2 2 diffusion term is expressed as x

y

y y n x n n x ui,j ) + D− (qxi,j+ 1 (N )ni,j+ 1 D+ ui,j ). ∇ · (qi,j Ni,j ∇uni,j ) = D− (q yi+ 1 ,j (N )ni+ 1 ,j D+ 2

2

2

2

For the convection term, the Engquist-Osher upwind scheme is used (see [39]) y EO n−1 n−1 n−1 x EO n−1 n−1 n−1 ∇ · (f (ui,j ), g(ui,j )) = D− f (ui,j , ui+1,j ) + D− g (ui,j , ui,j+1 ),

7

where n−1 n−1 (ui,j , ui+1,j )

1 1 n−1 n−1 = (f (ui,j ) + f (ui+1,j )) − 2 2

n−1 n−1 (ui,j , ui,j+1 )

1 1 n−1 n−1 = (g(ui,j ) + g(ui,j+1 )) − 2 2

f

EO

g

EO

Z

Z

un−1 i+1,j un−1 i,j un−1 i,j+1

un−1 i,j

|f ′ (ξ)|dξ, |g ′(ξ)|dξ.

The explicit formulas for f EO and g EO , for examples of f and g, are calculated in the Appendix A. 3.2. Fix point iteration Due to the implicit treatment of the diffusion term, the system that needs to be solved for each time level in the forward simulation is nonlinear. A fixpoint iteration is used to solve these resultant nonlinear systems. For the fixed time level, the nonlinear system is written on the form B(u)u = b, where B(u) is a nonsingular matrix depending nonlinearly on u, then we iterate B(uk−1)uk = b, until convergence, with u0 given. u0 is usually the solution from the previous time level. Let u⊤ = (un1,1, un1,2 , . . . , un1,n2−1 , un2,1 , un2,2 , . . . , un2,n2 −1 , . . . , unn1 −1,1 , unn1 −1,2 , . . . , unn1 −1,n2 −1 ),

then B is a block tridiagonal matrix. This block tridiagonal matrix equations can be solved conveniently with the chase method. This fix-point iteration method is only used to solve the forward problem (3.1). 3.3. Coarse grid for q From the above-mentioned method, it can be seen that the discretized permeability used in the numerical solution method of the forward problem is qi+ 1 ,j+ 1 , i = 0, 1, . . . , n1 − 1; j = 0, 1, . . . , n2 − 1. So the discretized 2 2 permeability vector q ∈ Rn1 ×n2 is needed. The initial and coarse grids are respectively denoted as Ωn1 ×n2 and Ωm1 ×m2 , where Ωm1 ×m2 is got by coarsening Ωn1 ×n2 , n1 × n2 and m1 × m2 are respectively the grid numbers of the initial and coarse grids. Similar to [40], 8

q ∈ Rm1 ×m2 is defined at Ωm1 ×m2 , and a constant prolongation is used when q ∈ Rn1 ×n2 is needed. The constant prolongation means that the entry of q at Ωn1 ×n2 is obtained by taking the value from the closest point of Ωm1 ×m2 . Actually, m1 × m2 is generally dependent on the number of available measurements. 3.4. Discretized model Eq. (1.9) can be discretized as unxd = kxnd ,

n = 1, 2, . . . , M;

d = 1, 2, . . . , D.

(3.2)

Eqs. (3.1) and (3.2) define a nonlinear vector-valued function F : Q → K, where Q and K denote vectors, which are respectively the compositions of qi,j and kxnd in the proper sequence. Let b kxnd denote the observations and form b in the same sequence as K. Especially, let the vector K Q⊤ = (q1,1 , q1,2 , . . . , q1,m2 , q2,1 , q2,2 , . . . , q2,m2 , . . . , qm1 ,1 , qm1 ,2 , . . . , qm1 ,m2 ),

K ⊤ = (kx1 1 , kx1 2 , . . . , kx1 D , kx2 1 , kx2 2 , . . . , kx2 D , . . . , kxM1 , kxM2 , . . . , kxMD ), b ⊤ = (b kxMD ). kxM2 , . . . , b kx2 D , . . . , b kxM1 , b kx2 2 , . . . , b kx1 D , b kx2 1 , b kx1 2 , . . . , b K kx1 1 , b

Then the inversion for Q is reduced to a nonlinear operator equation b F (Q) = K.

(3.3)

The observation points xd , d = 1, 2, . . . , D, are located at grid points in our numerical simulations. Nevertheless, in case the observation points are not on the grid points, the values of kxnd will be obtained using a spatial linear interpolation of un . 3.5. Regularized Gauss-Newton method As is well known, Eq. (3.3) is an ill-posed problem, so the Tikhonov regularization needs to be introduced. Instead of directly solving Eq. (3.3), the following optimal problem is considered n o 2 2 2 b min kF (Q) − Kk + β1 kL1 Qk + β2 kL2 Qk , (3.4)

where β1 , β2 are the regularization parameters, and L1 , L2 are the secondorder smooth matrices in the x− and y−direction, respectively (see [31]). 9

The regularized Gauss-Newton method can be used to find the minimum ⊤ −1 Qk+1 = Qk − [F ′ (Qk )⊤ F ′ (Qk ) + β1 L⊤ 1 L1 + β2 L2 L2 ] k ⊤ k b + β1 L⊤ × [F ′ (Qk )⊤ (F (Qk ) − K) 1 L1 Q + β2 L2 L2 Q ],

(3.5)

k = 0, 1, 2, . . .

It has a fast convergence speed and good stability. But like most iterative methods, it is only of local convergence. 4. Wavelet multiscale method Theoretically, the inverse problem of nonlinear convection-diffusion equation can effectively be solved by the general iterative methods, such as the regularized Gauss-Newton method, with a good initial estimate given. However, the direct application of such inversion methods to the real observed data has been disappointing when the initial estimate is far away from the true permeability or there is no prior information about permeability available. The reason for this is that numerous local minima in the objective function prevent convergence. As an emerging inversion strategy, wavelet multiscale method has explored the possibility of diminishing the problem of local minima by a decomposition of inverse problem by scale. When the inverse problem is decomposed into different scales from the largest to the smallest, the small scale components contain little varying features because at small scales the number of local minima is enormously reduced and those that remain are further apart from each other. Therefore, the global minimum at a small scale can more easily be found out than at a large scale by the general iterative methods. The global minimum at the smallest scale is such a good approximation of the global minimum at the second smallest scale that it can be regarded as a good initial estimate for the second smallest scale component. Then, the second smallest scale component can easily be solved by the general iterative methods with a good initial estimate given. Successively fining upward in this way can succeed in finding out the global minimum of the original inverse problem. As shown in Figure 1, there are three elements in wavelet multiscale method: restriction operator RO(J, J − 1), relaxation operator R(J − 1) and injection operator IO(J − 1, J). In Figure 1, S(J) represents the largest scale J on which the original inverse problem is defined. The restriction 10

operator RO(J, J − 1), which restricts the original inverse problem from S(J) to S(J − 1), is chosen as the wavelet decomposition algorithm RO(J, J − 1) : S(J) → S(J − 1). Once the problem has been restricted to the smaller scale S(J − 1), the global minimum at this scale can be found out by the relaxation operator R(J − 1). This paper uses the regularized Gauss-Newton method as the relaxation operator. Then, the global minimum at scale S(J − 1) is injected back up to the scale S(J) with the injection operator IO(J − 1, J), which is chosen as the wavelet reconstruction algorithm IO(J − 1, J) : S(J − 1) → S(J). The parts of the original inverse problem, which were not solved at scale S(J − 1), can be solved by the relaxation operator R(J). This completes the description of the two-scale wavelet multiscale solution for the inverse problem. The multi-scale wavelet multiscale method can be easily obtained by replacing the relaxation operator R(J − 1) with another two-scale wavelet multiscale solution recursively. R (J )

S (J ) RO(J ,J S (J

1)

-

IO J (

1)

-

R(J

J

)

1)

-

R ( 1)

S (1)

S (0 )

- 1,

I O(0, 1)

R O(1 ) ,0 R (0 )

Figure 1: Description of the wavelet multiscale method.

The object to be decomposed may be the source, the wave field, the permeability model, or a combination of the three. There are two basic methods to carry out the decomposition. The first method is to make wavelet transform on all the source, the wave field and the permeability model. The 11

increase in the scale leads to the great reduction in the dimension of the discretized inverse problem, especially for high-dimensional continuous problem. But the difficulties that confront us are that the treatment on the boundary and the direct wave is troublesome, and the high precision of forward computation is demanded. The second method is to make wavelet transformation only on the source and the wave field. This paper uses the second method, because it is very simple, and does not need the high precision of forward computation. The algorithm is summarized as follows: Algorithm 4.1. The wavelet multiscale method. (1) For two-dimensional nonlinear convection-diffusion equation, give the true permeability model Q∗ , the source function s, an initial permeability Q0 , b and the observed data K. b into different scales S(j), j = J, J −1, . . . , 0, using (2) Decompose s and K the compactly supported orthogonal wavelet “DB4” and Mallat algorithm [41], b j, j = namely the original problem (3.3) can be decomposed as F j (Q) = K ∗ 0 J, J − 1, . . . , 0. Set Q0 = Q , k = 0. b k with the (3) At scale S(k), calculate the solution Q∗k+1 of F k (P ) = K initial permeability Q∗k using the regularized Gauss-Newton method. (4) Let k = k + 1 and check if k > J, then the solution Q∗J+1 obtained at scale S(J) is the solution of original problem (3.3), else if k ≤ J, then go to step 3. 5. Numerical examples In the following numerical examples, we set Ω = [0, 1] × [0, 1],

u0 (x) = 0,

T = 0.06,

M = 30.

The source function is s(x, t) = g(t)

16 X i=1

i

δ(x − x ) − 4g(t − r)

4 X p=1

δ(x − xp ),

where g(t) = 1.6 sin(80πt) exp(−80π(t/2)), r ≈ 0.002, δ(·) is the Dirac δfunction, and xi , xp are the injector points and producer points, respectively. Here g(t) is chosen as an oscillatory function, which can effectively validate 12

the denoising ability of the wavelet multiscale method. In fact, a slowly decaying function or even a constant would be more realistic for the considered problem. For these functions, this method can also work, but the denoising ability is not obvious. Figure 2 shows the exact locations of injector points and producer points in simulation geometry, and we place 16 receivers in each injector point simultaneously, which means D = 16. In simulation, the grid cell size of Ωm1 ×m2 is (2∆x, 2∆y), so m1 = n21 , m2 = n22 .

Figure 2: Simulation geometry with the locations of injector points and producer points.

In the first three examples, ∆x = ∆y = 81 and the true permeability, q(x), is the piecewise constant  c q, x ∈ [0, 0.5] × [0, 0.5],    1c q2 , x ∈ [0, 0.5] × [0.5, 1], q(x) = c q , x ∈ [0.5, 1] × [0, 0.5],    3c q4 , x ∈ [0.5, 1] × [0.5, 1],

The results seem to be best when the regularization parameters are chosen as β1 = β2 = 0. As pointed out in [40], this is because the dimension is relatively 1 small in the first three examples. In the last three examples, ∆x = ∆y = 16 , −12 β1 = β2 = 10 and the complex permeability models will specifically be described in the each example. The flux function in the x−direction is an S-shaped Buckley-Leverett flux with gravitational effects f (u) =

u2 (1 − 5(1 − u2 )) , u2 + (1 − u)2 13

(5.1)

and the flux function in the y−direction is an S-shaped Buckley-Leverett flux g(u) =

u2 . u2 + (1 − u)2

(5.2)

The largest scale J of wavelet decomposition is 3 in all examples except for Example 5.1. The regularized Gauss-Newton method (3.5) is terminated when the modification is equal to or less than 0.01. Example 5.1. In the first example, the ground truth is qic = i, i = 1, . . . , 4, and we respectively use the wavelet multiscale method with J = 6 (WM6), wavelet multiscale method with J = 3 (WM3), and regularized GaussNewton method (RGN). The nonlinear diffusion function is N(u) = u2 +u+1, and the initial estimate Q0 is given as a constant. To compare differences among these methods, the iteration numbers when the parameters converge to the real value are compiled in Table 1. Table 1 The required iteration numbers by WM6, WM3 and RGN in Example 5.1. Number Initial value WM6 WM3 RGN 1 2.05 8 5 8 2 2.15 8 6 9 3 2.25 8 7 No convergence 4 2.35 8 8 No convergence 5 2.45 8 10 No convergence 6 2.55 10 No convergence No convergence We first compare WM3 with RGN. Looking at Table 1, from experiments 1 and 2, it is observed that WM3 is more efficient than RGN. Furthermore, from experiments 3-5, it can be seen that when the initial estimate is poor, there is no convergence by RGN, while for WM3, the convergence can be achieved. Then, WM6 is compared with WM3 and RGN. It can be seen from experiment 6 that as the scale of the wavelet multiscale method is increased, its convergence region gets wider. In addition, experiments 1-5 show that when the initial estimate is around the real value, the increase in the scale leads to the increase in the calculation amount; however, when the initial estimate is far away from the real value, the increase in the scale can lead to the decrease in the calculation amount. As a whole, it is concluded that the wavelet multiscale method has a wide convergence region and a fast convergence speed. 14

Example 5.2. In this example, q1c = q4c = 2.5, q2c = q3c = 4.5, and we use Algorithm 1. The nonlinear diffusion function is N(u) = u4 + u3 + u2 + u + 1, Q0 ≡ 3, and 5% Gaussian noise is added to the observed data. The true model and inversion result are shown in Figure 3. Example 5.3. In oil reservoirs the permeability usually has greatly large jumps, so we use Algorithm 1 with q1c = 8.5, q2c = 1.5, q3c = 2.5, q4c = 9.5 in this example. The nonlinear diffusion function is N(u) = u4 − u3 + u2 − u + 1, Q0 ≡ 3, and 5% Gaussian noise is added to the observed data. The true model and inversion result are shown in Figure 4. It can be seen from these two examples that the wavelet multiscale method is effective for the different N−functions, even if N is highly nonlinear. Moreover, in Example 5.3, the regularized Gauss-Newton method is not convergent when the initial estimate is chosen as an arbitrary constant. Hence, we conclude that the wavelet multiscale method is effective and even necessary in oil reservoir simulation.

4.5

5 4.5

4 4 3.5

3.5 3

3 2.5 2.5 8

2 8 6

8

6

6

4 2

4

2

2 0

8 6

4

4

2 0

0

(a) The true model.

0

(b) The inversion result.

Figure 3: The true model and inversion result with 5% Gaussian noise added in Example 5.2.

Example 5.4. In this example we consider a level-stratified medium containing two interfaces, which is depicted in Figure 5a. The nonlinear diffusion function is N(u) = u4 + u3 + u2 + u + 1, Q0 ≡ 3, and 1% Gaussian noise is added to the observed data. Figure 5b shows the inversion result obtained by Algorithm 1. Starting from the same initial estimate, RGN is not convergent. 15

10

10

8

8

6

6

4

4

2

2

0 8

0 8 6

6

8 6

4 2

4

2

2 0

8 6

4

4

2 0

0

(a) The true model.

0

(b) The inversion result.

Figure 4: The true model and inversion result with 5% Gaussian noise added in Example 5.3.

Example 5.5. In this example we consider the model of one anomalous body in a homogeneous medium with a permeability of 6.28, and the anomalous body has the permeability of 2.86. The nonlinear diffusion function is N(u) = u4 − u3 + u2 − u + 1, Q0 ≡ 3, and 1% Gaussian noise is added to the observed data. Figure 6 shows this model and the inversion result obtained by Algorithm 1. Starting from the same initial estimate, RGN is also convergent, but takes about 131.85% CPU time of the wavelet multiscale method. Example 5.6. In this example we consider the model of two anomalous bodies in a homogeneous medium with a permeability of 5.33. The anomalous bodies have the permeability of 2.25 and 7.01, respectively. The nonlinear diffusion function is N(u) = u2 + u + 1, Q0 ≡ 3, and 2% Gaussian noise is added to the observed data. This model and the inversion results at scale 2 and scale 3 obtained by Algorithm 1 are respectively shown in Figures 7 and 8. To compare the inversion results at different scales, the relative l2 errors vs. the scale numbers are plotted in Figure 9. It can be seen from these three examples that even though the model is large-scale and complicated, the wavelet multiscale method is still widely convergent, computationally efficient, and has the anti-noise ability. Moreover, we see from Figure 9 that in Example 5.6 the maximum relative l2 error occurs at scale 0, and the minimum relative l2 error occurs at scale 2, rather than scale 3. So it is easily concluded that the wavelet multiscale method can eliminate noise to a certain degree when the data is strongly noisy. The 16

reason for this is probably that wavelet has good de-noising performance.

5.5

5 2

2

5

4.5 4 4

6 8

4

4.5

6

4 3.5

8

3.5

3 10

10 3

2.5

12

12 2

2.5 14

14 1.5 2

16 2

4

6

8

10

12

14

16

16

1 2

(a) The true model.

4

6

8

10

12

14

16

(b) The inversion result.

Figure 5: The true model and inversion result with 1% Gaussian noise added in Example 5.4.

6. Conclusions This paper has studied the wavelet multiscale method for recovering the diffusion parameter in a nonlinear convection-diffusion equation, which arises as the saturation equation in the fractional flow formulation of the twophase porous media flow equations. The numerical simulations show that the proposed method is widely convergent, computationally efficient, and has the anti-noise and de-noising abilities. Furthermore, it is effective for the different N-functions, even if N is highly nonlinear, and also necessary in oil reservoir simulation due to the greatly large jumps in the permeability. It is interesting that the proposed method is without the loss of generality, and can be generalized to the other kinds of inverse problems. Acknowledgements The author is grateful to anonymous referees for their valuable comments. The author would like to thank Michael Ng and the typesetting team for their kind help. This work was supported by the National Natural Science Foundation of China (11271102), the Science and Technology Funds for the

17

7 2

6

2

4

5.5

4

6.5

6

6 5.5

6

5

8

5

8 4.5

4.5

10

10 4

4 12

12 3.5 3.5

14

14

3

16 2

4

6

8

10

12

14

3 2.5

16

16

2

(a) The true model.

4

6

8

10

12

14

16

(b) The inversion result.

Figure 6: The true model and inversion result with 1% Gaussian noise added in Example 5.5.

7 2

6.5

4

6 5.5

6

5 8 4.5 10 4 12

3.5 3

14

2.5

16 2

4

6

8

10

12

14

16

Figure 7: The true model in Example 5.6.

18

6.5

2

6.5

2

6

6 4

4 5.5

5.5 6

6 5

5 8 10 12 14

8

4.5 4

4

3.5

12

3.5

3

14

2

4

6

8

10

12

14

3 2.5

2.5

16

4.5

10

16

16

2

(a) The inversion result at scale 2

4

6

8

10

12

14

16

(b) The inversion result at scale 3

Figure 8: The inversion results with 2% Gaussian noise added in Example 5.6.

0.13 0.12

The relative l 2 error

0.11 0.1 0.09 0.08 0.07 0.06 0.05

0

0.5

1

1.5 2 The scale number

2.5

3

Figure 9: The relative l2 errors with respect to the scale numbers in Example 5.6.

19

Universities in He-bei Province of China (Z2015039), the Fundamental Research Funds for the Central Universities (N142303011), the PhD Foundation of Northeast University at Qinhuangdao (XNB2015002). Appendix A. The numerical flux function The explicit formula for the Engquist-Osher numerical flux functions f EO and g EO are given in this appendix. For simplicity, we only use one subscript index in the following calculations. Appendix A.1. Buckley-Leverett flux function with gravitational effects The flux function in the x−direction is an S-shaped Buckley-Leverett flux with gravitational effects, see Eq. (5.1). In order to calculate R ui+1function ′ |f (ξ)|dξ, the sign of f ′ in (0, 1) is studied first. f ′ has only one zero ui point in (0, 1), and this zero point is given by Thus, we see that

xzero ≈ 0.37.

 ξ ∈ (0, xzero),  < 0, ξ = xzero , f ′ (ξ) = 0,  > 0, ξ ∈ (xzero,1 ). In addition, we can calculate the integral  f (ui+1 ) − f (ui ), ui , ui+1 ≥ xzero ,   Z ui+1  f (u ) − f (u ), u i i+1 i , ui+1 < xzero , |f ′ (ξ)|dξ = f (ui ) + f (ui+1) − 2f (xzero), ui < xzero, ui+1 ≥ xzero ,  ui   2f (xzero ) − f (ui ) − f (ui+1), ui ≥ xzero , ui+1 < xzero ,

and by the definition of f EO , we then get  f (ui ),    f (ui+1), f EO (ui , ui+1) = f (xzero),    f (ui ) + f (ui+1 ) − f (xzero ),

ui , ui+1 ≥ xzero, ui , ui+1 < xzero , ui < xzero , ui+1 ≥ xzero, ui ≥ xzero , ui+1 < xzero.

Appendix A.2. Buckley-Leverett flux function The flux function in the y−direction is an S-shaped Buckley-Leverett flux function, see Eq. (5.2). Noticing that it follows

g ′ (u) ≥ 0,

u ∈ (0, 1),

g EO (ui , ui+1 ) = g(ui ). 20

References [1] M.S. Espedal, K.H. Karlsen, Numerical solution of reservoir flow models based on large time step operator splitting algorithm, in: A. Fasano (Eds.), Filtration in Porous Media and Industrial Applications: Lecture Notes in Mathematics, Springer, Berlin, 2000, pp. 9-77. [2] M.S. Gockenbach, A.A. Khan, An abstract framework for elliptic inverse problems: part 1. An output least-squares approach, Math. Mech. Solids. 12(2007)259-276. [3] M.S. Gockenbach, A.A. Khan, An abstract framework for elliptic inverse problems: part 2. An augmented Lagrangian approach, Math. Mech. Solids. 14(2009)517-539. [4] I. Knowles, Parameter identification for elliptic problems, J. Comput. Appl. Math. 131(2001)175-194. [5] Z. Chen, J. Zou, An augmented Lagrangian method for identifying discontinuous parameters in elliptic systems, SIAM J. Control Optim. 37(1999)892-910. [6] M. Hayek, P. Ackerer, E. Sonnendr¨ ucker, A new refinement indicator for adaptive parameterization: application to the estimation of the diffusion coefficient in an elliptic problem, J. Computat. Appl. Math. 224(2009)307-319. [7] T.F. Chan, X.C. Tai, Level set and total variation regularization for elliptic inverse problems with discontinuous coefficients, J. Comput. Phys. 193(2003)40-66. [8] K. Kunisch, G. Peichl, Estimation of a temporally and spatially varying diffusion coefficient in a parabolic system by an augmented Lagrangian technique, Numer. Math. 59(1991)473-509. [9] H.W. Engl, J. Zou, A new approach to convergence rate analysis of Tikhonov regularization for parameter identification in heat conduction, Inverse Probl. 16(2000)1907-1923. [10] T.K. Nilssen, X.C. Tai, Parameter estimation with the augmented Lagrangian method for a parabolic equation, J. Optim. Theory Appl. 124(2005)435-453. 21

[11] B. Gou, J. Zou, An augmented Lagrangian method for parameter identification in parabolic systems, J. Math. Anal. Appl. 263(2001)49-68. [12] Y.L. Keung, J. Zou, Numerical identifications of parameters in parabolic equations, Inverse Probl. 14(1998)83-100. [13] L. Li, W.Y. Liu, B. Han, Dynamical level set method for parameter identification of nonlinear parabolic distributed parameter systems, Commun. Nonlinear Sci. Numer. Simulat. 17(2012)2752-2765. [14] P.C. Shah, G.R. Gavalas, J.H. Seinfeld, Error analysis in history matching, the optimum level of parameterization, Soc. Pet. Eng. J. 6(1978)219228. [15] A.C. Reynolds, N. He, L. Chu, D. Oliver, Reparameterization techniques for generating reservoir descriptions conditioned to variograms and wellCtest pressure data, SPE 30588, in: Proceedings of the SPE Annual Technical Conference and Exhibition, Dallas, USA, 1995. [16] A.A. Grimstad, T. Mannseth, J.E. Nordtvedt, G. Nævdal, Reservoir characterization through scale splitting, in: Proceedings of the 7th European Conference on the Mathematics of Oil Recovery, Baveno, Italy, 2000. [17] M. Cuypers, O. Dubrule, P. Lamy, R. Bissell, Optimal choice of inversion parameters for history matching with the pilot point method, in: Proceedings of the 6th European Conference on the Mathematics of Oil Recovery, Peebles, UK, 1998. [18] Y. Efendiev, V. Ginting, T.Y. Hou, R. Ewing, Accurate multiscale finite element methods for two-phase flow simulations, J. Comput. Phys. 220(2006)155-174. [19] Y. Efendiev, T.Y. Hou, Multiscale finite element methods for porous media flows and their applications, Appl. Numer. Math. 57(2007)577596. [20] Y. Efendiev, T.Y. Hou, V. Ginting, Multiscale finite element methods for nonlinear problems and their applications, Commun. Math. Sci. 2(2004)553-589. 22

[21] Y. Efendiev, T.Y. Hou, Multiscale Finite Element Methods: Theory and Applications, Springer Science & Business Media, 2009. [22] P. Kolb, F. Collino, P. Lailly, Pre-stack inversion of a 1-D medium, Proc. IEEE 74(1986)498-508. [23] W. Hackbusch, Multigrid Methods and Applications, Spinger, Berlin, 1985. [24] S.R. Fulton, P.E. Ciesielski, W.H. Schubert, Multigrid methods for elliptic problems: A review, J. Atmos. Sci. 114(1986)943-959. [25] J.C. Ye, C.A. Bouman, K.J. Webb, R.P. Millane, Nonlinear multigrid algorithms for Bayesian optical diffusion tomography, IEEE T. Image Process. 10(2001)909-922. [26] L. Borcea, A nonlinear multigrid for imaging electrical conductivity and permittivity at low frequency, Inverse Probl. 17(2001)329-359. [27] S. Oh, A.B. Milstein, C.A. Bouman, K.J. Webb, A general framework for nonlinear multigrid inversion, IEEE T. Image Process. 14(2005)125-140. [28] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, 1992. [29] J. Liu, A multiresolution method for distributed parameter estimation, SIAM J. Sci. Comput. 14(1993)389-405. [30] H.S. Fu, B. Han, H.B. Liu, A wavelet multiscale iterative regularization method for the parameter estimation problems of partial differential equations, Neurocomputing 104(2013)138-145. [31] H.S. Fu, B. Han, A wavelet multiscale method for the inverse problems of a two-dimensional wave equation, Inverse Probl. Sci. Eng. 12(2004)643656. [32] H.S. Fu, B. Han, G.Q. Gai, A wavelet multiscale-homotopy method for the inverse problem of two-dimensional acoustic wave equation, Appl. Math. Comput. 190(2007)576-582. [33] J. Zhao, T. Liu, S. Liu, Identification of space-dependent permeability in nonlinear diffusion equation from interior measurements using wavelet multiscale method, Inverse Probl. Sci. Eng. 22(2014)507-529. 23

[34] L. Ding, B. Han, J.Q. Liu, A wavelet multiscale method for inversion of Maxwell equations, Appl. Math. Mech. 30(2009)1035-1044. [35] X.M. Zhang, K.A. Liu, J.Q. Liu, The wavelet multiscale method for inversion of porosity in the fluid-saturated porous media, Appl. Math. Comput. 180(2006)419-427. [36] Y. He, B. Han, A wavelet adaptive-homotopy method for inverse problem in the fluid-saturated porous media, Appl. Math. Comput. 208(2009)189-196. [37] J. Lei, S. Liu, Z. Li, M. Sun, X. Wang, A multi-scale image reconstruction algorithm for electrical capacitance tomography, Appl. Math. Model. 35(2011)2585-2606. [38] F. Dubot, Y. Favennec, B. Rousseau, D.R. Rousse, A wavelet multiscale method for the inverse problem of diffuse optical tomography, J. Comput. Appl. Math. 289(2015)267-281. [39] B. Engquist, S. Osher, One-sided difference approximations for nonlinear conservation laws, Math. Comput. 36(1984)321-351. [40] T.K. Nilssen, K.H. Karlsen, T. Mannseth, X.C. Tai, Identification of diffusion parameters in a nonlinear convection-diffusion equation using the augmented Lagrangian method, Comput. Geosci. 13(2009)317-329. [41] S.G. Mallat, A theory for multiresolution signal decomposition: the wavelet representation, IEEE Trans. Pattern Anal. Mach. Intell. 11(1989)674-693.

24