Numerical Method for a Filtration Model Involving a Nonlinear Partial Integro-Differential Equation

In this paper, we propose an efficient numerical method for solving an initial boundary value problem for a coupled system of equations consisting of a nonlinear parabolic partial integrodifferential equation and an elliptic equation with a nonlinear term. This problem has an important applied significance in petroleum engineering and finds application in modeling two-phase nonequilibrium fluid flows in a porous medium with a generalized nonequilibrium law. The construction of the numerical method is based on employing the finite element method in the spatial direction and the finite difference approximation to the time derivative. Newton’s method and the secondorder approximation formula are applied for the treatment of nonlinear terms. The stability and convergence of the discrete scheme as well as the convergence of the iterative process is rigorously proven. Numerical tests are conducted to confirm the theoretical analysis. The constructed method is applied to study the two-phase nonequilibrium flow of an incompressible fluid in a porous medium. In addition, we present two examples of models allowing for prediction of the behavior of a fluid flow in a porous medium that are reduced to solving the nonlinear integro-differential equations studied in the paper.


Introduction
Fluid flows through porous media are of fundamental significance to a wide range of natural and industrial processes in groundwater hydrology, enhanced oil recovery, climate engineering and many others. These processes are generally described by a system of partial differential equations. However, in certain circumstances, these equations may not provide an adequate and accurate representation of the fluid flow process due to neglecting important properties of the porous medium [1][2][3].
For example, a number of authors have conceived that the dynamics of fluid flow in a porous medium are influenced not only by the current state of the process but also by its previous states [4,5]. In particular, a number of studies have shown significant discrepancies between experimental results on wells and simulation results based on classical models [1,6]. These hereditary properties of the process, known as memory, can appear when modeling non-local transfer of reaction and contaminants, thermal conductivity [4,[7][8][9][10], and many others.
The memory effects can be taken into account by the inclusion of a non-trivial memory term in the model equation, which is often represented as a convolution of the function or its derivative with some kernel [11]. Thus, a fluid flow model in a porous media taking into account memory effects can be described by partial integro-differential equations.
Another example of the application of nonlinear partial integro-differential equations are models describing multi-phase nonequilibrium fluid flow in a porous medium. One of these models was proposed in [12,13], which generalized the classical nonequilibrium model of G. Barenblatt [14]. The idea behind this generalization lies in the assumption that all flows are nonequilibrium and differ only in the degree of nonequilibrium. In this regard, the authors of [12,13] introduced a function for each flow varying between 0 and 1 that shows to what extent the flow is non-equilibrium. Moreover, the limit values, 0 and 1, correspond to completely equilibrium and completely non-equilibrium cases.
The nonlinear integro-differential equation in the model studied in [12] describes the distribution of effective water saturation. The consideration of nonequilibrium effects appears to be necessary at all stages in the development of oil fields since the pressure drop versus time obtained during laboratory studies of porous media in order to determine the relative phase permeability functions significantly differ from the theoretical curves calculated in the framework of the classical filtration theory [15]. The effect of nonequilibrium can be significant: the time of saturation establishment in the conditions of oil fields can reach the order of a year [14].
The above two examples clearly demonstrate that nonlinear integro-differential equations make it possible to describe extremely important processes occurring in porous media.
Nonlinear integro-differential equations are known to admit an analytical solution only in a few particular cases due to their complexity. Over the last several decades, many techniques have been developed for the approximate solution of these equations, which include, for example, the combination of the radial basis functions and finite difference method [16], space-time spectral collocation method [17], Galerkin method with a nonclassical H 1 -projection [18].
Considering the importance of the subject and motivated by the model proposed in [12,13], we devote our study to an initial boundary problem for a nonlinear system of three coupled equations. The first equation of the system is a non-linear convectiondiffusion-reaction equation with an additional time-dependent integral term, which, in [12], describes the dynamics of the effective water saturation through a porous medium. The second and third equations together constitute an elliptic equation containing a nonlinear term and describe the dynamics of the pressure field.
The problem under consideration has been numerically solved in a few works. The paper [12] introduced an explicit difference scheme to numerically implement the integrodifferential equation, whereas the elliptic equation was solved using the successive overrelaxation iterative method. In contrast, in [39], an implicit finite difference scheme was utilized to solve the partial integro-differential equation, which was then implemented by the alternate directions method. The elliptic equation was solved using the iterative procedure of establishement method. The stability and convergence analysis of the constructed method has been discussed.
Although the numerical schemes constructed in [12,39] made it possible to assess the impact of nonequilibrium effects on the fluid flow within the framework of the constructed model, in our opinion, insufficient attention was paid to the handling of nonlinear terms which were treated as known functions, and their values were taken from previous time layers. This approach is known to lower the convergence order with respect to the time variable.
Unfortunately, the literature review did not reveal any other works aimed at solving the filtration problem under study. However there are a few alternative efficient approaches of solving nonlinear partial integro-differential equations similar to the first equation of the system under consideration. For example, the papers [4,8,40] have utilized a two-grid finite element or block-centered finite difference method to the nonlinear non-Fickian flow model such that a rough approximation was obtained on a coarse grid, and then the corresponding linearized problem was solved on a fine grid.
In [41], this equation was reduced to a system of ordinary differential equations by virtue of a combination of group preserving scheme and spectral meshless radial point interpolation, which was then solved by the method of lines. There are also known implementations of the non-linear non-Fickian flow model by the finite volume element method [42], mixed finite element method [43], expanded H 1 -Galerkin mixed finite element method [44].
The authors of [45] proposed a numerical method based on the matrix collocation method descended from the use of Laguerre and Taylor polynomials. In [46], the Kirchhoff transformation of the dependent variable was applied, and then the solution was obtained using Galerkin approximations.
Our work uses a different approach. We constructed a numerical method based on the use of a finite difference approximation in time and a finite element approximation in the spatial direction. To achieve a global second convergence order in time, we employ different approximations of the time derivative for the first and subsequent time layers. The integral term is approximated using the Lagrange interpolation formula. For treatment of one nonlinear term, Newton's method is applied, and the remaining nonlinear terms are linearized using a second-order difference approximation.
Secondly, we rigorously prove the stability and convergence of the fully discrete scheme as well as the convergence of Newton's iterative process. Further, the theoretical convergence order is confirmed by numerical tests conducted for a problem with a known exact solution. We also present a comparative analysis with the results obtained by the finite difference method proposed in [12].
Thirdly, we apply the constructed numerical method to modeling a two-phase nonequilibrium flow of an incompressible fluid in a porous medium according to the model proposed in [12]. The numerical scheme has shown itself to be stable for a wide range of time step values.
Fourthly, as a supplement, we consider two examples of fluid flow problems in porous media, which are reduced to solving the nonlinear integro-differential equation considered in the paper. The first problem is concerned with a special case of nonequilibrium filtration referred to as the countercurrent filtration problem. The second problem deals with the filtration process in a fractured porous medium.
The main contribution of our paper is an approach that make it possible to effectively account for nonlinear terms, thus, maintaining a balance between computational complexity and the achievement of the global second-order convergence.
The paper is structured as follows. The next section presents the problem statement and defines the assumptions under which the problem is solved. Section 2.2 is devoted to the construction of the fully discrete scheme and Newton's iterative method. In Sections 2.3 and 2.4, the stability and convergence of the fully discrete schemes are discussed. Section 2.5 determines the conditions under which the quadratic convergence of Newton's method is achieved. The remainder of the paper is aimed at providing numerical tests to show the effectiveness and applicability of the proposed method. First, we confirm the theoretical convergence order on the example of a problem with a known exact solution in Section 3.1. Finally, we show the applicability of the proposed method to solving three filtration problems in Sections 3.2 and 4 .
Although the boundary condition of the first kind for η is considered in the analysis, the results obtained can be extended to other kinds of boundary conditions with minor modifications. Further, we assume that η d ≡ 0 for simplicity of presentation.
The following assumptions are used throughout the paper. Assumption 1. Problem 1 has a unique solution that has a sufficient number of derivatives required for the analysis.

Derivation of a Fully Discrete Scheme
To construct a numerical scheme, we first discretize the time interval J by points t n = nτ, n = 0, 1, ..., N such that Nτ = T, where τ > 0 is the time discretization parameter, and denote f n (x) = f (x, t n ). In addition, we use the notation u n and p n due to the implicit temporal dependence.
At first, we approximate the time derivative and the integral term in (2). To achieve the global second order in temporal variable, we employ the following approximations of the time derivative at t = t n : , ξ n ∈ (t n−2 , t n ). Next, to approximate the third term in (2), represent the integral as below: where r n 3 = O τ 2 . Thus, the third term at t = t n in (2) is approximated as follows: Consider the differences Note that r n 7 = O τ 2 . Furthermore, Assumption 2 implies that r 1 6 = O τ 2 . A similar assertion is valid concerning the components of the vectors r n 5 and r 1 4 due to Assumption 4. Let us rewrite Problem 2 with the use of the obtained approximations.
.., N satisfying the following identities for any w ∈ W, v ∈ U, q ∈ M: when n = 1: when n ≥ 2: Next, we introduce a quasi-uniform triangulation T h in Ω with a diameter h and let W h ⊂ W, U h ⊂ U, M h ⊂ M be certain finite element spaces. Further, truncating the errors in (3)-(8), we define a fully discrete finite element procedure.
when n ≥ 2: Equations (9) and (12) are in fact nonlinear. We apply Newton's method for their linearization. To this end, we iterate the above procedure for the iterative parameter s = 1, 2, ... and denote the value of the function η n h on the sth iteration by η n,s h . Then, the iterative method is defined as follows.
when n ≥ 2: Thus, Problem 1 is numerically solved in accordance with the following algorithm. First, the pair u 1 , p 1 is determined simultaneously from Equations (16) and (17) employing a known initial value η 0 h . Further, by taking η 1,0 h = η 0 h as an initial approximation, the iterative process for computing η 1,s h , s = 1, 2, ... is conducted by (15) For subsequent time layers, the calculation is performed in a similar manner. For n ≥ 2, the pair (u n , p n ) is determined simultaneously from Equations (19) and (20) The former statement will be rigorously proven in Section 2.4, and the convergence of Newton's method will be shown in Section 2.5.

Stability of the Scheme
In this section, we discuss the stability of the scheme (9)- (14). For the sake of completeness, let us formulate the discrete analogue of the Gronwall's lemma, which will be used several times in obtaining the main results.
where C is a constant independent of τ. (10) and (11) and (v h , q h ) = u n h , p n h in (13) and (14). Then, the inequalities (21) and (22) follow immediately by using Assumptions 2 and 4, the Cauchy inequality, as well as the well-known inequality valid for all positive a, b and ε.
be the solution of Problem 4. The following inequality holds for all τ ≤ τ 1 such that under Assumptions 2 and 3: where C is a constant independent of τ.
Proof. Let us take w h = η n h in (12): and estimate the terms B n 1,m , m = 1, 2, ..., 8. By a direct check, one can verify that B n 1,1 is bounded from below by the quantity Under Assumptions 2 and 3, the condition τ ≤ τ 1 where τ 1 is defined in (24), by using the inequality (23) and noting that e −x < 1 for x > 0, it is easy to obtain that where is an arbitrarily small number. Making use of the inequalities (26)- (32) in (25) and summing the resulting inequality over n from 2 to n, we obtain: Evaluating the remaining dot products using Assumption 2 and noticing that we have from (33) that Then, for all τ ≤ τ 1 , where τ 1 is defined in (24), we obtain By applying the discrete Gronwall's lemma, we arrive at Now, choose w h = η 1 h in (9): Estimating the terms B 2,1 , B 2,2 , ..., B 2,7 in (36) in a similar manner as in obtaining the inequalities (27)-(32), we conclude that, for all τ ≤ τ 1 satisfying the condition (24): By combining (35) and (37), we obtain the inequality Finally, by taking the square root of both parts of the last inequality, then using the fact |a i | and observing that 1 √ n ≥ τ T , we arrive at the assertion of the lemma.
Proof. Lemmas 2 and 3 imply that, under the specified conditions: Applying the discrete Gronwall's lemma to this inequality, we obtain the assertion of the theorem.

Let us introduce the projections
satisfying the approximation properties where m and l depend on the choice of the finite element spaces. Let us decompose the errors as follows: where C is a constant independent of the mesh parameters.
Proof. Take the difference of Equations (7) and (13) and Equations (8) and (14), then consider the sum of the resulting equations for v h = θ n u , q h = θ n p , and take into account the notations (42) and (43) and definitions (39) and (40) to have Estimate the terms in (46) by using Assumptions 2 and 4, and relations (42), (43) as follows: Using these inequalities in (46), we obtain (45). The inequality (44) is derived in a similar manner.
under Assumptions 1-3: where C is a constant depending on the norms of the solution of Problem 1 but independent of the mesh parameters.
Proof. Subtract (3) and (6) from (9) and (12), respectively. Then, take w h = θ n η and use the notation (41) to obtain: when n = 1: when n ≥ 2: Estimate the terms in (49) using Assumptions 2, 3, definitions (38) as follows: Represent B n 4,7 as shown below: then utilize the relations (41) and (42), Assumptions 1 and 2 to obtain The term B n 4,8 is estimated as follows: Using the obtained inequalities in (49) and summing the resulting inequality with respect to n from 2 to n implies: Using the relation (34) and noting that θ 0 η = 0 and we conclude from (50) that Then, for all τ ≤ τ 2 satisfying the condition (47), we have Applying the discrete Gronwall's lemma gives We now turn to estimating the terms in (48) using Assumption 2 and the definition (38): Rewrite B 3,5 as follows: then use (41) and (42), Assumptions 1 and 2, and note that θ 0 η = 0 to have .
Taking into account the obtained inequalities, we conclude from (48) that Under the condition τ ≤ τ 2 , where τ 2 is defined in (47), we obtain By combining (51) and (52), we have Taking the square root of both parts of the inequality, we arrive at the assertion of the lemma.
for all τ ≤ τ 2 such that under Assumptions 1-4, where C is a constant depending on the norms of the solution of Problem 1 but independent of the mesh parameters.
Proof. Lemmas 4 and 5 imply that, for all τ ≤ τ 2 : By applying the discrete Gronwall's lemma to the last inequality, we obtain the assertion of the theorem.

Convergence of Newton's Method
Let us determine the conditions ensuring the quadratic convergence of Newton's iterative process (15), (18) using the technique proposed in [49]. Theorem 3. The iterative process (15) and (18) converges quadratically under Assumption 3 and the condition Cτ for sufficiently small τ.
Proof. Subtract (15) and (18) from (9) and (12) to obtain where θ n,s η = η n h − η n,s h , c 1 = 1 and c n = 3 2 , n ≥ 2. Choose w h = θ n,s η and estimate the terms in (53). Let us dwell in more detail on the estimate of the last term on the left-hand side of (53). By using Assumption 3 and the inequality [49] valid for a differentiable function g with a Lipschitz-continuous derivative, we obtain .
Further, using the inequality θ n,s−1 whence, for sufficiently small τ, we conclude that This inequality implies the assertion of the theorem.

Verification of the Convergence Order
This section provides an analysis of a few numerical tests conducted to validate the proposed numerical method. The purpose of the first computational experiment was to determine the empirical convergence order and compare it with the theoretical convergence order predicted in Theorem 2.
Consider Problem 1 in Q T = Ω × J, where Ω = (0, 1) 2 , J = (0, 1] with the following input data: The exact solution of the problem is as follows: η(x, t) = t 3 sin πx 1 sin πx 2 , One can verify that Assumptions 2-4 are satisfied for the selected data. The following finite element spaces were chosen: where P l is the set of polynomials of degree at most l. Therefore, we expected the second convergence order both in temporal and spatial variables.
A partition containing M 1 = 300 nodes was introduced on each side of the square Ω, and a quasi-uniform triangulation of Ω was generated on its basis. The diameter of the triangulation obtained in this way was h ≈ 0.004714. The iterative process (15) and (18) was performed until η n,s h − η n,s−1 h < 10 −9 . In all numerical tests, this condition was achieved in two to three iterations.
The numbers τ 1 and τ 2 in Theorems 1 and 2 that impose a constraint on the time step at = 0.04 were approximately equal to 0.46875 and 0.11719, respectively. In particular, the stability and convergence of the scheme is predicted to be achieved for all τ ≤ 1/10 according to these theorems. By successively halving the time step in the range from τ = 1/10 to τ = 1/160, we calculated the corresponding errors E τ in the L 2 -norm and the convergence order by the formula r = (log E 2τ − log E τ )/ log 2.
The calculation results are presented in the second and third columns of Table 1. It follows from the table that the actual convergence order slightly exceeded 2.0, which fully confirms the temporal convergence order predicted in Theorem 2.
This problem was also solved by the numerical scheme proposed in [12] for comparison. To this end, a uniform grid is introduced in the square Ω with a step h 1 = h 2 = 1/M 2 with M 2 = 300. Conducting a computational experiment as described above, we obtained the actual order of convergence of this method. The calculation results are presented in the fourth and fifth columns of Table 1, which clearly show that this method converges with the first order in the temporal variable. Thus, it follows from a comparison of the results that the method proposed in this paper is more efficient for solving the problem. A similar computational experiment was conducted by reducing the parameters M 1 and M 2 by half at a fixed time step τ = 0.01 to verify the convergence order with respect to the spatial variable. The results of the experiment showed the second order of convergence for both methods. Therefore, the numerical tests performed fully confirmed the order of convergence predicted in Theorem 2.

Modeling the Two-Phase Non-Equilibrium Flow in Porous Media
We now consider a more realistic example. Let us consider the following equations describing non-equilibrium flow of a two-phase incompressible fluid in a porous medium [12,13]: where η and s are the effective and true water saturations, respectively; u is the volumetric velocity of the flow; p is the pressure; K 0 (x) and m are the absolute permeability and porosity of the medium; σ is the replacement time; ν is the degree of non-equilibrium; φ is the porosity; and k w and k o are the relative phase permeabilities corresponding to the phases of water w and oil o. Equation (57) extends the non-equilibrium law by G. Barenblatt [50] and is referred to as the generalized law of non-equilibrium in [12]. The essence of this generalization lies in the assumption that all flows are nonequilibrium and differ only in the degree of nonequilibrium. The function ν determines the degree of nonequilibrium and takes values between 0 and 1, where the limit values ν = 1 and ν = 0 correspond to a completely nonequilibrium flow and a completely equilibrium flow, respectively.
We adopt a few assumptions regarding physical data. In contrast to [12], we exclude the explicit dependence of the capillary pressure function p c on x assuming its dependence only on η. In addition, the porosity and replacement time are assumed to be constants. In this case, the coefficients in (54)-(57) are defined as follows: a(x, η) = − dp c dη where ρ w and ρ o are the water and oil densities, and g is the gravity vector. Further, we employ a slightly modified LET-type model [51] to determine the relative phase permeabilities: with, say, α = β = 2, and assume that the capillary pressure function is defined as where A and B are positive constants. In this case, The calculations were performed on a mesh with a diameter h ≈ 0.0174875. Establishing the exact values of the constants τ 1 and τ 2 in Theorems 1 and 2 was found to be quite difficult. However due to the smallness of K 0 entering F, the restriction on τ was insignificant. In practice, the numerical scheme has shown itself to be stable for a wide range of time step values. The numerical test was conducted for τ = 10 −3 until the flow became steady. The values of the obtained effective saturation η at different times are illustrated in Figure 1. Note that discussion of the simulation results is beyond of this paper's remit since the aim of this study is to show the applicability of the numerical method. This issue is discussed in detail in [39].

Discussion
Based on the theoretical analysis performed and the results of numerical experiments, the following conclusions can be drawn.
(1) The proposed computational method makes it possible to obtain an approximate solution to the filtration problem, which includes a nonlinear partial integro-differential equation. The results of computational experiments conducted for a model problem with different mesh configurations fully confirmed the results of the theoretical analysis.
We demonstrated that the empirical convergence order in the temporal variable in the L 2 -norm slightly exceeded 2.00, while the spatial convergence order was consistently close to 2.00 for the chosen finite element spaces. This result agrees well with the results of Theorem 2 and significantly improves the result obtained in [12].
(2) It should be noted that the construction of a numerical scheme of the second convergence order in the temporal variable relies heavily on the properties of the functions adopted in classical filtration models. In most two-phase relative phase permeability models (see, for example, [52]), the function k w and its derivative vanishes at 0. Moreover, it is often assumed that the water saturation is uniformly distributed over the domain at the initial time. These simplifying considerations have allowed us to formulate Assumption 2(b) and Assumption 4(a).
(3) The proposed approach based on the combined use of Newton's method and a second-order approximation made it possible to effectively handle nonlinear terms, thereby, maintaining a balance between computational complexity and achieving second-order convergence. Moreover, the results of solving the model problem, presented in Section 3.1, showed the quadratic convergence of Newton's method.
It should be noted that the results obtained have more of a methodological nature. However, the proposed approach can be easily generalized for the case when all the coefficients in Equation (1) depend on the sought solution by applying Newton's method or high-order approximations for their linearization.
(4) Note that the initial boundary value problem for Equation (1) is also of independent interest. Let us show the applicability of the proposed method on a particular case of the nonequilibrium filtration model when the total flow velocity u is equal to zero. This case is known as countercurrent capillary impregnation [53]. Then, Equation (56) becomes an identity, and Equation (54) reduces to the following form: The iterative method for solving this problem is to find η n,s h ∈ W h satisfying the identities (15) and (18). The stability and convergence of a fully discrete scheme follows from Lemmas 3 and 5 with u n h = 0, and the convergence of the iterative method follows from Theorem 3.
Furthermore, the constructed numerical method can be applied to solving a filtration problem in a fractured porous medium. This generalized model is formulated under the assumption that the fluid and medium parameters depend on the fractional derivative of the pressure. Making use of a non-local change of variables allows one to reduce the obtained pressure equation to the non-linear partial integro-differential Equation (1). For a more detailed analysis, we refer the reader to Appendix A.
(5) Since the work is mainly aimed at a comprehensive study of the numerical method for solving Problem 1, little attention is paid to the interpretation of the results of modeling two-phase nonequilibrium fluid flow in porous media. In our opinion, this issue deserves more detailed consideration, which will be the subject of a separate paper. (6) In subsequent works, the authors intend to focus on other aspects of the numerical solution of Problem 1. For example, the obtained results can be extended for the fractional differential generalizations of the filtration model by replacing the time derivatives with their fractional differential counterparts. Similar generalizations of well-known fluid dynamics models have been the subject of many works [54][55][56][57]. In addition, in subsequent papers, the authors intend to investigate the asymptotic stability of the proposed scheme (see, for example, [58][59][60]). Furthermore, the study of control problems with qualitative properties of controllability [61] is also promising.
The outcomes of this study can be used in solving other classes of problems involving integro-differential equations. The proposed approach also forms the basis for further research in this direction.
Thus, the problem of fluid flow in a fractured porous medium under study has been reduced to solving a nonlinear partial integro-differential equation studied in the paper. The numerical solution of the obtained equation is similar to that described in the Discussion section for Equation (58).