---
abstract: |
    We prove existence and uniqueness of a solution to the Cauchy problem corresponding to the equation $$\begin{cases}
    \partial_t u_{\eps,\delta} +\mathrm{div} {\mathfrak f}_{\eps,\delta}(\mx, u_{\eps,\delta})=\eps \Delta u_{\eps,\delta}+\delta(\eps) {\partial}_t \Delta u_{\eps,\delta}, 
    \ \ {\bf x} \in M, \ \ t\geq 0\\
    u|_{t=0}=u_0(\mx).
    \end{cases}$$ Here, $\mff_{\eps,\delta}$ and $u_0$ are smooth functions while $\eps$ and $\delta=\delta(\eps)$ are fixed constants. Assuming $\mff_{\eps,\delta} \to \mff \in L^p( {\mathbb{R}}^d\times {\mathbb{R}};{\mathbb{R}}^d)$ for some $1<p<\infty$, strongly as $\eps\to 0$, we prove that, under an appropriate relationship between $\eps$ and $\delta(\eps)$ depending on the regularity of the flux $\mff$, the sequence of solutions $(u_{\eps,\delta})$ strongly converges in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d)$ towards a solution to the conservation law $$\partial_t u +\mathrm{div} {\mathfrak f}(\mx, u)=0.$$ The main tools employed in the proof are the Leray-Schauder fixed point theorem for the first part and reduction to the kinetic formulation combined with recent results in the velocity averaging theory for the second.
address:
- 'Melanie Graf, University of Vienna, Faculty of Mathematics'
- 'Michael Kunzinger, University of Vienna, Faculty of Mathematics'
- 'Darko Mitrovic, Faculty of Mathematics, University of Montenegro'
- 'Djordjije Vujadinovic, Faculty of Mathematics, University of Montenegro'
author:
- 'M. Graf'
- 'M. Kunzinger'
- 'D. Mitrovic'
- 'Dj. Vujadinovic'
title: A vanishing dynamic capillarity limit equation with discontinuous flux
---

[ subjclass[35K65, 42B37, 76S99]{} ifundefined[subjclassname@2010]{}[ ]{}[ xpxpsubjclassname@2010 ]{}]{}

Introduction and Notation {#intro}
=========================

Flow in a two-phase porous medium is governed by the Darcy law [@bear]

$$\label{D}
q=-K(S) \left( \nabla p +\rho g e_d \right),$$

where $e_d=(0,\dots,0,1),$ is the direction of gravity. The quantity $S$ is the saturation, $p$ is the pressure, and (the vector) $q$ is the flow velocity of the wetting phase (usually water, while the non-wetting one is oil or a gas).

The Darcy law represents conservation of momentum and, in order to close the system, we also need the conservation of mass $$\label{mc}
{\partial}_t S+ {\mathrm{div}}q =0.$$ In the two-dimensional situation, we have three equations, given by and , while we have four unknowns (two velocity components, saturation and pressure). Therefore, usually one assumes a constitutive relation between the pressure $p$, the capillary pressure $P_c$ (equal to differences of pressures between wetting and non-wetting phases), and the saturation $S$. If it is assumed that the capillary pressure is (almost) constant, one can derive the Buckley-Leverett equation (a scalar conservation law derived in [@BL]): $$\label{BLe}
{\partial}_t S+{\partial}_x f(S)=0,$$ for $f(S)=\frac{S^2}{S^2+A(1-S)^2}$ and an appropriate constant $A$.

If we assume that $P_c$ is “static” (independent of the $t$-derivative of $S$) in the sense that $p_c=p_c(S)$ then we arrive at a parabolic perturbation of the Buckley-Leverett equation [@bear]. However, both of the models (the standard Buckley-Leverett or the one perturbed by a parabolic term) appear to give results inconsistent with certain fairly simple experiments [@DC]. Namely, if we take a thin tube filled with dry sand and dip it in water at a constant rate at one side of the tube and then measure the concentration of the water as a function of time, for certain dipping rates the concentration will not be monotonic. This phenomenon, called overshoot, had been noticed long ago and is a sub-effect of the fingering effect [@H]. Namely, the tips of the fingers that appear when water is penetrating into dry land have larger water concentrations than the body of the finger and, interestingly, the tips of the fingers almost do not change their shape or size. On the other hand, an entropy solution to the Buckley-Leverett equation (in the sense of Kruzhkov [@Kru] or more precisely Oleinik [@O] since we are dealing with a Riemann problem) or the parabolic perturbation of the Buckley-Leverett equation with the Riemann initial data $$\label{riemann}
S(0,x)=\begin{cases}
S_L, & x<0\\
0, & x>0
\end{cases}$$ is monotonic (see e.g. the introduction of [@vDPS]). Therefore, many attempts to explain the gap between the standard theory (provided by the Buckley-Leverett or Richards [@R] equations) have been put forth recently. A purely mathematical approach can be found in [@KMN], where $\delta$-type-solutions to are constructed and $\delta$-distributions appearing as a part of the solution are explained as an inadequacy of the model. In [@PRE], a model is suggested that involves the notion of non-local energy, which eventually leads to a fourth-order equation whose solution (with the Riemann initial data given above) has a shape corresponding to the experimental results from [@DC].

In [@vDPS] one can find an interesting approach explaining the discrepancy between theory and experiment described above, based on the dynamics capillarity concept introduced in [@GH1; @GH2], which has drawn a lot of interest (especially after the publication of [@vDPS]; see also [@AndrCances]). Namely, in [@GH1; @GH2] it was supposed that the capillary pressure depends not only on the saturation $S$, but also on the time derivative of the saturation (of the wetting phase): $$P_c=p_c(S)-\phi \tau \frac{{\partial}S}{{\partial}t}.$$ Taking this into account and proceeding along the lines of deriving the Buckley-Leverett equation, one reaches a nonlinear pseudo-parabolic equation ([@vDPS (1.16)]) which, after linearization of higher order terms reduces to $$\label{BL3}
{\partial}_t S+{\partial}_x f(S)=\eps {\partial}_{xx}S+\tau \delta {\partial}_{xx} {\partial}_t S,$$ where $\eps$ and $\delta$ ($\delta=\eps^2$ in [@vDPS]) are small parameters, while $\tau$ is a fixed constant. By analyzing possible traveling wave solutions $S(\frac{x-ct}{\eps})$, $S(-\infty)=S_L$, $S(+\infty)=S_R$, where the constant $c$ is given by the Rankine-Hugoniot condition $c=\frac{f(S_L)-f(S_R)}{S_L-S_R}$, the authors arrive at solutions constructed from elementary waves (i.e., solutions consisting of shock waves and rarefaction waves) to , that are non-standard (i.e., non-admissible in the sense of Oleinik). Moreover, for some values $S_L$, such a solution exhibits an overshoot-type phenomenon (see [@vDPS Figure 7]).

In this paper, we shall make a step forward in the sense that we shall consider a multi-dimensional generalization to with a flux $\mff=\mff(\mx,\lambda)$ explicitly depending on the space variable and we shall analyze the dynamics capillarity limit to as $\eps\to 0$ for arbitrary initial data.

The $\mx$-dependence of the flux means that we assume that the medium in which we consider the phenomenon is heterogeneous, i.e. that it has different properties at different points (e.g. since in some parts of the medium we have sand while in some others clay). Moreover, we shall assume that the flux is discontinuous with respect to the space variable $\mx$, which means that the medium experiences abrupt (discontinuous) changes in its properties (one can imagine that we have sand which is highly permeable adjacent to the clay layer, which is weakly permeable — permeability is discontinuous in such a medium). Let us remark here that, stipulated by different applications, evolutionary equations with discontinuous flux have attracted considerable attention recently. For a (non-exhaustive) selection of recent results, cf. [@AKR; @AM; @crasta; @kar3; @Pan] and references therein.

Explicitly, we consider the conservation law

$$\label{cl-1}
{\partial}_t u+{\mathrm{div}}\mff(\mx,u)=0,$$

where

-   $\mff = (f_i)_{i=1}^d = (\mx,\lambda)\mapsto \mff(\mx,\lambda)$, $f_i \in L^1({\mathbb{R}}^d\times {\mathbb{R}}) \cap L^p({\mathbb{R}}^d\times {\mathbb{R}})$ and ${\partial}_{\lambda} f_i \in L^p({\mathbb{R}}^d\times {\mathbb{R}})$ for some $1<p<\infty$ and all $i=1,\dots, d$;

-   ${\mathrm{div}}_{\mx} \mff(\mx,\xi) \in {\cal M}({\mathbb{R}}^d\times {\mathbb{R}})$, where ${\cal M}({\mathbb{R}}^d\times {\mathbb{R}})$ is the space of Radon measures, and there exists a constant $\beta>0$ and a finite Radon measure $\mu \in {\cal M}({\mathbb{R}}^d)$ such that (in the sense of measures) $$|{\mathrm{div}}_{\mx} \mff(\mx,\xi)|\leq \frac{\mu(\mx)}{1+|\lambda|^{1+\beta}};$$

-   $\|\sup\limits_{\lambda \in {\mathbb{R}}} |{\partial}_{\lambda}f_i(\cdot,\lambda)| \|_{L^1({\mathbb{R}}^d)} <\infty$ for all $i=1,\dots, d$;

We then perturb by (a linearized) vanishing dynamic capillarity limit:

$$\label{vc-1}
{\partial}_t u_{\eps}+{\mathrm{div}}\mff_\eps (\mx,u_\eps)=\eps  \Delta u_\eps +\delta (\eps)\Delta {\partial}_t u_\eps.$$

Here, $\mff_\eps := K_\eps \cdot  \mff \star \omega_{n(\eps )}$, where componentwise convolution $\mff \star \omega_{n(\eps )}$ is a regularization of $\mff$ and $K_\eps=K_\eps(\mx,\lambda)$ is a bounded family of compactly supported functions equal to one on the ball $B(0,1/\eps)\subset {\mathbb{R}}^{d+1}$ and such that $\|\nabla K_\eps\|_{L^\infty({\mathbb{R}}^{d+1})} \leq 1$.

More precisely, we suppose that $\omega\in C^\infty_c({\mathbb{R}}^{d+1})$ is of the form $\omega(\mx,\lambda)=\omega^{(1)}(\mx)$ $\omega^{(2)}(\lambda)$ with $\omega^{(1)}, \omega^{(2)}$ test functions with unit integral on ${\mathbb{R}}^d$ and ${\mathbb{R}}$, respectively, and $$\label{mollifying}
\omega_{n(\eps)}(\mx,\lambda) \equiv \omega^{(1)}_{n(\eps)}(\mx)\omega^{(2)}_{n(\eps)}(\lambda) =
n(\eps)^{-d} \omega^{(1)}(\mx/n(\eps))n(\eps)^{-1} \omega^{(2)}(\lambda/n(\eps))$$ with an appropriate $n(\eps)$ such that $n(\eps )\to 0$ as $\eps\to 0$ (the form of $n(\eps)$ will be precisely determined later). Note that $$K_\eps \mff_\eps  \to \mff  \text{ strongly in } L^p({\mathbb{R}}^d\times {\mathbb{R}}) \text{ as } \eps\to 0.$$ We supplement with the initial data $$\label{id}
u|_{t=0}=u_0(\mx)\in L^1({\mathbb{R}}^d) \cap L^\infty({\mathbb{R}}^d),$$ while we take $$\label{id-pert}
u_\eps\big|_{t=0}=u_{0}^\eps(\mx) \to u_0(\mx) \ \ \text{strongly in} \ \ L^p({\mathbb{R}}^d),$$ where $$\label{cond}
\|u_{0}^{\eps}\|_{L^{2}({\mathbb{R}}^{d})}+n(\eps) \|\nabla u_{0}^{\eps}\|_{L^{2}({\mathbb{R}}^{d})}+n(\eps)^{2}\|\nabla u_{0x_{i}}\|_{L^{2}({\mathbb{R}}^{d})}\leq C_{0}.$$ It is not difficult to see that $u_{0}^\eps=u_0\star \omega^{(1)}_{n(\eps)}$, where $\omega^{(1)}_{n(\eps)}$ is as above, satisfies . We shall, however, use as an assumption.

With regards to pseudo-parabolic equations, we have already explained their importance in the porous media theory. Besides a possible justification of the experimental results from [@DC], there is a confirmed application to the seepage of homogeneous fluids through a fissured rock [@2Y]. Also, they describe the unidirectional propagation of nonlinear, dispersive long waves [@3Y; @29Y] (where $u$ is typically the amplitude or velocity), or population dynamics [@22Y] (where $u$ represents the population density).

In [@26; @30Y], the authors investigated the initial-boundary value problem and the Cauchy problem for linear pseudo-parabolic equations and established the existence and uniqueness of solutions. As for the non-linear variants, one can find numerous results even for singular pseudo-parabolic equations and degenerate pseudo-parabolic equations (see e.g. [@CYW; @KNS; @Mi; @12Y; @17Y; @Pt] and the references therein). Together with the existence and uniqueness results, among the given references, one can also find properties of solutions, such as asymptotic behavior and regularity. However, we are not aware of any corresponding results for the Cauchy problem for an equation of type .

This problem will be considered in Section \[sec-2\]. We shall use an approach that is characteristic for the theory of wave equations [@Ach]: we shall apply the Fourier transform with respect to $\mx$, solve the ordinary differential equation so obtained with respect to $t$ and prove the following theorem (the constants $\eps$ and $\delta$ are omitted for simplicity):

\[unique\_sol\_par\] There exists a unique solution to $$\label{main-eq}
u_{t}+{\mathrm{div}}(\mff(\mx, u)) =\triangle u+\partial_{t}\triangle u, \ \ (t,x)\in [0,T)\times {\mathbb{R}}^d,$$ where $T>0$, $\mff \in C^\infty_c({\mathbb{R}}^d\times  {\mathbb{R}})$ is a bounded function, supplemented with (merely) the initial condition $$\label{id-nonpert}
u|_{t=0}=u_0(\mx) \in C_c^\infty({\mathbb{R}}^d),$$ which belongs to $L^2((0,T)\times {\mathbb{R}}^d)\cap C^\infty((0,T)\times {\mathbb{R}}^d)$.

In the next step (Section 3), we are going to let the perturbation parameter $\eps\to 0$ in . This type of problem for conservation laws (when there is not only vanishing viscosity, but also third or higher order perturbation) was first addressed in [@Sch] and received considerable attention after that. A thorough analysis of this kind of limit can be found in [@Lfl-K], where the theory of non-classical shocks for conservation laws was essentially initiated. The standard approach here is to rewrite the equation under consideration in the kinetic formulation, or using Young measures (which is essentially equivalent) and then applying velocity averaging results [@AMP; @HTz], compensated compactness (in one dimensional situations) [@coclitte1; @coclitte; @Sch], or Di Perna techniques involving Young measures [@Dpe], as done in [@Lfl-b] and many others. Again we note that this list of citations is far from complete. The problem in our case is the low regularity of the flux function $\mff$ (it can be discontinuous with respect to $\mx\in {\mathbb{R}}^d$) as well as the multidimensional character of our problem, which prevents us from using the Young measures and compensated compactness approach here (see item (iii) below and Remark \[39\] for the case of a regular flux). We also remark that, to the best of our knowledge, the only diffusion-dispersion type result for equations with discontinuous coefficients is [@HKM].

Moreover, most of the velocity averaging results are given in the case of a homogeneous flux [@DLM; @GSR; @PS; @TT] or a flux for which $p \geq 2$ [@Ger; @LM2]. We have recently proved [@LM5] the velocity averaging lemma in the case $p>1$ and this will enable us to prove the strong convergence of the sequence $(u_\eps)$ along a subsequence towards a weak solution to . The following statements are the main results of the paper:

-   If $\delta=o(\eps^2)$ and $$\label{neps}
    \frac{\sqrt{\delta}}{\eps n(\eps)^{d+2}}
    \to 0$$ as $\eps\to 0$, then the family $(u_\eps)$ is strongly precompact in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d)$ under a suitable non-degeneracy condition (see Definition \[def-non-deg\]);

