# Nonlinear wave interactions
in geochemical modeling.

A. C. Alvarez, J. Bruining, D. Marchesin

###### Abstract

This paper is concerned with the study of the main wave interactions in a
system of conservation laws in geochemical modeling. We study the modeling of
the chemical complexes on the rock surface. The presence of stable surface
complexes affects the relative permeability. We add terms representing surface
complexes to the accumulation function in the model presented in [14]. This
addition allows to take into account the interaction of ions with the rock
surface in the modeling of the oil recovery by the injection of carbonated
water. Compatibility hypotheses with the modeling are made on the coefficients
of the system to obtain meaningful solutions. We developed a Riemann solver
taking into account the complexity of the interactions and bifurcations of
nonlinear waves. Such bifurcations occur at the inflection and resonance
surfaces. We present the solution of a generalized eigenvalue problem in a
(n+1)-dimensional case, which allows the construction of rarefaction curves. A
method to find the discontinuous solutions is also presented. We find the
solution path for some examples.

Keywords— Carbonated water, bifurcations, geochemical surface complexes,
Riemann solution, interaction with minerals, conservation laws.

##  1 Introduction

Composition modeling deals with the transport of ions and chemical species in
porous media. Modeling of such processes is based on a system of non-linear
partial differential equations composed of Darcy’s law and mass conservation
equations. Several processes can be studied with this system aimed at
understanding and improving oil recovery (see [11]). In this sense the study
of the Riemann problem for such flows plays a major role in the field of
Mathematics as well as in Engineering.

The class of models described in this paper has first appeared in problems of
geochemical flow in porous media, see for example [5, 3]. In the present work
we focus on the mathematical features of the above mentioned physical model.
This system represents conservation laws for a set of $\displaystyle n+1$
independent chemical species. These chemical species depend on a set of
$\displaystyle m$ ion compound concentrations, which are assumed to have
reached thermodynamic equilibrium. This situation leads to a system of
$\displaystyle n+1$ nonlinear partial equations and $\displaystyle k$
algebraic constraint equations. In this paper, we study the Riemann problem
for this kind of system.

We focus on mathematical features of systems of conservation laws in the
manifold given by the constraints. To simplify the above problem and to obtain
an equivalent system with only $\displaystyle n+1$ conservation laws, a novel
strategy based on Gibbs rule that reduces the dimension of the system is
developed in [5, 14]. This procedure was used successfully in studying the
injection of low salinity water.

This situation reduces to a system of conservation laws of the type

$\displaystyle{\frac{\partial\textbf{G}(U)}{\partial t}+\frac{\partial
u\textbf{F}(U)}{\partial x}=0,}$ (1)

with state variables $\displaystyle U=(U_{1},\ldots,U_{n})$ belonging to phase
space $\displaystyle\Omega$, with $\displaystyle
U_{i}(x,t):\mathbb{R}\times[0,+\infty)\longrightarrow\mathbb{R}$, accumulation
functions
$\displaystyle\textbf{G}(U)=(G_{1}(U),\cdots,G_{n+1}(U))^{T}:\Omega\longrightarrow\mathbb{R}^{n+1}$,
flux functions
$\displaystyle\textbf{F}(U)=(F_{1}(U),\cdots,F_{n+1}(U))^{T}:\Omega\longrightarrow\mathbb{R}^{n+1}$
and Darcy velocity $\displaystyle
u=u(x,t):\mathbb{R}\times[0,+\infty)\longrightarrow[0,+\infty)$. The physical
and chemical details of this model are described in [1], while in this work we
focus on mathematical and numerical aspects.

System (1) is typical of two phase multicomponent models for geochemical one-
dimensional flow, see, e.g., [21]. In this paper, we focus on the system with
accumulation and flux functions given by

$G_{j}=\varphi\rho_{wj}s_{w}+\varphi\rho_{oj}s_{o}+(1-\varphi)\rho_{rj}\quad\text{
and }\quad F_{j}=\rho_{wj}f_{w}+\rho_{oj}f_{o},\quad j=1,\cdots,n+1$ (2)

where $\displaystyle\varphi$ is a given parameter, $\displaystyle
s_{w}+s_{o}=1$ and $\displaystyle f_{w}+f_{o}=1$. In order to simplify the
notation we use $\displaystyle s=s_{w}$, $\displaystyle f=f_{w}$,
$\displaystyle s_{o}=1-s$ and $\displaystyle f_{o}=1-f$.

In such a model the state variables are the composition variables
$\displaystyle y_{i}$, for $\displaystyle i=1,\cdots,n-1$, the saturation
variable $\displaystyle s_{w}$ and the Darcy velocity $\displaystyle u$. The
coefficients $\displaystyle\rho_{wj}$, $\displaystyle\rho_{oj}$ and
$\displaystyle\rho_{rj}$ depend on the compositions $\displaystyle
y=(y_{1},\cdots,y_{n-1})$. Here $\displaystyle f_{\alpha}$, which is called
the fractional flow of phase $\displaystyle\alpha$, for
$\displaystyle\alpha=w,o$, depends on $\displaystyle s$ and $\displaystyle
y=(y_{1},\cdots,y_{n-1})$ and are written as $\displaystyle
f_{\alpha}=\lambda_{\alpha}/(\lambda_{w}+\lambda_{o})$, for
$\displaystyle\alpha=w$ and $\displaystyle o$, in which
$\displaystyle\lambda_{\alpha}=k_{r\alpha}/\mu_{\alpha}$. The functions
$\displaystyle\lambda_{\alpha}$, $\displaystyle k_{r\alpha}$ and
$\displaystyle\mu_{\alpha}$ are called mobility, relative permeability and
viscosity of phase $\displaystyle\alpha$ which depend on the variables
$\displaystyle s$ and $\displaystyle y$. The subscript $\displaystyle w$ in
$\displaystyle\rho_{wj}$ indicates that these functions refer to the water
phase while the subscript $\displaystyle o$ for $\displaystyle\rho_{oj}$
indicates the oil phase. Finally, $\displaystyle\rho_{rj}$ represent the ion
interactions with the porous medium.

Here, we incorporate to the model presented in [14], the coefficients
$\displaystyle\rho_{rj}$ in the accumulation functions $\displaystyle G$,
which arise from the interaction of ions with the rock surfaces. Such chemical
interaction enables incorporating the formation of surface complexes, which
enhances the oil recovery by injection of low salinity carbonated water.
Therefore, when these coefficients are introduced new features appear in the
Riemann solutions associated to the system (1).

In our study we consider the physical state variables
$\displaystyle\Omega=\\{(s,y)=[0,1]\times\mathcal{K}\subset\mathbb{R}\times\mathbb{R}^{n-1}\\}$,
where $\displaystyle\mathcal{K}$ denotes the $\displaystyle n-1$-dimensional
hypercube. The physical model consists of $\displaystyle n-1$ species
concentrations and two additional variables $\displaystyle(s,u)$, leading to
$\displaystyle n+1$ unknowns $\displaystyle U=(s,y,u)$.

A Riemann problem consists of a Cauchy (initial value) problem governed by
equations of type (1) with initial data

