# Traveling wave solutions for non-Newtonian foam flow in porous media

Weslley da Silva Pereira  Department of Mathematics, Federal University of
Juiz de Fora, Rua José Lourenço Kelmer, s/n - São Pedro, Juiz de Fora,
36036900, MG, BR  Department of Mathematical and Statistical Sciences,
University of Colorado Denver, 1201 Larimer Street, Denver, 80204, CO, US
<EMAIL_ADDRESS>Department of Mathematics, Federal University of
Juiz de Fora, Rua José Lourenço Kelmer, s/n - São Pedro, Juiz de Fora,
36036900, MG, BR  Department of Mathematical and Statistical Sciences,
University of Colorado Denver, 1201 Larimer Street, Denver, 80204, CO, US
Grigori Chapiro  Department of Mathematics, Federal University of Juiz de
Fora, Rua José Lourenço Kelmer, s/n - São Pedro, Juiz de Fora, 36036900, MG,
BR<EMAIL_ADDRESS>Department of Mathematics, Federal University of Juiz
de Fora, Rua José Lourenço Kelmer, s/n - São Pedro, Juiz de Fora, 36036900,
MG, BR

###### Abstract

The injection and in-situ generation of foam in porous media successfully
control gas mobility and improve the fluids’ sweep efficiency inside porous
media. Mathematical models describing this problem use two phases, foamed gas
and fluid, and usually have a term for foam generation and destruction.
Moreover, the non-Newtonian foam behavior is frequently modeled using the
Hirasaki and Lawson’s formula for foamed gas viscosity. In this paper, we
detail how the traveling wave analysis can be used to estimate the propagation
profiles and velocity for a range of non-Newtonian foam models in porous media
at constant total superficial flow velocity. We reformulate Hirasaki and
Lawson’s formula in an explicit form allowing us to find traveling wave
solutions for the non-Newtonian Linear Kinetic model. Comparing the solution
with the one for the Newtonian version, allows us to analyze qualitatively and
quantitatively the rheology of the foam flow in porous media.

### Funding

This research was carried out in association with the ongoing R&D project
registered as ANP 20715-9, “Modelagem matemática e computacional de injeção de
espuma usada em recuperação avançada de petróleo” (Universidade Federal de
Juiz de Fora (UFJF) / Shell Brasil / ANP) - Mathematical and computational
modeling of foam injection as an enhanced oil recovery technique applied to
Brazil pre-salt reservoirs, sponsored by Shell Brasil under the ANP R&D levy
as “Compromisso de Investimentos com Pesquisa e Desenvolvimento”. This project
is carried out in partnership with Petrobras. The author G. Chapiro was
supported in part by CNPq grant 303245/2019-0.

##  1 Introduction

One technique to control gas mobility and improve fluid sweep efficiency in a
porous medium consists of injecting foaming agents, e.g., surfactants and
nanoparticles, in aqueous solutions that create bubbles in the gas phase.
Several studies point out the favorable application of this method in enhanced
oil recovery, e.g., (Casteel and Djabbarah, 1988, Grassia et al., 2014), acid
diversion during matrix stimulation, e.g., (Behenna, 1995), and contaminated
aquifer remediation, e.g., (Burman and Hall, 1986).

The behavior of foam is frequently reported as shear-thinning (Fernø et al.,
2016, Heller and Kuntamukkula, 1987, Marsden and Khan, 1966, Zitha, 2006),
which means that its apparent viscosity decreases with increased shear stress.
Hirasaki and Lawson (1985) derived a shear-thinning expression for the
apparent viscosity of foam flowing in capillary tubes and validated it through
laboratory experiments. This non-Newtonian formula is consistent with the
classical result of Bretherton (1961) using isolated bubbles and has
application in several foam propagation models, e.g., (Izadi and Kam, 2019,
Kam, 2008, Kovscek et al., 1995, Simjoo and Zitha, 2015). Hirasaki and
Lawson’s expression implicitly defines the apparent gas viscosity as it
depends on the gas velocity, and the gas velocity depends on the viscosity.
Other works report that foam behaves as a Newtonian fluid for particular
regimes. For example, foam flow can be classified as Newtonian under nearly-
constant capillary pressure during steady radial flow (Rossen, 1991).
Following this idea, some models consider foam flow to be a Newtonian fluid
(Ashoori et al., 2011, Lozano et al., 2021) and present results compatible
with experimental data (Zavala et al., 2021). Hirasaki and Lawson’s formula is
so important that numerous papers investigate the foam flow in porous media
only focusing on the apparent viscosity (or, equivalently, the mobility
reduction factor).

The studies reported above motivate the current discussion on when to consider
non-Newtonian behavior of the foam flow. Vassenden and Holt (1998) present and
validate a model for the transition of foam flow behavior from Newtonian to
shear-thinning according to the gas flow rate. Alvarez et al. (2001) conclude
the foam behavior is shear-thinning in the low-quality regime and shear-
thickening in high-quality regimes. Rossen (1991) highlights the change from
Newtonian to non-Newtonian behavior of uniform texture foam for changing
capillary pressure.

One may find several predictive mathematical models in the literature aiming
to represent foam propagation inside porous media. Local Equilibrium (LE) foam
models use algebraic expressions to compute the foam texture. However, LE
models can be inadequate when strong foam generation is not certain (Kam et
al., 2007). Mechanistic foam models, introduced by Falls et al. (1988), Patzek
(1988), describe the foam texture dynamically using partial differential
equations (PDE). They have the unique potential to describe both transient and
stationary conditions for foam flow as they track the foam texture in space
and time (Eide et al., 2020). Many works report mechanistic models to match
laboratory experiments successfully, e.g., in (Kovscek et al., 1995, Kam,
2008, Simjoo and Zitha, 2015). Besides that, some authors prefer these models
when extracting analytical estimates of the foam flow, e.g., (Almajid et al.,
2019, Eide et al., 2020, Simjoo and Zitha, 2015). Moreover, several works
suggest that, after a transient stage, foam travels inside the porous media
under invariant water saturation profiles, e.g., (de Vries and Wit, 1990,
Ettinger and Radke, 1992). One strategy to obtain solution profiles, employed
in (Izadi and Kam, 2019, Kam, 2008, Simjoo and Zitha, 2015), consists of
disregarding the capillary pressure gradient in each phase’s superficial
velocity and applying the classical method of characteristics. Another
strategy consists of seeking traveling waves solutions, which are profiles
invariant to translation in space (Volpert et al., 2000). They are especially
suited for analyzing nonlinear wave propagation problems at constant speed.
Ashoori et al. (2011), Lozano et al. (2021) used this procedure to find
traveling wave solutions for Newtonian foam models at fixed total superficial
velocity. Ashoori et al. (2011) also find traveling waves for the model from
Kam (2008) which does not consider capillary pressure derivatives.

The present work shows how to compute traveling wave solutions that connect
two equilibrium states of a generic mechanistic non-Newtonian foam model that
consider capillary pressure derivatives. We assume the total superficial
velocity is fixed, which is the case of predominant rectilinear flow. The non-
Newtonian foamed gas viscosity is modeled by the Hirasaki and Lawson’s
formula. We rewrite this formula in terms of water saturation, foam texture
and total superficial velocity to avoid the direct dependency on the gas
superficial velocity. The traveling wave profiles serve as basis to compare a
linear kinetic foam model using Newtonian and non-Newtonian gas viscosity
expressions.

We organize this paper as follows. Section 2 presents the general population-
based foam model. Section 3 shows the Hirasaki and Lawson apparent gas
viscosity and its reformulated version. Section 4 presents the procedure for
obtaining traveling wave solutions for the non-Newtonian foam flow model.
Section 5 applies the procedure to a linear kinetics model, and use the
traveling wave profiles to compare Newtonian and non-Newtonian models. Section
6 presents the discussions and conclusions.

##  2 The population-balance foam model

Consider the two-phase flow of gas-water solution with foaming agents in a
porous medium, where the flow occurs only in one direction, and all fluids are
incompressible. The last assumption is undesirable but necessary for the
mathematical analysis performed in this paper. The process is modeled by the
mass balance of each phase $i$ as

$\displaystyle\frac{\partial(\phi\,S_{i})}{\partial t}+\frac{\partial
u_{i}}{\partial x}=0\,,\qquad i\in\\{w,g\\}\,,$ (1)

where $t$ is the time variable, $x$ is the space variable, $\phi$ is the
porosity, $S_{i}$ and $u_{i}$ are the saturation and superficial velocity of
the phase $i$. The subscript ‘$w$’ represents the aqueous phase and ‘$g$’
represents the gas phase. The phase superficial velocities satisfy the
generalized Darcy’s laws (Chen and Ewing, 1997)

$\displaystyle u_{\text{w}}$
$\displaystyle=f_{\text{w}}\,\left(u+\lambda_{\text{g}}\,K\,\frac{\partial
P_{\text{c}}}{\partial x}\right)\,,\qquad
u_{\text{g}}=f_{\text{g}}\,\left(u-\lambda_{\text{w}}\,K\,\frac{\partial
P_{\text{c}}}{\partial x}\right)\,,$ (2)

with fractional fluxes $f_{i}$ and phase mobilities $\lambda_{\text{i}}$ given
by

$\displaystyle f_{i}$
$\displaystyle:=\frac{\lambda_{\text{i}}}{\lambda_{\text{w}}+\lambda_{\text{g}}}\,,\qquad\lambda_{\text{i}}:=\frac{k_{\text{ri}}}{\mu_{\text{i}}}\,,\qquad
i\in\\{w,g\\},$ (3)

where $K$ is the permeability, $u=u_{w}+u_{g}$ is the (constant) total
superficial velocity, $P_{\text{c}}$ is the capillary pressure, and for each
phase $i$, $k_{\text{ri}}$ is the relative permeability, and $\mu_{\text{i}}$
is the apparent phase viscosity. The balance of foam texture $n_{\text{f}}$ in
the gas phase is

$\frac{\partial(\phi\,S_{\text{g}}\,n_{\text{f}})}{\partial
t}+\frac{\partial(n_{\text{f}}\,u_{\text{g}})}{\partial
x}=\phi\,S_{\text{g}}\,R\,,$ (4)

where $R:=r_{\text{g}}-r_{\text{c}}$, and $r_{\text{g}}$ and $r_{\text{c}}$
are the rates of generation and coalescence of foam, respectively. It is
useful to define the dimensionless foam texture ($n_{\text{D}}$) as follows

$n_{\text{D}}:=\frac{n_{\text{f}}}{n_{\max}}\,,$ (5)

where $n_{\max}$ is the reference foam texture.

Some auxiliary physical quantities are equally important in this work. We
define the effective water saturation ($S_{\text{we}}$), with values in [0,1],
as

$\displaystyle
S_{\text{we}}:=\frac{S_{\text{w}}-S_{\text{wc}}}{1-S_{\text{wc}}-S_{\text{gr}}}\,,$
(6)

where $S_{\text{wc}}$ is the connate water saturation, and $S_{\text{gr}}$ is
the residual gas saturation. The total apparent viscosity ($\mu_{app}$) is

$\displaystyle\mu_{app}:=\frac{1}{\lambda_{\text{w}}+\lambda_{\text{g}}}\,,$
(7)

and the mobility reduction factor (MRF) is

$\displaystyle MRF:=\frac{\mu_{\text{g}}}{\mu_{\text{g}}^{0}}\,.$ (8)

To avoid problems in the nomenclature, we emphasize using the word “total”
when referring to the apparent viscosity $\mu_{app}$. Note that the apparent
total viscosity represents the equivalent viscosity of the fluid system in the
following Darcy’s law

$\displaystyle u=-\frac{K}{\mu_{app}}\nabla P\,,$ (9)

where $P$ is the global pressure (Antontsev, 1972). Finally, the phase-$i$
interstitial velocity is