-   If $\delta={\cal O}(\eps^2)$ and $\frac{\sqrt{\delta}}{\sqrt{\eps} n(\eps)^{d/2+2}}\to 0$ as $\eps\to 0$, and the $\lambda$-derivative of the flux ${\partial}_\lambda \mff$ is bounded in addition to condition (C3), then the family $(u_\eps)$ is strongly precompact in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d)$ under the non-degeneracy condition;

-   If $\delta=o(\eps^2)$ and $\frac{\sqrt{\delta}}{\sqrt{\eps} n(\eps)^{d/2+2}}\to 0$ as $\eps\to 0$, and if the flux $\mff \in C^1({\mathbb{R}}^d\times {\mathbb{R}})$, then the family $(u_\eps)$ converges strongly in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d)$ towards the entropy solution to , (no non-degeneracy condition is needed).

Existence and uniqueness of the solution to the pseudo-parabolic equation with  {#sec-2}
===============================================================================

Throughout this section we suppose that $\mff$ is smooth and compactly supported. Under this assumption, we want to show that with the initial condition has a unique solution. The strategy for solving this problem is to define a mapping ${\cal T}: L^2([0,T)\times {\mathbb{R}}^d) \to L^2([0,T)\times {\mathbb{R}}^d)$ such that for every $v\in L^2([0,T)\times {\mathbb{R}}^d)$ the function $$\label{T}
u={\cal T}(v)$$ represents a solution to $$\label{main-eq-1}
u_{t}+{\mathrm{div}}(\mff(\mx, v)) =\triangle u+\partial_{t}\triangle u, \ \ (t,x)\in [0,T)\times {\mathbb{R}}^d,$$ with the initial conditions . Generally, for $u_0\in L^2({\mathbb{R}}^d)$, we are seeking weak solutions $u$ of this initial value problem in the sense that, for any smooth test function $\varphi$ with compact support in $[0,T)$, we require $$\int_0^T\int_{{\mathbb{R}}^d} u\varphi_t + f(x,v){\mathrm{div}}\varphi - u\Delta\varphi + u\partial_t  \Delta\varphi\,dxdt 
+ \int_{{\mathbb{R}}^d} u_0(\varphi(0,x)+\Delta\varphi(0,x))\,dx=0.$$

Then, we shall prove that the mapping ${\cal T}$ possesses a fixed point, which will turn out to be the solution to , .

To this end, we need the following consequence of the Leray-Schauder fixed point theorem (cf. [@GT Th. 11.3]):

\[fixed\] Let ${\cal T}$ be a compact mapping of a Banach space $\mathcal{B}$ into itself and suppose that there exists a constant $C$ such that $$\label{cond-fp}
\|u\|_{{\cal B}}\leq C$$ for all $u\in {\cal B}$ and $\sigma \in [0,1]$ satisfying $u=\sigma {\cal T}u$. Then ${\cal T}$ has a fixed point, that is, ${\cal T}u=u$ for some $u\in \mathcal{B}$.

Finally, let us fix our conventions for Fourier transform and inverse Fourier transform. For $u\in L^1({\mathbb{R}}^d)\cap L^2({\mathbb{R}}^d)$, the Fourier transform of $u$ is $${\cal F}(u)(\mxi)=\hat{u}(\mxi)=\int_{{\mathbb{R}}^d} e^{- i \mx \cdot \mxi} u(\mx)\, d\mx$$ and the inverse Fourier transform is given by $${\cal F}^{-1}(u)(\mx)=\check{u}(x)= \frac{1}{(2\pi)^d}\int_{{\mathbb{R}}^d} e^{i \mx \cdot \mxi} u(\mxi)\, d\mxi,$$ where $\cdot$ denotes the scalar product and $i$ is the imaginary unit.

We begin by proving existence and uniqueness of solutions to , .

\[ex-v\] Let $v\in L^2([0,T]\times {\mathbb{R}}^d)$. Then for any $u_0\in L^2({\mathbb{R}}^d)$ there exists a unique solution $u\in L^2([0,T]\times {\mathbb{R}}^d)$ to with $u|_{t=0}=u_0$.

After applying the Fourier transform with respect to $\mx \in {\mathbb{R}}^d$ to , we obtain $${\partial}_t \hat{u}(t,\mxi)+i \mxi \cdot \widehat{\mff(\cdot,v(t,\cdot))}(\mxi)=-|\mxi|^2 \hat{u}(t,\mxi) - |\mxi|^2  {\partial}_t \hat{u}(t,\mxi),$$ where $|\mxi|=(\xi_1^2+\dots \xi_d^2)^{\frac{1}{2}}$. We rewrite the last equation in the form $$\frac{d \hat{u}(t,\mxi)}{dt} +\frac{|\mxi|^2}{{1+|\mxi|^2}}\hat{u}(t,\mxi) =-\frac{i \mxi \cdot \widehat{\mff(\cdot,v(t,\cdot))}(\mxi)}{1+|\mxi|^2}$$ It follows from this that any weak solution to has $t$-derivative in $L^2$ as well, hence in particular is continuous with respect to $t$ and therefore possesses a classical trace on $t=0$. Solving the ODE with the initial data $\hat{u}(0,\mxi)=\hat{u}_0(\mxi)$ we arrive at $$\begin{aligned}
\label{sol-v}
\hat{u}(t,\mxi)=e^{-\frac{|\mxi|^2 t}{1+|\mxi|^2}} \left( \hat{u}_0(\mxi)-\int_0^t\frac{i \mxi \cdot \hat{\mff}}{1+|\mxi|^2} e^{\frac{|\mxi|^2 t'}{1+|\mxi|^2}}\,dt'\right),\end{aligned}$$ where $\hat{\mff}=\widehat{\mff(\cdot,v(t,\cdot))}(\mxi)$. By finding the inverse Fourier transform here with respect to $\mxi \in {\mathbb{R}}^d$, we indeed obtain a weak solution, so combined with the above considerations both existence and uniqueness follow.

\[estimates-1\] For any function $v\in L^2([0,T]\times {\mathbb{R}}^d)$, the solution to , belongs to $H^1((0,T)\times {\mathbb{R}}^d)$.

Let $u$ be the solution to , defined in Lemma \[ex-v\]. We need to estimate ${\partial}_{x_j} u$, $j=1,\dots,d$ and ${\partial}_t u$, where ${\partial}_{x_j}$, $j=1,\dots, d$ and ${\partial}_t$ are the weak derivatives of $u$. We have according to the Plancherel theorem (below, we take the Fourier transform with respect to $\mx\in {\mathbb{R}}^d$) and the Cauchy-Schwartz inequality $$\label{est-x}
\begin{split}
\|{\partial}_{x_j} & u\|_{L^2([0,T]\times {\mathbb{R}}^d)} = \| i \xi_j \hat{u}\|_{L^2([0,T]\times {\mathbb{R}}^d)} \\
&\leq \| \xi_j \hat{u}_0 \|_{L^2([0,T]\times {\mathbb{R}}^d)}+\left\|\int_0^t\frac{i \xi_j \mxi\cdot \hat{\mff}}{1+|\mxi|^2}e^{\frac{|\mxi|^2 t'}{1+|\mxi|^2}}dt'\right\|_{L^2([0,T]\times {\mathbb{R}}^d)} 
\\& \leq \| \xi_j \hat{u}_0 \|_{L^2([0,T]\times {\mathbb{R}}^d)}+\left( \int_0^T \int_{{\mathbb{R}}^d} \left| \int_0^t \frac{i \xi_j \mxi\cdot \hat{\mff}}{1+|\mxi|^2}e^{\frac{|\mxi|^2 t'}{1+|\mxi|^2}}dt' \right|^2 d\mxi dt\right)^{1/2}
\\& \leq \| \xi_j \hat{u}_0 \|_{L^2([0,T]\times {\mathbb{R}}^d)}+\left( \int_0^T \int_{{\mathbb{R}}^d} t \int_0^t \frac{\xi^2_j |\mxi\cdot \hat{\mff}|^2}{(1+|\mxi|^2)^2}e^{\frac{2|\mxi|^2 t'}{1+|\mxi|^2}}dt' d\mxi dt\right)^{1/2}\\
&\leq \| \xi_j \hat{u}_0 \|_{L^2([0,T]\times {\mathbb{R}}^d)}+T e^{T} 
\left(\int_0^T\int_{{\mathbb{R}}^d} \frac{\xi_j^2 \, |\mxi|^2 \, |\hat{\mff}|^2}{(1+|\mxi|^2)^2} d\mxi dt\right)^{1/2}\\
&\leq  \| \xi_j \hat{u}_0 \|_{L^2([0,T]\times {\mathbb{R}}^d)}+T e^{T} \left(\int_0^T\int_{{\mathbb{R}}^d} |\mff(\mx,v(t,\mx))|^2 d\mx dt\right)^{1/2},
\end{split}$$where in the last step we used the Plancherel theorem again. From here, since $f$ is compactly supported with respect to $\mx \in {\mathbb{R}}^d$ and bounded, we conclude that ${\partial}_{x_j} u \in L^2([0,T]\times {\mathbb{R}}^d)$.

As for the $t$-derivative, we have $$\begin{aligned}
\label{est-t}
{\partial}_t \hat{u}=
-\frac{|\mxi |^2}{1+|\mxi|^2} e^{\frac{-|\mxi |^2 t}{1+|\mxi|^2}} \left(\hat{u}_0(\mxi)-\int_0^t\frac{i\mxi \cdot \hat{\mff}}{1+|\mxi|^2} dt' \right)+\frac{i\mxi \cdot \hat{\mff}}{1+|\mxi|^2}, \end{aligned}$$ and therefore, repeating the procedure giving us the estimates on ${\partial}_{x_j} u$, we get $$\|{\partial}_t u \|_{L^2([0,T]\times {\mathbb{R}}^d)}=\|{\partial}_t \hat{u} \|_{L^2([0,T]\times {\mathbb{R}}^d)} \leq C$$ for a constant $C\in {\mathbb{R}}$.

Based on this we can proceed to applying Theorem \[fixed\] to prove existence of a solution to , .

\[ex-un\] There exists a unique solution to , . This solution belongs to $L^2([0,T]\times {\mathbb{R}}^d)\cap C^\infty([0,T]\times {\mathbb{R}}^d)$.

Let us first take a sequence of balls $B_n=B(0,n)$ centered at $0$ and of radius $n$. Set $u_0^n(\mx)=u_0(\mx) \chi_{B_n}(\mx)$, where $\chi_{B_n}: {\mathbb{R}}^d\to [0,1]$ is a smooth regularization of the characteristic function of the set $B_n$ such that $\chi_{B_n}(\mx)=1$ for $\mx \in B(0,n-1)$ and $\chi_{B_n}(\mx)=0$ for $\mx \notin B(0,n)$. Then we define the mapping ${\cal T}_n: L^2([0,T]\times B_n) \to  L^2([0,T]\times B_n)$ by $${\cal T}_n (v) :=\chi_{B_n} {\cal T} (v\chi_{B_n}),$$ where ${\cal T}$ is the operator from with the initial data $u_0 \chi_{B_n}$ and we extend $v\chi_{B_n}$ by zero so it is defined on $[0,T]\times {\mathbb{R}}^d$. This map is continuous since ${\cal T}$ is continuous. In fact, by the proof of Theorem \[estimates-1\], $T$ is in fact continuous as a map from $L^2([0,T)\times {\mathbb{R}}^d)$ to $ H^1((0,T)\times {\mathbb{R}}^d)$. So we can use Theorem \[estimates-1\] and the Rellich theorem (keeping in mind that $B_n$ is a bounded set) to conclude that ${\cal T}_n$ is a compact mapping.

Let us check condition from Theorem \[fixed\]. If a function $u\in L^2([0,T]\times B_n)$ satisfies $u=\sigma {\cal T}_n u$, then the function $\tilde{u}:=\sigma {\cal T} (u \chi_{B_n})$ satisfies $\chi_{B_n} \tilde{u}(t,\mx)=u(t,\mx)$ for $\mx\in B_{n}$ and is a solution of the Cauchy problem $$\begin{split}
\tilde{u}_{t}+\sigma {\mathrm{div}}(\mff(\mx, u \chi_{B_n})) &=\triangle \tilde{u}+\partial_{t}\triangle \tilde{u}, \ \ (t,x)\in [0,T)\times {\mathbb{R}}^d\\
\tilde{u}(0,\mx) &= \sigma u_0(\mx) \chi_{B_n}(\mx)
\end{split}$$ If we apply Fourier transform with respect to $\mx\in {\mathbb{R}}^d$ here and solve the ODE so obtained, we arrive at an expression analogous to :

$$\begin{aligned}
\hat{\tilde{u}}(t,\mxi)=\sigma e^{-\frac{t|\mxi |^2}{1+|\mxi|^2}} \left( \widehat{u_0\chi_{B_n}}(\mxi)-
\int_0^t\frac{i \mxi \cdot \hat{\mff}}{1+|\mxi|^2} e^{\frac{t' |\mxi |^2}{1+|\mxi|^2}}\,dt'\right),\end{aligned}$$

where $\hat{\mff}={\cal F}(\mff(\mx,u(t,\mx)\chi_{B_n}(t,\mx))$ is the Fourier transform with respect to $\mx$. Now $ \| u \|_{L^2([0,T]\times B_n)} \leq   \| \tilde{u} \|_{L^2([0,T]\times {\mathbb{R}}^d)}$ and repeating the procedure from and , we conclude that holds for the mapping ${\cal T}_n$ (keeping in mind that $n$ is fixed).

Thus, for every $n\in {\mathbb N}$ there exists a function $u_n\in L^2([0,T]\times B_n)$ solving $$\begin{split}
{\partial}_t{u_n}+ {\mathrm{div}}(\mff(\mx, u_n)) &=\triangle u_n+\partial_{t}\triangle u_n, \ \ (t,x)\in [0,T)\times B_{n-1}\\
u_n(0,\mx)&=u_0(\mx), \ \ \mx\in B_{n-1}
\end{split}$$ Extending $u_n$ by zero outside of $B_n$, we thereby obtain a sequence $(u_n)$ that is bounded in $L^2([0,T]\times {\mathbb{R}}^d)$, and locally bounded in $H^1((0,T)\times {\mathbb{R}}^d)$. Thus, $(u_n)$ admits a subsequence that locally converges towards $u\in L^2([0,T]\times {\mathbb{R}}^d)$. It is clear that the function $u$ is the weak solution to , and thus it must belong at least to $H^1((0,T)\times {\mathbb{R}}^d)$.

Moreover, the proof of Theorem \[estimates-1\] with $v=u$ is a classical example of a bootstrapping procedure, which enables us to conclude that $u$ actually belongs to each $H^k((0,T)\times {\mathbb{R}}^d)$ ($k\in {\mathbb N}$), hence is smooth (due to smoothness of $u_0$ and $\mff$).

Now we turn to the proof of uniqueness. Assume that $u_1$ and $u_2$ are solutions to , supplemented with the initial conditions $u_1\big{|}_{t=0}=u_{10}$ and $u_2\big{|}_{t=0}=u_{20}$. After applying the procedure from Theorem \[estimates-1\] we conclude that

$$\begin{aligned}
\label{u-1}
\hat{u}_1(t,\mxi)-\hat{u}_2(t,\mxi)=e^{\frac{-t|\mxi |^2}{1+|\mxi|^2}} \left( \hat{u}_{10}(\mxi)-\hat{u}_{20}(\mxi) -\int_0^t \frac{i\mxi \cdot (\hat{\mff}_1-\hat{\mff}_2)}{1+|\mxi|^2}e^{\frac{t' |\mxi |^2}{1+|\mxi|^2}}dt'\right),\end{aligned}$$

where $$\hat{\mff}_1-\hat{\mff}_2={\cal F}(\mff(\mx,u_1(t,\mx))-\mff(\mx,u_2(t,\mx))).$$ Let us first estimate $$\begin{aligned}
\label{u-3}
\left\|\frac{\mxi\cdot (\hat{\mff}_1-\hat{\mff}_2)}{1+|\mxi|^2}\right\|_{L^2({\mathbb{R}}^d)} &\leq  \| \hat{\mff}_1-\hat{\mff}_2\|_{L^2({\mathbb{R}}^d)} 
=  \| \mff_1-\mff_2\|_{L^2({\mathbb{R}}^d)}  \\
&\leq  \|{\partial}_\lambda \mff\|_{\infty}\| \hat{u}_1-\hat{u}_2\|_{L^2({\mathbb{R}}^d)}. \nonumber\end{aligned}$$ Taking the $L^2({\mathbb{R}}^d)$-norm of we get, due to : $$\begin{aligned}
\label{u-4}
&\|\hat{u}_1(t,\cdot)-\hat{u}_2(t,\cdot)\|_{L^2({\mathbb{R}}^d)} \leq \\ & \|\hat{u}_{10}-\hat{u}_{20}\|_{L^2({\mathbb{R}}^d)} + \|{\partial}_\lambda \mff \|_{\infty}  \int_0^t\|\hat{u}_1(t',\cdot)-\hat{u}_2(t',\cdot)\|_{L^2({\mathbb{R}}^d)} e^{t'} dt'.
\nonumber\end{aligned}$$ Applying Gronwall’s inequality, we arrive at $$\|\hat{u}_1(t,\cdot)-\hat{u}_2(t,\cdot)\|_{L^2({\mathbb{R}}^d)} \leq \|\hat{u}_{10}-\hat{u}_{20}\|_{L^2({\mathbb{R}}^d)} 
e^{(e^t-1)\|{\partial}_\lambda \mff\|_{\infty}}.$$ From here, we see that $u_{10}=u_{20}$ entails $u_1=u_2$. This concludes the proof.

Vanishing capillarity limit
===========================

In this section, we inspect the vanishing capillarity limit of . Unlike the situation that is typical in the case of the diffusion-dispersion limit ([@Sch; @AMP; @Lfl-b] etc.), where the existence of a solution to the perturbed equation is assumed together with all necessary properties of the solution, we have proved in the previous section the existence of the functional sequence whose convergence we analyze.

By Theorem \[unique\_sol\_par\], under the assumptions (C1)–(C3) from Section \[intro\], the equation with initial data (satisfying ) possesses a unique solution $u_\eps$ in $L^2((0,T)\times {\mathbb{R}}^d)\cap C^\infty((0,T)\times {\mathbb{R}}^d)$. As announced in Section \[intro\], our main result then is as follows:

Under the above assumptions, if $(F_1,\dots, F_d)=\partial_\lambda \mff$ lies in $L^{\bar{p}'}({\mathbb{R}}^{d+1})$ for some $\bar{p}'\in \big(1,2+\frac{4}{d}\big]$ then: \[main-theorem\]

-   If $\delta=o(\eps^2)$ and $n(\eps)$ satisfies then the family $(u_\eps)$ contains a strongly convergent subsequence in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d)$ under the non-degeneracy condition given in Definition \[def-non-deg\] below.

-   If $\delta={\cal O}(\eps^2)$ and $\frac{\sqrt{\delta}}{\sqrt{\eps} n(\eps)^{d/2+2}}\to 0$ as $\eps\to 0$, and the $\lambda$-derivative of the flux ${\partial}_\lambda \mff$ is bounded in addition to condition (C3), then the family $(u_\eps)$ contains a strongly convergent subsequence in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d)$ under the non-degeneracy condition;

-   If $\delta=o(\eps^2)$ and $\frac{\sqrt{\delta}}{\sqrt{\eps} n(\eps)^{d/2+2}}\to 0$ as $\eps\to 0$, and if the flux $\mff\in C^1({\mathbb{R}}^d\times {\mathbb{R}})$, then the family $(u_\eps)$ converges strongly in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d)$ towards the entropy solution to , .

Below we give a complete proof of (i). Although non-trivial, the other two claims are more standard and we shall only comment on the necessary modifications of the argument for establishing (i).

To proceed, let us first introduce the non-degeneracy condition. Such a condition is standard in many results using velocity averaging lemmas (cf., e.g., [@LM5; @LPT; @Pan; @PS; @TT]).

\[def-non-deg\] The flux appearing in is called non-degenerate if for the function ${\bf F}=(F_1,\dots, F_d)=\partial_\lambda \mff$, where $\mff=\mff(\mx,\lambda)$ we have: for almost every $(t,\mx)\in {\mathbb{R}}_+^{d}$ and every $\mxi\in S^{d}$, where $S^d$ is the unit sphere in ${\mathbb{R}}^{d+1}$ the mapping

$$\label{non-deg} \lambda \mapsto  \left( \xi_0+
\sum\limits_{k=1}^d F_k(\mx,\lambda) \xi_k \right),$$

is not zero on any set of positive measure.

Under these assumptions, we proved the following assertion in [@LM5 Theorem 3.4].

\[va\] Assume that $h_n\rightharpoonup h$ weakly in ${L}_{loc}^{s}({\mathbb{R}}^{d+1})$ for some $s\geq 2$, where $h_n$ are weak solutions to

$$\label{lin-eq}
\begin{split}
\sum\limits_{k=1}^d{\partial}_{x_k}
\left(F_k(\mx,\lambda) h_n(\mx,\lambda)\right)={\partial}_\lambda^\kappa
G_n(\mx,\lambda)\,.
\end{split}$$

Here, $\kappa\in {\mathbb N}_0$ and $\lambda \in {\mathbb{R}}$. Let $\bar{p}\in (1,s)$ and let $\frac{1}{\bar{p}}+\frac{1}{\bar{p}'}=1$. We assume

-   $F_k\in { L}^{\bar p'}({\mathbb{R}}^{d+1})$ for $k=1,\dots,d, \ $

-   The sequence $(G_n)$ is strongly precompact in the space ${ W}_{loc}^{-1,r}({\mathbb{R}}^{d+1})$, where $r>1$ satisfies the relation $1+\frac{1}{s}=\frac{1}{\bar{p}}+\frac{1}{r}$.

Finally, assume that the non-degeneracy condition is satisfied.

Then, for any $\rho\in L^2_c({\mathbb{R}})$, there exists a subsequence $(h_{n_k})$ of $(h_n)$ such that $$\label{va_result}
\int_{{\mathbb{R}}}\rho(\lambda) h_{n_k}(\mx,\lambda)d\lambda \to  \int_{{\mathbb{R}}}\rho(\lambda) h(\mx,\lambda)d\lambda \ \ \text{ strongly in
$L^1_{loc}({\mathbb{R}}^d)$}$$ as $k\to\infty$.

Now, we proceed to the estimates sufficient to apply the mentioned functional analytic tools.

\[osnov\] The solution $u_{\eps}$ of satisfies the following a priori estimates: $$\begin{aligned}
\label{prva}
\|u_{\eps}(t,\,.\,)\|_{L^{2}({\mathbb{R}}^{d})}\leq \|u_{0}^{\eps}\|_{L^{2}({\mathbb{R}}^{d})}+\sqrt{\delta}\|\nabla u^\eps_{0}\|_{L^{2}({\mathbb{R}}^{d})}+\sqrt{2 t} \| {\mathrm{div}}\mff_\eps(\mx,\lambda) \|^{\frac{1}{2}}_{L^{1}({\mathbb{R}}^{d}\times {\mathbb{R}})},\end{aligned}$$ $$\label{druga}
\begin{split}
\sqrt{2\eps}\|\nabla u_{\eps}\|_{L^{2}([0,T],L^{2}({\mathbb{R}}^{d}))} \leq \|u_{0}^{\eps}\|_{L^{2}({\mathbb{R}}^{d})} +\sqrt{\delta} & 
\|\nabla u_{0}^{\eps}\|_{L^{2}({\mathbb{R}}^{d})}
\\ &+\sqrt{2 T} \| {\mathrm{div}}\mff_\eps(\mx,\lambda) \|^{\frac{1}{2}}_{L^{1}({\mathbb{R}}^{d}\times {\mathbb{R}})},
\end{split}$$ $$\label{treca}
\begin{split}
\sqrt{\delta} \|\nabla u_{\eps}(t,\,.\,) \|_{L^{2}({\mathbb{R}}^{d})}\leq \|u_{0}^{\eps}\|_{L^{2}({\mathbb{R}}^{d})}+\sqrt{\delta}  \|\nabla & u_{0}^{\eps}\|_{L^{2}({\mathbb{R}}^{d})}\\
&+\sqrt{2 t} \| {\mathrm{div}}\mff_\eps(\mx,\lambda) \|^{\frac{1}{2}}_{L^{1}({\mathbb{R}}^{d}\times {\mathbb{R}})}
\end{split}$$ $$\label{cetvrta}
\begin{split}
\eps \delta \|\nabla {\partial}_{x_i} u_{\eps} & \|^2_{L^{2}([0,T],L^{2}({\mathbb{R}}^{d}))} \leq 
2\delta \|{\partial}_{x_i}u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}+2\delta^2 \|\nabla {\partial}_{x_i} u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}  \\  
&+\frac{d\delta}{\eps^2} \|{\partial}_\lambda \mff_\eps\|^2_{\infty} ( \|u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}+\delta\|\nabla u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})} \\
& +2 T \| {\mathrm{div}}\mff_\eps(\mx,\lambda) \|_{L^{1}({\mathbb{R}}^{d}\times {\mathbb{R}})}) +\frac{4  T \delta }{\eps} \| {\mathrm{div}}\mff^\eps(\mx, \lambda)\|^2_{L^2({\mathbb{R}}^d, L^\infty ({\mathbb{R}}))} 
\end{split}$$ $$\label{peta}
\begin{split}
 \delta \|\nabla & {\partial}_t u_\eps\|_{L^2([0,T],L^2({\mathbb{R}}^d))} 
 \leq  \frac {\sqrt{\delta}  \sqrt{d}}{2\sqrt{\eps}} \|{\partial}_\lambda \mff_\eps\|_{\infty}  \|u^\eps_{0}\|_{L^{2}({\mathbb{R}}^{d})}\\  
 &\hspace*{-2em}+
 \sqrt{\delta}  \left(\frac{\sqrt{\delta}\sqrt{d}}{2\sqrt{\eps}} \|{\partial}_\lambda \mff_\eps\|_{\infty}
 +\frac{\sqrt{\eps}}{\sqrt{2}}\right)\|\nabla u^\eps_{0}\|_{L^{2}({\mathbb{R}}^{d})} \\
 &\hspace*{-2em} +\frac{\sqrt{\delta}  \sqrt{dT}}{\sqrt{2\eps}}  \|{\partial}_\lambda \mff_\eps\|_{\infty}
 \|{\mathrm{div}}\mff_\eps \|_{L^1({\mathbb{R}}^d\times {\mathbb{R}})}^{\frac{1}{2}}+\sqrt{\delta}   \sqrt{T} 
 \|{\mathrm{div}}\mff_\eps \|_{L^2({\mathbb{R}}^d,L^\infty({\mathbb{R}}))}.
\end{split}$$

Let us consider a family $(u_{\eps})$ of smooth solutions to . Multiplying the equation by $\eta'(u_\eps)$ (with $\eta(u_\eps)=u_\eps^{2})$, we have $$\begin{aligned}
\label{-2}
&{\partial}_t(u_\eps^{2})+2\sum_{j=1}^{d}q^\eps_{j}(\mx,u_{\eps})_{x_j}+2 \int_0^{u_{\eps}} {\mathrm{div}}\mff_\eps(\mx,\lambda) d\lambda 
\\&=2\eps\sum_{j=1}^{n}(u_\eps {\partial}_{x_j}u_{\eps})_{x_j}-2\eps|\nabla u_\eps|^{2}+ \delta(\triangle\eta(u_\eps))_{t}-2\delta \sum_{j=1}^{n}\eta({\partial}_{x_j}u_{\eps})_t-2\delta\triangle u_{\eps} {\partial}_t u_{\eps},
\nonumber\end{aligned}$$ where $\bf q_\eps$ is defined by (recall that $\mff_\eps=(f_1^\eps,\dots,f_d^\eps)$) $${\partial}_\lambda{\bf q}_\eps(\mx,\lambda)=({\partial}_\lambda F^\eps_1(\mx,\lambda), \dots, {\partial}_\lambda F^\eps_d(\mx,\lambda))
=(\lambda {\partial}_\lambda f^\eps_1(\mx,\lambda),\dots, \lambda {\partial}_\lambda f^\eps_d(\mx,\lambda)),$$ and normalized by the condition $q^\eps_{j}(\mx,0)=0$, $j=1,...,d$. Integrating over $\mx\in {\mathbb{R}}^d$ and $0\leq t\leq T$, and using integration by parts we get $$\begin{split}
\int_{{\mathbb{R}}^{d}}(|u_\eps(t) &|^{2}  -|u^\eps_{0}|^{2})dx+2 \int_0^t \int_{{\mathbb{R}}^d}\int_0^{u_{\eps}(\mx,t')} {\mathrm{div}}\mff_\eps(\mx,\lambda) \,d\lambda d\mx dt'
\\&=-2\eps\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla u_\eps|^{2}d\mx dt'+\delta\int_{{\mathbb{R}}^{d}}\triangle (\eta(u_\eps(t,\mx))-\eta(u^\eps_0(\mx))d\mx\\
&-2\delta \int_{{\mathbb{R}}^{d}}\sum_{j=1}^{d}(|{\partial}_{x_j} u_{\eps}(t,\mx)|^{2}-|{\partial}_{x_j} u^\eps_{0}(\mx)|^{2}(\mx))d\mx-2\delta\int_{0}^{t}\int_{{\mathbb{R}}^{d}}\triangle u_{\eps} {\partial}_t u_{\eps}d\mx dt'\\
&=-2\eps\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla u_\eps|^{2}d\mx dt'-2\delta \int_{{\mathbb{R}}^{d}}\sum_{j=1}^{d}(|{\partial}_{x_j}u_{\eps}(t,\mx)|^{2}-|{\partial}_{x_j}u^\eps_{0}|^{2}(\mx))d\mx\\
&+\delta\sum_{j=1}^{d}\int_{{\mathbb{R}}^{d}}\int_{0}^{t}\eta({\partial}_{x_j}u_{\eps})_{t}dtdx=-2\eps\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla u_\eps|^{2}\,d\mx dt'\\
&-\delta \int_{{\mathbb{R}}^{d}}\sum_{j=1}^{d}|{\partial}_{x_j} u_{\eps}(t,\mx)|^{2}-|{\partial}_{x_j} u^\eps_{0}(\mx)|^{2}(\mx))d\mx
\end{split}$$ i.e., $$\label{main-eq2}
\begin{split}
\int_{{\mathbb{R}}^{d}}|u_\eps(t)|^{2}dx+2\eps\int_{0}^{t} &\int_{{\mathbb{R}}^{d}}|\nabla u_\eps|^{2}dxdt'+\delta \sum_{j=1}^{d}\int_{{\mathbb{R}}^{d}}|{\partial}_{x_j}u_{\eps}|^{2}d\mx\\
&=\int_{{\mathbb{R}}^{d}}|u_{0}^\eps|^{2}dx+\delta \sum_{j=1}^{d}\int_{{\mathbb{R}}^{d}}|{\partial}_{x_j} u^\eps_{0}(\mx)|^{2}(x)dx\\
&\hspace*{3em}-2 \int_0^t \int_{{\mathbb{R}}^d}\int_0^{u_{\eps}(\mx,t')} {\mathrm{div}}\mff_\eps(\mx,\lambda) d\lambda d\mx dt'.
\end{split}$$ The relations , and are easy consequences of .

To show , let us now differentiate equation with respect to $x_{i}$ and then multiply it by $2{\partial}_{x_i}u_{\eps}.$ Since the right hand side of is linear in $u_\eps$, the calculation for that side will remain the same and after integrating the equation over $[0,t]\times {\mathbb{R}}^d$ and applying integration by parts, we get $$\begin{split}
&\int_{{\mathbb{R}}^{d}}(|{\partial}_{x_i}u_{\eps}(t,\mx)|^{2}-|{\partial}_{x_i}u^\eps_{0}(\mx)|^{2})d\mx-2\sum_{k=1}^{d}\int_{0}^{t}\int_{{\mathbb{R}}^{d}} {\partial}_\lambda f^\eps_{k}(\mx, u_\eps) \, {\partial}_{x_k} u_\eps \, {\partial}_{x_i}^2 u_{\eps}d\mx dt\\
&-2\int_0^t \int_{{\mathbb{R}}^d} {\partial}_{x_i}^2 u_\eps \, {\mathrm{div}}\mff^\eps(\mx, \lambda)\big{|}_{\lambda=u_\eps} d\mx dt =-2\eps\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla {\partial}_{x_i}u_\eps|^{2}dxds\\
&-\delta \int_{{\mathbb{R}}^{d}}\sum_{j=1}^{d}(|{\partial}_{x_{i}x_j}u_\eps(t,\mx )|^{2}-|{\partial}_{x_{i}x_j}u^\eps_{0}|^{2})d\mx,
\end{split}$$ i.e.,

$$\begin{split}
&\int_{{\mathbb{R}}^{d}}|{\partial}_{x_i} u_{\eps}(t,\mx)|^{2} d\mx+2\eps\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla {\partial}_{x_i}u_\eps |^{2}d\mx dt'+\delta \int_{{\mathbb{R}}^{d}} |\nabla {\partial}_{x_{i}} u_\eps(t,\mx)|^{2}d\mx\\
&=\int_{{\mathbb{R}}^{d}}| {\partial}_{x_i}u^\eps_{0}(\mx)|^{2}d\mx+\delta \int_{{\mathbb{R}}^{d}} |\nabla {\partial}_{x_{i}}u^\eps_{0}(\mx)|^{2}d\mx 
+2 \int_0^t \int_{{\mathbb{R}}^d} {\partial}_{x_i}^2 u_\eps \, {\mathrm{div}}\mff^\eps (\mx, \lambda)\big{|}_{\lambda=u_\eps} d\mx dt' \\ &\hspace*{6em}+2\sum_{k=1}^{d}\int_{0}^{t}\int_{{\mathbb{R}}^{d}}{\partial}_\lambda f^\eps_{k}(\mx,u_\eps) {\partial}_{x_k}u_\eps  {\partial}_{x_i}^2 u_{\eps}dxdt' \\
&\leq \|{\partial}_{x_i}u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}+\delta \|\nabla {\partial}_{x_i} u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}   
+ \frac{\eps}{2} \int_{0}^{t}\int_{{\mathbb{R}}^{d}} |{\partial}_{x_{i}}^2 u_{\eps}|^2 d\mx dt'
\\
&  \hspace*{6em} +\frac{2}{\eps} \| {\mathrm{div}}\mff^\eps(\mx, \lambda)\big{|}_{\lambda=u_\eps}\|^2_{L^2([0,T]\times {\mathbb{R}}^d)}
+\frac{d \|{\partial}_\lambda \mff_\eps\|^2_{\infty}}{\eps} \int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla u_{\eps}|^2 d\mx dt'\\
 & \hspace*{12em} +\eps  \int_{0}^{t}\int_{{\mathbb{R}}^{d}} |{\partial}_{x_{i}}^2 u_{\eps}|^2 d\mx dt'
 \\ 
& \leq  \|{\partial}_{x_i}u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}+\delta \|\nabla {\partial}_{x_i} u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}  +\frac{d \|{\partial}_\lambda \mff_\eps\|^2_{\infty}}{\eps} \|\nabla u_{\eps}\|^2_{L^2([0,T],L^2({\mathbb{R}}^d))}  \\ 
&\hspace*{5em}+ \frac{3\eps}{2}  \int_{0}^{t}\int_{{\mathbb{R}}^{d}} |\nabla {\partial}_{x_{i}} u_{\eps}|^2 d\mx dt' +\frac{2}{\eps} \| {\mathrm{div}}\mff^\eps(\mx, \lambda)\big{|}_{\lambda=u_\eps}\|^2_{L^2([0,T]\times {\mathbb{R}}^d)}
\label{est-d-2} 
\end{split}$$

where $\|{\partial}_\lambda \mff_\eps\|_{\infty}=\max_{1\leq j\leq n}\|{\partial}_{\lambda} f^\eps_{j}\|_{L^\infty({\mathbb{R}}^{d+1})}$. Multiplying the above inequality by $\delta$, we have

$$\begin{split}
&\delta \int_{{\mathbb{R}}^{d}}|{\partial}_{x_i}u_{\eps}(t,\mx)|^{2}d\mx +\frac{\eps \delta}{2} 
\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla {\partial}_{x_i} u_{\eps}|^{2}d\mx dt + \delta^2 \int_{{\mathbb{R}}^{d}}|\nabla {\partial}_{x_{i}}u_{\eps}(t,\mx)|^{2}d\mx\\
& \leq \delta \|{\partial}_{x_i}u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}+\delta^2 \|\nabla {\partial}_{x_i} u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}  +
\frac{d \delta}{\eps} \|{\partial}_\lambda \mff_\eps\|^2_{\infty} \|\nabla u_{\eps}\|^2_{L^2([0,T],L^2({\mathbb{R}}^d))} \\&
+\frac{2\delta }{\eps} \| {\mathrm{div}}\mff^\eps(\mx, \lambda)\big{|}_{\lambda=u_\eps}\|^2_{L^2([0,T]\times {\mathbb{R}}^d)} \\
& \leq \delta \|{\partial}_{x_i}u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}+\delta^2 \|\nabla {\partial}_{x_i} u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})} 
+\frac{2T \delta }{\eps} \| {\mathrm{div}}\mff^\eps(\mx, \lambda)\|^2_{L^2({\mathbb{R}}^d, L^\infty ({\mathbb{R}}))} \\
&+\frac{d\delta}{2\eps^2} 
\|{\partial}_\lambda \mff_\eps\|^2_{\infty} ( \|u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}+\delta\|\nabla u_{0}^{\eps}\|^2_{L^{2}({\mathbb{R}}^{d})}
+2 T \| {\mathrm{div}}\mff_\eps(\mx,\lambda) \|_{L^{1}({\mathbb{R}}^{d}\times {\mathbb{R}})}),
\end{split}$$

where in the last inequality we used equation . From here, immediately follows.

We next aim to prove inequality . To this end we multiply equation by ${\partial}_t u_{\eps}$ and and then proceed similarly to the above. Namely,

$$\label{main-eq11}
\begin{split}
({\partial}_t u_{\eps})^{2}+\sum_{j=1}^{d} {\partial}_t u_{\eps} \, {\partial}_\lambda f_j^\eps(\mx,u_\eps) {\partial}_{x_j}u_{\eps} + {\partial}_t u_{\eps} \, &{\mathrm{div}}\mff^\eps(\mx,\lambda)|_{\lambda=u_\eps} \\
&=\eps {\partial}_t u_{\eps}\triangle u_\eps+\delta {\partial}_t u_{\eps}\, \triangle {\partial}_t u_{\eps},
\end{split}$$

and we integrate equation with respect to time and space.

Integration by part gives

$$\begin{aligned}
\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|{\partial}_t u_{\eps}|^{2}&d\mx dt'+\sum_{j=1}^{d}\int_{0}^{t}\int_{{\mathbb{R}}^{d}}{\partial}_t u_{\eps}{\partial}_\lambda f_j^\eps (\mx, u_\eps) {\partial}_{x_j} u_{\eps}d\mx dt' \\
&+\int_{0}^{t}\int_{{\mathbb{R}}^{d}} {\partial}_t u_{\eps} \, {\mathrm{div}}\mff^\eps(\mx,\lambda)|_{\lambda=u_\eps} d\mx d\lambda \\
&=\eps\int_{0}^{t}\int_{{\mathbb{R}}^{d}}\triangle u_\eps \, {\partial}_t u_{\eps} d\mx dt'+\delta\int_{0}^{t}\int_{{\mathbb{R}}^{d}}{\partial}_t u_{{\eps}}\triangle {\partial}_t u_{\eps}d\mx dt'
\nonumber \\
&=-\sum_{j=1}^{d}\int_{0}^{t}\eps\int_{{\mathbb{R}}^{d}}{\partial}_{x_j} u_{\eps}{\partial}_{tx_j}u_{\eps} d\mx dt'-\delta\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla {\partial}_t u_{\eps}|^{2}d\mx dt',
\nonumber\end{aligned}$$

i.e.

$$\begin{split}
\label{est-51}
&\int_{0}^{t}\int_{{\mathbb{R}}^{d}} |{\partial}_t u_{\eps}|^{2}d\mx dt'+\delta\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla {\partial}_t u_{\eps}|^{2}d\mx dt'=-\eps\sum_{j=1}^{d}\int_{0}^{t}\int_{{\mathbb{R}}^{d}}{\partial}_{x_j} u_{\eps} {\partial}_{tx_j} u_{\eps}d\mx dt'\\
&- \sum_{j=1}^{d}\int_{0}^{t}\int_{{\mathbb{R}}^{d}}\left({\partial}_t u_{\eps}{\partial}_\lambda f_j^\eps (\mx, u_\eps) {\partial}_{x_j} u_{\eps}
- {\partial}_t u_{\eps} \, {\partial}_{x_j} f_j^\eps(\mx,\lambda)|_{\lambda=u_\eps} \right) d\mx dt'    \\
&\leq \| {\partial}_\lambda \mff \|_{\infty}\sum_{j=1}^{d}\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|{\partial}_t u_{\eps}| \, |{\partial}_{x_j} u_{\eps}|d\mx dt' \\
& \hspace*{1em}+ \| {\partial}_t u_\eps \|_{L^2([0,T]\times {\mathbb{R}}^d)} \sqrt{T} \|{\mathrm{div}}\mff_\eps \|_{L^2({\mathbb{R}}^d,L^\infty ({\mathbb{R}}))}  -\frac{\eps}{2}\sum_{j=1}^{d}\int_{0}^{t}\int_{{\mathbb{R}}^{d}}(|{\partial}_{x_j} u_{\eps}|^{2})_t d\mx dt'   \\
&\leq \frac{3}{4}\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|{\partial}_t u_{\eps}|^{2} d\mx dt'+\frac{d \|{\partial}_\lambda \mff_\eps\|_{\infty}^{2}}{2}\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla u_\eps|^{2}d\mx dt'  \\
&+T \|{\mathrm{div}}\mff_\eps \|^2_{L^2({\mathbb{R}}^d, L^\infty({\mathbb{R}}))}
-\frac{\eps}{2}\int_{{\mathbb{R}}^{d}}|\nabla u_\eps(t)|^{2}dx+\frac{\eps}{2}\int_{{\mathbb{R}}^{d}}|\nabla u^\eps_{0}|^{2}d\mx .
 \end{split}$$

By multiplying by $\delta$ and taking into account , we conclude $$\begin{aligned}
&\frac{\delta}{4}\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|{\partial}_t u_{\eps}|^{2} d\mx dt'+\frac{\eps \delta}{2}\int_{{\mathbb{R}}^{d}}|\nabla u_\eps (t)|^{2}d\mx+\delta^{2} \int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla {\partial}_t u_{\eps}|^{2}d\mx dt' \\
& \nonumber \leq \frac{\delta d\| {\partial}_\lambda \mff_\eps\|_{\infty}^{2}}{2}\int_{0}^{t}\int_{{\mathbb{R}}^{d}}|\nabla u_\eps|^{2}dxdt'+\frac{\delta\eps}{2}\int_{{\mathbb{R}}^{d}}|\nabla u^\eps_{0}|^{2} d\mx+\delta T 
\|{\mathrm{div}}\mff_\eps \|^2_{L^2({\mathbb{R}}^d, L^\infty ({\mathbb{R}}))} \nonumber \\
&\leq \frac {\delta d\|{\partial}_\lambda \mff_\eps\|_{\infty}^{2}}{4\eps}
\|u^\eps_{0}\|_{L^{2}({\mathbb{R}}^{d})}^{2}+
\left(\frac {\delta^2 d\|{\partial}_\lambda \mff_\eps\|_{\infty}^{2}}{4\eps}+\frac{\delta  \eps}{2}\right)\|\nabla u^\eps_{0}\|_{L^{2}({\mathbb{R}}^{d})}^{2} \nonumber \\
& +\frac {\delta d\|{\partial}_\lambda \mff_\eps\|_{\infty}^{2}T}{2\eps}  \|{\mathrm{div}}\mff_\eps \|_{L^1({\mathbb{R}}^d\times {\mathbb{R}})}+\delta T \|{\mathrm{div}}\mff_\eps \|^2_{L^2({\mathbb{R}}^d,L^\infty({\mathbb{R}}))},
\nonumber\end{aligned}$$which gives the inequality .

We shall now derive the kinetic formulation for and then show that we can apply the velocity averaging result on this form of the equation.

[**(Kinetic formulation of )**]{} Assume that the function $u_\eps$ satisfies . Then, the function $h_\eps(t,\mx,\lambda)={\rm sgn}(u_\eps(t,\mx)-\lambda)$ solves a linear PDE of the form $$\begin{aligned}
\label{kinetic}
&{\partial}_t h_\eps +{\mathrm{div}}(h_\eps {\partial}_\lambda \mff_\eps (\mx,\lambda) )-{\partial}_\lambda  (h_\eps {\mathrm{div}}\mff_\eps(\mx,\lambda))={\partial}_\lambda {\mathrm{div}}(G_1^\eps)+ {\partial}^2_\lambda G^\eps_2,\end{aligned}$$ where

-   ${\mathrm{div}}G_1^\eps \to 0$ in $W^{-1,r}_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d\times {\mathbb{R}})$ for any $r\in [1,2)$;

-   $(G^\eps_2)$ is bounded in the space of measures ${\cal M}_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})$ and thus strongly precompact in $W_{loc}^{-1,r}({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})$ for any $r\in [1,\frac{d+2}{d+1}]$.