$U(x,t=0)=\left\\{\begin{matrix}U_{L}&\mbox{if }x<0,\\\ U_{R}&\mbox{if
}x>0.\end{matrix}\right.$ (3)

Riemann solutions are obtained using rarefaction, shock and composite wave
curves going from the state $\displaystyle U_{L}$ to $\displaystyle U_{R}$
passing through a set of intermediate constant states. Rarefactions are
continuous self-similar solutions of (1), which are represented by

$U=\widehat{U}(\xi),\quad\text{with}\quad\xi=x/t.$ (4)

Substituting (4) into system (1) we obtain the generalized eigenvalue problem

$Ar=\lambda Br,\quad\text{where}\quad A={\partial F}/{\partial U},\quad
B={\partial G}/{\partial U}.$ (5)

The eigenvector $\displaystyle r$ is parallel to $\displaystyle
d\widehat{U}/d\xi$, so the rarefaction curves are tangent to the
characteristic field given by the normalized eigenvector $\displaystyle r$.

A shock wave is a traveling discontinuity in a (weak) solution of system (1)
given by

$U(x,t)=\left\\{\begin{matrix}U^{-}&\mbox{if }x<\sigma t,\\\ U^{+}&\mbox{if
}x>\sigma t,\end{matrix}\right.$ (6)

where $\displaystyle\sigma$ is a real constant called the shock speed.
Solution (6) is a piecewise constant weak solution for the Riemann problem
defined by (1) and (3) provided $\displaystyle U^{+}$ and $\displaystyle
U^{-}$ satisfy the Rankine-Hugoniot condition:

$F(U^{-})-F(U^{+})=\sigma(G(U^{-})-G(U^{+})).$ (7)

In this work we use the wave curve method (see [17, 18, 4]), a valuable tool
to obtain semianalytical solutions. More details about the wave curve method
can be found in [15, 4, 12, 9], with theoretical justification provided by
[17, 18, 10, 23, 2] and references therein. Here, we provide a method for
calculating and concatenating waves in a set of relevant examples. To do this,
we provide the waves and the bifurcation surfaces necessary to solve the
Riemann problem for the system of conservation laws (1). To study this system
in the (n+1)-dimensional case is a hard task in general, however due to the
structure of the flux and accumulation functions we can obtain some useful
information about rarefaction and shock waves.

In general the system (1) is not strictly hyperbolic, i.e., the generalized
eigenvalues in (5) must coincide on some hypersurface. Moreover, genuine
nonlinearity fails, i.e., there is a locus where the characteristic speed has
zero directional derivative along the vector field $\displaystyle r=r_{i}(U)$,
i.e. $\displaystyle\lambda^{\prime}_{i}(r)=\nabla\lambda_{i}\cdot r_{i}$ = 0.
The rarefaction curve generically stops at this locus. Thus the determination
of the inflection locus is relevant for the wave curve method.

Following Liu (see [17, 18]) to cross the inflection locus, the rarefactions
wave need to be attached to a shock. This corresponds to the construction of a
composite curve, which arises in state space during the construction of wave
curves for the solution of a Riemann problem for non-strictly hyperbolic
systems of conservation laws.

There is one wave in which the saturation changes (called the Buckley-Leverett
wave (B-L)). There are also $\displaystyle n-1$ chemical waves. Having
identified these waves we provide the analysis of the Riemann solutions in
certain relevant situations.

An in-house Riemann solver was developed to represent wave curves that satisfy
compatibility and admissibility criteria (see [2] for technical details of
implementation). In this solver, we automate the construction of the solution
paths taking into account the bifurcation surfaces.

The paper is organized as follows. Section 2 presents the fractional flow
function used here. In Section 3, we obtain the characteristic speeds and
characteristic vectors of the system. In Section 4 we obtain the bifurcation
surfaces crucial in determining the wave interactions. In Section 5, we obtain
the Hugoniot locus. This locus is the main ingredient to obtain the Riemann
solutions. In Section 7, we solve the Riemann problem for some examples. In
two numerical examples, we study the case when the coefficients
$\displaystyle\rho_{jr}$ are constant. In Section 8, we draw our conclusions.

##  2 Fractional flow functions

The fractional flows $\displaystyle f_{w}$ and $\displaystyle f_{o}$ for water
and oil in (2) are frequently taken as saturation-dependent functions defined
as follows. We denote $\displaystyle
s_{e}=(s-s_{wc})/(1-s_{wc}),\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ s\geq s_{wc}\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ s_{e}=0\leavevmode\nobreak\ \text{for}$
$\displaystyle s<s_{wc};k_{rw}(s)=s_{e}^{(2/\lambda+3)}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{ and }\leavevmode\nobreak\ \leavevmode\nobreak\
k_{ro}(s_{o})=(1-s_{e})^{2}(1-s_{e}^{(2/\lambda+1)}),$ ($\displaystyle
s_{o}=1-s$) ([6]). Here, the water viscosity is taken as
$\displaystyle\mu_{w}=0.001$ and the oil viscosity as
$\displaystyle\mu_{o}=0.002$ when they are constant. In other examples
$\displaystyle\mu_{w}$ and $\displaystyle\mu_{o}$ depend on $\displaystyle y$
also. The fractional flow functions for the water and oil phases are given by

$f_{w}(s)=\frac{k_{rw}(s)/\mu_{w}}{(k_{rw}(s)/\mu_{w}+k_{ro}(s_{o})/\mu_{o})},\quad\text{and}\quad
f_{o}(s)=1-f_{w}(s).$ (8)

The fractional flow $\displaystyle f_{w}$ depends on the chemical variables
via the oil viscosity. In order to take into account this dependence in the
analysis of the wave interaction, we assume that the fractional flow
$\displaystyle f=f_{w}$ is $\displaystyle S$-shaped (see [7]):

######  Assumption 1.

S-shape. We assume that the fractional flows $\displaystyle f(s,y)$ has the
following suitable properties: $\displaystyle f(0,y)=0$ and $\displaystyle
f(1,y)=1$, $\displaystyle\frac{\partial f(0,y)}{\partial s}=\frac{\partial
f(1,y)}{\partial s}=0$. For each $\displaystyle y$, $\displaystyle f$ has one
inflection point $\displaystyle s^{*}$ (that depends on $\displaystyle y$ and
possibly on viscosities)

$s^{*}=s^{*}(y)\quad\text{such
that}\quad\frac{\partial^{2}f(s^{*},y)}{\partial s^{2}}=0.$ (9)

For each $\displaystyle y$, there are two states

$s^{{\dagger},1}(y)<s^{*}(y)<s^{{\dagger},2}(y)\quad\text{such
that}\quad\frac{\partial f(s^{{\dagger},1},\cdot)}{\partial s}=\frac{\partial
f(s^{{\dagger},2},\cdot)}{\partial s}=1.$ (10)

For $\displaystyle 0<s<1$, $\displaystyle\displaystyle{\frac{\partial
f(s,\cdot)}{\partial s}>0}$. If $\displaystyle 0\leq s<s^{*}(y)$, then
$\displaystyle\displaystyle{\frac{\partial^{2}f(s,y)}{\partial s^{2}}>0}$, and
if $\displaystyle s^{*}(y)<s\leq 1$, then
$\displaystyle\displaystyle{\frac{\partial^{2}f(s,y)}{\partial s^{2}}<0}$.

##  3 Eigenvalues analysis and elementary waves

The system has $\displaystyle n+1$ equations with the unknown variables
$\displaystyle(s,y,u)$ with $\displaystyle y=(y_{1},y_{2},\cdots,y_{n-1})$
given by

$\displaystyle\displaystyle\frac{\partial}{\partial
t}(\varphi\rho_{wj}(y)s_{w}+\varphi\rho_{oj}(y)s_{o}+(1-\varphi)\rho_{rj}(y))+\frac{\partial}{\partial
x}(u(\rho_{wj}(y)f_{w}+\rho_{oj}(y)f_{o}))=0.$ (11)

We take $\displaystyle s=s_{w}$, $\displaystyle f=f_{w}$, $\displaystyle
s_{o}=1-s$ and $\displaystyle f_{o}=1-f_{w}$. The accumulation and flux
functions $\displaystyle G$ and $\displaystyle uF$ are given in terms of

$\displaystyle\displaystyle
G_{j}=\varphi\rho_{wj}(y)s_{w}+\varphi\rho_{oj}(y)s_{o}+(1-\varphi)\rho_{rj}(y)$
(12) $\displaystyle\displaystyle F_{j}=\rho_{wj}(y)f_{w}+\rho_{oj}(y)f_{o}$
(13)

The index $\displaystyle w$ (water) is often replaced by the index
$\displaystyle a$ (aqueous phase), and the index $\displaystyle
j=1,\cdots,n+1$ is used to denote chemical species.

The main object to obtain the Riemann solutions of hyperbolic systems are the
rarefaction and shock curves. The rarefactions are obtained from integral
curves of the line fields given by the eigenvectors of system $\displaystyle
A\vec{r}=B\lambda\vec{r}$, with the Jacobian matrices $\displaystyle B=DG$ and
$\displaystyle A=D(uF)$ of the accumulation and flux functions $\displaystyle
uF$ and $\displaystyle G$. For system $\displaystyle(\ref{eq1})$, the matrices
$\displaystyle B$ and $\displaystyle A$ are given by

$\displaystyle\displaystyle B_{i,1}=\varphi[\rho_{i}],\quad
B_{i,k+1}=\varphi\frac{\partial\rho_{wi}}{\partial
y_{k}}s+\varphi\frac{\partial\rho_{oi}}{\partial
y_{k}}s_{o}+(1-\varphi)\frac{\partial\rho_{ri}}{\partial y_{k}},\text{ and
}B_{i,n+1}=0.$ (14) $\displaystyle\displaystyle
A_{i,1}=u[\rho_{i}]\frac{\partial f}{\partial s},\quad
A_{i,k+1}=u\left(\frac{\partial\rho_{wi}}{\partial
y_{k}}f+\frac{\partial\rho_{oi}}{\partial y_{k}}f_{o}\right)\text{ and
}A_{i,n+1}=F_{i},$ (15)

in which $\displaystyle i=1,\cdots,n+1$, $\displaystyle k=1,n-1$ and

$[\rho_{i}]=\rho_{wi}-\rho_{oi}.$ (16)

From the Jacobian matrices, we obtain the eigenpairs, which we summarize in

######  Proposition 1.

The eigenpairs of the eigenvalue problem $\displaystyle A\vec{r}=\lambda
B\vec{r}$, where the matrices $\displaystyle B$ and $\displaystyle A$
represent the Jacobian of the accumulation and flux terms of system
$\displaystyle(\ref{eq1})$ are the Buckley-Leverett eigenpair
($\displaystyle\lambda_{s},\vec{r}_{s}$) given by

$\lambda_{s}=\frac{u}{\varphi}\frac{\partial f}{\partial s}\quad\text{ and
}\quad\vec{r}_{s}=(1,\cdots,0)^{T},$ (17)

together with the $\displaystyle n-1$ chemical eigenpairs
$\displaystyle(\lambda_{\Lambda_{i}},\vec{r}_{\Lambda_{i}})$ given by

$\lambda_{\Lambda_{i}}=\frac{u}{\varphi}\frac{f-\Lambda_{i}}{s-\Lambda_{i}},$
(18)

with $\displaystyle i=1,\cdots,n-1$. We obtain $\displaystyle\Lambda_{i}$ and
$\displaystyle\vec{v}_{i}$ as the solutions of the generalized eigenvalue
problem

$((\mathcal{A}+\mathcal{A}_{1})-\Lambda_{i}(\mathcal{B}+\mathcal{B}_{1}))\vec{v}_{i}=0,$
(19)

where the matrices $\displaystyle\mathcal{A},\mathcal{A}_{1},\mathcal{B}$ and
$\displaystyle\mathcal{B}_{1}$ are given in (36)-(39). The matrices
$\displaystyle\mathcal{A}$ and $\displaystyle\mathcal{B}$ depend on the
variable $\displaystyle y$ while $\displaystyle\mathcal{A}_{1}$ and
$\displaystyle\mathcal{B}_{1}$ depend on the variables $\displaystyle s,y$.
Moreover, we obtain
$\displaystyle\vec{r}_{\Lambda_{i}}=(r_{\Lambda_{i}}^{1},\cdots,r_{\Lambda_{i}}^{n+1})^{T}$
as

$\vec{r}_{\Lambda_{i}}=P\vec{v}_{i},\\\ $ (20)

where the matrix $\displaystyle P$ is defined in (47). Furthermore, the left
eigenvectors $\displaystyle l_{\Lambda_{i}}$ are given by
$\displaystyle(0,l_{\vec{v_{i}}},0)$ where $\displaystyle l_{\vec{v_{i}}}$ are
the left eigenvectors associated to (19), i.e.,

$l_{\vec{v_{i}}}((\mathcal{A}+\mathcal{A}_{1})-\Lambda_{i}(\mathcal{B}+\mathcal{B}_{1}))=0.$
(21)

Proof of Proposition

The idea of the proof consists of reducing the matrix
$\displaystyle\mathcal{G}$ by straightforward application of the Gauss
elimination procedure, i.e., transforming the matrix
$\displaystyle\mathcal{G}$ into an equivalent matrix and obtaining the
eigenvalues and eigenvectors easily.

To determine the eigenvalues we solve $\displaystyle det(A-\lambda B)=0$, and
to calculate eigenvectors we solve $\displaystyle(A-\lambda B)\vec{r}=0$,
where the matrices $\displaystyle A$ and $\displaystyle B$ are given in (15).
Let us denote by $\displaystyle\mathcal{G}=(A-\lambda B)$, and
$\displaystyle\mathcal{G}=(\mathcal{G}_{i,j})$ for $\displaystyle
i,j=1,\cdots,n+1$. Here we use an auxiliary index $\displaystyle k$ that
ranges from $\displaystyle 1$ to $\displaystyle n-1$ and we write
$\displaystyle\mathcal{G}_{i,j}$ as

$\mathcal{G}_{i,1}(\lambda)=[\rho_{i}]\xi_{1}(\lambda),\mathcal{G}_{i,k+1}(\lambda)=\frac{\partial\rho_{wi}}{\partial
y_{k}}\xi_{2}(\lambda)+\frac{\partial\rho_{oi}}{\partial
y_{k}}\xi_{3}(\lambda)-(1-\varphi)\frac{\partial\rho_{ri}}{\partial
y_{k}}\lambda\text{ , }\mathcal{G}_{i,n+1}=F_{i}.$ (22)

where we define the auxiliary functions

$\xi_{1}(\lambda)=\left(u\frac{\partial f}{\partial
s}-\varphi\lambda\right),\quad\xi_{2}(\lambda)=(uf-\varphi\lambda
s)\quad\text{ and }\quad\xi_{3}(\lambda)=(uf_{o}-\varphi\lambda s_{o}).$ (23)

To express calculations done in the next steps of the proof, it is useful to
define the following coefficients $\displaystyle\gamma_{ij}$,
$\displaystyle\varrho_{ij}$, $\displaystyle\nu_{i}$,
$\displaystyle\vartheta_{ij}$ and $\displaystyle\varsigma_{ij}$ for
$\displaystyle i=1,\cdots,n+1$, $\displaystyle j=1,n-1$ as

$\displaystyle\displaystyle\gamma_{ij}=\frac{\partial\rho_{wi}}{\partial
y_{j}}[\rho_{1}]-\frac{\partial\rho_{w1}}{\partial
y_{j}}[\rho_{i}],\quad\varrho_{ij}=\frac{\partial\rho_{oi}}{\partial
y_{j}}[\rho_{1}]-\frac{\partial\rho_{o1}}{\partial y_{j}}[\rho_{i}],$ (24)
$\displaystyle\displaystyle\pi_{ij}=(1-\varphi)(\frac{\partial\rho_{ri}}{\partial
y_{j}}[\rho_{1}]-\frac{\partial\rho_{r1}}{\partial
y_{j}}[\rho_{i}]),\quad\nu_{i}=[\rho_{1}]F_{i}-[\rho_{i}]F_{1},$ (25)
$\displaystyle\displaystyle\vartheta_{ij}=\gamma_{ij}\nu_{4}-\gamma_{4,j}\nu_{i}\quad\varsigma_{ij}=\varrho_{ij}\nu_{4}-\varrho_{4,j}\nu_{i},$
(26) $\displaystyle\displaystyle\tau_{ij}=\pi_{ij}\nu_{4}-\pi_{4,j}\nu_{i}.$
(27)

All these coefficients ((24)- (27)) depend only of the variables
$\displaystyle y$. Now, we are able to perform the Gaussian elimination

  * 1)

Substituting the $\displaystyle i$-th row of matrix
$\displaystyle(\ref{matrixnew})$, for $\displaystyle i=2,\cdots,n+1$, by the
sum of the first row of $\displaystyle(A-\lambda B)$ multiplied by
$\displaystyle-[\rho_{i}]$ with its $\displaystyle i$-th row multiplied by
$\displaystyle[\rho_{1}]$, we obtain, for $\displaystyle j=1,n-1$ and
$\displaystyle i=2,\cdots,n+1$

$\begin{pmatrix}\mathcal{G}_{1,1}(\lambda)&\mathcal{G}_{1,j+1}(\lambda)&\mathcal{G}_{1,n+1}\\\
\mathbb{O}&\mathcal{G}^{(1)}_{i,j+1}(\lambda)&\mathcal{G}^{(1)}_{i,n+1}\\\
\end{pmatrix}\vec{r}=0.$ (28)

Here $\displaystyle\mathbb{O}$ is the column vector of three zeros and block
matrices $\displaystyle(\mathcal{G}^{(1)}_{i,j+1})$ and
$\displaystyle(\mathcal{G}^{(1)}_{i,n+1})$ for $\displaystyle i=2,\cdots,n+1$
and $\displaystyle j=1,n-1$ are the block matrices with elements

$\mathcal{G}^{(1)}_{i,j+1}(\lambda)=\left(\gamma_{ij}\xi_{2}(\lambda)+\varrho_{ij}\xi_{3}(\lambda)-\pi_{ij}\lambda\right)_{1\leq
j\leq n-1}\quad\text{ and }\quad\mathcal{G}^{(1)}_{i,n+1}=\nu_{i},$ (29)

where $\displaystyle\gamma_{ij}$, $\displaystyle\varrho_{ij}$
$\displaystyle\pi_{ij}$ and $\displaystyle\nu_{i}$ are given by Eq. (24) and
(25).

Notice here that if $\displaystyle[\rho_{i}]$ given in
$\displaystyle(\ref{diferenca})$ is zero for some index $\displaystyle i$ at
some $\displaystyle y$, the corresponding entry in the first column is zero
and we do not need to perform Gaussian elimination. Moreover, if
$\displaystyle[\rho_{1}]$ is zero for some $\displaystyle y$, we interchange
the first row with another row to obtain a non-zero pivot.

  * 2)

