Papers in Physics, vol. 1, art. 010007 (2009)

Received: 14 October 2009, Accepted: 28 December 2009
Edited by: A. C. Mart́ı
Licence: Creative Commons Attribution 3.0
DOI: 10.4279/PIP.010007

www.papersinphysics.org

ISSN 1852-4249

Parametric study of the interface behavior between two immiscible
liquids flowing through a porous medium

Alejandro David Mariotti,1∗ Elena Brandaleze,2† Gustavo C. Buscaglia3‡

When two immiscible liquids that coexist inside a porous medium are drained through
an opening, a complex flow takes place in which the interface between the liquids moves,
tilts and bends. The interface profiles depend on the physical properties of the liquids
and on the velocity at which they are extracted. If the drainage flow rate, the liquids
volume fraction in the drainage flow and the physical properties of the liquids are known,
the interface angle in the immediate vicinity of the outlet (θ) can be determined. In this
work, we define four nondimensional parameters that rule the fluid dynamical problem
and, by means of a numerical parametric analysis, an equation to predict θ is developed.
The equation is verified through several numerical assessments in which the parameters
are modified simultaneously and arbitrarily. In addition, the qualitative influence of each
nondimensional parameter on the interface shape is reported.

I. Introduction

The fluid dynamics of the flow of two immiscible
liquids through a porous medium plays a key role
in several engineering processes. Usually, though
the interest is focused on the extraction of one of
the liquids, the simultaneous extraction of both liq-
uids is necessary. This is the case of oil production
and of ironmaking. The water injection method
used in oil production consists of injecting water
back into the reservoir, usually to increase pressure

∗E-mail: mariotti.david@gmail.com
†E-mail: ebrandaleze@frsn.utn.edu.ar
‡E-mail: gustavo.buscaglia@icmc.usp.br

1 Instituto Balseiro, 8400 San Carlos de Bariloche, Ar-
gentina.

2 Departamento de Metalurgia, Universidad Tecnológica
Nacional Facultad Regional San Nicolás, 2900 San
Nicolás, Argentina.

3 Instituto de Ciências Matemáticas e de Computação,
Universidade de São Paulo, 13560-970 São Carlos, Brasil.

and thereby stimulate production. Normally, just
a small percentage of the oil in a reservoir can be
extracted, but water injection increases that per-
centage and maintains the production rate of the
reservoir over a longer period of time. The water
displaces the oil from the reservoir and pushes it to-
wards an oil production well [1]. In the steel indus-
try, this multiphase phenomenon occurs inside the
blast furnace hearth, in which the porous medium
consists of coke particles. The slag and pig iron
are stratified in the hearth and, periodically, they
are drained through a lateral orifice. The under-
standing of this flow is crucial for the proper design
and management of the blast furnace hearth [2]. In
both examples above, when the liquids are drained,
a complex flow takes place in which the interface
between the liquids moves, tilts and bends.

Numerical simulation of multiphase flows in
porous media is focused mainly in upscaling meth-
ods, aimed at solving for large scale features of in-
terest in such a way as to model the effect of the
small scale features [3–5]. Other authors [6–8] use

010007-1



Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.

the numerical methods to model the complex mul-
tiphase flow that takes place at the pore scale.

In this work, we numerically study the macro-
scopic behavior of the interface between two im-
miscible liquids flowing through a porous medium
when they are drained through an opening. The
effect of gravity on this phenomenon is considered.
We define four nondimensional parameters that
rule the fluid dynamical problem and, by means
of a numerical parametric analysis, an equation to
predict the interface tilt in the vicinity of the orifice
(θ) is developed. The equation is verified through
several numerical cases where the parameters are
varied simultaneously and arbitrarily. In addition,
the qualitative influence of each non-dimensional
parameter on the interface shape is reported.

II. Parametric Study

The numerical studies in this work were carried out
by means of the program FLUENT 6.3.26. Differ-
ent models to simulate the two-dimensional para-
metric study were used. The volume of fluid (VOF)
method was chosen to treat the interface problem
[9]. The drag force in the porous medium was
modeled by means of the source term suggested by
Forchheimer [10]. The source term for the ith mo-
mentum equation is:

Si = −
(

µ

α
Vi +

1
2
ρC|

−→
V |Vi

)
. (1)

For the constants α and C in Eq. (1), we use the
values proposed by Ergun [10]:

α =
ε3d2

150(1 − ε)2
, (2)

C = 1.75
(1 − ε)

ε3d
, (3)

where ε is the porosity, d is the particle equivalent
diameter, ρ is the density, V is the velocity and µ
is the dynamic molecular viscosity.

