\documentclass[JEP,XML,SOM,Unicode,NoEqCountersInSection]{cedram}
\datereceived{2018-07-12}
\dateaccepted{2019-06-06}
\dateepreuves{2017-06-18}

\usepackage{cases}
\usepackage{mathrsfs}
\let\mathcal\mathscr
\let\cal\mathcal
\makeatletter
\renewenvironment{numcases}[1]%
{$$\numc@opts \setbox\z@\hbox
 {\advance\c@equation\@ne\def\@currentlabel{\p@equation\theequation}% local
 $\displaystyle {#1\null}\m@th$}%
 \numc@setsub 
 \setbox\tw@\vbox\bgroup
 \stepcounter{equation}\def\@currentlabel{\p@equation\theequation}% 
 \global\@eqnswtrue\m@th \everycr{}\tabskip\numc@left\let\\\@eqncr
 \halign to\dimen@ii \bgroup \kern\wd\z@ \kern10
 %14
 \p@ % assume width of brace
 \tabskip\z@skip \global\@eqcnt\@ne $\displaystyle{##}$\hfil 
 &\global\@eqcnt\tw@ \skip@10\p@ \advance\skip@\tw@\arraycolsep \hskip\skip@
 ##\unskip\hfil\tabskip\@centering% \unskip removes space if no explanations
 &\global\@eqcnt\thr@@\hbox to\z@\bgroup\hss##\egroup\tabskip\z@skip\cr
}{\@@eqncr \egroup % end \halign, which does not contain first column or brace
 \global\advance\c@equation\m@ne 
%Measure the natural width of the alignment
 \unskip\unpenalty\unskip\unpenalty \setbox\z@\lastbox % grab last line
 \nointerlineskip \copy\z@ % then put it back
 \global\dimen@i\wd\z@ 
 \setbox\z@\hbox{\hskip-\numc@left\unhbox\z@}% Measure its natural width
 \ifdim \wd\z@<\dimen@i \global\dimen@i\wd\z@ \fi
\egroup% end \vbox (box\tw@, box\z@ is restored to LHS)
\hbox to\dimen@ii{\m@th % assemble the whole equation
 \hskip\numc@left 
 \hbox to\dimen@i{$\displaystyle \box\z@ % parameter #1
 \dimen@\ht\tw@ \advance\dimen@\dp\tw@ % get size of brace
 \left\{\vcenter to\dimen@{\vfil}\right.\n@space % make brace
 $\hfil}\hskip\@centering % finished first part (filled whole line)
 \kern-\dimen@ii % backspace the full width
 $\vcenter{\box\tw@}$% overlay the alignment
 }% end the \hbox to\dimen@ii
\numc@resetsub
$$\global\@ignoretrue}

\newcommand{\rightharpoonupfill}{$\m@th \mathord-\mkern-6mu \cleaders\hbox{$\mkern-2mu \mathord- \mkern-2mu$}\hfill\mkern-6mu \mathord \rightharpoonup$}
\makeatother

\newdimen\lengtharrow
\newbox\exponentbox

\newcommand{\Rightharpoonupts}[1]%
{\setbox\exponentbox=\hbox{$#1$}
\lengtharrow=\wd\exponentbox
\ifdim \lengtharrow=0pt
		\mathrel{\smash{\mathop{\hbox to 6mm{\rightharpoonupfill}}}}
\else \ifdim\lengtharrow<6mm
\lengtharrow=6mm
\else \advance\lengtharrow by 1mm
\fi
		\overset{\raisebox{3pt}{$\scriptstyle#1$}}{\mathrel{\smash{\hbox to\lengtharrow{\rightharpoonupfill}}}}
	\fi}
	
\newcommand{\Rightharpoonupds}[1]%
{\setbox\exponentbox=\hbox{$#1$}
\lengtharrow=\wd\exponentbox
\ifdim \lengtharrow=0pt
		\mathrel{\smash{\mathop{\hbox to 6mm{\rightharpoonupfill}}}}
\else \ifdim\lengtharrow<6mm
\lengtharrow=6mm
\else \advance\lengtharrow by 1mm
\fi
		\overset{\raisebox{5pt}{$\scriptstyle#1$}}{\mathrel{\smash{\hbox to\lengtharrow{\rightharpoonupfill}}}}
	\fi}
	
\def\Rightharpoonup#1{\mathchoice{\Rightharpoonupds{#1}}{\Rightharpoonupts{#1}}{\Rightharpoonupts{#1}}{\Rightharpoonupts{#1}}}

\newcommand{\longrightharpoonup}{\Rightharpoonup{}}



\multlinegap0pt
\newenvironment{enumeratei}
{\bgroup\def\theenumi{\roman{enumi}}\def\theenumii{\arabic{enumii}}\begin{enumerate}}
{\end{enumerate}\egroup}
\newcommand{\RedefinitSymbole}[1]{%
\expandafter\let\csname old\string#1\endcsname=#1
\let#1=\relax
\newcommand{#1}{\csname old\string#1\endcsname\,}%
}
\RedefinitSymbole{\forall} \RedefinitSymbole{\exists}

\newcommand{\<}{\bigl\langle}
\renewcommand{\>}{\bigr\rangle}
\newcommand{\Psfrac}[2]{(\sfrac{#1}{#2})}
\newcommand{\psfrac}[2]{\sfrac{(#1)}{#2}}
\newcommand{\spfrac}[2]{\sfrac{#1}{(#2)}}
\newcommand{\pspfrac}[2]{\sfrac{(#1)}{(#2)}}

\DeclareMathOperator{\meas}{meas}
\DeclareMathOperator{\spt}{spt}

\newcommand{\intTTT}{\int_0^T\hspace*{-2.5mm}\int_{\TT}}

\newcommand{\labeleq}[1]{}
\newenvironment{subnumcasess}{\begin{equation*}\begin{cases}}{\end{cases}\end{equation*}}

\newtheorem{thm}{Theorem}[section]
\newtheorem{prop}[thm]{Proposition}
\newtheorem{lem}[thm]{Lemma}

\theoremstyle{definition}
\newtheorem{definition}[thm]{Definition}
\newtheorem{assumption}[thm]{Assumption}

\theoremstyle{remark}
\newtheorem{remark}{Remark}

\newcommand{\TT}{\mathbb{T}^3}

\DeclareMathOperator{\curl}{curl}
\DeclareMathOperator{\Div}{div}
\newcommand{\R}{\mathcal{R}}
\newcommand{\di}{\,\mathrm{d}}
\DeclareMathOperator{\D}{D}
\newcommand{\Dt}{\dfrac{\mathrm{d}}{\mathrm{d} t}}
\let\ov\overline

\let\ep\varepsilon
\newcommand{\phil}{\phi_\lambda}
\newcommand{\phip}{\phi_p}

\begin{document}
\frontmatter
\title{Compression effects in heterogeneous media}

\author[\initial{D.} \lastname{Bresch}]{\firstname{Didier} \lastname{Bresch}}
\address{LAMA UMR 5127 CNRS, Univ. Savoie Mont Blanc\\
Chambéry, France}
\email{didier.bresch@univ-smb.fr}

\author[\initial{S.} \lastname{Ne\v casová}]{\firstname{\v{S}árka} \lastname{Ne\v casová}}
\address{Institute of Mathematics, Academy of Sciences of the Czech Republic\\
Žitná 25,
CZ-115 67 Praha 1,
Czech Republic }
\email{matus@math.cas.cz}
\urladdr{http://www.math.cas.cz/homepage/main_page.php?id_membre=22}

\author[\initial{C.} \lastname{Perrin}]{\firstname{Charlotte} \lastname{Perrin}}
\address{Aix Marseille Univ., CNRS, Centrale Marseille, I2M\\ Marseille, France}
\email{charlotte.perrin@univ-amu.fr}
\urladdr{https://www.chperrin.fr/}

\thanks{Part of this work was carried out when D.\,Bresch was invited in Prague by the GA\v CR project P201-16-03230S and RVO 67985840.
D.\,Bresch is supported by the Fraise project, grant ANR-16-CE06-0011 of the French National Research Agency (ANR).
\v S.\,Ne\v casová is supported by GA\v CR project P201-16-03230S and RVO 67985840.
C.\,Perrin is supported by a PEPS project \og Jeunes chercheuses et jeunes chercheurs\fg.
This work is also supported by GACR project 19-04243S and by the SingFlows project, grant ANR-18-CE40-0027 of the French National Research Agency (ANR)}

\begin{abstract}
We study in this paper compression effects in heterogeneous media with maximal packing constraint. Starting from compressible Brinkman equations, where maximal packing is encoded in a singular pressure and a singular bulk viscosity, we show that the global weak solutions converge (up to a subsequence) to global weak solutions of the two-phase compressible/incompressible Brinkman equations with respect to a parameter $\varepsilon$ which measures effects close to the maximal packing value. Depending on the importance of the bulk viscosity with respect to the pressure in the dense regimes, memory effects are activated or not at the limit in the congested (incompressible) domain.
\end{abstract}

\keywords{Compressible Brinkman equations, maximal packing, singular limit, free boundary problem, memory effect}

\subjclass{35Q35, 35B25, 76T20}

\altkeywords{Équations de Brinkman compressibles, contrainte
d'entassement maximal, limite singulière, problème à frontière libre,
effet mémoire}

\alttitle{Effets de compression en milieux hétérogènes}

\begin{altabstract}
Nous étudions dans cet article des effets de compression dans des milieux hétérogènes soumis à une contrainte d'entassement maximal. Partant des équations de Brinkman compressibles où la contrainte maximale est prise en compte au sein d'une pression et d'une viscosité de volume toutes deux singulières, nous montrons que les solutions faibles globales convergent (à une sous-suite près) vers des solutions faibles globales d'un modèle biphasique de type compressible/incompressible quand le paramètre $\varepsilon$, mesurant l'intensité de la résistance à la compression au voisinage de l'entassement maximal, tend vers $0$. En fonction de la prédominance relative de la viscosité de volume par rapport à la pression dans les régimes denses, nous mettons en évidence l'activation d'effets de mémoire à la limite dans le domaine congestionné (incompressible).
\end{altabstract}
\maketitle
\vspace*{-\baselineskip}\enlargethispage{\baselineskip}
\tableofcontents
\mainmatter

\section*{Introduction}
We analyze in this paper macroscopic models for heterogeneous media like mixtures, suspensions or crowds, in dense regimes. These regimes exhibit interesting behaviors such as transition phases with congestion (also called jamming for granular flows) and non-local (in time and/or in space) effects which are both due to a physical packing constraint, that is the finite size of the microscopic components. At the macroscopic scale this packing constraint corresponds to a maximal density constraint $ \rho \leq \rho^*$. A very challenging issue in physics and mathematics is then to model and analyze the change of behavior in congested domains $\rho = \rho^*$ and close to a transition phase $\rho^*-\ep<\rho<\rho^*$.\enlargethispage{.5\baselineskip}%

Two different approaches are generally considered in the literature to model congestion phenomena at the macroscopic level.
The first one, usually called \emph{hard approach}, consists in coupling compressible dynamics in the free domain $\{\rho < \rho^*\}$, with incompressible dynamics in the congested domain $\{\rho = \rho^*\}$.
Associated to the incompressibility constraint on the velocity field, an additional potential (seen as the Lagrange multiplier) is activated in the congested regions.
The second one which, by opposition, is called \emph{soft approach}, prevents the apparition of congested phases by introducing in the compressible dynamics repulsive forces which become singular as $\rho$ approaches $\rho^*$.
These repulsive effects can be describe either in the pressure (constraint on the fluid at equilibrium) or in the bulk viscosity coefficient, which represents the resistance of the material to a compression.
The interested reader is referred on these two approaches to~\cite{maury2012} and Section~\ref{sec:literature} below for additional references.
An intuitive link can be made between the two approaches: if the scope of action of the repulsive forces tends to~$0$, one expects that the soft congestion model degenerates towards a hard congestion model.
We give in the Section~\ref{sec:literature} below some conjectures on this singular limit and recent results that have been obtained in this direction.
In particular, one interesting conjecture made initially by {Lefebvre-Lepot} and {Maury} in~\cite{lefebvre2011} is that a singular bulk viscosity would degenerate in the singular limit towards a (incompressible) pressure and would activate memory effects in the limit congested domain.

We want to investigate rigorously the link between soft and hard systems, by showing how the choice of the constitutive laws, the pressure and the bulk viscosity as functions of the density in the soft models, impacts the behavior of the limit hard system in congested regions assuming a constant shear viscosity.
More precisely, the main objective of this paper is to characterize the respective effects of singular pressure and bulk viscosity close to the maximal density constraint in order to understand when memory and pressure effects are activated on the limit hard congestion system.
To~that end, we consider the following three-dimensional soft congestion system (based on compressible Brinkman equations) in $(0,T)\times \TT$:\vspace*{-3pt}
\begin{subnumcases}{\label{eq:semi-stat-0}}
\partial_t \rho_\varepsilon + \Div (\rho_\varepsilon u_\varepsilon) = 0,
\label{eq:semi-stat-0-mass}\\
\nabla p_\ep(\rho_\varepsilon) - \nabla(\lambda_\ep(\rho_\varepsilon) \Div u_\varepsilon) - 2\Div(\mu \D(u_\varepsilon)) + r u_\varepsilon
= f,
\label{eq:semi-stat-0-mom}
\end{subnumcases}
where $\rho_\varepsilon$ is the density, satisfying the constraint\vspace*{-3pt}
\begin{equation}\label{constraints}
0 \le \rho_\varepsilon < 1 \quad\hbox{a.e. } (t,x) \in [0,T]\times \TT,
\end{equation}
and $u_\ep$ is the velocity vector field in the material.
The coefficients $p_\ep$ and $\lambda_\ep$ are respectively the pressure law and the bulk viscosity coefficient, defined in this paper~as
\begin{equation} \label{eq:singular-law}
p_\ep(\rho_\varepsilon) = \ep \Bigl(\dfrac{\rho_\varepsilon}{1-\rho_\varepsilon}\Bigr)^\gamma, \quad
\lambda_\ep(\rho_\varepsilon) = \ep \Bigl(\dfrac{\rho_\varepsilon}{1-\rho_\varepsilon}\Bigr)^\beta \quad \text{with } \gamma, \beta > 1,
\end{equation}
while the shear viscosity is assumed to be constant: $\mu>0$.
Finally, $ru_\ep$ with $r>0$ represents the drag and the right-hand term, $f$, is a given external force.
Initially $\rho_\varepsilon\vert_{t=0} = \rho_0^\varepsilon$ with
\begin{equation}\label{hyp:bound-rho0}
\begin{gathered}
0 \le \rho^0_\varepsilon \le 1- R_\varepsilon < 1 \hbox{ and } R_\varepsilon \to 0
\hbox{ when } \varepsilon \to 0, \\
\frac{1}{|\TT|}\int_{\TT} \rho^0_\varepsilon (x)\, dx\le M^0 <1.
\end{gathered}
\end{equation}
Let us encode the effect of the singular bulk viscosity through the following PDE equation
that may be obtained from the mass equation
\[
\partial_t (\Lambda_\varepsilon (\rho_\varepsilon))
+ \Div (\Lambda_\varepsilon(\rho_\varepsilon) u_\varepsilon)
= - \lambda_\varepsilon(\rho_\varepsilon) \Div u_\varepsilon,
\]
where
\begin{equation}
\label{constraint1}
\Lambda_\varepsilon(\rho_\varepsilon) = \rho_\varepsilon \int_0^{\rho_\varepsilon} \lambda_\varepsilon (\tau)/\tau^2 \, d\tau =
\rho_\varepsilon \Bigl[\frac{1}{\beta-1} \, \varepsilon^{(1+\gamma-\beta)/\gamma}(p_\varepsilon(\rho_\varepsilon))^{(\beta-1)/\gamma}\Bigr].
\end{equation}

The main objective now is to understand the asymptotic regime which may be obtained by letting $\varepsilon$ go to zero. This corresponds to the limit towards the hard approach explained previously. Let us assume that $(\rho_\varepsilon, u_\varepsilon, p_\varepsilon(\rho_\varepsilon), \Lambda_\varepsilon(\rho_\varepsilon))$ tends to $(\rho, u, p, \Lambda)$. Then we get the following system in $(0,T)\times \TT$:
\begin{subnumcases}{\label{LimitSystem1}}
\partial_t \rho + \Div (\rho u) = 0, \\
\nabla p - \nabla \Pi - 2\Div(\mu \D(u)) + ru = f, \\
0\le \rho \le 1 \hbox{ and } p \ge 0,
\end{subnumcases}
where
\begin{equation}\label{LimitSystem2}
\Pi = - (\partial_t \Lambda+ \Div (\Lambda u)) \quad\hbox{with } \Lambda \ge 0.
\end{equation}
We also get the following limit initial data
\begin{equation} \label{inicongestion}
\rho\vert_{t=0} = \rho^0 \in [0,1], \qquad \Lambda\vert_{t=0} = \Lambda^0
\qquad \hbox{ in }~ \TT.
\end{equation}
It remains now to close the limit system by deriving two constraints. One of these constraints will result from Equality \eqref{constraint1} depending on the sign of $1+\gamma - \beta$ appearing explicitly in the power of $\ep$. For the last constraint, different scenarios will be obtained using one of the two following relations
\begin{equation}\label{constraint2}
(1-\rho_\varepsilon) \, p_\varepsilon(\rho_\varepsilon) =\varepsilon^{1/\gamma} \rho_\varepsilon (p_\varepsilon(\rho_\varepsilon))^{(\gamma-1)/\gamma},
\end{equation}
or
\begin{equation} \label{constraint3}
(1-\rho_\varepsilon) \, \Lambda_\varepsilon(\rho_\varepsilon)
= c(\beta) \varepsilon^{1/(\beta-1)} \rho_\varepsilon^{\beta/(\beta-1)} (\Lambda_\varepsilon(\rho_\varepsilon))^{(\beta-2)/(\beta-1)}.
\end{equation}

More precisely, passing to the limit in \eqref{constraint1} and \eqref{constraint2}--\eqref{constraint3}, we find the following relations in addition to the system \eqref{LimitSystem1}--\eqref{inicongestion}:
\begin{itemize}
\item If $1+\gamma -\beta=0$ (memory and pressure effect):
\begin{equation}\label{MPE}
\rho p = (\beta-1) \Lambda \qquad \hbox{ and } \qquad (1-\rho) \, \Lambda=0.
\end{equation}
\item If $1+\gamma - \beta <0$ (memory but no pressure effect):
\begin{equation}\label{MPPE}
p = 0 \qquad \hbox{ and } \qquad
(1-\rho)\, \Lambda= 0.
\end{equation}
\item If $1+\gamma -\beta>0$ (pressure but no memory effect):
\begin{equation} \label{PPME}
\Lambda=0
\qquad \hbox{ and } \qquad (1-\rho) \, p = 0.
\end{equation}
\end{itemize}

Observe that this formal analysis could be generalized to more general pressure and bulk viscosity laws than~\eqref{eq:singular-law}, to take into account different (singular) possible behaviors close to the maximal constraint.
The key argument relies here in the comparison between the pressure $p_\ep$ and the coefficient $\Lambda_\ep$ in the vicinity of the maximum density.
Let us emphasize the fact that there is no consensus in physics around the order of singularity of these laws (see for instance~\cite{andreotti2013} or~\cite{coussot2005}).

Note that it is well known that the compressibility of a fluid may be encoded in the pressure and in the bulk viscosity.
Indeed, incompressible systems may be obtained by letting the Mach number $\mathrm{Ma}$, which appears in the dimensionless Navier-Stokes equations in front of the pressure $\Psfrac{1}{\mathrm{Ma^2}}\nabla p(\rho)$, go
to zero (see for instance the works of {Desjardins} et al.~\cite{desjardins1999}, {Lions, Masmoudi}~\cite{lions1998Mach}, {Feireisl, Novotn{\'y}}~\cite{feireisl2009}).
But the incompressible equations can be also obtain from a large bulk viscosity limit:
if in the bulk viscosity term $\nabla(\lambda^0\Div u)$ one lets $\lambda^0$ go to $+\infty$ then, formally, $\Div u$ should tend to $0$.
This result has been recently proved by {Danchin} and {Mucha} in~\cite{danchin2017}.

In our paper, the main novelty is to consider both singular pressure and singular bulk viscosity depending on the density which will encode incompressibility of the material at the maximal packing value $\rho^*=1$, assuming the shear viscosity to be constant.
Below this maximal packing value, the material remains compressible.
It is also interesting from the physical point of view to consider density dependent shear viscosities $\mu(\rho)$, this case has been treated in the two papers~\cite{perrin2016,perrin2017}.
The results are presented in the section state of the art (subsection II-ii).
It has to be noted that the mathematical tools and difficulties are in that case completely different from those used in the present paper.
Historically, studies on compressible Navier-Stokes system with (non-singular) density dependent bulk viscosity $\lambda(\rho)$ and constant shear viscosity $\mu >0$, start from the beautiful paper~\cite{vaigant1995} by {Kazhikov} and {Waigant} where they prove global existence of strong solutions in two dimensions with periodic boundary conditions and with no vacuum state if initially no vacuum exists. In their paper, the pressure is assumed $p(\rho) = a \rho^\gamma$, $\mu > 0$ and $\lambda(\rho)= \rho^\beta$ with $\beta >3$.
Following this result, {Perepelitsa} proved in~\cite{perepelitsa2006} the global existence of a weak solution with uniform lower and upper bounds on the density when the initial density is away from vacuum.
Finally, the hypothesis on the coefficient $\beta$ has been recently relaxed with possible vacuum state in~\cite{huli2016} and bounded domains have been considered in~\cite{ducomet2013}.
It would be interesting to investigate the problem for singular bulk viscosity and singular pressure laws for the 3D compressible Navier-Stokes equations but this is not the main objective of our paper. We focus here on 3D Brinkman equations where the total acceleration of the fluid is neglected. A typical application we have in mind is the modeling of flows in porous media.
Brinkman equations are a classical extension of the Darcy equation:
\[
u = -\nabla p + f,
\]
with additional viscous terms, here $\nabla(\lambda\Div u) + 2\Div(\mu \D(u))$.
In their incompressible version, these equations have been rigorously derived by {Allaire} in~\cite{allaire1991} by homogenization techniques from Navier-Stokes equations in a perforated domain.
His result has been then extended by {Desvillettes} et al.\ \cite{desvillettes2008} and {Mecherbet}, {Hillairet}~\cite{mecherbet2018}.
The recent study~\cite{Nasser2017} provides some new analysis and numerical results on these equations in the incompressible case.
These equations in the compressible setting may also apply in biology in tissue growth modeling or in petroleum problem occurring in compressible porous matrix. The interested reader is referred to the study of {Perthame}, {Vauchelet}~\cite{perthame2015}, {Nasser El Dine} et al.~\cite{Nasser2017}, \cite{nasser2019},
or {Énault}~\cite{Enault10} and the references therein.\\
From a mathematical point of view, of course, it could be interesting to study the compressible Navier-Stokes equations with singular pressure and bulk viscosity.
Both the estimates on the effective flux and the compactness arguments are then of course much more tricky to handle due to the additional acceleration term.
In view of the applications we have in mind, this case is beyond the scope of our paper.

The paper will be organized as follows: We will first present the main existence and convergence results,
then we will review mathematical studies that have been realized recently around the subject of congestion problems. In the second section, we present important mathematical
properties linked to the system under consideration and the truncated system we first study.
Passing to the limit with respect to the parameter of the truncation, $\delta$, we get the global existence of weak solutions for the original system~\eqref{eq:semi-stat-0} at $\ep$ fixed.
It will be then possible to pass to the limit with respect to $\varepsilon$ to recover the hard congestion system \eqref{LimitSystem1}--\eqref{LimitSystem2}
with two additional relations which will be, depending on the parameters $\gamma$ and $\beta$, given by \eqref{MPE} or \eqref{MPPE} or
\eqref{PPME}.
We will divide the study in two sections depending on the sign of $\gamma-\beta$ which correspond to the dominant pressure regime $\gamma > \beta$ (Section~\ref{sec:pressure}) or dominant bulk regime $\beta \ge \gamma$ (Section~\ref{sec:bulk}).

\section{Main results}
We first prove in the paper the existence of global weak solutions to the soft congestion system \eqref{eq:semi-stat-0} when the pressure and the bulk viscosity are defined by~\eqref{eq:singular-law}. For simplicity, we assume in addition that
\begin{equation}\label{hyp:f}
f \in L^2\big(0,T; (L^q(\TT))^3\big) \quad\hbox{with } q>3.
\end{equation}
Let us first give our definition of global weak solutions of \eqref{eq:semi-stat-0}--\eqref{hyp:bound-rho0} and \eqref{LimitSystem1}--\eqref{inicongestion}.
\begin{definition}[Weak solutions of the soft congestion system]
A pair $(\rho_\ep,u_\ep)$ is called a global (bounded energy renormalized) weak solution to \eqref{eq:semi-stat-0}--\eqref{hyp:bound-rho0} if for any $T > 0$, the following properties hold.
\begin{itemize}
\item $\rho_\ep \in \mathcal{C}([0,T];L^q(\TT)) \cap L^\infty((0,T) \times \TT)$ for all $q \in [1, +\infty)$;
\item
$u_\ep\in L^2(0,T;(H^1(\TT))^3)$;
\item $0 \leq \rho_\ep(t,x) < 1$ a.e. in $(0,T)\times \TT$;
\item $(\rho_\ep,u_\ep)$ satisfies \eqref{eq:semi-stat-0-mass}--\eqref{eq:semi-stat-0-mom}
in the weak sense:
\begin{multline*}\labeleq{weak-ep-mass}
\intTTT{\rho_{\ep} \partial_t \phi}
+ \intTTT{\rho_{\ep} u_{\ep} \cdot \nabla \phi} \\
= \int_{\TT}{\rho_\ep(T) \phi(T)} - \int_{\TT}{\rho_\ep^0 \phi(0)}
\qquad \forall \phi \in \mathcal{C}^1([0,T] \times \TT);
\\
\labeleq{weak-ep-div-mom}
\shoveleft{- \intTTT{p_{\ep}(\rho_{\ep})\Div \psi}
+ \intTTT{(2\mu+\lambda_{\ep}(\rho_{\ep}))
\Div u_{\ep} \, \Div \psi}} \\
+ \intTTT
\mu\curl u_{\ep} \cdot \curl \psi
+ r\intTTT{u_{\ep} \cdot \psi}
= \int_{\TT}{f \cdot \psi}
\quad \forall \psi \in \mathcal{C}^1([0,T] \times \TT).
\end{multline*} 
\item The \emph{renormalized continuity equation} holds
\begin{equation*}\labeleq{renormalized-ep}
\partial_t b(\rho_{\ep}) + \Div(b(\rho_{\ep}) u_{\ep})
+ \big(b_+'(\rho_{\ep})\rho_{\ep}
- b(\rho_{\ep})\big)\Div u_{\ep} = 0
\quad \text{in}~~\mathcal{D}'((0,T)\times \TT),
\end{equation*}
for any $b\in \mathcal{C}([0,+\infty))$, piecewise ${\cal C}^1$, where $b'_+$
denotes the right derivative of $b$.

\item The energy inequality holds
\begin{multline*}\labeleq{energy_ep-0}
\sup_{t\in[0,T]} \int_{\TT} H_{\ep}(\rho_{\ep})
+ \intTTT{(\tfrac{3}{2}\mu+\lambda_{\ep}(\rho_{\ep}))(\Div u_{\ep})^2 } + \frac{\mu}{2} \intTTT{ |{\curl (u_{\ep})}|^2}\\
+ r \intTTT{ |u_{\ep}|^2}
\leq \int_{\TT}{H_{\ep}(\rho_\ep^0)} + \frac{1+C}{2\mu}\, \|f\|^2_{L^2_{t,x}}.
\end{multline*}
with $C$ the constant linked to Poincaré-Wirtinger inequality
$$\|g - \<g\>\|^2_{L^2(\TT)} \le C \|\nabla g \|^2_{L^2(\TT)}, $$
and where
\begin{equation*}\labeleq{H_ep}
H_{\ep}(\rho) = \dfrac{\ep}{\gamma-1}\cdot \dfrac{\rho^{\gamma}}{(1-\rho)^{\gamma-1}}.
\end{equation*} 
\end{itemize}
\end{definition}

\begin{definition}[Weak solutions of the hard congestion system]
We say that $(\rho,u,p,\Lambda)$ is a global weak solution to \eqref{LimitSystem1}--\eqref{inicongestion} if for any $T > 0$, it satisfies
\begin{itemize}
\item the following regularities
\begin{align*}
\rho &\in \mathcal{C}([0,T];L^q(\TT)) \cap L^\infty((0,T) \times \TT)\quad \text{for all } q \in [1, +\infty)\\
u &\in L^2(0,T;(H^1(\TT))^3),\\
p &\in \mathcal{M}_+((0,T)\times \TT) \quad \text{and} \quad \Lambda \in L^\infty(0,T;L^2(\TT));
\end{align*}
\item $0 \leq \rho(t,x) \leq 1$ a.e. in $(0,T)\times \TT$;
\item $(\rho,u,p,\Lambda)$ satisfies equations \eqref{LimitSystem1}--\eqref{LimitSystem2} in the sense of distributions.
\end{itemize}
\end{definition}

We prove in this paper the following existence results

\pagebreak[2]
\begin{thm}{\label{thm-ep}}
Let $\rho_\ep^0$ satisfying condition \eqref{hyp:bound-rho0}.
Assume in addition that
\begin{itemize}
\item if $1< \beta< \gamma$
\begin{equation}\label{hyp:energy-t0-C1}
\int_{\TT}{\dfrac{\ep}{(1-\rho^0_\ep)^{\gamma-1}}}\leq E^0 <+\infty;
\end{equation}
\item if $1 < \gamma \leq \beta$
\begin{equation}\label{hyp:energy-t0-C2}
\int_{\TT}{\big(\Lambda_\ep(\rho_\ep^0)\big)^2}
\leq \Lambda^0 < +\infty.
\end{equation}
\end{itemize} 
Then there exist global (bounded energy renormalized) weak solutions of
\eqref{eq:semi-stat-0}--\eqref{hyp:bound-rho0}
at $\ep$ fixed.
\end{thm}

The next result justifies the formal derivation of system \eqref{LimitSystem1}--\eqref{inicongestion} respectively with relations \eqref{MPE} if $\gamma= \beta-1$,
with relations \eqref{MPPE} if $\gamma<\beta-1$, relations
\eqref{PPME} if $\gamma>\beta-1$.
More precisely

\begin{thm}{\label{thm-limit}} As $\ep \rightarrow 0$, there exists a subsequence $(\rho_\ep,u_\ep)$ of global weak solutions of \eqref{eq:semi-stat-0}--\eqref{hyp:bound-rho0} such that $(\rho_\varepsilon, u_\varepsilon, p(\rho_\varepsilon), \Lambda_\varepsilon(\rho_\varepsilon))$ converges weakly to $(\rho,u, p, \Lambda)$ a global weak solution of the hard congestion system \eqref{LimitSystem1}--\eqref{inicongestion} satisfying two algebraic relations which encode the competition between the singular pressure and bulk viscosity, namely
\begin{enumerate}\renewcommand{\theenumi}{\Roman{enumi}}
\item
Memory effect in the congested domain {\rm ($\gamma \le \beta -1$)}:
\begin{itemize}
\item For $\gamma < \beta - 1$ {\rm (no pressure effect)}:
$$p=0, \qquad (1-\rho)\Lambda = 0 \quad \text{a.e.}$$ 
\item For $\gamma= \beta - 1$ {\rm (pressure effect)}:
$$\rho p= (\beta-1) \Lambda, \qquad (1-\rho)\Lambda = 0 \quad \text{a.e.}$$
\end{itemize} 

\item
No memory effect in the congested domain {\rm (pressure effect)}
{\rm ($\gamma > \beta -1$)}:
$$\Lambda= 0 \quad \text{a.e.}, \qquad \spt p \subset \{\rho =1\}.$$ 
\end{enumerate}
\end{thm}

\begin{remark}
Our system can be seen as a ``semi-stationary'' version of the compressible Navier-Stokes equations with the additional friction term $ru$. If there is no friction in the equations, namely $r=0$, then in the periodic case, the velocity $u$ is defined by Equation \eqref{eq:semi-stat-0-mom} up to a constant.
Therefore, we would need to impose an additional constraint on the velocity, \eg
\[
\int_{\TT}{u \di x} = 0.
\]
Integrating in space the momentum equation \eqref{eq:semi-stat-0-mom}, we would also need the ``compatibility relation''
\[
\int_{\TT}{f \di x} = 0.
\]
Provided these additional constraints, our two results remain unchanged.
\end{remark}

\subsubsection*{Incompressible dynamics in the congested domain}
In addition to the limit equations~\eqref{LimitSystem1}--\eqref{inicongestion}, one could have add the incompressibility condition $\Div u = 0$ on the congested sets $\{\rho = 1\}$.
More precisely, we have the following lemma.
\begin{lem}[{\cite[Lem.\,2.1]{lions1999}}]
Let $u \in L^2(0,T;(H^1(\mathbb{T}))^3)$ and $\rho\in L^2((0,T)\times\mathbb{T}^3)$ such that
\[
\partial_t \rho+\Div(\rho u)=0 \quad \text{in} \ (0,T)\times\Omega,\qquad \rho(0)=\rho^0,
\]
then the following two assertions are equivalent
\begin{enumeratei}
\item
$\Div u = 0$ a.e on $\{\rho \geq 1\}$ and $0 \leq \rho^0 \leq 1$,
\item
$0 \leq \rho(t,x) \leq 1$ a.e. $(t,x)$.
\end{enumeratei}
\end{lem}

\begin{remark}
An important issue concerning the limit systems that we obtain is the regularity of the limit pressure $p$.
Through our approximation procedure, the limit pressure $p$, if it is not $0$, is a priori a non-negative measure.
If one is able to prove that $p \in L^1((0,T)\times \TT)$, it is thus possible to give a sense a.e. to the product $\rho p$ at the limit and then to the ``exclusion constraint''
\[
(1-\rho) p = 0,
\]
which is another way to express the activation of the pressure in the congested zones (see the system~\eqref{eq:BBCR} written below).
In fact, this is less the justification the exclusion constraint than the regularity of the pressure which is crucial in the mathematical understanding of partially congested flows.
\end{remark}

\section{Historical remarks}{\label{sec:literature}}
For reader's convenience, we present below the context of this study and give some historical remarks concerning limits from soft approaches to hard approaches for congestion problems.

\subsubsection*{\textup{I}\enspace Derivation from compressible Euler equations}
A first generic hard congestion model is derived in~\cite{bouchut2000} by {Bouchut} et al. from one-dimensional two-phase gas/liquid flows. The equations read
\begin{subnumcases}{\label{eq:BBCR}}
\partial_t\rho + \partial_x(\rho u) = 0, \\
\partial_t(\rho u) + \partial_x(\rho u^2) + \partial_x p = 0, \\
0 \leq \rho \leq 1, ~ (1-\rho)p= 0, ~ p \geq 0,
\end{subnumcases}
in which the constraint $(1-\rho)p = 0$, sometimes called ``exclusion constraint'', expresses the activation of the pressure $p$ in the congested phase where $\rho = 1$.
The pressure ensures that the maximal density constraint $\rho^*=1$ is not exceeded.
This system has been then studied theoretically by {Berthelin} in~\cite{berthelin2002,berthelin2017} who constructs global weak solutions by means of an approximation with sticky blocks (see~\cite{maury2017} for an associated numerical method). {Degond} et al. approximate numerically in~\cite{degond2011,degond2017} the solutions of~\eqref{eq:BBCR} with an appropriate discretization of the soft congestion system
\begin{subnumcases}{\label{eq:BBCR-soft}}
\partial_t\rho + \partial_x(\rho u) = 0, \\
\partial_t(\rho u) + \partial_x(\rho u^2) + \partial_x p_\ep(\rho) = 0, \\
p_\ep(\rho) = \ep \Bigl(\dfrac{\rho}{1-\rho}\Bigr)^\gamma, \quad \gamma > 1.
\end{subnumcases}
Although the rigorous derivation of Equations~\eqref{eq:BBCR} from~\eqref{eq:BBCR-soft} (\ie the limit $\ep \rightarrow 0$ in~\eqref{eq:BBCR-soft}) has not been proved theoretically, the authors obtain satisfactory numerical results thanks to smart treatment of the singular pressure $p_\ep$ for small $\ep$.
Let us also mention on the subject the study~\cite{bresch2017} which addresses the issue of the creation of congested zones in 1D and highlights the multi-scale nature of the problem.

\Subsubsection*{\textup{II}\enspace Derivation from compressible Navier-Stokes equations}

\subsubsection*{\textup{II\,(i)}\enspace Compressible Navier-Stokes equations with constant viscosities}
The first justification of the link between a soft congestion system and a hard congestion system is given in~\cite{bresch2014} for the one-dimensional case.
In~\cite{peza2015}, the existence of global weak solutions to the multi-dimensional viscous equations
\begin{subnumcases}{\label{eq:peza-soft}}
\partial_t \rho + \Div(\rho u) = 0, \\
\partial_t(\rho u) + \Div(\rho u \otimes u) + \nabla p_\ep(\rho) - \nabla(\lambda \Div u) - 2\Div(\mu \D(u)) = 0,\\
p_\ep(\rho) = \ep \Bigl(\dfrac{\rho}{1-\rho}\Bigr)^\gamma, \quad
\gamma > 3, \quad 2\mu + \lambda > 0,
\end{subnumcases}
is first proved for a fixed $\ep > 0$.
Then, the authors show the weak convergence of these solutions as $\ep \rightarrow 0$ toward global weak solutions of the viscous hard congestion system
\begin{subnumcases}{\label{eq:peza-hard}}
\partial_t \rho + \Div(\rho u) = 0, \\
\partial_t(\rho u) + \Div(\rho u \otimes u) + \nabla p - \nabla(\lambda \Div u) - 2\Div(\mu \D(u)) = 0,\\
0 \leq \rho \leq 1, ~ (1-\rho)p= 0, ~ p \geq 0.
\end{subnumcases}

\begin{remark}
Note that the condition $\gamma > 3$ was assumed in~\cite{peza2015} to prove the existence of global weak solutions to~\eqref{eq:peza-soft}.
Precisely, it was used to prove the equi-integrability of the approximate truncated pressure $p_{\ep,\delta}(\rho_{\ep,\delta})$ as $\delta \rightarrow 0$ (see details of the truncation process in the next Section).
It is possible in fact to improve the bound on $\gamma$ and show the existence for $\gamma > \sfrac{5}{2}$ as it has been done by {Feireisl} et al. in~\cite{feireisl2016}.
\end{remark}

\begin{remark} Originally, {Lions} and {Masmoudi} in~\cite{lions1999} have obtained the same viscous system from the compressible Navier-Stokes equations with constant viscosities and pressure $p(\rho) = a \rho^{\gamma_n}$ letting $\gamma_n \rightarrow +\infty$.
The same limit has been used by {Perthame} et al.~\cite{perthame2014} for tumor growth modeling on the basis of the porous medium equation instead of Navier-Stokes equations (see the study of {Vauchelet} and {Zatorska}~\cite{vauchelet2017} in the case of Navier-Stokes equations with additional source term in the mass equation).
In this context the singular limit leads to the Hele-Shaw equations, this problem is sometimes called in the literature ``mesa problem''.
\end{remark}

\Subsubsection*{\textup{II\,(ii)}\enspace Compressible Navier-Stokes equation with singular density dependent viscosities}
In the modeling of immersed granular flows this type of singular limit has enabled to prove in~\cite{perrin2016} the link between the suspension regime and the granular regime which was an open conjecture in physics (see~\cite{andreotti2013}).
Precisely, global weak solutions to the following suspension model
\begin{subnumcasess}{\labeleq{suspension}}
\partial_t \rho + \Div(\rho u) = 0, \\
\partial_t(\rho u) + \Div(\rho u \otimes u) + \nabla p_\ep(\rho) - \nabla(\lambda_\ep(\rho) \Div u) \nonumber \\
\hspace{4.65cm}{} - 2\Div(\mu_\ep(\rho) \D(u)) + r \rho |u|u= 0,
\end{subnumcasess}%
are proved to exist at $\ep > 0$ fixed for singular viscosities and pressure such that
\begin{equation*}
\mu_\ep(\rho) = \mu_0(\rho + p_\ep(\rho)), \quad p_\ep(\rho) = \ep \Bigl(\dfrac{\rho}{1-\rho}\Bigr)^\gamma, \quad \gamma > 1,
\end{equation*}
and $\lambda_\ep$ satisfying a specific relation with the shear viscosity (and thus with the pressure), namely
\begin{equation}\label{eq:mu_lambda}
\lambda_\ep(\rho) = 2(\rho \mu_\ep'(\rho)-\mu_\ep(\rho)).
\end{equation}
Under these hypothesis, the solutions are shown to converge to global weak solutions~of
\begin{subnumcases}{\label{eq:granular}}
\partial_t \rho + \Div(\rho u) = 0, \\
\partial_t(\rho u) + \Div(\rho u \otimes u) + \nabla p + \nabla \Pi\\
\hspace*{4.1cm} {}- 2\mu_0\Div((\rho + p) \D(u)) + r\rho|u|u = 0,\nonumber\\
\partial_t p + u \cdot \nabla p = \dfrac{\Pi}{2\mu_0}, \label{eq:memory}\\
0 \leq \rho \leq 1, ~ (1-\rho)p = 0,~ p \geq 0,
\end{subnumcases}
where the pressures $p$ and $\Pi$ are respectively the weak limits of $p_\ep(\rho_\ep)$ and $\lambda_\ep(\rho_\ep)\Div u_\ep$.
The important difference between~\eqref{eq:granular} and~\eqref{eq:peza-hard} is the activation of an additional equation~\eqref{eq:memory} linking the two pressures $p$ and $\Pi$.
It results from the relation~\eqref{eq:mu_lambda} that is imposed at $\ep$ fixed.
Indeed, the conservation of mass and~\eqref{eq:mu_lambda} yield (at least formally)
\[
\partial_t \mu_\ep(\rho_\ep) + \Div(\mu_\ep(\rho_\ep) u_\ep) = - \dfrac{1}{2}\,\lambda_\ep(\rho_\ep)\Div u_\ep,
\]
which gives at the limit~\eqref{eq:memory} due to the incompressibility constraint $\Div u = 0$ that is satisfied in the congested domain.
From a modeling point of view, Equation~\eqref{eq:memory} expresses some memory effects in the congested regions, effects that were first identified by {Lefebvre-Lepot} and {Maury} in a macroscopic 1D model for ``viscous contact''~\cite{lefebvre2011} (see also~\cite{lefebvre2009} for a microscopic approach).
From a mathematical point of view, this equation is necessary to close the system and relates $\Pi$, which can be seen as the Lagrange multiplier associated to the constraint $\Div u = 0$ in the congested domain, and $p$ called \emph{adhesion potential} which characterizes the memory effects.
This is thus the singularity of bulk viscosity $\lambda_\ep$ which is responsible for the activation of memory effects in~\eqref{eq:granular}.

In the present paper, we characterize precisely the respective effects of pressure and bulk viscosity.
At the limit on the hard congestion system, we cover in particular the two cases introduced in~\cite{peza2015} and~\cite{perrin2017} where pressure effects or memory effects are activated.

\section{Structural properties and approximate system}

This section is divided into three parts.
After introducing some important quantities, such as the effective flux, and deriving crucial properties linking the pressure and the bulk viscosity, we present an approximate truncated system which formally degenerates to our original singular system~\eqref{eq:semi-stat-0} as the cut-off parameter tends to~$0$.
The last part details how we can construct global weak solutions to the truncated system.

\subsubsection*{Structural properties, effective flux}
Let $F$ be the viscous effective flux defined as
\begin{equation*}\labeleq{df:eff_flux}
F = (2\mu + \lambda(\rho))\Div u - p(\rho),
\end{equation*}
and the function $\nu$ defined from the viscosity coefficients
\begin{equation}\label{df:nu}
\nu(\rho) = \dfrac{1}{2\mu + \lambda(\rho)}.
\end{equation}
We prove the following Lemma.
\begin{lem}
Let $(\rho,u)$ satisfying in the weak sense the equations
\begin{subnumcasess}{\labeleq{sys-0}}
\partial_t \rho + \Div (\rho u) = 0, \\
\nabla p(\rho) - \nabla(\lambda(\rho) \Div u) - 2\Div(\mu \D(u)) + ru = f,
\end{subnumcasess}%
and denote
\[
S:= (-\Delta)^{-1}\Div \big(f -r u\big),
\]
where $(-\Delta)^{-1}$ is the inverse operator of the Laplacian.
More precisely, if for a periodic function $g$ such that $\<g\> =0$, where $\<g\>:= (|\TT|)^{-1} \int_{\TT}{g(x) \di x}$, $h=(-\Delta)^{-1}g$ is the unique periodic solution of
\[
-\Delta h = g\quad \text{in}~~ \TT, \quad \<h\> = 0.
\]
Then the following relations hold
\begin{align} \label{eq:lambda_divu}
\<\lambda(\rho)\Div u\> &= \<p(\rho)\> - \dfrac{\int_{\TT}{(p(\rho)+ S)\nu(\rho)}}{\int_{\TT}{\nu(\rho)}},\\
\label{eq:F-0}
F&= S - \frac{\<(p(\rho) +S)\, \nu(\rho)\>}{\<\nu(\rho)\>}.
\end{align}
\end{lem}

\begin{proof}
Observe first that integration in space of the momentum equation yields
\begin{equation*}\labeleq{int-f-u}
r \int_{\TT}{u} = \int_{\TT}{f}.
\end{equation*}
Applying the $\Div$ operator to~\eqref{eq:semi-stat-0-mom} we obtain
\begin{equation*}
\Delta F = \Div (f-ru).
\end{equation*}
Then
\begin{equation*}
F = \<F\> + S,
\end{equation*}
\ie
\begin{equation}\label{eq:expr_lambdadivu}
(2\mu+ \lambda(\rho)) \Div u - p(\rho) = S + \< \lambda(\rho) \Div u - p(\rho) \>.
\end{equation}
Let us now characterize the mean value of the effective flux in terms of
the density: rewriting this equation as
\begin{equation*}
\Div u = \nu(\rho)\bigl(S + \<\lambda(\rho)\Div u- p(\rho) \> + p(\rho) \bigr)
\end{equation*}
and integrating in space, we arrive at~\eqref{eq:lambda_divu}.
Replacing this expression in \eqref{eq:expr_lambdadivu}, we finally get~\eqref{eq:F-0}.
\end{proof}

\subsubsection*{Approximate system}
We introduce now a cut-off parameter $\delta \leq \delta^0 \in (0,1)$ in order to truncate the singular laws $p_\varepsilon$ and $\nu_\varepsilon$.
Namely, we define the truncated laws
\begin{align}\label{eq:p_delta}
p_{\ep,\delta}(\rho) &=
\begin{cases}
\ep \sfrac{\rho^\gamma}{(1-\rho)^\gamma}
& \text{if } \rho \leq 1-\delta, \\
\ep \sfrac{\rho^\gamma}{\delta^\gamma}
\quad & \text{if } \rho > 1-\delta,
\end{cases}
\\[5pt]
\label{eq:lambda_delta} \qquad\text{and} \qquad
\lambda_{\ep,\delta}(\rho) q&=
\begin{cases}
\ep \sfrac{\rho^\beta}{(1-\rho)^\beta} \quad & \text{if } \rho \leq 1-\delta, \\
\ep \sfrac{\rho^\beta}{\delta^\beta} \quad & \text{if } \rho > 1-\delta,
\end{cases}
\end{align}
and consider the associated system
\begin{subnumcases}{\label{eq:semi-stat-delta}}
\partial_t \rho_{\varepsilon, \delta}
+ \Div (\rho_{\varepsilon, \delta} u_{\varepsilon, \delta}) = 0,
\label{eq:system-delta-mass}\\
\nabla p_{\ep,\delta}(\rho_{\ep,\delta}) - \nabla(\lambda_{\ep,\delta}(\rho_{\ep,\delta}) \Div u_{\ep,\delta}) - 2\Div(\mu \D(u_{\ep,\delta})) + r u_{\ep,\delta},
= f \label{eq:system-delta-momentum}
\end{subnumcases}
with initial density
\begin{equation*}
\rho_{\ep,\delta}\vert_{t=0} = \rho_\ep^0.
\end{equation*}
Let us first give some properties related to $\nu_{\varepsilon,\delta} (\rho_{\varepsilon, \delta})$.
\begin{lem}{\label{lem:nu_delta}}
Assume that $\rho^0_\ep$ satisfies~\eqref{hyp:bound-rho0}.
There exist $C_1,C_2 >0$ which do not depend on $\delta$ or $\ep$ such that
\[
\|\nu_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^\infty_{t,x}}\leq C_1,
\qquad \int_{\TT}{\nu_{\ep,\delta}(\rho_{\ep,\delta})} \geq C_2.
\]
\end{lem}

\begin{proof} By definition of $\nu_{\ep,\delta}$~\eqref{df:nu}, we directly get
\[
\nu_{\ep,\delta} = \dfrac{1}{2\mu + \lambda_{\ep,\delta}} \leq \dfrac{1}{2\mu} < +\infty.
\]
Under the assumption on the initial mass~\eqref{hyp:bound-rho0}${}_2$ (which does not depend on $\delta$ or $\ep$), we have
\[
M^0 |\TT|
\geq \int_{\TT}{\rho_{\ep,\delta}} \geq \int_{\TT}{\rho_{\ep,\delta} \,\mathbf{1}_{\{\rho_{\ep,\delta} \geq \psfrac{1+M^0}{2}\}} }
\geq \dfrac{1+M^0}{2} \meas\{\rho_{\ep,\delta} \geq \psfrac{1+M^0}{2}\}.
\]
Now, since $ M^0 < \psfrac{1+M^0}{2 }$, it follows that
\begin{equation*}
\meas \left\{ \rho_{\ep,\delta} \geq \psfrac{1+M^0}{2} \right\} \leq \dfrac{2M^0}{1+M^0}\,|\TT| < |\TT|,
\end{equation*}
and
\begin{equation*}
\meas \left\{ \rho_{\ep,\delta} < \psfrac{1+M^0}{2} \right\} \geq K > 0,
\end{equation*}
with $K = 1- \spfrac{2M^0}{1+M^0}$ independent of $\delta$ and $\ep$.
Then,
\begin{align*}
\int_{\TT}{\nu_{\ep,\delta}(\rho_{\ep,\delta})}
& = \int_{\TT}{\dfrac{1}{2\mu + \lambda_{\ep,\delta}(\rho_{\ep,\delta})} }\\
& \geq \int_{\TT}{\dfrac{1}{2\mu + \ep \sfrac{\rho_{\ep,\delta}^\beta}{(1-\rho_{\ep,\delta})^\beta}}\,\mathbf{1}_{\{ \rho_{\ep,\delta} < \psfrac{1+M^0}{2} \}} } \\
& \geq \int_{\TT}{\dfrac{1}{2\mu + \ep \left(\pspfrac{1+M^0}{1-M^0}\right)^\beta }\,\mathbf{1}_{\{ \rho_{\ep,\delta} < \psfrac{1+M^0}{2} \}} }\\
& \geq \dfrac{K}{2\mu + \left(\pspfrac{1+M^0}{1-M^0}\right)^\beta},
\end{align*}
where we have used the fact that $\lambda_{\ep,\delta}(\rho_{\ep,\delta})$ is bounded (uniformly in $\delta$ and $\ep$) when $\rho_{\ep,\delta}$ is far from the singularity.
\end{proof}

\begin{lem}{\label{lambdadivu}} Let us assume $u_{\ep,\delta} \in L^2(0,T,(H^1(\TT))^3)$ and $f \in L^2((0,T) \times \TT)$.
Then, we get for all $(p,q) \in [1,+\infty)^2$:
\[
p_{\varepsilon,\delta} (\rho_{\varepsilon,\delta}) \in L^p (0,T;L^q(\TT))
\Longrightarrow 
\lambda_{\varepsilon, \delta}
(\rho_{\varepsilon,\delta}) \Div u_{\varepsilon,\delta} \in L^{\min{(2,p)}}(0,T;L^{\min(2,q)}(\TT)).
\]
\end{lem}

\begin{proof} We come back to the formula
$$\lambda_{\varepsilon,\delta} (\rho_{\varepsilon,\delta}) \Div u_{\varepsilon,\delta} =
-2 \mu \Div u_{\varepsilon,\delta}
+ p_{\varepsilon,\delta} (\rho_{\varepsilon,\delta}) + S_{\ep,\delta}
- \frac{\<(p_{\varepsilon,\delta} (\rho_{\varepsilon,\delta}) + S_{\ep,\delta})\, \nu_{\varepsilon,\delta} (\rho_{\varepsilon,\delta})\>}{\<\nu_{\varepsilon,\delta} (\rho_{\varepsilon,\delta})\>}.
$$
It suffices now to use previous Lemma \ref{lem:nu_delta} to conclude.
\end{proof}

\subsubsection*{Existence of global solutions to the approximate system}

The next theorem states that one can construct global weak solutions to the truncated system at $\delta> 0$ fixed.
\begin{thm}\label{approx}
Let $0 < \delta < \delta^0 = 1-R_\ep$ and $\rho_{\ep,\delta}^0 = \rho_\ep^0$ satisfying~\eqref{hyp:bound-rho0}.
Let us assume $f \in L^2(0,T; L^{\bar{q}}(\TT))$ with $\bar{q}>3$.
Then, for all $T \in (0, +\infty)$, there exists a global weak solution $(\rho_{\ep,\delta},u_{\ep,\delta})$ to the truncated system~\eqref{eq:system-delta-mass}--\eqref{eq:system-delta-momentum}, \ie
\begin{enumerate}
\item $u_{\ep,\delta}\in L^2(0,T;(H^1(\TT))^3)$, $\rho_{\ep,\delta} \in \mathcal{C}([0,T];L^q(\TT)) \cap L^\infty((0,T) \times \TT)$ for all $q \in [1, +\infty)$;
\item $(\rho_{\ep,\delta},u_{\ep,\delta})$ satisfies \eqref{eq:system-delta-mass}--\eqref{eq:system-delta-momentum}
in the weak sense:
\begin{multline*}\labeleq{weak-delta-mass}
\intTTT{\rho_{\ep,\delta} \partial_t \phi}
+ \intTTT{\rho_{\ep,\delta} u_{\ep,\delta} \cdot \nabla \phi} \\[-5pt]
\qquad = \int_{\TT}{\rho_\ep(T) \phi(T)} - \int_{\TT}{\rho_\ep^0 \phi(0)}
\qquad \forall \phi \in \mathcal{C}^1([0,T] \times \TT);
\end{multline*} 
\begin{multline}\label{eq:weak-delta-div-mom}
- \intTTT{p_{\ep,\delta}(\rho_{\ep,\delta})\Div \psi}
+ \intTTT{(2\mu+\lambda_{\ep,\delta}(\rho_{\ep,\delta}))
\Div u_{\ep,\delta} \, \Div \psi} \\
+ \intTTT
\mu\curl u_{\ep,\delta} \cdot \curl \psi
+ r\intTTT{u_{\ep,\delta} \cdot \psi}\\
= \int_{\TT}{f \cdot \psi}
\qquad \forall \psi \in \mathcal{C}^1([0,T] \times \TT);
\end{multline} 
\item the \emph{renormalized continuity equation} holds
\begin{multline}\label{eq:renormalized-delta}
\partial_t b(\rho_{\ep,\delta}) + \Div(b(\rho_{\ep,\delta}) u_{\ep,\delta})
+ \big(b_+'(\rho_{\ep,\delta})\rho_{\ep,\delta}
- b(\rho_{\ep,\delta})\big)\Div u_{\ep,\delta} = 0
\\ \text{in}~~\mathcal{D}'((0,T)\times \TT),
\end{multline}
for any $b \in \mathcal{C}^0([0,+\infty))$, piecewise ${\cal C}^1$, where $b'_+$
denotes the right derivative of $b$;

\item the energy inequality holds
\begin{multline}\label{eq:energy_delta}
\sup_{t\in[0,T]} \int_{\TT} H_{\ep,\delta}(\rho_{\ep,\delta})
+ \intTTT{(\frac{3}{2}\mu+\lambda_{\ep,\delta}(\rho_{\ep,\delta}))(\Div u_{\ep,\delta})^2 }\\
+ \frac{\mu}{2} \intTTT{ |{\curl (u_{\ep,\delta})}|^2}
+ r \intTTT{ |u_{\ep,\delta}|^2}
\leq \int_{\TT}{H_{\ep,\delta}(\rho_\ep^0)} + \frac{1+C}{2\mu} \|f\|^2_{L^2_{t,x}}
\end{multline}
where
\begin{equation}\label{eq:H_delta}
H_{\ep,\delta}(\rho) =
\begin{cases}
\dfrac{\ep}{\gamma-1}\cdot \dfrac{\rho^{\gamma}}{(1-\rho)^{\gamma-1}} \quad & \text{if } \rho \leq 1-\delta, \\
\dfrac{\ep}{\gamma-1}\cdot \dfrac{\rho^\gamma}{\delta^\gamma} - \dfrac{\ep}{(\gamma - 1) \delta^{\gamma }}\,(1-\delta)^{\gamma} \rho \quad & \text{if } \rho > 1-\delta.
\end{cases}
\end{equation} 
\end{enumerate}
\end{thm}
Note that defining $\Lambda_{\ep,\delta}$ as
\begin{equation*}\labeleq{df:Lambda_delta}
\Lambda_{\ep,\delta}(\rho) =
\begin{cases}
\dfrac{\ep}{\beta-1}\cdot \dfrac{\rho^\beta}{(1-\rho)^{\beta-1}} \quad & \text{if } \rho \leq 1-\delta, \\
\dfrac{\ep}{\beta-1}\cdot \dfrac{\rho^\beta}{\delta^\beta} - \dfrac{\ep}{(\beta-1)\delta^{\beta}}\, (1-\delta)^\beta \rho \quad & \text{if } \rho > 1-\delta,
\end{cases}
\end{equation*}
we get the following renormalized continuity equation in $\mathcal{D}'((0,T)\times \TT)$
\begin{equation}\label{eq:Lambda_delta}
\partial_t \Lambda_{\ep,\delta}(\rho_{\ep,\delta}) + \Div (\Lambda_{\ep,\delta}(\rho_{\ep,\delta}) u_{\ep,\delta}) = - \lambda_{\ep,\delta}(\rho_{\ep,\delta})\Div u_{\ep,\delta}.
\end{equation}
The existence of global weak solutions to the approximate system, namely Theorem~\ref{approx}, follows from a standard procedure.
For reader's convenience, since our main goal is the study of the singular systems, we just present the idea of the proof.
The analysis is in fact very similar to the classical case with constant bulk viscosity treated in~\cite[Chap.\,8.2]{lions1998}.
We construct exactly in the same way the solutions by solving first the system for a regular initial data $\rho^0$ via a fixed point argument.
Then, for a general initial density $\rho^0 \in L^\infty(\TT)$, we regularize $\rho^0$ and prove that we can pass to the limit with respect to the parameter of the regularization.
Compactness arguments are needed to identify the limit quantities and in particular we need to prove the strong convergence of the sequence of densities.
The arguments to justify this strong convergence are non-standard, and different from the case with a constant bulk viscosity term, but we justify in details this point in Section~\ref{sec:delta1} for the limit $\delta \rightarrow 0$.
We refer to~\cite{lions1998} for more details.

Now we have our global weak solutions $(\rho_{\ep,\delta},u_{\ep,\delta})$, we want to pass to the limit with respect to $\delta$ (at $\ep$ fixed) to get global existence of weak solutions for the compressible singular systems.
It will be then possible to pass to the limit with respect to $\varepsilon$ to get the congestion systems.\\
We will divide the study in two sections depending on the sign of $\gamma-\beta$. First, we treat the dominant pressure regime $\gamma > \beta$ (Section~\ref{sec:pressure}), then the dominant bulk viscosity regime $\beta \ge \gamma$ (Section~\ref{sec:bulk}).

\section{Dominant pressure regime \texorpdfstring{$\gamma > \beta > 1$}{gb1}}\label{sec:pressure}

\Subsection{Existence of weak solutions at $\ep$ fixed}

\subsubsection{Uniform estimates with respect to $\delta$}

First of all, observe that if we consider $\delta$ small enough, namely $\delta < 1-R_\ep$ ($\ep$ is fixed), we ensure that initially
\[
H_{\ep,\delta}(\rho_{\ep,\delta}^0) = H_{\ep}(\rho_{\ep}^0) \text{ is bounded in } L^1(\TT).
\]
Thanks to Theorem~\ref{approx}, the solutions $(\rho_{\ep,\delta},u_{\ep,\delta})$ satisfy the energy estimate~\eqref{eq:energy_delta}, so that $(u_{\ep,\delta})_\delta$ is bounded in $L^2(0,T;(H^1(\TT))^3)$
and $(H_{\ep,\delta}(\rho_{\ep,\delta}))_\delta$ is bounded in $L^\infty(0,T;L^1(\TT))$.
In particular, the control of the internal energy $H_{\ep,\delta}$ leads easily to a control of the density\vspace*{-3pt}
\begin{equation*}\labeleq{bound-rd-Lgamma}
(\rho_{\ep,\delta})_\delta \text{ is bounded in } L^{\infty}(0,T;L^\gamma(\TT)).
\end{equation*}
In the following lemma, we improve the control of the density by using the singularity of $H_{\ep,\delta}$ with respect to $\delta$ for large values of the density.
\begin{lem}\label{lem:meas_largerho_delta}
Let $(\rho_{\varepsilon, \delta}, u_{\varepsilon,\delta})$ be a global weak solution
of the compressible Brinkman system. Then\vspace*{-3pt}
\begin{equation}\label{eq:est-rhodelta}
\sup_{t\in [0,T]} \meas \big\{ x\in \TT, \ \rho_{\ep,\delta}(t,x) \geq 1-\delta \big\} \leq C(\ep) \ \delta^{\gamma -1}.
\end{equation}
\end{lem}

\begin{proof} The energy $H_{\ep,\delta}$ being defined as~\eqref{eq:H_delta}, recalling that $\gamma>1$, we have\vspace*{-3pt}
\begin{multline} \label{eq:ineg-Hdelta}
C \geq \sup_{t\in [0,T]} \int_{\TT}{H_{\ep,\delta}(\rho_{\ep,\delta})}
\geq \int_{\TT}{H_{\ep,\delta}(\rho_{\ep,\delta}) \,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}} \\[-2pt]
\geq C\int_{\TT}{\dfrac{\ep}{\delta^\gamma}\bigl[
(\rho_{\varepsilon,\delta}^{\gamma-1} - (1-\delta)^{\gamma-1})
+ (1-\delta)^{\gamma-1}\delta \bigr] \rho_{\ep,\delta}
\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}} \\[-3pt]
\geq C\int_{\TT}{\dfrac{\ep}{\delta^{\gamma-1} }\bigl[
(1-\delta)^{\gamma} \bigr]
\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}},
\end{multline}
which ends the proof.
\end{proof}

\begin{lem} Let $(\rho_{\varepsilon, \delta}, u_{\varepsilon,\delta})$ be a global weak
of the compressible Brinkman system~\eqref{eq:semi-stat-delta} with $\gamma > \beta > 1$. Then\vspace*{-3pt}
\begin{multline*}
\|\Lambda_{\varepsilon,\delta}(\rho_{\ep,\delta})\|_{L^1((0,T)\times \TT)} +
\|p_{\varepsilon,\delta}(\rho_{\ep,\delta})\|_{L^1((0,T)\times \TT)}\\
+
\|\lambda_{\varepsilon,\delta}(\rho_{\ep,\delta}) \Div u_{\varepsilon, \delta}\|_{L^1((0,T)\times \TT)}
+ \|\lambda_{\varepsilon,\delta}(\rho_{\varepsilon,\delta})\|_{L^{\gamma/\beta}((0,T)\times \TT)}
\le C_\varepsilon,
\end{multline*}
where $C_\varepsilon$ does not depend on $\delta$.
\end{lem}

\noindent {\it Proof.}
Integrating in space~\eqref{eq:Lambda_delta} and using~\eqref{eq:lambda_divu}, we have
\begin{equation*}
\Dt \int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta})} + \int_{\TT}{p_{\ep,\delta}(\rho_{\ep,\delta})}
= |\TT|\,\dfrac{\int_{\TT}{\big(S_{\ep,\delta}+p_{\ep,\delta}(\rho_{\ep,\delta})\big)\nu_{\ep,\delta}(\rho_{\ep,\delta})}}{\int_{\TT}{\nu_{\ep,\delta}(\rho_{\ep,\delta})}}.
\end{equation*}
Using Lemma~\ref{lem:nu_delta}, we can bound the right-hand side
\begin{equation}\label{eq:dt-Lambda-delta}
\Dt \int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta})} + \int_{\TT}{p_{\ep,\delta}(\rho_{\ep,\delta})}
\leq \dfrac{C_1|\TT|}{C_2}\,\|S_{\ep,\delta}\|_{L^1} + \dfrac{|\TT|}{C_2}\int_{\TT}{p_{\ep,\delta}(\rho_{\ep,\delta})\nu_{\ep,\delta}(\rho_{\ep,\delta})},
\end{equation}
where
\[
\|S_{\ep,\delta}\|_{L^1} = \|(-\Delta)^{-1}\Div \big(f-r u_{\ep,\delta}\big) \|_{L^1}
\leq C \big(\|f\|_{L^2(0,T; L^{\bar{q}}(\TT))} + \|u_{\ep,\delta}\|_{L^2H^1}\big),
\]
and
\begin{multline*}
p_{\ep,\delta}(\rho_{\ep,\delta})\nu_{\ep,\delta}(\rho_{\ep,\delta})
= \dfrac{p_{\ep,\delta}(\rho_{\ep,\delta})}{2\mu + \lambda_{\ep,\delta}(\rho_{\ep,\delta})} \\
\leq \dfrac{p_{\ep,\delta}(\rho_{\ep,\delta})}{2\mu}\,\mathbf{1}_{\{\rho_{\ep,\delta} < M_0 \}}
+ \dfrac{p_{\ep,\delta}(\rho_{\ep,\delta})}{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}\,\mathbf{1}_{\{M_0 \leq \rho_{\ep,\delta} < 1-\delta \}}
+ \dfrac{p_{\ep,\delta}(\rho_{\ep,\delta})}{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}\,\mathbf{1}_{\{\rho_{\ep,\delta} \geq 1-\delta \}}.
\end{multline*}
The first term of the right-hand side is bounded since $\rho_{\ep,\delta}$
is far from $1$.
For the two other terms, which become singular as $\delta \rightarrow 0$, we ensure that
\begin{align*}
\dfrac{p_{\ep,\delta}(\rho_{\ep,\delta})}{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}\,\mathbf{1}_{\{M_0 \leq \rho_{\ep,\delta} < 1-\delta \}}
& \leq \dfrac{C}{(1-\rho_{\ep,\delta})^{\gamma-\beta}}\,\mathbf{1}_{\{M_0 \leq \rho_{\ep,\delta} < 1-\delta \}}\\
& \leq C(\ep)H_{\ep,\delta}(\rho_{\ep,\delta})\,\mathbf{1}_{\{M_0 \leq \rho_{\ep,\delta} < 1-\delta \}} \labeleq{estim-p-v0},\\
\dfrac{p_{\ep,\delta}(\rho_{\ep,\delta})}{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}\,\mathbf{1}_{\{\rho_{\ep,\delta} \geq 1-\delta \}}
& \leq C \dfrac{\rho_{\ep,\delta}^{\gamma-\beta}}{\delta^{\gamma-\beta}}\,\mathbf{1}_{\{\rho_{\ep,\delta} \geq 1-\delta \}} \\
& \leq C(\ep)\bigl[H_{\ep,\delta}(\rho_{\ep,\delta})
+ \rho\bigr]\,\mathbf{1}_{\{\rho_{\ep,\delta} \geq 1-\delta \}}
\end{align*}
since $\beta \in (1,\gamma)$ and thus $0 < \gamma - \beta < \gamma - 1$ (recall Definition~\eqref{eq:H_delta} of $H_{\ep,\delta}$).
Using now the control of $H_{\ep,\delta}$ and the fact that the total mass is constant, we deduce ($\ep$ is fixed here)
\[
\Dt \int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta})} + \int_{\TT}{p_{\ep,\delta}(\rho_{\ep,\delta})}
\leq C(\ep).
\]
To conclude, let us observe that, using that $\beta <\gamma$ and the initial conditions~\eqref{hyp:bound-rho0} and~\eqref{hyp:energy-t0-C1}, we have
\begin{equation*}
\int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta}(0,\cdot)) }
= \int_{\TT}{\Lambda_{\ep}(\rho_\ep^0) }
\leq C.
\end{equation*}
Hence, we get from integration in time of~\eqref{eq:dt-Lambda-delta} that
\begin{equation*}\labeleq{bound-p-delta}
\big(p_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta~ \text{is bounded in} ~L^1\big((0,T)\times \TT \big).
\end{equation*}
Coming back to Lemma \ref{lambdadivu}, we obtain
\begin{equation*}\label{bound-lambdadivu-delta}
\big(\lambda_{\ep,\delta}(\rho_{\ep,\delta}) \Div u_{\ep,\delta}\big)_\delta \text{ bounded in}~ L^1\big((0,T)\times \TT\big).
\end{equation*}
Note that from the pressure estimate, since $\gamma > \beta$, we deduce that
\begin{equation*}\label{bound-lambda-delta}
\big(\lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta ~ \text{is bounded in} ~L^{\sfrac{\gamma}{\beta}}\big((0,T)\times \TT \big).
\end{equation*}

These controls on the pressure and the bulk viscosity can now be used to prove a maximal bound on the density.
\begin{prop}\label{prop:bound_rhodelta}
The density $\rho_{\ep,\delta}$ is bounded in $L^\infty((0,T)\times \TT)$ uniformly with respect to the cut-off parameter $\delta$.
\end{prop}

\begin{proof}
Using the renormalized continuity equation \eqref{eq:renormalized-delta} with $b(\rho) = \rho^m$ and $m \in (1,+\infty)$, we get
\[
\partial_t \rho_{\ep,\delta}^m + \Div (u_{\ep,\delta} \rho_{\ep,\delta}^m)
= (1-m) \, \rho_{\ep,\delta}^m \, \Div u_{\ep,\delta} \quad \text{in } \mathcal{D}'((0,T)\times \TT)
\]
and therefore
\[
\frac{d}{dt} \int_{\TT} \rho_{\ep,\delta}^m = (1-m)\, \int_{\TT} \rho_{\ep,\delta}^m \Div u_{\ep,\delta}.
\]
Using~\eqref{eq:expr_lambdadivu} to replace $\Div u_{\ep,\delta}$, we get
\begin{multline*}
\frac{d}{dt} \int_{\TT}\rho_{\ep,\delta}^m + (m-1) \int_{\TT} \rho_{\ep,\delta}^m
\nu_{\ep,\delta}(\rho_{\ep,\delta}) p_{\ep,\delta}(\rho_{\ep,\delta})
= (1-m) \int_{\TT} \rho_{\ep,\delta}^m
\nu_{\ep,\delta}(\rho_{\ep,\delta}) S_{\ep,\delta}\\
+ (1-m) \int_{\TT} \dfrac{\rho_{\ep,\delta}^m \nu_{\ep,\delta}(\rho_{\ep,\delta})}{|\TT|} \int_{\TT}{ \big((\lambda_{\ep,\delta})(\rho_{\ep,\delta}) \Div u_{\ep,\delta}
- p_{\ep,\delta}(\rho_{\ep,\delta}) \big) }.
\end{multline*}
Now, thanks to~\eqref{hyp:f} ($\bar{q} > 3$), we ensure that
\[
\|S_{\varepsilon,\delta}\|_{L^\infty_x} = \|(-\Delta)^{-1}\Div (f-ru_{\ep,\delta})\|_{L^\infty_x} \leq C \big(\|f\|_{L^{\bar{q}}_x}+ \|u_{\ep,\delta}\|_{L^6_x} \big)
\]
and therefore (we recall that $\|\nu_{\ep,\delta}\|_{L^\infty_{t,x}} \leq C_1$)
\begin{align*}
\frac{d}{dt} \int_{\TT} \rho_{\varepsilon,\delta}^m
\leq 2 \, C_1 \, m \Big(\|S_{\ep,\delta}\|_{L^\infty_x}
+ \| \lambda_{\ep,\delta}(\rho_{\ep,\delta}) \Div u_{\ep,\delta} \|_{L^1_x}
+ \|p_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^1_x} \Big)
\int_{\TT} \rho_{\ep,\delta}^m.
\end{align*}
This gives
\begin{align*}
\left(\int_{\TT} \rho_{\varepsilon,\delta}^m(t)\right)^{1/m} \leq
\left(\int_{\TT} \rho_{\varepsilon,\delta}^m(0)\right)^{1/m} \exp (\overline C) \qquad \forall t \geq 0,
\end{align*} 
with
\[
\overline{C} = 2\, C_1 \exp\Bigl(\|S_{\ep,\delta}\|_{L^1_t L^\infty_x}
+ \| (\lambda_{\ep,\delta})(\rho_{\ep,\delta}) \Div u_{\ep,\delta} \|_{L^1_{t,x}}
+ \|p_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^1_{t,x}}\Bigr).
\]
Hence, passing to the limit with respect to $m$ and recalling that $0 \leq \rho_{\ep,\delta}(0) = \rho_\ep^0 < 1$, we get the uniform upper bound:
\[
\|\rho_{\ep,\delta}\|_{L^\infty_{t,x}} \leq C.\qedhere
\]
\end{proof}

We now improve a little bit the estimate on the pressure.
This will ensure that the weak limit of the pressure is more regular than a measure.
\begin{lem}{\label{lem:estim-pressure}}
The sequence $(p_{\ep,\delta}(\rho_{\ep,\delta}))_\delta$ is bounded in $L^{1+\theta}\big((0,T)\times \TT \big)$.
\end{lem}
\begin{proof}
Let us consider in~\eqref{eq:weak-delta-div-mom} the test function
\[
\psi = \nabla \Delta^{-1}\big((H_{\ep,\delta}(\rho_{\ep,\delta}))^\alpha - \<(H_{\ep,\delta}(\rho_{\ep,\delta}))^\alpha\> \big) \quad \text{with } \alpha = \dfrac{\gamma - \beta}{\gamma-1} \in (0,1).
\]
We have
\begin{multline*}
\intTTT{p_{\ep,\delta}(\rho_{\ep,\delta}) \left[\big(H_{\ep,\delta}(\rho_{\ep,\delta}) \big)^{\alpha} - \<(H_{\ep,\delta}(\rho_{\ep,\delta}))^\alpha\> \right] } \\
= \intTTT{\Big(2\mu + \lambda_{\ep,\delta}(\rho_{\ep,\delta})\Big)\Div u_{\ep,\delta} \left[\big(H_{\ep,\delta}(\rho_{\ep,\delta}) \big)^{\alpha} - \<(H_{\ep,\delta}(\rho_{\ep,\delta}))^\alpha\> \right]} \\
- \intTTT{\big(f-ru_{\ep,\delta}\big) \cdot \psi }
\end{multline*}
and using the controls resulting from the energy inequality, we obtain
\begin{multline*}
\intTTT{p_{\ep,\delta}(\rho_{\ep,\delta}) \left[\big(H_{\ep,\delta}(\rho_{\ep,\delta}) \big)^{\alpha} - \<(H_{\ep,\delta}(\rho_{\ep,\delta}))^\alpha\> \right] } \\
\leq C + \intTTT{\big(H_{\ep,\delta}(\rho_{\ep,\delta}) \big)^{\alpha}\lambda_{\ep,\delta}(\rho_{\ep,\delta})\,|{\Div(u_{\ep,\delta})}| }.
\end{multline*}
Since $\alpha < 1$ and $H_{\ep,\delta}(\rho_{\ep,\delta})$ is bounded in $L^\infty(0,T;L^1(\TT))$,
we control
\[
\biggl|\intTTT{p_{\ep,\delta}(\rho_{\ep,\delta}) \<(H_{\ep,\delta}(\rho_{\ep,\delta}))^\alpha\> }\biggr|
\leq \|H_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^\infty L^1}^{\alpha} \|p_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^1L^1}
\]
and we get then
\begin{align*}
\intTTT{p_{\ep,\delta}(\rho_{\ep,\delta}) \big(H_{\ep,\delta}(\rho_{\ep,\delta}) \big)^{\alpha}}
& \leq C + \|H_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^\infty L^1}^{\alpha} \|p_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^1L^1} \\
& \hspace*{1cm} + \intTTT{\big(H_{\ep,\delta}(\rho_{\ep,\delta}) \big)^{\alpha}\lambda_{\ep,\delta}(\rho_{\ep,\delta})\,|{\Div(u_{\ep,\delta})}| } \\
& \leq C + \intTTT{\big(H_{\ep,\delta}(\rho_{\ep,\delta}) \big)^{2\alpha}\lambda_{\ep,\delta}(\rho_{\ep,\delta}) \,\mathbf{1}_{\{ \rho_{\ep,\delta} \geq M^0 \}} }.
\end{align*}
In the right-hand side we have
\begin{align*}
\big(H_{\ep,\delta}(\rho_{\ep,\delta})\big)^{2\alpha}&\lambda_{\ep,\delta}(\rho_{\ep,\delta}) \,\mathbf{1}_{\{ \rho_{\ep,\delta} \geq M^0 \}}\\
&\leq \dfrac{1}{(\gamma-1)^{2\alpha}}\,\dfrac{\ep^{2\alpha + 1}}{(1-\rho_{\ep,\delta})^{2\alpha (\gamma-1) + \beta}} \ \rho_{\ep,\delta}^{2\alpha\gamma + \beta}\,\mathbf{1}_{\{ \rho_{\ep,\delta} \geq M^0 \}} \\
& \leq \dfrac{1}{(\gamma-1)^{\alpha}}\,\left(H_{\ep,\delta}(\rho_{\ep,\delta})\right)^{\alpha} p_{\ep,\delta}(\rho_{\ep,\delta}) \ \ep^\alpha \dfrac{\rho_{\ep,\delta}^{(\alpha-1) \gamma + \beta}}{(1-\rho_{\ep,\delta})^{\alpha(\gamma-1)+\beta-\gamma}}\,\mathbf{1}_{\{ \rho_{\ep,\delta} \geq M^0 \}}.
\end{align*}
For $\alpha = \pspfrac{\gamma - \beta}{\gamma-1}$, we have
\[
(\alpha-1)\gamma+\beta = \dfrac{\gamma-\beta}{\gamma-1} > 0,
\qquad \alpha(\gamma-1)+\beta-\gamma = 0.
\]
So, using the $L^\infty$ bound on $\rho_{\ep,\delta}$, we obtain
\begin{align*}
&\big(H_{\ep,\delta}(\rho_{\ep,\delta})\big)^{2\alpha}\lambda_{\ep,\delta}(\rho_{\ep,\delta}) \,\mathbf{1}_{\{ \rho_{\ep,\delta} \geq M^0 \}}
< \ C\ep^\alpha \left(H_{\ep,\delta}(\rho_{\ep,\delta})\right)^{\alpha} p_{\ep,\delta}(\rho_{\ep,\delta}) \,\mathbf{1}_{\{ \rho_{\ep,\delta} \geq M^0 \}},
\end{align*}
with $\ep< \ep^0$ small enough, so that this term is absorbed in the left-hand side of the previous inequality.
Finally
\begin{equation*}\labeleq{add_control_p_delta}
\left(H_{\ep,\delta}(\rho_{\ep,\delta})\right)^{\alpha} p_{\ep,\delta}(\rho_{\ep,\delta}) \text{ is bounded in } L^1\big((0,T)\times \TT\big).
\end{equation*} 
The previous bound allows us to show that $(p_{\ep,\delta}(\rho_{\ep,\delta}))_\delta$ is even bounded in $L^{1+\theta}((0,T)\times \TT)$ for some $\theta > 0$.
Indeed, we have for $\rho_{\varepsilon,\delta} \ge 1-\delta$ (see \eqref{eq:ineg-Hdelta}
and the expression of $p_{\varepsilon,\delta}(\rho_{\varepsilon,\delta})$ given in \eqref{eq:p_delta}):
\begin{align*}
H_{\varepsilon,\delta}(\rho_{\varepsilon,\delta}) 1_{\{\rho_{\varepsilon,\delta} \geq 1-\delta\}}
& \geq C \,\dfrac{\ep}{\delta^\gamma}\, \Bigl[
(\rho_{\varepsilon,\delta}^{\gamma-1} - (1-\delta)^{\gamma-1})
+ (1-\delta)^{\gamma-1}\delta \Bigr] \rho_{\ep,\delta}
\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}\\
& \geq C \,\dfrac{\ep}{\delta^\gamma}\, (1-\delta)^{\gamma-1}\delta \rho_{\ep,\delta}
\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}\\
& \geq C\,\frac{\varepsilon^{1-1/\gamma}}{\delta^{\gamma-1}}\,\delta \big(p_{\varepsilon,\delta}(\rho_{\varepsilon,\delta})\big)^{1/\gamma}\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}.
\end{align*}
Observing that
\[
\delta \, \big(p_{\varepsilon,\delta}(\rho_{\varepsilon,\delta})\big)^{1/\gamma}\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}
\ge \varepsilon^{1/\gamma}(1-\delta)
\ge \varepsilon^{1/\gamma}(1-\delta^0),
\]
and introducing $\eta = \min(1,\gamma-1)$, we get that
\begin{align*}
H_{\varepsilon,\delta}(\rho_{\varepsilon,\delta}) 1_{\{\rho_{\varepsilon,\delta} \geq 1-\delta\}}
& \geq C \,\frac{\varepsilon^{1-1/\gamma}}{\delta^{\gamma-1}}\,
\big[\delta \big(p_{\varepsilon,\delta}(\rho_{\varepsilon,\delta})\big)^{1/\gamma}\big]^{1-\eta}
\big[\delta \big(p_{\varepsilon,\delta}(\rho_{\varepsilon,\delta})\big)^{1/\gamma}\big]^{\eta}
\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}} \\
& \geq C\, \frac{\varepsilon^{1-\sfrac{\eta}{\gamma}}(1-\delta^0)^{1-\eta}}{\delta^{\gamma-1-\eta}}\, \big(p_{\varepsilon,\delta}(\rho_{\varepsilon,\delta})\big)^{\eta/\gamma}
\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}.
\end{align*}
Hence
\[
(H_{\varepsilon,\delta}(\rho_{\varepsilon, \delta}))^{\alpha}
p_{\varepsilon,\delta}(\rho_{\varepsilon,\delta})\,
\,\mathbf{1}_{\{\rho_{\varepsilon,\delta} \ge 1-\delta\} }
\ge C_{\varepsilon,\delta^0,\eta}\, \frac{1}{\delta^{\alpha(\gamma-1-\eta)}}\,
(p_{\varepsilon,\delta}(\rho_{\varepsilon,\delta}))^{1+\alpha\eta/\gamma}
\,\mathbf{1}_{\{\rho_{\varepsilon,\delta} \ge 1-\delta\} }.
\]
On the set $\{\rho_{\ep,\delta} \leq M^0\}$ the pressure is of course bounded, so it remains to consider the set $\{M^0 < \rho_{\varepsilon,\delta} < 1-\delta\}$ on which we have
\begin{align*}
H_{\varepsilon,\delta}(\rho_{\varepsilon, \delta})\,\mathbf{1}_{\{\rho_{\varepsilon,\delta} < 1-\delta\}}
& = \dfrac{\ep}{\gamma-1}\, \dfrac{\rho_{\ep,\delta}^\gamma}{(1-\rho_{\ep,\delta})^{\gamma}}\,\mathbf{1}_{\{M^0 < \rho_{\varepsilon,\delta} < 1-\delta\}} \\
& = \dfrac{\ep^{1/\gamma}}{\gamma-1}\, \rho_{\ep,\delta} \big(p_{\ep,\delta}(\rho_{\ep,\delta})\big)^{\psfrac{\gamma-1}{\gamma}} \,\mathbf{1}_{\{M^0 < \rho_{\varepsilon,\delta} < 1-\delta\}} \\
& \geq C_{\ep,M^0} \big(p_{\ep,\delta}(\rho_{\ep,\delta})\big)^{\psfrac{\gamma-1}{\gamma}} \,\mathbf{1}_{\{M^0 < \rho_{\varepsilon,\delta} < 1-\delta\}}.
\end{align*}
We obtain again a control of the pressure (we recall that $\gamma > 1$)
\[
(H_{\varepsilon,\delta}(\rho_{\varepsilon, \delta}))^{\alpha}
p_{\varepsilon,\delta}(\rho_{\varepsilon,\delta})\,
\,\mathbf{1}_{\{M^0 <\rho_{\varepsilon,\delta} < 1-\delta\} }
\geq C_{\varepsilon}\, (p_{\varepsilon,\delta}(\rho_{\varepsilon,\delta}))^{1+\alpha\psfrac{\gamma-1}{\gamma}}
\,\mathbf{1}_{\{M^0 <\rho_{\varepsilon,\delta} < 1-\delta\} }.
\]
This concludes the proof of Lemma \ref{lem:estim-pressure}.
\end{proof}

\Subsubsection{Limit $\delta \rightarrow 0$}{\label{sec:delta1}}

\subsubsection*{First convergence results}
Thanks to the estimates we have just derived, there exists a limit density $\rho_\ep$ such that
\begin{equation*}
\rho_{\ep,\delta} \longrightharpoonup \rho_\ep \quad \text{weakly-* in} ~ L^\infty\big((0,T)\times \TT \big)
\end{equation*}
and, passing to the limit $\delta \rightarrow 0$ in~\eqref{eq:est-rhodelta} we get
\begin{equation}\label{eq:bound_rho_ep}
0 \leq \rho_\ep(t,x) < 1 \quad \text{a.e. } (t,x) \in [0,T]\times \TT.
\end{equation}
In addition, there exists a limit velocity $u_\ep$ such that
\begin{equation*}
u_{\ep,\delta} \longrightharpoonup u_\ep \quad \text{weakly in } L^2(0,T;(H^1(\TT))^3)
\end{equation*}
and, due to the continuity equation, we have
\begin{equation*}
\rho_{\ep,\delta} \rightarrow \rho_\ep \quad \text{in } \mathcal{C}_{\text{weak}}([0,T],L^r(\TT)) \quad \forall r \in [1,+\infty).
\end{equation*}
To identify the weak limit of the nonlinear term $\rho_{\ep,\delta} u_{\ep,\delta}$, we use the next compensated compactness Lemma.
\begin{lem}[{\cite[Lem.\,5.1]{lions1998}}]\label{lem:compensated_compactness}
Let $(g_n)$, $(h_n)$ be two sequences converging respectively to $g$, $h$ in $L^{r_1}(0,T;L^{q_1}(\TT))$ and $L^{r_2}(0,T;L^{q_2}(\TT))$ where $1\leq r_1,r_2\leq \infty$ and \hbox{$\sfrac{1}{r_1} + \sfrac{1}{r_2} = \sfrac{1}{q_1} + \sfrac{1}{q_2}$}. Assume in addition that
\begin{enumerate}
\item $(\partial_t g_n)_n$ is bounded in $L^1(0,T;W^{-m,1}(\TT))$ for some $m$ independent of $n$;
\item $\|h_n\|_{L^1_tH^s_x}$ is bounded for some $s >0$.
\end{enumerate}
Then $g_n h_n$ converges to $gh$ weakly in $\mathcal{D}'((0,T)\times \TT)$.
\end{lem}

We apply the result to $g_{\delta} = \rho_{\ep,\delta}$, $h_{\delta} = u_{\ep,\delta}$: we ensure the control of $\partial_t \rho_{\ep,\delta}$ in $L^2\big(0,T;H^{-1}(\TT)\big)$ from the continuity equation, while $(\nabla u_{\ep,\delta})_\delta$ is bounded in $L^2\big((0,T)\times \TT\big)$ thanks to the energy inequality.
Hence
\begin{equation*}
\rho_{\ep,\delta} u_{\ep,\delta} \longrightharpoonup \rho_\ep u_\ep \quad \text{weakly-* in } L^\infty(0,T;L^6(\TT)).
\end{equation*}
With the estimates on the pressure we deduce
\begin{equation*}
p_{\ep,\delta}(\rho_{\ep,\delta}) \longrightharpoonup \ov{p_\ep(\rho)} \quad \text{weakly in } L^{1}\big((0,T) \times \TT\big),
\end{equation*}
where $\overline{h}$ denotes the weak limit of the sequence $(h_\delta)_\delta$.
Our next goal is to get the strong convergence of the density $\rho_{\ep,\delta}$ in order to identify the limit of the pressure and the bulk viscosity which are non-linear functions of the density.

\subsubsection*{Convergence a.e. of the density.}
Thanks to the bounds on $\rho_{\ep,\delta}, u_{\ep,\delta}$, we can pass to the limit in the sense of distributions in the renormalized continuity equation
\[
\partial_t \rho_{\ep,\delta}^2 + \Div\big(\rho_{\ep,\delta}^2 u_{\ep,\delta} \big)
= - \rho_{\ep,\delta}^2 \Div u_{\ep,\delta},
\]
which reads at the limit
\[
\partial_t \ov{\rho_{\ep}^2} + \Div\big(\ov{\rho_{\ep}^2} u_{\ep} \big)
= - \ov{\rho_{\ep}^2 \Div u_{\ep}}.
\]
On the other hand, the limit density $\rho_\ep \in L^\infty((0,T)\times \TT)$ satisfies the renormalized continuity equation
\[
\partial_t \rho_{\ep}^2 + \Div\big(\rho_{\ep}^2 u_{\ep} \big)
= - \rho_{\ep}^2 \Div u_{\ep}.
\]
Defining $\Psi := \ov{\rho_\ep^2} - \rho_\ep^2 \geq 0$, we have then
\begin{equation*}
\partial_t \Psi + \Div(\Psi u) = \rho_\ep^2 \Div u_\ep - \ov{\rho_\ep^2 \Div u_\ep} \qquad \text{in } \mathcal{D}'.
\end{equation*}
By replacing $\Div u_\ep$ by its expression in terms of effective flux and pressure, the previous equation can be rewritten as
\begin{equation}\label{eq:Psi-0}
\partial_t \Psi + \Div(\Psi u)
 = \rho_\ep^2\, \ov{\nu_\ep (\rho_\ep) F_\ep } - \ov{\rho_\ep^2 \nu_\ep (\rho_\ep) F_\ep } + \rho_\ep^2\, \ov{ p_\ep(\rho) \nu_\ep(\rho_\ep) } - \ov{ \rho_\ep^2 p_\ep(\rho) \nu_\ep(\rho_\ep) }.
\end{equation}

\begin{remark}
In the classical case of constant bulk and shear viscosities $\mu^0$, $\lambda^0$, the previous equation writes
\[
\partial_t \Psi + \Div(\Psi u)
= \dfrac{1}{2\mu^0 + \lambda^0} \bigl[\rho_\ep^2\, \ov{F_\ep } - \ov{\rho_\ep^2 F_\ep }
+ \rho_\ep^2\, \ov{ p(\rho)} - \ov{ \rho_\ep^2 p(\rho)}
\bigr],
\]
and one can prove some weak compactness property of the effective flux (see~\cite[Chap.\,5]{lions1998} or~\cite{novotny2004}).
This property ensures that
\[
\rho_\ep^2\, \ov{F_\ep } = \ov{\rho_\ep^2 F_\ep},
\]
so that
\[
\partial_t \Psi + \Div(\Psi u) = \dfrac{1}{2\mu^0 + \lambda^0}\,\bigl(\rho_\ep^2\, \ov{ p(\rho)} - \ov{ \rho_\ep^2 p(\rho)}\bigr).
\]
In the usual case where the pressure is a monotone (increasing) function, independent of the parameter $\delta$, then (see~\cite[Lem.\,3.35]{novotny2004})
\begin{equation}\label{eq:monotone-press}
\rho_\ep^2\, \ov{ p(\rho)} \leq \ov{ \rho_\ep^2 p(\rho)}
\end{equation}
and\vspace*{-3pt}
\[
\partial_t \Psi + \Div(\Psi u) \leq 0.
\]
We conclude by integrating over space
\[
\Dt \int_{\TT}\Psi \leq 0.
\]
Recall that by the convexity of the functional $s \mapsto s^2$ we have $\Psi \geq 0$.
Hence, if initially $\Psi(0,\cdot) = 0$, we obtain
\[
\Psi = 0 \quad \text{a.e.} ~(t,x).
\]
This ensures the strong convergence of $(\rho_\delta)_\delta$.
Note finally that this calculation has been extended by {Feireisl} in~\cite{feireisl2002} to non-monotone pressure that are increasing only from a critical density. In this case, one controls the part where the pressure is non-monotone in such way that a Gronwall inequality can be applied to recover at the end $\Psi = 0$ a.e.
We will see below that we will have to use with such kind of arguments to prove the strong convergence of the density in case of density dependent bulk viscosities.
We refer the reader to~\cite{bresch2018} for recent developments on more general non-monotone pressures.
\end{remark} 

Our study is original in two ways: first, we have here to deal with a density dependent bulk viscosity, secondly, the pressure (as well as the bulk viscosity) depends on the parameter of approximation $\delta$.
We begin with proving some similar weak compactness properties satisfied by the effective flux.

\begin{prop}\label{prop:compactness_F}
We ensure the two following properties
\begin{align}\label{prop:compactness_F_1}
\ov{\Bigl(\dfrac{\rho^2 F_\ep}{2\mu +\lambda_\ep(\rho)} \Bigr)} &= \ov{\Bigl(\dfrac{\rho^2}{2\mu +\lambda_\ep(\rho)}\Bigr)}~ \ov{F_\ep}\qquad \text{in } \mathcal{D}',\\[3pt]
\label{prop:compactness_F_2}
\ov{\Bigl(\dfrac{F_\ep}{2\mu +\lambda_\ep(\rho)}\Bigr)} &= \ov{\Bigl(\dfrac{1}{2\mu +\lambda_\ep(\rho)}\Bigr)}~ \ov{F_\ep}\qquad \text{in } \mathcal{D}',
\end{align}
where $\ov{g}$ denotes the weak limit of the sequence $(g_\delta)$.
\end{prop}

\begin{proof}
The proof of these properties follows again from Lemma~\ref{lem:compensated_compactness} with $h_{\ep,\delta} = F_{\ep,\delta}$ and $g_{\ep,\delta}^1 = \rho_{\ep,\delta}^2 \nu_{\ep,\delta}(\rho_{\ep,\delta})$, $g_{\ep,\delta}^2=\nu_{\ep,\delta}(\rho_{\ep,\delta})$.
Let us check that we control the time derivative of $g_{\ep,\delta}^1$ and $g_{\ep,\delta}^2$ that satisfy the renormalized continuity equations
\begin{equation}\label{eq:renorm_gi}
\partial_t g_{\ep,\delta}^i(\rho_{\ep,\delta}) + \Div\big(g_{\ep,\delta}^i(\rho_{\ep,\delta}) u_{\ep,\delta}\big) + \Big((g_{\ep,\delta}^i)'_+(\rho_{\ep,\delta}) \, \rho_{\ep,\delta}-g_{\ep,\delta}^i(\rho_{\ep,\delta}) \Big) \Div u_{\ep,\delta} = 0
\end{equation}
with
\[ (g_{\ep,\delta}^1)'_+(\rho) =
\dfrac{2\rho}{2\mu + \lambda_{\ep,\delta}(\rho)}-\dfrac{(\lambda_{\ep,\delta})'_+(\rho)}{(2\mu + \lambda_{\ep,\delta}(\rho))^2}, \qquad (g_{\ep,\delta}^2)'_+(\rho)
= -\dfrac{(\lambda_{\ep,\delta})'_+(\rho)}{(2\mu + \lambda_{\ep,\delta}(\rho))^2}.\]
Recall that
\[
(\lambda_{\ep,\delta})'_+(\rho) = \begin{cases}
\ep \beta \sfrac{\rho^{\beta-1}}{(1-\rho)^{\beta+1}} \quad & \text{if } \rho < 1-\delta, \\
\ep \beta \sfrac{\rho^{\beta-1}}{\delta^\beta}\quad & \text{if } \rho \geq 1-\delta, \end{cases}
\]
hence
\begin{equation}\label{eq:inegal-dep-ep}
\begin{split}
\dfrac{(\lambda_{\ep,\delta})'_+(\rho)}{(2\mu + \lambda_{\ep,\delta}(\rho))^2}
& \leq C \,\mathbf{1}_{\{ \rho_{\ep,\delta} < M^0 \}}
+ \dfrac{(\lambda_{\ep,\delta})'_+(\rho)}{(\lambda_{\ep,\delta}(\rho))^2}\,\mathbf{1}_{\{ M^0 \leq \rho_{\ep,\delta} < 1-\delta \}}\\[-5pt]
&\hspace*{5cm}+ \dfrac{(\lambda_{\ep,\delta})'_+(\rho)}{(\lambda_{\ep,\delta}(\rho))^2}\,\mathbf{1}_{\{ \rho_{\ep,\delta} \geq 1-\delta \}} \\
& \leq C + \dfrac{C}{\ep}(1-\rho_{\ep,\delta})^{\beta -1}\,\mathbf{1}_{\{ M^0 \leq \rho_{\ep,\delta} < 1-\delta \}}
+ \dfrac{C}{\ep}\delta^{\beta}\,\mathbf{1}_{\{ \rho_{\ep,\delta} \geq 1-\delta \}} \\
& \leq C(\ep).
\end{split}
\end{equation}
As a consequence, the
$g^i_{\ep,\delta}$ and $(g^i_{\ep,\delta})'$ are bounded in $L^2\big((0,T)\times \TT\big)$ uniformly with respect to $\delta$ (but not uniformly in $\ep$, {\it cf.} Remark~\ref{rmk:cvg-rhoep}).
By a Cauchy-Schwarz inequality
\begin{multline*}
\left\|\bigl((g_{\ep,\delta}^i)'_+(\rho_{\ep,\delta}) \, \rho_{\ep,\delta}-g_{\ep,\delta}^i(\rho_{\ep,\delta}) \bigr) \Div u_{\ep,\delta} \right\|_{L^1L^1} \\
\leq \|(g_{\ep,\delta}^i)'_+(\rho_{\ep,\delta}) \, \rho_{\ep,\delta}-g_{\ep,\delta}^i(\rho_{\ep,\delta}) \|_{L^2L^2}\, \|{\Div u_{\ep,\delta}}\|_{L^2L^2} \leq C
\end{multline*}
and, coming back to Equation~\eqref{eq:renorm_gi}, we get
\begin{equation*}
\partial_t g_{\ep,\delta}^1, \, \partial_t g_{\ep,\delta}^2 \text{ bounded in } L^2(0,T;W^{-1,1}(\TT)).
\end{equation*}
On the other hand, we have
\[
h_{\ep,\delta} = S_{\ep,\delta} + \<\lambda_{\ep,\delta}(\rho_{\ep,\delta})\Div u_{\ep,\delta} - p_{\ep,\delta}\>,
\]
which is bounded in $L^1(0,T;H^1(\TT))$ since
\[
\|\nabla S_{\ep,\delta}\|_{L^2L^2} \leq C \big(\|f\|_{L^2L^{\bar{q}}}+ \|u_{\ep,\delta}\|_{L^2H^1}\big).
\]
Applying Lemma~\ref{lem:compensated_compactness}, we arrive finally at~\eqref{prop:compactness_F_1}--\eqref{prop:compactness_F_2}.
\end{proof}

Thanks to these properties, Equation~\eqref{eq:Psi-0} rewrites
\begin{align*}\labeleq{Psi}
\partial_t \Psi + \Div(\Psi u)
& = \rho_\ep^2\, \ov{\nu_\ep (\rho_\ep) F_\ep } - \ov{\rho_\ep^2 \nu_\ep (\rho_\ep) F_\ep } + \rho_\ep^2\, \ov{ p_\ep(\rho) \nu_\ep(\rho_\ep) } - \ov{ \rho_\ep^2 p_\ep(\rho) \nu_\ep(\rho_\ep) } \\
& = \ov{ F_\ep } \, \bigl(\rho_\ep^2 \,\ov{\nu_\ep(\rho_\ep)} - \ov{\rho_\ep^2 \nu_\ep(\rho_\ep)} \bigr) + \rho_\ep^2\, \ov{ p_\ep(\rho) \nu_\ep(\rho_\ep) } - \ov{ \rho_\ep^2 p_\ep(\rho) \nu_\ep(\rho_\ep) }.
\end{align*}
Let $b_{\ep,\delta}(s) := p_{\ep,\delta}(s) \nu_{\ep,\delta}(s) = \spfrac{p_{\ep,\delta}(s)}{2\mu +\lambda_{\ep,\delta}(s)}$.
As we said before, our constitutive laws depend on the parameter $\delta$ which prevents us from writing directly~\eqref{eq:monotone-press} on~$b_{\ep,\delta}$.
Nevertheless, for all $\delta < \delta^0$, $b_{\ep,\delta^0}(\cdot) \leq b_{\ep,\delta}(\cdot)$ then
\begin{multline*}
\rho_\ep^2 \ \ov{\Bigl(\dfrac{p_\ep(\rho_\ep)}{2\mu +\lambda_\ep(\rho_\ep)} \Bigr)} - \liminf_{\delta \rightarrow 0} \dfrac{\rho_{\ep,\delta}^2\, p_{\ep,\delta}(\rho_{\ep,\delta})}{2\mu +\lambda_{\ep,\delta}(\rho_{\ep,\delta})}\\
\leq \rho_\ep^2 \ \ov{\Bigl(\dfrac{p_\ep(\rho_\ep)}{2\mu +\lambda_\ep(\rho_\ep)} \Bigr)} - \liminf_{\delta \rightarrow 0} \dfrac{\rho_{\ep,\delta}^2\, p_{\ep,\delta^0}(\rho_{\ep,\delta})}{2\mu +\lambda_{\ep,\delta^0}(\rho_{\ep,\delta})}.
\end{multline*}
Since moreover $\gamma > \beta$, the function $b_{\ep,\delta^0}$ is increasing for any $\delta^0 > 0$ and it follows that (see for instance~\cite[Lem.\,3.35]{novotny2004})
\begin{multline*}
\rho_\ep^2 \ \ov{\Bigl(\dfrac{p_\ep(\rho_\ep)}{2\mu +\lambda_\ep(\rho_\ep)} \Bigr)} - \liminf_{\delta \rightarrow 0} \dfrac{\rho_{\ep,\delta}^2\, p_{\ep,\delta}(\rho_{\ep,\delta})}{2\mu +\lambda_{\ep,\delta}(\rho_{\ep,\delta})}\\
\leq \rho_\ep^2 \biggl[ \ov{\Bigl(\dfrac{p_\ep(\rho_\ep)}{2\mu +\lambda_\ep(\rho_\ep)} \Bigr)} - \ov{\Bigl(\dfrac{ p_{\ep,\delta^0}(\rho_{\ep})}{2\mu +\lambda_{\ep,\delta^0}(\rho_{\ep})}\Bigr)}\biggr] = 0,
\end{multline*}
where the last equality holds due to the strong convergence in $L^1\big((0,T)\times \TT\big)$ of $\ov{b_{\ep,\delta^0}(\rho)}$ to $\ov{b_{\ep}(\rho)}$ as $\delta^0 \rightarrow 0$ (this follows from the equi-integrability of the singular pressure and then from the equi-integrability of $b_\delta$).
We get then by integration in space
\begin{align*}\labeleq{Psi2}
\Dt \int_{\TT}{\Psi}~ \leq ~ \int_{\TT}{ |\ov{ F_\ep }|\, \bigl| \rho_\ep^2 \,\ov{\nu_\ep(\rho_\ep)} - \ov{\rho_\ep^2 \nu_\ep(\rho_\ep)} \bigr|}
\end{align*}
with
\begin{align*}
\big|\rho_\ep^2 \nu_{\ep,\delta}(\rho_{\ep,\delta}) - \rho_{\ep,\delta}^2 \nu_{\ep,\delta}(\rho_{\ep,\delta})\big|
& = \nu_{\ep,\delta}(\rho_{\ep,\delta}) \ \big|\rho_{\ep}^2 - \rho_{\ep,\delta}^2 \big|
\leq C \ \big|\rho_{\ep}^2 - \rho_{\ep,\delta}^2 \big|
\end{align*}
thanks to Lemma~\ref{lem:nu_delta}.
Letting $\delta \rightarrow 0$ we get
\begin{align*}
\int_{\TT}{ |\ov{ F_\ep }|\, \bigl| \rho_\ep^2 \,\ov{\nu_\ep(\rho_\ep)} - \ov{\rho_\ep^2 \nu_\ep(\rho_\ep)} \bigr|}
& \leq C \int_{\TT}{ |\ov{ F_\ep }|\, \big|\rho_{\ep}^2 - \ov{\rho_\ep^2} \big|} \\
& \leq C \int_{\TT}{ |\ov{ F_\ep }|\, \big(\ov{\rho_{\ep}^2} - \rho_\ep^2 \big)}
= C \int_{\TT}{|\ov{ F_\ep }|\Psi}.
\end{align*}
Since $|\ov{F_\ep}| \in L^1(0,T;L^\infty(\TT))$ and since $\displaystyle \int_{\TT}{\Psi(0,\cdot)} = 0$ (the initial data does not depend on $\delta$),
we conclude by Gronwall's Lemma that
\[
\Psi = 0 \quad \text{a.e.}
\]
As a consequence, there exists a subsequence (still denoted $\rho_{\ep,\delta}$) such that
\begin{equation*}
\rho_{\ep,\delta} \to \rho_\ep \quad \text{strongly in } L^q\big((0,T) \times \TT\big) \quad\forall q \in [1, +\infty).
\end{equation*}

\subsubsection*{Limit of the singular laws}
From the strong convergence of $\rho_{\ep,\delta}$ and the fact that $\rho_\ep < 1$ a.e. (see ~\eqref{eq:bound_rho_ep}) the following convergences hold
\begin{align*}
p_{\ep,\delta}(\rho_{\ep,\delta}) &\to p_\ep(\rho_\ep) \quad \text{strongly in } L^{1}\big((0,T)\times \TT\big),\\
\sqrt{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}&\to \sqrt{\lambda_\ep(\rho_\ep)} \quad \text{strongly in } L^{\sfrac{2\gamma}{\beta}}\big((0,T)\times \TT\big).
\end{align*}

\subsubsection*{Limit in the weak formulation of the equations.}
We can now pass to the limit in the weak formulation of the mass and momentum equations.
The only delicate term to deal with is the bulk viscosity term $\lambda_{\ep,\delta}(\rho_{\ep,_\delta}) \Div u_{\ep,\delta}$.
We use the strong convergence of $\sqrt{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}$ towards $\sqrt{\lambda_\ep(\rho_\ep)}$ in $L^2\big((0,T)\times \TT\big)$ combined with the weak convergence in $L^2\big((0,T)\times \TT\big)$ of $\Div u_{\ep,\delta}$ to get
\[ \sqrt{\lambda_{\ep,\delta}(\rho_{\ep,\delta})} \Div u_{\ep,\delta} \to \sqrt{\lambda_\ep(\rho_\ep)} \Div u_\ep \quad \text{weakly in } L^{1}((0,T)\times \TT). \]
Since $\sqrt{\lambda_{\ep,\delta}(\rho_{\ep,\delta})} \Div u_{\ep,\delta}$ also converges weakly in $L^2((0,T)\times \TT)$ (from the energy estimate), we deduce that the previous convergence holds in $L^2((0,T)\times \TT)$.
Using once again the strong convergence of $\sqrt{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}$ we obtain the weak convergence of the whole bulk viscosity term
\begin{multline*}
\lambda_{\ep,\delta}(\rho_{\ep,\delta}) \Div u_{\ep,\delta} =\sqrt{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}\sqrt{\lambda_{\ep,\delta}(\rho_{\ep,\delta})} \Div u_{\ep,\delta} \longrightharpoonup \lambda_\ep(\rho_\ep) \Div u_\ep \\ \text{in } L^{1}((0,T)\times \TT).
\end{multline*}
\noindent
Finally, the limit $(\rho_\ep, u_\ep)$ is a global weak solution of the system~
\begin{subnumcasess}{\labeleq{semi-stat-1-ep}}
\partial_t \rho_\ep + \Div (\rho_\ep u_\ep) = 0, \labeleq{semi-stat-1-ep-mass}\\
\nabla p_\ep(\rho_\ep) - \nabla(\lambda_\ep(\rho_\ep) \Div u_\ep) - 2\Div(\mu \D(u_\ep)) + r u_\ep = f, \labeleq{semi-stat-1-ep-mom} \\
0 \leq \rho_\ep < 1 ~ \quad \text{a.e.} ~~ (0,T)\times \TT.
\end{subnumcasess}%
In addition, we have the energy inequality
\begin{multline}\label{eq:energy_ep}
\sup_{t\in [0,T]} \int_{\TT}{H_\ep(\rho_\ep)} + \intTTT{(\tfrac{3}{2}\mu + \lambda_\ep(\rho_\ep))\, |{\Div u_\ep}|^2 }+ \frac{\mu}{2} \intTTT |{\curl u_\ep}|^2 \\
+ r \intTTT |u_\ep|^2 \leq ~ \int_{\TT}{H_\ep(\rho^0_\ep)}
+ \frac{1+C}{2\mu}\, \|f\|^2_{L^2(0,T;L^{\bar{q}}(\TT))},
\end{multline}
where
\[
H_\ep(\rho) = \dfrac{\ep}{\gamma -1}\, \dfrac{\rho^\gamma}{(1-\rho)^{\gamma-1}}.
\]

\Subsection{Congestion limit $\ep \rightarrow 0$}

\subsubsection{Uniform estimates in $\ep$}
From the energy estimate~\eqref{eq:energy_ep}, we deduce the controls of different quantities uniformly with respect to the parameter $\ep$
\begin{align*}
& \big(H_\ep(\rho_\ep)\big)_\ep ~ \text{is bounded in} ~L^\infty(0,T;L^1(\TT)) \labeleq{bound_internalernergy}, \\
& \big(u_\ep\big)_\ep ~ \text{is bounded in} ~ L^2\big(0,T;(H^1(\TT))^2\big),
\labeleq{bound_u}\\
& \big(\sqrt{\lambda_\ep(\rho_\ep)}\,\Div u_\ep\big)_\ep ~ \text{is bounded in} ~ L^2\big((0,T)\times\TT\big).
\end{align*}

\Subsubsection*{Control of the pressure}
\begin{lem}
Let $(\rho_\ep,u_\ep)$ be a global weak solution to~\eqref{eq:semi-stat-0} with $\gamma > \beta > 1$.
Then, there exists a constant $C > 0$ independent of $\ep$ such that
\begin{equation*}
\|p_\ep(\rho_\ep)\|_{L^1((0,T)\times \TT)} \leq C.
\end{equation*}
\end{lem}

\begin{proof} 
Let us consider in the weak formulation of the momentum equation the test function
\[
\psi = \nabla \Delta^{-1}(\rho_{\ep} - \<\rho_{\ep}\>).
\]
The resulting equation writes
\begin{multline*}
\intTTT{p_{\ep}(\rho_{\ep}) \big(\rho_{\ep} - \<\rho_{\ep}\>\big) } \\[-5pt]
= \intTTT{\big(2\mu+ \lambda_{\ep}(\rho_{\ep})\big)\Div(u_{\ep}) \big(\rho_{\ep} - \<\rho_{\ep}\>\big) }
- \intTTT{(f-ru_\ep) \cdot \psi}.
\end{multline*} 
Using the controls resulting from the energy inequality and the maximal bound on the density $\rho_\ep$, we get then
\begin{align*}
\intTTT p_{\ep}(\rho_{\ep}) &\big(\rho_{\ep} - \<\rho_{\ep}\>\big)\\[-5pt]
& ~ \leq C(\|\psi\|_{H^1}) + \intTTT{\sqrt{(2\mu+\lambda_{\ep}(\rho_{\ep}))} \sqrt{(2\mu+\lambda_{\ep}(\rho_{\ep}))}\,|{\Div(u_{\ep})}| } \\
& ~ \leq C + \dfrac{1}{2\eta}\intTTT{(2\mu+\lambda_{\ep}(\rho_{\ep}))\,|{\Div(u_{\ep})}|^2}
+ \dfrac{\eta}{2} \intTTT{(2\mu +\lambda_{\ep}(\rho_{\ep}))} \\
& ~ \leq C + \dfrac{\eta}{2}\intTTT{\lambda_{\ep}(\rho_{\ep})}
\end{align*}
for some $\eta>0$ determined below.
The integrals of the singular terms can be split into two parts depending on the value of $\rho_\ep$.
Recall that $\<\rho_\ep\> = \<\rho_\ep^0\>$ is far from $1$ uniformly in $\ep$ (\cf \eqref{hyp:bound-rho0}), and the functions $p_\ep(\rho_\ep)$, $\lambda_{\ep}(\rho_\ep)$ are bounded on the domain $ \{ \rho_{\ep} < {(1+\<\rho_{\ep}\>)}/{2}\}$.
Therefore:
\begin{align*}
\intTTT p_{\ep}(\rho_{\ep})& \big(\rho_{\ep} - \<\rho_{\ep}\>\big) \,\mathbf{1}_{\{ \rho_{\ep} \geq \psfrac{1+\langle\rho_{\ep}\rangle}{2} \} } \\
&\leq C - \intTTT{p_{\ep}(\rho_{\ep}) \big(\rho_{\ep} - \<\rho_{\ep}\>\big) \,\mathbf{1}_{\{ \rho_{\ep} < \psfrac{1+\langle\rho_{\ep}\rangle}{2} \} } }\\
&\hspace*{.75cm} + \dfrac{1}{2}\intTTT{\lambda_{\ep}(\rho_{\ep}) \,\mathbf{1}_{\{ \rho_{\ep} < \psfrac{1+\langle\rho_{\ep}\rangle}{2}\} } } + \dfrac{\eta}{2}\intTTT{\lambda_{\ep}(\rho_{\ep}) \,\mathbf{1}_{\{ \rho_{\ep} \geq \psfrac{1+\langle\rho_{\ep}\rangle}{2}\} } } \\
& ~ \leq C + \dfrac{\eta}{2}\intTTT{\lambda_{\ep}(\rho_{\ep}) \,\mathbf{1}_{\{ \rho_{\ep} \geq \psfrac{1+\langle\rho_{\ep}\rangle}{2}\} } }.
\end{align*}
Since $\beta < \gamma$, we have in addition that
\begin{align*}
\lambda_{\ep}(\rho_{\ep}) \,\mathbf{1}_{\{ \rho_{\ep} \geq \psfrac{1+\langle\rho_{\ep}\rangle}{2}\} }
& = \Bigl(\dfrac{1+\<\rho_{\ep}\>}{1-\<\rho_{\ep}\>}\Bigr)^{\beta-\gamma}\ep \Bigl(\dfrac{\rho_{\ep,\delta}}{1-\rho_{\ep,\delta}}\Bigr)^{\gamma} \,\mathbf{1}_{\{ \rho_{\ep} \geq \psfrac{1+\langle\rho_{\ep}\rangle}{2}\}} \\
& \leq p_{\ep}(\rho_{\ep}) \,\mathbf{1}_{\{ \rho_{\ep} \geq \psfrac{1+\langle\rho_{\ep}\rangle}{2}\} }.
\end{align*}
Now, with $\eta = \psfrac{1-M^0}{2}$,
\[
\dfrac{1-\<\rho_{\ep}\>}{2} \intTTT{p_{\ep}(\rho_{\ep})\,\mathbf{1}_{\{ \rho_{\ep} \geq \psfrac{1+\langle\rho_{\ep}\rangle}{2} \} } }\leq C + \dfrac{1-M^0}{4} \intTTT{p_{\ep}(\rho_{\ep}) \,\mathbf{1}_{\{ \rho_{\ep} \geq \psfrac{1+\langle\rho_{\ep}\rangle}{2}\} } },
\]
and recalling that we have
\[
\dfrac{1-\<\rho_{\ep}\>}{2} = \dfrac{1-\<\rho_{\ep}^0\>}{2} \geq \dfrac{1-M^0}{2}
\]
due our assumption~\eqref{hyp:bound-rho0}${}_2$,
we control the integral of the pressure on $\{ \rho_{\ep} \geq \psfrac{1+\langle\rho_{\ep}\rangle}{2}\}$
\[
\dfrac{1-M^0}{4}
\intTTT{p_{\ep}(\rho_{\ep})\,\mathbf{1}_{\{ \rho_{\ep} \geq \psfrac{1+\langle\rho_{\ep}\rangle}{2} \} } } \leq C.
\]
Finally, since the pressure is bounded on $\{ \rho_{\ep} < \psfrac{1+\langle\rho_{\ep}\rangle}{2}\}$ (far from the singularity), we get
\begin{equation*}\labeleq{bound-pep-L1}
\big(p_{\ep}(\rho_{\ep})\big)_\ep~ \text{bounded in} ~L^1\big((0,T)\times \TT \big)
\end{equation*}
which ends the proof.
\end{proof} 

\subsubsection{Limit $\ep\rightarrow 0$}

With the previous uniform bounds we deduce that there exist a limit density $\rho$ and a limit velocity $u$ such that
\begin{align*}
\rho_{\ep} &\longrightharpoonup \rho \quad \text{weakly-* in} ~ L^\infty\big((0,T)\times \TT \big), \quad \rho(t,x) \leq 1 ~\text{a.e.},\\
u_{\ep} &\longrightharpoonup u \quad \text{weakly in } L^2(0,T;(H^1(\TT))^3),
\end{align*}
and we pass to the limit in the nonlinear term $\rho_\ep u_\ep$ thanks to the compensated-compactness lemma~\ref{lem:compensated_compactness}:
\begin{equation*}
\rho_{\ep} u_{\ep} \longrightharpoonup \rho u \quad \text{weakly-* in } L^\infty(0,T;L^6(\TT)).
\end{equation*}
The boundedness of $(p_\ep(\rho_\ep))_\ep$ in $L^1\big((0,T)\times \TT\big)$ yields
\begin{equation*}
p_\ep(\rho_\ep) \to p \quad \text{weakly in } \mathcal{M}_+\big((0,T) \times \TT\big).
\end{equation*}
On the other hand, we have
\begin{align*}
&\lambda_\ep(\rho_\ep) \leq C \ep \,\mathbf{1}_{\{\rho_\ep \leq M^0\}}
+ C \ep^{1- \sfrac{\beta}{\gamma}} \big(p_\ep(\rho_\ep)\big)^{\sfrac{\beta}{\gamma}} \,\mathbf{1}_{\{\rho_\ep > M^0\}}
\end{align*}
with $\beta < \gamma$, so that $\lambda_\ep(\rho_\ep)$ converges strongly to $0$ in $L^{\sfrac{\gamma}{\beta}}((0,T)\times \TT)$.
Hence, $\sqrt{\lambda_{\ep}(\rho_\ep)}$ converges strongly to $0$ in $L^2\big((0,T)\times \TT\big)$ and
\begin{equation*}
\lambda_\ep(\rho_\ep) \Div u_\ep \longrightharpoonup \ov{\lambda_\ep \Div u} = 0 \quad \text{in}~ \mathcal{D}'.
\end{equation*}
Finally, it remains for us to prove that $\spt p \subset \{\rho=1\}$.
First, let us assume that there exists a set $A \subset \{\rho < 1 \}$ with non-zero measure such that $\lim_{\ep \rightarrow 0} \rho_\ep(t,x) = 1$ for $(t,x) \in A$.
Then, since $\rho_\ep$ converges weakly to $\rho$ in $L^q((0,T)\times \TT)$ for all $q \in [1,+\infty)$, we have
\[
\intTTT{\rho_\ep \,\mathbf{1}_A} \to \intTTT{\rho \,\mathbf{1}_A},
\]
which leads to a contradiction since
\[
\intTTT{\rho_\ep \,\mathbf{1}_A} \to |A|
\quad \text{while} \quad \intTTT{\rho \,\mathbf{1}_A} < |A|.
\]
Hence, there does not exist such $A$.
Let us assume now that there exists a set $B \subset \{\rho <1 \}$ such that $p(t,x) > 0$ for $(t,x) \in B$. Then, the weak convergence of $p_\ep(\rho_\ep)$ to~$p$ gives, for any $\Psi \in \mathcal{C}^\infty_c(\overline{B})$, $\Psi \geq 0$,
\[
\intTTT{p_\ep(\rho_\ep) \Psi} \to \<p,\Psi\> > 0.
\]
In view of the definition of $p_\ep$, we infer that there exists a subset $A \subset B$ with non-zero measure such that $\lim_{\ep \rightarrow 0} \rho_\ep(t,x) = 1$ for $(t,x) \in A$ (otherwise we would have $p_\ep(\rho_\ep) \rightarrow 0$ a.e. on $B$).
By the previous argument such a set $A$ does not exist, and we conclude that $|B|=0$.
Therefore
\begin{equation}\label{eq:spt-p}
p = 0 \quad \text{a.e. on } \{\rho < 1 \}.
\end{equation}
The limit system reads in this case
\begin{subnumcasess}{\labeleq{pb_limit}}
\partial_t \rho + \Div(\rho u) = 0, \\
\nabla p -2\Div(\mu \D(u)) + ru = f, \\
0 \leq \rho \leq 1, \quad \spt p \subset \{ \rho = 1\}, \quad p \geq 0,
\end{subnumcasess}%
where the memory effects are never activated.
\begin{remark}{\label{rmk:cvg-rhoep}}
Note that, unlike the previous limit $\delta\rightarrow 0$, we do not ensure the strong convergence of $(\rho_\ep)_\ep$ to $\rho$.
The problem relies in the lack of uniform estimates (see Equation \eqref{eq:inegal-dep-ep}) which prevents the derivation of the weak compactness properties on~$F_\ep$.
Nevertheless, as explained before, we are able to identify the weak limit of the nonlinear term $\rho_\ep u_\ep$, and to pass to the limit in the mass and momentum equations without the strong convergence of $(\rho_\ep)_\ep$.
\end{remark}

\section{Dominant bulk viscosity regime \texorpdfstring{$1<\gamma\leq \beta$}{1gb}}{\label{sec:bulk}}

Let us now consider the case where $\beta \geq \gamma$.
If the approach proceeds formally in the same way as before (regularization of the system by truncation of the singular laws and study of the behavior as $\ep \rightarrow 0$), we have here to adapt the arguments to get the appropriate uniform controls in $\delta$ and $\ep$.
For that purpose, we shall distinguish in the estimates three cases: $\gamma < \beta-1$, $\gamma = \beta -1$ and $\beta-1 \leq \gamma < \beta$ that correspond to the sub-cases presented in Theorem~\ref{thm-limit}.
In these three cases, we are not able to control the bulk viscosity coefficient $\lambda_\ep$ from the pressure $p_\ep$.
Nevertheless, we will see that Equation~\eqref{eq:expr_lambdadivu} enables to pass to the limit $\ep \rightarrow 0$ in the two cases $\gamma \in (\beta-1, \beta]$ (no memory effects at the limit), $\gamma \leq \beta-1$ (memory effects).

\Subsection{Case $\beta-1 < \gamma \leq \beta$, $\gamma > 1$}

\subsubsection{Existence of weak solutions at $\ep$ fixed}
We consider the same regularized system~\eqref{eq:semi-stat-delta} with the truncated pressure~\eqref{eq:p_delta} and bulk viscosity~\eqref{eq:lambda_delta} as in the previous section.
Recall that we ensure from Lemma~\ref{lem:nu_delta} the following properties on $\nu$:
\begin{equation*}
\exists\, C_1, C_2 > 0 \quad \text{s.t.} \quad \|\nu_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^\infty_{t,x}}\leq C_1,
\quad \int_{\TT}{\nu_{\ep,\delta}(\rho_{\ep,\delta})} \geq C_2.
\end{equation*}
\begin{lem}\label{lem:controls-1} Let $(\rho_{\varepsilon, \delta}, u_{\varepsilon,\delta})$ be a global weak solution of the compressible Brinkman system~\eqref{eq:semi-stat-delta} with $\gamma \leq \beta$. Then,
\[
\|\Lambda_{\varepsilon,\delta}(\rho_{\ep,\delta})\|_{L^1((0,T)\times \TT)}
+\|p_{\varepsilon,\delta}(\rho_{\ep,\delta})\|_{L^1((0,T)\times \TT)}
+\|\lambda_{\varepsilon,\delta}(\rho_{\ep,\delta}) \Div u_{\varepsilon,\delta}\|_{L^1((0,T)\times \TT)}
\leq C,
\]
where $C$ does not depend on $\delta$ or $\ep$.
\end{lem}
\begin{proof} Starting again from the equation on $\Lambda_{\ep,\delta}$~\eqref{eq:Lambda_delta}, using~\eqref{eq:lambda_divu} to replace $\lambda_{\ep,\delta} \Div u_{\ep,\delta}$, and integrating in space we have
\begin{multline*}{\labeleq{control-Lambda}}
\dfrac{\di}{\di t} \int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta})} + \int_{\TT}{p_{\ep,\delta}(\rho_{\ep,\delta})}
\leq C \Big(\|S_{\ep,\delta}\|_{L^1} + \int_{\TT}{p_{\ep,\delta}(\rho_{\ep,\delta}) \nu_{\ep,\delta}(\rho_{\ep,\delta}) \,\mathbf{1}_{\{\rho_{\ep,\delta} \leq M^0\}}}\\
\shoveright{+ \int_{\TT}{p_{\ep,\delta}(\rho_{\ep,\delta}) \nu_{\ep,\delta}(\rho_{\ep,\delta}) \,\mathbf{1}_{\{M^0 < \rho_{\ep,\delta} \}}} \Big)}\\
\leq C + C \int_{\TT}{\,\mathbf{1}_{\{\rho_{\ep,\delta} \leq M^0\}}}
+ \int_{\TT}{ \dfrac{C}{(1-\rho_{\ep,\delta})^{\gamma - \beta}}\,\mathbf{1}_{\{M^0 < \rho_{\ep,\delta} \leq 1-\delta \}}}
+ \int_{\TT}{ \dfrac{C}{\delta^{\gamma - \beta}}\,\mathbf{1}_{\{\rho_{\ep,\delta} > 1-\delta\}} },
\end{multline*}
for some constant $C$ independent of $\delta, \ep$.
Now, since $\gamma < \beta$, we can bound uniformly in $\delta,\ep$ the right-hand side and therefore
\begin{align*}
\big(\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta &\text{ is bounded in } L^\infty\big(0,T;L^1(\TT)\big),\\
\big(p_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta &\text{ is bounded in } L^1\big((0,T)\times \TT\big).
\end{align*}
We have as a byproduct (see Lemma~\ref{lambdadivu})
\begin{equation*}\labeleq{lambda-delta-2}
\big(\lambda_{\ep,\delta}(\rho_{\ep,\delta}) \Div u_{\ep,\delta}\big)_\delta ~ \text{bounded in}~ L^1\big((0,T)\times \TT\big).\qedhere
\end{equation*}
\end{proof}

\subsubsection*{Uniform upper bound on the density}
Proposition~\ref{prop:bound_rhodelta} still holds in the case $\gamma < \beta$: there exists $C$ which does not depend on $\delta$ such that
\begin{equation}\label{eq:bound_rhodelta_2}
\|\rho_{\ep,\delta}\|_{L^\infty_{t,x}} \leq \bar{\rho} < \infty.
\end{equation}

\begin{lem}\label{lem:bound-rhod-Lambda}
For any $t \in [0,T]$, we have
\begin{equation}\label{eq:est-rho-delta-2}
\meas \big\{ x\in \TT, \ \rho_{\ep,\delta}(t,x) \geq 1-\delta \big\} \leq C(\ep) \ \delta^{\beta-1}.
\end{equation}
\end{lem}

\begin{proof} 
Since initially we assume~\eqref{hyp:energy-t0-C2}, we can mimic the proof of~\eqref{eq:est-rhodelta} by using $\Lambda_{\ep,\delta}$ instead of $H_{\ep,\delta}$:
\begin{align*}
C &\geq \sup_{t\in [0,T]} \int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta})}
 \geq \int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta}) \,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}} \\