Assuming that $\displaystyle\nu_{n+1}$ is non-zero, we substitute the
$\displaystyle i$-th row of $\displaystyle(\ref{matrix3})$,i .e.,
($\displaystyle\mathcal{G}^{(1)}_{i,j+1}(\lambda)$, $\displaystyle
j=1,\cdots,n$), for $\displaystyle i=2,n$, by the sum of the
$\displaystyle(n+1)$-th row of $\displaystyle(\ref{matrix3})$,
($\displaystyle\mathcal{G}^{(1)}_{n+1,j+1}(\lambda)$, $\displaystyle
j=1,\cdots,n$) multiplied by $\displaystyle-\nu_{i}$ with the $\displaystyle
i$-th row of $\displaystyle(\ref{matrix3})$ multiplied by
$\displaystyle\nu_{n+1}$, Notice that the last column from $\displaystyle
i=2,\cdots,n$ becomes $\displaystyle\nu_{n+1}(-\nu_{i})-\nu_{n+1}\nu_{i}=0$,
and the first column for $\displaystyle i=2,\cdots,n$ remains as zero. Thus,
we obtain for $\displaystyle j=1,\cdots,n-1$ and $\displaystyle
i=2,\cdots,n-1$

$\begin{pmatrix}\mathcal{G}_{1,1}(\lambda)&\mathcal{G}_{1,j+1}(\lambda)&\mathcal{G}_{1,n+1}\\\
\mathbb{O}&\mathcal{G}^{(2)}_{i,j+1}(\lambda)&\mathbb{O}\\\
0&\mathcal{G}^{(1)}_{n+1,j+1}(\lambda)&\nu_{n+1}\\\ \end{pmatrix}\vec{r}=0.$
(30)

Here $\displaystyle\mathbb{O}$ is the column vector of $\displaystyle n-1$
zeros. Also, $\displaystyle\mathcal{G}^{(2)}_{i,j+1}(\lambda)$ and
$\displaystyle\mathcal{G}^{(1)}_{n+1,j+1}(\lambda)$ for $\displaystyle
i=2,\cdots,n$ and $\displaystyle j=1,\cdots,n-1$ are the block matrices given
by:

$\mathcal{G}^{(2)}_{i,j+1}(\lambda)=\vartheta_{ij}\xi_{2}(\lambda)+\varsigma_{ij}\xi_{3}(\lambda)-\tau_{ij}\lambda,\leavevmode\nobreak\
\\\
\mathcal{G}^{(1)}_{n+1,j+1}(\lambda)=\gamma_{n+1,j}\xi_{2}(\lambda)+\varrho_{n+1,j}\xi_{3}(\lambda)-\pi_{n+1,j}\lambda,$
(31)

where $\displaystyle\vartheta_{ij}$, $\displaystyle\varsigma_{ij}$ and
$\displaystyle\tau_{ij}$ are given by (26) and $\displaystyle(\ref{gamma})$,
for $\displaystyle i=1,\cdots,n+1$, $\displaystyle j=1,n-1$.

  * 3)

Since $\displaystyle\mathcal{G}$ is a block matrix, it is useful to define the
matrix $\displaystyle\left(G^{(2)}_{l,r}\right)$ for $\displaystyle
l,r=1,\cdots,n-1$ from matrix
$\displaystyle\left(\mathcal{G}^{(2)}_{i,j+1}\right)$, with rows from
$\displaystyle 2$ to $\displaystyle n$ and columns from $\displaystyle 2$ to
$\displaystyle n$.

$G^{(2)}_{l,r}=\mathcal{G}^{(2)}_{l+1,r+1}(\lambda)\quad\text{ for
}l,r=1,\cdots,n-1.$ (32)

From $\displaystyle(\ref{pms})$ and $\displaystyle\eqref{gge}$, we have
$\displaystyle\det(A-\lambda
B)=\mathcal{G}_{1,1}(\lambda)\det({G}^{(2)}_{l,r})\nu_{n+1}$

From $\displaystyle\mathcal{G}_{1,1}(\lambda)=0$ we get
$\displaystyle\xi_{1}(\lambda)=0$, therefore $\displaystyle(\ref{xis}.a)$
holds. To determine the corresponding eigenvalue, we substitute
$\displaystyle\lambda$ by $\displaystyle\lambda_{s}$ into
$\displaystyle(\ref{pms})$ and the statement $\displaystyle(\ref{xis}.b)$
holds. Finally, the eigenpair $\displaystyle(\lambda_{s},\vec{r}_{s})$ given
by $\displaystyle(\ref{lambdasa})$ is obtained. In this eigenpair, only
saturation varies and we call this wave a saturation wave, or Buckley-Leverett
type wave.

The other eigenvalues are obtained by solving

$\det({G}^{(2)}_{l,r})=0,$ (33)

where the matrix $\displaystyle{G}^{(2)}_{l,r}$ is given by (32).

To determine the other eigenvalues, we substitute $\displaystyle\lambda$ by
$\displaystyle\lambda_{\Lambda_{i}}$ given by $\displaystyle(\ref{eig2a})$
into $\displaystyle(\ref{gge})$. Notice that

$\xi_{2}(\lambda_{\Lambda_{i}})=uf-\varphi\lambda_{\Lambda_{i}}s=u\frac{f(s-\Lambda_{i})-(f-\Lambda_{i})s}{s-\Lambda_{i}}=u\Lambda_{i}\frac{s-f}{s-\Lambda_{i}},$
(34)
$\xi_{3}(\lambda_{\Lambda_{i}})=uf_{o}-\varphi\lambda_{\Lambda_{i}}s_{o}=u(1-f)-\lambda_{\Lambda_{i}}(1-s)=u(1-\Lambda_{i})\frac{s-f}{s-\Lambda_{i}}.$
(35)

Substituting $\displaystyle\xi_{2}(\lambda_{\Lambda_{i}})$ and
$\displaystyle\xi_{3}(\lambda_{\Lambda_{i}})$ given by (34)-(35) into (32),
the matrix $\displaystyle{G}^{(2)}_{l,r}$ can be rewritten as

$\displaystyle\displaystyle{G}^{(2)}_{l,r}=\frac{\lambda_{\Lambda_{i}}(s-f)}{f-\Lambda_{i}}\begin{bmatrix}\vartheta_{2,1}+\Lambda_{i}(\vartheta_{2,1}-\varsigma_{2,1})&\dots&\vartheta_{2,n-1}+\Lambda_{i}(\vartheta_{2,n-1}-\,\varsigma_{2,n-1})\\\
\vdots&\ddots&\vdots\\\
\vartheta_{n,1}+\Lambda_{i}(\vartheta_{n,1}-\,\varsigma_{n,1})&\dots&\vartheta_{n,n-1}+\Lambda_{i}(\vartheta_{n,n-1}-\,\varsigma_{n,n-1})\end{bmatrix}$
$\displaystyle\displaystyle-\frac{\lambda_{\Lambda_{i}}(s-f)}{f-\Lambda_{i}}\begin{bmatrix}\frac{f-\Lambda_{i}}{s-f}\tau_{2,1}&\cdots&\frac{f-\Lambda_{i}}{s-f}\tau_{2,n-1}\\\
\vdots&\ddots&\vdots\\\
\frac{f-\Lambda_{i}}{s-f}\tau_{n,1}&\cdots&\frac{f-\Lambda_{i}}{s-f}\tau_{n,n-1}\\\
\end{bmatrix}.$

Let us denote

$\mathcal{A}=\begin{bmatrix}\vartheta_{2,1}&\dots&\vartheta_{2,n-1}\\\
\vdots&\ddots&\vdots\\\ \vartheta_{n,1}&\dots&\vartheta_{n,n-1}\end{bmatrix},$
(36)
$\mathcal{B}=-\begin{bmatrix}(\vartheta_{2,1}-\varsigma_{2,1})&\dots&(\vartheta_{2,n-1}-\,\varsigma_{2,n-1})\\\
\vdots&\ddots&\vdots\\\
(\vartheta_{n,1}-\,\varsigma_{n,1})&\dots&\vartheta_{n,n-1}-\,\varsigma_{n,n-1})\end{bmatrix},$
(37)

$\mathcal{A}_{1}=-\begin{bmatrix}\frac{f}{s-f}\tau_{2,1}&\cdots&\frac{f}{s-f}\tau_{2,n-1}\\\
\vdots&\ddots&\dots\\\
\frac{f}{s-f}\tau_{n,1}&\dots&\frac{f}{s-f}\tau_{n,n-1}\\\ \end{bmatrix},$
(38)

and

$\mathcal{B}_{1}=-\begin{bmatrix}\frac{-1}{s-f}\tau_{2,1}&\cdots&\frac{-1}{s-f}\tau_{2,n-1}\\\
\vdots&\ddots&\dots\\\
\frac{-1}{s-f}\tau_{n,1}&\dots&\frac{-1}{s-f}\tau_{n,n-1}\\\ \end{bmatrix}.$
(39)

Let us denote
$\displaystyle\varUpsilon=(\lambda_{\Lambda_{i}}(s-f)/(f-\Lambda_{i}))$, we
have

${G}^{(2)}_{l,r}=\varUpsilon((\mathcal{A}-\Lambda_{i}\mathcal{B})+(\mathcal{A}_{1}-\Lambda_{i}\mathcal{B}_{1})),$
(40)

where the matrices $\displaystyle\mathcal{A}$ and $\displaystyle\mathcal{B}$
depend only on the variables $\displaystyle y$, while the matrices
$\displaystyle\mathcal{A}_{1}$ and $\displaystyle\mathcal{B}_{1}$ depend on
$\displaystyle y$ and $\displaystyle s$.

Thus to calculate the eigenvalues $\displaystyle\lambda_{\Lambda_{i}}$ it is
sufficient to determine $\displaystyle\Lambda_{i}$ using the generalized
eigenvalue problem given in (19).

  * 4)

The right eigenvectors $\displaystyle\vec{r}_{\Lambda_{i}}$ related to
$\displaystyle\lambda_{\Lambda_{i}}$ are obtained by substituting
$\displaystyle\lambda$ by $\displaystyle\lambda_{\Lambda_{i}}$ in (30).
Because of the structure of the matrix in (30), we can split the calculation
of $\displaystyle\vec{r}_{\Lambda_{i}}$. First, we obtain the coordinates
$\displaystyle(r_{\Lambda_{i}}^{2},\ldots,r_{\Lambda_{i}}^{n})$ of eigenvector
$\displaystyle\vec{r}_{\Lambda_{i}}$. Let us define the auxiliary vector
$\displaystyle\vec{v}_{i}$ of size $\displaystyle n-1$, which solves

$\left({G}^{(2)}_{l,r}\right)\vec{v}_{i}=0.$ (41)

Using $\displaystyle(\ref{ueq1})$-$\displaystyle(\ref{ueq1n})$, we can see
that $\displaystyle(\ref{pmsnok})$, after simplifications is written in the
form $\displaystyle(\ref{eigeqs2})$. Second, the coordinates $\displaystyle
r_{\Lambda_{i}}^{n+1}$ and $\displaystyle r^{1}_{\Lambda_{i}}$ are calculated
by solving the first and the last equations of $\displaystyle(\ref{pms})$ as

