Analysis of the interface between two immiscible fluids in porous media is considered as a classic problem with numerous industrial applications. Interface between two fluids could be unstable. In addition to defining dimensionless equations for a generalized geometry, this study developed equations for special cases. For the problem in one dimension, linear steady state analysis was done and it is indicated that steady state condition is maintained as long as the density of the upper fluid is less than the density of the lower fluid. By using the "K-Method", an analytical solution of the Boussinesq equation is presented and the results for constant and linear flux are obtained and discussed. Considering a Power-law flux at the origin and dimensional analysis, the drilling mud invasion into a permeable aquifer problems becomes a nonlinear boundary value problem that have been investigated in this paper. The obtained results show that the K-Method has a good ability for investigation of B-L problems such as this solved problem.

Keywords: interface, immiscible fluids, porous media, kourosh method, boussinesq equation

Evolution of the interface between two immiscible fluids is one of the most complicated problems in the fluid dynamics. The major complicity of this problem is due to difficulty of an explicitly closed that describe the interface movement. Parameter h denote the interface position is to deduce the close differential equation for the function $h\left(x, y, t\right)$ , $x$ and $y$ are the coordinates in horizontal plate, $t$ is the time variable and $h\left(x,y,t\right)$ is the thickness of the lower fluid.

The behavior of the stratified pair "gas-oil" or "oil-water" natural gas-oil reservoir, can be describe by the thickness of the lower fluid, on the modeling of this function h and then for an optimized recovery, it is absolutely necessary to have an appreciate modeling.¹ Whitham² has considered the case that the viscous force for both fluids can be neglected. Also the case when one of the two fluids is inviscid (like as in water-oil flow investigation in Shallow water theory or ground water flow) have been investigated by Dagan,³ Bear⁴ and Barenblatt.⁵ In ground water hydrology, Boussinesq models are involved to describe radial flow from or to well.^6,7 Axisymmetric flow is a feature of ground water systems subject to pumping.^8,9 Invasion of drilling fluid into a permeable bed is an example of involving pumping problems.^10,12 Dussan12 has been investigated the axisymmetric Boussinesq model by similarity variable approach. An early investigation of drilling mud invasion into a permeable aquifer has been provided by Doll.¹³ Calugaru1 proposed a generalized model for description of evolution of the interface between two immiscible fluids in porous media.

Tang¹⁴ studied transient groundwater flow in an unconfined aquifer subject to a constant water variation at the sloping water-land boundary. To characterize the transient groundwater flow, they presented a novel approximate solution to the 1-D Boussinesq equation. They applied the proposed method to various hydrological problems and showed that it can achieve desirable precisions, even in the cases with strong nonlinearity. Bartlet¹⁵ introduced a class of solutions of the nonlinear Boussinesq equation with source/sink terms. They applied their new solution to sloping aquifers and analytical results capture hysteresis between the groundwater level and groundwater flow rate as a function of hillslope characteristics. Lu¹⁶ suggested a homotopy analysis method (HAM) to solve the generalized Boussinesq equation. Due to the two-degree approximate solution of the variable coefficient Boussinesq equation, they showed that the homotopy perturbation method is effective to solve the variable solution equations. Bansal¹⁷ developed a new analytical solution of 2-dimensional linearized Boussinesq equation for approximation of subsurface seepage flow in confined and unconfined aquifers under varying hydrological conditions. He showed that the vertical flow through the base of the aquifer is an important factor in the determination of groundwater mound and cone of depression. Telyakovskiy¹⁸ modeled water injection at a single well in an unconfined aquifer by the Boussinesq equation with cylindrical symmetry. By introducing similarity variables, they reduced the original problem to a boundary-value problem for an ordinary differential equation. Their approximate solution incorporated both a singular part to model the behavior near the well and a polynomial part to model the behavior in the far field.