Before we start, notice that we require merely local convergence (see conditions (i) and (ii) of the lemma), which means that we can always choose $\eps$ small enough so that $K_\eps \equiv 1$ on the fixed relatively compact set on which we derive the estimates.

Therefore, it is enough to prove the lemma under the simplifying assumption $$\mff_\eps=\mff\star\omega_{n(\eps)}.$$ Also, in order to avoid proliferation of symbols, we will notationally adhere to integration over the entire space ${\mathbb{R}}^d\times {\mathbb{R}}$, with the implicit understanding we are in fact integrating over the support of the corresponding test function (see ).

Then, note that for any $\eta \in C^3({\mathbb{R}})$ such that $\eta'\in C^2_c({\mathbb{R}})$ we have (recall that $h_\eps(t,\mx,\lambda)={\rm sgn}(u_\eps(t,\mx)-\lambda)$) $$\label{kinetic-1}
\eta(u_\eps)  = \frac{1}{2} \int_{\mathbb{R}}\eta'(\lambda) h_\eps(t,\mx,\lambda)  d\lambda.$$ Moreover, for any function $f\in L^1({\mathbb{R}})$ $$\label{kinetic-2}
\int_{-\infty}^{u_\eps} f(\lambda) d\lambda=\frac{1}{2}\int_{\mathbb{R}}f(\lambda )  h_\eps(t,\mx,\lambda)  d\lambda +\frac{1}{2}\int_{\mathbb{R}}f(\lambda) d\lambda.$$ To proceed, we multiply equation by $\eta'(u_\eps)$. We get

