Integral Transforms Solution of Immiscible Two-Phase Flow in Porous Media with Capillary Pressure Effects

1. Introduction

Secondary oil recovery plays a critical role in hydrocarbon extraction, aiming to enhance recovery beyond the capabilities of primary methods. After the initial production phase, which relies on the reservoir’s natural energy, secondary recovery techniques, such as water and-or gas injection, are employed to sustain or restore reservoir pressure. These methods not only improve extraction efficiency and extend the productive lifespan of oil fields but also contribute to reducing the environmental footprint of oil exploration and production while optimizing the economic viability of recovery operations [1]. In this context, for more than half a century, waterflooding has served as a reliable method of maintaining reservoir pressure and increasing oil production efficiency [2,3,4,5,6,7].

Given the importance of advanced oil recovery methods, laboratory scale models are critical for understanding not only the dynamics of multiphase flows in heterogeneous porous media, but also the properties of the reservoir rock, to properly plan water floods [8,9]. Coreflood experiments, for example, allow researchers to simulate and study flow behavior under controlled conditions, revealing important information about the mechanisms of oil displacement and the effectiveness of various recovery techniques [10]. These experiments provide important data for inverse problems, allowing the determination of intrinsic properties of the plug subjected to fluid flow, such as relative permeabilities [11,12,13,14,15,16,17].

A significant body of research has focused on numerically predicting water saturation dynamics and oil recovery in porous media, aiming to elucidate the complexities of multiphase flow in both multidimensional [18,19,20] and one-dimensional frameworks [21,22,23]. One-dimensional models, such as the Buckley-Leverett formulation and its variations [12,24,25,26], have been widely used to study two-phase flow in porous media. These simplified formulations effectively capture the fundamental physics of the problem while offering cost effective solutions, making them invaluable for understanding flow behavior in complex systems.

However, since these solutions neglect the influence of capillary pressure, they fail to accurately represent flows where this effect plays a significant role. Incorporating such terms introduces a strong nonlinear behavior into the problem. Nonlinear convective and diffusive effects and-or nonlinear source terms are frequently encountered in transport phenomena modeling, often resulting in nonlinear partial differential equations. These nonlinear formulations are generally unsolvable using traditional analytical methods [27,28,29,30]. Instead, they typically require approximate numerical approaches [31,32,33,34] or the application of mathematical techniques that involve dependent or independent variable transformations in specific formulations or approximate analytical solutions of limited scope [35,36,37,38].

As a result, numerical approaches have become the dominant choice for tackling nonlinear problems, supported by a range of effective algorithms and schemes. Nonetheless, there is a growing need to further improve the efficiency of computationally intensive processes, particularly in applications involving nonlinear formulations. These challenges are especially relevant in tasks such as inverse problems, optimization and uncertainty quantification, due to the usually required large number of runs of the direct problem solution, where the demand for scalability and computational feasibility continues to rise [39].

In response to the limitations of classical exact or approximate analytical solution techniques, the development of hybrid approaches that combine analytical and numerical techniques has gained significant traction. Among these, one particularly notable and mature methodology is the Generalized Integral Transform Technique (GITT). Initially introduced by Cotta and Ozisik [40] and subsequently extended through numerous contributions [41,42,43,44,45,46,47,48,49,50,51,52,53], GITT generalized the classical integral transform method originally developed for linear diffusion problems [27,28,29]. This meshless technique demonstrates the effective combination of the accuracy and computational efficiency of analytical solutions with the flexibility and generality of numerical solutions to tackle intricate engineering challenges, yielding error controlled solutions with moderate computational expense, allowing seamless incorporation with different tasks that demand extensive direct model evaluations.

The application of the Generalized Integral Transform Technique (GITT) to nonlinear problems is based on the idea that nonlinear partial differential equations can be reformulated to confine the nonlinear terms within the source terms, both in the governing equations and boundary conditions. This reformulation preserves linear characteristic coefficients, which define the associated eigenvalue problem. Through eigenfunction expansion and the integral transformation process, the reformulated nonlinear partial differential equation is transformed into a system of coupled, nonlinear, infinite transformed ordinary differential equations. This transformed ODEs system is then numerically solved to determine the transformed potentials, facilitating the analysis of complex nonlinear dynamics in the original problem.

