The Cauchy Problem for a non Strictly Hyperbolic $3\times 3$ System of Conservation Laws Arising in Polymer Flooding

We consider a simple model for polymer flooding in two phase flow in rough media

$$\begin{cases}\partial_{t}s+\partial_{x}f(s,c,k)&=~{}0,\\ \partial_{t}[cs]+\partial_{x}[cf(s,c,k)]&=~{}0,\\ \partial_{t}k&=~{}0,\end{cases}$$

(1.1)

associated with the initial data

$$\left(s,c,k\right)(0,x)=\left(\bar{s},\bar{c},\bar{k}\right)(x),\quad x\in\mathbb{R}.$$

(1.2)

Here, the unknown vector is $(s,c,k)$, where $s$ is the saturation of water phase, $c$ is the fraction of polymer dissolved in water, and $k$ denotes the permeability of the porous media. We see that $k$ does not change in time, $k(t,x)=\bar{k}(x)$, and the initial data $\bar{k}(\cdot)$ might be discontinuous.

We neglect both the adsorption of polymers in the porous media and the gravitational force, where the solution to the Riemann problem becomes more complex. For such Riemann solvers, see [JohWin, MR3663611] for the effect of the adsorption term, and [WSCauchy] for the addition of the gravitational force term. In particular, when the adsorption effect is included, the $c$ family described below would no longer be linearly degenerate, while adding the gravitational force term, the $s$ waves described below could have negative speed. Both effects would disrupt the carefully designed wave front tracking algorithm we use to prove the main theorem.

The conserved quantities and their fluxes are given by, respectively

$$\mathbf{G}\left(s,c,k\right)=\begin{pmatrix}s\\ cs\\ k\end{pmatrix},\qquad\mathbf{F}\left(s,c,k\right)=\begin{pmatrix}f\left(s,c,k\right)\\ cf\left(s,c,k\right)\\ 0\end{pmatrix}.$$

Denoting the three families as the $s$, $c$ and $k$ family, we have the following 3 eigenvalues as functions of the variables $\left(\sigma,\gamma,\kappa\right)$ in the $\left(s,c,k\right)$ space.

$$\lambda_{s}=\partial_{\sigma}f\left(\sigma,\gamma,\kappa\right),\qquad\lambda_{c}=\frac{f\left(\sigma,\gamma,\kappa\right)}{\sigma},\qquad\lambda_{k}=0,$$

and the three corresponding right eigenvectors (in the $\left(s,c,k\right)$ space):

$$r_{s}=\begin{pmatrix}1\\ 0\\ 0\end{pmatrix},\qquad r_{c}=\begin{pmatrix}-\partial_{\gamma}f\left(\sigma,\gamma,\kappa\right)\\ \partial_{\sigma}f\left(\sigma,\gamma,\kappa\right)-\frac{f\left(\sigma,\gamma,\kappa\right)}{\sigma}\\ 0\end{pmatrix},\qquad r_{k}=\begin{pmatrix}-\partial_{\kappa}f\left(\sigma,\gamma,\kappa\right)\\ 0\\ \partial_{\sigma}f\left(\sigma,\gamma,\kappa\right)\end{pmatrix}.$$

A straight computation shows that both the $c$ and $k$ families are linearly degenerate. Furthermore, there exist regions in the domain such that $\lambda_{s}=\lambda_{c}$ or $\lambda_{s}=\lambda_{c}=\lambda_{k}$, where the system is parabolic degenerate.

The flux function $f\left(\sigma,\gamma,\kappa\right)$ has the following properties. For any given $(\gamma,\kappa)$, the mapping $\sigma\mapsto f\left(\sigma,\gamma,\kappa\right)$ is the well-known S-shaped Buckley-Leverett function [BL] with a single inflection point, see Fig. 1. To be specific, we have

$$f(\sigma,\gamma,\kappa)\in[0,1],\qquad\partial_{\sigma}f(\sigma,\gamma,\kappa)\geq 0,\qquad\mbox{for all }(\sigma,\gamma,\kappa),$$

and, for all $(\gamma,\kappa)$,