& \geq C\int_{\TT}{\dfrac{\ep}{\delta^\beta}\,\bigl[
(\rho_{\varepsilon,\delta}^{\beta-1} - (1-\delta)^{\beta-1})
+ (1-\delta)^{\beta-1}\delta \bigr] \rho_{\ep,\delta}
\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}} \\
& \geq C\int_{\TT}{\dfrac{\ep}{\delta^{\beta-1} }\,\bigl[
(1-\delta)^{\beta} \bigr]
\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}}.\qedhere
\end{align*}
\end{proof}

As explained in the previous section, passing to the limit in the weak formulation of the equations requires additional estimates on the singular laws.

\begin{lem}\label{lem:controls-2}Let $(\rho_{\varepsilon, \delta}, u_{\varepsilon,\delta})$ be a global weak solution of the compressible Brinkman system~\eqref{eq:semi-stat-delta} with $\gamma \in (\beta-1,\beta]$.
If initially~\eqref{hyp:energy-t0-C2} is satisfied, then there exist two constants $C^1>0$ independent of $\delta,\ep$, and $C^2_\ep> 0$ independent of $\delta$, such that
\begin{align*}
\|\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^\infty(0,T; L^2(\TT))}
&\leq C^1,\\
\|p_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^{1+\psfrac{\beta-1}{\gamma}}((0,T)\times \TT)} + \|\lambda_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^{\psfrac{\beta + \gamma -1}{\beta}}((0,T)\times \TT)}
&\leq C^2_\ep.
\end{align*}
\end{lem}