$\displaystyle\displaystyle
r_{\Lambda_{i}}^{n+1}=-\frac{\displaystyle\sum_{j=2}^{n}\mathcal{G}^{(1)}_{n+1,j}(\lambda_{\Lambda_{i}})v_{i}^{j-1}}{\mathcal{G}_{n+1,n+1}},$
(42) $\displaystyle\displaystyle
r^{1}_{\Lambda_{i}}=-\frac{\mathcal{G}_{1,n+1}r_{n+1}+\displaystyle\sum_{j=2}^{n}\mathcal{G}_{1,j}(\lambda_{\Lambda_{i}})v_{i}^{j-1}}{\mathcal{G}_{1,1}(\lambda_{\Lambda_{i}})}.$
(43)

and the coordinates $\displaystyle r_{\Lambda_{i}}^{j}$ for $\displaystyle
j=2,\cdots,n$ are given by

$r^{j}_{\Lambda_{i}}=v^{j-1}_{i},$ (44)

where $\displaystyle v_{i}$ is given in (41).

Here the matrices $\displaystyle\mathcal{G}$ and
$\displaystyle\mathcal{G}^{(1)}_{n+1,j}$ evaluated at
$\displaystyle\lambda=\lambda_{\Lambda_{i}}$ are given in (22) and (29).

Let us denote $\displaystyle a=(a_{1},\cdots,a_{n-1})$ and $\displaystyle
b=(b_{1},\cdots,b_{n-1})$ with

$a_{j}=-\mathcal{G}^{(1)}_{n+1,j}(\lambda_{\Lambda_{i}})/\mathcal{G}_{n+1,n+1},$
(45)

and

$b_{j}=\mathcal{G}_{1,n+1}\mathcal{G}^{(1)}_{n+1,j}(\lambda_{\Lambda_{i}})/(\mathcal{G}_{n+1,n+1}\mathcal{G}_{1,1})-\mathcal{G}_{1,j}(\lambda_{\Lambda_{i}})/\mathcal{G}_{1,1}.$
(46)

Let us denote the matrix $\displaystyle P=(P_{i,j})$ with $\displaystyle n+1$
rows and $\displaystyle n-1$ columns as

$P_{i,j}=\delta_{i,1}a_{j}+\delta_{i,n+1}b_{j}+\delta_{i,j+1},$ (47)

with $\displaystyle i=1,\ldots,n+1$, $\displaystyle j=1\ldots,n-1$,
$\displaystyle\delta_{i,j}=1$ if $\displaystyle i=j$, and
$\displaystyle\delta_{i,j}=0$ otherwise.

Using (47) and (42)- (43) we obtain (20).$\displaystyle\square$

The integral curves $\displaystyle\mathcal{W}^{s}$ and
$\displaystyle\mathcal{W}^{\Lambda_{i}}$ associated to
$\displaystyle\vec{r}_{s}$ and $\displaystyle\vec{r}_{\Lambda_{i}}$ are
obtained by integrating the ODEs

$\frac{d(s,y,u)}{d\xi}=\vec{r}_{s},$ (48)

and

$\frac{d(s,y,u)}{d\xi}=\vec{r}_{\Lambda_{i}}.$ (49)

Finally, we conclude that in the Buckley-Leverett wave only the saturation
changes. Moreover, it is possible to verify the following

######  Corollary 2.

If $\displaystyle\lambda_{\Lambda}=u{f}/{\varphi s}$
($\displaystyle\Lambda_{i}=0$) or
$\displaystyle\lambda_{\Lambda}=u{(1-f)}/\varphi(1-s)$
($\displaystyle\Lambda_{i}=1$), then the field associated to
$\displaystyle\lambda_{\Lambda}$ is linearly degenerate.

##  4 Bifurcation surfaces

###  4.1 Resonance

Bifurcation loci are used to divide the phase space in subregions in which the
sequences of waves for the Riemann solutions are the same. To study these loci
we take the flux functions $\displaystyle f$ given in 2.

For each $\displaystyle\Lambda_{i}$ there exist a coincidence locus, which is
denoted by $\displaystyle\Gamma_{\Lambda_{i}}$; it occurs where the
eigenvalues $\displaystyle\lambda_{s}$ and
$\displaystyle\lambda_{\Lambda_{i}}$ coincide, i.e., at the zero of

$\mathcal{G}_{\Lambda_{i}}(s,y)=\frac{\partial f(s,y)}{\partial
s}-\left(\frac{f(s,y)-\Lambda_{i}(s,y)}{s-\Lambda_{i}(s,y)}\right).$ (50)

A relevant situation occurs when the system of conservation laws in (11) has
constant coefficients $\displaystyle\rho_{rj}$ with $\displaystyle
j=1,\cdots,n+1$, which we summarize in the corollary

######  Corollary 3.

In the system of conservation laws (11), if the coefficients
$\displaystyle\rho_{rj}$ with $\displaystyle j=1,\cdots,n+1$ are constant then
the eigenvalues $\displaystyle\Lambda_{i}$ in (19) depend only on the
variables $\displaystyle y$.

If any coefficient $\displaystyle\rho_{rj}$ are not constant then there occur
significant changes in the behavior of the bifurcation loci. In particular
notice that we cannot provide a characterization of the inflection locus
$\displaystyle\mathcal{G}_{\Lambda_{i}}(s,y)=0$ in the nonconstant case.
However making some assumptions on the coefficients of the system of
conservation laws (11), we have that each $\displaystyle\Gamma_{\Lambda_{i}}$
consists of the union of two disconnected surfaces.

######  Lemma 4.

Assume that the coefficients $\displaystyle\rho_{rj}$ with $\displaystyle
j=1,\cdots,n+1$ in (11) are constant. Assume that either
$\displaystyle\Lambda_{i}<0$ or $\displaystyle\Lambda_{i}>1$. We conclude
$\displaystyle\Gamma_{\Lambda_{i}}=\Gamma_{1,\Lambda_{i}}\cup\Gamma_{2,\Lambda_{i}}$,
with

$\displaystyle\displaystyle\Gamma_{1,\Lambda_{i}}=\\{(s_{\Lambda_{i}}^{1},y):\mathcal{G}_{\Lambda_{i}}(s_{\Lambda_{i}}^{1}(y),y)=0,\quad
s_{\Lambda_{i}}^{1}(y)<s^{*}\\},$ (51)
$\displaystyle\displaystyle\Gamma_{2,\Lambda_{i}}=\\{(s_{\Lambda_{i}}^{2},y):\mathcal{G}_{\Lambda_{i}}(s_{\Lambda_{i}}^{2}(y),y)=0,\quad
s_{\Lambda_{i}}^{2}(y)>s^{*}\\},$ (52)

where $\displaystyle s^{*}$ is given in Assumption S-shaped (1).

Proof: Since $\displaystyle 0<s<1$ and $\displaystyle\Lambda_{i}<0$ or
$\displaystyle\Lambda_{i}>1$ the eigenvalues
$\displaystyle\lambda_{\Lambda_{i}}$ are well defined. Notice that
$\displaystyle\frac{\partial\mathcal{G}_{\Lambda_{i}}(s,y)}{\partial
s}=\frac{\partial^{2}f(s,y)}{\partial
s^{2}}+\frac{\lambda_{\Lambda_{i}}-\lambda_{s}}{s-\Lambda_{i}}$. Using the
S-Shape assumption, we have
$\displaystyle\frac{\partial\mathcal{G}_{\Lambda_{i}}(s,y)}{\partial
s}\not=0$, thus by the implicit function Theorem (IFT) there exist
$\displaystyle s_{\Lambda_{i}}^{1}(y)$ such that
$\displaystyle\mathcal{G}_{\Lambda_{i}}(s_{\Lambda_{i}}^{1}(y),y)=0$.
Moreover, by the S-Shape assumption there exist $\displaystyle
s_{\Lambda_{i}}^{2}(y)$ satisfying Eq. $\displaystyle(\ref{gam1})$ of Lemma
4.$\displaystyle\Box$

The coincidence between two different $\displaystyle\lambda_{\Lambda_{i}}$
occurs in some special cases.

######  Lemma 5.

At the states where
$\displaystyle{\lambda}_{\Lambda_{i}}={\lambda}_{\Lambda_{j}}$, we have
$\displaystyle f=s$ or $\displaystyle\Lambda_{i}=\Lambda_{j}$.

###  4.2 Inflection loci

Other structures appearing in this model are the inflection loci. The
inflections are, generically, co-dimension 1 structures, in which the monotone
increase of characteristic speed fails, i.e.,
$\displaystyle\nabla\lambda\cdot\vec{r}=0$.

In this model, we have $\displaystyle n-1$ fields. For the field
$\displaystyle\lambda_{s}$, it is easy to see that
$\displaystyle\nabla\lambda_{s}\cdot\vec{r}_{s}=\frac{1}{\varphi}\displaystyle{\frac{\partial^{2}f}{\partial
s^{2}}}$, thus the inflection locus consists of the states $\displaystyle
s^{*}=s^{*}(y)$, satisfying
$\displaystyle\displaystyle{\frac{\partial^{2}f}{\partial
s^{2}}(s^{*},\cdot)=0}$. We denote the inflection for the field
$\displaystyle(\lambda_{s},\vec{r}_{s})$ as $\displaystyle\mathcal{I}_{s}$:

$\mathcal{I}_{s}=\left\\{(s^{*},y)\text{ : }\frac{\partial^{2}f}{\partial
s^{2}}(s^{*},y)=0\right\\}.$ (53)

From IFT, if $\displaystyle\displaystyle{\frac{\partial^{3}f}{\partial
s^{3}}(s^{*},y)\neq 0}$ for $\displaystyle(s^{*},y)\in\mathcal{I}_{s}$, then
$\displaystyle\mathcal{I}_{s}$ is a smooth hypersurface with co-dimension 1 in
$\displaystyle\Omega$. Moreover, if $\displaystyle f$ does not depend on
$\displaystyle y$, then $\displaystyle\mathcal{I}_{s}$ is a hyperplane in
$\displaystyle\Omega$ with constant $\displaystyle s$.

######  Proposition 6.

The inflection locus $\displaystyle I_{\Lambda_{i}}$ of field
$\displaystyle(\lambda_{\Lambda_{i}},\vec{r}_{\Lambda_{i}})$ satisfies

$\nabla\lambda_{\Lambda_{i}}\cdot\vec{r}_{\Lambda_{i}}=\frac{1}{\varphi}\frac{s-f}{(s-\Lambda_{i})^{2}}(\lambda_{s}-{\color[rgb]{0,0,0}{\lambda_{\Lambda_{i}}}})\mathcal{H}_{i}(s,y),\quad\text{
where }\quad\mathcal{H}_{i}=\vec{w_{i}}\cdot\vec{v_{i}},$ (54)

with $\displaystyle\vec{w_{i}}=(\vec{w_{i}})_{j}$, for $\displaystyle
j=1,\cdots,n-1$, and $\displaystyle(\vec{w_{i}})_{j}$ given by
$\displaystyle(\ref{wee})$. Here $\displaystyle\vec{v_{i}}$ is a right
eigenvector associated to generalized eigenvalue $\displaystyle\Lambda_{i}$ in
(19).

Proof of Proposition 6

Let $\displaystyle\lambda_{\Lambda_{i}}$ be given by
$\displaystyle(\ref{eig2a})$. Calculating
$\displaystyle\nabla\lambda_{\Lambda_{i}}$, we obtain

$\displaystyle\displaystyle\nabla\lambda_{\Lambda_{i}}=\frac{1}{\varphi}\left(\omega,\varpi,\frac{\lambda_{\Lambda_{i}}}{u}\right),\text{
where }$
$\displaystyle\displaystyle\omega=\frac{1}{s-\Lambda_{i}}\left(\left(\lambda_{s}-\lambda_{\Lambda_{i}}\right)+u\frac{\partial\Lambda_{i}}{\partial
s}(\lambda_{\Lambda_{i}}-1)\right),$ (55)
$\displaystyle\displaystyle\varpi_{j}=-\frac{u}{(s-\Lambda_{i})^{2}}(f-s)\frac{\partial\Lambda_{i}}{\partial
y_{j}}+\frac{\partial f}{\partial y_{j}}\frac{1}{s-\Lambda_{i}},\quad\text{
for }j=1,\cdots,n-1.$ (56)