$$\begin{split}&f(0,\gamma,\kappa)=0,\qquad\quad~{}f(1,\gamma,\kappa)=1,\\ &\partial_{\sigma}f(0,\gamma,\kappa)=0,\qquad~{}\partial_{\sigma}f(1,\gamma,\kappa)=0,\\ &\partial_{\sigma\sigma}f(0,\gamma,\kappa)>0,\qquad\partial_{\sigma\sigma}f(1,\gamma,\kappa)<0.\end{split}$$

(1.3)

Remark that conditions (1.3) guarantee that the eigenvalues and the eigenvectors written above are well defined (can be extended) when $\sigma=0$. For each given $(\gamma,\kappa)$, there exists a unique $\sigma^{*}\left(\gamma,\kappa\right)\in\left]0,1\right[$ such that

$$\partial_{\sigma\sigma}f(\sigma^{*}\left(\gamma,\kappa\right),\gamma,\kappa)=0.$$

A detailed analysis of the wave properties for this system is carried out in [WS3x3], with the following highlights:

• •

  $k$ waves are the slowest with speed $0$. Both $f$ and $c$ are continuous across any $k$ wave;

• •

  $c$ waves travel with non negative speed. Both $\frac{f}{s}$ and $k$ are continuous across any $c$ wave;

• •

  $s$ waves travel with positive speed. Both $c$ and $k$ are continuous across any $s$ wave.

In [WS3x3], the global Riemann solver is constructed. Here we give a brief summary. Let $(s_{l},c_{l},k_{l})$ and $(s_{r},c_{r},k_{r})$ be the left and right state of a Riemann problem, respectively. In general, the solution of the Riemann problem consists of a $k$ wave, a $c$ wave and possibly some $s$ waves. They can be constructed as follows.

• •

  Let $(s_{m},c_{l},k_{r})$ denote the right state of the $k$ wave. The value $s_{m}$ is uniquely determined by the condition

  $$f(s_{m},c_{l},k_{r})=f(s_{l},c_{l},k_{l}).$$

• •

  For the remaining waves, we have $k\equiv k_{r}$ throughout. We then solve the Riemann problem for the $2\times 2$ sub-system

  $$\partial_{t}s+\partial_{x}f(s,c,k_{r})=0,\qquad\partial_{t}(cs)+\partial_{x}(cf(s,c,k_{r}))=0$$

  (1.4)

  with Riemann data $(s_{m},c_{l})$ and $(s_{r},c_{r})$ as the left and right states. The solution consists of waves with non-negative speed.

We now give a precise definition of weak solution to the Cauchy problem (1.1)–(1.2) and state the main theorem.

We emphasize the fact that $s=0$ is not excluded in our theorem since we do not make use of Lagrangian coordinates which would have required $s>0$. Indeed, under the hypotheses $s\geq s^{*}>0$, $k(t,x)=\text{const.}$, system (1.1) is equivalent to its Lagrangian formulation [Wagner]:

$$\begin{cases}\partial_{t}\left(\frac{1}{s}\right)-\partial_{y}\left(\frac{f\left(s,c,k\right)}{s}\right)&=0,\\ \partial_{t}c&=0,\\ k&=\text{const.},\end{cases}$$

(1.5)

where $y$ is the Lagrangian coordinate satisfying $\partial_{x}y=s$, $\partial_{t}y=-f(s,c,k)$. Therefore, under some additional hypotheses, the result in [BGS], which holds for scalar equations since it is based on the maximum principle for Hamilton–Jacobi equation, could be used to prove the existence of a unique vanishing viscosity solution to the first equation in (1.5) and hence a (in some sense) unique solution to system (1.5). However, since we consider the case where $s$ can become 0, the analysis in [BGS] cannot be applied. Instead, we need to solve the original non triangular $2\times 2$ system in Eulerian coordinates. Furthermore, since we consider rough permeability function $k$, the corresponding system in the Lagrangian coordinate is no longer triangular,

$$\begin{cases}\partial_{t}\left(\frac{1}{s}\right)-\partial_{y}\left(\frac{f\left(s,c,k\right)}{s}\right)&=0,\\ \partial_{t}c&=0,\\ \partial_{t}\left(\frac{k}{s}\right)-\partial_{y}\left(\frac{kf\left(s,c,k\right)}{s}\right)&=0.\end{cases}$$

Remark that the Lagrangian coordinates introduced in [BGS, Section 6] for (1.1) make the system triangular, but still require $s\geq s^{*}>0$ and moreover they mix time and space, therefore the Cauchy problem for (1.1) in Eulerian coordinates is not equivalent to the Cauchy problem in the coordinates introduced in [BGS, Section 6].

In this paper, the proof for the existence of solution is carried out by showing that wave front tracking approximate solutions are compact by a compensated compactness argument (see for instance [KarRasTad] for an application of compensated compactness to a $2\times 2$ bi–dimensional related model).

The remaining of the paper is organized as follows. In Section 2 the wave front tracking approximate solutions are constructed. In Section 3 we prove the necessary entropy estimates. Finally in Section 4 the compensated compactness argument is carried out to prove Theorem 1.1.

In this section we modify the algorithm constructed in [CocliteRisebro] and [WSCauchy] to adapt it to system (1.1). We define the functions (see Figure 1):

$$g\left(\sigma,\gamma,\kappa\right)=\frac{f\left(\sigma,\gamma,\kappa\right)}{\sigma},\qquad P\left(\sigma,\gamma,\kappa\right)=\int_{0}^{\sigma}\left|\partial_{\xi}g\left(\xi,\gamma,\kappa\right)\right|d\xi.$$

(2.6)

Since at $\sigma=0$ both $f$ and its derivative $\partial_{\sigma}f$ vanish, we define $g(0,\gamma,\kappa)=0$. Hypotheses (1.3) imply that $\partial_{\sigma}g\left(0,\gamma,\kappa\right)>0$ and that the function $\sigma\mapsto g\left(\sigma,\gamma,\kappa\right)$ has one single maximum point somewhere between the single inflexion point of $f$ and the point $\sigma=1$. The function $P$ is continuous with respect to its three variables and strictly increasing and invertible with respect to the variable $\sigma$.

Fixing initial data $\left(\bar{s},\bar{c},\bar{k}\right)$ satisfying the hypotheses of Theorem 1.1 and fixing an approximation parameter $\varepsilon>0$, we can construct piecewise constant approximate initial data $\left(\bar{s}^{\varepsilon},\bar{c}^{\varepsilon},\bar{k}^{\varepsilon}\right)$ with values in $\left[0,1\right]^{3}$ such that:

	$\displaystyle\left\|\bar{k}^{\varepsilon}-\bar{k}\right\|_{\mathbf{L}^{\infty}}<\varepsilon,$	$\displaystyle\operatorname{Tot.Var.}\bar{k}^{\varepsilon}\leq\operatorname{Tot.Var.}\bar{k},$
	$\displaystyle\left\|\bar{c}^{\varepsilon}-\bar{c}\right\|_{\mathbf{L}^{\infty}}<\varepsilon,$	$\displaystyle\operatorname{Tot.Var.}\bar{c}^{\varepsilon}\leq\operatorname{Tot.Var.}\bar{c},$
	$\displaystyle\left\|\bar{s}^{\varepsilon}-\bar{s}\right\|_{\mathbf{L}^{1}\left(\left(-\frac{1}{\varepsilon},\frac{1}{\varepsilon}\right),\mathbb{R}\right)}\leq\varepsilon,$			(2.7)

Let $\bar{x}_{1},\ldots,\bar{x}_{N}$ be the set of points in which $\bar{k}^{\varepsilon}$ has jumps such that

$$\bar{k}^{\varepsilon}\left(x\right)=k_{0}\chi_{]-\infty,\bar{x}_{1}]}+\sum_{i=1}^{N-1}k_{i}\chi_{]\bar{x}_{i},\bar{x}_{i+1}]}(x)+k_{N}\chi_{]\bar{x}_{N},+\infty[},$$