The problem of invasion by drilling fluid into a permeable empty bed solved by Li.¹⁹ They used the similarity variable and reduced the axisymmetric Boussinesq equation to a nonlinear ordinary differential equation. Their solution which was a sum of a singular term and a Taylor expansion at the wetting front was demonstrated to be highly accurate. Mortensen²⁰ considered the same problem²¹ in spherical coordinates with the prescribed power law point source boundary condition. They constructed an approximate similarity solution to a nonlinear diffusion equation in spherical coordinates. In this paper, first, the general dimension equation for interface position between two immiscible fluids in porous media has been presented. Then, the instability of this equation for one dimensional case by considering a self-similar solution and using the perturbation method has been investigated. And finally, an accurate solution to axisymmetric Boussinesq problem has been obtained. Also, in order to describe deformations of the interface separating two immiscible fluids in porous media, a system of partial differential equations which are true on either sides of the interface is required.

Theory

Figure 1 represents the generalized geometry of the problem in an orthogonal coordinate system. The initial position of the interface between two immiscible fluids is identified as initial assumed value of $t_o$ .When the flow is initiated by a driving force; changes of the interface are specified as a function of space $\overline{x}$ , $\overline{y}$ and time $\overline{t}$ .

Porous media is assumed homogeneous with constant porosity and isentropic. Permeability tensor for such media could be considered a diagonal matrix. By determining the characteristic length of $L$ for horizontal scale and characteristic height of $H$ , the dimensionless equation of changes of the surface could be written as follow [2]:

$\left(t_{el}\frac{\partial }{\partial t}-{\nabla }_{xy}\right)\left(h^2+\frac{{\delta }_1^2{t}_{r}t}{3K_z}{h}^2\frac{\partial h}{\partial t}+{\delta }_2{P}_r\right)=-t_r{t}_{er}\frac{\partial h}{\partial t}$ (1)

$\left(\frac{t_{el}}{{\mu }_r}\frac{\partial }{\partial t}-{\nabla }_{xy}\right)\left(-{\varnothing }^2+\frac{{\delta }_1^2{\rho }_r{t}_{r}t}{3K_z{\mu }_r{\delta }_o}{\varnothing }^2\frac{\partial \varnothing }{\partial t}+{\delta }_2{\delta }_o{\rho }_{r}P\right)=-\frac{{\rho }_r{\delta }_o{t}_r{t}_{el}}{{\mu }_r}\frac{\partial h}{\partial t}$ (2)

$\varnothing =-{\delta }_{o}h+{\delta }_o+1$ (3)

Where dimensionless parameters are specified as below:

$x=\frac{\overline{x}}{L}, y=\frac{\overline{y}}{L}, h=\frac{\overline{h}}{h_o}$

${\delta }_o=\frac{{\overline{h}}_0}{\overline{H}-\overline{h}o},{\delta }_1=\frac{{\overline{h}}_0}{L} , {\delta }_2=\frac{2{\overline{P}}_0}{\overline{\rho }g{\overline{h}}_0}$

$P_r=\frac{\overline{P}}{P_o}, {\rho }_r=\frac{{\overline{\rho }}_I}{{\overline{\rho }}_{II}}, {\mu }_r=\frac{{\overline{\mu }}_I}{{\overline{\mu }}_{II}} , K_r=\frac{{\overline{K}}_z}{{\overline{k}}_x}$ (4)

${\overline{t}}_{{}_{el}}=\frac{{\overline{\mu }}_L{\overline{L}}^2}{{\overline{K}}_x}, t_{gr}=2{\mu }_I\epsilon \frac{{\overline{L}}^2}{{\overline{K}}_x{\overline{\rho }}_{I}g{\overline{h}}_0},$ $t_{el}=\frac{{\overline{t}}_{el}}{{\overline{t}}_c}, t_r=\frac{2\epsilon }{\overline{\rho }g{\overline{h}}_0}=\frac{{\overline{t}}_{gr}}{{\overline{t}}_{el}}$

In these equations subscripts I and II stand for fluid one and fluid two, respectively. In addition, superscript "-" stands for vector quantities and subscript $r$ stands for the ratio of the two similar quantities. $t_{el}$ and $t_{gr}$ indicate the time of propagation of an elastic disturbance, respectively, and the time needed for the media to be filled with gravity flow of the fluid one and $t_c$ is the characteristic time.

For situations in which the time of propagation of the elastic wave, compared to $t_{gr}$ , is very short, general Equation 3 could be written as follows:

${\nabla }_{xy}\left(h+\frac{{\delta }_1^2}{3K_z}{h}^2\frac{\partial h}{\partial t}+{\delta }_{2}P\right)=\frac{\partial h}{\partial t}$ (4)