$\displaystyle v_{i}:=\frac{|u_{i}|}{\phi\,S_{i}}\,,\qquad i\in\\{w,g\\}.$
(10)

It is worth mentioning that several models exist that state $P_{\text{c}}$
depending exclusively on the fluid phase saturation, and many of them are
based on the classical works of Corey (1954), Brooks and Corey (1966), van
Genuchten (1980). The relative permeabilities $k_{\text{rw}}$ and
$k_{\text{rg}}$ are either derived directly from the capillary pressure or
fitted experimentally (see (Li and Horne, 2006) and references therein). Foam
generation and coalescence expressions $r_{\text{g}}$ and $r_{\text{c}}$ are
either inspired on microscopic mechanisms (Kam, 2008, Kovscek et al., 1995),
based on macroscopic observations (Ashoori et al., 2011, Simjoo and Zitha,
2015), or fitted using laboratory experiments (Thorat and Bruining, 2016). We
postpone the definitions of the expressions $k_{\text{rw}}(S_{\text{w}})$,
$k_{\text{rg}}(S_{\text{w}})$,
$r_{\text{g}}(S_{\text{w}},n_{\text{D}},u_{\text{w}})$,
$r_{\text{c}}(S_{\text{w}},n_{\text{D}},u_{\text{w}})$, and
$P_{\text{c}}(S_{\text{w}})$ to Section 5. Section 3 presents the formulas for
the apparent gas viscosity ($\mu_{\text{g}}$).

Some mechanistic (population-based) models use an expression for LE foam
texture, $n_{\text{D}}^{\text{LE}}$, that naturally arises from the
restriction $R=0$. For example:

  * •

The model from Simjoo and Zitha (2015) uses a constant value
$n_{\text{D}}^{\text{LE}}$, which is compatible with laboratory experiments.

  * •

In (Ashoori et al., 2011), $R=0$ if and only if
$n_{\text{D}}=n_{\text{D}}^{\text{LE}}(S_{\text{w}})=\tanh(400\,(S_{\text{w}}-0.37))$,
which approximates the Heaviside step function.

  * •

In (Chen et al., 2010), $n_{\text{D}}^{\text{LE}}$ is the only real solution
of

$\displaystyle(n_{\text{D}}^{\text{LE}})^{3}+\frac{k_{-1}v_{\text{g}}^{2/3}}{k_{1}^{0}v_{\text{w}}}n_{\text{D}}^{\text{LE}}-1=0\,,$
(11)

where $k_{-1}$ and $k_{1}^{0}$ are positive coefficients possibly depending on
$S_{\text{w}}$.

In the examples above, $0\leq n_{\text{D}}^{\text{LE}}\leq 1$. The expression
$n_{\text{D}}^{\text{LE}}$ from Kovscek et al. (1995) can be also written as a
function of $S_{\text{w}}$, $v_{\text{w}}$ and $v_{\text{g}}$ but it is not
constrained to the interval $[0,1]$. In (Kam, 2008),
$n_{\text{D}}^{\text{LE}}$ is also not constrained to the interval $[0,1]$ and
it is a function of $S_{\text{w}}$ and $\nabla P$ (See (Ashoori et al., 2011)
for a traveling wave analysis of the model in (Kam, 2008)).

######  Remark.

Some models make a distinction between trapped and flowing foam phases, e.g.,
(Izadi and Kam, 2019, Kovscek et al., 1995, Kharabaf and Yortsos, 1998). Other
studies report the importance of understanding the gas trapping mechanism to
model foam flow in porous media (Almajid et al., 2019, Jones et al., 2018).
Simjoo and Zitha (2015) observe no trapping, arguing that foam could be moving
slowly. On top of that, it is still not clear how to model trapped foam
fraction, although some models try to explain it (Kovscek et al., 1995, Cohen
et al., 1997, Tang and Kovscek, 2006). The analysis presented in this paper
may be applied to some trapped foam models, for instance, the model from
Kovscek et al. (1995) and its derivatives.

##  3 The Hirasaki and Lawson’s gas viscosity

The apparent gas viscosity in the presence of foam is frequently modeled by
the Hirasaki and Lawson’s formula

$\mu_{\text{g}}=\mu_{\text{g}}^{0}+\alpha\,\frac{n_{\text{f}}}{v_{\text{g}}^{1/3}}\,,$
(12)

when foam is considered a shear-thinning fluid, where $\mu_{\text{g}}^{0}$ is
the gas viscosity in the absence of foam. The proportionality constant
$\alpha$ depends on many factors, including the pore structure, the liquid
viscosity, and the gas-liquid surface tension (Hirasaki and Lawson, 1985). Due
to the inherent difficulty to obtain $\alpha$, many works use it as a fitting
parameter, e.g., (Chen et al., 2010, Eide et al., 2020, Kam, 2008, Kovscek et
al., 1995, Simjoo and Zitha, 2015, Thorat and Bruining, 2016). Since
$v_{\text{g}}$ itself depends on $\mu_{\text{g}}$ (see 2, 3 and 10), Eq. 12
defines $\mu_{\text{g}}$ implicitly, which hinders the analytical and
numerical analysis of the model.

The first contribution of this work is to rewrite 12 eliminating the direct
dependency on $v_{\text{g}}$. To do so, we replace (2), (3), and (10) into
(12) obtaining a cubic equation in the variable
$X=(k_{\text{rg}}+\lambda_{\text{w}}\,\mu_{\text{g}})^{1/3}$. The resulting
expression admits a unique positive root $X$, which can be written as

$\displaystyle\mu_{\text{g}}$
$\displaystyle=\mu_{\text{g}}^{0}+\frac{3\,A}{\lambda_{w}}\,\left(\sqrt[3]{B+\sqrt{B^{2}-A^{3}}}+\sqrt[3]{B-\sqrt{B^{2}-A^{3}}}\right)\,,$
(13)

where

$\displaystyle A$
$\displaystyle:=\frac{\alpha\,n_{\text{f}}\,\lambda_{\text{w}}}{3}\,\left(\frac{\phi\,(1-S_{\text{w}})}{k_{\text{rg}}\,\left|u-\lambda_{\text{w}}\,K\,\partial_{x}P_{\text{c}}\right|}\right)^{1/3}\quad\text{
and }\quad
B:=\frac{k_{\text{rg}}+\lambda_{\text{w}}\,\mu_{\text{g}}^{0}}{2}\,.$ (14)

Under LE conditions, the capillary pressure is constant. Therefore, term $A$
in (13) becomes

$\displaystyle A$
$\displaystyle:=\frac{\alpha\,n_{\max}\,n_{\text{D}}^{\text{LE}}\,\lambda_{\text{w}}}{3}\,\left(\frac{\phi\,(1-S_{\text{w}})}{k_{\text{rg}}\,|u|}\right)^{1/3}.$
(15)

Since $\lambda_{\text{w}}$ and $k_{\text{rg}}$ are functions of
$S_{\text{w}}$, and $\mu_{\text{g}}^{0}$, $\alpha$, $n_{\max}$, and $\phi$ are
constant parameters, we identify the apparent gas viscosity under LE
conditions as the function
$\mu_{\text{g}}^{LE}(S_{\text{w}},n_{\text{D}}^{\text{LE}},u)$.

##  4 Traveling wave solutions

Ashoori et al. (2011), Lozano et al. (2021) investigated traveling wave
solutions for a linear kinetic Newtonian foam model, where the apparent gas
viscosity ($\mu_{\text{g}}$) depends only on the foam texture
($n_{\text{D}}$). Below, we present the analysis for the generic non-Newtonian
model described in Section 2 with apparent viscosity formula detailed in
Section 3.

###  4.1 Traveling wave formulation

The traveling wave solutions of a system of evolutionary PDEs are the
solutions invariant to translation in space (see (Volpert et al., 2000) and
references therein). Mathematically speaking, they depend on a traveling
variable $\xi:=x-vt$, where $v$ is the constant wave velocity. For the system
(1), (2) and (4) with step function initial conditions