Considering that the subscript 1 and 2 represent
the fluid 1 and the fluid 2 respectively, three nondi-
mensional parameters were considered in the para-
metric study: viscosity ratio, µR = µ1/µ2, density

Figure 1: Sketch of the numerical 2D domain.

ratio, ρR = ρ1/ρ2, and nondimensional velocity,
VR = V0ρ2L/µ2; where V0 is the outlet velocity
and L a reference length.

i. Domain description

The numerical domain considered to carry out the
parametric study was a two-dimensional one com-
posed by the porous medium sub-domain and the
outlet sub-domain. The porous sub-domain is a
rectangle 10m wide and 10m tall. Inside of it, a
rigid, isotropic and homogeneous porous medium
was arranged. We use a porosity and particle di-
ameter of 0.32 and 0.006m, respectively.

For the outlet domain we use a rectangle 0.02m
wide and with a height L = 0.01m divided into
two equal parts and located at the center of one
of the lateral edges. The part located at the end
of the outlet domain is used to impose the outlet
velocity. Figure 1 shows a complete description of
the domain.

A quadrilateral mesh with 2.2 × 104 cells was
used, where the outlet sub-domain mesh consists
of 200 elements in all the cases studied.

As boundary conditions, on edge 1 we define a
zero gauge pressure condition normal to the bound-
ary and impose that only the fluid 1 can enter to the
domain through it. On edge 2 the boundary condi-
tion is the same as on edge 1 but the fluid consider
in this case is fluid 2. On edge 7 we impose a zero
gauge pressure normal to the boundary but in this

010007-2



Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.

Figure 2: Interface evolution when the initial posi-
tion is below the outlet level, without gravity.

case the fluids can only leave the domain. On the
other edges (edges 4, 5, 6, and 8) we impose a wall
condition where the normal and tangential velocity
is zero except for edge 3, at which the tangential
velocity is free and the stress tangential to the edge
is zero.

ii. Interface evolution

To illustrate how an interface reaches the stationary
position from an initially horizontal one, three sets
of curves were obtained.

Figure 2 shows the interface evolution for the
case without the gravity effect and the interface
initial position is below the outlet level. The inter-
face modifies its tilt to reach the exit and it changes
its shape to reach the stationary profile.

Figures 3 and 4 show the interface evolution
when the gravity is present but the interface initial
position is below and above the outlet level respec-
tively.

iii. Viscosity effect

One of the most important parameters to modify is
the dynamical viscosity of fluid 1. We maintain the
properties of the fluid 2 as the properties of water
(density 998Kg/m3, and dynamical viscosity 0.001
Pa.s) and the density of fluid 1 as the density of the
oil (850 Kg/m3). The dynamical viscosity of fluid
1 was varied from values smaller than those of fluid
2, to values much greater.

Figure 3: Interface evolution when the initial posi-
tion is below the outlet, with gravity.

Figure 4: Interface evolution when the initial posi-
tion is above the outlet, with gravity.

Two sets of curves were obtained, one consid-
ering the effect of gravity and the other without
considering it. Figure 5 shows the stationary in-
terface profiles for the different values of viscosity,
without gravity. A value of VR = 1.5 × 104 and
ρR = 1.17 were chosen. It is possible to observe
that, when fluid 1 has a viscosity higher than that
of fluid 2, the interface profile is above the outlet
and points downwards at the outlet. If fluid 1 has
a lower viscosity, the opposite happens.

When considering gravity, the value of VR was
changed to 1 × 105 (V0 = 10m/s), since for smaller
values the interface may not reach the outlet (this is

010007-3



Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.

Figure 5: Stationary interface profiles modifying
the fluid 1 viscosity without the gravity effect.

Figure 6: Stationary interface profiles modifying
the fluid 1 viscosity with the gravity effect.

later studied in Fig. 10). Figure 6 shows the curves
obtained for this situation, where the interface only
lies over the outlet level for the higher µR values.

iv. Outlet velocity effect

V0 is varied from a small value, similar to the porous
medium velocity (V0 = 0.2m/s or VR = 2000), to a
very large one (V0 = 50m/s or VR = 5×105). Main-
taining the properties of fluid 2 similar to those of
water, two sets of curves were obtained (µR > 1
and µR < 1), shown in Figs. 7 and 8, respectively.

When the effect of gravity was considered, two
additional sets of curves (Figs. 9 and 10) were ob-

Figure 7: Stationary interface profiles for several
values of VR, without gravity, for µR > 1 (µR =
35).

Figure 8: Stationary interface profiles for several
values of VR, without gravity, for µR < 1 (µR =
0.01).

tained.
Figure 7 shows the effect of VR when the viscosity

of fluid 1 is greater than that of fluid 2, without
gravity. It is possible to see that, as VR increases,
the interface tilt at the outlet is maximal for VR =
1 × 105.