$$\begin{aligned}
\label{kinetic-new}
&\partial_{t} \eta(u^{\eps})+{\mathrm{div}}\int_{-\infty}^{u_\eps} \eta'(\lambda) {\partial}_\lambda \mff_\eps(\mx,\lambda) d\lambda+\int_{-\infty}^{u_\eps} \eta''(\lambda){\mathrm{div}}\mff_\eps(\mx,\lambda) d\lambda\\&
=\sum_{j=1}^{d}\partial_{x_j}\left(\eps\eta'(u_{\eps}){\partial}_{x_j}u_{\eps}+\delta\eta'(u_{\eps}){\partial}_{t,x_{j}}u_{\eps}\right) \nonumber \\
&-\eps \eta''(u_{\eps})\sum_{j=1}^{d}|{\partial}_{{x_j}}u_\eps|^{2}-\delta\eta''(u_\eps)\sum_{j=1}^{d}{\partial}_{x_{j}}u_{\eps}{\partial}_{t,x_{j}}u_{\eps}.
\nonumber \end{aligned}$$

Now, we apply and and rewrite the latter equation in the form

$$\begin{aligned}
&\frac{1}{2}
\partial_{t} \int_{{\mathbb{R}}} \eta'(\lambda) h_\eps d\lambda +\frac{1}{2}{\mathrm{div}}\int_{\mathbb{R}}\eta'(\lambda) 
{\partial}_\lambda \mff_\eps(\mx,\lambda) h_\eps d\lambda+\frac{1}{2}\int_{\mathbb{R}}\eta''(\lambda){\mathrm{div}}\mff_\eps(\mx,\lambda) h_\eps d\lambda \\&
=\sum_{j=1}^{d}\partial_{x_j}\left(\eps \int_{\mathbb{R}}\eta''(\lambda) h_\eps d\lambda {\partial}_{x_j}u_{\eps}+\delta\int_{\mathbb{R}}\eta''(\lambda) h_\eps d\lambda {\partial}_{t,x_{j}}u_{\eps}\right)\\
&-\eps \int_{\mathbb{R}}\eta'''(\lambda) h_\eps d\lambda\sum_{j=1}^{d}|{\partial}_{{x_j}}u_\eps|^{2}-\delta\int_{\mathbb{R}}\eta'''(\lambda) h_\eps d\lambda\sum_{j=1}^{d}{\partial}_{x_{j}}u_{\eps}{\partial}_{t,x_{j}}u_{\eps}.\end{aligned}$$