Over the past three decades, following the introduction of GITT for nonlinear problems, this methodology has been progressively expanded to address a wide range of challenges. These extensions include applications to nonlinear models featuring variable physical properties, nonlinear source terms in both governing equations and boundary conditions, moving boundary problems, irregular domains, eigenvalue problems, nonlinear convection, diffusion or reaction-dissipation effects [41,42,43,45,46,47,48,49,51,54]. Recent applications of GITT that involve flow and transport in porous media can be cited such as pumping in confined, leaky, and unconfined aquifers [55], flow on fractured heterogeneous porous media [56], and flow in variably saturated porous media [57]. In reservoir engineering, GITT has been used in three-dimensional transient flow in heterogeneous reservoirs [58], two-dimensional energy balance for the transient temperature of sand faces in wells with mixed production [59], analysis of one-dimensional oil displacement by water injection in a core plug [60], and the two-dimensional convection-diffusion problem for tracer flow in a five-spot pattern [61,62].

Based on the discussion above and considering the potential of GITT and the importance of the addressed topic, this work presents a hybrid solution via GITT for immiscible two-phase flow in porous media with capillary pressure effects, which are highly relevant for reservoir studies. The formulation considers a two-phase flow in heterogeneous porous media, presenting linear and classical boundary conditions, which are independent of capillary pressure. In this demonstration, the effect of water injection at the rock’s boundary is expressed by a Dirichlet boundary condition, where the water saturation is at its maximum. The Kirchhoff transformation [28,37,63,64,65,66] is used to incorporate the nonlinearity of the diffusive operator originating from the capillary pressure into the transformed dependent variable to prepare the governing equation for the integral transformation via GITT. The test-case results here presented focus on the dynamics of the transformed and reconstructed solutions, highlighting convergence behavior, verification against an approximate numerical solution, and discusses the accuracy of low-order approximations.

2. Materials and Methods

2.1. Problem Formulation

Consider the two phase oil-water flow problem in a heterogeneous porous medium as shown in Eq. 1 [67]:

$$\varphi \left(X\right){\displaystyle \frac{\partial S_w}{\partial t}}+{\displaystyle \frac{u_t}{L}}{\displaystyle \frac{\partial f_w\left(S_w\right)}{\partial X}}={\displaystyle \frac{1}{L^2}}{\displaystyle \frac{\partial }{\partial X}}\left(D_c\left(S_w\right){\displaystyle \frac{\partial S_w}{\partial X}}\right),$$

(1)

subject to the following boundary conditions (Eq. 2a and 2b) and initial condition (Eq. 2c):

$$\begin{array}{cc}\hfill S_w(0,t)& =1-S_{or}\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill {\left.{\displaystyle \frac{\partial S_w(X,t)}{\partial X}}\right|}_{X=1}& =0\hfill \end{array}$$

(2b)

$$\begin{array}{cc}\hfill S_w(X,0)& =S_{wi}\hfill \end{array}$$

(2c)

In these, X represents the dimensionless spatial coordinate, such that the dimensional spatial coordinate x is defined as $x=XL$, with L being the length of the core, t the time, $S_w$ the water saturation such that $S_w\equiv S_w(x,t)$, $\varphi \left(X\right)$ the spatial porosity distribution, $u_t$ the injected flow velocity, and $f_w\left(S_w\right)$ the water fraction (Eq. 3) given as:

$$f_w\left(S_w\right)={\displaystyle \frac{{\displaystyle \frac{k{r}_w\left(S_w\right)}{{\mu }_w}}}{{\displaystyle \frac{k{r}_w\left(S_w\right)}{{\mu }_w}}+{\displaystyle \frac{k{r}_o\left(S_w\right)}{{\mu }_o}}}},$$

(3)

where $k{r}_w\left(S_w\right)$ and $k{r}_o\left(S_w\right)$ are the relative permeabilities of water and oil, respectively, and ${\mu }_w$ and ${\mu }_o$ are the dynamic viscosities of water and oil, respectively. The term $D_c\left(S_w\right)$, which is expressed as:

$$D_c\left(S_w\right)\equiv {\lambda }_o\left(S_w\right)f_w\left(S_w\right)\left(-{\displaystyle \frac{d{P}_c\left(S_w\right)}{d{S}_w}}\right),$$

(4)

is the nonlinear diffusion coefficient due to capillary pressure, with ${\lambda }_o\left(S_w\right)$ being the oil mobility and $P_c\left(S_w\right)$ the capillary pressure, and $S_{wi}$ and $S_{or}$ the residual water and oil saturation, respectively.

2.2. Formal Solution

Let us apply the Kirchhoff transform [27,28,66] to the dependent variable $S_w$, such that:

$$U(X,t)={\int }_{S_w^{*}}^{S_w}{\displaystyle \frac{D_c\left(S_w\right)}{D_c\left(S_w^{*}\right)}}d{S}_w,$$