On the other hand, Fig. 8 shows the interface
profiles when the viscosity of fluid 1 is smaller than
that of fluid 2. We observe that as VR increases the
interface tends to the horizontal position.

Figure 9 shows the different stationary interface

010007-4



Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.

Figure 9: Stationary interface profiles for several
values of VR, with gravity, for µR > 1 (µR = 35).

Figure 10: Stationary interface profiles for several
values of VR, with gravity, for µR < 1 (µR = 0.01).

positions when the gravity effect is present for µR >
1. The effect of gravity is quite significant, the
interface ascends but only for the highest value of
VR it lies above the outlet level.

Figure 10 shows the curves when the viscosity of
fluid 1 is lower than that of fluid 2 (µR < 1). The
behavior is different from that without gravity. In
fact, there exists a minimum outlet velocity below
which the interface does not reach the outlet.

v. Density effect

The effect of the density ratio on the interface pro-
file was studied in the presence of gravity. Keeping

Figure 11: Stationary interface profiles for several
values of the density of fluid 1, with VR = 1.54 and
µR = 35.

fluid 2 with the properties of water and the viscos-
ity of fluid 1 as 0.035Pa.s (µR = 35), the density
of fluid 1 was varied from its original value to one
three times smaller than that of fluid 2. In Fig. 11
it is seen that as ρR increases, the interface profile
ascends significantly, with a less significant change
in the tilt angle at the outlet.

III. Generic expression

From the study on the influence of each nondi-
mensional parameter on the interface behavior, an
equation that predicts the interface angle at the im-
mediate vicinity of the outlet (θ) was crafted. For
practical reasons, the cases where gravity is present
were considered to develop the equation.

In Sect. 2.2, it is possible to see that when the
nondimensional parameters ρR, µR and VR are con-
stant the interface changes its shape until it reaches
a stationary profile.

For this reason, a fourth nondimensional param-
eter is considered, the volume fraction of fluid 1 in
the outlet flow (VF).