If we rewrite the latter in the variational formulation, we get for every test function $\phi \in C_c^2({\mathbb{R}}^+\times {\mathbb{R}}^d)$ $$\begin{aligned}
\label{20}
&-\int_{{\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}}} \left( \partial_{t}\phi (t,\mx)  \eta'(\lambda) +  \nabla \phi \eta'(\lambda) \cdot  {\partial}_\lambda \mff_\eps(\mx,\lambda)\right) h_\eps dt d\mx d\lambda\\
&+\int_{{\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}}}  \eta''(\lambda) \, \phi \, {\mathrm{div}}\mff_\eps (\mx ,\lambda) h_\eps dt d\mx  d\lambda
\nonumber
\\&
=-2\sum_{j=1}^{d}\int_{{\mathbb{R}}+\times {\mathbb{R}}^d \times {\mathbb{R}}}h_\eps \partial_{x_j}\phi\eta''(\lambda)  \left(\eps  {\partial}_{x_j}u_{\eps}+\delta {\partial}_{t,x_{j}}u_{\eps}\right) dt d\mx d\lambda\nonumber\\
&- 2\int_{{\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}}}\phi(t,\mx) \eta'''(\lambda) h_\eps \left( \eps \sum_{j=1}^{d}|{\partial}_{{x_j}}u_\eps|^{2}-\delta \sum_{j=1}^{d}{\partial}_{x_{j}}u_{\eps}{\partial}_{t,x_{j}}u_{\eps} \right) dt d\mx d\lambda\nonumber\\
&=:-\int_{{\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}}} \left( \Gamma_{1}^{\eps}+\Gamma_{2}^{\eps}\right) \cdot \nabla \phi(t,\mx) \eta''(\lambda)dt d\mx d\lambda\nonumber\\
& \hspace*{8em} -\int_{{\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}}} \left(\Gamma_{3}^{\eps}+\Gamma_{4}^{\eps}\right) \phi(t,\mx) \eta'''(\lambda) dt d\mx d\lambda.\nonumber\end{aligned}$$ Now, since $\eta'\in C^2_c({\mathbb{R}})$ is arbitrary, if we put $$\label{G}
G_1^\eps =\Gamma_{1}^{\eps}+\Gamma_{2}^{\eps} \ \ {\rm and} \ \ G_2^\eps=\Gamma_{3}^{\eps}+\Gamma_{4}^{\eps}$$ we arrive at .