and let $\bar{y}_{1},\ldots,\bar{y}_{M}$ be the set of points in which $\bar{c}^{\varepsilon}$ has jumps such that

$$\bar{c}^{\varepsilon}\left(x\right)=c_{0}\chi_{]-\infty,\bar{y}_{1}]}+\sum_{j=1}^{M-1}c_{j}\chi_{]\bar{y}_{j},\bar{y}_{j+1}]}(x)+c_{M}\chi_{]\bar{y}_{M},+\infty[}.$$

Without loss of generality, we can suppose that no $\bar{y}_{j}$ coincides with any $\bar{x}_{i}$. Define the constant

$$L=\left\lceil\frac{1}{\varepsilon}\sup_{\gamma,\kappa}\left\|f\left(\cdot,\gamma,\kappa\right)\right\|_{\mathbf{C}^{2}}\right\rceil\cdot\left(N+M\right),$$

where $\lceil\alpha\rceil$ denotes the least integer greater than or equal to the real number $\alpha$. In the following we denote by $\land$ the logical operator and. We consider the following finite sets of possible values for the function $g$:

	$\displaystyle\mathcal{G}_{0}^{1}$	$\displaystyle=\left\{g\mid g=g\left(\bar{s}^{\varepsilon}(x),\bar{c}^{\varepsilon}(x),\bar{k}^{\varepsilon}(x)\right),\quad x\in\mathbb{R}\right\},$
	$\displaystyle\mathcal{G}_{0}^{2}$	$\displaystyle=\left\{g\mid g=g\left(\frac{\ell}{L},c_{j},k_{i}\right),\quad i=0,\ldots,N,\;j=0,\ldots,M,\;\ell=0,\ldots,L\right\},$
	$\displaystyle\mathcal{G}_{0}^{3}$	$\displaystyle=\left\{g\mid g=\max_{0\leq\sigma\leq 1}g\left(\sigma,c_{j},k_{i}\right),\quad i=0,\ldots,N,\;j=0,\ldots,M\right\},$
	$\displaystyle\mathcal{G}_{0}^{4}$	$\displaystyle=\Big{\{}g\mid g=g\left(\sigma,c_{j},k_{i}\right)=g\left(\sigma,c_{j^{*}},k_{i^{*}}\right)\land g_{s}\left(\sigma,c_{j},k_{i}\right)\cdot g_{s}\left(\sigma,c_{j^{*}},k_{i^{*}}\right)<0,$
		$\displaystyle\qquad\qquad\text{ for some }i,i^{*}=0,\ldots,N,\;j,j^{*}=0,\ldots,M,\;\sigma\in[0,1]\Big{\}},$
	$\displaystyle\mathcal{G}_{0}$	$\displaystyle=\mathcal{G}_{0}^{1}\cup\mathcal{G}_{0}^{2}\cup\mathcal{G}_{0}^{3}\cup\mathcal{G}_{0}^{4}.$