(5)

where $S_w^{*}$ is a reference water saturation. Then:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial U}{\partial S_w}}={\displaystyle \frac{D_c\left(S_w\right)}{D_c\left(S_w^{*}\right)}}\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial S_w}{\partial U}}={\displaystyle \frac{D_c\left(S_w^{*}\right)}{D_c\left(S_w\right)}},\hfill \end{array}$$

(6b)

one obtains:

$$\begin{array}{cc}\hfill \varphi \left(X\right){\displaystyle \frac{\partial S_w}{\partial U}}{\displaystyle \frac{\partial U}{\partial t}}+{\displaystyle \frac{u_t}{L}}{\displaystyle \frac{\partial f_w}{\partial X}}& ={\displaystyle \frac{1}{L^2}}{\displaystyle \frac{\partial }{\partial X}}\left(D_c\left(S_w\right){\displaystyle \frac{\partial S_w}{\partial U}}{\displaystyle \frac{\partial U}{\partial X}}\right)\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill \varphi \left(X\right){\displaystyle \frac{\partial U}{\partial t}}+{\displaystyle \frac{u_t}{L}}{\displaystyle \frac{d{f}_w}{d{S}_w}}{\displaystyle \frac{\partial U}{\partial X}}& ={\displaystyle \frac{D_c\left(S_w\right)}{L^2}}{\displaystyle \frac{{\partial }^{2}U}{\partial X^2}}\hfill \end{array}$$

(7b)

The transformed boundary and initial conditions are given by:

$$\begin{array}{cc}\hfill U(0,t)& =U_{or},\hfill \end{array}$$

(8a)

$$\begin{array}{cc}\hfill {\left.{\displaystyle \frac{dU(X,t)}{dX}}\right|}_{X=1}& =0,\hfill \end{array}$$

(8b)

$$\begin{array}{cc}\hfill U(X,0)& =U_{wi},\hfill \end{array}$$

(8c)

where:

$$\begin{array}{cc}\hfill U_{or}& \equiv {\int }_{S_w^{*}}^{1-S_{or}}{\displaystyle \frac{D_c\left(S_w\right)}{D_c\left(S_w^{*}\right)}}d{S}_w\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill U_{wi}& \equiv {\int }_{S_w^{*}}^{S_{wi}}{\displaystyle \frac{D_c\left(S_w\right)}{D_c\left(S_w^{*}\right)}}d{S}_w\hfill \end{array}$$

(9b)

A linear filter $U_F(X;t)$ is used to homogenize the boundary condition at $X=0$, according to the expression:

$$U(X,t)=U_F(X;t)+U^{*}(X,t),$$

(10)

in which $U^{*}(X,t)$ stands for the filtered potential, in the form:

$$U_F(X,t)=aX+b$$

(11)

The selected eigenvalue problem is based on the homogeneous version of the above filtered problem, such that:

$$\begin{array}{c}\hfill {\displaystyle \frac{d^2{\Psi }_i\left(X\right)}{d{X}^2}}+{\mu }_i^2{\Psi }_i\left(X\right)=0,\end{array}$$

(12a)

$$\begin{array}{c}\hfill {\Psi }_i\left(0\right)=0,\end{array}$$

(12b)

$$\begin{array}{c}\hfill {\left.{\displaystyle \frac{d{\Psi }_i\left(X\right)}{dX}}\right|}_{X=1}=0,\end{array}$$

(12c)

which can be solved analytically. In the eigenvalue problem, the eigenfunctions are normalized using the following normalization integrals:

$${\int }_0^1{\tilde{\Psi }}_i\left(X\right){\tilde{\Psi }}_j\left(X\right)dX={\delta }_{i,j},$$

(13)

where i and j are integers, ${\delta }_{i,j}$ is the Kronecker delta function, and the normalized eigenfunctions are written as:

$${\tilde{\Psi }}_i\left(X\right)={\displaystyle \frac{{\Psi }_i\left(X\right)}{N_i^{1/2}}},N_i={\int }_0^1{\Psi }_i^2\left(X\right)dX,$$

(14)

where $N_i$ is the normalization integral.

Based on the orthogonality relations, the transform-inverse pair is defined for the dependent variable:

$$\begin{array}{cc}\hfill {\overline{U}}_i^{*}\left(t\right)& ={\int }_0^1{U}^{*}(X,t){\tilde{\Psi }}_i\left(X\right)dX\left(\mathrm{Transform}\right)\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill U^{*}(X,t)& =\sum _{i=1}^{\infty }{\overline{U}}_i^{*}\left(t\right){\tilde{\Psi }}_i\left(X\right)\left(\mathrm{Inverse}\right)\hfill \end{array}$$