Using (40) and $\displaystyle\vec{r}_{\Lambda_{i}}$ given in
$\displaystyle(\ref{eiget1})$, we calculate
$\displaystyle\nabla\lambda_{\Lambda_{i}}\cdot\vec{r}_{\Lambda_{i}}$. After
some calculations we obtain $\displaystyle(\ref{fint})$, for which
$\displaystyle(\vec{w_{i}})_{j}$, for $\displaystyle j=1,\cdots,n-1$ are given
by

$\displaystyle\displaystyle(\vec{w_{i}})_{j}$
$\displaystyle\displaystyle=\omega(s-\Lambda_{i})^{2}a_{j}\delta_{1,j}+\frac{(s-\Lambda_{i})^{2}\lambda_{\Lambda_{i}}b_{j}}{u}\delta_{n+1,j}+(s-\Lambda_{i})^{2}\sum_{k=2}^{n}\delta_{k,j}\varpi_{j}.$
(57)

where $\displaystyle a_{j}$ and $\displaystyle b_{j}$ are given in (45) and
(46) $\displaystyle\square$.

Let us define the hypersurfaces $\displaystyle\mathcal{C}_{s}$,
$\displaystyle\Gamma_{\Lambda_{i}}$ and
$\displaystyle\mathcal{J}_{\mathcal{H}_{i}}$ as the loci where $\displaystyle
s=f$, $\displaystyle\lambda_{s}-\lambda_{\Lambda_{i}}=0$ and
$\displaystyle\mathcal{H}_{i}(s,y)=0$, respectively. From Proposition 6 the
inflection locus $\displaystyle\mathcal{I}_{\Lambda_{i}}$ is given by the
union of above hypersurfaces. A similar decomposition of the inflection
surface was obtained by [22].

##  5 Discontinuous solutions

Discontinuous solutions of systems of conservation laws
$\displaystyle(\ref{eq1})$ satisfy the Rankine-Hugoniot condition, i.e., for
given left and right states $\displaystyle(s^{-},y^{-},u^{-})$ and
$\displaystyle(s^{+},y^{+},u^{+})$ we have