We start the front tracking algorithm from the region $x<\bar{x}_{1}$. For this purpose we define the allowed values for $s$ in that region:

$$\mathcal{S}_{0,j}=\left\{\sigma\mid g\left(\sigma,c_{j},k_{0}\right)\in\mathcal{G}_{0}\right\}.$$

We call $f^{0,j}(\sigma)$ the linear interpolation of the map $\sigma\mapsto f\left(\sigma,c_{j},k_{0}\right)$ according to the points in $\mathcal{S}_{0,j}$. Observe that, since we have included $\mathcal{G}_{0}^{2}$ in $\mathcal{G}_{0}$, the set $\left\{\frac{\ell}{L}\right\}_{\ell=0}^{L}$ is included in $\mathcal{S}_{0,j}$ for every $j$, hence we have the uniform estimate

$$\begin{cases}\left|f\left(\sigma,c_{j},k_{0}\right)-f^{0,j}(\sigma)\right|\leq\varepsilon,\\ \left|\partial_{\sigma}f\left(\sigma,c_{j},k_{0}\right)-\partial_{\sigma}f^{0,j}(\sigma)\right|\leq\frac{\varepsilon}{N+M},\end{cases}\quad\text{ for all }\sigma\in\left[0,1\right],\ j=0,\ldots,M.$$