(15b)

If one applies the integral transform to Eq. 7b by integrating all terms in X from 0 to 1 after multiplying by the eigenfunction of the eigenvalue problem $\left({\int }_0^1\square {\tilde{\Psi }}_i\left(X\right)dX\right)$, one obtains:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\overline{U}}_i^{*}\left(t\right)}{dt}}& =\sum _{j=1}^{\infty }{B}_{ij}{\overline{U}}_j^{*}\left(t\right)+\sum _{j=1}^{\infty }{C}_{ij}{\overline{U}}_j^{*}\left(t\right),\hfill \end{array}$$

(16a)

$$\begin{array}{cc}\hfill {\overline{U}}_j^{*}\left(0\right)& ={\overline{U}}_{or,j}^{*},\hfill \end{array}$$

(16b)

where:

$$\begin{array}{cc}\hfill B_{ij}& =-{\displaystyle \frac{u_t}{L}}{\int }_0^1{\displaystyle \frac{1}{\varphi \left(X\right)}}{\displaystyle \frac{d{f}_w\left(S_w\right)}{d{S}_w}}{\tilde{\Psi }}_i\left(X\right){\tilde{\Psi }}_j{\left(X\right)}^{\prime }dX,\hfill \end{array}$$

(17a)

$$\begin{array}{cc}\hfill C_{ij}& =-{\displaystyle \frac{{\mu }_j^2}{L^2}}{\int }_0^1{\displaystyle \frac{1}{\varphi \left(X\right)}}{D}_c\left(S_w\right){\tilde{\Psi }}_i\left(X\right){\tilde{\Psi }}_j\left(X\right)dX,\hfill \end{array}$$

(17b)

$$\begin{array}{cc}\hfill {\overline{U}}_{or,j}^{*}& ={\int }_0^1{U}^{*}(X,0){\tilde{\Psi }}_j\left(X\right)dX,\hfill \end{array}$$

(17c)

$$\begin{array}{cc}\hfill U^{*}(X,0)& =U_{wi}-U_F(X,0),\hfill \end{array}$$

(17d)

where the matrices $B_{ij}$ and $C_{ij}$ are nonlinear and thus time-dependent. Once the system presented on Eq. 16a is numerically solved for the transformed potentials along the time variable, the inverse formula (Eq. 15b) is recalled to reconstruct the potential U, which is finally inverted to provide the water saturation distribuition, $S_w$, making use of the nonlinear relation (Eq. 5).

Depending on the functional form of the coefficient $D_c\left(S_w\right)$, the Kirchoff transformation, Eq. 5, shall require numerical integration. In this case, it is cost-effective to construct an interpolating function that provides a straightforward transformation within the domain of validity of the dependent variable $S_w$. As for the inversion, again for some functional forms of $D_c$, it may result in an analytical expression, otherwise a nonlinear algebraic equation roots algorithm needs to be employed. Again, an interpolating function construction shall be a cost-effective solution for the inversion process.

2.3. Test-Case

For the application of the solution presented in this work, a test case was selected based on the study by [68], specifically Case 3, in which $p_c\left(S_w\right)\ne 0$. The capillary pressure is provided in the form of a data table (blue dots), which was interpolated using 5^th order splines (black curve), as shown in Figure 1, for both the raw data and the interpolated curve within the domain of validity ($S_w<S_{wi}$ and $S_w>1-S_{or}$). The choice of 5^th order splines led to the least numerical oscillations in the derivative ${\displaystyle \frac{d{P}_c}{d{S}_w}}$ estimates, but artificial smoothing can also be implemented at the cost of loosing the exact interpolation for $P_c$ in all the data points.

The extrapolation beyond the original data domain, i.e., the enforcement of constant values for $P_c\left(S_w\right)$ at the extremes outside the valid range, is required because, during the numerical solution process, the potential function $U(X,t)$ may assume values slightly outside the dependent variable domain, due to lower truncation orders relatively to the time value being considered. A similar approach was adopted by [21] and is discussed in the GitHub repository by [69]. If $P_c\left(S_w\right)$ were to be estimated through a direct extrapolation of the fitted curve, unrealistic and physically inconsistent values might arise, resulting in amplification of the numerical error in the solution. Once the interpolation function for $P_c\left(S_w\right)$ is constructed, the corresponding diffusion coefficient $D_c\left(S_w\right)$ is computed according to Eq. 4.