$u^{+}F_{i}((s^{+},y^{+})-u^{-}F_{i}(s^{-},y^{-})=\sigma(G_{i}(s^{+},y^{+})-G_{i}(s^{-},y^{-})),$
(58)

with $\displaystyle i=1,\dots,n+1$, where $\displaystyle\sigma$ is the
propagation speed of the discontinuity. Equation (58) can be rewritten as

$\displaystyle\displaystyle\Phi_{i}.[\sigma,u^{+},u^{-}]=0,$ (59)

where $\displaystyle\Phi_{i}=(\Phi_{i1},\Phi_{i2},\Phi_{i3})$, with
$\displaystyle\Phi_{i1}=s^{+}(\rho_{wi}^{+}-\rho_{oi}^{+})+\rho_{oi}^{+}-(s^{-}(\rho_{wi}^{-}-\rho_{oi}^{-})+\rho_{oi}^{-})+(\rho_{ri}^{+}-\rho_{ri}^{+})$,
$\displaystyle\Phi_{i2}=-(f^{+}(\rho_{wi}^{+}-\rho_{oi}^{+})+\rho_{oi}^{+})$,
$\displaystyle\Phi_{i3}=(f^{-}(\rho_{wi}^{-}-\rho_{oi}^{-})+\rho_{oi}^{-})$
and $\displaystyle\rho^{+}=\rho(y^{+})$, $\displaystyle\rho^{-}=\rho(y^{-})$,
$\displaystyle f^{+}=f(s^{+},y^{+})$, $\displaystyle f^{-}=f(s^{-},y^{-})$
with $\displaystyle\rho=(\rho_{wi},\rho_{oi},\rho_{ri})$.

For each fixed state $\displaystyle(s^{-},y^{-})$, the Hugoniot locus
$\displaystyle\mathcal{HL}(s^{-},y^{-})$ consists of all states
$\displaystyle(s^{+},y^{+})$ satisfying $\displaystyle(\ref{rank1})$. In [14]
it is proven that the Hugoniot locus satisfies

$\mathcal{HL}(s^{-},y^{-})=\left\\{(s^{+},y^{+}):det(\Phi_{i}^{T},\Phi_{k}^{T},\Phi_{j}^{T})=0\right\\},$
(60)

for all combinations of distinct indices
$\displaystyle\\{i,k,j\\}\in\\{1,\ldots,n+1\\}$. Also, we verify that instead
of considering all index combination it is sufficient to reduce system (60) to
$\displaystyle n-1$ linearly independent equations. This result can summarized
as follows: let $\displaystyle i_{1}$ and $\displaystyle i_{2}$
$\displaystyle\in\\{1,2,\cdots,n+1\\}$ be two indices such that
$\displaystyle\Phi_{i_{1}}$ and $\displaystyle\Phi_{i2}$ are linearly
independent, then Equation (60) reduces to

$\mathcal{HL}(s^{-},y^{-})=\left\\{(s^{+},y^{+}):det(\Phi_{k}^{T},\Phi_{i_{1}}^{T},\Phi_{i_{2}}^{T})=0\right\\},$
(61)

for values of $\displaystyle k\in\\{1,\ldots,n+1\\}$ distinct from
$\displaystyle i_{1}$ and $\displaystyle i_{2}$. The curve given by (61) is
the intersection of $\displaystyle n-1$ independent surfaces in the
$\displaystyle n$ dimensional space $\displaystyle(s,y)$.

Assuming that there exist two indices $\displaystyle i,j$ with $\displaystyle
i\neq j$ such that the denominator is different from zero for given states
$\displaystyle(s^{-},y^{-})$ and $\displaystyle(s^{+},y^{+})$ then the values
of $\displaystyle u^{+}$ can be obtained from (58) as

$\frac{u^{+}}{u^{-}}=\frac{F_{i}(s^{-},y^{-})[G_{j}]-F_{j}(s^{-},y^{-})[G_{i}]}{F_{i}(s^{+},y^{+})[G_{j}]-F_{j}(s^{+},y^{+})[G_{i}]},$
(62)

where $\displaystyle[G_{i}]=G_{i}(s^{+},y^{+})-G_{i}(s^{-},y^{-})$ and
$\displaystyle[G_{j}]=G_{j}(s^{+},y^{+})-G_{j}(s^{-},y^{-})$.

A particular branch of the Hugoniot locus can be characterized in the
following

######  Lemma 7.

For fixed $\displaystyle(s^{-},y^{-})$, a branch of
$\displaystyle\mathcal{HL}(s^{-},y^{-})$ consists of the states of the form
$\displaystyle(s,y)\in\Omega$, with $\displaystyle s$ variable and
$\displaystyle y=y^{-}$.

This branch is called Buckley-Leverett or saturation branch.

##  6 Elementary waves and interactions

The features of the coefficients $\displaystyle\rho_{rj}$
$\displaystyle(j=1,\cdots,n+1)$ in the accumulation function have implications
on the quantitative behavior of the Riemann solutions. If such coefficients
are constant in the phase space then the matrices
$\displaystyle\mathcal{A}_{1}$ and $\displaystyle\mathcal{B}_{1}$ defined in
(38) and (39) satisfy $\displaystyle\mathcal{A}_{1}\equiv 0$ and
$\displaystyle\mathcal{B}_{1}\equiv 0$. This fact implies that the matrix
$\displaystyle{G}^{(2)}_{l,r}$ in (40) and the corresponding eigenvalues and
eigenvectors of equation (19) depend only on the variables $\displaystyle y$.
In [14] this case is studied and the existence of the splitting in the
variables $\displaystyle u$ and $\displaystyle s$ from $\displaystyle y$ to
obtain the chemical rarefaction and the Hugoniot locus is shown. This
splitting is extremely useful to obtain the wave curves.

In general to consider the coefficients $\displaystyle\rho_{rj}$
$\displaystyle(j=1,\cdots,n+1)$ as non constant functions in the accumulation
terms in (13) has significant implications in the structure of the solution.
In such a case it is no longer possible to solve the Riemann problem in terms
of the variables $\displaystyle y$ alone. Thus, considering surface complexes
increases difficulties because more wave interactions appear in the
construction of wave curves.

However, in general in the state variables $\displaystyle(s,y,u)$, the
variable $\displaystyle u$ can be split from the variables
$\displaystyle(s,y)$. Let us take $\displaystyle\Lambda$ as any of the
eigenvalues $\displaystyle\Lambda_{i}$ ($\displaystyle i=1,\cdots,n-1$) given
in (19). Using the Proposition 1 and equation $\displaystyle(\ref{integral})$,
we obtain along the integral curves $\displaystyle\mathcal{W}^{\Lambda}$ the
unknown variables $\displaystyle(s,y)$ independently of $\displaystyle u$ by
solving the system of differential equations:

$\frac{ds}{d\xi}=r_{\Lambda}^{1}(s,y),\;\frac{dy_{i}}{d\xi}=r_{\Lambda}^{i}(s,y),\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ i=2,\cdots,n,\leavevmode\nobreak\
\leavevmode\nobreak\ \frac{du}{d\xi}=ug_{\Lambda}^{n+1}(s,y),$ (63)

with $\displaystyle\xi=x/t$,
$\displaystyle\xi^{-}=\lambda(s^{-},y^{-},u^{-})$, $\displaystyle s(0)=s^{-}$,
$\displaystyle y(0)=y^{-}$, $\displaystyle u(0)=u^{-}$ and $\displaystyle
g_{\Lambda}^{n+1}(s,y)=r_{\Lambda}^{n+1}(s,y,u)/u$, where $\displaystyle
r_{\Lambda}^{i}$ ($\displaystyle i=1,\cdots,n+1$) are given in (42). Here the
superscript index is used to indicate the coordinates of the vector
$\displaystyle\vec{r}_{\Lambda}$ given in Eq. $\displaystyle(\ref{eig2a})$ of
Proposition 1. After obtaining $\displaystyle(s(\xi),y(\xi))$, we use the
expression from the last equation of $\displaystyle(\ref{equvv})$ yielding

$u=u^{-}exp({\gamma(\xi)}),\leavevmode\nobreak\ \leavevmode\nobreak\
\gamma(\xi)=\int_{\xi^{-}}^{\xi}{g_{\Lambda}^{n+1}}(s(\eta),y(\eta))d\eta,$
(64)

where $\displaystyle\xi=x/t$, $\displaystyle\xi^{-}=\lambda(s^{-},y^{-})$ and
$\displaystyle u^{-}$ is the initial value of $\displaystyle u$ on the
rarefaction wave, i.e. $\displaystyle u=u^{-}$ at $\displaystyle\xi=\xi^{-}$.
Furthermore, in [13], the authors proved that it is possible to project the
Rankine-Hugoniot locus from a state $\displaystyle(s^{-},y^{-},u^{-})$ in the
$\displaystyle\Omega=\\{(s,y)\\}$ space. In this sense, we can solve the
Riemann problem in the projection space $\displaystyle(s,y)$, which reduces
the number of calculations.

From formula (64), we can deduce that if the coefficients of system (11) are
such that the function $\displaystyle\gamma(\xi)$ assumes small values then
the values of the variable $\displaystyle u$ is small. This case was observed
in experiment studies, see [1].

We consider the interaction of the waves. There are two cases: namely, waves
in the same family, and waves in different families. The resonance condition
leads to the construction of curves formed by structures from distinct
families.

Here, we have $\displaystyle n-1$ kinds of shocks, one associated to
$\displaystyle\lambda_{s}$, which we denote by $\displaystyle\mathcal{H}_{s}$
and the other associated to $\displaystyle\lambda_{\Lambda_{i}}$
($\displaystyle i=1,\cdots,n-1$), i.e.,
$\displaystyle\mathcal{H}_{\Lambda_{i}}$. Similarly to
$\displaystyle\mathcal{R}_{s}$, the curve $\displaystyle\mathcal{H}_{s}$ is a
straight line parallel to the $\displaystyle s$ axis. On the other hand, for
states in the neighborhood of $\displaystyle(s-,y-,u-)$, the chemical shocks
$\displaystyle\mathcal{H}_{\Lambda_{i}}$ have structures similar to the
rarefaction $\displaystyle\mathcal{R}_{\Lambda_{i}}$. However, for states far
from the left state, we do not identify the behavior of Rankine-Hugoniot
locus.

In summary, the elementary waves consist of

  * •

The B-L saturation wave $\displaystyle\mathcal{R}_{s}$ where only the
saturation varies. This wave curve stops at the inflection locus
$\displaystyle I_{s}$ given in (53) and bifurcates at surface
$\displaystyle\lambda_{s}=\lambda_{\Lambda_{i}}$ with $\displaystyle
i=1,\ldots,n-1$,

  * •

The shock wave curve $\displaystyle\mathcal{H}_{s}$ consists of the locus
where only the saturation varies,

  * •

The saturation wave curves $\displaystyle\mathcal{R}_{\Lambda_{i}}$
($\displaystyle i=1,\ldots,n-1$) associated to the couple
$\displaystyle(\lambda_{\Lambda_{i}},r_{\Lambda_{i}})$. These wave curves stop
at the inflection loci $\displaystyle I_{\Lambda_{i}}$ given in (54) and
bifurcate at surfaces $\displaystyle\Gamma_{ij}$:
$\displaystyle\lambda_{\Lambda_{j}}=\lambda_{\Lambda_{i}}$ with $\displaystyle
i,j=1,\ldots,n$ $\displaystyle i\neq j$ or
$\displaystyle\lambda_{s}=\lambda_{\Lambda_{i}}$ with $\displaystyle
i=1,\ldots,n$,

  * •

The composite wave curve $\displaystyle
C_{s}=\mathcal{R}_{s}\cup\mathcal{H}_{s}$ consists of a characteristic B-L
shock,

  * •

The composite wave curves $\displaystyle
C_{\Lambda_{i}}=\mathcal{R}_{\Lambda_{i}}\cup\mathcal{H}_{\Lambda_{i}}$, where
$\displaystyle\mathcal{H}_{\Lambda_{i}}$ is a characteristic shock associated
to the branch of the family $\displaystyle\lambda_{\Lambda_{i}}$ with
$\displaystyle i=1,\ldots,n-1$,

  * •

The constant state $\displaystyle\mathcal{C}$.

When the Riemann problem is solved the collision of waves of the same and
different families may happen. Thus the possible collisions can be summarized
as

  * •

The collision of two shocks:
$\displaystyle\mathcal{H}_{s}\mathcal{H}_{\Lambda_{i}}$,
$\displaystyle\mathcal{H}_{\Lambda_{i}}\mathcal{H}_{\Lambda_{j}}$ with
$\displaystyle i\neq j$, $\displaystyle i,j=1,\cdots,n$ ,

  * •

The collision of two rarefaction waves:
$\displaystyle\mathcal{R}_{s}\mathcal{R}_{\Lambda_{i}}$,
$\displaystyle\mathcal{R}_{\Lambda_{i}}\mathcal{R}_{\Lambda_{j}}$ with
$\displaystyle i\neq j$, $\displaystyle i,j=1,\cdots,n-1$,

  * •

The collision of rarefaction waves and a shock:
$\displaystyle\mathcal{R}_{s}\mathcal{H}_{\Lambda_{i}}$,
$\displaystyle\mathcal{R}_{\Lambda_{i}}\mathcal{H}_{\Lambda_{j}}$ with
$\displaystyle i\neq j$, $\displaystyle i,j=1,\cdots,n-1$,

  * •

The collision of a shock and rarefaction waves:
$\displaystyle\mathcal{H}_{s}\mathcal{R}_{\Lambda_{i}}$,
$\displaystyle\mathcal{H}_{\Lambda_{i}}\mathcal{R}_{\Lambda_{j}}$ with
$\displaystyle i\neq j$, $\displaystyle i,j=1,\cdots,n-1$.

Overtaking of elementary waves from the same family

  * •

i-shock overtakes another i-shock
$\displaystyle\mathcal{H}_{\Lambda_{i}}\mathcal{H}_{\Lambda_{i}}$,

  * •

s-shock overtakes another s-shock
$\displaystyle\mathcal{H}_{s}\mathcal{H}_{s}$,

  * •

s-rarefaction wave overtakes i-shock
$\displaystyle\mathcal{H}_{\Lambda_{i}}\mathcal{R}_{s}$,

  * •

s-shock overtakes i-rarefaction wave
$\displaystyle\mathcal{R}_{\Lambda_{i}}\mathcal{H}_{s}$,

  * •

i-shock wave overtakes i-rarefaction
$\displaystyle\mathcal{R}_{\Lambda_{i}}\mathcal{H}_{\Lambda_{i}}$,

  * •

s-shock wave overtakes s-rarefaction
$\displaystyle\mathcal{R}_{s}\mathcal{H}_{s}$,

  * •

i-shock overtakes i-rarefaction wave
$\displaystyle\mathcal{R}_{\Lambda_{i}}\mathcal{H}_{\Lambda_{i}}$,

  * •

s-rarefaction wave overtakes another s-rarefaction wave
$\displaystyle\mathcal{R}_{s}\mathcal{R}_{s}$.

##  7 The Riemann Problem: Structure of waves

The wave curves and bifurcation loci studied in the previous Sections are the
main ingredients used to obtain the solution of $\displaystyle(\ref{eqt})$ for
given initial data

$\begin{cases}(s,y,u)_{L}&\text{ for }x<0,\\\ (s,y,\cdot)_{R}&\text{ for
}x>0.\end{cases}$ (65)

For the system 1 it is necessary to provide the initial data $\displaystyle u$
either on the right or left side. The missing value can be found from the
other state and the wave sequences in the solution.

The Riemann solution consists of a sequence of elementary waves $\displaystyle
w_{k}$ (shocks and rarefactions) for $\displaystyle k=1,2,\cdots,m$ and
constant states $\displaystyle U_{k}$ for $\displaystyle k=1,2,\cdots,p$, in
which $\displaystyle p$ is not known a priori. When the system is strictly
hyperbolic and each field is genuinely nonlinear or linearly degenerate, for a
system of $\displaystyle n$ equations, we can prove that we have a sequence of
at most $\displaystyle n$ waves, see [8]. However, for the general case (with
loss of strict hyperbolicity or loss of genuine nonlinearity) it is not easy
to determine a priori the number of waves in the solution. Here, we show that
the system exhibits several nonlinear effects such as loss of hyperbolicity.

In any case, these waves are written as

$\mathcal{U}_{L}\equiv\mathcal{U}_{0}\overset{w_{1}}{\longrightarrow}\mathcal{U}_{1}\overset{w_{2}}{\longrightarrow}\cdots\overset{w_{m}}{\longrightarrow}\mathcal{U}_{m}\equiv\mathcal{U}_{R},$
(66)

where $\displaystyle\mathcal{U}=(s,y,u)$. In the Riemann solution it is
necessary that the waves have increasing speed, the so called geometrical
compatibility. Sometimes, this geometrical compatibility is sufficient to
furnish existence and uniqueness for the solution. Moreover, this condition is
used to select the physical sequence of waves for the Riemann solution.

With the previous considerations, we study the behavior of waves and their
interactions. In that direction, we study the Riemann solution for a
particular problem, where we are able to obtain some interesting behavior of
the wave curves, when these wave curves cross the inflection or coincidence
loci.

In the previous Sections we find bifurcations, eigenvalues and eigenvectors
for system (11). In summary, the eigenpairs for states $\displaystyle(s,y)$,
$\displaystyle y=(y_{1},\cdots,y_{n-1})$ are: a saturation eigenpair of form
$\displaystyle\displaystyle{\lambda_{s}=\frac{u}{\varphi}\frac{\partial
f}{\partial s}\quad\vec{r}_{s}=(1,0,\cdots,0)}$ and $\displaystyle n-1$
compositional eigenvalues of form
$\displaystyle\displaystyle{\lambda_{\Lambda}=\frac{u}{\varphi}\frac{f-\Lambda}{s-\Lambda}}$.
We will also use the fact that each eigenvector generates its own rarefaction
curves, which are denoted $\displaystyle\mathcal{R}_{s}$ for the field
$\displaystyle\vec{r}_{s}$ and $\displaystyle\mathcal{R}_{\Lambda_{i}}$ for
the field $\displaystyle\vec{r}_{\Lambda_{i}}$.

For the rarefaction $\displaystyle\mathcal{R}_{s}$ only the saturation
changes. The rarefactions associated to
$\displaystyle\mathcal{R}_{\Lambda_{i}}$ are in general transverse to
$\displaystyle\mathcal{R}_{s}$.

In this model, we assume there is a coincidence surface between eigenvalues of
different families, i.e., there is $\displaystyle(s,y)$ for which
$\displaystyle\lambda_{s}=\lambda_{\Lambda_{i}}$. This system loses strict
hyperbolicity on this coincidence surface.

###  7.1 Application

In this section, we present numerical examples for a system of four
conservation laws having the same structure as system (1). We give an overview
of the main parameters of the model, a full description of which can be found
in [1]. The physical model includes oil recovery improvement due to the
injection of carbonated water. The model also includes enhancement of the
solubility of carbon dioxide in the water due to low salinity. We have four
state variables. Two of these variables are the water saturation
$\displaystyle s$ and the Darcy velocity $\displaystyle u$. The two other
state variables are associated to the water phase, namely the hydrogen ion
concentration $\displaystyle[H^{+}]$ (denoted here by $\displaystyle y_{2}$)
and the chloride concentration $\displaystyle[Cl^{-}]$ (denoted here by
$\displaystyle y_{1}$).

The system has four state variables and four equations. These equations
include the molar concentrations of the total hydrogen
$\displaystyle(\rho_{w1})$, the total chloride $\displaystyle(\rho_{w2})$, the
total hydrogen minus twice the total oxygen $\displaystyle(\rho_{w3})$ with
carbon dioxide $\displaystyle(\rho_{o3})$ and decane
$\displaystyle(\rho_{o4})$ (see [1]). They also include the two coefficients
$\displaystyle\rho_{r3}$ and $\displaystyle\rho_{r4}$ in the accumulation term
representing the formation of chemical complexes on the rock surface.

Disregarding the diffusive terms as in [1], the accumulation $\displaystyle
G=(G_{1},G_{2},G_{3},G_{4})$ and flux functions $\displaystyle F=u\widehat{F}$
(with
$\displaystyle\widehat{F}=(\widehat{F}_{1},\widehat{F}_{2},\widehat{F}_{3},\widehat{F}_{4})$)
are written as:

$\displaystyle\displaystyle G$ $\displaystyle\displaystyle=(\varphi
s\rho_{w1},\varphi
s\rho_{w2},s\varphi\rho_{w3}+4\varphi(1-s)\rho_{o3}+(1-\varphi)\rho_{r3},(1-s)\varphi\rho_{o4}+(1-\varphi)\rho_{r4}),$
(67) $\displaystyle\displaystyle\widehat{F}$
$\displaystyle\displaystyle=(f\rho_{w1},f\rho_{w2},f\rho_{w3}+4(1-f)\rho_{o3},(1-f)\rho_{o4}).$
(68)

Here, $\displaystyle\varphi$ is the porosity. The coefficient functions
$\displaystyle\rho_{w1}$, $\displaystyle\rho_{w3}$, $\displaystyle\rho_{o3}$,
$\displaystyle\rho_{w2}$ and $\displaystyle\rho_{o4}$ are obtained through the
PHREEQC program (see [19, 20]). We use the coefficients given in [1]. These
coefficients depend on the variables $\displaystyle y=(y_{1},y_{2})$ only.
Here, we consider several possibilities for the coefficients
$\displaystyle\rho_{r3}$ and $\displaystyle\rho_{r4}$.

The fractional flow function $\displaystyle f$ has the $\displaystyle S$-shape
with properties described in Assumption 1. Here, we use the function presented
in Section 2.

We denote by $\displaystyle\lambda_{k}(A)$, with $\displaystyle
k=s,\Lambda_{1},\Lambda_{2}$ the characteristic speeds at state $\displaystyle
A$. Also we denote by $\displaystyle\sigma(A,B)$ the shock speed from state
$\displaystyle A$ to state $\displaystyle B$.

In order to construct the Riemann solution in a region where
$\displaystyle\lambda_{\Lambda_{1}}<\lambda_{s}<\lambda_{\Lambda_{2}}$, the
geometrical compatibility condition requires that we start with a rarefaction
wave $\displaystyle\mathcal{R}_{\Lambda_{1}}$, followed by either a constant
state or a shock wave (see below), which must have a speed larger than the
speed of the previous wave. On the other hand, if
$\displaystyle\lambda_{\Lambda_{2}}>\lambda_{\Lambda_{1}}>\lambda_{s}$, we
start with a saturation wave; either a rarefaction
$\displaystyle\mathcal{R}_{s}$ or a shock $\displaystyle H_{s}$.

Finally, we can initiate building the solution of the Riemann problem starting
from the injection conditions of the reservoir (left state) until arriving at
the initial conditions of the reservoir (right state).

An in-house Riemann solver was developed to calculate and represent wave
curves, which satisfy the compatibility and admissibility criteria (see [2]
for technical implementation details). In this solver, we automated the
construction of slow and fast solution paths taking into account the
bifurcation structures. Rarefaction curves are obtained by solving ODE
systems, one for B-L rarefactions wave curves given in (48) and the other for
chemical wave curves in (49). Moreover, shock waves are obtained by finding
the intersection of the two surfaces constructed using (61) for $\displaystyle
n=4$. Additional composition wave curves are necessary for crossing the
inflection surfaces. When two eigenvalues coincide the rarefaction waves
collide and numerical difficulties appear for the continuation of the wave
curve. In this case the desingularizing algorithm developed in [2] is used to
allow wave curves to cross resonance surfaces.

####  7.1.1 Riemann solutions when the coefficients $\displaystyle\rho_{rj}$
are constant.

Inflection and Coincidence Surfaces

Rarefaction wave solutions suffer structural modifications along the loci
where the eigenvalues coincide. The coincidence surface has four connected
parts given as the sets where speeds satisfy either
$\displaystyle\lambda_{s}=\lambda_{\Lambda_{1}}$ (two parts) or
$\displaystyle\lambda_{\Lambda_{2}}=\lambda_{\Lambda_{1}}$ (one part) or
$\displaystyle\lambda_{s}=\lambda_{\Lambda_{2}}$ (one part) (see Figure 1).

Figure 1: Coincidence Surfaces
$\displaystyle\lambda_{\Lambda_{i}}=\lambda_{s}$, i=1,2 and
$\displaystyle\lambda_{\Lambda_{i}}=\lambda_{\Lambda_{j}}$ for $\displaystyle
i\neq j$. Each of these surfaces have almost constant saturation.

From Proposition 6, we verify that the inflection surface consists of four
connected parts: one part for $\displaystyle\lambda_{s}$, two parts for
$\displaystyle\lambda_{\Lambda_{1}}$ and one part for
$\displaystyle\lambda_{\Lambda_{2}}$ (see Figure 2).

Figure 2: Inflection surfaces for three eigenpairs
$\displaystyle(\lambda_{s},\vec{r}_{s})$,
$\displaystyle(\lambda_{\Lambda_{1}},\vec{r}_{\Lambda_{1}})$ and
$\displaystyle(\lambda_{\Lambda_{2}},\vec{r}_{\Lambda_{2}})$.

Riemann Solutions

In the general case, despite building elementary waves involves simple
constructions, solving the Riemann problem is a difficult task due to the
large number of combinations of elementary waves as well as to the presence of
several bifurcation surfaces. However, when the coefficients
$\displaystyle\rho_{jr}$, in the accumulation function in (67), either are
constant or have small derivatives relative to the other coefficients some
useful simplifications occur.

One of such simplifications is that the field associated to the eigenvalue
$\displaystyle\lambda_{\Lambda_{2}}=({u}/{\varphi})({f}/{s})$ is linearly
degenerate (see Corollary 2) and as a consequence the rarefaction and shock
curves coincide. Moreover, curves associated to
$\displaystyle\lambda_{\Lambda_{1}}$ and $\displaystyle\lambda_{s}$ are
invariant in the planes $\displaystyle y_{1}=(y_{1})_{C}$ with constant
$\displaystyle(y_{1})_{C}$, while curves associated to
$\displaystyle\lambda_{\Lambda_{2}}$ are transversal to these planes. This
permits the displacement in the state space in the direction of the variable
$\displaystyle y_{1}$.

Let us to construct the solution for the Riemann-Goursat problem associated to
the system of conservation laws $\displaystyle(\ref{eq1})$ with flux and
accumulation functions given in (67)-(68), that is the solution of these
equations with piecewise constant initial data

$\left\\{\begin{array}[c]{ll}L=(s_{L},(y_{1})_{L},(y_{2})_{L},u_{L})&\textbf{if}\hskip
5.69046ptx<0,\\\ R=(s_{R},(y_{1})_{R},(y_{2})_{R},\;\cdot\;)&\textbf{if}\hskip
5.69046ptx>0.\end{array}\right.$ (69)

To do that, we consider two examples taking into account the elementary waves
and the bifurcation surfaces described above. In the first example the left
state $\displaystyle L$ and the right state $\displaystyle R$ have different
values for the variable $\displaystyle y_{1}$, i.e., $\displaystyle
y_{1}=C_{1}$ and $\displaystyle y_{1}=C_{2}$, with constant $\displaystyle
C_{1}$ and $\displaystyle C_{2}$. However, in the second example the left and
right state have the same value for the variable $\displaystyle y_{1}$.
Because the variables $\displaystyle(s,y_{1},y_{2})$ can be split from the
variable $\displaystyle u$, we are able to represent the Riemann solutions in
the projected state space $\displaystyle(s,y_{1},y_{2})$ in
$\displaystyle\\{[0,1]\times[0,0.4]\times[2.7,4]\\}$ (see [14]). In the two
examples we fixed $\displaystyle\rho_{3r}=0.01$ and
$\displaystyle\rho_{4r}=0.02$.

In both examples the left state $\displaystyle L$ belongs to the region where
$\displaystyle\nabla\Lambda_{s}\cdot\vec{r}_{s}<0$ and $\displaystyle f>s$ and
the right state $\displaystyle R$ belongs to the region where
$\displaystyle\nabla\lambda_{s}\cdot\vec{r}_{s}>0$ and $\displaystyle f<s$.
These restrictions imply that to cross the surface
$\displaystyle\nabla\lambda_{s}\cdot\vec{r}_{s}=0$, we need only type B-L
shocks. This is because $\displaystyle\nabla\lambda_{s}\cdot\vec{r}_{s}=0$ is
a barrier for the B-L rarefactions. Moreover, the surface $\displaystyle f=s$
constitutes a barrier for the remaining rarefaction curves (see [14]). This
fact can be verified by simple inspection of the formula for the infection
surface given in (54).

In the first example, we take the left state $\displaystyle
L=(s_{L},(y_{1})_{L},(y_{2})_{L})=(0.99,0.03,2.74)$ belonging to a region
where the condition
$\displaystyle\lambda_{s}<\lambda_{\Lambda_{1}}<\lambda_{\Lambda_{2}}$ holds.
We take the right state given by $\displaystyle
R=(s_{R},(y_{1})_{R},(y_{2})_{R})=(0.27,0.3,2.74)$. In this case left and
right states belong to distinct planes i.e. $\displaystyle
y_{1}=(y_{1})_{L}=0.03$ and $\displaystyle y_{1}=(y_{1})_{R}=0.3$.

The solution of the Riemann problem can be obtained using the following wave
curves. First, we take a saturation curve $\displaystyle\mathcal{R}_{s}(L)$
from $\displaystyle L$ to the state $\displaystyle A$ where
$\displaystyle\lambda_{s}(A)=\lambda_{\Lambda_{1}}(A)$ (See Figure 4). From
the state $\displaystyle A$, we take a rarefaction curve
$\displaystyle\mathcal{R}_{\Lambda_{1}}(A)$ connected to the physical
boundary.

From here the continuation of curve
$\displaystyle\mathcal{R}_{\Lambda_{1}}(A)$ requires the determination of a
point belonging to both curves $\displaystyle\mathcal{R}_{\Lambda_{2}}(B^{*})$
and $\displaystyle\mathcal{H}_{s}(R)$ defined below. Since seeking such a
point in the three-dimensional space is a numerically unstable procedure, the
construction of an auxiliary surface $\displaystyle\varUpsilon$ is necessary
(see Figure 3). Thus in this example, we construct a surface
$\displaystyle\Upsilon$ made up of rarefaction curves
$\displaystyle\mathcal{R}_{\Lambda_{2}}(B^{*})$ starting at all points
$\displaystyle B^{*}$ on the rarefaction
$\displaystyle\mathcal{R}_{\Lambda_{1}}(A)$. In this way, such curves form a
surface which connects the rarefaction curve
$\displaystyle\mathcal{R}_{\Lambda_{1}}(A)$ with the plane
$\displaystyle\mathcal{R}$ ($\displaystyle y_{1}=0.3$), determining the
intersection curve $\displaystyle\Phi_{c}$.

Figure 3: Here we show the auxiliary surface $\displaystyle\Upsilon$ made up
of rarefaction curves of type $\displaystyle\mathcal{R}_{\Lambda_{2}}(B^{*})$
starting at all points $\displaystyle B^{*}$ on the rarefaction
$\displaystyle\mathcal{R}_{\Lambda_{1}}(A)$ (shortened
$\displaystyle\mathcal{R}_{\Lambda_{1}}$ in the Figure). The plane
$\displaystyle\mathcal{R}$ corresponds to $\displaystyle y_{2}=0.3$. Here
$\displaystyle\Phi_{c}$ represents the intersection curve between the surface
$\displaystyle\Upsilon$ and the plane $\displaystyle\mathcal{R}$.

Now, we take a backward B-L shock curve out of the right state $\displaystyle
R=(0.27,0.3,2.74)$ on which the characteristic speed decreases. Afterwards, we
look for the intersection state $\displaystyle C$ of the auxiliary surface
$\displaystyle\Upsilon$ with the backward shock curve
$\displaystyle\mathcal{H}_{s}(R)$ from $\displaystyle R$ to $\displaystyle C$
(see Figure 4). The existence of the state $\displaystyle C$ is guaranteed by
the facts that in the curve $\displaystyle\mathcal{H}_{s}(R)$ only the
saturation variable changes and the auxiliary surface $\displaystyle\Upsilon$
has an intersection curve $\displaystyle\Phi_{c}$ with the plane
$\displaystyle\mathcal{R}$ ($\displaystyle y_{1}=0.3$), which passes through
the point with the same coordinate $\displaystyle y_{2}=2.74$ of the right
state $\displaystyle R$.

Now, we look for the rarefaction $\displaystyle\mathcal{R}_{\Lambda_{2}}(B)$
on the surface $\displaystyle\Upsilon$ from the point $\displaystyle B$ to
$\displaystyle C$ (See Figure 4), where the state $\displaystyle B$ is chosen
such that the curve $\displaystyle\mathcal{R}_{\Lambda_{2}}(B)$ crosses the
plane $\displaystyle\mathcal{R}$, $\displaystyle y_{1}=0.3$ at the point
$\displaystyle C$ with coordinate $\displaystyle y_{2}=2.74$. Here, we use the
fact that the curve $\displaystyle\mathcal{R}_{\Lambda_{2}}(B)$ is a contact
wave, therefore the characteristic speed along this curve remains constant
equal to $\displaystyle\lambda_{\Lambda_{1}}(B)$. This fact guarantees
compliance to the monotone increase of the wave speed.

Figure 4: Schematics of Riemann solution corresponding to left (L) and right
states (R). The solution is obtained as a concatenation of four wave curves
separated by two constant states: the first consists of a rarefaction
saturation wave $\displaystyle\mathcal{R}_{s}(L)$ followed by a rarefaction
curve $\displaystyle\mathcal{R}_{\Lambda_{1}}(A)$ (shortened
$\displaystyle\mathcal{R}_{\Lambda_{1}}$ in the Figure), the second is a
rarefaction curve $\displaystyle\mathcal{R}_{\Lambda_{2}}(B)$ (shortened
$\displaystyle\mathcal{R}_{\Lambda_{2}}$ in the Figure) and the third is an
admissible BL-shock $\displaystyle\mathcal{H}_{s}(C)$. The coordinates of each
point are $\displaystyle L=(0.99,0.03,2.74)$, $\displaystyle
A=(0.77,0.03,2.74)$, $\displaystyle B=(0.70,0.03,2.77)$, $\displaystyle
C=(0.70,0.3,2.74)$ and $\displaystyle R=(0.15,0.03,2.74)$. Denoting
$\displaystyle\lambda=(\Lambda_{s},\lambda_{\Lambda_{1}},\lambda_{\Lambda_{2}})$,
we have $\displaystyle\lambda(L)=(1.15\times 10^{-5},2.66\times
10^{-5},2.73\times 10^{-5})$, $\displaystyle\lambda(A)=(2.14\times
10^{-5},2.14\times 10^{-5},3.36\times 10^{-5})$,
$\displaystyle\lambda(B)=(4.76\times 10^{-5},2.47\times 10^{-5},3.36\times
10^{-5})$, $\displaystyle\lambda(C)=(5.26\times 10^{-5},2.45\times
10^{-5},3.42\times 10^{-5})$, $\displaystyle\lambda(R)=(2.55\times
10^{-6},2.58\times 10^{-7},3.04\times 10^{-5})$ and
$\displaystyle\sigma(C,L)=5.26\times 10^{-5})$.

Since at the state $\displaystyle C$ the wave speed satisfies
$\displaystyle\lambda_{\Lambda_{2}}(C)=\sigma(C,R)$, we have at this point a
characteristic BL-shock $\displaystyle\mathcal{H}_{s}(C)$ from state
$\displaystyle C$ to $\displaystyle R$. Finally, the solution is given by the
wave sequence $\displaystyle
L\xrightarrow{\mathcal{R}_{s}}(A)\xrightarrow{\mathcal{R}_{\Lambda_{1}}}(B)\xrightarrow{\mathcal{R}_{\Lambda_{2}}}(C)\xrightarrow{\mathcal{H}_{s}}(R)$.

Similarly, it can be verified that the wave sequences $\displaystyle
L\xrightarrow{\mathcal{R}_{s}}(A)\xrightarrow{\mathcal{R}_{\Lambda_{1}}}(B)\xrightarrow{\mathcal{H}_{s}}(R^{*})$
solve the Riemann problem from the left state $\displaystyle L$ to all right
states $\displaystyle R^{*}$ given by $\displaystyle R^{*}=(s,0.3,2.74)$ where
$\displaystyle s\in(0.15,0.26)$ and the shock speeds satisfy
$\displaystyle\sigma(C,R^{*})\in(4.13\times 10^{-5},5.11\times 10^{-5})$. We
verify that $\displaystyle\sigma(C,R^{*})>\lambda_{\Lambda_{2}}(C)$ holds for
all states $\displaystyle R^{*}$. Moreover, the shocks from state
$\displaystyle C$ to $\displaystyle R^{*}$ are 3-Lax shocks (see [15, 16]). We
verify that at each intermediate state of the wave curves the geometrical
compatibility condition holds. We also check the possibility of the existence
of other combinations of curves from the left to right states, but these cases
do not satisfy such condition. However, we do not have a theoretical result
which guarantees that the solution found constitutes the unique solution for
this problem for chosen $\displaystyle L$ and $\displaystyle R$.

In the second example, we take the left state $\displaystyle
L=(s_{L},(y_{1})_{L},(y_{2})_{L})=(0.99,0.03,2.74)$ satisfying
$\displaystyle\lambda_{s}(L)<\lambda_{\Lambda_{1}}(L)<\lambda_{\Lambda_{2}}(L)$.
Assuming that the right state is given by $\displaystyle
R=(s_{R},(y_{1})_{R},(y_{2})_{R})=(0.33,0.03,2.85)$, with $\displaystyle L$
and $\displaystyle R$ belonging to the same plane
$\displaystyle\\{y_{1}=0.03\\}$, the solution of the Riemann problem can be
obtained with the following sequence of curves. First, from $\displaystyle L$
to the state $\displaystyle A$ where
$\displaystyle\lambda_{s}(A)=\lambda_{\Lambda_{1}}(A)$, we take the saturation
wave $\displaystyle R_{s}(L)$ (see Figure 5). From the state $\displaystyle
A$, we take the rarefaction $\displaystyle\mathcal{R}_{\Lambda_{1}}(A)$
connecting to the physical boundary.

To reach the right state $\displaystyle R$ from here, it is useful to
construct an extension of the rarefaction $\displaystyle R_{\Lambda_{1}}(A)$.
Such extension curve is defined as those points $\displaystyle X$ belonging to
a shock between states $\displaystyle X$ and $\displaystyle B$ in
$\displaystyle R_{\Lambda_{1}}(A)$ such that
$\displaystyle\sigma(X,B)=\lambda_{s}(B)$.

In this way, we find an intermediate state $\displaystyle B$ such that
$\displaystyle\sigma(B,R)=\lambda_{s}(B)$ which consists of a characteristic
B-L shock to the point $\displaystyle R$. Finally, the solution is given by
the wave sequence $\displaystyle
L\xrightarrow{\mathcal{R}_{s}}(A)\xrightarrow{\mathcal{R}_{\Lambda_{1}}}(B)\xrightarrow{\mathcal{H}_{s}}(R)$.
We verify that at the intermediate states $\displaystyle A$ and $\displaystyle
B$ the compatibility condition of increasing speed is satisfied.

Figure 5: Schematics of Riemann solution corresponding to left (L) and right
states (R). The solution is obtained from a concatenation of three wave curves
separated by two constant states: the first consists of a rarefaction
saturation curve $\displaystyle R_{s}(L)$ (shortened
$\displaystyle\mathcal{R}_{s}$ in the Figure) followed by a compatible
rarefaction curve $\displaystyle\mathcal{R}_{\Lambda_{1}}(A)$ (shortened
$\displaystyle\mathcal{R}_{\Lambda_{1}}$ in the Figure), the second one is the
rarefaction curve $\displaystyle\mathcal{R}_{\Lambda_{1}}(B)$ and the third
one is an admissible BL-shock $\displaystyle\mathcal{H}_{s}(C)$. The
coordinates of each point are $\displaystyle L=(0.99,0.03,2.74)$,
$\displaystyle A=(0.77,0.03,2.74)$, $\displaystyle B=(0.68,0.03,2.85)$ and
$\displaystyle R=(0.33,0.03,2.85)$. Let us denote
$\displaystyle\lambda=(\lambda_{s},\lambda_{\Lambda_{1}},\lambda_{\Lambda_{2}})$,
we have $\displaystyle\lambda(L)=(1.15\times 10^{-5},2.66\times
10^{-5},2.73.15\times 10^{-5})$, $\displaystyle\lambda(A)=(2.14\times
10^{-5},2.14\times 10^{-5},3.36\times 10^{-5})$,
$\displaystyle\lambda(B)=(6.01\times 10^{-5},2.70\times 10^{-5},3.28\times
10^{-5})$, $\displaystyle\lambda(R)=(7.96\times 10^{-6},1.02\times
10^{-6},2.7\times 10^{-5})$ and $\displaystyle\sigma(C,L)=6.01\times
10^{-5})$.