\begin{proof}
Taking $b(s) = \big(\Lambda_{\ep,\delta}\big)(s)^2$ in \eqref{eq:renormalized-delta} we get
\begin{equation}\label{eq:eta-Lambda}
\begin{split}
\partial_t\big(\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)^2& + \Div\big(\big(\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)^2 u_{\ep,\delta} \big)\\
& = - \Big[\rho_{\ep,\delta} (\Lambda_{\ep,\delta})'_+(\rho_{\ep,\delta})\Lambda_{\ep,\delta}(\rho_{\ep,\delta})
+ \Lambda_{\ep,\delta}(\rho_{\ep,\delta}) \lambda_{\ep,\delta}(\rho_{\ep,\delta})
\Big] \Div u_{\ep,\delta} \\
& = - \Bigl[\frac{\rho_{\ep,\delta} (\Lambda_{\ep,\delta})'_+(\rho_{\ep,\delta})}{\lambda_{\ep,\delta}(\rho_{\ep,\delta})} + 1\Bigr]
\Lambda_{\ep,\delta}(\rho_{\ep,\delta})
\lambda_{\ep,\delta}(\rho_{\ep,\delta})\Div u_{\varepsilon, \delta},
\end{split}
\end{equation}
where
\[
(\Lambda_{\ep,\delta})'_+(\rho) =
\begin{cases}
\dfrac{\ep}{(1-\rho)^\beta}\, \dfrac{(\beta-\rho)\rho^{\beta-1}}{\beta-1} \quad&\text{if } \rho < 1-\delta, \\[10pt]
\dfrac{\ep}{\delta^\beta}\, \dfrac{\beta\rho^{\beta-1} - (1-\delta)^\beta}{\beta-1} \quad&\text{if } \rho \geq 1-\delta,
\end{cases}
\]
which is such that
\[
0 \leq \dfrac{\rho_{\ep,\delta} \ (\Lambda_{\ep,\delta})'_+(\rho_{\ep,\delta})}{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}
\leq \dfrac{\beta}{\beta-1}.
\]
Recall now that
\begin{multline*}
\lambda_{\ep,\delta}(\rho_{\ep,\delta})\Div u_{\ep,\delta}
= p_{\ep,\delta} (\rho_{\ep,\delta})
+ S
- 2\mu \, \Div u_{\ep,\delta} \\
+ \dfrac{1}{|\TT|} \int_{\TT}{\big(\lambda_{\ep,\delta}(\rho_{\ep,\delta})\Div u_{\ep,\delta} - p_{\ep,\delta} (\rho_{\ep,\delta}) \big)}
\end{multline*}
that we can replace in~\eqref{eq:eta-Lambda}.
Integrating next in space, we obtain (recall that$\beta > 1$)
\begin{align*}\labeleq{eta-2}
\Dt \int_{\TT}&\big(\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)^{2}
+ \int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\, p_{\ep,\delta} (\rho_{\ep,\delta})} \\
& \le
\dfrac{2\beta -1}{\beta-1}\int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\, |S_{\ep,\delta}}|
+ 2\mu\,\dfrac{2\beta -1}{\beta-1} \int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\,
| {\Div u_{\ep,\delta}}|} \\
&
\hspace*{1cm}+ \dfrac{2\beta -1}{\beta-1} \dfrac{1}{|\TT|}\left(\int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta})}\right) \biggl|\left(\int_{\TT}{\big(\lambda_{\ep,\delta}(\rho_{\ep,\delta})\Div u_{\ep,\delta} - p_{\ep,\delta} (\rho_{\ep,\delta})\big)} \right) \biggr| \\
& = I_1 + I_2 + I_3.
\end{align*}
By a Cauchy-Schwarz inequality, we have
\begin{align*}
I_1 + I_2
& \leq C (\|S_{\ep,\delta}\|_{L^2}^2 + \|{\Div u_{\ep,\delta}}\|_{L^2}^2) + \int_{\TT}{\big(\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)^{2}}
\end{align*}
while, due to the previous estimates,
\begin{align*}
\int_0^T I_3
& \leq C \|\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^\infty L^1}
\big(\|\lambda_{\ep,\delta}(\rho_{\ep,\delta})\Div u_{\ep,\delta}\|_{L^1L^1} + \|p_{\ep,\delta} (\rho_{\ep,\delta})\|_{L^1L^1} \big).
\end{align*}
We conclude then with a Gronwall inequality (recall that initially we assume~\eqref{hyp:energy-t0-C2}) and get
\begin{equation*}
\sup_{t \in [0,T]}\int_{\TT}{\big(\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)^{2}}
+ \int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta}) p_{\ep,\delta} (\rho_{\ep,\delta})}
\leq C + C\int_{\TT}{\big(\Lambda_{\ep,\delta}(\rho_{\ep}^0)\big)^{2}}.
\end{equation*}
With the control of $\Lambda_{\ep,\delta}(\rho_{\ep,\delta}) p_{\ep,\delta}(\rho_{\ep,\delta})$ in $L^1\big((0,T)\times \TT\big)$, we deduce that
\[
\Bigl\|\dfrac{1}{(1-\rho_{\ep,\delta})^{\gamma + \beta-1} }\,\mathbf{1}_{\{ \rho_{\ep,\delta} \leq 1-\delta \}} \Bigr\|_{L^1} \leq C(\ep),
\]
and
\[
\Bigl\|\dfrac{\rho_{\ep,\delta}^{\gamma + \beta -1}}{\delta^{\gamma+ \beta-1}} \,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}} \Bigr\|_{L^1} \leq C(\ep),
\]
since (see Lemma~\ref{lem:bound-rhod-Lambda} and~\eqref{eq:bound_rhodelta_2})
\[
\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\,\mathbf{1}_{\{ \rho_{\ep,\delta} > 1-\delta \}}
\geq \dfrac{\ep (1-\delta)^\beta}{\delta^{\beta-1}}
\quad \text{and} \quad
\|\rho_{\ep,\delta}\|_{L^\infty} \leq C.
\]
Hence,
\begin{gather}\labeleq{pdelta2}
\big(p_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta ~ \text{is bounded in} ~ L^{1 + \psfrac{\beta-1}{\gamma}}((0,T)\times \TT),\notag\\
\label{eq:lambdadelta2}
\big(\lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta ~ \text{is bounded in} ~ L^{\psfrac{\gamma + \beta - 1}{\beta}}((0,T)\times \TT).
\end{gather}
We emphasize the fact that these controls are not uniform in $\ep$.
This is due to the fact that the pressure $p_\ep$ and the bulk viscosity $\lambda_\ep$ are ``more singular'' than $\Lambda_\ep$ close to $1$ (since $\gamma, \beta > \beta-1$).
\end{proof}

\begin{remark}
We assumed in this case that $\gamma > 1$ which enables, thanks to~\eqref{eq:lambdadelta2}, to bound $\lambda_{\ep,\delta}(\rho_{\ep,\delta})$ in $L^1((0,T) \times \TT)$.
\end{remark}

\subsubsection*{Limit $\delta \rightarrow 0$}
We can pass to the limit in the weak formulation of the equation except in the non-linear terms.
As in the previous section, we need to prove the strong convergence of $\rho_{\ep,\delta}$.
Note first that the results of Proposition~\ref{prop:compactness_F} still hold:
\begin{align*}\labeleq{prop:compactness_F_1-2}
\ov{\Bigl(\dfrac{\rho^2 F_\ep}{2\mu +\lambda_\ep(\rho)} \Bigr)} &= \ov{\Bigl(\dfrac{\rho^2}{2\mu +\lambda_\ep(\rho)} \Bigr)}~ \ov{F_\ep}\qquad \text{in } \mathcal{D}',\\\labeleq{prop:compactness_F_2-2}
\ov{\Bigl(\dfrac{F_\ep}{2\mu +\lambda_\ep(\rho)} \Bigr)} &= \ov{\Bigl(\dfrac{1}{2\mu +\lambda_\ep(\rho)} \Bigr)}~ \ov{F_\ep}\qquad \text{in } \mathcal{D}'.
\end{align*}
Hence, for $\Psi = \ov{\rho_\ep^2} - \rho_\ep^2 \geq 0$,
\[\labeleq{Psi-2}
\partial_t \Psi + \Div(\Psi u)
= \ov{ F_\ep } \, \bigl(\rho_\ep^2 \,\ov{\nu_\ep(\rho_\ep)} - \ov{\rho_\ep^2 \nu_\ep(\rho_\ep)} \bigr) 
 + \rho_\ep^2\, \ov{ p_\ep(\rho_\ep) \nu_\ep(\rho_\ep) } - \ov{ \rho_\ep^2 p_\ep(\rho_\ep) \nu_\ep(\rho_\ep) }.
\]
where, in this case, the function $s \mapsto b_{\ep,\delta}(s) = p_{\ep,\delta}(s) \nu_{\ep,\delta}(s)$ is non-monotone.
In fact it is not a problem, because $b_{\ep,\delta}$ is now bounded and we can then treat it as the first part of the right-hand side.
We obtain similarly:
\begin{equation*}
\dfrac{\di}{\di t} \int_{\TT}{\Psi(t,\cdot)}
\leq C \int_{\TT}{(|\ov{F_\ep}|+1)\Psi(t,\cdot)},
\end{equation*}
with $|\ov{F_\ep}|\in L^1(0,T;L^\infty(\TT))$ which yields again $\Psi = 0$ by Gronwall's inequality.
The strong convergence of the density is thus preserved in this case.
In addition, due to~\eqref{eq:est-rho-delta-2}, we ensure that the limit density satisfies
\begin{equation}\label{eq:rho-ep-max}
0\leq \rho_{\ep} < 1 \quad \text{a.e.}
\end{equation}
We can then pass to the limit in all the terms of the equations.
In particular, for the bulk viscosity term, the strong convergence of $\rho_{\ep,\delta}$ and~\eqref{eq:rho-ep-max} imply that $\lambda_{\ep,\delta}(\rho_{\ep,\delta})$ converges strongly to $\lambda_\ep(\rho_\ep)$ in $L^{\psfrac{\gamma + \beta -1}{\beta}}((0,T)\times \TT)$.
As in the previous section we deduce that
\[
\lambda_{\ep,\delta}(\rho_{\ep,\delta})\Div u_{\ep,\delta}
= \sqrt{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}\,\sqrt{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}\,\Div u_{\ep,\delta}
\]
converges weakly to $\lambda_\ep(\rho_\ep)\Div u_\ep$ in $L^q((0,T)\times \TT)$ for some $q>1$.

Finally, we pass to the limit in~\eqref{eq:Lambda_delta} (recall $\Lambda_{\ep,\delta}(\rho_{\ep,\delta})$ is bounded in $L^\infty(0,T;L^2(\TT))$ so that $\Lambda_{\ep,\delta}(\rho_{\ep,\delta}) u_{\ep,\delta}$ is bounded in $L^2(0,T;(L^{3/2}(\TT))^3)$):
\begin{equation}\label{eq:Lambda_ep-2}
\partial_t \Lambda_\ep(\rho_\ep) + \Div\big(\Lambda_\ep(\rho_\ep) u_\ep \big) = -\lambda_\ep(\rho_\ep)\Div u_\ep.
\end{equation}

\subsubsection{Limit $\ep \rightarrow 0$}
At this stage, we control uniformly
\begin{align*}
&\big(\rho_\ep \big)_\ep \quad \text{in } L^\infty((0,T)\times \TT), \\
&\big(u_\ep \big)_\ep \quad \text{in } L^2(0,T;(H^1(\TT))^3), \\
&\big(\Lambda_\ep(\rho_\ep)\big)_\ep \quad \text{in } L^\infty(0,T; L^2(\TT)),\\
&\big(p_\ep(\rho_\ep)\big)_\ep \quad \text{in } L^1((0,T)\times \TT),\\
&\big(\lambda_\ep(\rho_\ep) \Div u_\ep\big)_\ep \quad \text{in } L^1((0,T)\times \TT).
\end{align*}
Therefore
\begin{equation*}
p_\ep(\rho_\ep) \to p \quad \text{in}~\mathcal{M}_+((0,T)\times \TT)
\end{equation*}
and
\begin{equation*}
\lambda_\ep(\rho_\ep) \Div u_\ep\to \Pi \quad \text{in}~\mathcal{M}((0,T)\times \TT).
\end{equation*}
On the other hand
\begin{equation*}
\Lambda_\ep(\rho_\ep) \leq C \ep \,\mathbf{1}_{\{\rho_\ep \leq M^0\}}
+ C \ep^{1- \psfrac{\beta-1}{\gamma}} \big(p_\ep(\rho_\ep) \big)^{\psfrac{\beta-1}{\gamma}} \,\mathbf{1}_{\{\rho_\ep > M^0\}},
\end{equation*}
with $\gamma > \beta-1$, so that $\Lambda_\ep(\rho_\ep)$ converges strongly to $0$ in $L^{\sfrac{\gamma}{\beta-1}}((0,T)\times \TT)$ (and thus in $L^2(0,T;L^{\sfrac{6}{5}}(\TT))$).
Passing to the limit in~\eqref{eq:Lambda_ep-2}, we have then
\[
\Pi = 0.
\]
As a consequence, we conclude that the triple $(\rho, u, p)$ is a global weak solution of
\begin{subnumcasess}{\labeleq{xxx3}}
\partial_t \rho + \Div (\rho u) = 0, \\
\nabla p - 2 \Div(\mu \D(u)) + ru = f, \\
0 \leq \rho \leq 1, ~\spt p \subset \{\rho = 1\},~ p \geq 0.
\end{subnumcasess}%

\subsection{Case $1 <\gamma < \beta - 1$}

In this case, we expect the activation of the memory effects in the limit $\ep \rightarrow 0$.

\subsubsection{Existence of weak solutions at $\ep$ fixed}
First we observe that Lemma~\ref{lem:controls-1} still holds in this case, so that\vspace*{-3pt}
\begin{equation*}
\big(\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta \text{ is bounded in } L^\infty\big(0,T;L^1(\TT)\big)
\end{equation*}
and, since $\gamma < \beta-1$,\vspace*{-3pt}
\[
\big(p_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta \text{ is bounded in } L^\infty(0,T;L^{\psfrac{\beta-1}{\gamma}}(\TT)).
\]
With this control of the pressure and the energy estimate, we deduce from Lemma~\ref{lambdadivu} that\vspace*{-3pt}
\begin{equation*}\labeleq{lambdadivu-3}
\big(\lambda_{\ep,\delta}(\rho_{\ep,\delta}) \Div u_{\ep,\delta}\big)_\delta ~ \text{is bounded in}~
L^2\big(0,T;L^{\min\{2,\psfrac{\beta-1}{\gamma}\}}(\TT)\big).
\end{equation*}

\subsubsection*{Uniform upper bound on the density}
The same arguments as before show that\vspace*{-3pt}
\[
\|\rho_{\ep,\delta}\|_{L^\infty_{t,x}} \leq \bar{\rho} < \infty,
\]
and\vspace*{-3pt}
\[
\meas \big\{ x\in \TT, \ \rho_{\ep,\delta}(t,x) \geq 1-\delta \big\} \leq C(\ep) \ \delta^{\beta-1}.
\]

\Subsubsection*{Additional integrability on $\Lambda_{\ep,\delta}$}
\begin{lem}\label{lem:controls-3}Let $(\rho_{\varepsilon, \delta}, u_{\varepsilon,\delta})$ be a global weak solution of the compressible Brinkman system~\eqref{eq:semi-stat-delta} with $1 <\gamma < \beta-1$.
If initially $\|\Lambda_\ep(\rho^0_\ep)\|_{L^2} \leq C$, then there exist $C^1>0$ independent of $\delta,\ep$, and $C^2_\ep> 0$ independent of $\delta$, such that
\begin{align*}
\|\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^\infty L^2} + \|p_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^\infty L^{\sfrac{2(\beta-1)}{\gamma}}}
&\leq C^1,\\
\|\lambda_{\ep,\delta}(\rho_{\ep,\delta})\|_{L^\infty L^{\sfrac{2(\beta-1)}{\beta}}}
&\leq C^2_\ep.
\end{align*}
\end{lem}

\begin{proof} We can reproduce exactly the same estimate on $\big(\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)^{2}$ by taking $b(s) = \big(\Lambda_{\ep,\delta}\big)(s)^2$ in \eqref{eq:renormalized-delta}.
We obtain\vspace*{-3pt}
\begin{equation*}
\sup_{t \in [0,T]}\int_{\TT}{\big(\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)^{2}}
+ \int_{\TT}{\Lambda_{\ep,\delta}(\rho_{\ep,\delta}) p_{\ep,\delta} (\rho_{\ep,\delta})}
\leq C + C\int_{\TT}{\big(\Lambda_{\ep,\delta}(\rho_{\ep}^0)\big)^{2}}.
\end{equation*}
Hence\vspace*{-3pt}
\begin{equation*}\labeleq{Lambdadelta3}
\big(\Lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta ~ \text{is bounded in} ~ L^\infty(0,T;L^2(\TT)),
\end{equation*}
and consequently\vspace*{-3pt}
\begin{gather*}\labeleq{pdelta3}
\big(p_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta ~ \text{is bounded in} ~ L^\infty(0,T;L^{\sfrac{2(\beta-1)}{\gamma}}(\TT)),
\\
\labeleq{lambdadelta3}
\big(\lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta ~ \text{is bounded in} ~ L^\infty(0,T;L^{\sfrac{2(\beta-1)}{\beta}}(\TT)).
\end{gather*}
This latter control on $\lambda_\ep(\rho_\ep)$ being not uniform in $\ep$.
\end{proof}

\begin{remark}\label{rk:control-lambda}
We assumed in this case that $1 <\gamma<\beta-1$, so that $\beta > 2$ and
\[
\dfrac{2(\beta-1)}{\beta} > 1.
\]
The bulk viscosity coefficient $\big(\lambda_{\ep,\delta}(\rho_{\ep,\delta})\big)_\delta$ is then bounded in $L^1((0,T) \times \TT)$.
\end{remark}

\noindent
\textit{Limit $\delta \rightarrow 0$.}
We can pass to the limit in the weak formulation of the equations exactly in the same way as before. We show first the strong convergence of the density and use then the control of $\sqrt{\lambda_{\ep,\delta}(\rho_{\ep,\delta})}$ in $L^2\big((0,T)\times \TT\big)$ to identify the limit of the bulk viscosity term:
\[
\lambda_{\ep,\delta}(\rho_{\ep,\delta})\Div u_{\ep,\delta}
\to
\lambda_\ep(\rho_\ep)\Div u_\ep\quad \text{weakly in}~ L^1((0,T)\times \TT).
\]

\subsubsection{Limit $\ep \rightarrow 0$}
Thanks to the uniform controls of the singular quantities, we ensure the following convergences
\begin{gather*}
\Lambda_\ep(\rho_\ep) \to \Lambda \quad \text{weakly-* in} ~ L^\infty(0,T;L^2(\TT)),\\
\lambda_\ep(\rho_\ep) \Div u_\ep \to \Pi \quad \text{weakly in} ~
\mathcal{M}((0,T)\times \TT).
\end{gather*}
Invoking similar arguments as those used to obtain \eqref{eq:spt-p}, we can get
\[
\spt \Lambda \subset \{ \rho = 1\},
\]
or, written differently (the product $\rho \Lambda$ making sense almost everywhere):
\[
(1-\rho) \Lambda = 0.
\]
On the other hand, we have $\gamma < \beta-1$, which yields
\begin{equation*}
p_\ep(\rho_\ep) \to 0 \quad\text{strongly in} ~ L^1((0,T)\times \TT).
\end{equation*}
Compared to the previous case, the limit $\Lambda$ is not $0$ and we close the system by passing to the limit in
\[
\partial_t \Lambda_\ep(\rho_\ep) + \Div(\Lambda_\ep(\rho_\ep) u_\ep) = - \lambda_\ep(\rho_\ep)\Div u_\ep.
\]
Here, $(\Lambda_\ep)_\ep$ is controlled in $L^\infty(0,T; L^2(\TT))$ and $\|\partial_t \Lambda_\ep\|_{L^1(W^{-1,1})}\leq C$, while $(u_\ep)_\ep$ is bounded in $L^2(0,T;(H^1(\TT))^3)$.
We can then pass to the limit in the product $\Lambda_\ep(\rho_\ep) u_\ep$ thanks to Lemma~\ref{lem:compensated_compactness} and get the limit equation
\[
\partial_t \Lambda + \Div(\Lambda u) = -\Pi \quad \text{in} ~ ~\mathcal{D}'.
\]
We have thus justified the activation of memory effects in the congestion limit.
The tuple $(\rho, u, \pi,\Lambda)$ is finally a global weak solution of
\begin{subnumcasess}{\labeleq{xxx}}
\partial_t \rho + \Div (\rho u) = 0, \\
-\nabla \Pi - 2 \Div(\mu \D(u)) + r u = f, \\
\partial_t \Lambda + \Div(\Lambda u) = -\Pi,\\
0 \leq \rho \leq 1, ~(1-\rho)\Lambda = 0,~ \Lambda \geq 0.
\end{subnumcasess}%

\subsection{Case $1<\gamma = \beta-1$}
The estimates remain mostly unchanged compared to the previous case $\gamma < \beta-1$.
Nevertheless, in the limit $\ep \rightarrow 0$, we do not have the convergence of $p_\ep(\rho_\ep)$ to $0$ anymore and $p_\ep(\rho_\ep)$ converges now to $p = (\beta-1) \Lambda$.
The~limit system then writes
\begin{subnumcasess}{\labeleq{xxx2}}
\partial_t \rho + \Div (\rho u) = 0, \\
(\beta-1)\nabla \Lambda - \nabla \Pi - 2 \Div(\mu \D(u)) + ru = f, \\
\partial_t \Lambda + \Div(\Lambda u) = -\Pi,\\
0 \leq \rho \leq 1, ~(1-\rho)\Lambda = 0,~ \Lambda \geq 0.
\end{subnumcasess}%

\backmatter
\bibliographystyle{jepplain+eid}
\bibliography{bresch-et-al}
\end{document}