Next, we prove (i) and (ii). To begin with, we show that ${\mathrm{div}}\Gamma_{1}^{\eps}\rightarrow 0\enspace\mbox{as} \enspace \eps\rightarrow 0^{+}$ in $H^{-1,r}_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d\times {\mathbb{R}})$, $r\leq 2$.

We have for any $\phi \in C_c^{1}({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})$ supported in the hypercube with side length $M$ centered at zero:

$$\begin{split}
|\left<{\mathrm{div}}\Gamma_{1}^{\eps},\phi \right>|&= \left|\sum_{j=1}^{d}\int_{0}^{T}\int_{{\mathbb{R}}^{d+1}}\eps h_\eps u_{x_j}^{\eps}\phi_{x_j}(t,\mx,
\lambda)d\mx d\lambda dt\right| \\
&\leq \eps M^{1/2} \| \phi\|_{H^1({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})} \, \| \nabla u^\eps \|_{L^2([0,T]\times {\mathbb{R}}^d)}.
\end{split}$$

To proceed further, we use . For estimating the term $\|{\mathrm{div}}(\mff_\varepsilon)\|_{L^1({\mathbb{R}}^{d+1})}$ appearing in we use approximations of the form $\mff_\eps :=\mff \star \omega_{n(\eps )}$ as in with $n(\eps)$ satisfying . Since $\delta=o(\eps^2)$, the $\mff_\eps$ will converge to $\mff $ as $\eps \to 0$ in $L^r_{loc}({\mathbb{R}}^d\times {\mathbb{R}})$ for any $1\le r \le p$. Using (C2) and (C3), we get

$$\label{divfepsL1}
 \begin{split}
\|{\mathrm{div}}(\mff_\varepsilon)\|_{L^1({\mathbb{R}}^{d+1})}&= \| |{\mathrm{div}}(\mff)|\star \omega_{n(\eps)}\|_{L^1({\mathbb{R}}^{d+1})} \\
&\leq 
\left\| \frac{\mu(\mx)}{1+|\lambda|^{1+\beta}}\star \omega_{n(\eps)}(\mx,\lambda) \right\|_{L^1({\mathbb{R}}^{d+1})} 
\\&\leq C \mu({\mathbb{R}}^{d}) \| \omega_{n(\eps)} \|_{L^1({\mathbb{R}}^{d+1})} \leq C,
\end{split}$$

where here and below, $C$ denotes a generic constant that may alter from line to line. Thus (keeping in mind and ), from we conclude $$\begin{aligned}
|\left<{\mathrm{div}}\Gamma_{1}^{\eps},\phi  \right>| \leq \sqrt{\eps} (n(\eps)^{-1}\sqrt{\delta}+1)C   
\|\phi\|_{H^{1}({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}}})\rightarrow 0^{+},\enspace as\enspace \eps\rightarrow 0^{+}\end{aligned}$$ by .

For the next estimate we have $$\label{divgammaeps2}
\begin{split}
|\left<{\mathrm{div}}\Gamma_{2}^{\eps},\phi \right>|\leq 2 M^{1/2}
\|\phi\|_{H_{0}^{1}([0,T]\times{\mathbb{R}}^{d+1})} \delta \|\nabla u^{\eps}_{t}\|_{L^{2}([0,T],L^{2}({\mathbb{R}}^{d}))}.
\end{split}$$

Comparing with , we see that it remains to estimate the terms $\|{\partial}_\lambda \mff_\eps\|_{\infty} $ and $\|{\mathrm{div}}(\mff_\varepsilon)\|_{L^2({\mathbb{R}}^{d},L^\infty ({\mathbb{R}}))}$. For the first one we have $$\label{palambdafepsinfty}
\begin{split}
\|{\partial}_\lambda &\mff_\eps \|_{\infty} = \sup_{i} \| ({\partial}_\lambda f_i)\star \omega_{n(\eps)}\|_{L^\infty({\mathbb{R}}^{d+1})}\\ 
&= \sup_{i} \| ({\partial}_\lambda f_i)\star \omega_{n(\eps)}\|_{L^\infty({\mathbb{R}}^{d}, L^\infty({\mathbb{R}}))} \\ &\leq \sup_{i} \| {\partial}_\lambda f_i\|_{L^1({\mathbb{R}}^{d}, L^\infty({\mathbb{R}}))}  \| \omega_{n(\eps)} \|_{L^\infty({\mathbb{R}}^{d+1})}  \leq C n(\eps)^{-d-1} 
\end{split}$$ by (C3). For the second one we obtain from (C2) and the Young inequality $$\label{divfepsL2}
\begin{split}
\|{\mathrm{div}}(\mff_\varepsilon)\|_{L^2({\mathbb{R}}^{d},L^\infty ({\mathbb{R}}))}&\leq 
\| \frac{\mu(\mx)}{1+|\lambda|^{1+\beta}}\star \omega_{n(\eps)}(\mx,\lambda) \|_{L^2({\mathbb{R}}^{d},L^\infty({\mathbb{R}}))}\\
& =  \| \mu\star \omega_{n(\eps)}^{(1)}(\mx)\frac{1}{1+|\lambda|^{1+\beta}}\star \omega^{(2)}_{n(\eps)}(\lambda) \|_{L^2({\mathbb{R}}^{d},L^\infty({\mathbb{R}}))} \\
&\leq C n(\eps)^{-1}\|\mu \star \omega^{(1)}_{n(\eps)} \|_{L^2({\mathbb{R}}^d)}  
\\ &\leq C n(\eps)^{-1}\mu({\mathbb{R}}^d) \| \omega^{(1)}_{n(\eps)} \|_{L^2({\mathbb{R}}^d)} \leq C 
n(\eps)^{-d/2-1}.
\end{split}$$ Combining , , , and we see that $$\begin{aligned}
\nonumber
|\langle{\mathrm{div}}& \Gamma_{2}^{\eps},\phi \rangle| 
\leq \\ 
& C \delta^{1/2} \left(\frac{1}{\sqrt{\eps} n(\eps)^{d+1}}
+ \Big(\frac{\sqrt{\delta}}{\sqrt{\eps}n(\eps)^{d+1}}+\sqrt{\eps}\Big)
\frac{1}{n(\eps)} + \frac{1}{\sqrt{\eps}n(\eps)^{d+1}}
+\frac{1}{n(\eps)^{d/2 + 1}} \right) \nonumber \\
&\cdot\|\phi\|_{H_{0}^{1}([0,T]\times {\mathbb{R}}^{d+1})}
\rightarrow 0^{+},\enspace as\enspace \eps\rightarrow 0^{+}.
\label{G2-new}\end{aligned}$$ by . This implies claim (i) of the current lemma. Next, using , , , and we estimate $$\label{G4-new}
\begin{split}
|\left<\Gamma_{4}^{\eps},\phi \right>|& \leq \delta\left|\sum_{j=1}^{d}\int_{0}^{T}\int_{[-M,M]^{d+1}}|u_{x_j}^{\eps}| \, |u_{t,x_{j}}^{\eps}| \, |\phi(t,\mx,\lambda)|d\mx dt d\lambda \right|\\
&\leq M \delta \|\phi\|_{\infty}\sum_{j=1}^{d}\int_{0}^{T}\int_{{\mathbb{R}}^{d}}|u_{x_j}^{\eps}||u_{t,x_{j}}^{\eps}|dx dt\\
&\leq M \delta \|\phi\|_{\infty} \|\nabla u^{\eps}\|_{L^{2}([0,T]\times {\mathbb{R}}^{d})} \|\nabla u_{t}^{\eps}\|_{L^{2}([0,T]\times {\mathbb{R}}^{d})}\\
& \le C \sqrt{\delta} \|\phi\|_\infty  \Big(\frac{1}{\sqrt{\eps}} +
\frac{{\sqrt{\delta}}}{\sqrt{\eps} n(\eps)}\Big) \Big[\frac{1}{\sqrt{\eps}n(\eps)^{d+1}}+ \\
& + \Big(\frac{\sqrt{\delta}}{\sqrt{\eps}n(\eps)^{d+1}}
+\sqrt{\eps}\Big)\frac{1}{n(\eps)} + 
\frac{1}{\sqrt{\eps}n(\eps)^{d+1}} +  \frac{1}{n(\eps)^{d/2+1}}
\Big],  
\end{split}$$ which goes to $0$ as $\eps\to 0+$ by .