A generic expression [(Eq. (4)], consisting of
three terms and containing 22 constants, was ad-
justed by trial and error until satisfactory agree-
ment with the numerical results was found.

010007-5



Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.

a -29.1 i 1.162 × 107 p 0.1
b -0.04 j -0.18 q 0.72
c 0.1 k -1.64 r 2763.1
d 0.45 l -1 s 0.4
e 206 m 0.63 t -0.6
f 0.035 n 12.74 u -1.82
g -0.05 o 0.085 v 3.2
h -1.26

Table 1: Constant values in the generic expression.

PT D ε 1/α C Resistance

A 0.005 0.3 10.88 × 107 1.81 × 104 Very high
B 0.006 0.32 5.88 × 107 1.21 × 104 High
C 0.02 0.2 3 × 107 1.75 × 104 Medium
D 0.02 0.25 1.35 × 107 8400 Low
E 0.05 0.17 8.41 × 106 1.18 × 104 Very low

Table 2: Porous medium types.

θ = αµbRV
c
Rρ

d
R

+ eµfRV
g
Rρ

h
R exp(−iµ

j
RV

k
R ρ

l
RV F

m)

+ nµoRV
p
Rρ

q
R exp(−rµ

s
RV

t
Rρ

u
RV F

v) (4)

Table 1 shows the values of the constants in the
generic expression.

i. Equation verification

To verify that the generic expression (4) predicts
the value of θ correctly when the parameters are
arbitrarily modified, 21 additional numerical cases
were simulated. These cases cover a wide range of
physical properties of the liquids and of the char-
acteristics of the porous medium (by means of the
coefficients 1/α and C).

The interface angle obtained from the generic
expression (θGE ) was compared with the interface
angle obtained from the simulations (θS ). Table 2
shows the five porous medium types (PT) that were
chosen.

Some cases (1, 2, 4, 5, 9, 11, 17-21) were chosen
based on the possible combinations of immiscible
liquids that can be manipulated in real situations.
For the remaining cases, the properties of the liq-
uids were fixed at arbitrary values (fictitious liq-
uids) so that a wide range of the nondimensional

Case Fluid 1 Fluid 2 µ1 µ2 ρ1 ρ2
1 Heavy

oil
Water
solu-
tion

0.4 0.005 850 998

2 Light
oil

Water
emul-
sion

0.012 0.06 850 998

3 – – 0.08 0.008 4000 4680

4 Kero-
sene

Water 24 × 10−4 0.001 780 998

5 Ace-
tone

Water 3.3 × 10−4 0.001 791 998

6 – – 0.003 0.01 400 720
7 – – 0.2 0.01 1500 2700
8 – – 0.3 0.004 3300 5940
9 Light

slag
Hot
pig
iron

0.02 0.001 2800 7000

10 – – 0.5 0.013 500 1250
11 Medium pig 0.4 0.005 2800 7000

slag iron
12-16 – – 0.08 0.008 4000 4680
17-21 Heavy

slag
pig
iron

0.4 0.005 2800 7000

Table 3: Cases description.

Case PT ρR µR VR VF

1 B 1.17 80 5 × 104 2.8
2 B 1.17 0.2 5 × 104 87.2
3 B 1.17 10 10 × 104 28.4
4 B 1.28 2.4 5 × 104 96.1
5 B 1.26 0.33 1 × 105 96.8
6 B 1.8 0.3 1 × 105 93.2
7 B 1.8 20 1 × 105 16.1
8 B 1.8 80 8 × 104 18.4
9 B 2.5 20 1 × 105 100.0
10 B 2.5 40 1 × 105 5.5
11 B 2.5 80 5 × 104 70.8
12 A 1.2 10.0 1 × 105 23.1
13 B 1.2 10.0 1 × 105 28.4
14 C 1.2 10.0 1 × 105 41.5
15 D 1.2 10.0 1 × 105 53.4
16 E 1.2 10.0 1 × 105 63.7
17 A 2.5 80.0 5 × 104 15.8
18 B 2.5 80.0 5 × 104 24.3
19 C 2.5 80.0 5 × 104 43.2
20 D 2.5 80.0 5 × 104 70.8
21 E 2.5 80.0 5 × 104 94.0

Table 4: Parameter values and PT for all cases
described in Table 3.

parameters was covered. Table 3 shows the descrip-
tion of each numerical case, while Table 4 shows
the corresponding nondimensional parameter val-
ues and PT.

010007-6



Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.

Figure 12: Comparison between the predictions of
the generic expression and the numerical result for
the 21 validation cases.

We define an error (e = 100|∆θ|/180) as the
percentage of the absolute value of the difference
between the interface angles (∆θ = θS − θGE ) di-
vided by the interface angle range (180◦). Figure
12 shows the comparison between the generic ex-
pression and the numerical cases. It is seen that
the generic expression (4) predicts the interface an-
gle, for the cases used in this study, with an error
smaller than 10%.

IV. Conclusions

A numerical study of the macroscopic interface
behavior between two immiscible liquids flowing
through a porous medium, when they are drained
through an opening, has been reported. Four
nondimensional parameters that rule the fluid-
dynamical problem were identified. Thereby, a nu-
merical parametric analysis was developed where
the qualitative observation of the resulting inter-
face profiles contributes to the understanding of the
effect of each parameter. In addition, a generic ex-
pression to predict the interface angle in the imme-
diate vicinity of the outlet opening (θ) was devel-
oped. To verify that the generic equation predicts
the value of θ correctly, 21 numerical cases with
widely different parameters were simulated. Con-
sidering that the cases encompass a large class of
liquids and porous media, the prediction of θ within
an error of 10% is considered satisfactory.

Acknowledgements - A. D. M. and E. B. are
grateful for the support from Metallurgical Depart-
ment and DEYTEMA (UTNFSRN). G. C. B. ac-
knowledges partial financial support from CNPq
and FAPESP (Brazil).

[1] W C Lyons, G J Plisga, Standard Handbook of
Petroleum & Natural Gas Engineering - 2nd
ed, Gulf Professional Publishing, Burlington
(2005).

[2] A K Biswas, Principles of blast furnace iron-
making: theory and practice, Cootha Publish-
ing House, Brisbane (1981).

[3] A Westhead, Upscaling for two-phase flows in
porous media, PhD thesis: California Institute
of Technology, Pasadena, California (2005).

[4] R E Ewing, The Mathematics of reservoir sim-
ulation, SIAM, Philadelphia (1983).

[5] M A Cardoso, L J Durlofsky, Linearized
reduced-order models for subsurface flow sim-
ulation, J. Comput. Phys. 229, 681 (2010).

[6] M J Blunt, Flow in porous media - pore-
network models and multiphase flow, Curr.
Opin. Colloid Interface Sci. 6, 197 (2001).

[7] Z Chen, G Huan, Y Ma, Computational meth-
ods for multiphase flows in porous media,
SIAM, Philadelphia (2006).

[8] Y Efendiev, T Houb, Multiscale finite element
methods for porous media flows and their ap-
plications, Appl. Num. Math. 57, 577 (2007).

[9] FLUENT 6.3 User’s Guide, Fluent Inc. (2006).

[10] J Bear, Dynamics of fluids in porous media,
Dover Publications Inc., New York (1988).

010007-7