${\nabla }_{xy}\left(-{\varnothing }^2+\frac{{\delta }_1{\rho }_r}{3K_z{\mu }_r{\delta }_o}{\varnothing }^2\frac{\partial \varnothing }{\partial t}+{\delta }_2{\delta }_o{\rho }_{r}P\right)=-\frac{{\rho }_r{\delta }_o}{{\mu }_r}\frac{\partial \varnothing }{\partial t}$ (5)

$\varnothing =-{\delta }_{o}h+{\delta }_o+1$

Time and considering the fact that $t_{el}\ll t_{gr}$ . Relative time of $t_r$ would be significantly less than 1. By substituting h for $\varnothing$ and eliminating P from the two equations, the system could be revised as a single equation of the coefficient h.

${\nabla }_{xy}\left(\frac{t_{el}}{{\mu }_r}\frac{\partial }{\partial t}-{\nabla }_{xy}\right)\left(h^2-\frac{2\left({\delta }_o+1\right)}{{\delta }_o+{\rho }_r}h+\frac{{\delta }_1^2{\rho }_r}{3K_z\left({\delta }_o+{\rho }_r\right)}\left[h^2+\frac{1}{{\mu }_r{\delta }_o}{\left(1+{\delta }_o-{\delta }_{o}h\right)}^2\right]\frac{\partial h}{\partial t}\right)$ $=\frac{{\rho }_r\left({\delta }_o+{\mu }_r\right)}{{\mu }_r\left({\delta }_o+{\rho }_r\right)}\frac{\partial h}{\partial t}$ (6)

By substituting ${\mu }_r\to \infty$ , ${\rho }_r\to \infty$ and ${\delta }_o\to 0$ in a problem containing one fluid with free boundary, Equation (6) could be written as follows:

${\nabla }_{xy}\left(h^2+\frac{{\delta }_1^2{\rho }_r}{3K_z\left({\delta }_o+{\rho }_r\right)}{h}^2\frac{\partial h}{\partial t}\right)=\frac{\partial h}{\partial t}$ (7)

One of the applications of Equation (7) is in cases where the height of the porous layer can be considered negligible in comparison to the horizontal length. In such a case ${\delta }_1\ll 1$ and as a result $h$ can be expressed as a series based on ${\delta }_1^2/3K_z$ powers.

In limit conditions, by ignoring the term $\left({\delta }_1^2{\rho }_r/3K_z\right)h^2\left(\partial h/\partial t\right)$ , the Equation (7) will turn into the Boussinesq equation.

${\nabla }_{xy}\left(h^2\right)=\frac{\partial h}{\partial t}$ (8)

Figure 2 Schematic of diffusion of drilling fluid.

Equation (8) is known as the classic equation of "Shallow Water Flow" theory in porous media. To model the diffusion of drilling fluid in porous media surrounding a circular well (Figure 2), Eq. (8)is given by:

$\frac{1}{r}\frac{\partial }{\partial r}\left(rh\frac{\partial h}{\partial r}\right)=\frac{\partial h}{\partial t}$ (9)

where $r=\overline{r}/r_o$ and $r_o$ is considered to be the characteristic radius. Equation (6) for porous layer with negligible thickness compared to horizontal characteristic length could be indicated as follows:

${\nabla }_{xy}\left(h^2-\frac{2\left({\delta }_o+1\right)}{\left({\delta }_o+{\rho }_r\right)}h\right)=\frac{\partial h}{\partial t}$ (10)

This equation points out the position of the interface between two immiscible fluids in porous media with insignificant height compared to its characteristic length, in $x$ and $y$ directions.

Figure 1 Generalized geometry of the problem in an orthogonal coordinate system.

Analyzing the equation and instability

For a one-dimensional flow, Equation (10) changes into a parabolic equation. By considering new independent variable $\eta =x/\sqrt{t}$ , following ordinary differential equation (ODE) could be substituted:

$\frac{d^2}{d{\eta }^2}\left(h^2-\frac{2\left({\delta }_o+1\right)}{{\delta }_o+{\rho }_r}h\right)=-\frac{h}{2}\frac{{\rho }_r\left({\delta }_o+{\mu }_r\right)}{{\mu }_r\left({\delta }_o+{\rho }_r\right)}\frac{\partial h}{\partial \eta }$ (11)