Finally, it remains to estimate $\Gamma_3^\eps$. Using , and we obtain $$\begin{split}
|\left<\Gamma_{3}^{\eps},\phi \right>|& \leq \eps  \left|\sum_{j=1}^{d}\int_{0}^{T}\int_{[-M,M]^{d+1}} |u_{x_j}^{\eps}(t,\mx)|^2 \phi(t,\mx,\lambda) \, d\mx dt d\lambda\right|\\
&\leq  \eps M \|\phi\|_{\infty} \|\nabla u^{\eps}\|^2_{L^{2}([0,T]\times {\mathbb{R}}^{d})} \\
&\leq C \|\phi\|_{\infty} (1+\frac{\delta}{n(\eps)^2}) \leq  C \|\phi\|_{\infty}.
\end{split}$$

From the estimates given above, we see that ${\mathrm{div}}\Gamma_1^\eps$ and ${\mathrm{div}}\Gamma_2^\eps$ converge strongly to zero in $H^{-1}_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})$ implying that they converge strongly in $W^{-1,r}_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})$, $r\in [1,2)$, as well.

As for $\Gamma_3^\eps$ and $\Gamma_4^\eps$, they are locally bounded in the space of Radon measures (which we denote it by ${\cal M}({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})$) and thus, they are precompact in $W^{-1,r}_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})$, $r\in \big[1,\frac{d+2}{d+1}\big]$ [@evans Th. 1.3.2]. This implies condition (ii) from the Lemma.

According to the previous theorem, we see that we can apply the velocity averaging lemma Theorem \[va\]. Indeed, we have:

Any solution $h_\eps$ of satisfies $$\label{for-va}
{\partial}_t h_\eps +{\mathrm{div}}(h_\eps {\partial}_\lambda \mff (\mx,\lambda) )={\partial}_\lambda {\mathrm{div}}(G_1^\eps)+ {\partial}^2_\lambda G^\eps_2 + {\mathrm{div}}(G^\eps_3)+{\partial}_\lambda  (h_\eps {\mathrm{div}}\mff_\eps(\mx,\lambda)),$$ where $(G_j^\eps)$, $j=1,3$, are strongly precompact in $L_{loc}^r({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})$ while $(G^\eps_2)$ and $(h_\eps {\mathrm{div}}\mff_\eps(\mx,\lambda))$ are strongly precompact in $W^{-1,r}_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})$, $r\in \big[1,\frac{d+2}{d+1}\big]$.

It is enough to rewrite in the form $$\begin{split}
& {\partial}_t h_\eps +{\mathrm{div}}(h_\eps {\partial}_\lambda \mff (\mx,\lambda) )
={\partial}_\lambda {\mathrm{div}}(G_1^\eps)+ {\partial}^2_\lambda G^\eps_2  \\& 
\hspace*{5em}
+{\mathrm{div}}(h_\eps ({\partial}_\lambda \mff (\mx,\lambda)-{\partial}_\lambda \mff_\eps (\mx,\lambda)))+{\partial}_\lambda  (h_\eps {\mathrm{div}}\mff_\eps(\mx,\lambda)),
\end{split}$$ and denote $$G_3^\eps=h_\eps ({\partial}_\lambda \mff (\mx,\lambda)-{\partial}_\lambda \mff_\eps (\mx,\lambda)).$$ Clearly, since $\mff_\eps=K_\eps\cdot\mff \star \omega_{n(\eps)}$, we have $G_3^\eps \to 0$ as $\eps\to 0$ in $L_{loc}^p({\mathbb{R}}^d\times {\mathbb{R}})$ (were $p$ is given in (C1)). On the other hand, according to condition (C2), we see that $(h_\eps {\mathrm{div}}\mff_\eps(\mx,\lambda)))$ is bounded in the space of measures and thus strongly precompact in $ W^{-1,r}_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}})$. This concludes the proof.

Now, we are ready to use Theorem \[va\]. In fact, we may choose $r\in \big(1,\frac{d+2}{d-1}\big]$ and $s\ge 2$ in such a way that $$1+\frac{1}{s} - \frac{1}{r} <1$$ and then set $\bar p := (1+\frac{1}{s} - \frac{1}{r})^{-1}$. Then fixing $\bar{p}'$ such that $\frac{1}{\bar p} + \frac{1}{\bar p'}=1$, the assumptions of Theorem \[main-theorem\] precisely allow to apply Theorem \[va\] (noting that $(h_\eps)$ is bounded in any $L^s_{loc}$, hence contains a weakly convergent subsequence). More precisely, one easily checks that any $\bar p\in \big[\frac{2d+4}{d+4},\infty\big)$, and thereby any $\bar p' \in \big(1,2+\frac{4}{d}\big]$ can be obtained in this way. This is the reason for the specific assumption in Theorem \[main-theorem\].

Consequently, for any sequence $\eps_n\searrow 0$ and any $\rho \in L^2({\mathbb{R}})$, setting $h_n:=h_{\eps_n}$ and $u_n:=u_{\eps_n}$, the sequence $\left(\int_{\mathbb{R}}\rho(\lambda) h_n(t,\mx,\lambda) d\lambda \right) $ is strongly precompact in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d)$. As we shall see in the next theorem, this implies strong convergence of the sequence $(u_n)$. We first need the following auxiliary result.

\[lconv-1\] Assume that the sequence $(u_n)$ is bounded in $L^p(\Omega)$, $\Omega \Subset {\mathbb{R}}^d$, $d\in {\mathbb N}$, for some $p\geq 1$. Define $$\Omega_n^l=\{ \mx\in\Omega:\, |u_{n}(\mx)| > l \}.$$ Then $$\label{conv1}
\lim\limits_{l\to \infty} \sup\limits_{n\in {\mathbb N}} {\rm
meas}(\Omega_n^l) =0.$$

Since $(u_{n})$ is bounded in $L^p(\Omega)$ and thus in $L^1(\Omega)$ as well (since $\Omega \Subset {\mathbb{R}}^m$), we have $$\begin{aligned}
&\sup\limits_{n\in {\mathbb N}} \int_{\Omega }|u_{n}(\mx)| d\mx \geq
\sup\limits_{n\in {\mathbb N}} \int_{\Omega_n^l} l d\mx\, \implies \frac{1}{l}\sup\limits_{k\in {\mathbb N}} \int_{\Omega }|u_{n}(\mx)| d\mx \geq
\sup\limits_{n\in {\mathbb N}} {\rm meas}(\Omega_n^l),\end{aligned}$$ implying after letting $l\to \infty$ here.

With the above notations we finally arrive at:

Under the non-degeneracy conditions , the sequence $(u_n=u_{\eps_n})$ of solutions to , strongly converges along a subsequence in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}}^d)$.

For $l\in {\mathbb N}$ we set $\rho(\lambda)=\chi_{(-l,l)}(\lambda)$, where $\chi_{(-l,l)}$ is the characteristic function of the interval $(-l,l)$. Then by Theorem \[va\] combined with a diagonalization argument there exists a common subsequence, again denoted by $(u_n)$, such that for any $l\in {\mathbb N}$: $$\begin{aligned}
\frac{1}{2}\int_{-l}^l h_n(t,\mx,\lambda) d\lambda&=u_n
\chi_{\{|u_n|\leq l\}}(t,\mx)+l \chi_{\{u_n>l\}}(t,\mx)-l \chi_{\{u_n<-l\}}(t,\mx)\nonumber\\
&=:T_l(u_n) \to u^l \label{3.8first}\end{aligned}$$ as $n\to \infty$ in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}})$. The operators $T_l$ are known as truncation operators [@DHM].

It is not difficult to prove that from here we can conclude about the convergence of $(u_n)$. First, we show that the sequence $(u^l)$ converges strongly in $L^1_{loc}({\mathbb{R}}^+\times {\mathbb{R}})$ as $l\to \infty$.

To this end, let $\Omega\Subset {\mathbb{R}}^+\times {\mathbb{R}}^d$. We claim that $$\label{uni_n}
\lim_l \sup_n \|{T_l(u_{n}) -  u_{n}}\|_{L^1(\Omega)} \to 0\,.$$ In fact, let $$\Omega_n^l=\{ (t,\mx)\in\Omega:\, |u_{n}(t,\mx)| > l \}.$$ Then since (by , , , and ) $(u_{n})$ is bounded in $L^2(\Omega)$, we have

$$\int_\Omega |u_{n}-T_l(u_{n})| dt d\mx
\leq \int_{\Omega_n^l} |u_{n}|dt d\mx  \leq {\rm meas}(\Omega_n^l)^{1/2} \, \|u_n\|_{L^2(\Omega)}\to 0$$ as $l\to \infty$, uniformly with respect to $n$ according to . This gives .

Next, we estimate $$\begin{aligned}
\label{ul} \|u^{l_1}-u^{l_2}\|_{L^1( \Omega)}
&\leq
\|u^{l_1}-T_{l_1}(u_n)\|_{L^1( \Omega)}+\|T_{l_1}(u_n)-u_n\|_{L^1( \Omega)}\nonumber\\
&+\|T_{l_2}(u_n)-u_n\|_{L^1( \Omega)}+\|T_{l_2}(u_n)-u^{l_2}\|_{L^1( \Omega)}\,,\nonumber\end{aligned}$$ which together with and implies that $(u^l)$ is a Cauchy sequence. Thus, there exists $u\in L^1(\Omega)$ such that $$\label{conv2}
u^l \to u \ \ {\rm in} \ \ L^1(\Omega).$$

Now it is not difficult to see that the entire sequence $(u_n)$ converges towards $u$ in $L^1(\Omega)$ as well. Namely, $$\begin{aligned}
   \|u_n-u\|_{L^1(\Omega)} \leq
\|u_n-T_l(u_n)\|_{L^1(\Omega)}+\|T_l(u_n)-u^l\|_{L^1(\Omega)}+\|u^l-u\|_{L^1(\Omega)},\end{aligned}$$ which by the definition of the functions $u^l$, in conjunction with and gives the claim.

We actually proved only item (i) of Theorem \[main-theorem\]. The other two items (item (ii) and item (iii)) can be proven by an adaptation of the proof for item (i). We provide a more precise explanation in the following remark.

\[39\] Derivation of (ii) and (iii) from Theorem \[main-theorem\].

If we additionally assume that $\|{\partial}_\lambda \mff\| \leq C <\infty$ then we can use in the case $\delta={\cal O}(\eps^2)$ and $\frac{\sqrt{\delta}}{\sqrt{\eps} n(\eps)^{d/2+2}}\to 0$ and the considerations thereafter will remain the same. This gives (ii). Indeed, we have to check whether the sequence of equations satisfies the conditions of Theorem \[va\]. To this end, we need to estimate $\Gamma^\eps_j$, $j=1,\dots,4$, appearing in . If we assume $\|{\partial}_\lambda \mff\| \leq C <\infty$ then we can omit and get instead of $$\begin{aligned}
\label{G2-newReg}
|\left<{\mathrm{div}}\Gamma_{2}^{\eps},\phi \right>| 
\leq C \delta^{1/2} \left(\frac{1}{\sqrt{\eps} }
+ \Big(\frac{\sqrt\delta}{\sqrt{\eps}}+\sqrt{\eps}\Big)\frac{1}{n(\eps)} + \frac{1}{\sqrt{\eps}}
+\frac{1}{n(\eps)^{d/2 +1}} \right) \\
\cdot\|\phi\|_{H_{0}^{1}([0,T]\times {\mathbb{R}}^{d+1})}
\rightarrow 0^{+},\enspace as\enspace \eps\rightarrow 0^{+}.
\nonumber\end{aligned}$$ Similarly, instead of , we have $$\label{G4-newReg}
\begin{split}
|\left<\Gamma_{4}^{\eps},\phi \right>|& \leq \delta\left|\sum_{j=1}^{d}\int_{0}^{T}
\int_{[-M,M]^{d+1}}|u_{x_j}^{\eps}| \, |u_{t,x_{j}}^{\eps}| \, |\phi(t,\mx,\lambda)|d
\mx dt d\lambda \right|\\
& \le C \sqrt{\delta} \|\phi\|_\infty \Big(\frac{1}
{\sqrt{\eps}} +
\frac{\sqrt{\delta}}{\sqrt{\eps} n(\eps)}\Big) \Big[\frac{1}{\sqrt{\eps}} 
+ \\
& + \Big(\frac{\sqrt{\delta}}{\sqrt{\eps}}
+\sqrt{\eps}\Big)\frac{1}{n(\eps)} + 
\frac{1}{\sqrt{\eps}}+  \frac{1}{n(\eps)^{d/2 +1}}
\Big]\leq C \|\phi\|_\infty 
\end{split}$$

If we have a regular flux i.e. $\mff \in C^1({\mathbb{R}}^d\times {\mathbb{R}})$ and the conditions $\delta=o(\eps^2)$ and $\frac{\sqrt{\delta}}{\sqrt{\eps} n(\eps)^{d/2+2}}\to 0$ as $\eps\to 0$, we can use Young measures in the way given in [@Lfl-b] to derive the convergence. Actually, in this case the diffusion given by the second order term in will dominate over the dynamic capillarity given by the third order term in and therefore, we will end up with the unique Kruzhkov admissible solution to with the corresponding initial data.