We solve all the Riemann problems in $x<\bar{x}_{1}$ at $t=0$ in the following way. Let $\bar{x}\in\left]\bar{y}_{j},\bar{y}_{j+1}\right[$ ($\bar{y}_{0}=-\infty$) be a jump in $\bar{s}^{\varepsilon}$. Here we take the entropic solution to the Riemann problem

$$\partial_{t}s+\partial_{x}f^{0,j}\left(s\right)=0,\qquad s\left(0,x\right)=\begin{cases}\bar{s}^{\varepsilon}\left(\bar{x}-\right)&\text{ for }x<\bar{x},\\ \bar{s}^{\varepsilon}\left(\bar{x}+\right)&\text{ for }x>\bar{x}.\end{cases}$$

Since $f^{0,j}$ is piecewise linear, the solution to the Riemann problem is piecewise constant and takes values in the set $\mathcal{S}_{0,j}$ of the kink points of $f^{0,j}$, moreover the entropy condition in [Bbook, Theorem 4.4] is satisfied. The same Riemann solver is used whenever at $t>0$ two $s$ waves interact in some region $\Omega_{0,j}$ defined below.

At the points $\bar{y}_{j}$, we solve the Riemann problem according to the minimum jump condition described in [WSCauchy] (see also [GimseRisebro]) that we briefly outline (see Fig 3).

Define

$$s^{L}=\bar{s}^{\varepsilon}\left(\bar{y}_{j}-\right),\quad c^{L}=\bar{c}^{\varepsilon}\left(\bar{y}_{j}-\right)=c_{j-1},\quad s^{R}=\bar{s}^{\varepsilon}\left(\bar{y}_{j}+\right),\quad c^{R}=\bar{c}^{\varepsilon}\left(\bar{y}_{j}+\right)=c_{j},$$

and the two auxiliary monotone functions (the first one non increasing and the second one non decreasing)

$$G^{L}\left(\sigma\right)=\begin{cases}\max\left\{g\left(\varsigma,c^{L},k_{0}\right)\mid\varsigma\in\left[\sigma,s^{L}\right]\right\},&\text{ for }\sigma\leq s^{L},\\[10.0pt] \min\left\{g\left(\varsigma,c^{L},k_{0}\right)\mid\varsigma\in\left[s^{L},\sigma\right]\right\},&\text{ for }\sigma\geq s^{L},\end{cases}$$

$$G^{R}\left(\sigma\right)=\begin{cases}\min\left\{g\left(\varsigma,c^{R},k_{0}\right)\mid\varsigma\in\left[\sigma,s^{R}\right]\right\},&\text{ for }\sigma\leq s^{R},\\[10.0pt] \max\left\{g\left(\varsigma,c^{R},k_{0}\right)\mid\varsigma\in\left[s^{R},\sigma\right]\right\},&\text{ for }\sigma\geq s^{R}.\end{cases}$$

Call $\gamma$ the unique level at which $G^{L}$ and $G^{R}$ intersect. Because of the hypotheses on $g$, $\gamma$ is equal to either $g\left(s^{L},c^{L},k_{0}\right)$, $g\left(s^{R},c^{R},k_{0}\right)$, a maximum of either $g\left(\cdot,c^{L},k_{0}\right)$ or $g\left(\cdot,c^{R},k_{0}\right)$, or a point in which these two function intersect with derivatives having opposite sign. In any case $\gamma\in\mathcal{G}_{0}$ holds. Define the closed intervals

$$I^{L}=\left[G^{L}\right]^{-1}\left(\left\{\gamma\right\}\right),\qquad I^{R}=\left[G^{R}\right]^{-1}\left(\left\{\gamma\right\}\right).$$

Finally call $s^{-}$ and $s^{+}$ respectively the unique projections of $s^{L}$ and $s^{R}$ on the closed strictly convex sets $I^{L}$ and $I^{R}$. It is not difficult to show that

$$\begin{array}[]{ll}\gamma=G^{L}\left(s^{-}\right)=g\left(s^{-},c^{L},k_{0}\right),&s^{L},s^{-}\in\mathcal{S}_{0,j-1},\\ \gamma=G^{R}\left(s^{+}\right)=g\left(s^{+},c^{R},k_{0}\right),&s^{R},\;s^{+}\in\mathcal{S}_{0,j}.\end{array}$$

Take any wave, with left and right states $s_{l},s_{r}\in\mathcal{S}_{0,j-1}$, of the entropic solution to the Riemann problem

$$\partial_{t}s+\partial_{x}f^{0,j-1}\left(s\right)=0,\qquad s\left(0,x\right)=\begin{cases}s^{L}&\text{ for }x<\bar{y}_{j},\\ s^{-}&\text{ for }x>\bar{y}_{j},\end{cases}$$

(2.8)

and suppose $s^{L}<s^{-}$. Then, $s^{L}\leq s_{l}<s_{r}\leq s^{-}$ and, because of the entropy condition [Bbook, Theorem 4.4], its speed satisfies ($f^{0,j-1}$ coincides with $f\left(\cdot,c^{L},k_{0}\right)$ on $\mathcal{S}_{0,j-1}$):

$$\begin{split}\lambda_{s}&=\frac{f\left(s_{r},c^{L},k_{0}\right)-f\left(s_{l},c^{L},k_{0}\right)}{s_{r}-s_{l}}\leq\frac{f\left(s^{-},c^{L},k_{0}\right)-f\left(s_{l},c^{L},k_{0}\right)}{s^{-}-s_{l}}\\ &\quad=g\left(s^{-},c^{L},k_{0}\right)+s_{l}\frac{g\left(s^{-},c^{L},k_{0}\right)-g\left(s_{l},c^{L},k_{0}\right)}{s^{-}-s_{l}}\\ &\quad=\lambda_{c}+s_{l}\frac{G^{L}\left(s^{-}\right)-g\left(s_{l},c^{L},k_{0}\right)}{s^{-}-s_{l}}\\ &\quad\leq\lambda_{c},\end{split}$$

with $\lambda_{c}=\gamma=g\left(s^{-},c^{L},k_{0}\right)=g\left(s^{+},c^{R},k_{0}\right)$ and where we used the definition of $G^{L}$. If instead $s^{-}<s^{L}$, then $s^{-}\leq s_{r}<s_{l}\leq s^{L}$ and, as in the previous computation, we have

$$\lambda_{s}\leq\lambda_{c}+s_{l}\frac{G^{L}\left(s^{-}\right)-g\left(s_{l},c^{L},k_{0}\right)}{s^{-}-s_{l}}\leq\lambda_{c},$$

Therefore, in any case, the solution to the Riemann problem (2.8) can be patched with a $c$ wave that travel with speed $\lambda_{c}$ and connect the left state $\left(s^{-},c^{L},k_{0}\right)$ to the right state $\left(s^{+},c^{R},k_{0}\right)$. Similar computations can be done at the right of the $c$ wave so that the complete solution includes a $c$ wave travelling with speed $\lambda_{c}$, possibly together some entropic $s$ waves to its left (solutions to $\partial_{t}s+\partial_{x}f^{0,j-1}\left(s\right)=0$) and some entropic $s$ waves to its right (solutions to $\partial_{t}s+\partial_{x}f^{0,j}\left(s\right)=0$). We also use this Riemann solver whenever, for $t>0$, a $c$ wave interact with one or more $s$ waves.

We point out the following properties of this Riemann solver that will be needed in the proof of the main theorem.

• •

  The $c$ wave satisfies Rankine-Hugoniot

  $$\begin{cases}f\left(s^{+},c^{R},k_{0}\right)-f\left(s^{-},c^{L},k_{0}\right)=\lambda_{c}\left(s^{+}-s^{-}\right),\\ c^{R}f\left(s^{+},c^{R},k_{0}\right)-c^{L}f\left(s^{-},c^{L},k_{0}\right)=\lambda_{c}\left(c^{R}s^{+}-c^{L}s^{-}\right).\\ \end{cases}$$

• •

  The $c$ wave is an “admissible” path as defined in [WSCauchy] and satisfies the following entropy condition:

  • –

    if $s^{-}<s^{+}$ there exists $s^{*}\in\left[s^{-},s^{+}\right]$ such that

    $$\begin{cases}g\left(\sigma,c^{L},k_{0}\right)\geq\lambda_{c}&\text{ for all }\sigma\in\left[s^{-},s^{*}\right],\\ g\left(\sigma,c^{R},k_{0}\right)\geq\lambda_{c}&\text{ for all }\sigma\in\left[s^{*},s^{+}\right].\end{cases}$$

    (2.9)

  • –

    if $s^{+}<s^{-}$ there exists $s^{*}\in\left[s^{+},s^{-}\right]$ such that

    $$\begin{cases}g\left(\sigma,c^{R},k_{0}\right)\leq\lambda_{c}&\text{ for all }\sigma\in\left[s^{+},s^{*}\right],\\ g\left(\sigma,c^{L},k_{0}\right)\leq\lambda_{c}&\text{ for all }\sigma\in\left[s^{*},s^{-}\right].\end{cases}$$

    (2.10)

We define the open regions

$$\Omega_{0,j}=\left\{(t,x)\in\left[0,+\infty\right)\times\mathbb{R}\mid x<\bar{x}_{1}\land y_{j}(t)<x<y_{j+1}(t)\right\},$$

and the flux

$$F^{\varepsilon}\left(t,x,\sigma\right)=f^{0,j}\left(\sigma\right),\text{ for }\left(t,x\right)\in\Omega_{0,j}.$$

The wave front tracking approximation $s^{\varepsilon}$ so constructed is an exact weak entropic solution to

$$\begin{cases}\partial_{t}s^{\varepsilon}+\partial_{x}\left[F^{\varepsilon}\left(t,x,s^{\varepsilon}\right)\right]=0,\\ \partial_{t}\left(c^{\varepsilon}s^{\varepsilon}\right)+\partial_{x}\left[c^{\varepsilon}F^{\varepsilon}\left(t,x,s^{\varepsilon}\right)\right]_{x}=0,\\ \end{cases}$$

in the region $x<\bar{x}_{1}$.

Since the $c$ family is linearly degenerate, $c$ waves will not interact with each other. Indeed, given two consecutive $c$ waves located respectively in $y_{j}(t)$ and $y_{j+1}(t)$, the first conservation law in the previous system implies

$$\begin{split}\frac{d}{dt}\int_{y_{j}(t)}^{y_{j+1}(t)}s^{\varepsilon}\left(t,x\right)\;dx&=\dot{y}_{j+1}\left(t\right)s^{\varepsilon}\left(t,y_{j+1}\left(t\right)-\right)-\dot{y}_{j}\left(t\right)s^{\varepsilon}\left(t,y_{j}\left(t\right)+\right)\\ &\quad-f^{0,j}\left(s^{\varepsilon}\left(t,y_{j+1}\left(t\right)-\right)\right)+f^{0,j}\left(s^{\varepsilon}\left(t,y_{j}\left(t\right)+\right)\right)=0\end{split}$$

since $\dot{y}=\lambda_{c}=\frac{f}{\sigma}$. Hence $c$ waves cannot interact with each other (if $s^{\varepsilon}=0$ between two $c$ waves, then they both must travel with zero speed and therefore even in this case they cannot interact).

Since any interaction with the $k$ wave located at $\bar{x}_{1}$ cannot give rise to waves entering the region $x<\bar{x}_{1}$, following [WSCauchy], the wave front tracking algorithm can be carried out for all times in that region. Observe that for a fixed $\varepsilon$ the total variation of the singular variable $P$ introduced in (2.6) is bounded. Since the grid $\mathcal{S}_{0,j}$ contains all the possible maximum points of $g\left(\cdot,c_{j},k_{0}\right)$, in the regions $\Omega_{0,j}$, $P\left(\sigma,c_{j},k_{0}\right)=\int_{0}^{\sigma}\left|\partial_{\xi}\frac{f^{0,j}\left(\xi\right)}{\xi}\right|d\xi$ for any $\sigma\in\mathcal{S}_{0,j}$. The variable $P$ is well behaved in the interplay between the two resonant waves $s$ and $c$. Unfortunately, this behavior is disrupted by the third family of waves, the $k$ waves, except in the case where very strong hypothesis are assumed on the flux as in [CocliteRisebro]. In fact, in [CocliteRisebro] it is assumed that the point of maximum of $g$ does not change with $k$ (actually in [CocliteRisebro] our $k$ waves correspond to discontinuities in time because of Lagrangian coordinates). Since such assumptions are not realistic for our model, we are not able to prove a bound on the total variation of $P$ uniformly in $\varepsilon$. Instead, we resolve this difficulty by applying a compensated compactness argument.

Up to now, all the values of $\left(s^{\varepsilon},c^{\varepsilon},k^{\varepsilon}\right)$ are determined for $x<\bar{x}_{1}$. Since $k^{\varepsilon}$ is constant in time, its value at the right of $\bar{x}_{1}$ is known. The jump conditions $\Delta f=\Delta c=0$ determine all the values of $\left(s^{\varepsilon},c^{\varepsilon},k^{\varepsilon}\right)$ to the right of $\bar{x}_{1}$. These values could introduce new values for the function $g$ that must be added to the grid, i.e.,

$$\mathcal{G}_{1}=\mathcal{G}_{0}\cup\left\{g\left(s^{\varepsilon}\left(t,\bar{x}_{1}+\right),c^{\varepsilon}\left(t,\bar{x}_{1}+\right),k^{\varepsilon}\left(\bar{x}_{1}+\right)\right)\mid\;t\geq 0\right\}.$$

This gives the new allowed values for $s$ for $\bar{x}_{1}<x<\bar{x}_{2}$:

$$\mathcal{S}_{1,j}=\left\{\sigma\mid g\left(\sigma,c_{j},k_{1}\right)\in\mathcal{G}_{1}\right\}.$$

Using these values we now build the corresponding approximations $f^{1,j}\left(\sigma\right)$ of the flux as the linear interpolation of $f\left(\sigma,c_{j},k_{1}\right)$ according to the points in $\mathcal{S}_{1,j}$. Then we solve, as before, all the Riemann problems at $t=0$, $\bar{x}_{1}<x<\bar{x}_{2}$ and $t\geq 0$, $x=\bar{x}_{1}$. As before $c$ waves cannot interact with each other so, by induction, we can carry out the wave front tracking algorithm on the semi plane $t\geq 0$.

We define the open regions ($\bar{x}_{0}=y_{0}(t)=-\infty$, $\bar{x}_{N+1}=y_{M+1}(t)=+\infty$)

$$\Omega_{i,j}=\left\{(t,x)\in\left[0,+\infty\right)\times\mathbb{R}\mid\bar{x}_{i}<x<\bar{x}_{i+1}\land y_{j}(t)<x<y_{j+1}(t)\right\}.$$

The wave front tracking approximations so obtained are weak entropic solutions to

$$\begin{cases}\partial_{t}s^{\varepsilon}+\partial_{x}\left[F^{\varepsilon}\left(t,x,s^{\varepsilon}\right)\right]=0,\\ \partial_{t}\left(c^{\varepsilon}s^{\varepsilon}\right)+\partial_{x}\left[c^{\varepsilon}F^{\varepsilon}\left(t,x,s^{\varepsilon}\right)\right]=0,\\ \partial_{t}k^{\varepsilon}=0,\\ \left(s^{\varepsilon},c^{\varepsilon},k^{\varepsilon}\right)(0,x)=\left(\bar{s}^{\varepsilon},\bar{c}^{\varepsilon},\bar{k}^{\varepsilon}\right)(x),\end{cases}$$

(2.11)

where the flux $F^{\varepsilon}$ is defined by

$$F^{\varepsilon}\left(t,x,\sigma\right)=f^{i,j}\left(\sigma\right),\text{ for }\left(t,x\right)\in\Omega_{i,j}.$$

For all $\left(t,x\right)\in\left[0,+\infty\right[\times\mathbb{R}\text{ and }\sigma\in\left[0,1\right]$, the flux satisfies the estimates

$$\begin{cases}\left|F^{\varepsilon}\left(t,x,\sigma\right)-f\left(\sigma,c^{\varepsilon}\left(t,x\right),k^{\varepsilon}\left(x\right)\right)\right|\leq\varepsilon,\\ \left|\partial_{\sigma}F^{\varepsilon}\left(t,x,\sigma\right)-\partial_{\sigma}f\left(\sigma,c^{\varepsilon}\left(t,x\right),k^{\varepsilon}\left(x\right)\right)\right|\leq\frac{\varepsilon}{N+M}.\end{cases}$$

(2.12)

We remark that, in any region $\Omega_{i,j}$, $s^{\varepsilon}$ is an entropic solution to the scalar conservation law

$$\partial_{t}s^{\varepsilon}+\partial_{x}\left[f^{i,j}\left(s^{\varepsilon}\right)\right]=0.$$