Table 1 presents the input data used for the test case. The relative permeabilities are defined according to the Brooks & Corey relationship [70], given by:

$$k{r}_w\left(S_w\right)=k{r}_w^{\max }{\left({\displaystyle \frac{S_w-S_{wi}}{1-S_{wi}-S_{or}}}\right)}^{n_w},k{r}_o\left(S_w\right)=k{r}_o^{\max }{\left({\displaystyle \frac{1-S_w-S_{or}}{1-S_{wi}-S_{or}}}\right)}^{n_o},$$

(18)

where $S_{wi}$ and $S_{or}$ represent the irreducible water and oil saturations, respectively, and $n_w$ and $n_o$ are real-valued exponents fitted to experimental observations. The maximum relative permeabilities of water and oil are denoted by $k{r}_w^{\max }$ and $k{r}_o^{\max }$, respectively.

Figure 2 provides a graphical representation of the Kirchhoff transformation and its inversion for this test-case. It should be pointed out that for this particular behavior of the coefficient $D_c$, in consequence of the $P_c$ behavior, while the transformation to the U domain is a fairly smooth function, the inversion to the $S_w$ domain is quite steep close to the domain limits, $S_{wi}$ and $1-S_{or}$. This requires a refinement of the interpolating points near these borders to ensure a proper accurate interpolation, otherwise a small deviation in the U value could cause an amplified deviation in the $S_w$ value and cause numerical weegles to appear in the representation of the water saturation field.

2.4. Shanks Transformation

The Shanks transformation [71,72] is a non-linear sequence-to-sequence acceleration method commonly used to improve the convergence of slowly convergent or even divergent series. Consider that each term of the eigenfuntion expansion represents one term $\left\{a_m\right\}$ in a given sequence. From this sequence, one can readily form the sequence of partial sums $\left\{A_n\right\}$, up to the truncation order N. The first-order Shanks transform is defined as :

$$S\left(A_n\right)={\displaystyle \frac{A_{n+1}{A}_{n-1}-A_n^2}{A_{n+1}-2A_n+A_{n-1}}}=A_{n+1}-{\displaystyle \frac{{(A_{n+1}-A_n)}^2}{A_{n+1}-2A_n+A_{n-1}}},$$

(19)

where $S\left(A_n\right)$ represents an extrapolated estimate of the sequence’s limiting value. This transformation effectively removes the leading-order error term under the assumption that the sequence behaves locally as a linear combination of exponentials, thus providing an accelerated approximation to the series limit. Further convergence acceleration may be obtained by repeating the use of the Shanks transformation, applied to the new sequence $\left\{S\left(A_n\right)\right\}$, thus providing a second-, third-, and higher-order Shanks transforms. In this work, the Shanks transformation is employed to accelerate, when needed, or at least confirm the convergence of the computed GITT results.

3. Results

3.1. Convergence Analysis

This section presents the convergence analysis of the proposed integral transform solution. The truncation order in the eigenfunction expansion is initially fixed at $N_T=12$, and the results are examined for different number of terms in the resulting series, $N\in \{2,4,6,8,10,12\}$ in order to inspect the convergence behavior of the proposed expansion. All numerical integrations for the transformed system matrices and initial conditions vector are performed using Gaussian quadrature with $N_g=80$ points, ensuring that integration errors remain sufficiently low and controlled below the desired four digits accuracy in the final results throughout the analysis.

To accelerate and further verify the expansion convergence, the first-level Shanks transformation is applied to the partial sums of each sequence. The Shanks method provides an accelerated estimate of the series limit, which also serves as an additional consistency check for the computed values. Given the favorable convergence behavior already observed for moderate values of N, the contribution of the Shanks transformation in accelerating convergence is useful mostly at lower time values, yet confirms the robustness of the solution procedure.

The numerical results are summarized in Table 2, Table 3 and Table 4, for three distinct time instants ($t\in \{300,600,1000\}s$). In all cases, the solution $U(X,t)$ is evaluated at five representative spatial positions ($X\in \{0.1,0.3,0.5,0.7,0.9\}$). A clear trend of convergence is observed as N increases, with full convergence to four significant digits at the highest value of N (± 1 in the last digit given), as confirmed by the Shanks transformation.

To complement the convergence study, a relative error analysis is conducted to quantify the accuracy of the truncated eigenfunction expansions. The objective is to assess how the partial sum of N terms compares to a reference solution obtained with $N_T=12$ terms. The relative residue estimate ${\epsilon }_{\mathrm{rel}}(X,t)$ is defined as:

$${\epsilon }_{\mathrm{rel}}=\left|{\displaystyle \frac{{\displaystyle \sum _{k=1}^N{U}_k(X,t)}}{{\displaystyle \sum _{k=1}^{N_T}{U}_k(X,t)}}}-1\right|.$$

(20)

The results presented on Figure 3, Figure 4 and Figure 5, for five different dimensionless positions and three different times, demonstrate that even for relatively low truncation levels, such as $N=6$, the relative residue remains below 1% at all relevant spatial and temporal positions. This highlights the efficiency of the proposed methodology in capturing the dominant dynamics of the system with a reduced computational cost.

The convergence tables presented so far refer to the values of $U(X,t)$ computed at selected spatial positions and time instants. However, from a physical point of view, one needs to analyze the behavior of the water saturation $S_w(X,t)$, as this variable directly represents the fluid distribution within the porous medium. For completeness, we therefore provide below the corresponding convergence tables for $S_w(X,t)$, evaluated at the same dimensionless spatial positions and time instants used previously for $U(X,t)$.

These results are obtained by post-processing the values of $U(X,t)$ through a numerical interpolation of the inverse Kirchhoff transform, but were also critically compared with the direct inversion through the nonlinear algebraic equation solver, without making use of the interpolated curve, and no difference has been observed with respect to the results shown below.

Table 5. Convergence of $S_w(X,t)$ for $t=300$ s and N up to 12.

N	$\mathit{X}=0.1$	$\mathit{X}=0.3$	$\mathit{X}=0.5$	$\mathit{X}=0.7$	$\mathit{X}=0.9$
2	0.5608	0.4669	0.4043	0.3421	0.2
4	0.5578	0.4665	0.4056	0.3414	0.2
6	0.5579	0.4665	0.4056	0.3413	0.2
8	0.5584	0.4660	0.4055	0.3431	0.2
10	0.5586	0.4662	0.4056	0.3426	0.2
12	0.5586	0.4662	0.4055	0.3424	0.2
Shanks	0.5586	0.4662	0.4055	0.3424	0.2

Table 5. Convergence of $S_w(X,t)$ for $t=300$ s and N up to 12.

N	$\mathit{X}=0.1$	$\mathit{X}=0.3$	$\mathit{X}=0.5$	$\mathit{X}=0.7$	$\mathit{X}=0.9$
2	0.5608	0.4669	0.4043	0.3421	0.2
4	0.5578	0.4665	0.4056	0.3414	0.2
6	0.5579	0.4665	0.4056	0.3413	0.2
8	0.5584	0.4660	0.4055	0.3431	0.2
10	0.5586	0.4662	0.4056	0.3426	0.2
12	0.5586	0.4662	0.4055	0.3424	0.2
Shanks	0.5586	0.4662	0.4055	0.3424	0.2

Table 6. Convergence of $S_w(X,t)$ for $t=600$ s and N up to 12.

N	$\mathit{X}=0.1$	$\mathit{X}=0.3$	$\mathit{X}=0.5$	$\mathit{X}=0.7$	$\mathit{X}=0.9$
2	0.5888	0.5096	0.4670	0.4405	0.4271
4	0.5876	0.5092	0.4674	0.4407	0.4269
6	0.5874	0.5093	0.4673	0.4408	0.4269
8	0.5873	0.5092	0.4673	0.4408	0.4269
10	0.5873	0.5092	0.4673	0.4408	0.4269
12	0.5873	0.5092	0.4673	0.4408	0.4269
Shanks	0.5873	0.5092	0.4673	0.4408	0.4269

Table 6. Convergence of $S_w(X,t)$ for $t=600$ s and N up to 12.

N	$\mathit{X}=0.1$	$\mathit{X}=0.3$	$\mathit{X}=0.5$	$\mathit{X}=0.7$	$\mathit{X}=0.9$
2	0.5888	0.5096	0.4670	0.4405	0.4271
4	0.5876	0.5092	0.4674	0.4407	0.4269
6	0.5874	0.5093	0.4673	0.4408	0.4269
8	0.5873	0.5092	0.4673	0.4408	0.4269
10	0.5873	0.5092	0.4673	0.4408	0.4269
12	0.5873	0.5092	0.4673	0.4408	0.4269
Shanks	0.5873	0.5092	0.4673	0.4408	0.4269

Table 7. Convergence of $S_w(X,t)$ for $t=1000$ s and N up to 12.