Indeed, it is not difficult to see that under the regularity assumptions on $\mff$ and $\delta=o(\eps^2)$, we have $$\label{regular}
\begin{split}
& -\int_{{\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}}} \left( \Gamma_{1}^{\eps}+\Gamma_{2}^{\eps}\right) \cdot \nabla \phi(t,\mx) \eta''(\lambda)dt d\mx d\lambda \to 0\\
&  -\int_{{\mathbb{R}}^+\times {\mathbb{R}}^d \times {\mathbb{R}}} \Gamma_{4}^{\eps}(t,\mx,\lambda) \phi(t,\mx) \eta'''(\lambda) dt d\mx d\lambda \to 0.
\end{split}$$ Now, denote by $\nu_{(t,\mx)}$ the Young measure (e.g. [@evans]) corresponding to a subsequence of the family $(u^\eps)$. If we let $\eps\to 0$ in along the subsequence for a convex entropy $\eta$, and take into account, we get in the sense of distributions (below, $q(\mx,\lambda)=\int_{-\infty}^\lambda \eta'(v) {\partial}_v \mff(\mx,v) dv$ is the entropy flux)

$$\partial_{t} \int \eta(\lambda)d\nu_{(t,\mx)}(\lambda) +{\mathrm{div}}\int q(\mx,\lambda) d\nu_{(t,\mx)}(\lambda)+\int \int_{-\infty}^{\lambda} \eta''(v){\mathrm{div}}\mff(\mx,v) dv d\nu_{(t,\mx)}(\lambda) \leq 0.$$

Now, we simply rely on the result from [@Dpe] (see also [@KLfl2; @Sze]) to conclude that the Young measure is unique and atomic i.e. of the form $\nu_{(t,\mx)}(\lambda)=\delta(\lambda-u(t,\mx))$ for the unique entropy admissible solution to the underlying conservation law .

[99]{}

B. Andreianov and C. Cancès, [*Vanishing capillarity solutions of Buckley-Leverett equation with gravity in two-rocks medium*]{}, Comput. Geosci., [**17**]{}(3) (2013), 551–572.

B. Andreianov, K. H. Karlsen and N. H. Risebro, [*A theory of $L^1$-dissipative solvers for scalar conservation laws with discontinuous flux*]{}, Arch. Ration. Mech. Anal., [**201**]{} (2011), 27–86.

B. Andreianov, D.  Mitrović, [*Entropy conditions for scalar conservation laws with discontinuous flux revisited*]{}, Ann. Inst. H. Poincaré Anal. Non Lin. [**32**]{} (2015), 1307–1335.

J.Aleksić, D.Mitrović, S.Pilipović, [*Hyperbolic conservation laws with vanishing nonlinear diffusion and linear dispersion in heterogeneous media*]{}, J. Evol. Equ. [**9**]{} (2009), 809–828.

G. Barenblat, I. Zheltov, I. Kochiva, [*Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks*]{}, J. Appl. Math. Mech. 24 (5) (1960) 1286–1303.

T.B. Benjamin, J.L. Bona, J.J. Mahony, [*Model equations for long waves in nonlinear dispersive systems*]{}, Philos. Trans. R. Soc. Lond. Ser. A 272 (1220) (1972) 47–78.

J. Bear, [*Dynamics of Fluids in Porous Media*]{}, American Elsevier Publishing Company, Inc., New York, 1972.

S. E.  Buckley and M. C.  Leverett, [*Mechanism of fluid displacements in sands*]{}, Transactions of the AIME [****]{}146 (1942), 107–116.

Y. Cao, J.X. Yin, C.P. Wang, [*Cauchy problems of semilinear pseudo-parabolic equations*]{}, J. Differential Equations [**246**]{} (2009), 4568–4590.

G. M. Coclite, L. di Ruvo, [*A note on convergence of the solutions of Benjamin-Bona-Mahony type equations*]{}, Nonlinear Anal. Real World Appl. [**40**]{} (2018), 64–81.

G. M. Coclite, L. di Ruvo, [*A singular limit problem for conservation laws related to the Rosenau-Korteweg-de Vries equation*]{}, J. Math. Pures Appl. (9) 107 (2017), 315–335.

G. Crasta, V.  De Cicco, G. De Philippis, F. Ghiraldin, [*Structure of solutions of multidimensional conservation laws with discontinuous flux and applications to uniqueness*]{}, Arch. Ration. Mech. Anal. [**221**]{} (2016), 961–985.

L.  C.  Cueto-Felgueroso, R.  Juanes, [*Stability analysis of a phase-field model of gravity-driven unsaturated flow through porous media*]{}, Physical Review E 79 (2009), 13 pages.

D. A.  DiCarlo, [*Modeling observed saturation overshoot with continuum additions to standard unsaturated theory*]{}, Advances in Water Resources [**28**]{} (2005), 1021–1027.

R.J. DiPerna, [*Measure-valued solutions to conservation laws*]{}, Arch. Rational Mech. Anal. [**88**]{} (1985), 223–270.

R. J. DiPerna, P. L. Lions, Y. Meyer, [*$L^p$-regularity of velocity averages*]{}, Ann. Inst. H. Poincaré Anal. Non Linéaire, [**8**]{} (1991), 271—287.

G. Dolzmann, N. Hungerbuhler, S. Müller, [*Nonlinear elliptic systems with measure valued right-hand side*]{}, Math. Zeitschrift, [**226**]{} (1997), 545–574.

C. J.  van Duijn, L. A.  Peletier and I. S.  Pop, [*A new class of entropy solutions of the Buckley- Leverett equation*]{} SIAM J. Math. Anal. [**39**]{} (2007), 507–536.

L. C. Evans, [*Weak Convergence Methods for Nonlinear Partial Differential Equations*]{}, Regional Conference Series in Mathematics, No. 74. Conference Board of the Mathematical Sciences, 1990.

P. Gérard, [*Microlocal Defect Measures*]{}, Comm. Partial Differential Equations [**16**]{} (1991), 1761–1794.

D. Gilbarg, N.S.Trudinger, Elliptic Partial Differential Equations of Second Order, Springer, 2001.

F. Golse, L. Saint-Raymond, [*Velocity averaging in $L^1$ for the transport equation*]{}, C. R. Acad. Sci. Paris Ser. I Math., [**334**]{} (2002), 557–562.

S.  M.  Hassanizadeh, W. G.  Gray, [*Thermodynamics basis of capillary pressure in porous media*]{}, Water Resour. Res. [**29**]{} (1993), 3389–3405.

S.  M.  Hassanizadeh, W. G.  Gray, Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries, Adv. Water Resor. [**13**]{} (1990), 179–186.

H.Holden, K.H.Karlsen, D.Mitrovic, [*Zero diffusion-dispersion-smoothing limits for a scalar conservation law with discontinuous flux function*]{}, Int. J. Differ. Equ. 2009, Art. ID 279818, 33 pp.

G.M.Homsy, [*Viscous fingering in porous media*]{}, Annu. Rev. Fluid Mech. [**19**]{} (1987), 271– 311.

S.Hwang, A.E.Tzavaras, [*Kinetic decomposition of approximate solutions to conservation laws: application to relaxation and diffusion-dispersion approximations*]{}, Comm. Partial Differential Equations [**27**]{} (2002), 1229–1254.

H. Kalisch, D. Mitrovic, J.M. Nordbotten, [*Rayleigh-Taylor instability of immiscible fluids in porous media*]{}, Continuum Mech. Thermodyn. [**28**]{} (2016), 721–731.

G. Karch, [*Asymptotic behaviour of solutions to some pesudoparabolic equations*]{}, Math. Methods Appl. Sci. [**20**]{}, (1997) 271–289.

K. H. Karlsen, N. H. Risebro and J. Towers, [*$L^1$-stability for entropy solutions of nonlinear degenerate parabolic connection-diffusion equations with discontinuous coefficients*]{}, Skr. K. Nor. Vid. Selsk, [**3**]{} (2003), 1–49.

E.I. Kaikina, P.I. Naumkin, I.A. Shishmarev, [*The Cauchy problem for a Sobolev type equation with power like nonlinearity*]{}, Izv. Math. [**69**]{} (2005), 59–111.

C. I. Kondo, P. G. LeFloch, [*Zero diffusion-disspersion limits for scalar conservation laws*]{}, SIAM J. Math. Anal., [**33**]{} (2002), 1320–1329.

C. Kondo and P.G. LeFloch, [*Measure-valued solutions and well-posedness of multi- dimensional conservation laws in a bounded domain*]{}, Portugal. Math. [**58**]{} (2001), 171–194.

S.Khomrutai, [*Global well-posedness and grow-up rate of solutions for a sublinear pseudoparabolic equation*]{}, J. Differential Equations [**260**]{} (2016), 3598–3657.

S. N. Kruzhkov, [*First order quasilinear equations in several independent variables*]{}, Mat. Sb. [**81**]{} (1970), 217-243.

K.H. Kwek, C.C. Qu, [*Alternative principle for pseudo-parabolic equations*]{}, Dynam. Systems Appl. [**5**]{} (1996) 211–217.

M. Lazar, D. Mitrović, [*Velocity averaging – a general framework*]{}, Dynamics of PDEs, [**9**]{} (2012), 239–260.

M. Lazar, D. Mitrović, [*Existence of solutions to a scalar conservation laws with a flux of low regularity*]{}, Electronic Journal of Differential Equations, Vol. 2016 (2016), No. 325, pp. 1–18.

P.-L. Lions, B. Perthame, E. Tadmor, [*A kinetic formulation of multidimensional scalar conservation laws and related equations*]{}, J. Amer. Math. Soc. [**7**]{} (1994), 169–191.

P.G.LeFloch, [*Hyperbolic systems of conservation laws. The theory of classical and nonclassical shock waves*]{}, Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, 2002. x+294 pp. ISBN: 3-7643-6687-7.

A.Mikelić, [*A global existence result for the equations describing unsaturated flow in porous media with dynamic capillary pressure*]{}, J. Differential Equations [**248**]{} (2010), 1561–1577

A. Novick-Cohen, R.L. Pego, [*Stable patterns in a viscous diffusion equation*]{}, Trans. Amer. Math. Soc. [**324**]{} (1991) 331–351.

O. A. Oleinik, [*Discontinuous solutions of nonlinear partial differential equations*]{}, Dokl. Akad. Nauk SSSR, [**95**]{} (1954), 451–455.

V. Padron, [*Effect of aggregation on population recovery modeled by a forward-backward pseudoparabolic equation*]{}, Trans. Amer. Math. Soc. [**356**]{} (2004), 2739–2756.

B. Perthame, P. Souganidis, [*A limiting case for velocity averaging*]{}, Ann. Sci. Ec. Norm. Sup. [**4**]{} (1998), 591–598.

E. Yu. Panov, [*Existence and strong pre-compactness properties for entropy solutions of a first-order quasilinear equation with discontinuous flux*]{}, Arch. Ration. Mech. Anal. [**195**]{} (2010), 643–673.

M. Ptashnyk, [*Degenerate quasilinear pseudoparabolic equations with memory terms and variational inequalities*]{}, Nonlinear Anal. TMA [**66**]{} (2007), 2653–2675.

L.  A.  Richards, [*Capillary conduction of liquids through porous mediums*]{}, Physics N.Y. [**1**]{} (1931), 318–333.

M. E. Schonbek, [*Convergence of solutions to nonlinear dispersive equations*]{}, Comm. Partial Differential Equations [**7**]{} (1982), 959–1000.

S. Selberg, A. Tesfahun, [*Low regularity well-posedness for some nonlinear Dirac equations in one space dimension*]{}, Differential Integral Equations 23 (2010), no. 3-4, 265–278.

R.E. Showalter, T.W. Ting, [*Pseudoparabolic partial differential equations*]{}, SIAM J. Math. Anal. 1 (1970) 1–26.

A. Szepessy, [*An existence result for scalar conservation laws using measure-valued solutions*]{}, Comm. Partial Differential Equations, [**14**]{} (1989), 1329–1350.

T. Tao, E. Tadmor, [*Velocity Averaging, Kinetic Formulations, and Regularizing Effects in Quasi-Linear Partial Differential Equations*]{}, Comm. Pure Appl. Math. [**60**]{} (2007), 1488–1521.

Taylor, M. E., [*Partial Differential Equations I. Basic Theory*]{} Second edition. Applied Mathematical Sciences, 115. Springer, New York, 2011.

T.W. Ting, [*Certain non-steady flows of second-order fluids*]{}, Arch. Ration. Mech. Anal. [**14**]{} (1963), 1–26.

T.W. Ting, [*Parabolic and pseudo-parabolic partial differential equations*]{}, J. Math. Soc. Japan [**21**]{} (1969), 440–453.