$(S_{\text{w}},n_{\text{D}})=\left\\{\begin{array}[]{cc}(S_{\text{w}}^{-},n_{\text{D}}^{-}),&x=0\\\\[5.0pt]
(S_{\text{w}}^{+},n_{\text{D}}^{+}),&x>0\end{array}\right.$ (16)

the traveling wave solution is represented by bounded and differentiable
functions
$\tilde{S}_{\text{w}},\tilde{n}_{\text{D}}:\mathbb{R}\rightarrow\mathbb{R}$
satisfying
$(\tilde{S}_{\text{w}},\tilde{n}_{\text{D}})(x-vt)=(S_{\text{w}},n_{\text{D}})(x,t)$
for a some $v\in\mathbb{R}$, all $x\in(0,L)$, $L>0$, and all
$t\in(0,+\infty)$. Moreover, the wave front
$(\tilde{S}_{\text{w}},\tilde{n}_{\text{D}})$ must satisfy the boundary
conditions (Volpert et al., 2000):

$\lim_{\xi\rightarrow\pm\infty}(\tilde{S}_{\text{w}},\tilde{n}_{\text{D}})(\xi)=(S_{\text{w}}^{\pm},n_{\text{D}}^{\pm})\in\mathbb{R}^{2}\,,\quad\text{and}\quad\lim_{\xi\rightarrow\pm\infty}d_{\xi}(\tilde{S}_{\text{w}},\tilde{n}_{\text{D}})(\xi)=0\,.$
(17)

The (left) state $(S_{\text{w}}^{-},n_{\text{D}}^{-})$ corresponds to the
inflow boundary conditions, and the (right) state
$(S_{\text{w}}^{+},n_{\text{D}}^{+})$ corresponds to the initial conditions
for the original PDE.

Starting from (1), (2) and (4), one may verify that the wave front
$(\tilde{S}_{\text{w}},\tilde{n}_{\text{D}})$ is a solution of the following
system of ordinary differential equations (ODEs):

$\left\\{\begin{aligned}
&d_{\xi}(-v\,\phi\,\tilde{S}_{\text{w}}+\tilde{u}_{\text{w}})=0\,,\\\
&d_{\xi}\left[\tilde{n}_{\text{D}}\,n_{\max}\,(-v\,\phi\,(1-\tilde{S}_{\text{w}})+(u-\tilde{u}_{\text{w}}))\right]=\phi\,(1-\tilde{S}_{\text{w}})\,\tilde{R},\end{aligned}\right.$
(18)

for all $\xi\in\mathbb{R}$, where
$\tilde{R}:=R\left(\,\tilde{S}_{\text{w}}\,,\;\tilde{n}_{\text{D}}\,,\;\tilde{u}_{\text{w}}\,\right)$,

$\displaystyle\tilde{u}_{\text{w}}=\frac{\lambda_{\text{w}}(\tilde{S}_{\text{w}})}{\lambda_{\text{w}}(\tilde{S}_{\text{w}})+\frac{k_{\text{rg}}(\tilde{S}_{\text{w}})}{\tilde{\mu}_{\text{g}}}}\,\left(u+\frac{k_{\text{rg}}(\tilde{S}_{\text{w}})}{\tilde{\mu}_{\text{g}}}\,K\,\frac{\partial
P_{\text{c}}}{\partial
S_{w}}(\tilde{S}_{\text{w}})\,\frac{\partial\tilde{S}_{\text{w}}}{\partial\xi}\right)\,,$
(19)

and $\tilde{\mu}_{\text{g}}$ is given below in 23. The solution
$(\tilde{S}_{\text{w}},\tilde{n}_{\text{D}})$ of (18) satisfying restrictions
(17) is known as traveling wave solution. We can use the limits in (17) to
conclude that there exists
$u_{\text{g}}^{+}:=u-\lim_{\xi\rightarrow+\infty}\tilde{u}_{\text{w}}$ such
that

$\displaystyle-v\,\phi\,(1-\tilde{S}_{\text{w}})+(u-\tilde{u}_{\text{w}})=-v\,\phi\,(1-S_{\text{w}}^{+})+u_{\text{g}}^{+}$
(20)

for all $\xi\in\mathbb{R}$. Thus, we can rewrite 18-19 as follows

$\displaystyle\left\\{\begin{aligned}
&\frac{-\lambda_{\text{w}}\,k_{\text{rg}}}{k_{\text{rg}}+\lambda_{\text{w}}\,\tilde{\mu}_{\text{g}}}\,K\,\frac{\partial
P_{\text{c}}}{\partial
S_{w}}\,d_{\xi}\tilde{S}_{\text{w}}=\left[u_{\text{g}}^{+}-\frac{k_{\text{rg}}\,u}{k_{\text{rg}}+\lambda_{\text{w}}\,\tilde{\mu}_{\text{g}}}\right]-v\,\phi\,(\tilde{S}_{\text{w}}-S_{\text{w}}^{+})\,,\\\
&n_{\max}\,(u_{\text{g}}^{+}-v\,\phi\,(1-S_{\text{w}}^{+}))\,d_{\xi}\tilde{n}_{\text{D}}=\phi\,(1-\tilde{S}_{\text{w}})\,\tilde{R}\,,\end{aligned}\right.$
(21)

Note that 20 also implies that, for all $\xi\in\mathbb{R}$,

$\widetilde{v_{\text{g}}}:=\frac{|u-\tilde{u}_{\text{w}}|}{\phi\,(1-\tilde{S}_{\text{w}})}=\left|v\frac{\left(S_{\text{w}}^{+}+\frac{u_{\text{g}}^{+}}{\phi
v}\right)-\tilde{S}_{\text{w}}}{1-\tilde{S}_{\text{w}}}\right|\,.$ (22)

Then, one may replace (22) into (12) to obtain

$\tilde{\mu}_{\text{g}}=\mu_{\text{g}}^{0}\left(1+\frac{\alpha\,n_{\max}}{\mu_{\text{g}}^{0}\,|v|^{\frac{1}{3}}}\,\left(\frac{1-\tilde{S}_{\text{w}}}{\big{|}S_{\text{w}}^{crit}-\tilde{S}_{\text{w}}\big{|}}\right)^{\frac{1}{3}}\,\tilde{n}_{\text{D}}\right)\,,$
(23)

where $S_{\text{w}}^{crit}:=S_{\text{w}}^{+}+u_{\text{g}}^{+}/(v\,\phi)$.
Equation (23) is an alternative version of (12) valid only for the traveling
wave profile.

Let us summarize these calculations. The traveling wave solution of the
problem (1)-(4) in the domain $(0,L)\times(0,+\infty)$, with initial condition
$(S_{\text{w}},n_{\text{D}})(x,0)=(S_{\text{w}}^{+},n_{\text{D}}^{+})$ and
boundary condition
$(S_{\text{w}},n_{\text{D}})(0,t)=(S_{\text{w}}^{-},n_{\text{D}}^{-})$ is
given by the solution of 21 which satisfies limit conditions (17). We discuss
how to obtain these limit conditions in Subsection 4.2.

The existence of the solution of (21) is not trivial. Lozano et al. (2021)
investigate the existence of solutions for 21 using Newtonian foam flow and
particular choice of $R=r_{\text{g}}-r_{\text{c}}$. They show that such
solutions exist in specific parameter regions, depending on the model and the
limit states. We do not perform such deep analysis here; however, we give
evidence that the solutions exist for the examples from Section 5.

###  4.2 Traveling wave velocity and equilibria

The usual procedure to determine the traveling wave velocity and equilibria is
to apply limits 17 to the ODE system 21. Let us assume
$k_{\text{rw}}(S_{\text{w}}^{\pm})$, $k_{\text{rg}}(S_{\text{w}}^{\pm})$, and
$d_{S_{\text{w}}}P_{\text{c}}(S_{\text{w}}^{\pm})$ are real numbers and that
$k_{\text{rg}}(S_{\text{w}}^{\pm})>0$. As
$\lim_{\xi\rightarrow\pm\infty}d_{\xi}\tilde{S}_{\text{w}}(\xi)=0$ from (17),
it follows that $\lim_{\xi\rightarrow\pm\infty}d_{x}P_{\text{c}}(\xi)=0$.
Therefore,
$\lim_{\xi\rightarrow\pm\infty}\tilde{\mu}_{\text{g}}(\xi)=\mu_{\text{g}}^{LE}(S_{\text{w}}^{\pm},n_{\text{D}}^{\pm},u)$.
Moreover, we use the fractional flux definition in 3, and apply the limits in
17 to 21 and 19 to obtain the algebraic system of nonlinear equations

$\displaystyle\frac{k_{\text{rg}}(S_{\text{w}}^{-})}{k_{\text{rg}}(S_{\text{w}}^{-})+\lambda_{\text{w}}(S_{\text{w}}^{-})\,\mu_{\text{g}}^{LE}(S_{\text{w}}^{-},n_{\text{D}}^{-},u)}=f_{\text{g}}^{-}\,,$
(24)
$\displaystyle\frac{k_{\text{rg}}(S_{\text{w}}^{+})}{k_{\text{rg}}(S_{\text{w}}^{+})+\lambda_{\text{w}}(S_{\text{w}}^{+})\,\mu_{\text{g}}^{LE}(S_{\text{w}}^{+},n_{\text{D}}^{+},u)}=\frac{u_{\text{g}}^{+}}{u}\,,$
(25)
$\displaystyle\left(u_{\text{g}}^{+}-f_{\text{g}}^{-}\,u\right)-v\,\phi\,(S_{\text{w}}^{-}-S_{\text{w}}^{+})=0\,,$
(26) $\displaystyle
R\left(\,S_{\text{w}}^{-}\,,\;n_{\text{D}}^{-}\,,\;(u-u_{\text{g}}^{+})+v\,\phi\,(S_{\text{w}}^{-}-S_{\text{w}}^{+})\,\right)=0\,,$
(27) $\displaystyle
R\left(\,S_{\text{w}}^{+}\,,\;n_{\text{D}}^{+}\,,\;u-u_{\text{g}}^{+}\,\right)=0\,.$
(28)

Most laboratory and field scale simulations use homogeneous initial saturation
profiles, and measure the injection gas fractional flow to predict foam
propagation. Therefore, we consider $u$, $S_{\text{w}}^{+}$ and
$f_{\text{g}}^{-}$ to be the known values and
$(S_{\text{w}}^{-},n_{\text{D}}^{-},n_{\text{D}}^{+},u_{\text{g}}^{+},v)$ as
the main unknowns in 24-28.

As discussed in Section 2, there is a particular case where the restrictions
$R=0$ and $\partial_{x}P_{\text{c}}=0$ result in a LE foam texture
$n_{\text{D}}=n_{\text{D}}^{\text{LE}}(S_{\text{w}},u_{\text{w}},u)$. The
models from (Ashoori et al., 2011, Chen et al., 2010, Kam, 2008, Kovscek et
al., 1995, Simjoo and Zitha, 2015), for instance, allow for such an expression
of local-equilibrium foam texture. For these models, one may use the following
strategy to solve 24-28 given $u$, $S_{\text{w}}^{+}$ and $f_{\text{g}}^{-}$:

  1. 1.

Solve 24 for $S_{\text{w}}^{-}$ using
$n_{\text{D}}^{-}=n_{\text{D}}^{\text{LE}}(S_{\text{w}}^{-},(1-f_{\text{g}}^{-})\,u,u)$;

  2. 2.

Solve 25 for $u_{\text{g}}^{+}$ using
$n_{\text{D}}^{+}=n_{\text{D}}^{\text{LE}}(S_{\text{w}}^{+},u-u_{\text{g}}^{+},u)$;

  3. 3.

Solve 26 for $v$.

This algorithm involves the solution of two uncoupled scalar nonlinear
equations, which can be done either analytically or applying a classical
numerical procedure, see (Mathews, 1992). After obtaining all constant
parameters, one must solve the ODE system (21) with limit conditions in 17
using usual ODE solvers for nonlinear systems to obtain the traveling wave
profile.

######  Remark.

The procedure to obtain the traveling wave limits and wave velocity described
above can still be used for initial gas at residual saturation. In this case,
the existence of traveling wave limits rely on the following limit

$\displaystyle\lim_{\xi\rightarrow+\infty}d_{\xi}\tilde{S}_{\text{w}}=\lim_{S_{\text{w}}\rightarrow
1-S_{\text{gr}}}\frac{u-v\,\phi\,(1-S_{\text{w}}-S_{\text{gr}})\frac{k_{\text{rg}}+\lambda_{\text{w}}\,\tilde{\mu}_{\text{g}}}{k_{\text{rg}}}}{\lambda_{\text{w}}\,K\,\frac{\partial
P_{\text{c}}}{\partial S_{w}}}=0\,.$

Ashoori et al. (2011) use a modified Brooks-Corey capillary pressure equation
obtained by multiplying the original formula by the term
$(1-S_{\text{w}}-S_{\text{gr}})^{c}$, $c=0.01$. Such modification does not
affect the phenomenological conclusions, but guarantees the limit above holds.
In the present paper, we use the model from Ashoori et al. (2011) with the
referred modification. We stress that one can propose other modifications to
the capillary pressure provided the limit above remains valid.

##  5 Application in a linear kinetics model

In this section, we apply the methodology from the previous sections to obtain
traveling wave solutions for the linear kinetics foam propagation model from
Ashoori et al. (2011). We use them to compare the influence of the non-
Newtonian apparent viscosity in the foam model. In this model, the relative
permeabilities, $k_{\text{rw}}(S_{\text{w}})$ and
$k_{\text{rg}}(S_{\text{w}})$, are generalized Brooks-Corey expressions,
fitted for the nitrogen-water flow in the Boise sandstone. The capillary
pressure has an equivalent unitary pore size distribution index, also inspired
in the work of Brooks and Corey (1966). Foam generation and coalescence are
controlled by expressions $r_{\text{g}}(S_{\text{w}})$ and
$r_{\text{c}}(n_{\text{D}})$ as follows

$\displaystyle\begin{aligned}
k_{\text{rw}}(S_{\text{w}})&:=0.2\,S_{\text{we}}^{4.2}\,,\\\
k_{\text{rg}}(S_{\text{w}})&:=0.94\,(1-S_{\text{we}})^{1.3}\,,\\\
P_{\text{c}}(S_{\text{w}})&:=15000\,\frac{0.022\,(1-S_{\text{w}}-S_{\text{gr}})^{0.01}}{S_{\text{w}}-S_{\text{wc}}}\,,\\\
r_{\text{g}}(S_{\text{w}})&:=K_{\text{c}}\,n_{\max}n_{\text{D}}^{\text{LE}}(S_{\text{w}})\,,\\\
r_{\text{c}}(n_{\text{f}})&:=K_{\text{c}}\,n_{\max}\,n_{\text{D}}\,,\\\
n_{\text{D}}^{\text{LE}}(S_{\text{w}})&:=\left\\{\begin{array}[]{ll}\tanh(400\,(S_{\text{w}}-0.37))&\text{
if $S_{\text{w}}>0.37$\,,}\\\ 0&\text{ if $S_{\text{w}}\leq 0.37$\,,}\\\
\end{array}\right.\end{aligned}$ (29)

where $K_{\text{c}}$ is the kinetic generation/coalescence rate parameter. The
coefficient 15000 stems from the Leverett J-function (Leverett, 1941) with a
gas-water surface tension of 0.03 $\mathrm{N}\mathrm{/}\mathrm{m}$. Table 1
shows the remaining parameter values from (Ashoori et al., 2011). To
summarize, we analyze the model described by (1)-(4), (12) and (29) with
parameters from Table 1.

Ashoori et al. (2011) studied the incompressible foam flow considering foam as
a Newtonian fluid. They used the following apparent gas viscosity formula

$\mu_{\text{g}}^{N}:=\mu_{\text{g}}^{0}\,MRF^{N}(n_{\text{D}})\,,\qquad
MRF^{N}(n_{\text{D}})=1+C_{MRF}\,n_{\text{D}}\,,$ (30)

where $C_{MRF}$ is a reference mobility reduction factor for the strongest
foam. Two well-known foam properties inspired the simplification: (1) foamed
gas mobility reduction is large and nearly constant at high water saturation,
and (2) there exists a limiting capillary pressure, where foam collapses
abruptly. Parameter $\alpha$ is the one that minimizes the mean squared
difference between the gas mobilities under equilibrium conditions, i.e.,

$\displaystyle\min_{\alpha}\int_{S_{\text{wc}}}^{1-S_{\text{gr}}}\left(\frac{k_{\text{rg}}(S_{\text{w}})}{\mu_{\text{g}}^{LE}(S_{\text{w}},n_{\text{D}}^{\text{LE}}(S_{\text{w}}),u)}-\frac{k_{\text{rg}}(S_{\text{w}})}{\mu_{\text{g}}^{0}\,MRF^{N}(n_{\text{D}}^{\text{LE}}(S_{\text{w}}))}\right)^{2}\,dS_{\text{w}}\,.$
(31)

for $u=$29.3$$\mu\mathrm{m}\mathrm{/}\mathrm{s}$$. This mapping generalizes
the one proposed by Zavala et al. (2021).

Table 1: Parameters for nitrogen and water flow in a Boise sandstone (Ashoori et al., 2011). Property | Symbol | Value
---|---|---
Newtonian gas viscosity coefficient | $C_{MRF}$ |  $18\,500$ [-]
Permeability | $K$ |  ${10}^{-12}$ $\mathrm{m}^{2}$
Kinetic parameter | $K_{\text{c}}$ | 200; 1; 0.01 $\mathrm{s}^{-1}$
Reference foam texture | $n_{\max}$ |  $8\text{x}{10}^{13}$ $\mathrm{m}^{-3}$
Residual gas saturation | $S_{\text{gr}}$ |  $0.18$ [-]
Residual (connate) water saturation | $S_{\text{wc}}$ |  $0.2$ [-]
Total superficial velocity | $u$ |  $2.93\text{x}{10}^{-5}$ $\mathrm{m}\mathrm{/}\mathrm{s}$
Non-Newtonian gas viscosity coefficient | $\alpha$ |  $2.44\text{x}{10}^{-16}$ $\mathrm{P}\mathrm{a}\,\mathrm{s}^{\frac{2}{3}}\,\mathrm{m}^{\frac{10}{3}}$
Water viscosity | $\mu_{\text{w}}$ |  ${10}^{-3}$ $\mathrm{P}\mathrm{a}\,\mathrm{s}$
Gas viscosity in the absence of foam | $\mu_{\text{g}}^{0}$ |  $2\text{x}{10}^{-05}$ $\mathrm{P}\mathrm{a}\,\mathrm{s}$
Porosity | $\phi$ |  $0.25$ [-]

In what follows, we compute traveling wave profiles for the linear kinetic
model using the non-Newtonian apparent gas viscosity 12. We compare those
profiles with the ones obtained using the Newtonian apparent viscosity 30. We
analyze two scenarios. The first one reproduces the example found in (Ashoori
et al., 2011). The second one show a case where the non-Newtonian foam model
brings significant information. We also show how apparent total viscosity
changes with flow velocity in the equilibria, which helps understanding the
connection between Newtonian and non-Newtonian foam models. Hereafter, the
superscript ‘$N$’ identifies the properties of the Newtonian and ‘$S$’ the
properties of the non-Newtonian models. Table 2 shows the equilibria computed
using 24-28 for the Newtonian and non-Newtonian models in the two scenarios.

Table 2: Parameter values at limit states and traveling wave velocities in two scenarios used to compare Newtonian and non-Newtonian models. We use $u=29.3$\mu\mathrm{m}\mathrm{/}\mathrm{s}$$. | Scenario 1 | Scenario 2
---|---|---
|  $f_{w}^{-}$ = 0.268; $S_{w}^{+}$ = 0.72 |  $f_{w}^{-}$ = 0.9; $S_{w}^{+}$ = 0.819
Parameter | Newtonian | non-Newt. | Newtonian | non-Newt.
$S_{\text{w}}^{-}$ [-] | $0.372\,004$ | $0.371\,940$ | $0.501\,128$ | $0.467\,238$
$n_{\text{D}}^{-}$ [-] | $0.664\,985$ | $0.650\,301$ | $1.0$ | $1.0$
$n_{\text{D}}^{+}$ [-] | $1.0$ | $1.0$ | $1.0$ | $1.0$
$u_{\text{g}}^{+}$ [$\mu\mathrm{m}\mathrm{/}\mathrm{s}$] | $7.250\,802\,325\,588\,347\text{x}{10}^{-02}$ | $6.111\,156\,621\,159\,01\text{x}{10}^{-03}$ | $8.781\,610\text{x}{10}^{-05}$ | $3.192\,801\text{x}{10}^{-07}$
$|\nabla P|^{-}$ [$\mathrm{M}\mathrm{P}\mathrm{a}\mathrm{/}\mathrm{m}$] | $8.565\,506$ | $8.579\,051$ | $2.737\,615$ | $4.520\,126$
$|\nabla P|^{+}$ [$\mathrm{M}\mathrm{P}\mathrm{a}\mathrm{/}\mathrm{m}$] | $3.059\,092\text{x}{10}^{-01}$ | $3.066\,042\text{x}{10}^{-01}$ | $1.474\,962\text{x}{10}^{-01}$ | $1.474\,966\text{x}{10}^{-01}$
$v$ [$\mathrm{m}\mathrm{/}\mathrm{s}$] | $2.456\,937\text{x}{10}^{-04}$ | $2.464\,111\text{x}{10}^{-04}$ | $3.686\,905\text{x}{10}^{-05}$ | $3.331\,801\text{x}{10}^{-05}$

###  5.1 Scenario 1

Scenario 1 corresponds to the flow data from the literature (Ashoori et al.,
2011) with the limit states $S_{\text{w}}^{+}=0.72$ and
$f_{\text{w}}^{-}=0.268$. The autonomous system (21) can be rewritten in the
standard form:

$\displaystyle d_{\xi}\left[\begin{array}[]{c}\tilde{S}_{\text{w}}\\\
\tilde{n}_{\text{D}}\end{array}\right]=F(\tilde{S}_{\text{w}},\tilde{n}_{\text{D}});$
where (34) $\displaystyle F(S_{\text{w}},n_{\text{D}})$
$\displaystyle:=\left[\begin{array}[]{c}\dfrac{\left(u_{\text{g}}^{+}-\frac{k_{\text{rg}}\,u}{k_{\text{rg}}+\lambda_{\text{w}}\,\mu}\right)-v\,\phi\,(S_{\text{w}}-S_{\text{w}}^{+})}{-\frac{k_{\text{rg}}}{k_{\text{rg}}+\lambda_{\text{w}}\,\mu}\,\lambda_{\text{w}}\,K\,\frac{\partial
P_{\text{c}}}{\partial S_{w}}}\\\ \\\
\dfrac{\phi\,(1-S_{\text{w}})\,R}{n_{\max}\,(u_{\text{g}}^{+}-v\,\phi\,(1-S_{\text{w}}^{+}))}\end{array}\right]$
(38) $\displaystyle\text{and}\quad\mu$
$\displaystyle:=\mu_{\text{g}}^{0}+\frac{\alpha\,n_{\max}}{v^{\frac{1}{3}}}\,\left(\frac{1-S_{\text{w}}}{\big{|}S_{\text{w}}^{crit}-S_{\text{w}}\big{|}}\right)^{\frac{1}{3}}\,n_{\text{D}},$
(39)

One standard approach for proving the existence of the solution connecting two
equilibria in an autonomous system is based on the analysis of the ODE’s phase
space, as done by Lozano et al. (2021). We repeat some steps of their
analysis.

We start by proving that:

  1. 1.

$F(S_{\text{w}},n_{\text{D}})=0$ if and only if
$n_{\text{D}}=n_{\text{D}}^{\text{LE}}(S_{\text{w}})$ and
$u_{\text{g}}^{+}-k_{\text{rg}}/(k_{\text{rg}}+\lambda_{\text{w}}\,\mu)\,u=v\,\phi\,(S_{\text{w}}-S_{\text{w}}^{+})$;

  2. 2.

The curve
$u_{\text{g}}^{+}-k_{\text{rg}}\,u/(k_{\text{rg}}+\lambda_{\text{w}}\,\mu)$
restricted to $n_{\text{D}}=n_{\text{D}}^{\text{LE}}(S_{\text{w}})$ intersects
the straight line $v\,\phi\,(S_{\text{w}}-S_{\text{w}}^{+})$ in at most two
values of $S_{\text{w}}\in(S_{\text{wc}},1-S_{\text{gr}})$.

Since
$F(S_{\text{w}}^{+},n_{\text{D}}^{+})=F(S_{\text{w}}^{-},n_{\text{D}}^{-})=0$,
those are the only two roots of $F$ in
$(S_{\text{wc}},1-S_{\text{gr}})\times[0,1]$.

Second, let us comment on the dynamics in the neighborhood of the limit states
of Scenario 1 in Table 2. In two-dimensional phase space, this analysis
reduces to look at the real part’s sign of the eigenvalues of the Jacobian
matrix $(\partial_{S_{\text{w}}}F,\partial_{n_{\text{D}}}F)$; positive signs
indicate source directions, and negative signs indicate sink directions. Since
$F(S_{\text{w}},n_{\text{D}})=0$ only over the equilibria, the orbit (solution
of 34) exits if there is at least one source direction in the left state and
at least one sink direction in the right state. Table 3 shows the eigenvalues
for states $-$ and $+$ using the non-Newtonian model. In all cases, the states
$-$ have eigenvalues with positive real parts (sources), and states $+$ have
eigenvalues’ real parts with opposite signs (saddle points), which proves
existence of an orbit, i.e., traveling wave solution.

Table 3: Eigenvalues of the Jacobian matrix $(\partial_{S_{\text{w}}}F,\partial_{n_{\text{D}}}F)$ at the equilibrium states $-$ and $+$ of Scenario 1. The classification follows a typical nomenclature (Guckenheimer et al., 1984). $K_{\text{c}}$ [$\mathrm{1}\mathrm{/}\mathrm{s}$] | State |  $\lambda_{1}$ [$\mathrm{1}\mathrm{/}\mathrm{m}$] |  $\lambda_{2}$ [$\mathrm{1}\mathrm{/}\mathrm{m}$] | Classification
---|---|---|---|---
200 | $-$ | $108\,891.085\,247\,048\,88$ | $1\,715\,138.395\,853\,513\,4$ | Source
200 | $+$ | $-1\,019\,924.184\,787\,762\,2$ | $809\,184.986\,264\,093\,5$ | Saddle
1 | $-$ | $1091+2854\text{$\mathrm{i}$}$ | $1091-2854\text{$\mathrm{i}$}$ | Spiral source
1 | $+$ | $-1\,019\,924.184\,787\,762\,2$ | $4045.924\,931\,320\,468$ | Saddle
0.01 | $-$ | $12\,087.064\,939\,222\,36$ | $772.574\,989\,099\,824\,7$ | Source
0.01 | $+$ | $-1\,019\,924.184\,787\,762\,2$ | $40.459\,249\,313\,204\,68$ | Saddle

Table 2 shows that the relative difference for the states $n_{\text{D}}^{-}$
is lower than 3%, and the values of $u_{\text{g}}^{+}$ are approximately zero.
The remaining parameter values vary less than 1% between the models. Note that
the Newtonian model’s values in scenario 1 coincide with those from (Ashoori
et al., 2011). The solution profiles for $S_{\text{w}}$ and $n_{\text{D}}$ for
the non-Newtonian case are very similar to the ones of the Newtonian model
found in (Ashoori et al., 2011). Figure 1 shows the MRF and $f_{w}$ profiles.
Although the MRF we observe in the left panel of Fig. 1 a significant
qualitative and quantitative change in MRF, this does not affect fractional
flux curves as depicted in the right panel of the same figure.

Figure 1: Mobility Reduction Factor and $f_{w}$ profiles in the non-Newtonian
(filled lines) and Newtonian (dashed lines) models using $K_{\text{c}}=200$
$\mathrm{s}^{-1}$ in Scenario 1.

###  5.2 Scenario 2

The second scenario explores higher differences between the Newtonian and non-
Newtonian behaviors by taking a higher initial water saturation,
$S_{\text{w}}^{+}=0.819$, and injecting a mixture with lower foam quality,
$f_{\text{w}}^{-}=0.9$. This scenario, usually denoted by the term low-quality
regime, was not based on experimental data from the literature; it was built
to enlarge the differences between Newtonian and non-Newtonian regimes
describing foam flow in porous media. Scenario 2 addresses a fundamental
scientific discussion on understanding the non-Newtonian nature of the foam
displacement. The analysis of the existence of traveling wave solutions for
Scenario 2 is analogous to the one from Section 5.1 and, therefore, we omit
it.

Table 2 shows the equilibria for the Newtonian and non-Newtonian models in
Scenario 2. We note that the Newtonian traveling wave velocity $v$, water
saturation $S_{\text{w}}^{-}$, and pressure gradient $|\nabla P|^{-}$ differ
from the corresponding values of the non-Newtonian values 11%, 7%, and 39%,
respectively. Figure 2 shows the traveling wave profiles of
$\tilde{S}_{\text{w}}$ and $\tilde{n}_{\text{D}}$, and the pressure gradient
over the traveling wave for the non-Newtonian foam model using
$K_{\text{c}}=200$ $\mathrm{s}^{-1}$. Besides the difference in the limit
states, the non-Newtonian foam model profiles are sharper than the ones in the
Newtonian model. Figure 3 shows the MRF and $f_{w}$ over the traveling wave.
We use 23 to compute $MRF^{S}$. As in scenario 1, MRF is qualitatively and
quantitatively different between the models. In this case, however, the
fractional flux profile changes between the models leading to different water
saturation profiles, see left panel of Fig. 2.

Figure 2: Wave profiles in the non-Newtonian (filled lines) and Newtonian
(dashed lines) models using $K_{\text{c}}=200$ $\mathrm{s}^{-1}$ in Scenario
2. Left: $\tilde{S}_{\text{w}}(\xi)$ and $\tilde{n}_{\text{D}}(\xi)$. Right:
pressure gradient obtained using the formulas from Section 2.

Figure 3: Mobility Reduction Factor and $f_{w}$ profiles in the non-Newtonian
(filled lines) and Newtonian (dashed lines) models using $K_{\text{c}}=200$
$\mathrm{s}^{-1}$ in Scenario 2.

###  5.3 Traveling wave equilibria for different flow velocities

Figure 4 shows isolines of the superficial gas velocity $u_{\text{g}}$ as a
function of the foam quality $f_{\text{g}}$ and the total apparent viscosity
$\mu_{app}$ defined in (7). All quantities are taken considering LE conditions
and
$u\in[0.1$\mu\mathrm{m}\mathrm{/}\mathrm{s}$,100$\mu\mathrm{m}\mathrm{/}\mathrm{s}$]$,
a range that covers laboratory and field experimental values. One can see that
the lower is the superficial gas velocity the stronger is the foam (higher
$\mu_{app}$). In agreement with (Eide et al., 2020, Gassara et al., 2020), the
high-quality regime region (after the peak of the strongest foam) shrinks when
the superficial gas velocity increases. In the Newtonian model, total apparent
viscosity and fractional fluxes depend only on $S_{\text{w}}$, so there is a
single curve for each $C_{MRF}$. We plot the curves $C_{MRF}=$37\,000$$ and
$C_{MRF}=$185\,000$$ in Fig. 4 to show that one should vary $C_{MRF}$ so that
the Newtonian model better represent the non-Newtonian gas viscosity
accordingly to the flow regime.

Figure 4: Apparent viscosity $\mu_{app}$ as a function of $u_{\text{g}}$ and
$f_{\text{g}}$, all using LE foam, for the non-Newtonian foam model. The
dashed and dotted lines represent Newtonian total apparent viscosity with a
given $C_{MRF}$.

######  Remark.

Scenario 2 evidences some difference between the traveling wave solutions of
non-Newtonian and Newtonian models. The parameter values in this scenario
assume the injection of a mixture with low gas concentration. Although it is
uncommon, such configuration occurs in industrial application, e.g.,
surfactant alternated gas (SAG) injection in the Enhanced Oil Recovery (EOR).
During this procedure, pure water (containing surfactant) and pure gas are
injected into the reservoir periodically resulting in alternated low and high
gas concentrations inside the porous medium.

######  Remark.

We used the explicit Runge-Kutta method RK5(4) (Shampine, 1986) to solve (21)
in the variable $\eta:=-\xi$ on $(0,\eta_{F})$ with initial conditions
$(\tilde{S}_{\text{w}},\tilde{n}_{\text{D}})(0)=(S_{\text{w}}^{+}-10^{-6},n_{\text{D}}^{+})$
and time steps bounded by $\eta_{F}/1000$. This change of variables, $\xi$ to
$-\eta$, helps the numerical solver to find the correct solution.

##  6 Discussion and Conclusions

In this paper, we detailed the procedure of obtaining traveling wave solutions
for a range of models describing propagation of foam in porous media. For
that, we reformulate the Hirasaki and Lawson’s apparent viscosity in terms of
the foam model’s unknowns $S_{\text{w}}$, $n_{\text{D}}$ and $u$. We restrict
the analysis to the case of constant total superficial velocity $u$ as in
previous works. It is worth mentioning that the reformulated apparent
viscosity formula can be used in other contexts. One example is to reconstruct
fractional fluxes directly from the experimental data.

In order to evidence the potential of the proposed methodology, we find the
traveling wave solutions for the linear kinetic model from Ashoori et al.
(2011). We use the solutions to quantify the importance of considering the
non-Newtonian apparent viscosity in the model. We proposed a mapping procedure
to map one model in the other. This procedure involves the solution of a
minimization problem for the mean square difference between gas mobilities in
the two models. We conclude that:

  * •

For the cases reproduced from Ashoori et al. (2011), the non-Newtonian
viscosity brings no new relevant information. This is expected in such a high-
quality regime, as observed in several other works, e.g., (Cheng et al., 2000,
Salazar Castillo et al., 2020, Valdez et al., 2021).

  * •

In the low-quality example, we identify a foam flow regime with stronger non-
Newtonian fluid behavior, and obtain more significant discrepancies between
the two models. Our example showed a relative difference of 39% in the
pressure gradient and 11% in the traveling wave velocity. Note that estimating
the traveling wave velocity is crucial to predict the gas breakthrough during
an enhanced recovery process in porous media. It is worth mentioning that we
achieve such contrast even using the same constant total superficial velocity
from Scenario 1.

  * •

The mobility reduction factor gives significantly less information on the foam
flow than the fractional flux curves. In the two scenarios we analyze, the
mobility reduction factors between Newtonian and non-Newtonian presented
qualitative and quantitative differences. However, the water saturation and
foam texture wave profiles coincide qualitatively in the first scenario,
while, in the second one, they don’t.

  * •

For different total flow velocities and quality regimes, one may still use the
Newtonian gas viscosity formula, although the parameters should be adjusted so
that the shear-thinning behavior of foam can be accurately represented.

## References

  * Almajid et al. (2019) M. M. Almajid, N. Nazari, and A. R. Kovscek.  Modeling Steady-State Foam Flow: Hysteresis and Backward Front Movement.  _Energy Fuels_ , 33(11):11353–11363, 2019.  ISSN 1520-5029.  doi: 10.1021/acs.energyfuels.9b01842.
  * Alvarez et al. (2001) J. M. Alvarez, H. J. Rivas, and W. R. Rossen.  Unified model for steady-state foam behavior at high and low foam qualities.  _SPE J._ , 6(3):325–333, Sept. 2001.  ISSN 1086-055X.  doi: 10.2118/74141-PA.
  * Antontsev (1972) S. N. Antontsev.  On the solvability of boundary value problems for degenerate two-phase porous flow equations.  _Dinamika Splosnoi Sredy Vyp_ , 10:28–53, 1972.
  * Ashoori et al. (2011) E. Ashoori, D. Marchesin, and W. R. Rossen.  Roles of transient and local equilibrium foam behavior in porous media: Traveling wave.  _Colloids Surf., A_ , 377(1-3):228–242, Mar. 2011\.  ISSN 0927-7757.  doi: 10.1016/j.colsurfa.2010.12.042.
  * Behenna (1995) F. R. Behenna.  Acid Diversion from an Undamaged to a Damaged Core Using Multiple Foam Slugs.  In _SPE - European Formation Damage Conference_ , pages 441–452, 1995\.  ISBN 9781555634476.  doi: 10.2118/30121-ms.
  * Bretherton (1961) F. P. Bretherton.  The motion of long bubbles in tubes.  _J. Fluid Mech._ , 10(2):166–188, Mar. 1961.  ISSN 1469-7645.  doi: 10.1017/S0022112061000160.
  * Brooks and Corey (1966) R. H. Brooks and A. T. Corey.  Properties of Porous Media Affecting Fluid Flow.  _J. Irrig. Drain. Div._ , 92(2):61–88, 1966.  ISSN 0044-7978.
  * Burman and Hall (1986) J. W. Burman and B. E. Hall.  Foam As a Diverting Technique for Matrix Sandstone Stimulation.  In _SPE of AIME_ , Apr. 1986.  doi: 10.2118/15575-ms.
  * Casteel and Djabbarah (1988) J. F. Casteel and N. F. Djabbarah.  Sweep improvement in the CO2 flooding by use of foaming agents.  _SPE Reserv. Eng._ , 3(4):1186–1192, Nov. 1988\.  ISSN 0885-9248.  doi: 10.2118/14392-PA.
  * Chen et al. (2010) Q. Chen, M. G. Gerritsen, and A. R. Kovscek.  Modeling foam displacement with the local-equilibrium approximation: Theory and experimental verification.  _SPE J._ , 15(1):171–183, Mar. 2010.  ISSN 1086-055X.  doi: 10.2118/116735-pa.
  * Chen and Ewing (1997) Z. Chen and R. E. Ewing.  Fully discrete finite element analysis of multiphase flow in groundwater hydrology.  _SIAM J. Numer. Anal._ , 34(6):2228–2253, Dec. 1997.  ISSN 0036-1429.  doi: 10.1137/S0036142995290063.
  * Cheng et al. (2000) L. Cheng, A. B. Reme, D. Shan, D. A. Coombe, and W. R. Rossen.  Simulating Foam Processes at High and Low Foam Qualities.  In _SPE/DOE improved oil recovery symposium_. Society of Petroleum Engineers, 2000.  doi: 10.2118/59287-ms.
  * Cohen et al. (1997) D. Cohen, T. W. Patzek, and C. J. Radke.  Onset of Mobilization and the Fraction of Trapped Foam in Porous Media.  _Transp. Porous Media_ , 28(3):253–284, 1997\.  ISSN 0169-3913.  doi: 10.1023/A:1006552320036.
  * Corey (1954) A. T. Corey.  The Interrelation Between Gas and Oil Relative Permeabilites.  _Producers Monthly_ , 19(1):38–41, 1954.
  * de Vries and Wit (1990) A. S. de Vries and K. Wit.  Rheology of gas/water foam in the quality range relevant to steam foam.  _SPE Reserv. Eng._ , 5(2):185–192, 1990.  ISSN 0885-9248.  doi: 10.2118/18075-pa.
  * Eide et al. (2020) Ø. Eide, M. Fernø, S. Bryant, A. Kovscek, and J. Gauteplass.  Population-balance modeling of CO2 foam for CCUS using nanoparticles.  _J. Nat. Gas Sci. Eng._ , 80:103378, Aug. 2020.  ISSN 1875-5100.  doi: 10.1016/j.jngse.2020.103378.
  * Ettinger and Radke (1992) R. A. Ettinger and C. J. Radke.  Influence of texture on steady foam flow in Berea sandstone.  _SPE Reserv. Eng._ , 7(1):83–90, Feb. 1992.  ISSN 0885-9248.  doi: 10.2118/19688-PA.
  * Falls et al. (1988) A. H. Falls, G. J. Hirasaki, T. W. Patzek, D. A. Gauglitz, D. D. Miller, and T. Ratulowski.  Development of a mechanistic foam simulator: The population balance and generation by snap-off.  _SPE Reserv. Eng._ , 3(3):884–892, 1988.  ISSN 0885-9248.  doi: 10.2118/14961-PA.
  * Fernø et al. (2016) M. A. Fernø, G. Jarand, P. Monrawee, H. Åsmund, G. Arne, A. R. Kovscek, and G. Hirasaki.  Experimental study of foam generation, sweep efficiency, and flow in a fracture network.  _SPE J._ , 21(4):1140–1150, Aug. 2016.  ISSN 1086-055X.  doi: 10.2118/170840-pa.
  * Gassara et al. (2020) O. Gassara, F. Douarche, B. Braconnier, and B. Bourbiaux.  Calibrating and Scaling Semi-empirical Foam Flow Models for the Assessment of Foam-Based EOR Processes (in Heterogeneous Reservoirs).  _Transp. Porous Media_ , 131(1):193–221, 2020\.  ISSN 1573-1634.  doi: 10.1007/s11242-018-01223-5.
  * Grassia et al. (2014) P. Grassia, E. Mas-Hernández, N. Shokri, S. J. Cox, G. Mishuris, and W. R. Rossen.  Analysis of a model for foam improved oil recovery.  _J. Fluid Mech._ , 751:346–405, 2014.  doi: 10.1017/jfm.2014.287.
  * Guckenheimer et al. (1984) J. Guckenheimer, P. Holmes, and M. Slemrod.  Nonlinear Oscillations Dynamical Systems, and Bifurcations of Vector Fields.  _J. Appl. Mech._ , 51(4):947–947, Dec. 1984.  ISSN 0021-8936.  doi: 10.1115/1.3167759.
  * Heller and Kuntamukkula (1987) J. P. Heller and M. S. Kuntamukkula.  Critical Review of the Foam Rheology Literature.  _Ind. Eng. Chem. Res._ , 26(2):318–325, Feb. 1987\.  ISSN 1520-5045.  doi: 10.1021/ie00062a023.
  * Hirasaki and Lawson (1985) G. J. Hirasaki and J. B. Lawson.  Mechanisms of Foam Flow in Porous Media: Apparent Viscosity in Smooth Capillaries.  _SPE J._ , 25(2):176–190, 1985.  ISSN 0197-7520.  doi: 10.2118/12129-pa.
  * Izadi and Kam (2019) M. Izadi and S. I. Kam.  Bubble-population-balance modeling for supercritical carbon dioxide foam enhanced-oil-recovery processes: From pore-scale to core-scale and field-scale events.  _SPE Reserv. Eval. Eng._ , 22(4):1467–1480, Nov. 2019.  ISSN 1094-6470.  doi: 10.2118/191202-PA.
  * Jones et al. (2018) S. A. Jones, N. Getrouw, and S. Vincent-Bonnieu.  Foam flow in a model porous medium: II. the effect of trapped gas.  _Soft Matter_ , 14(18):3497–3503, 2018.  ISSN 1744-6848.  doi: 10.1039/c7sm02458d.
  * Kam (2008) S. I. Kam.  Improved mechanistic foam simulation with foam catastrophe theory.  _Colloids Surf., A_ , 318(1-3):62–77, 2008.  ISSN 0927-7757.  doi: 10.1016/j.colsurfa.2007.12.017.
  * Kam et al. (2007) S. I. Kam, Q. P. Nguyen, Q. Li, and W. R. Rossen.  Dynamic simulations with an improved model for foam generation.  _SPE J._ , 12(1):35–48, Mar. 2007.  ISSN 1086-055X.  doi: 10.2118/90938-PA.
  * Kharabaf and Yortsos (1998) H. Kharabaf and Y. C. Yortsos.  Pore Network Model for Foam Formation and Propagation in Porous Media.  _SPE J._ , 3(1):42–53, Mar. 1998.  ISSN 1086-055X.  doi: 10.2118/36663-pa.
  * Kovscek et al. (1995) A. R. Kovscek, T. W. Patzek, and C. J. Radke.  A mechanistic population balance model for transient and steady-state foam flow in Boise sandstone.  _Chem. Eng. Sci._ , 50(23):3783–3799, 1995.  ISSN 0009-2509.  doi: 10.1016/0009-2509(95)00199-F.
  * Leverett (1941) M. C. Leverett.  Capillary Behavior in Porous Solids.  _Trans. AIME_ , 142(01):152–169, 1941.  ISSN 0081-1696.  doi: 10.2118/941152-g.
  * Li and Horne (2006) K. Li and R. N. Horne.  Comparison of methods to calculate relative permeability from capillary pressure in consolidated water-wet porous media.  _Water Resour. Res._ , 42(6):1–9, June 2006.  ISSN 0043-1397.  doi: 10.1029/2005WR004482.
  * Lozano et al. (2021) L. F. Lozano, R. Q. Zavala, and G. Chapiro.  Mathematical properties of the foam flow in porous media.  _Comput. Geosci._ , 25(1):515–527, Feb. 2021\.  ISSN 1420-0597.  doi: 10.1007/s10596-020-10020-3.
  * Marsden and Khan (1966) S. S. Marsden and S. A. Khan.  The Flow of Foam Through Short Porous Media And Apparent Viscosity Measurements.  _SPE J._ , 6(01):17–25, Mar. 1966.  ISSN 0197-7520.  doi: 10.2118/1319-pa.
  * Mathews (1992) J. H. Mathews.  _Numerical Methods For Mathematics, Science, and Engineering_.  Prentice-Hall International, Englewood Cliffs, N.J, 1992.  ISBN 0136249906.  doi: 10.2307/2282127.
  * Patzek (1988) T. W. Patzek.  Description of Foam Flow in Porous Media by the Population Balance Method.  In _ACS Symposium Series_ , volume 373, pages 326–341. ACS, July 1988.  doi: 10.1021/bk-1988-0373.ch016.
  * Rossen (1991) W. R. Rossen.  Rheology of foam in porous media at the ”limiting capillary pressure.  In _6th European Symposium on IOR_. EAGE, 1991.  doi: 10.3997/2214-4609.201411236.
  * Salazar Castillo et al. (2020) R. O. Salazar Castillo, S. F. Ter Haar, C. G. Ponners, M. Bos, and W. Rossen.  Fractional-Flow Theory for Non-Newtonian Surfactant-Alternating-Gas Foam Processes.  _Transp. Porous Media_ , 131(2):399–426, jan 2020\.  ISSN 0169-3913.  doi: 10.1007/s11242-019-01351-6.
  * Shampine (1986) L. F. Shampine.  Some practical Runge-Kutta formulas.  _Math. Comp._ , 46(173):135–135, Jan. 1986.  ISSN 0025-5718.  doi: 10.1090/s0025-5718-1986-0815836-3.
  * Simjoo and Zitha (2015) M. Simjoo and P. L. J. Zitha.  Modeling of Foam Flow Using Stochastic Bubble Population Model and Experimental Validation.  _Transp. Porous Media_ , 107(3):799–820, Apr. 2015.  ISSN 0169-3913.  doi: 10.1007/s11242-015-0468-y.
  * Tang and Kovscek (2006) G. Q. Tang and A. R. Kovscek.  Trapped gas fraction during steady-state foam flow.  _Transp. Porous Media_ , 65(2):287–307, Nov. 2006\.  ISSN 0169-3913.  doi: 10.1007/s11242-005-6093-4.
  * Thorat and Bruining (2016) R. Thorat and H. Bruining.  Foam Flow Experiments. I. Estimation of the Bubble generation-Coalescence Function.  _Transp. Porous Media_ , 112(1):53–76, 2016.  ISSN 1573-1634.  doi: 10.1007/s11242-016-0632-z.
  * Valdez et al. (2021) A. R. Valdez, B. M. Rocha, J. M. F. Façanha, A. V. O. de Souza, A. Perez-Gramatges, G. Chapiro, and R. W. Santos.  Foam-assisted water–gas flow parameters: From core-flood experiment to uncertainty quantification and sensitivity analysis.  _Transp. Porous Media_ , pages 1–21, 2021.
  * van Genuchten (1980) M. T. van Genuchten.  A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils.  _Soil Sci. Soc. Am. J._ , 44(5):892–898, 1980\.  ISSN 1435-0661.  doi: 10.2136/sssaj1980.03615995004400050002x.
  * Vassenden and Holt (1998) F. Vassenden and T. Holt.  Experimental Foundation for Relative Permeability Modeling of Foam.  In _SPE/DOE IOR Symposium_. SPE, 1998.  doi: 10.2118/39660-ms.
  * Volpert et al. (2000) A. I. Volpert, V. A. Volpert, and V. A. Volpert.  _Translations of Mathematical Monographs- Traveling Wave Solutions of Parabolic Systems_ , volume 140.  AMS, 2000.  ISBN 0821897578.
  * Zavala et al. (2021) R. Q. Zavala, L. F. Lozano, G. Chapiro, and P. L. J. Zitha.  Analytical solution for the population-balance model describing foam displacement.  _Transport in Porous Media_ , -:1–17, 2021.  doi: 10.1007/s11242-021-01589-z.
  * Zitha (2006) P. L. J. Zitha.  A new stochastic bubble population model for foam in porous media.  In _Proceedings - SPE Symposium on IOR_ , volume 1, pages 116–129. SPE, 2006.  doi: 10.2118/98976-ms.


# Universal Free-Fall Law for Liquid Jets under Fully Developed Injection
Conditions

M. Beneitez  DAMTP, Centre for Mathematical Sciences, Wilberforce Road,
Cambridge CB3 0WA, United Kingdom  D. Moreno-Boza  Grupo de Mecánica de
Fluidos, Departamento de Ingeniería Térmica y de Fluidos, Universidad Carlos
III de Madrid. Avda. de la Universidad 30, 28911, Leganés, Madrid, Spain.  A.
Sevilla  Corresponding author<EMAIL_ADDRESS>Grupo de Mecánica de
Fluidos, Departamento de Ingeniería Térmica y de Fluidos, Universidad Carlos
III de Madrid. Avda. de la Universidad 30, 28911, Leganés, Madrid, Spain.

###### Abstract

We show that vertical slender jets of liquid injected in air with a fully-
developed outlet velocity profile have a universal shape in the common case in
which the viscous force is much smaller than the gravitational force. The
theory of ideal flows with vorticity provides an analytical solution that,
under negligible surface tension forces, predicts
$R_{j}(Z)=[(1+Z/4)^{1/2}-(Z/4)^{1/2}]^{1/2}$, where $R_{j}$ is the jet radius
scaled with the injector radius and $Z$ is the vertical distance scaled with
the gravitational length, $l_{g}=u_{o}^{2}/2g$, where $u_{o}$ is the mean
velocity at the injector outlet and $g$ is the gravitational acceleration. In
contrast with Mariotte’s law, $R_{j}=(1+Z)^{-1/4}$, previously reported
experiments employing long injectors collapse almost perfectly onto the new
solution.

††preprint: Revised version submitted to Phys. Rev. Lett.

_Introduction_.– The structure, dynamics and rupture of freely falling liquid
jets in air have been the subject of scientific inquiry for centuries [1, 2].
Although early insightful observations appear in Da Vinci’s Codex Leicester
[3], the first quantitative theory of falling liquid columns is due to
Mariotte [4], who made use of Torricelli’s free-fall law
$u_{j}(z)=(u_{o}^{2}+2gz)^{1/2}$, where $u_{j}(z)$ is the axial velocity of
the liquid, assumed to be uniform in the radial direction, $u_{o}$ is the
injection velocity, $z$ is the vertical distance measured downwards from the
release height, and $g$ is the gravitational acceleration [5]. When combined
with the liquid continuity equation, $u_{j}(z)r_{j}(z)^{2}=u_{o}r_{o}^{2}$,
where $r_{o}$ and $r_{j}(z)$ are the injector and jet radii, respectively, the
dimensionless jet radius $R_{j}(z)=r_{j}(z)/r_{o}$ obeys Mariotte’s law,
$R_{j}=(1+2gz/u_{o}^{2})^{-1/4}$. The convective acceleration of the liquid
due to gravity, $\bm{v}\cdot\bm{\nabla}\bm{v}\sim u_{o}^{2}/l_{g}\sim g$
defines a gravitational length $l_{g}=u_{o}^{2}/2g$ such that

$R_{j}(Z)=(1+Z)^{-1/4},$ (1)

in terms of the rescaled height $Z=z/l_{g}=z/(Fr\,r_{o})$, where
$Fr=l_{g}/r_{o}=u_{o}^{2}/(2gr_{o})$ is the Froude number.

Mariotte’s law (1) provides a highly accurate description of the asymptotic
structure of free liquid jets for $Z\gg 1$, a limit in which all the forces
other than liquid inertia and gravity become negligibly small [6]. These
additional forces include, among others, surface tension [7], viscous stresses
[8, 9], aerodynamic forces [10, 11], and the radial diffusion of liquid
momentum that takes place when the initial velocity profile is non-uniform
[12, 13, 14, 15, 16, 17]. Here we focus on the latter viscous relaxation
effect, which is particularly important in actual engineering flows due to the
fact that most liquid injection devices involve the development of viscous
boundary layers at the inner injector wall, with the fully developed
Poiseuille profile as the canonical example of non-uniform exit conditions.
The downstream diffusion of the outlet vorticity in the free-jet region
depends on the Reynolds number, $Re=u_{o}r_{o}/\nu$, where $\nu$ is the
kinematic viscosity of the liquid. When $Re\lesssim 1$ the velocity profile
becomes uniform within a distance $l_{\nu}\sim r_{o}$ from the outlet. In
contrast, when $Re\gg 1$, the viscous development length $l_{\nu}\sim
u_{o}r_{o}^{2}/\nu=Re\,r_{o}\gg r_{o}$, as deduced by balancing the convective
acceleration, $\bm{v}\cdot\bm{\nabla}\bm{v}\sim u_{o}^{2}/l_{\nu}$, with the
radial diffusion of momentum, $\nu\bm{\nabla}^{2}\bm{v}\sim\nu
u_{o}/r_{o}^{2}$. The structure of the relaxing free jet is thus controlled by
the ratio $l_{\nu}/l_{g}\sim 2gr_{o}^{2}/(\nu u_{o})=G=Re/Fr$. The gravity
number $G$ allows to distinguish two relevant limiting cases: when $G\ll 1$,
the viscous relaxation takes place under almost buoyancy-free conditions [12,
14], and the gravitational thinning occurs far downstream, where the velocity
profile is almost uniform. In the opposite limit $G\gg 1$, the viscous
diffusion of momentum is negligible along the initial gravitational thinning
region.

Here we provide a first-principles description of the jet structure in the
frequently realized limit $G\gg 1$ assuming, for simplicity, the slenderness
hypothesis $Fr\gg 1$. To that end, advantage is taken of the fact that viscous
forces can be neglected, in a first approximation, within the gravitational
stretching region $0<z\lesssim l_{g}\ll l_{\nu}$. Under these conditions, the
theory of ideal flows with vorticity [18] is the appropriate framework to
describe the structure of the falling jet.

Figure 1: Sample numerical jet with parameter definitions and cylindrical
system of coordinates for $Re=400,Fr=We=10$. Accompanying axial velocity
profiles are represented for selected axial locations.

_Theory of inviscid vortical jets_.– Consider a liquid jet emerging from a
circular tube of radius $r_{o}$, and injected downwards into a passive gaseous
atmosphere at constant pressure $p_{a}$ under the presence of a gravitational
acceleration $g$. An example of the resulting liquid-air interface, of
constant surface tension coefficient $\sigma$, can be observed in Fig. 1,
where the main variables used in the analysis are also indicated. We describe
the jet with use made of a cylindrical coordinate system $(z,r,\theta)$ with
origin at the center of the injector exit and the axial coordinate $z$
pointing in the direction of gravity. The liquid velocity field,
$\bm{v}=u(r,z)\bm{e}_{z}+v(r,z)\bm{e}_{r}$, is assumed incompressible,
inviscid and axisymmetric, and thus satisfies the momentum equation

$\bm{\nabla}\left(\dfrac{u^{2}+v^{2}}{2}+\dfrac{p}{\rho}-gz\right)=\bm{v}\wedge\left(\bm{\nabla}\wedge\bm{v}\right),$
(2)

where $p$ is the liquid pressure and $|\bm{v}|^{2}/2=(u^{2}+v^{2})/2$ is the
kinetic energy per unit mass. The projection of Eq. (2) along the jet
streamlines yields Bernoulli’s equation

$\dfrac{u^{2}+v^{2}}{2}+\dfrac{p}{\rho}-gz=C_{l},$ (3)

indicating that the total pressure remains constant along each streamline, but
is different for different streamlines. Indeed, in contrast with the case of
an irrotational flow, here the value of $C_{l}$ depends on the streamline due
to the fact that $\bm{\nabla}\wedge\bm{v}\neq 0$. The continuity equation
$\partial_{z}(ru)+\partial_{r}(rv)=0$ provides the estimate $v_{c}\sim
u_{o}r_{o}/l_{g}=u_{o}/Fr\ll u_{o}$ for the characteristic radial velocity,
whereby the kinetic energy $|\bm{v}|^{2}/2\simeq u^{2}/2$ with small
$O(Fr^{-2})$ relative errors. In addition, it is useful to introduce a stream
function $\psi$ defined by $ru=\partial_{r}\psi$ and $rv=-\partial_{z}\psi$,
satisfying the continuity equation, such that Eq. (3) reads

$U^{2}(\Psi,Z)=U^{2}(\Psi,0)+Z+P(\Psi,0)-P(\Psi,Z),$ (4)

in its simplified dimensionless form, where $U=u/u_{o}=R^{-1}{\rm d}\Psi/{\rm
d}R$, $\Psi=\psi/(u_{o}r_{o}^{2})$, $P=2p/(\rho u_{o}^{2})$, $R=r/r_{o}$ and
$Z=z/l_{g}$. To close the model, an appropriate expression for the pressure
$P$ is needed. To that end, use is made of the normal stress balance at the
free surface, $p-p_{a}=\sigma\mathcal{C}$, where $\mathcal{C}$ is twice the
mean curvature of the free surface. The slenderness assumption $Fr\gg 1$
allows to simplify the curvature, $\mathcal{C}=r_{j}^{-1}$, and therefore the
normal stress balance may be reduced to $p=p_{a}+\sigma/r_{j}$, with relative
errors $O(Fr^{-2})\ll 1$. In addition, the radial pressure variations inside
the liquid, $\Delta_{r}p$, can also be neglected with the same relative error.
Indeed, the radial momentum equation provides the estimate
$\Delta_{r}p\sim\rho u_{o}v_{c}r_{o}/l_{g}\sim\rho u_{o}^{2}\,Fr^{-2}\ll\rho
u_{o}^{2}$. With these simplifications, the non-dimensional pressure
difference between the jet exit, $Z=0$, and a generic downstream station $Z$,
reads $P(0)-P(Z)=We^{-1}(1-1/R_{j})$, where $We=\rho u_{o}^{2}r_{o}/(2\sigma)$
is the Weber number, which is the only parameter governing the structure of
the vortical ideal jet.

Assuming that the outlet velocity profile is fully developed 111The theory can
be readily applied to arbitrary outlet velocity profiles but must, in general,
be integrated numerically., $U(R,0)=2(1-R^{2})$, the stream function at $Z=0$
is obtained by radial integration as $\Psi=R^{2}(1-R^{2}/2)$, which can be
inverted to yield $R^{2}=1-(1-2\Psi)^{1/2}$. It follows then that
$U^{2}(\Psi,0)=4(1-2\Psi)$, whereby Eq. (4) reduces to the ordinary
differential equation 222In the non-slender case with $Fr\lesssim 1$ the
radial and axial velocities are comparable, and the capillary pressure must
include the contribution of the axial curvature. The mathematical model can
then be formulated as an elliptic free-boundary problem whose analysis is much
more involved than the slender approximation employed herein.

$\dfrac{1}{R^{2}}\left(\frac{\mathrm{d}\Psi}{\mathrm{d}R}\right)^{2}=4(1-2\Psi)+F(Z),$
(5)

for the function $\Psi(R;Z,R_{j},We)$, where
$F(Z)=Z-\frac{1}{We}\,\frac{1-R_{j}}{R_{j}}$, and $Z,\,R_{j}(Z)$ play the role
of parameters in the slender approximation employed here. The integration of
Eq. (5) by separation of variables yields

$\int_{0}^{1/2}\frac{{\rm d}\Psi}{\sqrt{4(1-2\Psi)+F(Z)}}=\frac{R_{j}^{2}}{2}$
(6)

where application has been made of the symmetry boundary condition at the
axis, $\Psi(R=0)=0$, and of the conservation of the flow rate,
$\Psi(R=R_{j})=1/2$. A straightforward integration of Eq. (6) finally provides
the octic

$R_{j}^{8}-(2+Z+We^{-1})R_{j}^{4}+We^{-1}R_{j}^{3}+1=0,$ (7)

which, for $We\gg 1$, has the explicit solution

$R_{j}(Z)=\sqrt{\sqrt{1+Z/4}-\sqrt{Z/4}},$ (8)

to be compared with Mariotte’s law (1). The far-field shape of the jet is
found by expanding Eq. (8), to yield $R_{j}=Z^{-1/4}-Z^{-5/4}/2+O(Z^{-9/4})$
for $Z\gg 1$, while Eq. (1) has the large-$Z$ expansion
$Z^{-1/4}-Z^{-5/4}/4+O(Z^{-9/4})$. Thus, both solutions converge to the same
asymptotic free-fall law $R_{j}\to Z^{-1/4}$ as $Z\to\infty$. In contrast, the
shapes predicted near the injector are completely different. Indeed, Eq. (8)
provides $R_{j}=1-Z^{1/2}/4+Z/32+O(Z^{3/2})$ for $Z\ll 1$, while Mariotte’s
law (1) yields $1-Z/4+5Z^{2}/32+O(Z^{3})$. It is thus deduced that the initial
contraction of the jet is much more abrupt for the jet with initial Poiseuille
profile than for the uniform jet, as a consequence of the smaller velocity of
the fluid particles near the free surface.

Figure 2: (a) and (b) Comparison of the jet shapes given by Mariotte’s law
(1), and the new solution (8), with experimental data for $G\gg 1$ and
$We\gtrsim 1$ extracted from Duda and Vrentas [12] (DV67), Philippe and
Dumargue [16] (PD91) and González-Mendizábal  _et al._ [21] (GM87). The free
surface from a steady numerical simulation of a jet for $Re=400$, $Fr=We=10$
is also presented for comparison. (c) Downstream evolution of the jet velocity
at the axis (solid lines) and at the free surface (dashed lines) obtained from
numerical simulations at different values of $Re$, together with the solution
of the ideal vortical flow (4) and Mariotte’s free-fall law (1). (d) Influence
of capillarity on the shape of the jet for several $We$ for $Re=400$, $Fr=10$,
showing that surface-tension effects are negligible for $We\gtrsim 3$.

_Effect of surface tension_.– Weak surface-tension effects can be rationalized
by means of a perturbation expansion of Eq. (7) in powers of $We^{-1}$, to
yield
$R_{j}(Z;We)=[(1+Z/4)^{1/2}-(Z/4)^{1/2}]^{1/2}\left[1+We^{-1}f(Z)\right]$. The
function $f(Z)$ has a maximum value of 0.032, whereby the relative errors
introduced by Eq. (8) compared with Eq. (7) are smaller than 1% for $We\gtrsim
3$. Similarly, surface tension effects can also be introduced in Mariotte’s
solution (1) by considering the capillary pressure jump in the slender
approximation, to yield $U_{j}^{2}=1+Z-We^{-1}(1-R_{j})/R_{j}$, leading to the
quartic $(1+Z+We^{-1})R_{j}^{4}-We^{-1}R_{j}^{3}-1=0$ which has a large-$We$
expansion $R_{j}=(1+Z)^{-1/4}[1+We^{-1}s(Z)]$, where the function $s(Z)$ has a
single maximum $s(Z=175/81)=27/1024\approx 0.0264$, so that the relative error
introduced by Eq. (1) is smaller than 1% for $We\gtrsim 2$. These conclusions
are confirmed by the results of Fig. 2(d), where full Navier-Stokes
simulations are compared with Eq. (8) for several values of $We$.

_Numerical simulations_.– The steady axisymmetric Navier-Stokes equations were
integrated numerically using Comsol. Finite elements are employed to
discretize the deformable domain, whose displacement field is computed by
solving a Laplace equation within an arbitrary Lagrangian-Eulerian (ALE)
framework. Conversely to what has been presented above, for this section we
will take $r_{o}$ as the characteristic length scale for both axial and radial
directions so that the dimensionless equations of motion read

$\bm{\nabla}\cdot\bm{v}=0,\quad\bm{v}\cdot\bm{\nabla}\bm{v}=-\bm{\nabla}P/2+\bm{\nabla}\cdot\bar{\bar{\tau}}+(2Fr)^{-1}\,\bm{e}_{z},$
(9)

where $\bar{\bar{\tau}}=(\bm{\nabla}\bm{v}+\bm{\nabla}\bm{v}^{T})/Re$ is the
viscous stress tensor, and the necessary boundary conditions are
$U(-L/r_{o},R)-2(1-R^{2})=V(-L/r_{o},R)=0$ (Poiseuille flow prescribed at the
upstream boundary of an injector of length $L$), $U(Z,1)=V(Z,1)=0$ (no-slip at
the injector inner wall) for $Z\leq 0$, $V(Z,0)=0$ (radial symmetry), and
$-(P/2)\bm{n}+\bar{\bar{\tau}}\cdot\bm{n}=-(2We)^{-1}\bm{n}\bm{\nabla}\cdot\bm{n}$
(stress balance) at the free surface $R=R_{j}(Z)$, which needs to be
determined as part of the solution, where $\bm{n}$ is the outwards normal to
the interface. To that end, the vertical displacement of the mesh at the free-
surface boundary is left as a Lagrange multiplier that enforces the kinematic
condition $\bm{v}\cdot\bm{n}=0$ or, equivalently,
$U\,\partial{R_{j}}/\partial{Z}=V$. A stress-free boundary condition at $Z\gg
1$ completes the description of the numerical set-up. The solution is found
employing a root-finding algorithm with a prescribed normalized tolerance not
greater than $10^{-4}$. The prescription of Eq. (8) as initial seed for the
free-surface shape was seen to facilitate the convergence greatly. Several
other initial guesses were also employed to assess numerical independence. An
example of the numerical jet so obtained is presented in Fig. 1 for $Re=400$,
$Fr=We=10$.

_Comparison with experiments_.– Due to the fundamental and applied interest of
free liquid jets, a large number of experimental studies can be found in the
literature covering many different aspects of the jet structure, instability
and break-up [1, 2]. Among these studies, several ones were devoted to the
generation of slender gravity-dominated jets emerging from long injectors
fulfilling the hypotheses behind our analysis. Figure 2(a,b) shows the free-
surface shape extracted from several available experiments (symbols), a
Navier-Stokes simulation at $Re=400$ and $Fr=We=10$, together with the
classical free-fall solution given by Eq. (1) and our parameter-free, leading-
order solution (8). In all cases the experimental jet shapes are much closer
to our free-fall law (8) than to Mariotte’s equation (1), most notably in the
region close to the injector outlet, as illustrated in Fig. 2(b). The slight
deviations from Eq. (8) observed in Fig. 2 can be attributed to the finite
values of $G$ in the experiments, leading to a non-negligible radial diffusion
of momentum across the jet core (see e.g. the star symbols extracted from Ref.
[16]). The downstream evolution of the axial velocity at the jet centerline
and at the free surface provided by the numerical simulation is represented in
Fig. 2(c) for several values of $Re$, together with the predictions of (4) for
the ideal vortical flow, as well as Torricelli’s free-fall law. These results
clearly show that the Navier-Stokes flow converges to the ideal flow solution
as the Reynolds number is increased.

_Concluding remarks_.– Our new analytical solution (8) for the steady
structure of slender laminar liquid jets with fully-developed outlet velocity
profiles provides much better agreement with experiments than the routinely
used Mariotte’s law (1). This comparison unambiguously demonstrates that,
while Mariotte’s solution fails to describe liquid jets with vorticity at
large Reynolds numbers, the new theory predicts a universal shape function
that collapses all the experiments onto a single master curve.

The hypotheses behind the theory, namely $Fr\gg 1$ and $Re\gg Fr$, clearly
need some justification. To that end, it proves convenient to introduce the
Morton number $Mo=g\mu^{4}/\rho\sigma^{3}$, and the Bond number $Bo=\rho
gr_{o}^{2}/\sigma$ as auxiliary dimensionless parameters, in terms of which
$Fr=We/Bo$ and $G=Re/Fr=\sqrt{2}Mo^{-1/4}Bo^{5/4}We^{-1/2}$. Large values of
the Froude number require that $Bo\ll We$, a condition that is accomplished
for injectors of sufficiently small diameter in the jetting regime, which
requires $We\gtrsim 1$ for $Bo\ll 1$ [22, 23]. In contrast, the condition
$G\gg 1$ depends crucially on the viscosity of the working liquid. For
instance, the value of $Mo\approx 2.5\times 10^{-11}$ in the case of water at
room temperature, in which case $G\approx 631Bo^{5/4}We^{-1/2}$, and the
validity condition writes $We\ll 4\times 10^{5}Bo^{5/2}$, which is
accomplished under the Rayleigh break-up regime for any realistic injector
diameter [24]. Although similar conclusions hold for any low-viscosity liquid,
the validity of our theory breaks down for liquids of sufficiently large
viscosity, for which either full Navier-Stokes simulations or simplified one-
dimensional descriptions may be used instead [1, 2].

###### Acknowledgements.

Useful discussions with A. Martínez-Calvo and M. Rubio-Rubio are gratefully
acknowledged. The authors thank the Spanish MCIU-Agencia Estatal de
Investigación through project PID2020-115655GB-C22, partly financed through
FEDER European funds.

## References

  * Eggers [1997] J. Eggers, Nonlinear dynamics and breakup of free surface flows, Rev. Mod. Phys. 69, 865 (1997).
  * Eggers and Villermaux [2008] J. Eggers and E. Villermaux, Physics of liquid jets, Rep. Prog. Phys. 71, 036601 (2008).
  * Da Vinci [1980] L. Da Vinci,  _Codex Leicester_ (London: Christie’s, 1980).
  * Mariotte [1686] E. Mariotte,  _Traite du Mouvement des Eaux et Des Autres Corps Fluides_ (Paris: E. Michallet, 1686).
  * Torricelli [1644] E. Torricelli,  _Opera geometrica_ (Florentiae, typis A. Masse & L. de Landis, 1644).
  * Senchenko and Bohr [2005] S. Senchenko and T. Bohr, Shape and stability of a viscous thread, Phys. Rev. E 71, 056301 (2005).
  * Geer and Strikwerda [1983] J. F. Geer and J. C. Strikwerda, Vertical slender jets with surface tension, J. Fluid Mech. 135, 155 (1983).
  * Eggers and Dupont [1994] J. Eggers and T. F. Dupont, Drop formation in a one-dimensional approximation of the Navier-Stokes equation, J. Fluid Mech. 262, 205 (1994).
  * García and Castellanos [1994] F. J. García and A. Castellanos, One-dimensional models for slender axisymmetric viscous liquid jets, Phys. Fluids 6, 2676 (1994).
  * Weber [1931] C. Weber, On the breakdown of a fluid jet, Z. Angew. Math. Mech. 11, 136 (1931).
  * Sterling and Sleicher [1975] A. M. Sterling and C. A. Sleicher, The instability of capillary jets, J. Fluid Mech. 68, 477 (1975).
  * Duda and Vrentas [1967] J. L. Duda and J. S. Vrentas, Fluid mechanics of laminar liquid jets, Chem. Eng. Sci. 22, 855 (1967).
  * Goren [1966] S. L. Goren, Development of the boundary layer at a free surface from a uniform shear flow, J. Fluid Mech. 25, 87 (1966).
  * Goren and Wronski [1966] S. L. Goren and S. Wronski, The shape of low-speed capillary jets of Newtonian liquids, J. Fluid Mech. 25, 185 (1966).
  * Og̃uz [1998] H. Og̃uz, On the relaxation of laminar jets at high Reynolds numbers, Phys. Fluids 10, 361 (1998).
  * Philippe and Dumargue [1991] C. Philippe and P. Dumargue, Etude de l’etablissement d’un jet liquide laminaire emergeant d’une conduite cylindrique verticale semi-infinie et soumis a l’influence de la gravite, ZAMP 42, 227 (1991).
  * Sevilla [2011] A. Sevilla, The effect of viscous relaxation on the spatiotemporal stability of capillary jets, J. Fluid Mech. 684, 204 (2011).
  * Batchelor [2000] G. K. Batchelor,  _An introduction to fluid dynamics_ (Cambridge University Press, 2000).
  * Note [1] The theory can be readily applied to arbitrary outlet velocity profiles but must, in general, be integrated numerically.
  * Note [2] In the non-slender case with $Fr\lesssim 1$ the radial and axial velocities are comparable, and the capillary pressure must include the contribution of the axial curvature. The mathematical model can then be formulated as an elliptic free-boundary problem whose analysis is much more involved than the slender approximation employed herein.
  * González-Mendizábal  _et al._ [1987] D. González-Mendizábal, C. Olivera-Fuentes, and J. M. Guzmán, Hydrodynamics of laminar liquid jets. Experimental study and comparison with two models, Chem. Eng. Comm. 56, 117 (1987).
  * Clanet and Lasheras [1999] C. Clanet and J. C. Lasheras, Transition from dripping to jetting, J. Fluid Mech. 383, 307 (1999).
  * Rubio-Rubio  _et al._ [2013] M. Rubio-Rubio, A. Sevilla, and J. M. Gordillo, On the thinnest steady threads obtained by gravitational stretching of capillary jets, J. Fluid Mech. 729, 471 (2013).
  * Gordillo  _et al._ [2001] J. M. Gordillo, M. Pérez-Saborid, and A. M. Ganán-Calvo, Linear stability of co-flowing liquid–gas jets, J. Fluid Mech. 448, 23 (2001).