Other calculations were made using nonconstant coefficients
$\displaystyle\rho_{jr}$ depending on the variables $\displaystyle y_{1}$ and
$\displaystyle y_{2}$. In such examples changes of the intermediate states are
observed but there are no changes with the respect to the constant
coefficients case in the structure of the wave curves. However, these
calculations are not sufficient to elucidate the issue of structural stability
of the solutions. This issue will be studied in future work.

##  8 Conclusion

In this paper the eigenvalue analysis of a system of $\displaystyle n+1$
conservation laws is performed. A method to build the discontinuous solutions
and the bifurcation surfaces is presented as well. These elementary waves and
surfaces are the basis for the subsequent analysis of the solution paths. It
is noteworthy that these solution paths are useful for verification of
numerical methods.

An overview of possible wave collisions is presented. Using the elementary
waves and the bifurcation surfaces some Riemann solutions are found. Due to
the complexity of the problem in general, some hypotheses on the coefficients
are made to obtain algorithms that allow easy application of the wave curve
method. Structural stability of the solution of this kind of system is an open
question.

## Acknowledgement

The authors thank Sergio Pilotto as well as support from CAPES under grant
88881.156518/2017-01, CAPES/NUFFIC grant 88887.156517/2017-00, CNPq under
grants 405366/2021-3, 306566/2019-2, FAPERJ under grants E- 26/210.738/2014,
E-26/202.764/2017, E-26/201.159/2021.