The first point about Equation (11) is attained by investigating the linear stability. By considering $h=\overline{h}+\epsilon {\delta }_h$ where $\overline{h}$ a non-disturbance response of the equation is, it could simply be given as:

$2\left(h-\gamma \right){\delta }_h^{"}=-\frac{1}{2}\eta {\delta }_h^{\text{'}}$ (12)

$H^{″}=-\frac{\eta }{4\left(1-\gamma \right)}{\delta }_h^{\text{'}}$ (13)

The domain of disturbance and is equal to $\left({\delta }_o+1\right)/\left({\delta }_o+{\rho }_r\right)$ . Therefore, to limit the amount of ${\delta }_h$ for infinite time period, $\gamma$ should possess measures below one, in other words, ${\rho }_{II}$ should be less than ${\rho }_I$ , ( ${\rho }_{II}\le {\rho }_I$ ). As a result, if the upper fluid is heavier than the lower fluid, instability phenomenon occurs.

Similarity solution of boussinesq equation

By introducing the similarity variable $\eta =r/t\frac{1+\alpha }{4}$ and relation $h=t\frac{\alpha -1}{2}f(\eta )N^2$ ,⁵ Equation (9) could be written as follows:

$\left(\frac{\alpha +1}{4}\right){\eta }^2\frac{df}{d\eta }-\frac{\left(\alpha -1\right)}{2}f+\frac{1}{\eta }\frac{d}{d\eta }\left(\eta f\frac{df}{d\eta }\right)$ (14)

A power-law flux has been considered at the origin where the power is related to the parameter $\alpha$ [5]:

$q=2\pi rh{\left(\frac{\partial h}{\partial r}\right)}_{r=0}=2\pi t^{\alpha -1}{\left[\eta f{f}_{\eta }\right]}_{\eta =0}{N}^4$ (15)

Dussan¹² disclosed that state $\alpha =1.2$ indicates filtrate invasion state in special physical problem. $\alpha =1$ Discusses the problem with constant flux, while $\alpha =2$ is linearly dependent to the increase in flux with respect to time. With respect to front position $\left(r_o\right)$ , $N$ is given by:⁵

$N=\frac{r_o}{t^{\frac{\alpha +1}{4}}}$ (16)

In this paper the "K-method" has been applied to solve the Equation (14) by adding and removing the term $\left(\alpha +1/2\right)\eta f$ to this equation. The new form of equation would be obtained as follows:²¹

$\left(\frac{\alpha +1}{4}\right){({\eta }^{2}f)}^{\text{'}}-{\left(\eta f{f}^{\prime }\right)}^{\text{'}}=\alpha \eta f$ (17)

${\left[\left(\eta f\right)f^{\prime }+\frac{\alpha +1}{4}\eta \left(\eta f\right)\right]}_1^{\eta }=\underset{1}{\overset{\eta }{{\displaystyle \int }}}\alpha \eta fd\eta$ (18)

And as for ${\eta }_1=1$ , f equals 0:

$\{\begin{array}{c}{f}^{\text{'}}+\left(\frac{\alpha +1}{4}\right)\eta =\frac{K}{\eta f}{{\displaystyle \int }}^{\text{}}\alpha \left(\eta f\right)d\eta \\ K\underset{}{\to }1 \end{array}$ (19)

Where $K$ is the K-method parameter and possesses measures between 0 and 1. For $K\to 1$ , response of the Equation (19) and (15) would be the same. With regard to the series $1/\eta f$ , it could be indicated that:

$\frac{1}{\eta f}=\frac{1}{\eta f_o}-K\frac{f_1}{\eta f_o^2}+K^2\frac{f_1^2-f_o{f}_2}{\eta f_o^3}+\dots$ $\frac{K\alpha }{\eta f}\underset{1}{\overset{\eta }{{\displaystyle \int }}}\eta fd\eta =K\left(\frac{1}{\eta f_o}\right)\underset{1}{\overset{\eta }{{\displaystyle \int }}}\alpha \eta f_{o}d\eta +K^2\left[\frac{1}{\eta f_o}\underset{1}{\overset{\eta }{{\displaystyle \int }}}\alpha \eta f_{1}d\eta -\frac{f_1}{\eta f_o^2}\underset{1}{\overset{\eta }{{\displaystyle \int }}}\alpha \eta f_{o}d\eta \right]+\dots$ (20)

or

$\{\begin{array}{l}{f}_o^{\text{'}}+\left(\frac{\alpha +1}{4}\right)\eta =0 f_1\left(1\right)=0 \\ f_o\left(1\right)=0f_1^{\text{'}}=\left(\frac{1}{\eta f_o}\right)\underset{1}{\overset{\eta }{{\displaystyle \int }}}\alpha \eta f_{o}d\eta f_1\left(1\right)=0 \\ f_2^{\text{'}}=\frac{1}{\eta f_o}\underset{1}{\overset{\eta }{{\displaystyle \int }}}\alpha \eta f_{1}d\eta -\frac{f_1}{\eta f_o^2}\underset{1}{\overset{\eta }{{\displaystyle \int }}}\alpha \eta f_{o}d\eta f_2\left(1\right)=0 \end{array}$ (21)

By implementing simple mathematical operations, results of the Equation (21) could be defined as three following statements:

$f_o\left(\eta \right)=-\frac{1}{8}\left(1+\alpha \right)\left({\eta }^2-1\right)$

$f_1\left(\eta \right)=\frac{1}{8}\alpha \left({\eta }^2-1\right)-\frac{\alpha }{4}\log x$ (22)

$f_2\left(\eta \right)=\frac{{\pi }^2{\alpha }^2}{24\left(1+\alpha \right)}$ $+{\alpha }^2\left[\left(-\log x-\frac{1}{2}\log {\left(x\right)}^2+\log x\log \left(1+x\right)-poly\log \left(2,2x\right)+poly\log \left(2,-x\right)\right)\right]/\left(2+2\alpha \right)$

Regarding the Equation (15) the amount of flux accumulation could be written as follows:

$\underset{0}{\overset{t}{{\displaystyle \int }}}\alpha dt=C{t}^{\alpha }$ (23)

Where C is the constant standing for flux accumulation and equals:

$C=2\pi N^4\underset{0}{\overset{1}{{\displaystyle \int }}}\eta fd\eta$ (24)

Thus, characteristic radius $r_o$ and interface between the two fluids could be indicated with respect to $C$ as follows:

$r_o^4=\frac{C{t}^{\alpha +1}}{2\pi {{\displaystyle \int }}_0^1\eta fd\eta }$ (25)

$h=\left(\frac{C}{2\pi {{\displaystyle \int }}_0^1\eta d\eta }\right)t^{\alpha }f\left(\eta \right)$ (26)

Considering the importance of the amount of $\underset{0}{\overset{1}{{\displaystyle \int }}}\eta fd\eta$ in assessing variables of the problem and relations obtained from Equation (21), this integral could be solved with respect to $\alpha$ .

$\underset{0}{\overset{1}{{\displaystyle \int }}}\eta fd\eta =\frac{1}{32}+\alpha \frac{1+\left(4-{\pi }^2/3\right)}{16\left(1+\alpha \right)}$ (27)

This study generally represents the function $f\left(\eta \right)/\left(1+\alpha \right)$ as the solution of Equation (14). Figure 3 compares the obtained results with both the numerical results and the results obtained by Li et al¹⁹ in $\alpha =1$ . Also, in Figure 4, the $f\left(\eta \right)/\left(1+\alpha \right)$ is presented with regards to the various values of $\eta$ and $\alpha$ . In view of the significance of amount ${{\displaystyle \int }}^{\text{}}fd{\eta }^2/1+\alpha$ in evaluating the position of the interface, Figure 5 compares assessed amount of this variable in this study with results obtained from numerical solutions.

Figure 3 Compares the obtained results with numerical results and the results.
By19 in$\alpha =1$ .

Figure 4 The values of $f\left(\eta \right)/\left(1+\alpha \right)$ respect with to the various values of η and α.

Figure 5 Comparing of the obtained values of $\left(\underset{0}{\overset{1}{{\displaystyle \int }}}fd{\eta }^2 \right)/\left(1+\text{α}\right)$ with results obtained from numerical solutions.

General model for interface between two immiscible fluids in porous media especially for case that elastic perturbation are propagating very faster than gravity perturbation have been investigated. Stability analysis of the Boussinesq equation as governing equation in this case has been done. The result showed that gravity is unstable if the lower fluid is lighter. Also the Boussinesq equation for dynamic movement of interface between two immiscible fluids has been solved by the K-Method. With a Power-law flux at the origin and dimensional analysis, the drilling mud invasion into a permeable aquifer problems becomes a nonlinear boundary value problem that have been investigated in this paper. The obtained results show that the K-Method has a good ability for investigation of B-L problems such as this solved problem.

None.

The author declares that there is no conflict of interest.

1. Calugaru CI, Calugaru DG, Crolet JM, et al. Evolution of fluid-fluid interface in porous media as the model of gas-oil fields. Electronic Journal of Differential Equations. 2003;2003(72):1–13.
2. Whitham. Linear and Nonlinear Waves. New York: John Wiley& Sons; 1974. 638 p.
3. Dagan G. second-order theory of shallow free surface flow in porous media. QJ Mech Appl Math. 1967;20(4):517–526.
4. Bear J. Dynamics of Fluid in Porous Media. American Elsevier New York: American Elsevier; 1972. 764 p.
5. Barenblatt GI, Entov VM, Ryzhik VM. Theory of Fluid Flows through Natural Rocks. Dordrecht: Kluwer Academic Publishers; 1990. 395 p.
6. Childs E. Drainage of groundwater resting on a sloping bed. Water Res journal. 1971;7(5):1256–1263.
7. Aravin VI, numerov SN. Theory of fluid flow in under formable porous media. Isr Program for Science Translation Jerusalem; 1965. 511 p.
8. Verhoest Niko EC, Troch PA. Some analytical solutions of the linearized Boussinesq equation with recharge for a sloping aquifer. Water Resources research. 2000;36(3):793–800.
9. Pistiner A. Radial Oil–Water Transport toward an Extraction Well, Transp. Porous Med. 2018;124(2):479–493.
10. Teo HT, Jeng DS, Seymour BR, Barry DA. A new analytical solution for water table fluctuations in coastal aquifers with sloping beaches. Advances In Water Resources. 2003;26(12):1239–1247.
11. Gomes AFC, Marinho JLG, Oliveira LMTM. Numeric study of a drilling fluid leak in a rock formation, permeability aspects. Brazilian Journal of Petroleum and Gas. 2016;10(4):221–232.
12. Dussan E, Auzerias FM. Buoyancy-induced flow in porous media generated near a drilled oil well. The accumulation of filtrate at a horizontal impermeable boundary. J Fluid Mech. 1993;254:283–311.
13. Doll HG. Filtrate invasion in highly permeable sands. USA: Petrol Eng; 1955. 115 p.
14. Tang Y, Jiang Q, Zhou C. Approximate analytical solution to the Boussinesq equation with a sloping water-land boundary. Water Resour. Res. 2016;52(4):2529–2550.
15. Bartlett MS, Porporato A. A Class of Exact Solutions of the Boussinesq Equation for Horizontal and Sloping Aquifers. Water Resources Research. 2018;54(2):767–778.
16. Lu D, Shen J, Cheng Y. The approximate solutions of nonlinear Boussinesq equation. Journal of Physics: Conference Series. 2016;11(2):1–9.
17. Bansal RK. Approximate Analytical Solution of Boussinesq Equation in Homogeneous Medium with Leaky Base. Applications and Applied Mathematics: An International Journal. 2016;11(1):184–191.
18. Telyakovskiy AS, Kurita S, Allen MB. Polynomial-based approximate solutions to the Boussinesq equation near a well. Advances in Water Resources. 2016;96:68–73.
19. Li L, Lockington DA, Parlange MB, et al. Similarity solution of axisymmetric flow in porous media. Advances in water resources. 2005;28(10):1076–1082.
20. Mortensen J, Sayaka O, Jean-Yves P, et al. Approximate similarity solution to a nonlinear diffusion equation with spherical symmetry. International Journal of Numerical Analysis & Modeling. 2012;9(1):105–114.
21. Shahnazari MR. A novel homotopy perturbation method: Koroush method for a thermal boundary layer in a saturated porous medium. International Journal of Engineering (IJE). 2012;25(1):59–64.