N	$\mathit{X}=0.1$	$\mathit{X}=0.3$	$\mathit{X}=0.5$	$\mathit{X}=0.7$	$\mathit{X}=0.9$
2	0.6318	0.5687	0.5380	0.5213	0.5139
4	0.6304	0.5683	0.5384	0.5215	0.5137
6	0.6302	0.5684	0.5383	0.5215	0.5137
8	0.6301	0.5684	0.5384	0.5215	0.5137
10	0.6301	0.5684	0.5384	0.5215	0.5137
12	0.6301	0.5684	0.5384	0.5215	0.5137
Shanks	0.6301	0.5684	0.5384	0.5215	0.5137

Table 7. Convergence of $S_w(X,t)$ for $t=1000$ s and N up to 12.

N	$\mathit{X}=0.1$	$\mathit{X}=0.3$	$\mathit{X}=0.5$	$\mathit{X}=0.7$	$\mathit{X}=0.9$
2	0.6318	0.5687	0.5380	0.5213	0.5139
4	0.6304	0.5683	0.5384	0.5215	0.5137
6	0.6302	0.5684	0.5383	0.5215	0.5137
8	0.6301	0.5684	0.5384	0.5215	0.5137
10	0.6301	0.5684	0.5384	0.5215	0.5137
12	0.6301	0.5684	0.5384	0.5215	0.5137
Shanks	0.6301	0.5684	0.5384	0.5215	0.5137

In all examined time instants, the values of $S_w$ exhibit a clear convergence trend as N increases. Convergence to three significant digits can be observed for truncation orders as low as $N=4$ even at the lowest value of time. The application of the Shanks transformation further confirms this observation.

3.2. Verification Analysis

Figure 6 illustrates the spatial evolution of $U(X,t)$ at selected time instants, namely $t\in \{100,250,400,550,700,1000\}$ s. These time points were chosen to encompass both the initial transient regime and the post-breakthrough phase of the saturation front. The results show that the solution maintains a smooth and well-behaved structure throughout the spatial domain for all times considered. Only for the lowest time value ($t=100s$) it can be noticed, after careful inspection, that the transition region after the saturation front might require a few more terms in the expansion.

Figure 7 presents a verification analysis performed by comparing the proposed solution with a numerical solution obtained using NDSolve, Mathematica’s built-in PDE solver. Due to numerical instabilities encountered when solving the full nonlinear PDE directly, it was necessary to benefit from the Kirchoff transformation as well. Even then, NDSolve was not able to handle the nonlinearity due to the capillary pressure, and it was necessary to introduce a simplification in the diffusive term, the nonlinear coefficient $D_c\left(S_w\right)$ was evaluated using the known solution obtained by the present method itself. This semi-coupled strategy yielded an approximate numerical result that still serves as a meaningful result for verification purposes.

The comparison revealed a good agreement between the proposed solution and the approximate numerical solution with NDSolve, particularly at early and intermediate times. As expected, a gradual loss of accuracy is observed at larger times, which is attributable to the approximation introduced and propagated in the numerical solution by NDSolve. in the numerical model. Nevertheless, the level of agreement remains sufficient to verify the GITT-based solution.

3.3. Low Order Truncation Solution

In the context of inverse problems and other computationally intensive tasks that require multiple evaluations of the direct problem solution, the use of low-order truncation expansions is particularly attractive, in combination with the Approximation Error Model (AEM) methodology [73,74], as it significantly reduces computational cost, without compromising the overall precision in the estimations. In this section, we investigate the behavior of the solution when the number of terms is limited to $N=6$.

It is important to emphasize that when the eigenfunction expansion is restricted to low order truncation such as $N=6$, it is inherently assumed that all transformed components $U_k(X,t)$ for $k>6$ are identically zero. On the other hand, one may also extract the truncated expansions for instance up to $k=6$ from a computed solution with $N_T=12$, in which case the contributions from higher-order modes ($k=7$ to 12) are present in the solution space but excluded from the partial sum. As such, these two approaches yield different results, and this distinction becomes relevant when evaluating the implications of model truncation in inverse analyses, combined with the AEM approach.

To illustrate the behavior of a low truncation expansion, we compare three sets of results for three representative time instants: (i) direct computations with increasing N up to 6 (with $N_T=6$); (ii) the Shanks-transformed value at $N_T=6$; and (iii) the truncated expansion at $N=6$ extracted from the $N_T=12$ truncated expansion. Table 8, Table 9, and Table 10 report the values of $S_w(X,t)$ for each time value $(t=300,600,1000\mathrm{s})$ at five spatial positions.