## References

  * [1] A.C. Alvarez, T. Blom, W.J. Lambert, J. Bruining, and D. Marchesin.  Analytical and numerical validation of a model for flooding by saline carbonated water.  Journal of Petroleum Science and Engineering, 167:900–917, 2018\.
  * [2] A.C. Alvarez, G.T. Goedert, and D. Marchesin.  Resonance in rarefaction and shock curves: local analysis and numerics of the continuation method.  Journal of Hyperbolic Differential Equations, 17(04):639–676, 2020\.
  * [3] A.C. Alvarez, W.J. Lambert, J. Bruining, and D. Marchesin.  The Riemann solution for carbonated waterflooding.  In ECMOR XV-15th European Conference on the Mathematics of Oil Recovery, 2016.
  * [4] A. V. Azevedo and D. Marchesin.  Multiple viscous solutions for systems of conservation laws.  Transactions of the American Mathematical Society, 347(8):3061–3077, 1995.
  * [5] T. Blom, A.C. Alvarez, W.J. Lambert, D. Marchesin, and J. Bruining.  Low salinity carbonated waterflooding.  In ECMOR XV-15th European Conference on the Mathematics of Oil Recovery, 2016.
  * [6] R. H. Brooks and A. T. Corey.  Hydraulic properties of porous media and their relationship to drainage design.  PhD thesis, Colorado State University. Libraries, 1963.
  * [7] Se. E. Buckley and MCi. Leverett.  Mechanism of fluid displacement in sands.  Transactions of the AIME, 146(01):107–116, 1942.
  * [8] C. M. Dafermos.  Hyperbolic Conservation Laws in Continuum Physics, volume 3.  Springer, 2005.
  * [9] W. Dahmen, S. Müller, and A. Voß.  Riemann problem for the Euler equation with non-convex equation of state including phase transitions.  In Analysis and Numerics for Conservation Laws, pages 137–162. Springer, 2005.
  * [10] E. L. Isaacson, D. Marchesin, C. F. Palmeira, and B. J. Plohr.  A global formalism for nonlinear waves in conservation laws.  Communications in Mathematical Physics, 146(3):505–552, 1992.
  * [11] Larry W Lake, Steven Lawrence Bryant, and Aura N Araque-Martinez.  Geochemistry and Fluid Flow, volume 7.  Gulf Professional Publishing, 2002.
  * [12] W. Lambert, D. Marchesin, and J. Bruining.  The Riemann solution for the injection of steam and nitrogen in a porous medium.  Transport in Porous Media, 81(3):505–526, 2010.
  * [13] W. J. Lambert and D. Marchesin.  The Riemann problem for multiphase flows in porous media with mass transfer between phases.  Journal of Hyperbolic Differential Equations, 6(04):725–751, 2009\.
  * [14] W.J. Lambert, A.C. Alvarez, V. Matos, D. Marchesin, and J. Bruining.  Nonlinear wave analysis of geochemical injection for multicomponent two phase flow in porous media.  Journal of Differential Equations, 266(1):406–454, 2019.
  * [15] P. D. Lax.  Hyperbolic systems of conservation laws ii.  Communications on Pure and Applied Mathematics, 10(4):537–566, 1957.
  * [16] Peter D. Lax.  Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves.  SIAM, 1973.
  * [17] T. P. Liu.  The Riemann problem for general 2$\displaystyle\times$ 2 conservation laws.  Transactions of the American Mathematical Society, 199:89–112, 1974\.
  * [18] T. P. Liu.  The Riemann problem for general systems of conservation laws.  Journal of Differential Equations, 18(1):218–234, 1975.
  * [19] D. L. Parkhurst and C.A.J. Appelo.  Description of input and examples for PHREEQC version 3: a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations.  Technical report, US Geological Survey, 2013.
  * [20] D. L. Parkhurst, C.A.J. Appelo, et al.  User’s guide to PHREEQC (version 2): A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations.  1999\.
  * [21] G. A. et. al. Pope.  The application of fractional flow theory to enhanced oil recovery.  Society of Petroleum Engineers Journal, 20(03):191–205, 1980.
  * [22] H. Wahanik.  Thermal Effects in the Injection of CO2 in Deep Underground Aquifers.  PhD thesis, PhD thesis, IMPA. Brazil, 2011.
  * [23] B. Wendroff.  The Riemann problem for materials with nonconvex equations of state I: Isentropic flow.  Journal of Mathematical Analysis and Applications, 38(2):454–466, 1972.