The comparison reveals that the direct $N=6$ solutions are remarkably close to the Shanks accelerated estimates, which confirms that the expansion with $N_T=6$ is essentially converged. Nevertheless, when comparing with the truncated sum at $N=6$ of the expansion with $N_T=12$ terms, the deviations are more noticeable, though always below 2.5% from these tables, but for accurate estimates should be accounted for by the approximation error model when implementing the inverse problem analysis.

This observation supports the use of low-order truncations, such as $N=6$ in the present example, as efficient and reliable. Therefore, the results suggest that a truncated solution with $N_T=6$ can serve as a robust surrogate for the full model in inverse formulations, enabling efficient parameter estimation while preserving the essential physical features of the system.

In sequence, we analyze the influence of the numerical integration accuracy on the final solution, by varying the number of Gaussian quadrature points $N_g$ used in the evaluation of the transformed terms $U_k(X,t)$. Although all previous results were computed with $N_g=80$ to eliminate integration-induced errors from the convergence study, it is relevant to assess the effect of reducing $N_g$, particularly in the context of inverse problems, when computational savings are always at a premium.

Table 11, Table 12 and Table 13 present the values of $S_w(X,t)$ obtained for $N_T=6$ and $N_g\in \{20,40,60,80\}$, for the same spatial positions and time instants previously analyzed. The Shanks transformation results computed with $N_g=80$ are included for comparison.

The results indicate that the solution converges rapidly with respect to $N_g$, and that $N_g=40$ already provides values very close to the reference with $N_g=80$ and $N_T=6$. At $N_g=40$, the results are already converged to ± 1 in the fourth significant digit, indicating that the numerical integration error becomes irrelevant to within the desired accuracy in the final results.

For applications that require a large number of the direct problem evaluations, such as most parameter estimation in inverse formulations, this reduction in integration cost can lead to substantial performance gains.

4. Discussion

The results presented in the previous sections confirm the robustness, efficiency, and accuracy of the proposed solution for immiscible two-phase flow in porous media based on the Generalized Integral Transform Technique (GITT). The convergence analysis demonstrated that a relatively low number of terms in the truncated eigenfunction expansion ($N=6$) suffices to achieve converged and precise solutions for both the Kirchhoff-transformed variable $U(X,t)$ and the physically meaningful saturation field $S_w(X,t)$. The use of the Kirchhoff transform is handy in partially incorporating the nonlinearity due to the capillary pressure effect into the dependent variable and thus reducing the stiffness in the transformed ODE system numerical solution. Besides, as previously observed in [39], it offers an inherent convergence acceleration bonus, as more information on the original problem nonlinear diffusive operator is aggregated to the dependent variable. These findings are particularly relevant in the context of inverse problems, where computational tractability is a key requirement.

The agreement between the GITT-based solution and the approximate numerical one provided by NDSolve provides verification of the implemented simulation. Despite the simplification that had to be adopted in the numerical simulation, where the nonlinear coefficient in the diffusive term was evaluated from the hybrid GITT solution itself, the match remains consistently good across the spatial and temporal domains. Minor discrepancies at longer times are attributed to the error propagation in the approximate numerical procedure rather than flaws in the proposed methodology.

From a computational standpoint, the analysis of numerical integration effects revealed that Gaussian quadrature with as few as $N_g=40$ points already yields negligible error to the final water saturation results at the lower truncation order or $N_T=6$, indicating that the method is not only accurate but also computationally efficient. This supports the potential use of the present formulation in repeated simulations, such as those encountered in parameter estimation or uncertainty quantification tasks.

Moreover, the truncated low-order reduced model was shown to retain the essential features of the full solution while offering significant reductions in computational cost. The comparison between the truncated models and the higher-order solution $(N_T=12)$ confirms that even the low order truncated expansion preserved accuracy always below 2.5%, which is certainly acceptable for combining with the AEM for inverse problem analysis. This suggests that the approach can serve as a reliable reduced-order model in inverse formulations, where repeated forward evaluations are required.

In the broader context of modeling nonlinear flow in porous media, the present methodology offers a valuable alternative to conventional numerical solvers. It combines the interpretability and structure of analytical approaches with the flexibility of numerical integration, enabling hybrid strategies that are particularly well-suited for data assimilation, control, inverse problems, and optimization scenarios.

Future research may extend the current formulation to include more complex boundary conditions and multi-dimensional geometries. The flexibility of the integral transform approach also allows for its integration with Bayesian inference frameworks, further expanding its applicability in modern reservoir characterization workflows.