\section{Introduction} \label{sec.introl}
The modeling of multiphase flow through porous media is an important task in many engineering applications like for example geothermal systems, oil reservoir engineering, ground-water hydrology, and thermal energy storage.  In unsaturated soils the flow can be essentially described by a two-phase flow of two immiscible fluids: air and water. A prominent example for the physical processes involved in this kind of flow is the infiltration of contaminants into upper layers of unsaturated soil 
where the air present in the soil is taken to be at a constant pressure. 
This situation can be modeled by the nonisothermal Richards equation  where the water-phase  consists of a mixture of $N$ chemical components undergoing cross-diffusion effects. The model is based on equations representing mass, momentum and energy balances, and are coupled to constitutive relations derived in a thermodynamically consistent way starting from the Helmholtz free energy. 
Compared to the main difficulties related to the mathematical analysis of the similar model considered by Amaziane, Jurak, Pankratov and Piatnitski in \cite{AJPP17} (coupling, strong nonlinearity, degeneracy of the diffusion term in the saturation equation), in this article we additionally deal with cross-diffusion terms which appear in the diffusion and in the thermal flux due to the presence of chemical components in the water-phase. 
 Note also that in \cite{AJPP17} the authors first regularized the phase pressures before they were using the time discretization to get a sequence of elliptic problems on which they could apply the fixed point argument,
while in our approach, we apply the techniques of Buli\v cek, J\"ungel, Pokorn\'y and Zamponi \cite{BJPZ20} taking into account the relation between the capillary pressure function and the interfacial boundary energy \eqref{def.Eint} introduced by Jurak, Koldoba, Konyukhov and Pankratov in \cite{JKKP19}.  


 The mathematical analysis of multiphase flows in porous media started due to their importance in applications in geological sciences and reservoir simulations. The first results on the existence of weak solutions of a simplified system describing a nonisothermal two-phase flow in porous media was obtained by Bocharov and Monakhov in \cite{BoMo88.1, BoMo88.2}. In \cite{LiSun10} the authors proved the global existence of weak solutions for the nonisothermal, one-dimensional multicomponent heat-air-vapor transport model in porous textile materials. Existence of global weak solutions for Richards' model arising from the heat and moisture flow through a partially saturated porous medium was obtained by Bene\v s and Pa\v zanin in \cite{BePaz17}. The existence result for a nonisothermal immiscible incompressible two-phase flow in heterogeneus porous media was given by Amaziane, Jurak, Pankratov and Piatnitski in \cite{AJPP17} where the authors used the concept of the nonisothermal global pressure introduced by Bocharov and Monakhov in \cite{BoMo88.1} as the crucial mathematical  tool to obtain a priori estimates and  compactness results. Note that all these results were obtained for the single-species case only.
 
The mathematical existence analyis of cross-diffusion systems has been investigated extensively in the last two decades using entropy methods (see \textit{e.g.} \cite{BSW12, Jue16, Jue15}), but the focus on how to model cross-diffusion effects (for instance in fluid mixtures) in a thermodynamically consistent way attracted the mathematical community's attention only recently, see \textit{e.g.} \cite{BH15, BJPZ20, BPZ17, DMZ20, Drey20, FHKM20, HJ21, JMZ18}. For seminal works on the thermodynamic theory of fluid mixtures we want to mention \textit{e.g.} \cite{BD15, Gio99}. 
In general, cross diffusion means that the concentration gradient of one species induces a flux of another species.
 Mathematically, cross diffusion is usually described by strongly coupled parabolic systems of reaction-diffusion type, where the diffusion matrix is non-diagonal and in general neither symmetric nor positive definite, but the system possesses a (formal) gradient-flow structure which helps to perform the analysis without the use of a maximum principle, which is in general not available for such strongly coupled systems. Instead, the key tool in these entropy methods \cite{BSW12, Jue16, Jue15} is to perform a change of variables into so called \textit{entropy variables}, which turn out to have also a physical meaning in terms of the chemical potentials of the system. 
 
 We want to give a brief overview of some of the works performed so far in the direction of thermodynamically consistent modelling of physical systems with cross-diffusion terms. 
 In \cite{BH15}, Buli\v{c}ek and Havrda studied an incompressible ﬂuid mixture, where the diffusion of all species and the heat
ﬂux were driven by the gradient of the chemical potentials and of the temperature, allowing for cross effects in such a way that the first and second laws of thermodynamics were still satisfied.  
 In \cite{Drey20}, Dreyer, Druet, Gajewski and Guhlke performed the existence analysis for an improved Nernst-Planck-Poisson model of compressible isothermal electrolytes with cross-diffusion phenomena due to mass conservation. In \cite{BPZ17} an incompressible, multicomponent, heat conducting, electrically charged fluid model with cross-diffusion effects was studied. In \cite{JMZ18} a single-phase fluid mixture in porous media was described and analyzed by using cross-diffusion equations. In \cite{DMZ20} the authors studied an isothermal, unsaturated two-phase flow mixture with dynamic capillary pressure in porous media with cross-diffusion effects, where the assumption of constant temperature simplified the analysis. A recent breakthrough was made by Buli\v{c}ek, J\"{u}ngel, Pokorn\'y and Zamponi in \cite{BJPZ20} by studying a steady compressible thermodynamically consistent Navier-Stokes-Fourier system for chemically reacting fluids with cross-diffusion effects. Very recent results were performed in \cite{HJ21}, where the Maxwell-Stefan-Fourier equations in Fick-Onsager form were investigated and in \cite{FHKM20}, where a thermodynamically consistent reaction-cross-diffusion system coupled to an equation for the heat transfer was studied. The aim of this paper, namely to study in a thermodynamically consistent way the nonisothermal immiscible compressible two-phase unsaturated flow in porous media with cross diffusion, is completely new, to our best knowledge. 

 
 This paper is organized as follows. 
 In the next section we start with an introduction of the main physical quantities, the presentation of the governing equations of the model of interest, followed by the constitutive relations and main mathematical hypotheses. In subsection \ref{SolConc} we present the concept of weak solutions called variational entropy solution. Section \ref{sec.apriori} is devoted to the derivation of the entropy balance equation which is the starting point for obtaining a priori estimates for smooth solutions. In Section \ref{sec.weakstab}, we deal with the weak sequential stability for smooth solutions. In Section \ref{sec.appr}, we introduce the approximate problem using different regularization parameters. The existence of solutions to the approximate system is shown using the Leray-Schauder fixed point argument. Finally, using the compactness results established in Section \ref{sec.weakstab}, we perform the limit as the time-discretization and the regularization parameters go to zero, obtaining finally a solution of the degenerate system. We conclude the article with an Appendix where we present some auxiliary technical results.


\subsection{The physical quantities and main assumptions}\label{sec.phys}
We consider a nonisothermal Richards flow process with cross-diffusion effects in a porous reservoir $\Omega \subset \R^3$, which is a bounded Lipshitz domain, with a nondeformable skeleton. The time interval of interest is $(0,\mathcal{T})$ and $\mathcal{Q} = \Omega \times (0,\mathcal{T})$.
The indexes $w,n,s$ correspond to the wetting (water) phase, the non-wetting (air) phase, and the skeleton. \\
{\sf The temperature:} $T = T(x,t)$ stands for the temperature in $\Omega$. We assume that the temperature is locally in a heat equilibrium state, i.e. $T_w = T_n = T_s = T$.\\
{\sf The phase pressures:} by $p_w = p(x,t)$ we denote the water pressure while $p_n = p_{at}$ represents the constant air-pressure.\\
{\sf The water mass density:} $\rho_w = \rho(T,p)$. We assume that the water phase is a mixture of $N$ different components. 
We denote the mass density of the $i-$th water component by $\rho_i$, for $i=1,\ldots,N$. It holds that $\rho = \sum_{i=1}^N \rho_i$. Moreover, let be $\vec\rho = (\rho_1,\ldots,\rho_N)$.\\
{\sf The saturation functions:} By $S_w = S(x,t)$ we denote the water saturation in $\Omega$.\\
 {\sf The energy densities of the water-phase  and skeleton energy:} By $E_w$ we denote the energy density of the water phase.  We assume that the skeleton energy density $E_s$ depends on the temperature only, i.e. $E_s = E_s(T)$. 
   Following \cite{JKKP19} the total energy density $E_f$ is given by
\begin{align}
\label{def.Ef}
E_f = E_w S  + E_{int}(S).
\end{align}
By denoting with $e$ the specific internal energy, we have $E_w = (\rho e)$.\\

{\sf The capillary pressure:}  The constitutive equation relating the capillary pressure $P_c$ to the water-phase saturation $S$ is classically given as an algebraic relationship between $P_c$ and $S$: $P_c(S) = p_n - p_w = p_{at} - p$. A detailed discussion of this relationship can be found e.g. in \cite{Bear68}. The relationship between $P_c$ and $S$ has been generalized on the basis of thermodynamical arguments by Gray and Hassanizadeh \cite{HassGrey93}. They derived an extended relationship, which, in the case of Richards flow and after scaling, can be written in the form    
 \begin{equation}
 p = -P_c(S) + \partial_t f(S),   
 \label{DCP}
 \end{equation}
 where  $f$ denotes the dynamic capillary pressure function.\\
{\sf The skeleton and the interfacial boundary energies:}
      Following \cite{JKKP19, Coussy04}, by $E_{int} = E_{int}(S)$ we denote the energy of the interphase density. The capillary pressure function is related to the interfacial boundary energy by
      \[  P_c(S) = - \frac{\partial E_{int}}{\partial S}.\]
  It follows that, up to a constant, the interface energy density (i.e.~the energy density of the boundary between the two phases) reads as
  \begin{align}
\label{def.Eint}
E_{int}(S) &= \int_{S}^1 P_c(\xi)d\xi.
\end{align}
{\sf The viscosity function:} By $\mu_w = \mu(T)$ we denote the viscosity of water in $\Omega$.\\
{\sf The porosity function:} By $\Phi$ we denote the porosity of the domain $\Omega$. We assume that the porosity $\Phi : \Omega\to [0,1]$ is a Lebesgue-measurable function such that
\begin{equation}
\label{Phi.pos}
\mbox{ess}\inf_{\Omega}\Phi > 0,\qquad 
\mbox{ess}\sup_{\Omega}\Phi < 1.
\end{equation}
{\sf The relative permeability function:} By $k_{r,w} = k_r(S)$ we denote the relative permeability of the water-phase in $\Omega$. We assume that $k_r$ is a continuous function that satisfies: $(i)$ $0 \leq k_r \leq 1$ on $\R$; $(ii)$ $k_r(S) = 0$ for $S \leq 0$ and $k_r(S) = 1$ for $S \geq 1$.\\
{\sf The mobility functions:}  The water mobility function $\lambda_w = \lambda(S,T)$ is defined by
        \begin{equation}
	\label{mobility}
\lambda(S,T) = \frac{k_{r}(S)}{\mu(T)},
\end{equation}
where $\mu(T)$ is the viscosity of the water phase.\\
{\sf The water-phase velocity:}   Following \cite{Bear68}, the generalized Darcy law gives
\begin{align}
	\label{vow}
\veb{v} = -K\lambda(S,T)\nabla p,
\end{align}
where $\lambda$ is the mobility function, and $0<K<+\infty$ is the constant absolute permeability of the domain $\Omega$.\\
{\sf The chemical potentials:} By $\mu_1,\ldots,\mu_N$ we denote the chemical potentials of the contaminants in the water-phase.\\
{\sf The diffusion and the thermal flux:}  
Following thermodynamics of irreversible processes (TIP) in linear approximations \cite{GrootMazur}, we have the following 
expressions for the diffusion and the thermal fluxes
(provided that the external forces are neglected):  \\        
  (i) the diffusion flux  of the $i$th species is denoted by  $\veb{J}_i$ and is given by
        \begin{align}
  \veb{J}_i = L_{i0} \nabla \Big(\frac{1}{T}\Big) - \sum_{j=1}^{N} L_{ij} \nabla \Big(\frac{\mu_j}{T}\Big), \qquad i=1,\ldots,N,  \label{flux.J}
\end{align}
where $L_{ij} = L_{ij}(\vec{\rho},T)$ are the diffusion coefficients (mobilities) which form the mobility matrix $(L_{ij})_{i,j=1,\ldots,N}$.\\
  (ii) the heat flux is denoted by $\veb{q}$ and consists of Fourier's law and the molecular diffusion term
         \begin{align}
  \veb{q}  = -\kappa(T) \nabla T + \sum_{j=1}^{N} L_{0{j}} \nabla \Big( \frac{\mu_j}{T} \Big),
  \label{flux.q}
\end{align}
where $\kappa(T) > 0$ is the heat conductivity given by
\begin{equation}
 \kappa(T) = L_{00} T^{-2}.  \label{heat_cond}
\end{equation}  

\subsection{Model equations, constitutive relations and hypotheses}\label{sec.model}

  Let $\veb{v}$ be the water velocity and $\veb{v}_i$ be the velocity of the $i-$th water component, for $i=1,\ldots,N$.
A water mixture (barycentric) velocity $\veb{v}$ is defined as 
\[ \rho \veb{v} = \sum_{i=1}^N \rho_i \veb{v}_i. \] 
Each component of the water mixture satisfies the mass conservation law, i.e. 
\begin{equation}\label{eq.comp}
\Phi\frac{\partial}{\partial t} (S \rho_i) + \dv(  \rho_i \veb{v}_i ) = {r_i},\quad i=1,\ldots,N,
\end{equation}
{where $r_1,\ldots,r_N$ are reaction terms modeling e.g.~chemical reactions between the components in the water mixture. We assume that $\sum_{i=1}^N r_i = 0$ (total mass conservation).}

By summing these equations we get the total mass conservation of the water-phase:
\begin{equation*}
\Phi\frac{\partial}{\partial t} (S \rho) + \dv(  \rho\veb{v} ) = 0.
  
\end{equation*}
The transport of each water-component is divided into convective part ($\rho_i \veb{v}$)  and diffusive part
($\veb{J}_i$) by dividing the component flux    
\begin{equation*}
  \rho_i \veb{v}_i = \rho_i \veb{v} + \veb{J}_i,
\end{equation*}
where the diffusive flux is defined as $\veb{J}_i=\rho_i (\veb{v}_i - \veb{v})$. Obviously, 
$  \sum_{i=1}^N \veb{J}_i = 0.  $

The model consists of the mass conservation equation \eqref{eq.comp} for each water-component phase $i$, the energy conservation equation taken from \cite{Kav95} and the equation for the water saturation \eqref{DCP}:
\begin{align}
\Phi  \frac{\partial}{\partial t} (S \rho_i) + \dv(  \rho_i\veb{v} + \veb{J}_i ) &= {r_i},\quad i=1,\ldots,N,  \label{cons_w} \\
  \frac{\partial}{\partial t}\big (
\Phi E_f + (1-\Phi) E_s \big) 
+ \dv \big( ((\rho e) + p) \veb{v}  + \veb{q} \big) & = 0,   \label{cons_e} \\
\partial_t f(S) + P_c(S) + p &= 0. \label{PC.Sw}
\end{align}
We assume the following initial conditions
\begin{align}
	\label{ic.1}
	&\rho_i(\cdot,0) = \rho_i^{in},\quad
	S(\cdot,0) = S^{in},\quad
	T(\cdot,0) = T^{in}\quad\mbox{in }\Omega,~~ i=1,\ldots,N ,
\end{align}
where the initial data are Lebesgue measurable functions such that 
\begin{align}
	\label{hp.data}
\rho_i^{in}>0,\quad 0\leq S^{in}\leq 1,
\quad T^{in}>0,\quad\mbox{a.e.~in }\Omega,\\
\label{hp.data.2}
S^{in}\rho^{in}, \quad 
\Phi E_f^{in} + (1-\Phi)E_s^{in},\quad
\Phi S^{in} (\rho\eta)^{in} + (1-\Phi)(\rho\eta)_s^{in}\in L^1(\Omega).
\end{align}
 Equations \eqref{cons_w}-\eqref{PC.Sw} are solved in a bounded domain $\Omega \subset \R^3$ and are supplemented with complete slip boundary conditions for the velocity, and Robin boundary conditions for the diffusion fluxes and the heat flux on $\partial \Omega \times (0,\infty)$, $i=1,\ldots,N$:
\begin{align}\label{bc}
&\veb v \cdot \nu = 0,\quad
 \veb J_i\cdot\nu = \sum_{k=1}^Nb_{ik} 
 \left( \frac{\mu_k}{T} - \frac{\mu_{0,k}}{T_0} \right),\quad 
\veb q \cdot \nu = \alpha (T-T_0),
\end{align}
where $\alpha>0$, $T_0>0$, $\mu_{0,1},\ldots,\mu_{0,N}\in\R$ 
are scalar constants, while
$(b_{ij})_{i,j=1,\ldots N}\in\R^{N\times N}$ is a constant, symmetric, positive semidefinite matrix such that
\begin{align}
\label{hp.b}
\sum_{i=1}^N b_{ij} = 0\qquad j=1,\ldots,N.
\end{align}
{\bf Constitutive relations.}
We assume that the thermodynamic quantities appearing in \eqref{cons_w} and \eqref{cons_e}
are induced by the Helmholtz free energy density $\rho\Psi$ in a thermodynamically consistent way. Precisely, we assume
that the free energy density $\rho\Psi$ decomposes into a sum of water- and skeleton-related contributions
\begin{align}
&\rho \Psi  = (\rho \Psi)_w + (\rho\Psi)_s , \label{Helm}\\ \nonumber
& (\rho \Psi)_w = (\rho \Psi)_w(\vec\rho, T),\quad
(\rho\Psi)_s = (\rho\Psi)_s(T).
\end{align}
Following \cite[Remark 1.2]{BJPZ20}, we take the following expression for the functions $(\rho \Psi)_w, (\rho\Psi)_s$, under the simplifying assumption that the molar masses of each water-component phase are all coincident:
\begin{align}\label{Helm.w}
(\rho\Psi)_w(\vec\rho,T) &= T\sum_{i=1}^N\rho_i\log\rho_i + \rho^\gamma - c_w\rho T\log T  + p_{at},\\
\label{Helm.s}
(\rho\Psi)_s(T) &= T -c_s T\log T .
\end{align}
Here $c_w$, $c_s>0$ are (scaled) heat capacities (of water and skeleton, respectively), $p_{at}>0$ is the atmospheric pressure, while 
\begin{equation}
\label{gamma}
\gamma > 2 .
\end{equation}
As a consequence,
the chemical potentials, pressure, internal energy density and entropy density for each phase are given by the following relations.\smallskip\\
{\sf Chemical potentials:}
 \begin{align}
 \mu_i  = \frac{\partial(\rho \Psi)_w}{\partial \rho_i}
 = T(\log\rho_i + 1) + \gamma\rho^{\gamma-1} - c_w T\log T,
 \quad i=1,\ldots,N. \label{chp.1}
\end{align}
{\sf Phase pressure:}
\begin{align}
  p  = -(\rho \Psi)_w  + \sum_{i=1}^N \rho_i \mu_i
  = T\rho + (\gamma-1)\rho^\gamma - p_{at}. \label{thp.1}
\end{align}
{\sf Phase energy density: } 
\begin{align}
 (\rho e)  = (\rho \Psi)_w - T \frac{\partial(\rho \Psi)_w}{\partial T} = \rho^\gamma + c_w\rho T + p_{at}.
 \label{eng.1}
\end{align}
{\sf Phase entropy densities:}
\begin{align}\label{entropy}
(\rho \eta) = -\frac{\partial(\rho \Psi)_w}{\partial T} = 
-\sum_{i=1}^N\rho_i\log\rho_i + c_w\rho(\log T + 1).
\end{align}
We point out that the Gibbs-Duhem relations hold
\begin{align}
	\label{GibbsDuhem}
	T d(\rho\eta) = d(\rho e) - \sum_{i=1}^N\mu_i d\rho_i ,
\end{align}
where the differential $d$ stands here for $\partial_t$ or $\partial_{x_i}$, $i=1,2,3$.\smallskip\\
{\sf Skeleton entropy and energy:}
\begin{align}
(\rho \eta)_s &= -\frac{\partial(\rho \Psi)_s}{\partial T} =
c_s - 1 + c_s\log T,
\label{Skel.Entr}\\
E_s &= (\rho \Psi)_s - T \frac{\partial(\rho \Psi)_s}{\partial T}
= c_s T. \label{Skel.Ener}
\end{align}
{\bf Hypotheses.} 
In this subsection we collect the main mathematical hypotheses which we imposed throughout this article.\\
{\sf Phase velocity.}
We also assume that $\mu : \R_+\to\R_+$ is a continuous, uniformly positive and bounded function, while $k_r\in C^0([0,1])$ is an increasing, nonnegative function such that
 
\begin{equation}
\label{hp.kr}
\exists \alpha_r \in\left(\frac{2}{\gamma},\,
\frac{14}{3} \right),\quad
\exists k_r^*>0:\quad
\lim_{s\to 0} s^{-\alpha_r}k_r(s) = k_r^* .
\end{equation}
{\sf Capillary pressure.}
We assume that $P_c , f : (0,1)\to\R$ are $C^1$ strictly decreasing functions such that
\begin{align}
\label{hp.Pcf}
\exists c_{f}>0,~~ \frac{\gamma}{\gamma-1}\leq q<2 :\quad 
\inf_{0<s<1/2}\frac{|P_c'(s)|k_r(s)^{\frac{q}{2(q-1)}}}{
	|f'(s)|} \geq c_{f},\\
\label{hp.krf}
\lim_{s\to 1}f(s) = -\infty,\quad
\exists c_f'>0:\quad 
\left|\frac{d}{ds}\sqrt{k_r(s)}\right|\leq c_f'
\left|f'(s)\right|,\quad 0<s<1,\\
\label{hp.Pcf2}
f(0)>0,\quad P_c(S)>0,\quad
\exists s_0>0:\quad
\inf_{s\in (0,s_0)}\frac{P_c(s)}{f(s)} > \frac{p_{at}}{f(0)},\\
\label{hp.Pc.bound}
\exists c_{p}>0,~~ k_p\geq 0:\quad
\lim_{s\to 0^+} s^{k_p} P_c(s) = c_p ,\\
\label{lb.dfeps}
\inf_{(0,1)}|f'| > 0,\qquad 
\int_0^{1/2}|f'(u)| |\log u| du < \infty .
\end{align}
 

{\sf Heat conductivity.}
We assume that the heat conductivity $\kappa(T)$ satisfies the following assumption
\begin{align}
	\label{kappa}
\exists \kappa_1, \kappa_2>0:\quad
\kappa_1(1+T^\beta)\leq\kappa(T)\leq \kappa_2(1+T^\beta),
\end{align}
where the exponent $\beta>0$ satisfies
\begin{equation}
\label{hp.beta}
\beta\geq \frac{q}{2-q},\quad \beta > \frac{\gamma}{\gamma-2},\quad \beta\geq\frac{4}{3}.
\end{equation}
We point out that assumptions \eqref{lb.dfeps} and $\beta\geq 4/3$ are only needed in the construction of the approximate scheme.
A possible way to consider the case $\beta<4/3$ might be to adapt the (more involved) stationary scheme presented in \cite{BJPZ20} to the time-dependent system \eqref{cons_w}--\eqref{PC.Sw}. However, we prefer here to restrict ourself to the case $\beta\geq 4/3$ for the sake of simplicity.
{\sf Mobilities.} The Onsager-Casimir 
reciprocity relations imply the symmetry of the mobilities $L_{ij}=L_{ji}$, $L_{i0}=L_{0i}$ for $i,j=1,\ldots,N$,
and the second law of thermodynamics requires 
the positive semidefiniteness of the matrix 
$(L_{ij}) \in \mathbb{R}^{N\times N}$ as well as 
$L_{00}\geq 0$. 

Furthermore, we assume that $L_{ij} = L_{ji} = \tilde L_{ij}(\vec\rho,T)$,
$L_{i0} = L_{0i} = \tilde L_{i0}(\vec\rho,T)$ for $i,j=1,\ldots,N$, where
$\tilde L_{ij}, \tilde L_{i0}\in C^0(\R^N_+\times\R_+)$ such that
\begin{equation}
\label{hp.Ltilde}
|\tilde L_{ij}(\vec\rho,T)| + \frac{1}{T}|\tilde L_{i0}(\vec\rho,T)|
\leq C,\quad i,j=1,\ldots,N.
\end{equation}
Moreover, we assume that there exists $C, C'>0$ such that
\begin{align}
		\label{Ass.L}
		C |\Pi^N \vec{u} |^2\leq 
	\sum_{i,j=1}^N \tilde L_{ij}(\vec\rho,T) u_i u_j \leq 
	C' |\Pi^N \vec{u} |^2 , \quad 
	\forall \vec{u}\in\R^N,~~ \vec\rho\in\R_+^N,~~ T>0,
\end{align}  
where $\Pi^N = \Ib - \veb{1} \otimes \veb{1}/N$, is the orthogonal projector on $\sspan \{\veb{1}\}^\top$.
Thus,
\begin{align}
	\label{Pi}
(\Pi^N \vec{u})_i = u_i - \frac{1}{N} \sum_{j=1}^N u_j,\quad i=1,\ldots,N.
\end{align}
We also point out that, being $(L_{ij})_{i,j=1,\ldots,N}$ symmetric and positive semidefinite, \eqref{Ass.L} implies that 
$\sum_{i=1}^N L_{ij}=0$ for $j=1,\ldots,N$.\smallskip\\
{\sf Reaction terms.} For $i=1,\ldots,N$ we assume $r_i = \tilde r_i(\rho,T,\Pi^N\vec\mu/T)$ for some continuous function $\tilde r_i : \R_+\times (0,\infty)\times\R^N\to\R$ such that, for every $\rho\geq 0$, $T>0$, $\vec\zeta\in\R^N$,
\begin{align}
	\label{hp.r}
&	\sum_{j=1}^N \tilde r_j(\rho,T,\vec\zeta) = 0,\\
\label{hp.r.2}
&\exists C_0, C_1, C_2>0,~~
 a>2:\quad
\begin{cases} \displaystyle
	\sum_{j=1}^N \tilde r_j(\rho,T,\vec\zeta)\zeta_j 
	\leq C_0 - C_1\left|\Pi^N\vec\zeta\right|^{a},\smallskip\\
	\sum_{j=1}^N |\tilde r_j(\rho,T,\vec\zeta)|\leq 
	C_2\left(1 + \left|\Pi^N\vec\zeta\right|^{a-1}\right),
\end{cases}\\
\label{hp.r.3}
&\forall \rho > 0, T>0,
\mbox{ the mapping }
 \vec\zeta\in\R^N\mapsto -\sum_{j=1}^N \tilde r_j(\rho,T,\vec\zeta)\zeta_j \in\R
\mbox{ is convex,}\\
\label{hp.r.4}
&\exists\omega\in C^0(\R_+),~~\omega(0)=0:\quad\
\left|\sum_{j=1}^N (\tilde r_j(\rho,T,\vec\zeta)
-\tilde r_j(\rho',T',\vec\zeta))
\zeta_j\right|\leq\omega(|\rho-\rho'|+|T-T'|)|\vec\zeta|^a\\
\nonumber
&\hspace{11.3cm}\forall \rho,\rho',T,T'>0,~~ \vec\zeta\in\R^N.
\end{align}
\subsection{Solution concept, main result and the methodology of the proof}\label{SolConc}

Here we define the concept of a weak solution to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}, \eqref{bc} and state the main result of the paper. Due to the fact that some terms from the total energy balance are non-integrable in the classical weak sense, we consider so-called variational entropy solutions (following the terminology in \cite{NovPok11}), i.e. solutions which fulfill \eqref{cons_w} weakly, the integrated form of the total energy balance \eqref{cons_e}, and the weak formulation of the entropy balance equation \eqref{EBE}. 

\begin{definition} (Variational entropy solution.) \label{def.weaksol}
	A variational entropy solution to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}, \eqref{bc} is a Lebesgue-measurable function 
	$(\vec\rho,T,S): \Omega\times\R_+
	\to \R^N_+\times\R_+\times [0,1]$ such that, for every $\mathcal{T}>0$ and $i=1,\ldots,N$
\begin{align*}
P_c(S), \quad p,\quad \rho^\gamma \in L^1(Q_\mathcal{T}),\\
S\rho^\gamma,\quad  
S(\Phi E_f + (1-\Phi) E_s) \in L^\infty(0,\mathcal{T}; L^1(\Omega)),\\
\log T,\quad T^{\beta/2} \in L^2(0,\mathcal{T}; H^1(\Omega)), \quad
T \in L^\infty(0,\mathcal{T}; L^1(\Omega))\cap 
L^{ \beta + \frac{2}{3} }(Q_\mathcal{T}),\\
\Pi^N \big( \vec\mu/T \big) \in L^2(0,\mathcal{T}; H^1(\Omega))\cap L^a(Q_\mathcal{T}),\quad
\sqrt{\frac{\lambda(S,T)}{T}}\nabla p \in L^2(Q_\mathcal{T}),\\
\frac{1}{T},\quad
\sum_{i,j=1}^N b_{ij}
\left(\Pi^N \frac{\vec\mu}{T} \right)_i
\left(\Pi^N \frac{\vec\mu}{T} \right)_j
\in L^1(\pa\Omega\times (0,\mathcal{T})),\\
T^{-1} f'(S) \big( \partial_t S \big)^2\in L^1(Q_\mathcal{T}),\quad
f(S)\in L^\infty(0,\mathcal{T}; W^{1,q}(\Omega)) ,
\end{align*}
for some suitable exponent
$a > 2$ and $\beta$, $q$ given by \eqref{hp.beta}, \eqref{hp.Pcf} (respectively),
the weak formulation of the partial mass balances \eqref{cons_w} holds
	\begin{align}
		\label{weak.1}
		&\int_0^\mathcal{T} \langle \Phi\pa_t (S\rho_i) , \varphi\rangle dt -
		\int_0^\mathcal{T}\int_\Omega(\rho_i \veb{v} + \veb{J}_i)\cdot\nabla\varphi dx dt\\
		\nonumber
		&\qquad +\int_0^\mathcal{T}\int_{\pa\Omega}\sum_{k=1}^Nb_{ik}\left(\frac{\mu_k}{T} - \frac{\mu_{0,k}}{T_0}\right)\varphi d\sigma dt = 
		\int_0^\mathcal{T}\int_\Omega r_i\varphi dx dt,\\
		\nonumber
		& i=1,\ldots,N,\qquad
		\forall \varphi\in C^1(\overline{Q_\mathcal{T}}),
	\end{align}
as well as the integrated energy balance
\begin{align}\label{weak.2}
	&\int_\Omega\big (\Phi E_f(t) + (1-\Phi) E_s(t) \big)dx
	+ \alpha\int_0^{t}\int_{\pa\Omega} (T-T_0)d\sigma dt \\
	&\qquad = 
	\int_\Omega\big (\Phi E_f^{in} + (1-\Phi) E_s^{in} \big)dx,
	\quad\forall t \in [0,\mathcal{T}], \nonumber
\end{align}
	eq.~\eqref{PC.Sw} holds a.e.~$x\in\Omega$, $t>0$, 
	the constitutive relations listed in subsection \ref{sec.model} 
	are satisfied a.e.~$x\in\Omega$, $t>0$, 
	the entropy balance equation holds 
	\begin{align}\label{weak.EBE}
		&\int_0^\mathcal{T}\langle\partial_t\left[\Phi S(\rho\eta) + (1-\Phi)(\rho \eta)_s \right],\varphi\rangle dt \\ \nonumber
		&\qquad -\int_0^\mathcal{T}\int_\Omega\Big( (\rho \eta) \veb{v} - \sum_{i=1}^N \frac{\mu_i}{T} \veb{J}_i + \frac{\veb{q}}{T} \Big)\cdot\nabla\varphi dx dt  \\
		& =\int_0^\mathcal{T}\int_\Omega\left( 
		\sum_{i,j=1}^N L_{ij} \nabla \big( \frac{\mu_i}{T} \big) \cdot \nabla \big( \frac{\mu_j}{T} \big) + L_{00} \big| \nabla \big( \frac{1}{T} \big) \big|^2\right)\varphi dx dt \nonumber\\ 
		&\qquad + \int_0^\mathcal{T}\left( 
		\frac{K}{T}\lambda(S,T)|\nabla p|^2
		-\Phi \frac{1}{T} f'(S) \big( \partial_t S \big)^2
		-\sum_{i=1}^N r_i\frac{\mu_i}{T}
		\right)  \varphi dx dt\nonumber	\\
		&\qquad +\int_0^{\mathcal{T}}\int_{\pa\Omega}\left(\alpha\frac{T_0 - T}{T}
		+\sum_{i,j=1}^N b_{ij}\frac{\mu_i}{T}\left(\frac{\mu_j}{T} - \frac{\mu_{0,j}}{T_0}\right)
		\right)\varphi d\sigma dt  + \langle\xi , \varphi\rangle,
		\nonumber
	\end{align}
	for every $\varphi\in C^1(\overline{Q_\mathcal{T}})$, $\varphi\geq 0$ a.e.~in $Q_\mathcal{T}$, for a suitable choice of $\xi\in \mathcal M(\overline{Q_\mathcal{T}})$, $\xi\geq 0$ (nonnegative Radon measure), 
	the saturation balance \eqref{PC.Sw} holds a.e.~in $Q_\mathcal{T}$, and	
	the initial conditions \eqref{ic.1} hold in the sense
	\begin{align}
		\label{weak.ic}
		&\lim_{t\to 0}\int_\Omega \Phi S(t)\rho_i(t)\varphi dx =
		\int_\Omega \Phi S^{in}\rho^{in}_i\varphi dx,\qquad
		i=1, \ldots , N,\\
		\label{weak.ic.1}
		&\lim_{t\to 0}\int_\Omega 
		(\Phi E_f(t) + (1-\Phi)E_s(t))\varphi dx =
		\int_\Omega (\Phi E_f^{in} + (1-\Phi)E_s^{in})\varphi dx,\\
		\label{weak.ic.2}
		&\lim_{t\to 0}S(t) = S^{in}\quad\mbox{strongly in $L^1(\Omega)$,}
	\end{align}
	for all $\varphi\in C^{1}(\overline\Omega)$.	
\end{definition}

The main result of this paper is the following theorem.

\begin{theorem} (Large-data existence of solutions.)\label{Thm1}\\
Let hypotheses \eqref{hp.kr}--\eqref{hp.r.2} hold. Then there exists a variational entropy solution triplet $(\vec\rho,T,S)$ to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}, \eqref{bc}.
\end{theorem}

\subsection{The methodology of the proof of Theorem \ref{Thm1}} 
In the proof of Theorem \ref{Thm1} we mostly follow the approach given in \cite{BJPZ20}. 
The main difficulties in the analysis of the system are the combined presence of cross-diffusion and degenerate terms, as well as the fact that the system is time dependent. The cross diffusion means a lack of parabolicity (in the standard sense) for the system, while the presence of degeneracies in some terms prevent us from obtaining useful gradient estimates. In fact,
because of the degeneracy of the relative permeability, we do not have uniform estimates for the gradient of the pressure or the gradient of the mass densities, which causes difficulties in the proof of the compactness result. Furthermore,
the strong nonlinearities in the equations imply that a clever use of the Div-Curl Lemma was needed in order to show compactness of  all the involved terms.
Moreover, in comparison to \cite{BJPZ20}, we do not have equations for the linear momentum, which is instead explicitly defined via Darcy's law. As a consequence, higher integrability estimates {\em \`{a}-la-}Bogovsky for the mass density \cite[Lemma 3.3]{BJPZ20} are unavailable. On the other hand, a higher integrability bound for the pressure on domains separated from the initial time (i.e.~$\Omega\times [\epsilon,\mathcal{T}]$, $\epsilon>0$ arbitrary) is derived from the equation for the capillary pressure, and is crucially employed in the compactness argument. Furthermore, while in \cite{BJPZ20} the stationary case is considered, the system studied here is time dependent, which clearly makes the analysis quite more complicated. Indeed, a consequence of dealing with the nonstationary case is the impossibility of obtaining a weak formulation for the energy balance (we only get an equation for the {\em integrated} energy balance) due to the lack of control of terms in the energy flux. 
To overcome the lack of parabolicity in the equations due to cross diffusion,
we start by proving the existence of solutions to an approximate system, which is uniformly elliptic with respect to a special set of variables (``entropy variables''), i.e.~chemical potentials and the logarithm of the temperature. Then we derive suitable a priori estimates from the entropy balance equation, and finally pass to the limit of vanishing approximation parameters. For the construction of a solution to the approximate problem, 
we employ a time semi-discretization (with timestep $\tau$) as well as lower-order and higher-order regularizations with parameters $\varepsilon > 0$, $\delta>0$ (respectively) of the mass conservation equation for each water-component phase and the conservation of the energy equation.
These levels of approximation allow us to
obtain $H^2(\Omega)$ regularity for the approximated, time-discretized solutions. In particular,
regularizations of the cross-diffusion operator and the source terms yield $H^2(\Omega)$ bounds for the chemical potentials, which imply that the (approximated) partial mass densities are bounded and uniformly positive. 
A further regularization of the free energy and the capillary pressure are required in order to retain higher integrability bounds for the mass densities and temperature, as well as uniform positivity of the saturation, after taking the time-continuous limit $\tau\to 0$ and the higher-order regularization limit $\delta\to 0$.
We point out that, while making the analysis more involved for the aforementioned reasons, the cross-diffusion terms yield in the end an $L^2$ bound for the relative chemical potentials, whose consequence is the following property: at any point of the considered space-time domain, either the total density vanishes (i.e.~there is a vacuum) or each partial mass density is positive. 
In order to prove the strong convergence of the sequence of approximated total density, partial densities and the temperature, we cannot employ the Aubin-Lions Lemma because of the high nonlinearities present in the equations; instead, 
a smart application of the Div-Curl lemma (a standard tool in fluid dynamics) is required in order to complete the compactness argument. 
Finally, we point out that the addition of a nonstationary term in the expression for the capillary pressure \eqref{PC.Sw} is mainly motivated by the need of a bound for the time derivative of (a function of) the saturation, which is crucial in the compactness argument.
 

\section{Entropy balance equation and a priori estimates for smooth solutions}\label{sec.apriori}
 In this section we will derive suitable a priori estimates for smooth solutions. The starting point is the following result, which shows that the entropy balance equation can be derived from the mass and energy balance equations \eqref{cons_w} and \eqref{cons_e}.
\begin{lemma}\label{lem.EBE} (Entropy balance equation.)\\
Let $(\vec\rho,T,S): \Omega\times\R_+
	\to \R^N_+\times\R_+\times [0,1]$ be a smooth solution to  \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}, \eqref{bc}. Then the entropy balance equation holds 
 \begin{align}\label{EBE}
&\frac{\partial}{\partial t}\left[\Phi S(\rho\eta) + (1-\Phi)(\rho \eta)_s \right] + \dv\Big( (\rho \eta) \veb{v} - \sum_{i=1}^N \frac{\mu_i}{T} \veb{J}_i + \frac{\veb{q}}{T} \Big)   \\
& = \frac{K}{T}\lambda(S,T)|\nabla p|^2 +
  \sum_{i,j=1}^N L_{ij} \nabla \big( \frac{\mu_i}{T} \big) \cdot \nabla \big( \frac{\mu_j}{T} \big) + L_{00} \big| \nabla \big( \frac{1}{T} \big) \big|^
- \Phi \frac{1}{T} f'(S) \big( \partial_t S \big)^2
-\sum_{i=1}^N r_i\frac{\mu_i}{T}. \nonumber
\end{align}
\end{lemma} 
 

\begin{proof}
The Gibbs-Duhem relations \eqref{GibbsDuhem} allow us to write
\begin{align}
 \partial_t (\rho \eta) &= - \sum_{i=1}^N \frac{\mu_i}{T}  \partial_t \rho_i +  \frac{1}{T}  \partial_t (\rho e) . \label{Eq.1}
\end{align} 
Next, from \eqref{Skel.Entr}, \eqref{Skel.Ener} 
it follows
\begin{align}\label{Eq.s}
 \partial_t (\rho \eta)_s =   \frac{1}{T} \partial_t E_s. 
\end{align} 
We also point out that, summing \eqref{thp.1} and \eqref{eng.1} and employing \eqref{chp.1} and \eqref{entropy} leads to
\begin{align}
	\label{rhoeta}
p + (\rho e) - \sum_{i=1}^N\rho_i\mu_i = T(\rho\eta) .
\end{align}

Let us multiply \eqref{cons_w} by $-\mu_i/T$ and sum from $i=1,\ldots,N$, multiply \eqref{cons_e} by $1/T$ and then sum the obtained equations. 
 In this way we get
 \begin{align}
&\Phi\left( 
-\sum_{i=1}^N\frac{\mu_i}{T}\pa_t(\rho_i S) +
\frac{1}{T}\pa_t E_f \right)
+(1-\Phi)\frac{1}{T}\pa_t E_s
\nonumber\\
&\qquad - \sum_{i=1}^N \frac{\mu_i}{T}\dv\big( \rho_i \veb{v} + \veb{J}_i  \big) + \frac{1}{T} \dv \big( (\rho e + p) \veb{v} + \veb{q} \big) = 
-\sum_{i=1}^N r_i\frac{\mu_i}{T}.
 \label{Eq.3}
 \end{align}
Let us now consider the first term on the left-hand side of \eqref{Eq.3}. From \eqref{def.Ef}, \eqref{def.Eint} we deduce
\begin{align*}
& -\sum_{i=1}^N\frac{\mu_i}{T}\pa_t(\rho_i S) +
\frac{1}{T}\pa_t E_f\\ 
&
= S\left( -\sum_{i=1}^N\frac{\mu_i}{T}\pa_t\rho_i
+\frac{1}{T}\pa_t (\rho e) \right)
+\left( -\sum_{i=1}^N\frac{\mu_i}{T}\rho_i
+\frac{1}{T}(\rho e) \right)\pa_t S 
- \frac{1}{T}P_c(S)\pa_t S .
\end{align*}
Using \eqref{Eq.1}, \eqref{rhoeta} and \eqref{PC.Sw} 
from the equation above we get
\begin{align*}
&-\sum_{i=1}^N\frac{\mu_i}{T}\pa_t(\rho_i S) +
\frac{1}{T}\pa_t E_f = \partial_t \big( S(\rho \eta) \big) + \frac{1}{T}f'(S)(\partial_t S)^2.
\end{align*}
Thanks to the above-written identity and \eqref{Eq.s},
equation \eqref{Eq.3} becomes
\begin{align*}
	&\Phi\pa_t \big( S(\rho\eta) \big)
	+(1-\Phi)\pa_t(\rho\eta)_s + \frac{1}{T} \Phi
	f'(S) (\partial_t S)^2
	\nonumber\\
	& - \sum_{i=1}^N \frac{\mu_i}{T}\dv\big( \rho_i \veb{v} + \veb{J}_i  \big) 
	 + \frac{1}{T} \dv \big( (\rho e + p) \veb{v} + \veb{q} \big) = -\sum_{i=1}^N r_i\frac{\mu_i}{T},
\end{align*}
which can be equivalently rewritten as
\begin{align}
	\nonumber
\Phi \partial_t & \big( S(\rho\eta) \big) + (1-\Phi) \partial_t (\rho \eta)_s 
+ \frac{1}{T}\Phi f'(S) (\partial_t S)^2\\
&+ \veb{v} \cdot \Big( -\sum_{i=1}^N \frac{\mu_i}{T} \nabla \rho_i + \frac{1}{T} \nabla (\rho e) 
+\frac{1}{T}\nabla p \Big) \nonumber \\
& +( \dv \veb{v}) \Big( -\sum_{i=1}^N \frac{\mu_i}{T} \rho_i + \frac{(\rho e)}{T} + \frac{p}{T} \Big) 
+ \dv \Big( -\sum_{i=1}^N \frac{\mu_i}{T} \veb{J}_i + \frac{\veb{q}}{T} \Big) \nonumber \\
 & = -\sum_{i=1}^N  \nabla \Big(  \frac{\mu_i}{T}\Big) \cdot \veb{J}_i + \nabla \Big( \frac{1}{T} \Big) \cdot \veb{q} -\sum_{i=1}^N r_i\frac{\mu_i}{T}.
\label{Eq.4}
\end{align}
Due to \eqref{GibbsDuhem} one has
\begin{align}
 \nabla (\rho \eta) &= -\sum_{i=1}^N \frac{\mu_i}{T} \nabla \rho_i + \frac{1}{T} \nabla(\rho e) ,
\label{Eq.5}
\end{align}
while from \eqref{flux.J} and \eqref{flux.q} it follows
\begin{align}
-\sum_{i=1}^N \nabla \Big( \frac{\mu_i}{T} \Big) \cdot \veb{J}_i & + \nabla \Big( \frac{1}{T} \Big) \cdot \veb{q}
  \nonumber \\ & = \sum_{i,j=1}^N L_{ij} \nabla \big( \frac{\mu_i}{T} \big) \cdot \nabla \big( \frac{\mu_j}{T} \big) + L_{00} \big| \nabla \big( \frac{1}{T} \big) \big|^2 \geq 0.
\label{Eq.9}
\end{align}
By inserting \eqref{Eq.5}-\eqref{Eq.9} into \eqref{Eq.4} and making use of \eqref{rhoeta} we conclude 
 \begin{align*}
\Phi \partial_t \big( & S(\rho\eta) \big) 
+ (1-\Phi) \partial_t (\rho \eta)_s + \veb{v} \cdot \nabla(\rho \eta)
 + ( \dv \veb{v}) (\rho \eta)\\
& + \frac{1}{T}\veb{v}\cdot\nabla p
+ \dv \Big( -\sum_{i=1}^N \frac{\mu_i}{T} \veb{J}_i + \frac{\veb{q}}{T} \Big)  \nonumber \\
 &
  = \sum_{i,j=1}^N L_{ij} \nabla \big( \frac{\mu_i}{T} \big) \cdot \nabla \big( \frac{\mu_j}{T} \big) + L_{00} \big| \nabla \big( \frac{1}{T} \big) \big|^2 
   - \frac{1}{T}\Phi f'(S) \big( \partial_t S \big)^2
   -\sum_{i=1}^N r_i\frac{\mu_i}{T},
\end{align*}
which, thanks to \eqref{vow}, yields the entropy balance equation \eqref{EBE}.
\end{proof}
\begin{remark}
The right-hand side of \eqref{EBE} is nonnegative
thanks to the assumptions that  $(L_{ij}) \in \R^{N \times N}$ is a positive semidefinite matrix, $L_{00}\geq 0$, the mobility function $\lambda$ in \eqref{mobility} is nonnegative, $f'(S)\leq 0$ and $-\sum_{i=1}^N r_i\frac{\mu_i}{T}\geq 0$. 
Therefore, the second law of thermodynamics holds true
 \begin{align}
\partial_t\left[\Phi S(\rho\eta) + (1-\Phi)(\rho \eta)_s \right] + \dv\Big( (\rho \eta) \veb{v} - \sum_{i=1}^N \frac{\mu_i}{T} \veb{J}_i + \frac{\veb{q}}{T} \Big) \geq 0.
\label{EBI}
\end{align}
\end{remark}
 The main a priori estimates are given by the following proposition.
\begin{lemma}\label{lem.ap.est}
Any smooth solution to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}--\eqref{bc} satisfies
\begin{align}
\|P_c(S)\|_{L^1(Q_\mathcal{T})} + 
\|p\|_{L^1(Q_\mathcal{T})} +
\|\rho\|_{L^\gamma(Q_\mathcal{T})} &\leq C,\label{Est.p}\\
\|S^{1/\gamma}\rho\|_{L^\infty(0,\mathcal{T}; L^\gamma(\Omega))} + 
\|S p\|_{L^\infty(0,\mathcal{T}; L^1(\Omega))} + 
\|\Phi E_f + (1-\Phi) E_s\|_{L^\infty(0,\mathcal{T}; L^1(\Omega))} &\leq C,
\label{Est.E.LinfL1}\\
\|T\|_{L^\infty(0,\mathcal{T}; L^1(\Omega))} + 
 \|\log T \|_{L^2(0,\mathcal{T}; H^1(\Omega))} 
 + \|T^{\beta/2} \|_{L^2(0,\mathcal{T}; H^1(\Omega))} &\leq C, \label{Est.A}\\
 \label{Est.T}
 \|T\|_{L^{ \beta + \frac{2}{3} }(Q_\mathcal{T})} &\leq C,\\
 \|\Pi^N \big( \vec\mu/T \big)\|_{L^2(0,\mathcal{T}; H^1(\Omega))} 
 +\|\Pi^N \big( \vec\mu/T \big)\|_{L^a(Q_\mathcal{T})} &\leq C,
\label{Est.B.2}\\
\left\|\sqrt{\frac{\lambda(S,T)}{T}}\nabla p\right\|_{L^2(Q_\mathcal{T})}
&\leq C,\label{AJPP.1}\\
\|T^{-1}\|_{L^1(\pa\Omega\times (0,\mathcal{T}))} & \leq C,
\label{T.boundary}\\
\label{mu.boundary}
\left\| \sum_{i,j=1}^N b_{ij}\frac{\mu_i}{T}\frac{\mu_j}{T} \right\|_{L^1(\pa\Omega\times (0,\mathcal{T}))} &\leq C,\\
\|T^{-1} f'(S) \big( \partial_t S \big)^2\|_{L^1(Q_\mathcal{T})} +
\|f(S)\|_{L^\infty(0,\mathcal{T}; W^{1,q}(\Omega))} &\leq C, \label{Est.S}
\end{align}
where here and in the following, $C > 0$ denotes a generic constant depending only on the given data and (possibly) on $\mathcal{T}$.
\end{lemma}

\begin{proof}
We start by integrating the energy balance \eqref{cons_e} in space and time:
\begin{align*}
\int_\Omega\big (
\Phi E_f + (1-\Phi) E_s \big)\vert_{t=0}^{t=t_1}dx
+ \int_0^{t_1}\int_{\pa\Omega}
\big( ((\rho e) + p) \veb{v}  + \veb{q} \big)\cdot\nu d\sigma dt = 0,
\quad\forall t_1 \in (0,\mathcal{T}).
\end{align*}
From \eqref{bc} it follows
\begin{align*}
\int_\Omega\big (
\Phi E_f + (1-\Phi) E_s \big)\vert_{t=0}^{t=t_1}dx
+ \alpha\int_0^{t_1}\int_{\pa\Omega} (T-T_0)
 d\sigma dt = 0,\quad\forall t_1 \in (0,\mathcal{T}),
\end{align*}
and so
\begin{equation}
\label{est.totE}
\sup_{t\in [0,\mathcal{T}]}\int_\Omega\big (
\Phi E_f + (1-\Phi) E_s \big)dx
+ \alpha\int_0^\mathcal{T}\int_{\pa\Omega} T d\sigma dt\leq C .
\end{equation}
The $L^\infty(0,\mathcal{T}; L^\gamma(\Omega))$ bound for $S^{1/\gamma}\rho$
and the $L^\infty(0,\mathcal{T}; L^1(\Omega))$ bound for $Sp$
are derived from the constitutive relations 
\eqref{def.Ef}, \eqref{thp.1}, \eqref{eng.1} and estimate \eqref{est.totE}.
Therefore \eqref{Est.E.LinfL1} holds. As a byproduct,
from \eqref{Skel.Ener}, \eqref{Phi.pos}, \eqref{Est.E.LinfL1}
it follows
\begin{align}\label{T.LinfL1}
	\|T\|_{L^\infty(0,\mathcal{T}; L^1(\Omega))}\leq C.
\end{align}
Let $\zeta : [0,1]\to [0,1]$ be a $C^1([0,1])$ 
cutoff function such that 
$\zeta = 1$ on $[0,1/3]$, $\zeta = 0$ on $[2/3,1]$, and
$\zeta$ is decreasing in $(1/3,2/3)$.
Let us multiply \eqref{PC.Sw} by $\zeta(S)$ and define
$f_1(s) = \int_0^s\zeta(u)f'(u)du$ for $0\leq s\leq 1$.
Notice that $f_1$ is bounded by construction and the assumptions on $\zeta$. We obtain
\begin{align*}
\pa_t f_1(S) + \zeta(S)P_c(S) + \zeta(S)p = 0.
\end{align*}
Integrating the above identity leads to
\begin{align*}
\int_\Omega f_1(S(t))dx + \int_0^t\int_\Omega(\zeta(S)P_c(S) + \zeta(S)p)dx dt' = 
\int_\Omega f_1(S^{in})dx,\quad t>0.
\end{align*}
Since $p\geq -p_{at}$ a.e.~in $Q_\mathcal{T}$ and $P_c$ is a nonnegative function while $f_1$ is bounded, we obtain
that $\zeta(S)P_c(S)$ is bounded in $L^1(Q_\mathcal{T})$. Since 
$(1-\zeta)P_c$ is bounded (because $P_c$ is nonnegative, decreasing, and $1-\zeta$ vanishes near $0$), we infer
\begin{align}\label{int.PcS.ub}
\int_0^t\int_\Omega P_c(S)dx dt\leq C.
\end{align}
On the other hand, multiplying \eqref{PC.Sw} times $1-\zeta(S)$ and defining 
$f_2(s) = \int_0^s(1-\zeta(u))f'(u)du$ yields
\begin{align*}
\pa_t f_2(S) + (1-\zeta(S))P_c(S) + (1-\zeta(S))p = 0.
\end{align*}
Integrating the above identity leads to
\begin{align*}
-\int_\Omega f_2(S(t))dx =
\int_0^t\int_\Omega((1-\zeta(S))P_c(S) + 
(1-\zeta(S))p)dx dt' - 
\int_\Omega f_2(S^{in})dx,\quad t>0.
\end{align*}
We know that $(1-\zeta)P_c$ is bounded.
Moreover, the properties of $\zeta$, the constitutive relations \eqref{def.Ef}, \eqref{thp.1}, \eqref{eng.1} and \eqref{Est.E.LinfL1} imply that $(1-\zeta(S))p$ is bounded in $L^1(Q_\mathcal{T})$. It follows that
\begin{align*}
-\int_\Omega f_2(S(t))dx\leq C,\quad t>0.
\end{align*}
Since $-f_2(s) = -f_2(2/3)+f(2/3) - f(s)$ for $s>2/3$ and the fact that $f$ is bounded in $[0,2/3]$ (because it is decreasing, smooth in $(0,1)$, and $f(0)<\infty$), we conclude
\begin{align}\label{int.fS.lb}
-\int_\Omega f(S(t))dx\leq C,\quad t>0.
\end{align}
Integrating \eqref{PC.Sw} in $Q_\mathcal{T}$ and employing \eqref{int.PcS.ub}, \eqref{int.fS.lb} and the fact that $p\geq -p_{at}$ a.e.~in $Q_\mathcal{T}$ yields \eqref{Est.p}.

Let us now integrate the entropy balance equation \eqref{EBE} with respect to space and time (from $t=0$ to $t=\mathcal{T}$) and employ \eqref{ic.1}--\eqref{bc}. We get
 \begin{align}
  \int_0^{\mathcal{T}} \int_\Omega 
  &\Big( \sum_{i,j = 1}^N  L_{ij} \nabla \big(  \frac{\mu_i}{T}\big) \cdot \nabla \big(   \frac{\mu_j}{T} \big) + \kappa(T)|\nabla \log T|^2   \Big) dx dt \nonumber \\
  & + \int_0^{\mathcal{T}}\int_\Omega
  \frac{K}{T}\lambda(S,T)|\nabla p|^2 dx dt
  -\int_0^\mathcal{T}\int_\Omega\sum_{i=1}^N r_i\frac{\mu_i}{T}dx dt\nonumber\\
  & +\int_0^{\mathcal{T}}\int_{\pa\Omega}\left(\alpha\frac{T_0 - T}{T}
  +\sum_{i,j=1}^N b_{ij}\frac{\mu_i}{T}\left(\frac{\mu_j}{T} - \frac{\mu_{0,j}}{T_0}\right)
   \right)d\sigma dt
  \nonumber\\
  & = \int_\Omega \Big(  \Phi S(\rho\eta)  
  + (1-\Phi)(\rho \eta)_s   \Big) dx\vert_{t=0}^{t=\mathcal{T}} ,
 \label{EBE.1}
 \end{align}
 where $\kappa(T) $ is given by \eqref{heat_cond}.
Let us consider the right-hand side of \eqref{EBE.1}. 
From \eqref{entropy} and \eqref{Skel.Entr} we deduce
\begin{align*}
\int_\Omega \Big(  \Phi S(\rho\eta)  
+ (1-\Phi)(\rho \eta)_s   \Big) dx\vert_{t=0}^{t=\mathcal{T}} \leq 
C_1 +
C\int_\Omega\left( 
\Phi S \rho (\log T)_+ + (1-\Phi)(\log T)_+
 \right)  dx\vert_{t=\mathcal{T}} ,
\end{align*}
where $C_1$ depends on the initial datum and $(u)_+:=\max \{u,0\}$ denotes the positive part.
Since \eqref{def.Ef}, \eqref{eng.1}, \eqref{Skel.Ener} hold, 
we deduce
\begin{align*}
\int_\Omega \Big(  \Phi S(\rho\eta) 
+ (1-\Phi)(\rho \eta)_s   \Big) dx\vert_{t=0}^{t=\mathcal{T}} \leq 
C_1 +
C\int_\Omega \left(\Phi E_f + (1-\Phi)E_s \right) dx\vert_{t=\mathcal{T}} .
\end{align*}
From \eqref{est.totE} we get
\begin{equation}
\label{est.eta}
\int_\Omega \Big(  \Phi S(\rho\eta) 
+ (1-\Phi)(\rho \eta)_s   \Big) dx\vert_{t=0}^{t=\mathcal{T}} \leq C.
\end{equation}
Let us now consider the boundary term in \eqref{EBE.1}. Since the matrix $(b_{ij})_{i,j=1}^N$ is positive semidefinite and \eqref{hp.b} holds, it follows (via Cauchy-Schwartz)
\begin{align*}
\sum_{i,j=1}^N b_{ij}\frac{\mu_i}{T}\left(\frac{\mu_j}{T} - \frac{\mu_{0,j}}{T_0}\right) &\geq 
\frac{1}{2}\sum_{i,j=1}^N b_{ij}\frac{\mu_i}{T}\frac{\mu_j}{T} 
-\frac{1}{2}\sum_{i,j=1}^N b_{ij}\frac{\mu_{0,i}}{T_0}\frac{\mu_{0,j}}{T_0},
\end{align*}
which yields \eqref{mu.boundary}. 
Furthermore, \eqref{hp.r.2} implies
\begin{equation*}
-\int_0^\mathcal{T}\int_\Omega\sum_{i=1}^N r_i\frac{\mu_i}{T}dx dt
\geq C_1\|\Pi^N(\frac{\vec\mu}{T})\|_{L^a(Q_\mathcal{T})}^a - C_0,
\end{equation*}
which yields the $L^a(Q_\mathcal{T})$ bound for $\Pi^N\vec\mu/\mathcal{T}$.
On the other hand, \eqref{Ass.L} implies that 
\begin{align}
 &\sum_{i,j=1}^N L_{ij} \nabla \big( \frac{\mu_i}{T} \big) \cdot
  \nabla \big( \frac{\mu_j}{T} \big)  = \sum_{i,j=1}^N L_{ij} \sum_{k=1}^N \partial_k \big( \frac{\mu_i}{T} \big)  \partial_k \big( \frac{\mu_i}{T} \big) \nonumber \\
   & = \sum_{k=1}^N \Big(\sum_{i,j=1}^N L_{ij} \partial_k \big( \frac{\mu_i}{T} \big)  \partial_k \big(\frac{\mu_j}{T} \big) \Big) \geq C \sum_{k=1}^N |\Pi^N \partial_k \big( \frac{\vec\mu}{T}  \big)|^2 \nonumber \\
   & = C \sum_{k=1}^N \sum_{i=1}^N \Big( \Pi^N \partial_k \big(   \frac{\vec\mu}{T}\big)  \Big)_i^2 
    = C \sum_{i=1}^N |\nabla \Pi^N \big( \frac{\vec\mu}{T} \big)_i |^2
     = C |\nabla \Pi^N \big( \frac{\vec\mu}{T}  \big) |^2.
\end{align}
In this way from \eqref{EBE.1} we get the estimate
\begin{align*}
 \| \nabla \Pi^N \big( \frac{\vec\mu}{T} \big) \|_{L^2(Q_T)} \leq C,
\end{align*}
which, together with the $L^a(Q_\mathcal{T})$ bound for $\Pi^N\vec\mu/\mathcal{T}$, yields \eqref{Est.B.2}.

Clearly \eqref{T.boundary} follows immediately from \eqref{EBE.1} and the previous bounds for the boundary integrals and \eqref{est.eta}.
On the other hand, from \eqref{kappa} it follows
\begin{align} 
	\kappa(T)|\nabla \log T|^2 & \geq \kappa_1(1+T^\beta) |\nabla \log T|^2 \geq c(|\nabla \log T|^2 + |\nabla T^{\beta/2}|^2 ).  \label{Eq.15}
\end{align}
Therefore \eqref{T.boundary}, \eqref{T.LinfL1}, \eqref{EBE.1}, \eqref{Eq.15} lead to the uniform bounds for $T$, $\log T$ in 
$L^\infty(0,T; L^1(\Omega))$, $L^2(0,T; H^1(\Omega))$ (respectively) contained in \eqref{Est.A}. In order to prove the uniform bound for $T^{\beta/2}$ in $L^2(0,T; H^1(\Omega))$ (which will yield \eqref{Est.A}), it is necessary to distinguish two cases.\\
{\em Case $1< \beta\leq 2$.} In this case $T^{\beta/2}$ is bounded in $L^\infty(0,T; L^1(\Omega))$ thanks to \eqref{T.LinfL1}, so Poincar\'e's Lemma yields the statement.\\
{\em Case $\beta>2$.} In this case we observe that
\begin{align*}
	\nabla T = \chf{(0,1)}(T) T \nabla\log T + \chf{[1,\infty)}(T)\frac{2}{\beta}T^{1-\beta/2}\nabla T^{\beta/2}
\end{align*}
is bounded in $L^2(Q_T)$ given the uniform bounds for $\nabla\log T$ and $\nabla T^{\beta/2}$ in $L^2(Q_T)$. Given that \eqref{T.LinfL1} holds, it follows via Poincar\'e's Lemma that $T$ is bounded in $L^2(0,T; H^1(\Omega))$ and therefore (via Sobolev embedding) in $L^2(0,T; L^6(\Omega))$. We apply now an iterative interpolation argument to prove that $T^{\beta/2}$ is bounded in $L^2(0,T; L^1(\Omega))$; the statement will then follow via Poincar\'e's Lemma. Assume that, for some $p\geq 2$,
\begin{align}\label{T.LpL3p}
\|T\|_{L^p(0,T; L^{3p}(\Omega))}\leq C .
\end{align}
Since \eqref{T.LinfL1} holds, a straightforward interpolation argument yields
\begin{align*}
\|T\|_{L^{2r}(0,T; L^{r}(\Omega))}\leq C,\quad r = \frac{5}{6} + \frac{p}{2}.
\end{align*}
This means that $T^r$ is bounded in $L^2(0,T; L^1(\Omega))$. If $\beta\leq 2r$ then we deduce (easy interpolation with \eqref{T.LinfL1}) that $T^{\beta/2}$ is bounded in $L^2(0,T; L^1(\Omega))$, which is what we want. Let us therefore assume that $\beta>2r$. It follows from the known bounds for $\nabla\log T$, $\nabla T^{\beta/2}$ in $L^2(Q_T)$ that 
$\nabla T^{r}$ is bounded in $L^2(Q_T)$. Since $T^r$ is bounded in $L^2(0,T; L^1(\Omega))$, Poincar\'e's Lemma yields that $T^r$ is bounded in $L^2(0,T; H^1(\Omega))$, thus Sobolev's embedding implies that $T^r$ is bounded in $L^2(0,T; L^6(\Omega))$, that is, \eqref{T.LpL3p} holds with $p$ replaced by $2r = \frac{5}{3} + p$. Since \eqref{T.LpL3p} holds true with $p=2$, an easy induction argument implies that \eqref{T.LpL3p}
holds for every $p$ having the form $p = 2 + \frac{5}{3}n$, with $n\in\N$ arbitrary, as long as $2 + \frac{5}{3}n < \frac{\beta}{2}$. Since $2 + \frac{5}{3}n\to\infty$ as $n\to\infty$, iterating this argument allows us to eventually deduce that $T^{\beta/2}$ is bounded in $L^2(0,T; L^1(\Omega))$.
Bound \eqref{AJPP.1} follows immediately from \eqref{EBE.1}.

Let us now show \eqref{Est.S}. 
The bound for the first term in \eqref{Est.S} follows immediately from \eqref{EBE.1}; let us try to bound the second term.
Take the gradient of \eqref{PC.Sw}, and multipy the resulting equation times $|\nabla f(S)|^{q-2}\nabla f(S)$ and integrate in $\Omega\times [0,t_1]$, where $t_1\in [0,\mathcal{T}]$ is generic. It follows
\begin{align*}
&\frac{1}{q}\int_\Omega |\nabla f(S)|^q dx\Big\vert_{t=0}^{t=t_1} 
+ \int_0^{t_1}\int_\Omega |f'(S)|^{q-2} f'(S)P_c'(S) |\nabla S|^q dx dt\\
&\qquad
= -\int_0^{t_1}\int_\Omega |\nabla f(S)|^{q-2}\nabla f(S)\cdot\nabla p dx dt .
\end{align*}
Young's inequality and \eqref{AJPP.1} lead to
\begin{align*}
&\frac{1}{q}\int_\Omega |\nabla f(S)|^q dx\Big\vert_{t=0}^{t=t_1} 
+ \int_0^{t_1}\int_\Omega |f'(S)|^{q-2}f'(S)P_c'(S) |\nabla S|^q dx dt\\
&\quad\leq \frac{\eps}{2}\int_0^{t_1}\int_\Omega\frac{T}{\lambda} 
|\nabla f(S)|^{2(q-1)} dx dt + 
\frac{1}{2\eps}
\int_0^{t_1}\int_\Omega\frac{\lambda}{T}|\nabla p|^2 dx dt\\
&\quad\leq \frac{\eps (q-1)}{q}\int_0^{t_1}\int_\Omega\lambda^{-q/(2q-2)} |\nabla f(S)|^q dx dt + 
\frac{\eps (2-q)}{2q}\int_0^{t_1}\int_\Omega T^{q/(2-q)} dx dt 
+ C\eps^{-1} .
\end{align*}
Thanks to \eqref{mobility} and the boundedness of $\mu$, it follows
$1/\lambda = \mu(T)/k_r(S)\leq C/k_r(S)$. Furthermore \eqref{hp.beta}, 
\eqref{Est.A} yield that $\|T^{q/(2-q)}\|_{L^1(Q_\mathcal{T})}\leq C$.
From these estimates and the fact that both $P_c$ and $f$ are nonincreasing, we deduce
\begin{align*}
&\frac{1}{q}\int_\Omega |\nabla f(S)|^q dx\Big\vert_{t=0}^{t=t_1} 
+ \int_0^{t_1}\int_\Omega |f'(S)|^{q-1} |P_c'(S)| |\nabla S|^q dx dt\\
&\qquad\leq C\eps\int_0^{t_1}\int_\Omega
k_r(S)^{-q/(2q-2)}|f'(S)|^q |\nabla S|^q dx dt + C\eps^{-1} .
\end{align*}
Since $k_r$ is uniformely positive in $[1/2,1]$, it follows
\begin{align*}
&\frac{1}{q}\int_\Omega |\nabla f(S)|^q dx\Big\vert_{t=0}^{t=t_1} 
+ \int_0^{t_1}\int_\Omega |f'(S)|^{q-1} |P_c'(S)| |\nabla S|^q dx dt\\
&\quad\leq 
C\eps\int_\Omega |\nabla f(S(t_1))|^q dx\\
&\qquad + C\eps\int_{Q_{t_1}\cap\{ S<1/2 \} }
k_r(S)^{-q/(2q-2)}|f'(S)|^q |\nabla S|^q dx dt + C\eps^{-1} .
\end{align*}
Applying \eqref{hp.Pcf} with $\eps< c_{f}/qC$ 
leads to
\begin{align*}
\int_\Omega |\nabla f(S)|^q dx\Big\vert_{t=0}^{t=t_1} 
	\leq C + 
	C\int_\Omega |\nabla f(S(t_1))|^q dx ,\quad \forall t_1\in [0,\mathcal{T}],
\end{align*}
which via a Gronwall argument yields a bound for $\nabla f(S)$ in $L^\infty(0,\mathcal{T}; L^q(\Omega))$. However, since \eqref{int.fS.lb} holds and $f$ is upper bounded, \eqref{Est.S} follows via Poincar\'e's Lemma.
This finishes the proof of Lemma~\ref{lem.ap.est}.
\end{proof}

The following three propositions (Prop. \ref{prop.5}, Prop. \ref{prop.6}, Prop. \ref{prop.7}) give  a priori estimates needed for the existence result.
Their proofs are rather technical and they are given in the Appendix.
The first one gives estimates for the terms in \eqref{cons_w}, \eqref{cons_e}.
\begin{proposition}\label{prop.5}
Any smooth solution to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}--\eqref{bc} satisfies
\begin{align}
& \|\rho\veb v\|_{L^{\frac{2\beta}{\beta+1}}(0,\mathcal{T}; L^{\frac{6\beta\gamma}{\gamma + 6\beta + 3\beta\gamma}}(\Omega))}\leq C. \label{Est.rhov} \\
& \|\veb J_i\|_{L^2(Q_\mathcal{T})}\leq C,\quad i=1,\ldots,N. \label{Est.J} \\
& \|\rho_i\veb v + \veb J_i\|_{L^m(Q_\mathcal{T})} \leq C,\quad i=1,\ldots,N. \label{Est.flux.1} \\
 & \|\pa_t(\Phi S\rho_i)\|_{L^m(0,\mathcal{T}; W^{-1,m}(\Omega))} \leq C,
\label{Est.Srhot} \\
& \|\pa_t f(S)\|_{L^1(Q_\mathcal{T})}\leq C. \label{Est.St} \\
& \|\pa_t F(S)\|_{L^\infty(0,\mathcal{T}; L^1(\Omega))}\leq C,\quad
F(s) \equiv -\int_s^1 s_1 f'(s_1)ds_1, \label{Est.FSt}
\end{align}
where $\gamma$ and $\beta$ were given by \eqref{gamma}, \eqref{hp.beta} and
$$
m\equiv\min\left\{ \frac{2\beta}{\beta+1}, \frac{6\beta\gamma}{\gamma + 6\beta + 3\beta\gamma} \right\} > 1.
$$
\end{proposition}


The following proposition gives the bound of the total entropy $\Phi S(\rho\eta) + (1-\Phi)(\rho\eta)_s$ and for the entropy flux.
\begin{proposition}\label{prop.6}
Any smooth solution to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}--\eqref{bc} satisfies
\begin{equation}
\label{Est.entr}
\|\Phi S(\rho\eta) + (1-\Phi)(\rho\eta)_s\|_{L^{\frac{2\gamma}{\gamma+2}}(Q_\mathcal{T})} \leq C.
\end{equation}
\begin{equation}
\label{Est.entrflux}
\exists s>1:\quad 
\left\|(\rho\eta)\veb v -\sum_{i=1}^N\frac{\mu_i}{T}\veb J_i +
\frac{\veb q}{T}\right\|_{L^s(Q_\mathcal{T})} \leq C .
\end{equation}
\end{proposition}


Finally, we derive the gradient bound for $\rho$, the bound for $\log \rho_i/\rho$ and the bound for the reaction terms $r_i$, $i=1,\ldots,N$.
\begin{proposition}\label{prop.7}
Any smooth solution to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}--\eqref{bc} satisfies
\begin{equation}
\label{Est.narho}
\|\nabla [\sqrt{k_r(S)} G(\rho^\gamma)]\|_{L^{a_2}(0,\mathcal{T}; L^q(\Omega))}\leq C[G],\qquad \forall G\in W^{1,\infty}(\R_+).
\end{equation}
\begin{equation}
\label{est.relchp}
\left\| \log\frac{\rho_i}{\rho} \right\|_{L^2(Q_\mathcal{T})}\leq C,\qquad
i=1,\ldots,N.
\end{equation}
\begin{equation}
\label{Est.r}
\|r_i\|_{L^{\frac{a}{a-1}}(Q_\mathcal{T})} \leq C,\quad i=1,\ldots,N.
\end{equation}
\end{proposition}

\section{Weak sequential stability for smooth solutions}\label{sec.weakstab}
 The concept of weak sequential stability in the sense of Feireisl \cite{Feir04}, together with weak compactness results, constitute an important step in the global existence analysis.  Generally speaking, weak sequential stability means that, given a sequence of (smooth) solutions to a PDE system of interest, there exists a subsequence which converges to a (weak) solution to the problem. Weak stability implies that the limiting solution is a solution to the original system.

In our case, the sequence of solutions solves an approximate and regularized system of equations, obtained
from the original one by a semidiscretization in time,
and the index of the sequence is related to the approximation and regularization parameters. In order to prove weak compactness, the Div-Curl lemma \cite[Prop. 3.3]{FeirNov09}, developed by Murat \cite{Murat78} and Tartar \cite{Tartar75}, has been used, which represents an efficient tool for handling compactness in nonlinear problems, where the classical Rellich--Kondraschev argument is not applicable.

In this section we assume that a sequence $(\vec{\rho}^{(n)},T^{(n)},S^{(n)})$ of variational entropy  solutions to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}, \eqref{bc} exists, and show that $(\vec{\rho}^{(n)},T^{(n)},S^{(n)})$ converges (up to subsequences) to another variational entropy solution $(\vec{\rho},T,S)$ of \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}, \eqref{bc} as $n\to \infty$.

\begin{theorem}[Weak stability]\label{thm.ws}
Let $(\vec{\rho}^{(n)},T^{(n)},S^{(n)})$ be a sequence of variational entropy solutions to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}, \eqref{bc} according to Def.~\ref{def.weaksol}.
Assume furthermore that $\rho^{(n)}_i, T^{(n)}, S^{(n)} > 0$ a.e.~in $Q_\mathcal{T}$ for $i=1,\ldots,N$, $n\in\N$. 
Then $(\vec{\rho}^{(n)},T^{(n)},S^{(n)})$ is (up to subsequences) strongly convergent in $L^1(Q_\mathcal{T})$ (as $n\to \infty$) to some function triplet
$(\vec{\rho},T,S)$ which is a variational entropy solution to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1}, \eqref{bc} according to Def.~\ref{def.weaksol}.	
\end{theorem}

We shall prove Theorem \ref{thm.ws} in several steps described in the following subsections.
\subsection{Convergence results}
  In the following lemma we will collect the strong convergence results, $\forall\eps>0$, needed for passing to the limit in \eqref{weak.1}--\eqref{weak.2}--\eqref{weak.EBE}, when $n \to \infty$. 
\begin{lemma}\label{str.convergences}
The following strong convergence results hold:
\begin{equation}
\label{S.cnv}
S^{(n)}\to S\quad\mbox{strongly in }L^r(Q_\mathcal{T}),~~\forall r<\infty,
\end{equation}
where 
\begin{equation}
\label{S.pos}
 S>0 \quad \textrm{a.e. in }\quad Q_\mathcal{T}.
\end{equation}
\begin{equation}
	\label{p.cnv.Qeps}
p^{(n)}\to p\quad\mbox{strongly in $L^{4/3}(\Omega\times [\eps,\mathcal{T}])$ for every $\eps>0$,} 
\end{equation}
\begin{equation}
\label{p.cnv}
S^{(n)} p^{(n)}\to S p\quad\mbox{strongly in }L^1(Q_\mathcal{T}).
\end{equation}
\begin{equation}
\label{rho.cnv}
\rho^{(n)}\to \rho\quad\mbox{strongly in }L^{\gamma-\eps}(Q_\mathcal{T}),
\end{equation}
\begin{equation}
\label{rhoi.cnv}
\rho^{(n)}_i\to\rho_i\quad\mbox{strongly in }L^{\gamma-\eps}(Q_\mathcal{T}),\quad
i=1,\ldots,N,
\end{equation}
\begin{align}
	\label{cnv.aux}
S^{(n)}\rho^{(n)}_i\log\rho^{(n)}_i\to S\rho_i\log\rho_i,\quad
S^{(n)}\rho^{(n)}_i\to S\rho_i,\quad i=1,\ldots,N,\\
\nonumber
\mbox{strongly in }L^{\gamma-\eps}(Q_\mathcal{T}),
\end{align}
\begin{align}
	\label{T.cnv}
\logT^{(n)}\to\log T\quad\mbox{strongly in }L^{2-\eps}(0,\mathcal{T}; L^{6-\eps}(\Omega)),\\
\label{T.cnv.2}
T^{(n)}\to T\quad\mbox{strongly in }L^{\beta + \frac{2}{3} - \eps}(Q_\mathcal{T}).
\end{align}
\end{lemma}
\begin{proof}
This proof will be divided into 4 steps. In Step 1 we will prove \eqref{S.cnv}, \eqref{S.pos} and \eqref{p.cnv}. In Step 2 we will prove \eqref{rho.cnv}. Step 3 will be about proving \eqref{rhoi.cnv} and \eqref{cnv.aux}, and finally Step 4 proves \eqref{T.cnv} and \eqref{T.cnv.2}.   
Let us define preliminarly the following vector fields:
\begin{equation}
	\label{Un.i}
U_i^{(n)} = (\Phi S^{(n)}\rho^{(n)}_i, \rho^{(n)}_i\mathbf{v}^{(n)} + \mathbf{J}^{(n)}_i),\qquad
i=1,\ldots,N,
\end{equation}
\begin{equation}
\label{Un}
U^{(n)} = \sum_{i=1}^N U^{(n)}_i = (\Phi S^{(n)}\rho^{(n)}, \rho^{(n)}\mathbf{v}^{(n)}),
\end{equation}
\begin{equation}
\label{Vn}
V^{(n)}[G] = (\sqrt{k_r(S^{(n)})}G( (\rho^{(n)})^{\gamma}), 0,0,0),\qquad
G\in W^{1,\infty}(\R_+),
\end{equation}
\begin{equation}
\label{Zn}
Z^{(n)}[G] = ( G(\Pi^N\vec\mu^{(n)}/T^{(n)}),0,0,0),\qquad
G\in W^{1,\infty}(\R^N),
\end{equation}
\begin{equation}
\label{Wn}
W^{(n)} = (\PhiS^{(n)}(\rho\eta)^{(n)} + (1-\Phi)(\rho\eta)_s^{(n)},  (\rho\eta)^{(n)} \veb{v}^{(n)} - \sum_{i=1}^N \frac{\mu_i^{(n)}}{T^{(n)}} \veb{J}_i^{(n)} + \frac{\veb{q}^{(n)}}{T^{(n)}} ),
\end{equation}
\begin{equation}
\label{Yn}
Y^{(n)}[G] = ( G(T^{(n)}),0,0,0),\qquad G\in W^{1,\infty}(\R_+).
\end{equation}


{\em \textbf{Step 1: Strong convergence, a.e. positivity of saturation and convergence of $S^{(n)} p^{(n)}$.}} \\
In this step we show \eqref{S.cnv}--\eqref{p.cnv}.
From \eqref{Est.S}, \eqref{Est.St} it follows
\begin{align*}
\|\pa_t f(S^{(n)})\|_{L^1(0,\mathcal{T}; L^1(\Omega))} + 
\|f(S^{(n)})\|_{L^\infty(0,\mathcal{T}; W^{1,q}(\Omega))}\leq C.
\end{align*}
Since $W^{1,q}(\Omega)\hookrightarrow L^q(\Omega)$ compactly and 
$L^q(\Omega)\hookrightarrow L^1(\Omega)$ continuously, the Aubin-Lions Lemma
allows us to deduce that, up to subsequences,
\begin{align*}
f(S^{(n)})\quad\mbox{is strongly convergent in $L^q(Q_\mathcal{T})$ as $n\to \infty$.}
\end{align*}
In particular $f(S^{(n)})$ is a.e.~convergent in $Q_\mathcal{T}$. Since $f$ is strictly decreasing (and a fortiori one-to-one), it follows that $S^{(n)}$
is a.e.~convergent in $Q_\mathcal{T}$. Thanks to the uniform $L^\infty(Q_\mathcal{T})$ bounds for $S^{(n)}$, we deduce \eqref{S.cnv}.

We now prove \eqref{S.pos}. Since $p^{(n)}\geq -p_{at}$
and \eqref{hp.Pcf2} holds, it follows
\begin{align*}
\pa_t f(S^{(n)}) + \Gamma f(S^{(n)}) \leq p_{at},\quad
\Gamma\equiv \inf_{0<s<s_0}\frac{P_c(s)}{f(s)}>\frac{p_{at}}{f(0)}.
\end{align*}
A Gronwall argument yields
\begin{align*}
f(S^{(n)}(t))\leq f(S^{in})e^{-\Gamma t} + \frac{p_{at}}{\Gamma}(1-e^{-\Gamma t})
,\quad t>0,~~ \mbox{a.e.~in }\Omega .
\end{align*}
Relation \eqref{S.cnv} allows us to take the limit $n \to\infty$ in 
the above inequality and obtain
\begin{align}\label{ub.fS}
f(S(t))\leq f(S^{in})e^{-\Gamma t} + \frac{p_{at}}{\Gamma}(1-e^{-\Gamma t}) 
,\quad t>0,~~ \mbox{a.e.~in }\Omega .
\end{align}
Consider the set $E := \{ x\in\Omega~:~ f(S^{in}(x))\geq p_{at}/\Gamma \}$. Clearly $\{S^{in}=0\}\subset E$ since 
$\Gamma>p_{at}/f(0)$. Moreover
\begin{align*}
f(S(t))\leq \frac{p_{at}}{\Gamma},\quad t>0,\quad\mbox{a.e.~in }\Omega\backslash E ,
\end{align*}
which means that 
\begin{align*}
S(t)\geq f^{-1}\left( \frac{p_{at}}{\lambda} \right) > 0,\quad t>0,\quad\mbox{a.e.~in }\Omega\backslash E .
\end{align*}
On the other hand, the right-hand side of \eqref{ub.fS} is decreasing in time for $x\in E$, so
\begin{align*}
\forall\eps>0,\quad
f(S(t))\leq f(S^{in})e^{-\lambda \eps} + \frac{p_{at}}{\lambda}(1-e^{-\lambda \eps}),
\quad t\geq\eps,~~ \mbox{a.e.~in }E .
\end{align*}
But $f(S^{in}(x))\leq f(0)$, so
\begin{align*}
\forall\eps>0,\quad
S(t)\geq \tilde s_\eps\equiv f^{-1}\left( 
f(0)e^{-\lambda \eps} + \frac{p_{at}}{\lambda}(1-e^{-\lambda \eps}) \right), \quad t\geq\eps,~~ \mbox{a.e.~in }E .
\end{align*}
Clearly $\tilde s_\eps>0$ for every $\eps>0$. This fact, combined with the uniform positivity of $S$ on $(\Omega\backslash E)\times [0,\mathcal{T}]$, allows us to deduce
\begin{equation}
\label{lb.S}
\forall\eps>0\quad\exists s_\eps>0:\quad
S(t)\geq s_\eps,\quad t\geq \eps,\quad\mbox{a.e.~in }\Omega.
\end{equation}
In particular $S>0$ a.e.~in $Q_\mathcal{T}$ and thus \eqref{S.pos} holds.
Furthermore, it is immediate to see from the argument above that \eqref{lb.S} also holds for every term of the sequence $S^{(n)}$, $n\in\N$. 
\medskip\\
The strong convergence of the saturation \eqref{S.cnv} and its positivity \eqref{S.pos} will help us to  show \eqref{p.cnv}.
Let us define $Q_\mathcal{T}^\eps = \Omega\times [\eps,\mathcal{T}]$.
From \eqref{mobility}, \eqref{AJPP.1}, \eqref{lb.S} 
as well as the uniform boundedness of $\mu$ it follows
\begin{align*}
\forall\eps>0,\quad\exists C_\eps>0:\quad
\|(T^{(n)})^{-1/2}\nabla p^{(n)}\|_{L^2(Q_\mathcal{T}^\eps)}\leq C_\eps .
\end{align*}
From the above estimate, \eqref{Est.T} and the fact that $\beta>2$ we obtain via H\" older's inequality
\begin{align*}
\forall\eps>0,\quad\exists C_\eps>0:\quad
\|\nabla p^{(n)}\|_{L^{16/11}(Q_\mathcal{T}^\eps)}\leq 
\|(T^{(n)})^{1/2}\|_{L^{16/3}(Q_\mathcal{T}^\eps)}
\|(T^{(n)})^{-1/2}\nabla p^{(n)}\|_{L^{2}(Q_\mathcal{T}^\eps)} \leq 
C_\eps .
\end{align*}
Via Poincar\'e Lemma we get
\begin{align*}
\forall\eps>0,\quad\exists C_\eps>0:\quad
\left\|p^{(n)} - |\Omega|^{-1}\int_\Omega p^{(n)} dx\right\|_{L^{16/11}(Q_\mathcal{T}^\eps)}\leq C_\eps .	
\end{align*}
On the other hand, \eqref{lb.S}, \eqref{Est.E.LinfL1} yield
\begin{align*}
\left|\int_\Omega p^{(n)}(t) dx\right|\leq \frac{1}{s_\eps}
\int_{\Omega}S^{(n)}(t) |p^{(n)}(t)|dx\leq \frac{C}{s_\eps},\quad
t\geq \eps,
\end{align*}
which means
\begin{align}\label{Est.p2.a}
\forall\eps>0,\quad\exists C_\eps>0:\quad
\|p^{(n)}\|_{L^{16/11}(Q_\mathcal{T}^\eps)}\leq C_\eps .
\end{align}
Since $p^{(n)}\to p$ a.e.~in $Q_\mathcal{T}$, 
the Vitali-De la Vall\'ee Poussin Theorem \cite[Corollary 10.2]{FeirNov09} implies \eqref{p.cnv.Qeps}. 
Given that \eqref{S.cnv}, \eqref{p.cnv.Qeps} hold, it follows
\begin{align*}
S^{(n)} p^{(n)}\to S p\quad\mbox{strongly in }L^{1}(Q_\mathcal{T}^\eps),\quad\forall\eps>0.
\end{align*}
However, 
\begin{align*}
\|S^{(n)} p^{(n)}\|_{L^{1}(0,\eps; L^1(\Omega))}\leq
\eps\|S^{(n)} p^{(n)}\|_{L^{\infty}(0,\eps; L^1(\Omega))}\leq
C \eps,
\end{align*}
and a similar estimate holds for $S p$ (e.g.~via Fatou's Lemma).
It follows that
\begin{align*}
&\limsup_{n\to\infty}\|S^{(n)} p^{(n)} - S p\|_{L^1(Q_\mathcal{T})}\\
&\leq \limsup_{n\to\infty}\|S^{(n)} p^{(n)} - S p\|_{L^1(Q_\mathcal{T}^\eps)}
+ \limsup_{n\to\infty}\|S^{(n)} p^{(n)} - S p\|_{L^1(Q_\mathcal{T}\backslash Q_\mathcal{T}^\eps)}\\
&\leq \sup_{n\in\N}\left( 
\|S^{(n)} p^{(n)}\|_{L^{1}(0,\eps; L^1(\Omega))} + 
\|S p\|_{L^{1}(0,\eps; L^1(\Omega))}
\right)\leq C\eps .
\end{align*}
This means that \eqref{p.cnv} holds.
\medskip
\noindent

{\em \textbf{Step 2: Strong convergence of the total density.}} In this step we show \eqref{rho.cnv}.\\
Since \eqref{Est.p} holds, it follows that (up to subsequences)
$$
\rho^{(n)}\rightharpoonup\rho\quad\mbox{weakly in }L^\gamma(Q_\mathcal{T}).
$$
Let us consider the vector fields $U^{(n)}$, $V^{(n)}[G]$ defined in \eqref{Un}, \eqref{Vn}, with $G\in W^{1,\infty}(\R_+)$ arbitrary. Thanks to \eqref{Est.E.LinfL1}, \eqref{Est.flux.1}, we deduce that $U^{(n)}$ is bounded in $L^{m}(Q_\mathcal{T})$ with $m>1$, while $V^{(n)}[G]$ is trivially bounded in $L^\infty(Q_\mathcal{T})$. On the other hand,
summing \eqref{cons_w} in $i=1,\ldots,N$ yields $\dv_{(t,x)}U^{(n)}=0$,
while the antisymmetric part of the Jacobian of $V^{(n)}[G]$, which we denote by $\curl_{(t,x)} V^{(n)}[G]=\nabla V^{(n)}[G]-\nabla^T V^{(n)}[G]$, satisfies
$$
|\curl_{(t,x)} V^{(n)}[G]|\leq C |\nabla [\sqrt{k_r(S^{(n)})}G( (\rho^{(n)})^{\gamma})]|
$$
and is therefore bounded in $L^1(Q_\mathcal{T})$ thanks to \eqref{Est.narho}.
In particular both $\dv_{(t,x)}U^{(n)}$, $\curl_{(t,x)} V^{(n)}[G]$ are relatively compact in $W^{-1,z}(Q_\mathcal{T})$ for some $z>1$. Therefore we can apply the Div-Curl Lemma \cite[Prop. 3.3]{FeirNov09} and deduce
\begin{align*}
\overline{U^{(n)}\cdot V^{(n)}[G]} = 
\overline{U^{(n)}}\cdot \overline{V^{(n)}[G]},
\end{align*}
where $\overline{F(U^{(n)})}$ denotes a weak $L^1$-limit of the sequence $(F(U^{(n)}))_{n \in \mathbb{N}}$, meaning that $\overline{F(U^{(n)})} \in L^1(\Omega)$ is defined by
$$ \int_{\Omega}\overline{F(U^{(n)})}\phi\,dx := \lim_{n\to\infty} \int_{\Omega}F(U^{(n)})\phi\,dx \quad \textnormal{for all}~~\phi \in L^{\infty}(\Omega). $$
This is equivalent to
\begin{align*}
\overline{\PhiS^{(n)}\rho^{(n)} \sqrt{k_r(S^{(n)})}G( (\rho^{(n)})^{\gamma})} = 
\overline{\PhiS^{(n)}\rho^{(n)}}\,
\overline{\sqrt{k_r(S^{(n)})}G( (\rho^{(n)})^{\gamma})} .
\end{align*}
However, since $\Phi>0$ and does not depend on $n$, while $S^{(n)}\to S$
strongly in $L^p(Q_\mathcal{T})$ for every $p<\infty$, it follows
\begin{align*}
	S\sqrt{k_r(S)}\,
	\overline{\rho^{(n)} G( (\rho^{(n)})^{\gamma})} = 
	S\sqrt{k_r(S)}\,
	\rho\,\overline{G( (\rho^{(n)})^{\gamma})} .
\end{align*}
It holds $S\sqrt{k_r(S)}>0$ because $S>0$ a.e.~in $Q_\mathcal{T}$.
We deduce
\begin{align}\label{rGr}
\overline{\rho^{(n)} G( (\rho^{(n)})^{\gamma})} = 
\rho\,\overline{G( (\rho^{(n)})^{\gamma})} \quad
\mbox{on }Q_\mathcal{T} , \quad
\forall G\in W^{1,\infty}(\R_+).
\end{align}
Let us now choose $G(s) = G_k(s)\equiv \min\{ (1+s)^{1/\gamma}, (1+k)^{1/\gamma}  \}$,
for $s\geq 0$, $k\in\N$. Define also $G_\infty(s) = (1+s)^{1/\gamma}$.
The weak lower semicontinuity of the $L^1$ norm and Cauchy-Schwartz inequality allow us to estimate
\begin{align*}
&\|
\overline{\rho^{(n)} [G_k( (\rho^{(n)})^{\gamma})-G_\infty( (\rho^{(n)})^{\gamma})] }
\|_{L^1(Q_\mathcal{T})}\\ 
&\quad\leq\liminf_{n\to\infty}\|
\rho^{(n)} [G_k( (\rho^{(n)})^{\gamma})-G_\infty( (\rho^{(n)})^{\gamma})] 
\|_{L^1(Q_\mathcal{T})}\\ 
&\quad\leq\sup_{n\in\N}\|\rho^{(n)}\|_{L^2(Q_\mathcal{T})}
\|G_k( (\rho^{(n)})^{\gamma})-G_\infty( (\rho^{(n)})^{\gamma})\|_{L^2(Q_\mathcal{T})}.
\end{align*}
Thanks to \eqref{Est.p}, it follows (remember that $\gamma>2$):
\begin{align*}
&\|\overline{\rho^{(n)} [G_k( (\rho^{(n)})^{\gamma})-G_\infty( (\rho^{(n)})^{\gamma})] }\|_{L^1(Q_\mathcal{T})}^2\\
&\quad\leq C\sup_{n\in\N}\int_{Q_\mathcal{T}\cap\{\rho^{(n)} > k^{1/\gamma}\}}
(1+(\rho^{(n)})^\gamma)^{2/\gamma}dx dt\\
&\quad\leq C k^{2/\gamma-1}\sup_{n\in\N}\int_{Q_\mathcal{T}\cap\{\rho^{(n)} > k^{1/\gamma}\}}
(\rho^{(n)})^{\gamma-2}(1+(\rho^{(n)})^\gamma)^{2/\gamma}dx dt\\
&\quad\leq C k^{2/\gamma-1}\sup_{n\in\N}\|\rho^{(n)}\|_{L^\gamma(Q_\mathcal{T})}^\gamma
\leq C k^{2/\gamma-1} .
\end{align*}
In a similar way one shows that
\begin{align*}
\|\rho\,
\overline{G_k( (\rho^{(n)})^{\gamma})-G_\infty( (\rho^{(n)})^{\gamma}) }
\|_{L^1(Q_\mathcal{T})}^2\leq C k^{2/\gamma-1}.
\end{align*}
From the above inequalities and \eqref{rGr} with $G=G_k$ it follows
\begin{align*}
&\|\overline{\rho^{(n)} G_\infty( (\rho^{(n)})^{\gamma})} - 
\rho\,\overline{G_\infty( (\rho^{(n)})^{\gamma})}\|_{L^1(Q_\mathcal{T})}\\
&\leq \|
\overline{\rho^{(n)} [G_k( (\rho^{(n)})^{\gamma})-G_\infty( (\rho^{(n)})^{\gamma})] }
\|_{L^1(Q_\mathcal{T})} \\
&\qquad + \|\rho\,
\overline{G_k( (\rho^{(n)})^{\gamma})-G_\infty( (\rho^{(n)})^{\gamma}) }
\|_{L^1(Q_\mathcal{T})}\leq
C k^{1/\gamma-1/2}\to 0\quad\mbox{as }k\to\infty,
\end{align*}
implying that (recall that $G_\infty(s) = (1+s)^{1/\gamma}$)
\begin{align}\label{id.rho}
	\overline{\rho^{(n)} (1 + (\rho^{(n)})^{\gamma})^{1/\gamma}} = 
	\rho\,\overline{(1 + (\rho^{(n)})^{\gamma})^{1/\gamma}} \quad
	\mbox{on }Q_\mathcal{T} .
\end{align}
Since $s\in\R_+\mapsto (1+s^\gamma)^{1/\gamma}\in\R_+$ is strictly increasing and strictly convex, we conclude \cite[Thr.~10.19, Thr.~10.20]{FeirNov09} that, up to subsequences,
$\rho^{(n)}$ is a.e.~convergent in $Q_\mathcal{T}$. Since $\rho^{(n)}$ is bounded in $L^\gamma(Q_\mathcal{T})$, Vitali - De la Vall\'ee Poussin theorem implies that \eqref{rho.cnv} holds. 

\medskip
\noindent
{\em \textbf{Step 3: Strong convergence of partial densities.}}\\
In this step we show \eqref{rhoi.cnv} and \eqref{cnv.aux}.
For every $i=1,\ldots,N$,
the vector field $U_i^{(n)}$ defined in \eqref{Un.i} is bounded in 
$L^m(Q_\mathcal{T})$ thanks to \eqref{Est.flux.1}, while its time-space divergence $\dv_{(t,x)}U_i^{(n)} = r_i^{(n)}$ is bounded in $L^{a/(a-1)}(Q_\mathcal{T})$ due to \eqref{Est.r}, and a fortiori relatively compact in $W^{-1,p}(Q_\mathcal{T})$ for some $p>1$.
On the other hand, for every $G\in W^{1,\infty}(\R^N)$,
the vector field $Z^{(n)}[G]$ defined in \eqref{Zn} is trivially bounded in $L^\infty(Q_\mathcal{T})$,
while \eqref{Est.B.2} implies
\begin{align*}
\|\curl_{(t,x)}Z^{(n)}[G]\|_{L^2(Q_\mathcal{T})}\leq C[G]\|\nabla \Pi^N\vec\mu^{(n)}/T^{(n)} \|_{L^2(Q_\mathcal{T})}\leq C[G].
\end{align*}
In particular $\curl_{(t,x)}Z^{(n)}[G]$ is relatively compact in $H^{-1}(Q_\mathcal{T})$.
Therefore we can apply the Div-Curl Lemma with the vector fields
$U_i^{(n)}$, $Z^{(n)}$ and obtain
\begin{align*}
\overline{U_i^{(n)}\cdot Z^{(n)}} = 
\overline{U_i^{(n)}}\cdot\overline{Z^{(n)}},\qquad i=1,\ldots,N,
\end{align*}
which means
\begin{align*}
\overline{\PhiS^{(n)}\rho^{(n)}_i G(\Pi^N\vec\mu^{(n)}/T^{(n)})} = 
\overline{\PhiS^{(n)}\rho^{(n)}_i}\,
\overline{G(\Pi^N\vec\mu^{(n)}/T^{(n)})},\qquad i=1,\ldots,N.
\end{align*}
Once again, the strong convergence of $S^{(n)}$ and the weak convergence of $\rho^{(n)}_i$, $G(\Pi^N\vec\mu^{(n)}/T^{(n)})$ (as well as the positivity of $\Phi$) imply
\begin{align}\label{riGmuT}
\overline{\rho^{(n)}_i G(\Pi^N\vec\mu^{(n)}/T^{(n)})} = 
\rho_i\,
\overline{G(\Pi^N\vec\mu^{(n)}/T^{(n)})},\quad i=1,\ldots,N,\quad\mbox{on }Q_\mathcal{T} ,
\end{align}
for every $G\in W^{1,\infty}(\R^N)$.
However, from \eqref{chp.1}, \eqref{Pi} it follows
\begin{equation}
\label{relchp}
(\Pi^N\vec\mu^{(n)}/T^{(n)})_i = \log\rho^{(n)}_i - \frac{1}{N}\sum_{j=1}^N\log\rho^{(n)}_j ,\quad i=1,\ldots,N,
\end{equation}
so \eqref{riGmuT} becomes
\begin{align*
	\overline{\rho^{(n)}_i G\left(\log\rho^{(n)}_i - \frac{1}{N}\sum_{j=1}^N\log\rho^{(n)}_j\right)} = 
	\rho_i\,
	\overline{G\left(\log\rho^{(n)}_i - \frac{1}{N}\sum_{j=1}^N\log\rho^{(n)}_j\right)},\quad i=1,\ldots,N,\quad\mbox{on }Q_\mathcal{T} .
\end{align*}
Let $k\in\N$ arbitrary. Define
$ Q_{\mathcal{T},k}\equiv\{ (x,t)\in Q_\mathcal{T} ~:~ \rho(x,t)\geq 1/k \}$.
Since $\rho^{(n)}\to\rho$ strongly in $L^1(Q_\mathcal{T})$, by Egorov-Severini's theorem $\rho^{(n)}\to\rho$ almost uniformly in $Q_{\mathcal{T},k}$, that is
\begin{align*}
\forall\epsilon>0\quad\exists E_{\epsilon,k}\subset Q_{\mathcal{T},k}:\quad
\rho^{(n)}\to\rho ~~ \mbox{strongly in }L^\infty(E_{\epsilon,k}),\quad
|Q_{\mathcal{T},k}\backslash E_\epsilon|<\epsilon .
\end{align*}
We can therefore assume w.l.o.g.~that $\rho^{(n)}\geq 1/2k$ a.e.~in $E_{\epsilon,k}$, $n\in\N$. Defining $\sigma_i^{(n)} = \rho^{(n)}_i/\rho^{(n)}$ on $E_{\epsilon,k}$, for $i=1,\ldots,N$, 
allows us to write
\begin{align*
\sum_{i=1}^N\overline{\rho^{(n)}_i G\left(\log\sigma^{(n)}_i - \frac{1}{N}\sum_{j=1}^N\log\sigma^{(n)}_j\right)} = 
\sum_{i=1}^N\rho_i\,
	\overline{G\left(\log\sigma^{(n)}_i - \frac{1}{N}\sum_{j=1}^N\log\sigma^{(n)}_j\right)}
	\quad\mbox{on }E_{\epsilon,k} .
\end{align*}

Fix $i\in\{1,\ldots,N\}$ generic. By choosing $G(\veb u) = G_M(u_i) = \min\{ u_i, M \}$ in the above identity, with $M\in\N$ generic, exploiting \eqref{Est.B.2} and proceeding like in the proof of \eqref{id.rho} one finds out that
\begin{align*}
\sum_{i=1}^N\overline{\rho^{(n)}_i \left(
	\log\sigma^{(n)}_i - \frac{1}{N}\sum_{j=1}^N\log\sigma^{(n)}_j\right)} = 
\sum_{i=1}^N\rho_i\,
\overline{\left(
	\log\sigma^{(n)}_i - \frac{1}{N}\sum_{j=1}^N\log\sigma^{(n)}_j\right)}\quad
\mbox{on }E_{\epsilon,k},
\end{align*}
which, thanks to \eqref{rho.cnv} and \eqref{est.relchp}, leads to
\begin{align*}
\sum_{i=1}^N\overline{\rho^{(n)}_i \log\sigma^{(n)}_i} = 
\sum_{i=1}^N\rho_i\,
\overline{\log\sigma^{(n)}_i}\quad \mbox{on }E_{\epsilon,k}.
\end{align*}
Once again, \eqref{rho.cnv} and \eqref{est.relchp} imply
\begin{align*}
	\sum_{i=1}^N\overline{\sigma^{(n)}_i \log\sigma^{(n)}_i} = 
	\sum_{i=1}^N\overline{\sigma^{(n)}_i}\,
	\overline{\log\sigma^{(n)}_i}\quad \mbox{on }E_{\epsilon,k} .
\end{align*}
The fact that $\log$ is strictly monotone and strictly concave allows us to conclude that (up to subsequences)

$\sigma^{(n)}_i = \rho^{(n)}_i/\rho^{(n)}$ is a.e.~convergent on $E_{\epsilon,k}$,
which, together with \eqref{rho.cnv}, implies that 
$\rho^{(n)}_i\to \rho_i$ a.e.~in $E_{\epsilon,k}$.
Since $|Q_{\mathcal{T},k}\backslash E_{\epsilon,k}|<\epsilon$ and \eqref{Est.p} holds, we easily deduce that 
$\rho^{(n)}_i\to \rho_i$ strongly in $L^1(Q_{\mathcal{T},k})$ for every $k\in\N$.
Since $0\leq\rho\leq 1/k$ on $Q_\mathcal{T}\backslash Q_{\mathcal{T},k}$, it follows that 
\begin{align*}
\limsup_{n\to\infty}\int_{Q_\mathcal{T}}|\rho^{(n)}_i - \rho_i|dx dt &\leq
\limsup_{n\to\infty}\int_{Q_{\mathcal{T},k}}|\rho^{(n)}_i - \rho_i|dx dt +
\limsup_{n\in\N}\int_{Q_\mathcal{T}\backslash Q_{\mathcal{T},k}}|\rho^{(n)}_i - \rho_i|dx dt\\
&\leq 2 \int_{Q_\mathcal{T}\backslash Q_{\mathcal{T},k}}\rho\, dx dt\leq
\frac{2}{k}|Q_\mathcal{T}|,
\end{align*}
and therefore $\rho^{(n)}_i\to \rho_i$ strongly in $L^1(Q_{\mathcal{T}})$,
for $i=1,\ldots,N$.

Together with \eqref{Est.p}, we conclude \eqref{rhoi.cnv}. Moreover, \eqref{cnv.aux} holds as it is a straightforward consequence of \eqref{S.cnv}, \eqref{rhoi.cnv} and Vitali - De la Vall\'ee Poussin theorem.

{\em \textbf{Step 4: Strong convergence of the temperature.}} \\
Finally, we show \eqref{T.cnv} and \eqref{T.cnv.2}.
Let us consider the vector fields $W^{(n)}$, $Y^{(n)}$ defined in \eqref{Wn}, \eqref{Yn}. Thanks to \eqref{Est.entr}, \eqref{Est.entrflux}
we deduce that $W^{(n)}$ is bounded in $L^s(Q_\mathcal{T})$ for some $s>1$,
while $Y^{(n)}[G]$ is bounded in $L^\infty(Q_\mathcal{T})$ for every $G\in W^{1,\infty}(\R_+)$. On the other hand, the time-space divergence
$\dv_{(t,x)}W^{(n)}$ is the right-hand side of \eqref{EBE}, which is bounded in $L^1(Q_\mathcal{T})$ thanks to Lemma \ref{lem.ap.est}, while
$$
|\curl_{(t,x)}Y^{(n)}[G]|\leq C[G]|\nabla T|
$$
is bounded in 
$L^1(Q_\mathcal{T})$ due to \eqref{Est.A}. In order to see this, we distinguish two cases:\medskip\\
{\em Case $1<\beta\leq 2$.} In this case it holds 
$\nabla T^{(n)} = \frac{2}{\beta}(T^{(n)})^{1-\beta/2}\nabla (T^{(n)})^{\beta/2}$, which, given that $0\leq 1-\beta/2<1/2$, implies that $\nabla T^{(n)}$ is bounded in $L^1(Q_T)$ via the uniform bounds for $T^{(n)}$ and $\nabla (T^{(n)})^{\beta/2}$ in $L^\infty(0,T; L^1(\Omega))$ and $L^2(Q_T)$, respectively.\medskip\\
{\em Case $\beta>2$.} In this case we observe that
\begin{align*}
\nabla T^{(n)} = \chf{(0,1)}(T^{(n)}) T^{(n)} \nabla\logT^{(n)} + \chf{[1,\infty)}(T^{(n)})\frac{2}{\beta}(T^{(n)})^{1-\beta/2}\nabla (T^{(n)})^{\beta/2}
\end{align*}
is bounded in $L^2(Q_T)$ given the uniform bounds for $\nabla\logT^{(n)}$ and $\nabla (T^{(n)})^{\beta/2}$ in $L^2(Q_T)$.
In particular both 
$\dv_{(t,x)}W^{(n)}$, $\curl_{(t,x)}Y^{(n)}[G]$ are relatively compact in $W^{-1,r}(Q_\mathcal{T})$ for some $r>1$. From the Div-Curl Lemma it follows
\begin{align*}
\overline{W^{(n)}\cdot Y^{(n)}[G]} = 
\overline{W^{(n)}}\cdot \overline{Y^{(n)}[G]}
\end{align*}
which means
\begin{align*}
\overline{(\PhiS^{(n)}(\rho\eta)^{(n)} + (1-\Phi)(\rho\eta)_s^{(n)})G(T^{(n)})}
=
\overline{(\PhiS^{(n)}(\rho\eta)^{(n)} + (1-\Phi)(\rho\eta)_s^{(n)})}
\overline{G(T^{(n)})}\\
\mbox{a.e~in }Q_\mathcal{T},\quad\forall G\in W^{1,\infty}(\R_+).
\end{align*}
From the definitions \eqref{entropy}, \eqref{Skel.Entr},
bound \eqref{Est.A} and strong convergence relations \eqref{cnv.aux} 
we deduce
\begin{align*}
(c_w\Phi S \rho + c_s(1-\Phi))\overline{(\logT^{(n)})G(T^{(n)})} = 
(c_w\Phi S \rho + c_s(1-\Phi))\overline{(\logT^{(n)})}~
\overline{G(T^{(n)})},
\end{align*}
which, thanks to \eqref{Phi.pos}, yields
\begin{equation*}
\overline{(\logT^{(n)})G(T^{(n)})} = \overline{(\logT^{(n)})}~ \overline{G(T^{(n)})}\quad
\mbox{a.e.~in }Q_\mathcal{T},\quad
\forall G\in W^{1,\infty}(\R_+).
\end{equation*}
Since \eqref{Est.A} holds, by arguing in a similar way as the derivation of \eqref{id.rho} one obtains
\begin{equation}
\label{TGT}
\overline{(\logT^{(n)})T^{(n)}} = \overline{(\logT^{(n)})}~ \overline{T^{(n)}}
\quad\mbox{a.e.~in }Q_\mathcal{T} .
\end{equation}
Once again, the strict monotonicity and strict convexity of $\log$ allows us to conclude that $T^{(n)}$ is (up to subsequences) a.e.~convergent in $Q_\mathcal{T}$. From this fact, \eqref{Est.A}, \eqref{Est.T} and Sobolev's embedding $H^1(\Omega)\hookrightarrow L^6(\Omega)$ it follows
\eqref{T.cnv} and \eqref{T.cnv.2}.
\end{proof}
\subsection{Passing to the limit in the equations \eqref{weak.1}--\eqref{weak.EBE}, when $n \to \infty$.} In this subsection we will use the previous convergence results for passing to the limit in the variational entropy formulation.

\medskip

{\em Limit in equation \eqref{weak.1}.}
From \eqref{Est.rhov}, \eqref{Est.J}, 
\eqref{Est.Srhot}, \eqref{cnv.aux} it follows
\begin{align*}
\Phi\pa_t(S^{(n)}\rho^{(n)}_i)\rightharpoonup^*\Phi\pa_t(S\rho_i)\quad
\mbox{weakly* in }L^m(0,\mathcal{T}; W^{-1,m}(\Omega)),\\
\veb J_i^{(n)}\rightharpoonup \veb J_i\quad\mbox{weakly in }L^2(Q_\mathcal{T}),
\end{align*}
for $i=1,\ldots,N$.
Thanks to \eqref{vow}, \eqref{AJPP.1}, $(\lambda(S^{(n)},T^{(n)})T^{(n)})^{-1/2}\mathbf{v}^{(n)}$ is bounded in $L^2(Q_\mathcal{T})$;
however, eq.~\eqref{mobility}, the boundedness of $k_r\in C^0([0,1])$ and the uniform positivity of $\mu$ (see assumptions on phase velocity) imply that $\lambda(S^{(n)},T^{(n)}) = k_r(S^{(n)})/\mu(T^{(n)})$ is bounded in $L^\infty(Q_T)$, yielding that
$(T^{(n)})^{-1/2}\mathbf{v}^{(n)}$ is bounded in $L^2(Q_\mathcal{T})$. So, from \eqref{Est.T} we deduce that
$\mathbf{v}^{(n)}$ is bounded in $L^{\frac{6\beta + 4}{3\beta + 5}}(Q_\mathcal{T})$.
As a consequence,
\begin{align*}
\mathbf{v}^{(n)}\rightharpoonup \veb v\quad\mbox{weakly in }L^{\frac{6\beta + 4}{3\beta + 5}}(Q_\mathcal{T}).
\end{align*}
Since \eqref{rhoi.cnv} holds and (thanks to \eqref{hp.beta})
$\frac{1}{\gamma} + \frac{3\beta + 5}{6\beta + 4} < 1$, it follows
\begin{align*}
\rho^{(n)}_i\veb v^{(n)}\rightharpoonup \rho_i\veb v\quad\mbox{weakly in 
$L^s(Q_\mathcal{T})$ for some $s>1$.}
\end{align*}
We must now identify $\veb v$.
We start by pointing out that \eqref{lb.S} easily holds for $S^{(n)}$, too (just remake the same computations).
Therefore \eqref{vow}, \eqref{mobility} lead to
\begin{align*}
-K\nablap^{(n)}\rightharpoonup \frac{\veb v}{\lambda(S,T)}
\quad\mbox{weakly in }L^{\frac{6\beta + 4}{3\beta + 5}-\eps_1}(\Omega\times [\eps_2,\mathcal{T}]),\quad\forall\eps_1,\eps_2>0.
\end{align*}
From \eqref{p.cnv.Qeps} it follows
\begin{align*}
-K\lambda(S,T)\nabla p = \veb v\quad\mbox{a.e.~in }\Omega\times [\eps_2,\mathcal{T}],\quad\forall\eps_2>0.
\end{align*}
Being $\eps_2>0$ arbitrary and $\veb v\in L^{\frac{6\beta + 4}{3\beta + 5}}(Q_\mathcal{T})$, we deduce that \eqref{vow} holds a.e.~in $Q_\mathcal{T}$. As a consequence of this fact and  \eqref{AJPP.1}, we also deduce
\begin{equation}
\label{wcnv.nap.1}
\frac{\lambda(S^{(n)},T^{(n)})^{1/2}}{(T^{(n)})^{1/2}}\nabla p^{(n)}
\rightharpoonup
\frac{\lambda(S,T)^{1/2}}{T^{1/2}}\nabla p\quad\mbox{weakly in }L^2(Q_\mathcal{T}).
\end{equation} 
We will now identify the limits of the diffusion fluxes $\veb J_1,\ldots,\veb J_N$.
From \eqref{Est.B.2} it follows
\begin{align}\label{muT.cnv}
\Pi^N\vec\mu^{(n)}/T^{(n)} \rightharpoonup\vec\zeta
\quad\mbox{weakly in }L^2(0,\mathcal{T}; H^1(\Omega)),
\end{align}
which, thanks to \eqref{hp.Ltilde}, \eqref{Ass.L}, yields, for $i=1,\ldots,N$:
\begin{align}\label{lim.muT}
\sum_{j=1}^N \tilde L_{ij}(\vec\rho^{(n)},T^{(n)})\nabla(\mu_j^{(n)}/T^{(n)}) \rightharpoonup
\sum_{j=1}^N \tilde L_{ij}(\vec\rho,T)\nabla\zeta_j
\quad\mbox{weakly in }L^2(Q_\mathcal{T}).
\end{align}
In order to identify $\vec\zeta$, we point out that 
\begin{equation}
\label{rho.0}
\mbox{For a.e.~}(x,t)\in Q_\mathcal{T},\quad
\rho(x,t)>0\quad\Rightarrow\quad\min_{i=1,\ldots,N}\rho_i(x,t)>0.
\end{equation}
This is a consequence of \eqref{est.relchp}. In fact, since $\rho^{(n)}_i$ is a.e.~convergent in $Q_\mathcal{T}$, Fatou's Lemma implies
\begin{align*}
\int_{Q_\mathcal{T}\cap\{\rho>0\}}\log(\rho/\rho_i)dx\leq
\liminf_{n\to\infty}
\int_{Q_\mathcal{T}\cap\{\rho>0\}}\log(\rho^{(n)}/\rho^{(n)}_i)dx\leq C,
\end{align*}
for $i=1,\ldots,N$. This means that $\rho_i>0$ on $Q_\mathcal{T}\cap\{\rho>0\}$ for $i=1,\ldots,N$.

Thanks to the a.e.~convergence of $\rho^{(n)}_i$ and \eqref{rho.0}, we deduce
\begin{align}\label{lim.muT.2}
\Pi^N\frac{\vec\mu^{(n)}}{T^{(n)}}\to 
\vec\zeta = \Pi^N\frac{\vec\mu}{T}\quad\mbox{a.e.~on }Q_\mathcal{T}\cap\{\rho>0\} ,\\
\nonumber
\left(\Pi^N\frac{\vec\mu}{T}\right)_i = \log\left( \frac{\rho_i}{\rho} \right)\quad\mbox{a.e.~in }Q_\mathcal{T}\cap\{\rho>0\} ,\quad i=1,\ldots,N.
\end{align}
From \eqref{hp.Ltilde}, \eqref{Est.A}, \eqref{rhoi.cnv}, \eqref{T.cnv} we deduce that,
for $i=1,\ldots,N$,
\begin{align}
\label{Li0T.cnv}
\frac{\tilde L_{i0}(\vec{\rho}^{(n)},T^{(n)})}{T^{(n)}}\to
\frac{\tilde L_{i0}(\vec\rho,T)}{T}\quad\mbox{strongly in }
L^p(Q_\mathcal{T}),\quad \forall p<\infty,\\
\nonumber
\nabla\logT^{(n)}\rightharpoonup\nabla\log T\quad\mbox{weakly in }L^2(Q_\mathcal{T}).
\end{align}
Since $\frac{\tilde L_{i0}(\vec{\rho}^{(n)},T^{(n)})}{T^{(n)}}\nabla\logT^{(n)}$ is bounded in $L^2(Q_\mathcal{T})$, we deduce (up to subsequences)
\begin{align*}
\frac{\tilde L_{i0}(\vec{\rho}^{(n)},T^{(n)})}{T^{(n)}}\nabla\logT^{(n)}\rightharpoonup
\frac{\tilde L_{i0}(\vec\rho,T)}{T}\nabla\log T\quad\mbox{weakly in }
L^2(Q_\mathcal{T}),
\end{align*}
for $i=1,\ldots,N$, which means
\begin{align}\label{lim.cross.1}
\tilde L_{i0}(\vec{\rho}^{(n)},T^{(n)})\nabla\frac{1}{T^{(n)}}\rightharpoonup
\tilde L_{i0}(\vec\rho,T)\nabla\frac{1}{T}\quad\mbox{weakly in }
L^2(Q_\mathcal{T}).
\end{align}
From \eqref{lim.muT}, \eqref{lim.muT.2}, \eqref{lim.cross.1} we obtain
\begin{align}
\label{lim.J}
\veb{J}_i &= L_{i0} \nabla \Big(\frac{1}{T}\Big) 
- \sum_{k=1}^{N} L_{ik} \nabla\zeta_k, \quad \mbox{a.e.~on }Q_\mathcal{T} ,\\
L_{ij} &= \tilde L_{ij}(\vec\rho,T),\quad
L_{i0} = \tilde L_{i0}(\vec\rho,T)\quad
\mbox{a.e.~on }Q_\mathcal{T} .\nonumber
\end{align}
for $i,j=1,\ldots,N$. Let us now focus on the reaction terms.
From \eqref{Est.r} it follows that
\begin{align*}
r_i^{(n)}\rightharpoonup r_i\quad\mbox{weakly in }L^{a/(a-1)}(Q_\mathcal{T}),
\quad i=1,\ldots,N.
\end{align*}
However, since $r_i^{(n)} = \tilde r_i(\rho^{(n)}, T^{(n)}, \Pi^N\vec\mu^{(n)}/T^{(n)})$, the continuity of $\tilde r_i$ as well as
\eqref{rho.cnv}, \eqref{T.cnv}, \eqref{lim.muT.2} we deduce
\begin{equation}
\label{lim.r}
r_i = \tilde r_i(\rho, T, \Pi^N\vec\mu/T)\quad\mbox{a.e.~in }Q_\mathcal{T}\cap\{\rho>0\}.
\end{equation}
Finally, we point out that the continuous Sobolev embedding $H^{1}(\Omega)\hookrightarrow H^{1/2}(\pa\Omega)$ yields the convergence of the boundary integral:
\begin{align*}
\sum_{k=1}^N b_{ik}\left(\frac{\mu^{(n)}_k}{T^{(n)}} - \frac{\mu_{0,k}}{T_0}  \right)\rightharpoonup
\sum_{k=1}^N b_{ik}\left(\zeta_k - \frac{\mu_{0,k}}{T_0}\right)
\quad\mbox{weakly in }L^2(\pa\Omega\times (0,\mathcal{T})).
\end{align*}
This shows that \eqref{weak.1} holds.

\medskip

{\em Limit in equation \eqref{weak.EBE}.}
Let us now turn our attention to $\veb q^{(n)}$. 
From \eqref{Est.qT} it follows 
\begin{align}\label{q.cnv.1}
\frac{\veb q^{(n)}}{T^{(n)}}\rightharpoonup\tilde{\veb q}\quad
\mbox{weakly in }L^{\frac{2+3\beta}{1+3\beta}}(Q_\mathcal{T}).
\end{align}
On the other hand, from \eqref{flux.q} it follows
\begin{align*}
\frac{\veb q^{(n)}}{T^{(n)}} = 
-\frac{\kappa(T^{(n)})}{T^{(n)}}\nablaT^{(n)} 
+ \sum_{j=1}^N\frac{\tilde L_{0j}(\vec{\rho}^{(n)},T^{(n)})}{T^{(n)}}\nabla\frac{\mu_j^{(n)}}{T^{(n)}} .
\end{align*}
From \eqref{hp.Ltilde}, \eqref{muT.cnv}, \eqref{Li0T.cnv} and the boundedness of $\sum_{j=1}^N\frac{\tilde L_{0j}(\vec{\rho}^{(n)},T^{(n)})}{T^{(n)}}\nabla\frac{\mu_j^{(n)}}{T^{(n)}}$
in $L^2(Q_\mathcal{T})$ we deduce
\begin{align*}
\sum_{j=1}^N
\frac{\tilde L_{0j}(\vec{\rho}^{(n)},T^{(n)})}{T^{(n)}}\nabla\frac{\mu_j^{(n)}}{T^{(n)}}
\rightharpoonup
\sum_{j=1}^N\frac{\tilde L_{0j}(\vec\rho,T)}{T}\nabla\zeta_j
\quad\mbox{weakly in }L^2(Q_\mathcal{T}).
\end{align*}
On the other hand, \eqref{Est.kappaT}, \eqref{T.cnv} imply
\begin{align*}
\frac{\kappa(T^{(n)})}{T^{(n)}}\nablaT^{(n)} \rightharpoonup
\frac{\kappa(T)}{T}\nabla T \quad\mbox{weakly in }L^{\frac{2+3\beta}{1+3\beta}}(Q_\mathcal{T}) .
\end{align*}
We conclude 
\begin{equation}
\label{qtilde}
\tilde{\veb q} = 
-\frac{\kappa(T)}{T}\nabla T 
+ \sum_{j=1}^N\frac{\tilde L_{0j}(\vec\rho,T)}{T}\nabla\zeta_j 
=\frac{\veb q}{T},\quad
\mbox{a.e.~in }Q_\mathcal{T}.
\end{equation}
From \eqref{Est.muJ}, \eqref{lim.muT.2} it follows
\begin{align*}
-\sum_{i=1}^N\frac{\mu_i^{(n)}}{T^{(n)}}\veb J_i^{(n)}\rightharpoonup
\Xi\quad\mbox{weakly in }L^{\frac{2a}{a+2}}(Q_\mathcal{T}),\quad
\Xi = -\sum_{i=1}^N\frac{\mu_i}{T}\veb J_i\quad\mbox{a.e.~in }Q_\mathcal{T}\cap\{\rho>0\}.
\end{align*}
Finally, consider the quantity
\begin{align*}
\xi^{(n)}
&=\frac{K}{T^{(n)}}\lambda(S^{(n)},T^{(n)})|\nabla p^{(n)}|^2 +
\sum_{i,j=1}^N L_{ij}^{(n)} \nabla \big( \frac{\mu_i^{(n)}}{T^{(n)}} \big) \cdot \nabla \big( \frac{\mu_j^{(n)}}{T^{(n)}} \big) + L_{00}^{(n)} \big| \nabla \big( \frac{1}{T^{(n)}} \big) \big|^2\nonumber\\ 
&\qquad - \Phi \frac{1}{T^{(n)}} f'(S^{(n)}) \big( \partial_t S^{(n)} \big)^2
-\sum_{i=1}^N r_i^{(n)}\frac{\mu_i^{(n)}}{T^{(n)}}\\
&-\frac{K}{T}\lambda(S,T)|\nabla p|^2 
-\sum_{i,j=1}^N L_{ij} \nabla \big( \frac{\mu_i}{T} \big) \cdot \nabla \big( \frac{\mu_j}{T} \big) - L_{00} \big| \nabla \big( \frac{1}{T} \big) \big|^2\nonumber\\ 
&\qquad +\Phi \frac{1}{T} f'(S) \big( \partial_t S \big)^2
+\sum_{i=1}^N r_i\frac{\mu_i}{T}.
\end{align*}
It is clear that $\xi^{(n)}$ is bounded in $L^1(Q_\mathcal{T})$ and therefore
also in $\mathcal M(\overline{Q_\mathcal{T}})$. Therefore 
$\xi^{(n)}\rightharpoonup^*\xi$ weakly* in $\mathcal M(\overline{Q_\mathcal{T}})$. 


We will prove that $\xi$ is a nonnegative measure.
Let $\varphi\in C^0(\overline{Q_\mathcal{T}})$, $\varphi\geq 0$
in $Q_\mathcal{T}$.	
From \eqref{wcnv.nap.1} and the weak lower semicontinuity of 
the $L^2(Q_\mathcal{T})$ norm we deduce
\begin{align}\label{Fatou.1}
\liminf_{n\to\infty}\int_{Q_\mathcal{T}}\frac{K}{T^{(n)}}\lambda(S^{(n)},T^{(n)})|\nabla p^{(n)}|^2 \varphi dx dt \geq
\int_{Q_\mathcal{T}}\frac{K}{T}\lambda(S,T)|\nabla p|^2 
\varphi dx dt.
\end{align}
Furthermore, from \eqref{Est.S} it follows that
\begin{align*}
\frac{1}{(T^{(n)})^{1/2}} \sqrt{-f'(S^{(n)})} \partial_t S^{(n)}
\rightharpoonup
\omega\quad\mbox{weakly in }L^2(Q_\mathcal{T}).
\end{align*}
Since \eqref{T.cnv.2} holds we deduce
\begin{align*}
\sqrt{-f'(S^{(n)})} \partial_t S^{(n)}
\rightharpoonup
\sqrt{T}\omega\quad\mbox{weakly in }L^1(Q_\mathcal{T}).
\end{align*}
Testing $\sqrt{-f'(S^{(n)})} \partial_t S^{(n)}$ against an arbitrary test function, integrating by parts in time and 
taking the limit $n\to\infty$ allows us to easily determine the value of $\omega$, yielding 
\begin{align*}
\frac{1}{(T^{(n)})^{1/2}} \sqrt{-f'(S^{(n)})} \partial_t S^{(n)}
\rightharpoonup
\frac{1}{T^{1/2}} \sqrt{-f'(S)} \partial_t S
\quad\mbox{weakly in }L^2(Q_\mathcal{T}).
\end{align*}
Once again, the above relation and the weak lower semicontinuity of the $L^2(Q_\mathcal{T})$ norm lead to
\begin{align}\label{Fatou.2}
\liminf_{n\to\infty}\int_{Q_\mathcal{T}}
- \Phi \frac{1}{T^{(n)}} f'(S^{(n)}) \big( \partial_t S^{(n)} \big)^2 
\varphi dx dt \geq
\int_{Q_\mathcal{T}} 
-\Phi \frac{1}{T} f'(S) \big( \partial_t S \big)^2 
\varphi dx dt .
\end{align}
Finally, Lemma \ref{lem.Fatou.gen} and assumptions \eqref{hp.r.3}, \eqref{hp.r.4} imply that
\begin{align}
\label{Fatou.3}
\liminf_{n\to\infty}\int_{Q_\mathcal{T}}
\sum_{i,j=1}^N L_{ij}^{(n)} \nabla \big( \frac{\mu_i^{(n)}}{T^{(n)}} \big) \cdot \nabla \big( \frac{\mu_j^{(n)}}{T^{(n)}} \big) \varphi dx dt &\geq 
\int_{Q_\mathcal{T}}
\sum_{i,j=1}^N L_{ij} \nabla \big( \frac{\mu_i}{T} \big) \cdot \nabla \big( \frac{\mu_j}{T} \big) \varphi dx dt,\\
\label{Fatou.4}
\liminf_{n\to\infty}\int_{Q_\mathcal{T}}
L_{00}^{(n)} \big| \nabla \big( \frac{1}{T^{(n)}} \big) \big|^2
\varphi dx dt &\geq \int_{Q_\mathcal{T}}
L_{00} \big| \nabla \big( \frac{1}{T} \big) \big|^2
\varphi dx dt,\\
\label{Fatou.5}
\liminf_{n\to\infty}\int_{Q_\mathcal{T}}
-\sum_{i=1}^N r_i^{(n)}\frac{\mu_i^{(n)}}{T^{(n)}} \varphi dx dt
&\geq\int_{Q_\mathcal{T}}
-\sum_{i=1}^N r_i\frac{\mu_i}{T} \varphi dx dt.
\end{align}
From \eqref{Fatou.1}--\eqref{Fatou.5} and the definition of $\xi$ we conclude
\begin{align*}
\langle\xi,\varphi\rangle = 
 \lim_{n\to\infty}\langle\xi^{(n)},\varphi\rangle\geq 0.
\end{align*}
This means that $\xi$ is nonnegative. Therefore \eqref{weak.EBE} holds.

\medskip

{\em Limit in equation \eqref{weak.2}.}
We will now prove that \eqref{weak.2} holds in the limit $n\to\infty$. From \eqref{Est.A}, \eqref{T.cnv} it follows that
$T^{(n)}\rightharpoonup T$ weakly in $L^2(0,\mathcal{T}; H^1(\Omega))$.
Thanks to the Sobolev embedding $H^1(\Omega)\hookrightarrow H^{1/2}(\pa\Omega)$ this implies that
$T^{(n)}\rightharpoonup T$ weakly in $L^2(0,\mathcal{T}; H^{1/2}(\pa\Omega))$. In particular
\begin{align*}
\int_0^t\int_{\pa\Omega}(T^{(n)} - T_0)d\sigma dt'\to
\int_0^t\int_{\pa\Omega}(T - T_0)d\sigma dt',\quad n\to\infty .
\end{align*}
Furthermore $E_s(T^{(n)}) = c_sT^{(n)}\to c_s T = E_s(T)$ strongly in $L^2(Q_\mathcal{T})$ thanks to 
\eqref{Skel.Ener}, \eqref{T.cnv.2}, while the term (recall \eqref{def.Eint})
\begin{align*}
E_{int}(S^{(n)}) = \int_{S^{(n)}}^1 P_c(s) ds	
\end{align*}
can be estimated from \eqref{hp.Pc.bound} and \eqref{Est.p} 
through a similar argument as in \eqref{Est.SPc}, implying that
\begin{align*}
E_{int}(S^{(n)})\to E_{int}(S)\quad\mbox{strongly in }L^1(Q_\mathcal{T}).
\end{align*}
Therefore we must now only show that the term (recall \eqref{eng.1})
\begin{align*}
S^{(n)} (\rho e)^{(n)} = S^{(n)} (\rho^{(n)})^\gamma + c_w\rho^{(n)} T^{(n)} + p_{at}
\end{align*}
is strongly convergent in $L^1(Q_\mathcal{T})$. However, from \eqref{thp.1} we obtain
\begin{align*}
(\gamma-1)S^{(n)} (\rho e)^{(n)} - S^{(n)} p^{(n)} = 
[c_w(\gamma-1)-1]S^{(n)}\rho^{(n)}T^{(n)} + \gamma p_{at} ,
\end{align*}
which, thanks to \eqref{S.cnv}, \eqref{rho.cnv}, \eqref{T.cnv}, \eqref{p.cnv}, leads to
\begin{align*}
(\gamma-1)S^{(n)} (\rho e)^{(n)} \to S p + 
[c_w(\gamma-1)-1]S\rho T + \gamma p_{at} = 
(\gamma-1)S (\rho e)\\
\mbox{strongly in }L^1(Q_\mathcal{T}).
\end{align*}
We conclude that \eqref{weak.2} holds.
\medskip

 {\em Limit in equation \eqref{PC.Sw}. We will now show that the saturation balance equation \eqref{PC.Sw} holds in the limit $n\to\infty$. }

We have by assumption
\begin{align*}
	\pa_t f(S^{(n)}) + P_c(S^{(n)}) + p^{(n)} = 0,\quad t>0,\quad\mbox{a.e.~in }\Omega .
\end{align*}
Multiplying the equation times $S^{(n)}$ leads to
\begin{align*}
	\pa_t F(S^{(n)}) + S^{(n)} P_c(S^{(n)}) + S^{(n)} p^{(n)} = 0,\quad t>0,\quad\mbox{a.e.~in }\Omega ,
\end{align*}
where $F$ is like in \eqref{Est.FSt}.
Integrating the above equation against a test function $\varphi\in C^1_c(Q_\mathcal{T})$ yields
\begin{align}\label{eq.F}
	-\int_0^\mathcal{T}\int_\Omega F(S^{(n)})\pa_t\varphi dx dt 
	+ \int_0^\mathcal{T}\int_\Omega \left(S^{(n)} P_c(S^{(n)}) + S^{(n)} p^{(n)} \right)\varphi dx dt = 0.
\end{align}
Since $F(S) = -\int_S^1 s_1 f'(s_1)d s_1 = -f(1) + Sf(S) + 
\int_S^1 f(s_1)d s_1$ and $f\in C^0([0,1))$, then also $F\in C^0([0,1))$. Therefore from \eqref{S.cnv} we deduce
$F(S^{(n)})\to F(S)$ strongly in $L^r(Q_\mathcal{T})$ for every $r<\infty$.
On the other hand, \eqref{hp.Pc.bound} and \eqref{Est.p} imply
\begin{align}\label{Est.SPc}
	\begin{cases}
		\|S^{(n)} P_c(S^{(n)})\|_{L^{\infty}(Q_\mathcal{T})} \leq
		C & k_p\leq 1,\\
		\|S^{(n)} P_c(S^{(n)})\|_{L^{\frac{k_p}{k_p-1}}(Q_\mathcal{T})} \leq
		C\| P_c(S^{(n)})^{1-1/k_p} \|_{L^{\frac{k_p}{k_p-1}}(Q_\mathcal{T})} 
		\leq C & k_p>1
	\end{cases} .
\end{align}
Since $z_p>1$ and \eqref{S.cnv} holds, we deduce that 
$S^{(n)} P_c(S^{(n)})\to S P_c(S)$ strongly in $L^1(Q_\mathcal{T})$.
Using \eqref{S.cnv} and \eqref{p.cnv} we can take the limit $n\to\infty$ in \eqref{eq.F} and get
\begin{align*}
	-\int_0^\mathcal{T}\int_\Omega F(S)\pa_t\varphi dx dt 
	+ \int_0^\mathcal{T}\int_\Omega \left(S P_c(S) + S p \right)\varphi dx dt = 0,\quad\forall\varphi\in C^1_c(Q_\mathcal{T}),
\end{align*}
which is the weak formulation of
\begin{align*}
	\pa_t F(S) + S P_c(S) + S p = 0,\quad t>0,\quad\mbox{a.e.~in }\Omega .
\end{align*}
Since $F'(S) = S f'(S)$ and $S>0$ a.e.~in $Q_\mathcal{T}$, dividing the above equation times $S$ yields \eqref{PC.Sw}.






This finishes the proof of Theorem \ref{thm.ws}.







\section{The approximate scheme}\label{sec.appr}
In this section we build a sequence of approximated solutions to \eqref{cons_w}--\eqref{PC.Sw}, \eqref{ic.1},
\eqref{bc}.

Let $K_1>0$, $K_2>0$, $K_3>0$ be generic but fixed positive constants satisfying
\begin{align}
\label{hp.K123}
\begin{cases}
K_1 > \frac{6}{5}\gamma & \alpha_r \leq \frac{4}{3}\\
\frac{6}{5}\gamma < K_1 < \frac{3\alpha_r - 2}{3\alpha_r - 4}\gamma  & \alpha_r > \frac{4}{3}
\end{cases},\qquad
3 < K_2 < 3\beta ,\qquad
K_3 > \frac{5 K_1 + 6\gamma}{5 K_1 - 6\gamma}.
\end{align}
We point out that, for $\alpha_r>4/3$, it holds $\frac{3\alpha_r - 2}{3\alpha_r - 4}>\frac{6}{5}$ thanks to \eqref{hp.kr}. 
\subsection{Discretization and regularization}
We introduce the following parameters of the approximate scheme:
\begin{itemize}
\item[(i)]  $\eps> 0$ lower-order regularization of \eqref{cons_w}-\eqref{cons_e},
as well as the regularization in the free energy, the source term, and the capillary pressure,
 \item[(ii)]  $\delta > 0$ higher-order regularization in \eqref{cons_w}-\eqref{cons_e},
 \item[(iii)] $\tau > 0$ time-step in the implicit Euler discretization.
\end{itemize}

{\bf Regularization.}
We introduce the "regularized" free energy for the approximate system:
\begin{equation}
\label{free.eps}
(\rho\Psi)_{w,\eps} = (\rho\Psi)_w + \eps\rho^{K_1},\quad
(\rho\Psi)_{s,\eps} = (\rho\Psi)_{s} - \eps T^{K_2}.
\end{equation}
The definitions of the thermodynamics quantities $\mu_i$, $p$, $\rho e$, $(\rho\eta)_s$, $E_s$ change accordingly:
\begin{equation}
\label{mu.eps}
\mu_{i,\eps} = \frac{\pa (\rho\Psi)_{w,\eps}}{\pa\rho_i} =
\frac{\pa (\rho\Psi)_{w}}{\pa\rho_i} + \eps K_1\rho^{K_1-1},
\end{equation}
\begin{equation}
\label{p.eps}
p_\eps = -(\rho\Psi)_{w,\eps} + \sum_{i=1}^N\rho_i\mu_{i,\eps} = 
-(\rho\Psi)_{w} + \sum_{i=1}^N\rho_i\mu_{i} 
+\eps (K_1 - 1)\rho^{K_1},
\end{equation}
\begin{equation}
\label{rhoe.eps}
(\rho e)_{\eps} = (\rho\Psi)_{w,\eps} - T\frac{\pa(\rho\Psi)_{w,\eps}}{\pa T} =
(\rho\Psi)_{w} - T\frac{\pa(\rho\Psi)_{w}}{\pa T} 
+\eps\rho^{K_1},
\end{equation}
\begin{equation}
\label{rhoetas.eps}
(\rho\eta)_{s,\eps} = -\frac{\pa(\rho\Psi)_s}{\pa T} 
+\eps K_2 T^{K_2-1},
\end{equation}
\begin{equation}
\label{Es.eps}
E_{s,\eps} = (\rho\Psi)_{s,\eps} - T\frac{\pa(\rho\Psi)_{s,\eps}}{\pa T}=
(\rho\Psi)_{s} - T\frac{\pa(\rho\Psi)_{s}}{\pa T}
+ \eps(K_2-1)T^{K_2}.
\end{equation}
However, in what follows, we will denote the regularized quantities $\mu_{i,\eps}$, $p_\eps$, $(\rho e)_\eps$,
$(\rho\eta)_{s,\eps}$, $E_{s,\eps}$ with
$\mu_i$, $p$, $\rho e$, $(\rho\eta)_{s}$, $E_{s}$
to avoid a too cumbersome notation.
Also, notice that the phase entropy $\rho\eta$ remains unchanged. We also point out that the regularized skeleton entropy $(\rho\eta)_s$ is still a concave function of the regularized skeleton energy $E_s$ (easy computations).

We regularize the source term $r_i$:
\begin{align}
	\label{ri.eps}
r_{i,\eps} = r_i - \eps\left|\frac{\mu_i}{T}\right|^{a-2}\frac{\mu_i}{T},\quad
i=1,\ldots,N.
\end{align}
We also regularize the dynamic capillary pressure term $P_c$ in the following way:
\begin{align}\label{def.Pceps}
& P_{c,\eps}(s) =
\begin{cases}
 P_c(s) & \mbox{if }k_p>0\\
 P_c(s) - \eps\log s & \mbox{if }k_p=0
\end{cases},\quad
E_{f,\eps}^k = (\rho e)_{\eps}^k S^k 
- \int_{1/2}^{S^k} P_{c,\eps}(\xi)d\xi,
\end{align}
where the parameter $k_p\geq 0$ is defined in \eqref{hp.Pc.bound}.

Finally, we truncate the initial value for the saturation:
\begin{align}\label{Sin.eps}
S^{in,\eps} = \max\{\eps, \min\{S^{in},1-\eps\}\} .
\end{align}
{\bf Discretization.} Fix $\mathcal{T} > 0$. For $L \in \N$ we define $\tau = \mathcal{T}/L$, $t_k = \tau k$ $(k = 0, \ldots,L)$.
Consider the implicit Euler time discretization in \eqref{cons_w}--\eqref{PC.Sw} with regularizing terms:
\begin{align}
&\int_\Omega\tau^{-1}(S^k\rho_i^k - S^{k-1}\rho_i^{k-1})
\Phi\varphi dx
-\int_\Omega( \rho_i^k\veb v^k + \veb J_i^k )\cdot\nabla\varphi dx 
 -\int_\Omega r_{i,\eps}^k\varphi dx \nonumber \\
 & +\int_{\pa\Omega}\sum_{\ell=1}^N b_{i\ell}\left( \frac{\mu_\ell^k}{T^k}
- \frac{\mu_{0,\ell}}{T_0} \right)\varphi ds = - \eps R_{\eps,1}(\mu_i^k/T^k, \varphi) - \delta \tilde{R}_{\delta,1}(\mu_i^k/T^k,\varphi),  \label{app.1} 
\end{align}
\begin{align}
&\int_\Omega\tau^{-1}( \Phi(E_{f,\eps}^k - E_{f,\eps}^{k-1}) + 
(1-\Phi)(E_s^k - E_s^{k-1}) )\psi dx \nonumber \\
& -\int_\Omega( ((\rho e)^k + p^k)\veb v^k + \veb q^k )\cdot\nabla\psi dx 
 +\alpha\int_{\pa\Omega}(T^k - T_0)\varphi ds \nonumber \\
& =-\eps \Big(  R_{\eps,2}(T^k,\log T^k , \psi)  +  R_{\eps,3}(T^k,\log T^k,\psi) \Big)  \nonumber\\
& - \delta \Big( \tilde{R}_{\delta,2}(T^k, \log T^k, \psi)  +  \tilde{R}_{\delta,3}(T^k, \log T^k, \psi)   \Big),  \label{app.2}
 \end{align}
 \begin{align}
 \frac{f(S^k)-f(S^{k-1})}{\tau} +
P_{c,\eps}(S^k) + p^k= 0, \label{app.3}
 \end{align}
for every $\varphi, \psi \in C^2(\overline\Omega)$, where we defined:
\begin{align}
 R_{\eps,1}(\mu_i^k/T^k, \varphi) & = \int_\Omega\left(
\nabla\frac{\mu_i^k}{T^k}\cdot\nabla\varphi + 
\frac{\mu_i^k}{T^k}\varphi\right)dx, \label{Reps1} \\
 R_{\eps,2}(T^k,\log T^k , \psi)  & = \int_\Omega((1+T^k)\nabla\log T^k\cdot\nabla\varphi + (1+(T^k)^{-K_3})(\log T^k)\psi ) dx, \label{Reps2}\\
 R_{\eps,3}(T^k,\log T^k,\psi) & = \int_\Omega (T^k)^{-K_3} |\nabla\log T^k|^{K_3-1}\nabla\log T^k\cdot\nabla\psi dx, \label{Reps3}\\
 \tilde{R}_{\delta,1}(\mu_i^k/T^k,\varphi)  & = \int_\Omega D^2\frac{\mu_i^k}{T^k} : D^2\varphi dx, \label{Rdelta1}\\
 \tilde{R}_{\delta,2}(T^k, \log T^k, \psi)  & = \int_\Omega (1+T^k) D^2\log T^k : D^2\psi dx, \label{Rdelta2}\\
 \tilde{R}_{\delta,3}(T^k, \log T^k, \psi) & =  \int_\Omega (1 + T^k)|\nabla\log T^k|^2\nabla\log T^k\cdot\nabla\psi) dx. \label{Rdelta3}
\end{align}

We will first prove the well-posedness of the approximate system and then perform the limits $\delta\to 0$, $\tau\to 0$, $\eps\to 0$ {\em in this order.}


\subsection{Existence of solutions to the approximate system.} 
We formulate \eqref{app.1}-\eqref{app.2}-\eqref{app.3} as a fixed-point problem for an operator $F : X \times [0,1] \to X$ defined by means of a linearized problem.

\subsubsection{Reformulation as a fixed-point problem.}
Define the spaces
\[   V =  W^{1,4}(\Omega; \R^N) \times  W^{1,4}(\Omega),\quad
   V_0 = H^2(\Omega; \R^N) \times  H^2(\Omega),\]
and variables 
$$\vec z = \frac{\vec\mu}{T},\quad w = \log T . $$ 
For arbitrary given $(\vec z^*, w^*)\in V \hookrightarrow L^\infty(\Omega; \R^{N}) \times L^\infty(\Omega)$ (remember $\Omega\subset\R^3$), define 
$T^* = \exp(w^*)$, $\vec\mu^* = T^*\vec z^*$, and the quantities $\vec\rho^*$, $\veb v^*$ etc.~accordingly to the constitutive relations.

For $(\vec z^*, w^*)\in V$ and $\sigma\in [0,1]$, let us consider the linearized problem
\begin{align}
\label{lin.1}
&\sigma\left\{ \int_\Omega\tau^{-1}(S\rho_i^* - S^{k-1}\rho_i^{k-1})
\Phi\varphi dx
-\int_\Omega( \rho_i^*\veb v^* + \veb J_i^* )\cdot\nabla\varphi dx\right.\\ \nonumber
&\qquad\left.
-\int_\Omega r_{i,\eps}^*\varphi dx
+\int_{\pa\Omega}\sum_{\ell=1}^N b_{i\ell}\left( \frac{\mu_\ell^*}{T^*}
- \frac{\mu_{0,\ell}}{T_0} \right)\varphi ds
\right\}\\ \nonumber
&\quad = -\eps   R_{\eps,1}(z_i,\varphi) - \delta R_{\delta,1}(z_i, \varphi),
\quad\varphi\in H^2(\Omega),\\
\label{lin.2}
&\sigma\left\{
\int_\Omega\tau^{-1}( \Phi(E_{f,\eps}^* - E_{f,\eps}^{k-1}) + 
(1-\Phi)(E_s^* - E_s^{k-1}) )\psi dx\right. \\ \nonumber
&\left.
-\int_\Omega( ((\rho e)^* + p^*)\veb v^* + \veb q^* )\cdot\nabla\psi dx
+\alpha\int_{\pa\Omega}(T^* - T_0)\psi ds\right\}\\ \nonumber
& = -\eps \Big( R_{\eps,2}(T^*, w, \psi) + R_{\eps,3}(T^*,w, \psi)  \Big) - 
  \delta \Big(  R_{\delta,2}(T^*, w, \psi) + R_{\delta,3}(T^*,w, \psi)   \Big), \\
  \nonumber
&\qquad\psi\in H^2(\Omega),
\end{align}
coupled with the nonlinear algebraic equation
\begin{equation}
\label{lin.3}
\tau^{-1}(f(S)-f(S^{k-1})) + \sigma\left( P_{c,\eps}(S) + p^* \right) = 0.
\end{equation}
The strict monotonicity of $f$, $P_{c,\eps}$ as well as the fact that 
\begin{align*}
\lim_{s\to 0^+} P_{c,\eps}(s) = \infty,\quad
\lim_{s\to 1^-} f(s) = -\infty
\end{align*}
imply that \eqref{lin.3} has a unique solution $S : \Omega\to (0,1)$ which is Lebesgue-measurable. Moreover, 
the Lax-Milgram Lemma allows us to state that system \eqref{lin.1}, \eqref{lin.2} has a unique solution $(\vec z, w)\in V_0\hookrightarrow V$. We consider therefore the mapping
\begin{align*}
F: ((\vec z^*, w^*),\sigma)\in V\times [0,1]\mapsto (\vec z, w)\in V.
\end{align*}
It holds that $F(\cdot,0) \equiv 0$. Furthermore $F$ is continuous (standard argument) and bounded as an operator 
$V \times [0,1]\to V_0$ (just choose $\varphi = z_i$, $\psi=w$ in \eqref{lin.1}, \eqref{lin.2}, respectively, and estimate the left-hand sides of the resulting equations in a straightforward way), which means (thanks to the compact Sobolev embedding $V_0\hookrightarrow V$) that $F$ is compact.
Now we wish to prove that the set
\begin{align}\label{set.fp}
\mathcal F_\sigma\equiv
\left\{ (\vec z, w)\in V~ : ~
F((\vec z,w),\sigma) = (\vec z,w) \right\}
\end{align}
is bounded in $V$ uniformly w.r.t.~$\sigma\in [0,1]$. This will allow us to apply the Leray-Schauder fixed point theorem and deduce the existence of a fix point for $F(\cdot,1)$, that is, a solution to the approximate system \eqref{app.1}--\eqref{app.3}.

\subsubsection{Global entropy and energy balance for approximate solutions.}
Let $(\vec z, w)\in\mathcal F_\sigma$. Recall that 
$T = \exp(w)$, $\vec\mu = T\vec z$. Choose $\varphi = z_i = \mu_i/T$, $\psi = -1/T$ in \eqref{lin.1}, \eqref{lin.2}, respectively (also keep in mind that $(\vec z, w) = (\vec z^*, w^*)$, so we omit the $*$ sign everywhere). It follows (the index $i$ is summed from $1$ to $N$)
\begin{align}\label{app.ei.1}
\sigma &\left\{ \tau^{-1}\int_\Omega\left(  
\Phi (S\rho_i - S^{k-1}\rho_i^{k-1})\frac{\mu_i}{T} 
-( \Phi(E_{f,\eps} - E_{f,\eps}^{k-1}) + (1-\Phi)(E_s - E_s^{k-1}) )\frac{1}{T}
\right) dx\right. \\ \nonumber
&\left. -\int_\Omega( \rho_i\veb v + \veb J_i )\cdot\nabla\frac{\mu_i}{T} dx
+\int_\Omega( ((\rho e) + p)\veb v + \veb q )\cdot
\nabla\frac{1}{T} dx
+\alpha\int_{\pa\Omega}\frac{T_0 - T}{T} ds
\right.\\ \nonumber
&\left.
-\int_\Omega r_{i,\eps}\frac{\mu_i}{T} dx
+\int_{\pa\Omega}\sum_{\ell=1}^N b_{i\ell}\left( \frac{\mu_\ell}{T}
- \frac{\mu_{0,\ell}}{T_0} \right)\frac{\mu_i}{T} ds
\right\}\\ \nonumber
= & -\eps\int_\Omega\left(
|\nabla \frac{\mu_i}{T}|^2 + |\frac{\mu_i}{T}|^2
\right)dx
-\delta\int_\Omega |D^2 \frac{\mu_i}{T}|^2 dx\\
\nonumber
&-\eps\int_\Omega((1+T)\frac{|\nabla\log T|^2}{T} -
(1+T^{-K_3})\frac{\log T}{T})dx \\ \nonumber
& -\eps\int_\Omega |T|^{-1-K_3} |\nabla\log T|^{K_3+1}dx\\
\nonumber
& -\delta\int_\Omega( (1+T)
D^2 \log T : D^2\left( -\frac{1}{T} \right) + 
(T^{-1} + 1)|\nabla\log T|^4)dx .
\end{align}
We consider firstly the terms coming from the regularization and from the discrete time derivatives.
Note that the Gibbs-Duhem relations \eqref{GibbsDuhem} give:
\[ \frac{\pa(-\rho\eta)}{\pa\rho_i} =  \frac{\mu_i}{T}, \quad
  \frac{\pa(-\rho\eta)}{\pa (\rho e)} = -\frac{1}{T}.  \]
In this way we get:
\begin{align*}
&\Phi (S\rho_i - S^{k-1}\rho_i^{k-1})\frac{\mu_i}{T} 
-( \Phi(E_{f,\eps} - E_{f,\eps}^{k-1}) + 
(1-\Phi)(E_s - E_s^{k-1}) )\frac{1}{T} \\
&= S^{k-1}\Phi\left\{  (\rho_i - \rho_i^{k-1})\frac{\pa(-\rho\eta)}{\pa\rho_i} 
+ ((\rho e) - (\rho e)^{k-1})\frac{\pa(-\rho\eta)}{\pa (\rho e)}\right\} \\
&\quad + \Phi (S - S^{k-1})\left\{
\rho_i\frac{\mu_i}{T}
-\frac{(\rho e)}{T}\right\}\\
&\quad 
+ \frac{\Phi}{T}\int_{S^{k-1}}^S P_{c,\eps}(\xi)d\xi
-(1-\Phi)c_s\frac{T-T^{k-1}}{T} \\
&\quad 
-(1-\Phi)\eps(K_2-1)\frac{T^{K_2} - (T^{k-1})^{K_2} }{T}.
\end{align*}
The concavity of the functions
\begin{align*}
(\vec\rho,\rho e)\mapsto\rho\eta,\quad
T\mapsto\log T,\quad
T\mapsto T^{1-1/K_2},\quad
S\mapsto\int_{1/2}^S P_{c,\eps}(\xi)d\xi
\end{align*}
and \eqref{rhoeta} lead to
\begin{align*}
&\Phi (S\rho_i - S^{k-1}\rho_i^{k-1})\frac{\mu_i}{T} 
-( \Phi(E_{f,\eps} - E_{f,\eps}^{k-1}) + 
(1-\Phi)(E_s - E_s^{k-1}) )\frac{1}{T} \\
&\geq 
 \Phi S^{k-1}(\rho\eta)^{k-1} + (1-\Phi)
[c_s\log T^{k-1} + \eps K_2 (T^{k-1})^{K_2-1}]
-\Phi S(\rho\eta) \\
&\quad - (1-\Phi)
[c_s\log T + \eps K_2 T^{K_2-1} ]
+ \Phi (S-S^{k-1})\frac{P_{c,\eps}(S) + p}{T} .
\end{align*}
Multiplying the above inequality times $\sigma\tau^{-1}$ and employing
\eqref{lin.3} and \eqref{Skel.Entr} allows us to obtain
\begin{align}\label{d.eta}
&\frac{\sigma}{\tau}\left\{ \Phi (S\rho_i - S^{k-1}\rho_i^{k-1})\frac{\mu_i}{T} 
-( \Phi(E_{f,\eps} - E_{f,\eps}^{k-1}) + 
(1-\Phi)(E_s - E_s^{k-1}) )\frac{1}{T} \right\}\\
\nonumber
&\geq -\frac{\sigma}{\tau}  \Big\{
\Phi \Big( S(\rho\eta)   -S^{k-1}(\rho\eta)^{k-1}  \Big) + (1-\Phi) \Big(
     (\rho\eta)_s   -(\rho\eta)_s^{k-1} \Big)
\Big\}\\ \nonumber
&\qquad - \frac{\Phi}{T} \tau^{-2}(S-S^{k-1})(f(S)-f(S^{k-1})) 
 -  \frac{\sigma}{\tau} (1-\Phi) \eps K_2 [ T^{K_2 - 1} - (T^{k-1})^{K_2-1}    ]. 
\end{align}
Next, using Young's inequality, direct calculations give the following estimate:
\begin{align*}
\delta\int_\Omega( (1+T)
D^2 \log T : D^2\left( -\frac{1}{T} \right) + 
(T^{-1} + 1)|\nabla\log T|^4)dx\\
\geq c\delta\int_\Omega\left(
|D^2\log T|^2 
+ (T^{-1} + 1)|\nabla\log T|^4\right)dx
\end{align*}
for some constant $c>0$.

Therefore, from \eqref{app.ei.1}, using \eqref{Eq.9} and the discrete entropy balance equation, we obtain the global entropy balance for the approximate solutions:
\begin{align} \label{New.1}
&\sigma\int_\Omega\left\{ \sum_{i,j = 1}^N  L_{ij} \nabla \big(  \frac{\mu_i}{T}\big) \cdot \nabla \big(   \frac{\mu_j}{T} \big) + \frac{\kappa(T)|\nabla T|^2}{T^2}
+ \frac{K}{T}\lambda(S,T)|\nabla p|^2
-\vec r\cdot\frac{\vec\mu}{T}   \right\} dx  \\
& -\int_\Omega\Phi\tau^{-2}(S-S^{k-1})(f(S)-f(S^{k-1}))dx
\nonumber\\
& +\sigma\int_{\pa\Omega}\left(\alpha\frac{T_0 - T}{T}
+\sum_{i,j=1}^N b_{ij}\frac{\mu_i}{T}\left(\frac{\mu_j}{T} - \frac{\mu_{0,j}}{T_0}\right)\right)ds
\nonumber\\
&+\eps\int_\Omega\left(
|\nabla \frac{\mu_i}{T}|^2 + |\frac{\mu_i}{T}|^2
\right)dx
+\delta\int_\Omega |D^2 \frac{\mu_i}{T}|^2 dx\nonumber\\
&+\eps\int_\Omega((1+T)\frac{|\nabla\log T|^2}{T} -
(1+T^{-K_3})\frac{\log T}{T})dx \nonumber \\
&+\eps\int_\Omega |T|^{-1-K_3} |\nabla\log T|^{K_3+1}dx
\nonumber\\
& +c\delta\int_\Omega\left(
|D^2\log T|^2 
+ (T^{-1} + 1)|\nabla\log T|^4\right)dx\nonumber\\
& \leq \sigma\tau^{-1}\int_\Omega \Big(\Phi (S(\rho\eta) - S^{k-1}(\rho\eta)^{k-1})
+ (1-\Phi)((\rho \eta)_{s,\eps} - (\rho \eta)_{s,\eps}^{k-1} )  \Big) dx .
\nonumber
\end{align}
On the other hand, choosing $\psi=1$ in \eqref{lin.2} leads to the global energy balance inequality for the approximate solutions:
\begin{align} \label{New.2}
&\sigma\left\{
\int_\Omega\tau^{-1}( \Phi(E_{f,\eps} - E_{f,\eps}^{k-1}) + 
(1-\Phi)(E_s - E_s^{k-1}) ) dx
+\alpha\int_{\pa\Omega}(T - T_0) ds\right\}\\ \nonumber
&\qquad = -\eps\int_\Omega (1+T^{-K_3})\log T dx .
\end{align}
Summing \eqref{New.1} and \eqref{New.2} yields the global entropy-energy inequality for the approximate solutions:
\begin{align}\label{app.est}
&\sigma\tau^{-1}
\int_\Omega\left[ \Phi(E_{f,\eps} - S(\rho\eta)) + 
(1-\Phi)(E_s - (\rho\eta)_{s,\eps}) \right] dx\\ \nonumber
&+\sigma\int_\Omega\left\{ \sum_{i,j = 1}^N  L_{ij} \nabla \big(  \frac{\mu_i}{T}\big) \cdot \nabla \big(   \frac{\mu_j}{T} \big) + \frac{\kappa(T)|\nabla T|^2}{T^2}
+ \frac{K}{T}\lambda(S,T)|\nabla p|^2
-\vec r\cdot\frac{\vec\mu}{T}   \right\} dx  \\
& -\int_\Omega \frac{\Phi}{T}\tau^{-2}(S-S^{k-1})(f(S)-f(S^{k-1}))dx
\nonumber\\
& +\sigma\int_{\pa\Omega}\left(\alpha
\left(1 - \frac{1}{T}\right)(T-T_0)
+\sum_{i,j=1}^N b_{ij}\frac{\mu_i}{T}\left(\frac{\mu_j}{T} - \frac{\mu_{0,j}}{T_0}\right)\right)ds
\nonumber\\
&+\eps\int_\Omega\left(
\left|\nabla \frac{\vec\mu}{T}\right|^2 + \left|\frac{\vec\mu}{T}\right|^2
\right)dx
+\delta\int_\Omega \left|D^2 \frac{\vec\mu}{T}\right|^2 dx\nonumber\\
&+\eps\int_\Omega\left((1+T)\frac{|\nabla\log T|^2}{T} 
+(1+T^{-K_3})\left(1-\frac{1}{T}\right)\log T \right)dx \nonumber \\
&+\eps\int_\Omega |T|^{-1-K_3} |\nabla\log T|^{K_3+1}dx
\nonumber\\
& +c\delta\int_\Omega\left(
|D^2\log T|^2 
+ (T^{-1} + 1)|\nabla\log T|^4\right)dx\nonumber\\
& \leq \sigma\tau^{-1}\int_\Omega \left[\Phi 
(E_{f,\eps}^{k-1} - S^{k-1}(\rho\eta)^{k-1})
+ (1-\Phi)(E_s^{k-1} - (\rho \eta)_{s,\eps}^{k-1} ) \right] dx .
\nonumber
\end{align}
Furthermore, since $\int_{s_\eps^*}^S P_{c,\eps}(\xi)d\xi\leq C$ thanks to the definition on $P_{c,\eps}$, 
we deduce from the definition of $E_{f,\eps}$ and \eqref{eng.1}, \eqref{entropy} that 
\begin{align*}
E_{f,\eps} - S(\rho\eta) &= 
S(\rho^\gamma + c_w\rho T + p_{at})
- \int_{s_\eps^*}^S P_{c,\eps}(\xi)d\xi
+ S \eps\rho^{K_1}\\
&+S\sum_{i=1}^N\rho_i\log\rho_i - c_w S\rho (\log T + 1)\\
&\geq -C + S\left(
\rho^\gamma + c_w\rho\left( T - \log T \right)
+ \sum_{i=1}^N\rho_i(\log\rho_i - c_w)
\right)+ S \eps\rho^{K_1}\\
&\geq S\rho^\gamma + S \eps\rho^{K_1} - C ,
\end{align*}
while \eqref{Skel.Entr}, \eqref{Skel.Ener} imply
\begin{align*}
E_s - (\rho\eta)_s &= c_s(T - \log T) + 1 - c_s
+\eps T^{K_2-1}( (K_2-1)T - K_2 ) \\
&\geq C(T + |\log T| + \eps T^{K_2}) - C'.
\end{align*}
It follows
\begin{align}\label{lb.ESeta}
\Phi(E_{f,\eps} &- S(\rho\eta)) + 
(1-\Phi)(E_s - (\rho\eta)_s) \\ \nonumber
&\geq 
c_1\Phi S\left(\rho^\gamma + \eps\rho^{K_1}\right)
+ c_2(1-\Phi)(T + |\log T| + \eps T^{K_2}) - C .
\end{align}
We deduce from \eqref{app.est}, \eqref{lb.ESeta} that 
(recall $w=\log T$, $\vec z=\vec\mu/T$)
\begin{align*}
\|\vec z\|_{H^2(\Omega;\R^N)} + 
\|D^2 w\|_{L^2(\Omega)} + 
\|\nabla w\|_{L^2(\Omega)} + \int_\Omega w (1-e^{-w}) dx
\leq C(\eps,\delta,\tau).
\end{align*}
Since $C>0$ exists such that $|w|\leq C(w(1-e^{-w})+1)$ for every $w\in\R$, we obtain through the Poincar\'e Lemma
\begin{align*}
\|\vec z\|_{H^2(\Omega;\R^N)} + \|w\|_{H^2(\Omega)}\leq C(\eps,\delta,\tau).
\end{align*}

\subsubsection{Conclusion: existence for the approximate problem.}
We conclude that the set in \eqref{set.fp} is bounded in $V_0$ (and a fortiori in 
$V$) uniformely w.r.t.~$\sigma\in [0,1]$. Leray-Schauder's fixed point theorem allows us to conclude that $(\vec z, w)\in V_0$
exists such that $F((\vec z, w),1)=(\vec z, w)$, i.e.~the approximate system has a solution $(\vec\mu^k,T^k,S^k)$ such that
$\vec\mu^k/T^k = \vec z\in H^2(\Omega;\R^N)$,
$\log T^k = w\in H^2(\Omega)$, $S\in L^\infty(\Omega)$,
$0<S<1$ a.e.~in $\Omega$. In particular 
$\vec\mu^k\in H^2(\Omega;\R^N)$, $T^k\in H^2(\Omega)$.
Furthermore \eqref{app.est} holds with $\sigma=1$.
We conclude this part by stating the following technical result, whose proof can be found in the Appendix.
\begin{lemma}\label{lem.estSk}
The following estimates for $S_k$ hold:
\begin{align}
\label{Sceps}
& S^k \geq \eps > 0\mbox{ a.e.~in }\Omega,~~k\geq 0,~~\mbox{provided that }0<\eps<f^{-1}(p_{at}/\lambda_0),\\
\label{est.f.k}
& \left|\int_\Omega f(S^k)dx\right| \leq 
\left|\int_\Omega f(S^{k-1})dx\right|
+\tau C(\eps)\int_\Omega( 1 + |p^k| )dx ,\quad k\geq 0,\\
\label{est.naf.k}
&(1-C\tau)
\int_\Omega (1+|\nabla f(S^k)|^2)^{2/3} dx + 
\tau \int_\Omega \frac{|f'(S^k)P_{c,\eps}'(S^k)| |\nabla S^k|^2}{
	(1+|\nabla f(S^{k})|^2)^{1/3}} dx
\\ \nonumber
&\qquad\leq \int_\Omega (1+|\nabla f(S^{k-1})|^2)^{2/3} dx
+ \tau\int_\Omega |\nabla p^k|^{4/3} dx,\quad k\geq 0 .
\end{align}
\end{lemma}



\subsection{Limit $\delta\to 0$.}
We denote all the variables with a $\delta$ apex to put their dependence on $\delta$ in evidence.

From \eqref{app.est} it follows
\begin{align*}
\sqrt{\eps}\|\frac{\vec\mu^{(\delta)}}{T^{(\delta)}}
\|_{H^1(\Omega)} 
+ \|\log T^{(\delta)}\|_{H^1(\Omega)}
+ \|(T^{(\delta)})^{\frac{\beta}{2}}\|_{H^1(\Omega)}
\leq C(\tau).
\end{align*}
The compact Sobolev's embedding $H^1(\Omega)\hookrightarrow L^{q}(\Omega)$ for every $q<6$ implies that 
$\vec\mu^{(\delta)}/T^{(\delta)}$, $\log T^{(\delta)}$,
$(T^{(\delta)})^{\beta/2}$ are strongly convergent in $L^{q}(\Omega)$
for every $q<6$, and in particular also a.e.~in $\Omega$.
From \eqref{chp.1} it follows, for $i=1,\ldots,N$:
\begin{align}\label{rhomu.1}
\log\rho_i^{(\delta)} + \frac{\gamma (\rho^{(\delta)})^{\gamma-1}}{T^{(\delta)}} 
+\eps K_1(\rho^{(\delta)})^{K_1-1}
= \frac{\mu_i^{(\delta)}}{T^{(\delta)}} + c_w\log T^{(\delta)} - 1.
\end{align}
Taking the exponential of both sides in \eqref{rhomu.1} and summing in $i=1,\ldots,N$ yields
\begin{align*}
\rho^{(\delta)}\exp\left(
\frac{\gamma (\rho^{(\delta)})^{\gamma-1}}{T^{(\delta)}} 
+\eps K_1(\rho^{(\delta)})^{K_1-1}
\right) = \sum_{i=1}^N\exp\left(
\frac{\mu_i^{(\delta)}}{T^{(\delta)}} + c_w\log T^{(\delta)} - 1
\right).
\end{align*}
The right-hand side of the above identity is a.e.~convergent in $\Omega$ to an a.e.~positive function, while $T^{(\delta)}\to T>0$ a.e.~in $\Omega$. By the implicit function theorem $\rho^{(\delta)}$ is a.e.~convergent in $\Omega$ to an a.e.~positive function $\rho$. 
From this fact and
\eqref{rhomu.1} it follows that $\rho_i^{(\delta)}$ is a.e.~convergent in $\Omega$ to $\rho_i>0$, for $i=1,\ldots,N$.

From \eqref{app.est}, \eqref{lb.ESeta}, \eqref{Sceps} we also deduce
\begin{equation}
\label{app.rho.K1}
\|\rho^{(\delta)}\|_{L^{K_1}(\Omega)}\leq C(\eps,\tau).
\end{equation}
In particular $\rho^{(\delta)}_i\to\rho_i$ strongly in $L^q(\Omega)$, for every $q<K_1$, $i=1,\ldots,N$.

From \eqref{app.est} and \eqref{mobility} (as well as the boundedness of $\mu$) it follows that 
$$
\left\|\sqrt\frac{k_r(S^{(\delta)})}{T^{(\delta)}} \nabla p^{(\delta)}
\right\|_{L^2(\Omega)}\leq C(\tau) .
$$
However, \eqref{Sceps} and the $L^6(\Omega)$ bound for $(T^{(\delta)})^{\beta/2}$ yield
\begin{align*}
\|\nabla p^{(\delta)}\|_{L^{\frac{6\beta}{1+3\beta}}(\Omega)}
&\leq
\|\sqrt{T^{(\delta)}}\|_{L^{6\beta}(\Omega)}
\left\| \frac{1}{\sqrt{T^{(\delta)}}}\nabla p^{(\delta)} \right\|_{L^2(\Omega)}\\
&\leq C(\eps) \|\sqrt{T^{(\delta)}}\|_{L^{6\beta}(\Omega)}
\left\|\sqrt\frac{k_r(S^{(\delta)})}{T^{(\delta)}} \nabla p^{(\delta)}
\right\|_{L^2(\Omega)}\leq C(\eps,\tau).
\end{align*}
The above estimate, together with the bound for $p^{(\delta)}$ in $L^1(\Omega)$ which comes from the $L^1(\Omega)$ estimate for $E_{f,\eps}^{(\delta)}$ and relation 
$p^{(\delta)}\leq C E_{f,\eps}^{(\delta)} + C' $
implies (via Poincar\'e's Lemma)
\begin{align}\label{est.p.1}
\|p^{(\delta)}\|_{W^{1,\frac{6\beta}{1+3\beta}}(\Omega)}\leq C(\eps,\tau).
\end{align}
Therefore via Sobolev's embedding
\begin{align}\label{est.p.2}
\|p^{(\delta)}\|_{L^{\frac{6\beta}{1+\beta}}(\Omega)}\leq C(\eps,\tau),\\ \nonumber
p^{(\delta)}\to p\quad \mbox{strongy in }L^q(\Omega),~~\forall q<\frac{6\beta}{1+\beta}.
\end{align}
Furthermore the trivial bound $(\rho^{(\delta)})^\gamma\leq C(p^{(\delta)}+1)$ and the a.e.~convergence of $\rho_i^{(\delta)}$ imply
\begin{align*}
\rho_i^{(\delta)}\to\rho_i\quad\mbox{strongly in }L^q(\Omega),\quad
\forall q<\frac{6\beta\gamma}{1+3\beta},\quad i=1,\ldots,N.
\end{align*}
The constitutive relations between all the variables (see subsection \ref{sec.model}) hold because $\rho_i^{(\delta)}$, $T^{(\delta)}$ are a.e.~convergent in $\Omega$ towards a.e.~positive limits.

Finally, \eqref{app.3}, the strong convergence of $p^{(\delta)}$ and the strict monotonicity of $f$, $P_{c,\eps}$ imply that $S^{(\delta)}$ is a.e.~convergent in $\Omega$ towards a limit $S$
which satisfies $0<\eps\leq S < 1$ a.e.~in $\Omega$ due to \eqref{Sceps} and $\lim_{s\to 1}f(s)=-\infty$. Furthermore $S^{(\delta)}$ is also strongly convergent to $S$ in $L^q(\Omega)$ for every $q<\infty$.

It is immediate to see that
the $\delta-$depending regularization terms in \eqref{app.1}, \eqref{app.2} easily tend to 0 when $\delta\to 0$, the only exception being the term 
\begin{align*}
\delta\int_\Omega T^{(\delta)} |\nabla\log T^{(\delta)}|^2 \nabla\log T^{(\delta)}\cdot\nabla\psi dx
\end{align*}
which appears in \eqref{app.2}.
Since $\beta\geq 4/3$ by assumption, it follows
\begin{align*}
\delta\int_\Omega T^{(\delta)} |\nabla\log T^{(\delta)}|^3 dx &\leq \delta\|T^{(\delta)}\|_{L^{4}(\Omega)}
\|\nabla\log T^{(\delta)}\|_{L^{4}(\Omega)}^3\\
&\leq C\delta\|T^{(\delta)}\|_{L^{3\beta}(\Omega)}
\|\nabla\log T^{(\delta)}\|_{L^{4}(\Omega)}^3 .
\end{align*}
Since $T^{(\delta)}$ is bounded in $L^{3\beta}(\Omega)$ while
$\delta^{1/4}\nabla\log T^{(\delta)}$
is bounded in $L^4(\Omega)$, we conclude
\begin{align*}
\delta\int_\Omega T^{(\delta)}|\nabla\log T^{(\delta)}|^3 dx\leq C\delta^{1/4}.
\end{align*}
Let us now check whether the terms in the energy equation are bounded in $L^q(\Omega)$ for some $q>1$. 
From \eqref{est.p.1}, \eqref{vow} and the boundedness of $\lambda(S^{(\delta)},T^{(\delta)})$ it follows
\begin{equation}
\label{est.v.1}
\|\veb v^{(\delta)}\|_{L^{\frac{6\beta}{1+3\beta}}(\Omega)}\leq C,
\end{equation}
and so
\begin{align*}
\veb v^{(\delta)}\rightharpoonup \veb v\quad
\mbox{weakly in }L^{\frac{6\beta}{1+3\beta}}(\Omega).
\end{align*}
Furthermore \eqref{est.p.2} and the strong convergence of $S^{(\delta)}$, $T^{(\delta)}$ imply
\begin{align*}
\nabla p^{(\delta)}\rightharpoonup\nabla p\quad
\mbox{weakly in }L^{\frac{6\beta}{1+3\beta}}(\Omega),\\
\lambda(S^{(\delta)},T^{(\delta)})\to\lambda(S,T)\quad\mbox{strongly in }L^q(\Omega),\quad\forall q<\infty.
\end{align*}
It follows $\veb v = -K\lambda(S,T)\nabla p$.

From \eqref{est.p.2}, \eqref{est.v.1} we deduce
\begin{align*}
p^{(\delta)}\veb v^{(\delta)}\rightharpoonup p\veb v
\quad\mbox{weakly in }L^{\frac{3\beta}{1+2\beta}}(\Omega).
\end{align*}
Notice $\frac{3\beta}{1+2\beta}>1$ since $\beta>1$.
Also, since $(\rho e)^{(\delta)}\leq C p^{(\delta)} + C'$,
it follows
\begin{align*}
(\rho e)^{(\delta)}\veb v^{(\delta)}
\rightharpoonup (\rho e)\veb v
\quad\mbox{weakly in }L^{\frac{3\beta}{1+2\beta}}(\Omega).
\end{align*}
Let us consider the heat flux $\veb q^{(\delta)}$.
Since
\begin{align*}
\veb q^{(\delta)} = -\kappa(T^{(\delta)})\nabla T^{(\delta)} + \sum_{j=1}^N L_{0j}\nabla\frac{\mu_j^{(\delta)}}{T^{(\delta)}} 
\end{align*}
and $\mu_j^{(\delta)}/T^{(\delta)} \rightharpoonup \mu_j/T$ weakly in $L^2(\Omega)$,
for $j=1,\ldots,N$, one only needs to make sure that 
$\kappa(T^{(\delta)})\nabla T^{(\delta)}$ converges (weakly) to the correct limit. From 
the $\delta-$uniform bound for $\nabla (T^{(\delta)})^{K_4/2}$ in $L^2(\Omega)$ it follows
\begin{align*}
\|\kappa(T^{(\delta)})\nabla T^{(\delta)}\|_{L^{\frac{3\beta}{2\beta+1}}(\Omega)}
&\leq C\| (T^{(\delta)})^{\beta/2 + 1}\nabla (T^{(\delta)})^{\beta/2} \|_{L^{\frac{3\beta}{2\beta+1}  }(\Omega)}\\
&\leq
C\|(T^{(\delta)})^{\beta/2 + 1}\|_{L^{\frac{6\beta}{\beta + 2}}(\Omega)}
\|\nabla (T^{(\delta)})^{\beta/2}\|_{L^2(\Omega)}\\
&\leq C\|T^{(\delta)}\|_{
	L^{3\beta}(\Omega)}^{\beta/2 + 1}
\|\nabla (T^{(\delta)})^{\beta/2}\|_{L^2(\Omega)}
\leq C(\eps,\tau).
\end{align*}
Therefore
\begin{align*}
\kappa(T^{(\delta)})\nabla T^{(\delta)}\rightharpoonup
\kappa(T)\nabla T\quad\mbox{weakly in }L^{\frac{3\beta}{2\beta+1}}(\Omega).
\end{align*}
It is immediate to see that
$E_{f,\eps}^{(\delta)}$ is strongly convergent (in $L^1(\Omega)$) towards the correct limit.
At this point taking the limit $\delta\to 0$ in \eqref{app.1}--\eqref{app.3} is straightforward. The resulting limit equations are just like \eqref{app.1}--\eqref{app.3} but without the regularizing terms in $\delta$.

\subsection{Limit $\tau\to 0$.}
Let us introduce the new notation
$$\rho_i^{(\tau)}(x,t) = \rho_i^{in}(x)\chf{\{0\}}(t) + 
\sum_{k=1}^N\rho_i^k(x)\chf{ ((k-1)\tau,k\tau] }(t),\quad
i=1,\ldots,N, $$
and the same for the other variables. We introduce also the discrete time derivative
\begin{align*}
D_\tau f(t) = \tau^{-1}( f(t) - f(t-\tau) ),\quad t\in [\tau,\mathcal{T}],
\end{align*}
for any function $f : [0,\mathcal{T}]\to X$ ($X$ is any vector space on $\R$). We can rewrite \eqref{app.1}--\eqref{app.3} (after the limit $\delta\to 0$) as
\begin{align}
\label{app.1b}
&\int_0^\mathcal{T}\int_\Omega D_\tau(S^{(\tau)}\rho^{(\tau)}_i)
\Phi\varphi dx dt
-\int_0^\mathcal{T}\int_\Omega( \rho^{(\tau)}_i\mathbf{v}^{(\tau)} + \mathbf{J}^{(\tau)}_i )\cdot\nabla\varphi dx dt\\ \nonumber
&\qquad
-\int_0^\mathcal{T}\int_\Omega r_{i,\eps}^{(\tau)}\varphi dx dt
+\int_0^\mathcal{T}\int_{\pa\Omega}\sum_{\ell=1}^N b_{i\ell}\left( \frac{\mu^{(\tau)}_\ell}{T^{(\tau)}}
- \frac{\mu_{0,\ell}}{T_0} \right)\varphi ds dt\\ \nonumber
&\quad = -\eps\int_0^\mathcal{T}\int_\Omega\left(
\nabla\frac{\mu^{(\tau)}_i}{T^{(\tau)}}\cdot\nabla\varphi + 
\frac{\mu^{(\tau)}_i}{T^{(\tau)}}\varphi\right)dx dt,
\end{align}
\begin{align}
\label{app.2b}
&\int_0^\mathcal{T}\int_\Omega D_\tau( \Phi E_{f,\eps}^{(\tau)} + 
(1-\Phi)E_s^{(\tau)} )\psi dx dt\\ \nonumber
&\quad
-\int_0^\mathcal{T}\int_\Omega( ((\rho e)^{(\tau)} + p^{(\tau)})\mathbf{v}^{(\tau)} + \veb q^{(\tau)} )\cdot\nabla\psi dx dt
+\alpha\int_0^\mathcal{T}\int_{\pa\Omega}(T^{(\tau)} - T_0)\psi ds dt\\ \nonumber
& = -\eps\int_0^\mathcal{T}\int_\Omega((1+T^{(\tau)})\nabla\logT^{(\tau)}\cdot\nabla\psi + (1+(T^{(\tau)})^{-K_3})(\logT^{(\tau)})\psi)dx dt\\
\nonumber
&\quad -\eps\int_0^\mathcal{T}\int_\Omega (T^{(\tau)})^{-K_3} |\nabla\log T^{(\tau)}|^{K_3-1}\nabla\logT^{(\tau)}\cdot\nabla\psi dxdt,
\end{align}
\begin{equation}
\label{app.3b}
D_\tau f(S^{(\tau)}) + P_{c,\eps}(S^{(\tau)}) + p^{(\tau)} = 0.
\end{equation}
Estimate \eqref{app.est} can be rewritten as 
\begin{align}\label{app.est.2}
&\sup_{[0,\mathcal{T}]}\int_\Omega\left[ \Phi(E_{f,\eps}^{(\tau)} - S^{(\tau)}(\rho\eta)^{(\tau)}) + 
(1-\Phi)(E_s^{(\tau)} - (\rho\eta)_s^{(\tau)}) \right] dx\\ \nonumber
&+\int_0^\mathcal{T}\int_\Omega\left\{\sum_{i,j = 1}^N  
L_{ij}\nabla \big(  \frac{\mu^{(\tau)}_i}{T^{(\tau)}}\big) \cdot \nabla \big(\frac{\mu^{(\tau)}_j}{T^{(\tau)}}\big) 
+ \frac{\kappa(T^{(\tau)})|\nablaT^{(\tau)}|^2}{|T^{(\tau)}|^2}\right.\\
\nonumber
&\left. + \frac{K}{T^{(\tau)}}\lambda(S^{(\tau)},T^{(\tau)})|\nablap^{(\tau)}|^2
-\vec r^{(\tau)}\cdot\frac{\vec\mu^{(\tau)}}{T^{(\tau)}}
-\Phi (D_\tauS^{(\tau)})(D_\tau f(S^{(\tau)}))
\right\}dx  \\
& +\int_0^\mathcal{T}\int_{\pa\Omega}\left[\alpha
\left(1 - \frac{1}{T^{(\tau)}}\right)(T^{(\tau)}-T_0)
+\sum_{i,j=1}^N b_{ij}\frac{\mu_i^{(\tau)}}{T^{(\tau)}}\left(\frac{\mu_j^{(\tau)}}{T^{(\tau)}} - \frac{\mu_{0,j}}{T_0}\right)\right]ds dt
\nonumber\\
&+\eps\int_0^\mathcal{T}\int_\Omega\left(
\left|\nabla \frac{\vec\mu^{(\tau)}}{T^{(\tau)}}\right|^2 + \left|\frac{\vec\mu^{(\tau)}}{T^{(\tau)}}\right|^2
+(1+T^{(\tau)})\frac{|\nabla\log T^{(\tau)}|^2}{T^{(\tau)}} 
\right.\nonumber\\
&\left. \qquad
+\left|\nabla\frac{1}{T^{(\tau)}}\right|^{K_3+1}
+(1+(T^{(\tau)})^{-K_3})
\left(1-\frac{1}{T^{(\tau)}}\right)\log T^{(\tau)} \right)dxdt \nonumber \\
& \leq \int_\Omega \left[\Phi 
(E_{f,\eps}^{in} - S^{in,\eps}(\rho\eta)^{in})
+ (1-\Phi)(E_s^{in} - (\rho \eta)_s^{in} ) \right] dx .
\nonumber
\end{align}
The a priori estimates derived in Section \ref{sec.apriori} are satisfied by the approximate solution. Moreover, the latter satisfies also the following bounds, which follow from \eqref{app.est.2}:
\begin{align}\label{est.nablamutau}
\sqrt{\eps}\|\mu^{(\tau)}_i/T^{(\tau)}\|_{L^2(0,\mathcal{T}; H^1(\Omega))}
+\eps^{1/a}\|\mu^{(\tau)}_i/T^{(\tau)}\|_{L^a(Q_\mathcal{T})}
\leq C,\quad i=1,\ldots,N.
\end{align}
Let us now find a better estimate for $T^{(\tau)}$. The 
$L^\infty(0,\mathcal{T}; L^1(\Omega))$ bound for $E_{s,\eps}$ (see \eqref{app.est.2}) as well as \eqref{Es.eps} lead to
\begin{equation}
\label{est.TK2.eps}
\|T^{(\tau)}\|_{L^\infty(0,\mathcal{T}; L^{K_2}(\Omega))}\leq C\eps^{-1/K_2},
\end{equation}
while the $L^2(Q_\mathcal{T})$ bound for $\nabla (T^{(\tau)})^{\beta/2}$ and the Sobolev embedding $H^1\hookrightarrow L^6$ imply
\begin{align*}
\|T^{(\tau)}\|_{L^{\beta}(0,\mathcal{T}; L^{3\beta}(\Omega))}\leq C.
\end{align*}
The previous two estimates yield by interpolation
(we also employed the fact that $\beta + (2/3)K_2>\beta+2$, which follows from \eqref{hp.K123})
\begin{align}\label{est.Ttau}
\|T^{(\tau)}\|_{L^{\beta+2}(Q_\mathcal{T})}\leq C
\|T^{(\tau)}\|_{L^{\beta + (2/3)K_2 }(Q_\mathcal{T})} &\leq
C \|T^{(\tau)}\|_{L^\infty(0,\mathcal{T}; L^{K_2}(\Omega))}^{
\frac{(2/3)K_2}{\beta + (2/3)K_2}}
\|T^{(\tau)}\|_{L^{\beta}(0,\mathcal{T}; L^{3\beta}(\Omega))}^{
\frac{\beta}{\beta + (2/3)K_2}}\\ \nonumber
&\leq C \eps^{-\frac{(2/3)}{\beta + (2/3)K_2}}.
\end{align}
Furthermore, from \eqref{app.est.2} and the fact that
\begin{align*}
(T-1)\log T\geq 0\quad\forall T>0,\qquad
\inf_{0<T<1/2}(T-1)\log T = \frac{1}{2}\log 2 > 0,
\end{align*}
we also deduce
\begin{equation}
\label{est.Ttau.2}
\eps\int_{Q_\mathcal{T}\cap\{T^{(\tau)} < 1\}}
(T^{(\tau)})^{-K_3-1}\log\frac{1}{T^{(\tau)}}dx dt + 
\eps^{\frac{1}{1+K_3}}
\left\| \frac{1}{T^{(\tau)}} \right\|_{L^{1+K_3}(0,\mathcal{T}; W^{1,1+K_3}(\Omega))}\leq C .
\end{equation}
Let us now find an estimate for the pressure.
From \eqref{AJPP.1}, 
\eqref{Sceps}, \eqref{est.TK2.eps} it follows
\begin{align*}
\|\nablap^{(\tau)}\|_{L^{2}(0,\mathcal{T}; L^{\frac{2K_2}{K_2+1}}(\Omega) )}
&\leq C(\eps)\|\sqrt{T^{(\tau)}}\|_{L^{\infty}(0,\mathcal{T}; L^{2K_2}(\Omega))}
\left\|\sqrt\frac{k_r(S^{(\tau)})}{T^{(\tau)}}
\nabla p^{(\tau)}\right\|_{L^2(Q_\mathcal{T})}\leq C(\eps),
\end{align*}
so via the $L^\infty(0,T; L^1(\Omega))$ bound for $p^{(\tau)}$ in \eqref{Est.p} and Poincar\'e's Lemma
\begin{equation}
\label{est.naptau.K2}
\|p^{(\tau)}\|_{L^{2}(0,\mathcal{T}; W^{1,\frac{2K_2}{K_2+1}}(\Omega) )}\leq
C(\eps).
\end{equation}
Sobolev's embedding $W^{1,\frac{2K_2}{K_2+1} }(\Omega)\hookrightarrow L^{\frac{6K_2}{K_2+3}}(\Omega)$
implies
\begin{align*}
\|p^{(\tau)}\|_{L^2(0,T; L^{\frac{6K_2}{K_2+3}}(\Omega))}\leq C(\eps).
\end{align*}
The estimate above together with the $L^\infty(0,\mathcal{T}; L^1(\Omega))$ bound for $p^{(\tau)}$ coming from \eqref{app.est.2} 
leads (by interpolation) to
\begin{align}
\label{est.ptau.K2}
\|p^{(\tau)}\|_{L^{2/\theta}(0,\mathcal{T}; L^\xi(\Omega))}\leq C(\eps),
\quad \frac{1}{\xi} = 1-\theta + \theta\frac{K_2+3}{6K_2},
\quad \frac{3(K_2+1)}{5K_2-3}<\theta<1.
\end{align}
We point out that 
thanks to \eqref{hp.K123} it holds that $K_2>3$, which implies
$\frac{3(K_2+1)}{5K_2-3}<1$. Therefore one can choose
$\frac{3(K_2+1)}{5K_2-3}< \theta < 1$ in \eqref{est.ptau.K2}.

Another estimate for $p^{(\tau)}$ can be found by combining \eqref{AJPP.1} and \eqref{est.Ttau}. Noticing that 
$\beta+2\geq 3$ and \eqref{Sceps}, \eqref{hp.kr} hold, one finds
\begin{align*}
\|\nablap^{(\tau)}\|_{L^{3/2}(Q_\mathcal{T})}
&\leq C\eps^{-\alpha_r}\|\sqrt{T^{(\tau)}}\|_{L^{6}(Q_\mathcal{T})}
\left\|\sqrt\frac{k_r(S^{(\tau)})}{T^{(\tau)}}
\nabla p^{(\tau)}\right\|_{L^2(Q_\mathcal{T})}
\leq C \eps^{-\alpha_r} ,
\end{align*}
so via the $L^\infty(0,T; L^1(\Omega))$ bound for $p^{(\tau)}$ in \eqref{Est.p} and Poincar\'e's Lemma
\begin{equation}
\label{est.ptau.1}
\|p^{(\tau)}\|_{L^{\frac{3}{2}}(0,\mathcal{T};
	W^{1,\frac{3}{2}}(\Omega))}\leq C \eps^{-\alpha_r} .
\end{equation}
Sobolev's embedding $W^{1,3/2}(\Omega)\hookrightarrow L^{3}(\Omega)$ yields
\begin{align*}
\|p^{(\tau)}\|_{L^{3/2}(0,\mathcal{T};
	L^{3}(\Omega)}\leq C \eps^{-\alpha_r},
\end{align*}
which, together with the $L^\infty(0,\mathcal{T}; L^1(\Omega))$ bound for $p^{(\tau)}$ coming from \eqref{Est.p}, allows us to deduce by interpolation that
\begin{equation}
\label{est.ptau}
\|p^{(\tau)}\|_{L^{2}(Q_\mathcal{T})}\leq 
\|p^{(\tau)}\|_{L^\infty(0,\mathcal{T}; L^1(\Omega))}^{1/4}
\|p^{(\tau)}\|_{L^{3/2}(0,\mathcal{T}; L^3(\Omega))}^{3/4}
\leq C\eps^{-\frac{3}{4}\alpha_r} .
\end{equation}
Inequality \eqref{est.naf.k} implies in the new notation (upon summation in $k$)
\begin{align*}
\int_\Omega & (1+|\nabla f(S^{(\tau)}(t))|^2)^{2/3} dx + 
\int_0^t\int_\Omega \frac{|f'(S^{(\tau)})P_{c,\eps}'(S^{(\tau)})| |\nablaS^{(\tau)}|^2}{
	(1+|\nabla f(S^{(\tau)})|^2)^{1/3}} dx dt'
\\ \nonumber
&\leq \int_\Omega (1+|\nabla f(S^{in,\eps})|^2)^{2/3} dx
+ \int_0^t\int_\Omega |\nablap^{(\tau)}|^{4/3} dx dt'\\
&\qquad + C\int_0^t\int_\Omega (1+|\nabla f(S^{(\tau)})|^2)^{2/3} dx dt ,\qquad t\in (0,\mathcal{T}] .
\end{align*}
Gronwall's inequality (in its discrete version) and  \eqref{est.ptau.1} 
lead to
\begin{align}\label{est.naf}
\|\nabla f(S^{(\tau)})\|_{L^\infty(0,\mathcal{T}; L^{\frac{4}{3} }(\Omega))}^{4/3}
+ \int_0^\mathcal{T}\int_\Omega \frac{|f'(S^{(\tau)})P_{c,\eps}'(S^{(\tau)})| |\nablaS^{(\tau)}|^2}{
	(1+|\nabla f(S^{(\tau)})|^2)^{1/3}} dx dt
\leq C(\eps).
\end{align}
Thanks to \eqref{lb.dfeps}, the above estimate implies
\begin{equation}
\label{est.naStau}
\|\nabla S^{(\tau)}\|_{L^\infty(0,\mathcal{T}; L^{4/3}(\Omega))}\leq C(\eps).
\end{equation}
From \eqref{est.f.k}, \eqref{Est.p} as well as the fact that $f$ is upper bounded, one easily deduces that 
\begin{align}\label{est.f}
\|f(S^{(\tau)})\|_{L^\infty(0,\mathcal{T}; L^{1}(\Omega))}
\leq C(\eps).
\end{align}
From \eqref{est.naf}, \eqref{est.f} and Poincar\'e inequality
we get
\begin{equation}
\label{est.fStau.W}
\|f(S^{(\tau)})\|_{L^\infty(0,\mathcal{T}; W^{1,4/3}(\Omega))}\leq
C(\eps).
\end{equation}
Furthermore, thanks to \eqref{Sceps} 
and the fact that $P_{c,\eps}$ is nonnegative, decreasing and continuous in $(0,1)$, we immediately deduce that
\begin{equation}
\label{est.PcStau}
\|P_{c,\eps}(S^{(\tau)})\|_{L^{\infty}(Q_\mathcal{T})}\leq C(\eps).
\end{equation}
From \eqref{app.3b}, \eqref{est.ptau}, \eqref{est.PcStau}, 
it follows
\begin{equation}
\label{est.Df}
\|D_\tau f(S^{(\tau)})\|_{L^2(Q_\mathcal{T})}\leq C(\eps).
\end{equation}
Bounds \eqref{est.fStau.W}, \eqref{est.Df} allow us to apply Aubin-Lions Lemma in the version of \cite[Thr.~1]{DreJue12}
and obtain (up to subsequences)
\begin{equation}
\label{cnv.fStau}
f(S^{(\tau)})\to f^*\qquad\mbox{strongly in $L^1(Q_\mathcal{T})$, as }\tau\to 0.
\end{equation}
Being $f$ invertible, this implies that
$S^{(\tau)}$ is a.e.~convergent in $Q_\mathcal{T}$. The uniform $L^\infty$ bounds for $S^{(\tau)}$ yield
\begin{equation}
\label{cnv.Stau}
S^{(\tau)}\to S\quad\mbox{strongly in }L^q(Q_\mathcal{T}),~~\forall q<\infty.
\end{equation}
Given \eqref{Sceps}, \eqref{cnv.fStau}, \eqref{cnv.Stau} and the fact that 
$\lim_{s\to 1}f(s)=-\infty$, Fatou's Lemma allows us to state
\begin{equation}
\label{Seps.bounds}
\eps\leq S < 1\qquad\mbox{a.e.~in }Q_\mathcal{T} .
\end{equation}
From \eqref{p.eps}, \eqref{est.ptau} it follows 
\begin{align}\label{est.rhotau}
\|\rho^{(\tau)}\|_{L^{2 K_1}(Q_\mathcal{T})}\leq C \eps^{-\frac{3\alpha_r}{4K_1}} .
\end{align}
Let us now find an estimate for  $\nabla\rho_i$, $i=1,\ldots,N$.
Differentiating \eqref{rhomu.1} in $x$ leads to
\begin{align}\label{na.rhotau}
\nabla\log\rho^{(\tau)}_i 
+ \gamma\frac{\nabla (\rho^{(\tau)})^{\gamma-1}}{T^{(\tau)}}
- \gamma \frac{(\rho^{(\tau)})^{\gamma-1}}{T^{(\tau)}}\nabla\logT^{(\tau)}
+\eps K_1\nabla(\rho^{(\tau)})^{K_1-1}\\ \nonumber
= \nabla\frac{\mu^{(\tau)}_i}{T^{(\tau)}} + c_w\nabla\logT^{(\tau)},\quad
i=1,\ldots,N .
\end{align}
Easy computations allow us to equivalently rewrite the above identity as
\begin{align*}
\sum_{j=1}^N\mathcal M_{ij}\nabla\rho^{(\tau)}_j = \mathcal F_i,\qquad
i=1,\ldots,N,\\
\mathcal M_{ij} = \frac{\delta_{ij}}{\rho^{(\tau)}_i} + 
\frac{\gamma(\gamma-1)}{T^{(\tau)}}(\rho^{(\tau)})^{\gamma-2}
+\eps K_1 (K_1-1)(\rho^{(\tau)})^{K_1-2},\\
\mathcal F_i = \nabla\frac{\mu^{(\tau)}_i}{T^{(\tau)}} + 
\left(c_w + \frac{\gamma (\rho^{(\tau)})^{\gamma-1}}{T^{(\tau)}} \right)
\nabla\logT^{(\tau)} .
\end{align*}
Since $\sum_{i,j=1}^N\mathcal M_{ij}v_i v_j\geq\sum_{i=1}^N \frac{v_i^2}{\rho^{(\tau)}_i}$, we deduce
\begin{align*}
4\sum_{i=1}^N |\nabla\sqrt{\rho^{(\tau)}_i}|^2\leq 
\sum_{i,j=1}^N\mathcal M_{ij}\nabla\rho^{(\tau)}_i\cdot\nabla\rho^{(\tau)}_j = 
\sum_{i=1}^N\mathcal F_i\cdot\nabla\rho^{(\tau)}_i\\
\leq
2\left( \sum_{i=1}^N |\mathcal F_i|^2\rho_i \right)^{1/2}
\left( \sum_{i=1}^N |\nabla\sqrt{\rho^{(\tau)}_i}|^2 \right)^{1/2},
\end{align*}
which implies
\begin{align*}
\sum_{i=1}^N |\nabla\sqrt{\rho^{(\tau)}_i}| &\leq
C \sum_{i=1}^N |\mathcal F_i|\sqrt{\rho^{(\tau)}_i}\\
&\leq C\sqrt{\rho^{(\tau)}}\left|\nabla\frac{\vec\mu^{(\tau)}}{T^{(\tau)}}
\right| + 
C\left(1 + \frac{(\rho^{(\tau)})^{\gamma-1/2}}{T^{(\tau)}}\right)|\nabla\logT^{(\tau)}|.
\end{align*}
Multiplying the above estimate times $\sqrt{\rho^{(\tau)}}$ leads to
\begin{align}
\label{ub.narho}
\sum_{i=1}^N |\nabla\rho^{(\tau)}_i| &\leq
C\rho^{(\tau)}\left|\nabla\frac{\vec\mu^{(\tau)}}{T^{(\tau)}}
\right| + 
C\left(1 + \frac{(\rho^{(\tau)})^{\gamma}}{T^{(\tau)}}\right)|\nabla\logT^{(\tau)}|
\end{align}
From \eqref{hp.K123}, \eqref{est.rhotau}, \eqref{est.Ttau.2} it follows
\begin{align*}
\left\|
\rho^{(\tau)}\nabla\frac{\vec\mu^{(\tau)}}{T^{(\tau)}}
\right\|_{L^1(Q_\mathcal{T})}\leq 
\left\|\rho^{(\tau)}\right\|_{L^2(Q_\mathcal{T})}
\left\|\nabla\frac{\vec\mu^{(\tau)}}{T^{(\tau)}}
\right\|_{L^2(Q_\mathcal{T})}\leq C(\eps),
\end{align*}
as well as
\begin{align*}
\left\|
\frac{(\rho^{(\tau)})^{\gamma}}{T^{(\tau)}}\nabla\logT^{(\tau)}
\right\|_{L^1(Q_\mathcal{T})} &\leq
C\left\|(\rho^{(\tau)})^{\gamma}\right\|_{L^{ 
		\frac{5K_1}{3\gamma} }(Q_\mathcal{T})}
\left\|\frac{1}{T^{(\tau)}}
\right\|_{L^{1+K_3}(Q_\mathcal{T})}
\left\|\nabla\logT^{(\tau)}
\right\|_{L^2(Q_\mathcal{T})}\\
&\leq C(\eps).
\end{align*}
We deduce
\begin{equation}
\label{est.narhotau}
\|\nabla\rho^{(\tau)}_i\|_{L^1(Q_\mathcal{T})}\leq C(\eps),\quad i=1,\ldots,N.
\end{equation}
From \eqref{gamma} and \eqref{hp.K123}, \eqref{est.rhotau} it follows $\|\rho^{(\tau)}_i\|_{L^4(Q_\mathcal{T})}\leq C(\eps)$ for $i=1,\ldots,N$. From this fact, the uniform $L^\infty$ bounds for $S^{(\tau)}$, as well as \eqref{est.naStau}, \eqref{est.narhotau}, we conclude
$$
\|\nabla (S^{(\tau)}\rho^{(\tau)}_i)\|_{L^1(Q_\mathcal{T})}\leq C(\eps),\qquad
i=1,\ldots,N,
$$
which, thanks to the uniform $L^1$ bound for $S^{(\tau)}\rho^{(\tau)}$,
yields
\begin{equation}
\label{est.Srhotau}
\|S^{(\tau)}\rho^{(\tau)}_i\|_{L^1(0,\mathcal{T}; W^{1,1}(\Omega))}\leq C(\eps),\qquad
i=1,\ldots,N.
\end{equation}
An estimate for $D_\tau (S^{(\tau)}\rho^{(\tau)}_i)$,
$i=1,\ldots,N$, can be derived from \eqref{app.1b} in a rather straightforward way, just like we did in the previous section (devoted to the limit $\delta\to 0$), although now one has to keep in mind the dependency of every quantity on time. 
The following bounds are easily derived from 
\eqref{hp.Ltilde}, \eqref{hp.r.2}, \eqref{Est.A}, \eqref{Est.B.2}:
\begin{align*}
\|\veb J_i^{(\tau)}\|_{L^2(Q_\mathcal{T})} + \|r_{i,\eps}^{(\tau)}\|_{L^{a/(a-1)}(Q_\mathcal{T})} \leq C,\qquad i=1,\ldots,N.
\end{align*}
On the other hand, from \eqref{vow}, \eqref{est.ptau.1}, \eqref{est.rhotau} and the fact that $K_1>\frac{6}{5}\gamma>\frac{12}{5}$ (see \eqref{hp.K123}, \eqref{vow})
it follows
\begin{align*}
\exists q>1:\quad
\|\rho^{(\tau)}_i\veb v^{(\tau)}\|_{L^{q}(Q_\mathcal{T})}\leq C
\|\rho^{(\tau)}\|_{L^{\frac{5}{3}K_1}(Q_\mathcal{T})}
\|\nablap^{(\tau)}\|_{L^{\frac{4}{3}}(Q_\mathcal{T})}\leq C(\eps).
\end{align*}
We deduce
\begin{align*}
\|D_\tau(S^{(\tau)}\rho^{(\tau)}_i\Phi)\|_{L^q(0,\mathcal{T}; W^{1,q}(\Omega)')}\leq C(\eps),\quad i=1,\ldots,N.
\end{align*}
Since $\Phi\in L^\infty(0,\mathcal{T}; W^{1,\infty}(\Omega))$, from the estimate above and \eqref{est.Srhotau} we conclude via Aubin-Lions Lemma that (up to subsequences)
$S^{(\tau)}\rho^{(\tau)}_i\Phi$ is strongly convergent in $L^1(Q_\mathcal{T})$, for $i=1,\ldots,N$. Given the uniform $L^\infty$ bounds for $S^{(\tau)}$
and \eqref{est.rhotau}, the convergence holds in $L^{z}(Q_\mathcal{T})$
for every $z<\frac{5}{3}K_1$. Furthermore, since $\Phi$ does not depend on $\tau$ and is uniformly positive, while $S^{(\tau)}$ is strongly convergent in $L^q(Q_\mathcal{T})$ for every $q<\infty$ and satisfies \eqref{Sceps}, we conclude
\begin{equation}
\label{cnv.rhotau}
\rho^{(\tau)}_i\to\rho_i\qquad\mbox{strongly in }L^z(Q_\mathcal{T}),\quad\forall z<\frac{5}{3}K_1,\quad i=1,\ldots,N.
\end{equation}
We wish to point out that $\rho_i>0$ a.e.~in $Q_\mathcal{T}$ for $i=1,\ldots,N$ thanks to Fatou's Lemma and \eqref{est.nablamutau}.

Now we show that the temperature $T^{(\tau)}$ is strongly convergent.
We begin by estimating the energy flux $\mathcal J_E^{(\tau)}$ given by
\begin{align*}
\mathcal J_E^{(\tau)} = ((\rho e)^{(\tau)} + p^{(\tau)})\veb v^{(\tau)}
+\veb q^{(\tau)} 
-\eps(1+T^{(\tau)})\nabla\logT^{(\tau)}
-\eps \left|\nabla\frac{1}{T^{(\tau)}}\right|^{K_3-1}
\nabla\frac{1}{T^{(\tau)}} .
\end{align*}
Given the definitions of $p^{(\tau)}$, $(\rho e)^{(\tau)}$, $\veb v^{(\tau)}$, it holds
\begin{align*}
|((\rho e)^{(\tau)} + p^{(\tau)})\veb v^{(\tau)}|\leq
C(p^{(\tau)} + 2p_{at})|\nablap^{(\tau)}|
\end{align*}
From \eqref{est.naptau.K2}, \eqref{est.ptau.K2} it follows
\begin{align*}
\|p^{(\tau)}|\nablap^{(\tau)}|\|_{L^{\frac{2}{1+\theta}}(0,\mathcal{T};
	L^{\xi_2}(\Omega))}\leq C(\eps),\\
\frac{1}{\xi_2} = \frac{1}{\xi} + \frac{K_2+1}{2K_2} = 
1-\theta + \theta\frac{K_2+3}{6K_2} + \frac{K_2+1}{2K_2}.
\end{align*}
Since $\frac{3(K_2+1)}{5K_2-3}<\theta<1$, it follows
that $\xi_2>1$, $2/(1+\theta)>1$. Therefore 
\begin{align*}
\exists q>1:\quad \|p^{(\tau)}|\nablap^{(\tau)}|\|_{L^q(Q_\mathcal{T})}\leq C(\eps).
\end{align*}
As a consequence
\begin{align}
\label{est.pv.tau}
\exists q>1:\quad 
\|((\rho e)^{(\tau)} + p^{(\tau)})\veb v^{(\tau)}\|_{L^q(Q_\mathcal{T})}\leq C(\eps).
\end{align}
Let us now consider the heat flux
\begin{align*}
\veb q^{(\tau)} = -\kappa(T^{(\tau)})\nablaT^{(\tau)} 
+ \sum_{j=1}^N L_{0j}^{(\tau)}\nabla\frac{\mu_j^{(\tau)}}{T^{(\tau)}}.
\end{align*}
It is clear that 
\begin{align*}
\left\|\sum_{j=1}^N L_{0j}^{(\tau)}\nabla\frac{\mu_j^{(\tau)}}{T^{(\tau)}}\right\|_{L^2(Q_\mathcal{T})}\leq C.
\end{align*}
On the other hand, from \eqref{kappa} one finds
\begin{align*}
|\kappa(T^{(\tau)})\nablaT^{(\tau)}| 
&\leq C(1 + (T^{(\tau)})^\beta)|\nablaT^{(\tau)}|\\
&\leq CT^{(\tau)}|\nabla\logT^{(\tau)}|+
C (T^{(\tau)})^{1+\beta/2}|\nabla (T^{(\tau)})^{\beta/2}|.
\end{align*}
From \eqref{hp.K123}, \eqref{est.Ttau} it follows that
$T^{(\tau)}$, $(T^{(\tau)})^{1+\beta/2}$ are bounded in $L^{2+\omega}(Q_\mathcal{T})$ for some $\omega>0$ uniformely in $\tau$.
Therefore from \eqref{Est.A} one finds
\begin{align}
\label{est.kappa.tau}
\exists q>1:\quad \|\kappa(T^{(\tau)})\nablaT^{(\tau)}\|_{L^q(Q_\mathcal{T})}\leq C(\eps).
\end{align}
From \eqref{est.Ttau.2} it follows immediately
\begin{equation}
\label{est.AAA}
\eps\left\| \left|\nabla\frac{1}{T^{(\tau)}}\right|^{K_3-1}
\nabla\frac{1}{T^{(\tau)}}
\right\|_{L^{\frac{K_3+1}{K_3}}(Q_\mathcal{T})}\leq C\eps^{\frac{1}{K_3+1}}.
\end{equation}
Since the exponents $q$ in \eqref{est.pv.tau}, \eqref{est.kappa.tau} do not depend on $K_3$,
we conclude that, if $K_3$ is large enough, the energy flux $\mathcal J_E^{(\tau)}$ can be estimated as follows:
\begin{equation}
\label{est.JE.tau}
\|\mathcal J_E^{(\tau)}\|_{L^{\frac{K_3+1}{K_3}}(Q_\mathcal{T})}\leq C(\eps).
\end{equation}
An estimate for the total energy $\mathcal E^{(\tau)} = \Phi E_{f,\eps}^{(\tau)} + (1-\Phi)E_s^{(\tau)}$ is now required.
Since
\begin{align}\label{def.E}
\mathcal E^{(\tau)} &= 
\Phi\left( 
S^{(\tau)}( (\rho^{(\tau)})^\gamma + c_w\rho^{(\tau)}T^{(\tau)} + p_{at}
+ \eps (\rho^{(\tau)})^{K_1} )
- \int_{1/2}^{S^{(\tau)}}P_{c,\eps}(\xi)d\xi
\right)\\
&+ (1-\Phi)\left( 
c_sT^{(\tau)} + \eps(K_2-1)(T^{(\tau)})^{K_2}
\right),\nonumber
\end{align}
if follows from \eqref{Seps.bounds} (also thanks to \eqref{gamma} and Young's inequality)
\begin{align*}
|\mathcal E^{(\tau)}|\leq C(\eps)\left(1 + (\rho^{(\tau)})^{K_1} 
+ (T^{(\tau)})^{K_2}\right).
\end{align*}
From \eqref{hp.K123}, \eqref{est.rhotau}, \eqref{est.Ttau}, 
it follows that 
\begin{align*}
\exists q>1:\quad \|(\rho^{(\tau)})^{K_1} + (T^{(\tau)})^{K_2}\|_{L^q(Q_\mathcal{T})}\leq C(\eps).
\end{align*}
It follows that $\mathcal E^{(\tau)}$ is bounded in 
$L^q(Q_\mathcal{T})$ uniformely in $\tau$ for some $q>1$ which does not depend on $K_3$. Therefore choosing $K_3$ large enough yields
\begin{equation}
\label{est.E.tau}
\|\mathcal E^{(\tau)}\|_{L^{\frac{K_3+1}{K_3}}(Q_\mathcal{T})}\leq C(\eps).
\end{equation}
We aim to prove strong convergence of the temperature via the Div-Curl Lemma. For this reason we find it more convenient to work with time-continuous functions. With this idea in mind, we define now a piecewise linear interpolation operator for the time variable. Given any $u : [0,\mathcal{T}]\to X$ (where $X$ is a generic Banach space), define
\begin{align*}
\mathscr L_\tau[u](t) &= u(0)\chf{ \{0\} }(t) + 
\sum_{k=1}^N\left(
u(t_{k-1})\frac{t_k - t}{\tau} + 
u(t_k)\frac{t-t_{k-1}}{\tau}
\right)\chf{(t_{k-1},t_k] }(t),\\
t_k &= k\tau,\quad k=0,\ldots,N,\quad N = T/\tau.
\end{align*}
We point out that, for every function $u : [0,\mathcal{T}]\to X$,
the following identity holds:
\begin{align}\label{Ldt}
\frac{d}{dt}\mathscr{L}_\tau[u] = D_\tau u^{(\tau)},
\end{align}
where $u^{(\tau)}$ is the piecewise constant approximation of $u$ and $D_\tau$ is the usual discrete time derivative operator.

Let us now define the vector fields
\begin{align*}
\mathcal U^{(\tau)} = ( \mathscr{L}_\tau[\mathcal E^{(\tau)}],
\mathcal J^{(\tau)} ),\qquad
\mathcal V^{(\tau)}[G] = (G(1/T^{(\tau)}),0,0,0),
\end{align*}
where $G\in W^{1,\infty}(\R_+)$ is arbitrary.
From \eqref{est.JE.tau}, \eqref{est.E.tau} and the definition of $\mathscr{L}_\tau$ it follows
\begin{align*}
\|\mathcal U^{(\tau)}  \|_{L^q(Q_\mathcal{T})}\leq C(\eps).
\end{align*}
On the other hand, 
eq.~\eqref{app.2b} implies the following identity:
\begin{align}\label{eq.E}
\int_0^\mathcal{T}\int_\Omega & 
D_\tau(\mathcal E^{(\tau)} )\psi dx dt - 
\int_0^\mathcal{T}\int_\Omega\mathcal J_E^{(\tau)}\cdot\nabla\psi dx dt \\
\nonumber
&= -\eps\int_0^\mathcal{T}\int_\Omega (1+(T^{(\tau)})^{-K_3})(\logT^{(\tau)})\psi dx dt\\ \nonumber
&\qquad\forall\psi\in L^{K_3+1}(0,\mathcal{T}; W_0^{1,K_3+1}(\Omega)),
\end{align}
which, thanks to \eqref{Ldt}, implies
\begin{align*}
\left\langle\dv_{(t,x)}\mathcal U^{(\tau)}, \psi\right\rangle
= -\eps\int_0^\mathcal{T}\int_\Omega (1+(T^{(\tau)})^{-K_3})(\logT^{(\tau)})\psi dx dt\quad\forall\psi\in C^1_c(Q_\mathcal{T}),
\end{align*}
which, thanks to \eqref{est.Ttau.2}, yields
\begin{align*}
\|\dv_{(t,x)}\mathcal U^{(\tau)}\|_{L^1(Q_\mathcal{T})}\leq C.
\end{align*}
In particular, $\dv_{(t,x)}\mathcal U^{(\tau)}$ is relatively compact in $W^{-1,r}(Q_\mathcal{T})$ for some $r>1$.
On the other hand, while $\|\mathcal V^{(\tau)}[G]\|_{L^\infty(Q_\mathcal{T})}\leq C[G]$, it holds
\begin{align*}
|\curl_{(t,x)}\mathcal V^{(\tau)}[G]|\leq C[G]|\nabla(1/T^{(\tau)})|,
\end{align*} 
which, together with \eqref{est.Ttau.2}, leads to
\begin{align*}
\|\curl_{(t,x)}\mathcal V^{(\tau)}[G]\|_{L^{K_3+1}(Q_\mathcal{T})}\leq C(\eps).
\end{align*}
In particular $\curl_{(t,x)}\mathcal V[G]$ is relatively compact in $W^{-1,r}(Q_\mathcal{T})$ for some $r>1$.

At this point we are in the condition of applying the Div-Curl lemma and deduce
\begin{align*}
\overline{\mathcal{U}^{(\tau)}\cdot\mathcal{V}^{(\tau)}[G]}
=\overline{\mathcal{U}^{(\tau)}}\, \cdot\,
\overline{\mathcal{V}^{(\tau)}[G]}
\end{align*}
that is
\begin{align*}
\overline{\mathscr{L}_\tau[\mathcal E^{(\tau)}]G(1/T^{(\tau)})}
=
\overline{\mathscr{L}_\tau[\mathcal E^{(\tau)}]}\,\,
\overline{G(1/T^{(\tau)})}\quad\mbox{a.e.~}Q_\mathcal{T} .
\end{align*}
We now wish to get rid of the piecewise linear interpolation operator $\mathcal L_\tau$. Let $\psi\in C^1_c(Q_\mathcal{T})$ arbitrary and define $\omega^{(\tau)} = \mathcal E^{(\tau)} - \mathscr{L}_\tau[\mathcal E^{(\tau)}]$.
From the above identity it follows
\begin{align}\label{aa.1}
&\int_{Q_\mathcal{T}} (\overline{\mathcal E^{(\tau)}G(1/T^{(\tau)})} - 
\overline{\mathcal E^{(\tau)}}\,\,
\overline{G(1/T^{(\tau)})})\psi  dx dt\\ 
\nonumber
& = 
\int_{Q_\mathcal{T}}( \overline{(\mathcal E^{(\tau)} - \mathscr{L}_\tau[\mathcal E^{(\tau)}])
	G(1/T^{(\tau)})})\psi dx dt\\ \nonumber
	& - 
\int_{Q_\mathcal{T}}(\overline{(\mathcal E^{(\tau)} - \mathscr{L}_\tau[\mathcal E^{(\tau)}])}\,\,\overline{
	G(1/T^{(\tau)})})\psi dx dt \\ \nonumber
&=\lim_{\tau,\sigma\to 0}
\frac{1}{2}\int_{Q_\mathcal{T}}(\omega^{(\tau)}-\omega^{(\sigma)})
(G(1/T^{(\tau)}) - G(1/T^{(\sigma)}))\psi dx dt \\ \nonumber
&\leq \lim_{\tau,\sigma\to 0}\frac{1}{2}
\| \omega^{(\tau)}-\omega^{(\sigma)} \|_{L^{\frac{K_3+1}{K_3}}(0,\mathcal{T}; W^{-1,\frac{K_3+1}{K_3}}(\Omega))}\\
\nonumber
&\qquad\times
\| (G(1/T^{(\tau)}) - G(1/T^{(\sigma)}))\psi\|_{
	L^{K_3+1}(0,\mathcal{T}; W^{1,K_3+1}(\Omega))}\\
\nonumber
&\leq C\lim_{\tau\to 0}\|\omega^{(\tau)}\|_{L^{\frac{K_3+1}{K_3}}(0,\mathcal{T}; W^{-1,\frac{K_3+1}{K_3}}(\Omega))}
\|G(1/T^{(\tau)})\psi\|_{
	L^{K_3+1}(0,\mathcal{T}; W^{1,K_3+1}(\Omega))}.
\end{align}
Since $G\in W^{1,\infty}(\R_+)$ and \eqref{est.Ttau.2} holds,
it follows
\begin{align}\label{aa.2}
\|G(1/T^{(\tau)})\psi\|_{
	L^{K_3+1}(0,\mathcal{T}; W^{1,K_3+1}(\Omega))}
\leq C(\eps)\|\psi\|_{
	L^{K_3+1}(0,\mathcal{T}; W^{1,K_3+1}(\Omega))
\cap L^\infty(Q_\mathcal{T})}.
\end{align}
On the other hand, since
\begin{align*}
\omega^{(\tau)}(t) &= \mathcal E^{(\tau)}(t) - \mathscr{L}_\tau[\mathcal E^{(\tau)}](t)\\ 
&= \sum_{k=1}^N\frac{t_k-t}{\tau}(\mathcal E^{(\tau)}(t_{k-1}) - \mathcal E^{(\tau)}(t_k))
\chf{(t_{k-1},t_k] }(t) =-\tau\sum_{k=1}^N\frac{t_k-t}{\tau}
(D_\tau\mathcal{E}^{(\tau)})(t_k)\chf{(t_{k-1},t_k] }(t),
\end{align*}
therefore
\begin{align*}
&\|\omega^{(\tau)}\|_{L^{\frac{K_3+1}{K_3}}(0,\mathcal{T}; W^{-1,\frac{K_3+1}{K_3}}(\Omega))}\\
&\qquad\leq
\tau\sum_{k=1}^N\left\| \frac{t_k - \cdot }{\tau}
\chf{(t_{k-1},t_k] } \right\|_{L^{\frac{K_3+1}{K_3}}(0,\mathcal{T})}
\|D_\tau\mathcal{E}^{(\tau)}(t_k)\|_{
	W^{-1,\frac{K_3+1}{K_3}}(\Omega)}\\
&\qquad\leq
\tau^{1 + \frac{K_3}{K_3+1} }\sum_{k=1}^N
\|D_\tau\mathcal{E}^{(\tau)}(t_k)\|_{
	W^{-1,\frac{K_3+1}{K_3}}(\Omega)} = \tau^{\frac{K_3}{K_3+1} }
\|D_\tau\mathcal{E}^{(\tau)}\|_{L^1(0,\mathcal{T};
	W^{-1,\frac{K_3+1}{K_3}}(\Omega))}.
\end{align*}
From \eqref{est.Ttau.2}, \eqref{est.JE.tau}, \eqref{eq.E} we deduce
\begin{align*}
\|\omega^{(\tau)}\|_{L^{\frac{K_3+1}{K_3}}(0,\mathcal{T}; W^{-1,\frac{K_3+1}{K_3}}(\Omega))}\leq C \tau^{\frac{K_3}{K_3+1} },
\end{align*}
which, together with \eqref{aa.1}, \eqref{aa.2} leads to
\begin{align*}
\int_{Q_\mathcal{T}} (\overline{\mathcal E^{(\tau)}G(1/T^{(\tau)})} - 
\overline{\mathcal E^{(\tau)}}\,\,
\overline{G(1/T^{(\tau)})})\psi  dx dt\leq 0\qquad\forall\psi\in C^1_c(Q_\mathcal{T}),
\end{align*}
which easily implies
\begin{align}
\label{aa.3}
\overline{\mathcal E^{(\tau)}G(1/T^{(\tau)})} - 
\overline{\mathcal E^{(\tau)}}\,\,\overline{G(1/T^{(\tau)})} = 0
\qquad\mbox{a.e.~in }Q_\mathcal{T} .
\end{align}
From the \eqref{def.E}
and the strong convergence of $S^{(\tau)}$, $\rho^{(\tau)}$ it follows
\begin{align}\label{aa.4}
&(c_w\Phi\rho + c_s(1-\Phi))\left(
\overline{T^{(\tau)} G(1/T^{(\tau)})} 
- \overline{T^{(\tau)}}\,\,\overline{G(1/T^{(\tau)})}
\right)\\ \nonumber
&+ \eps(K_2-1)(1-\Phi)\left(
\overline{ (T^{(\tau)})^{K_2}G(1/T^{(\tau)}) } - 
\overline{ (T^{(\tau)})^{K_2} }\,\,\overline{G(1/T^{(\tau)})}
\right) = 0\\ \nonumber
&\mbox{a.e.~in }Q_\mathcal{T} .
\end{align}
Choose now $G(s) = \arctan(s)$, $s\geq 0$. The resulting function
$s\in\R_+\mapsto G(1/s)\in\R_+$ is strictly decreasing and strictly convex. From \cite[Thr.~10.19]{FeirNov09} it follows
\begin{align*}
\overline{T^{(\tau)} G(1/T^{(\tau)})} 
- \overline{T^{(\tau)}}\,\,\overline{G(1/T^{(\tau)})}\leq 0\quad
\mbox{a.e.~in }Q_\mathcal{T},\\
\overline{ (T^{(\tau)})^{K_2}G(1/T^{(\tau)}) } - 
\overline{ (T^{(\tau)})^{K_2} }\,\,\overline{G(1/T^{(\tau)})}\leq 0
\quad\mbox{a.e.~in }Q_\mathcal{T},
\end{align*}
which, together with \eqref{aa.4}, yields
\begin{align*}
\overline{T^{(\tau)} G(1/T^{(\tau)})} 
- \overline{T^{(\tau)}}\,\,\overline{G(1/T^{(\tau)})} =
\overline{ (T^{(\tau)})^{K_2}G(1/T^{(\tau)}) } - 
\overline{ (T^{(\tau)})^{K_2} }\,\,\overline{G(1/T^{(\tau)})} = 0\\
\mbox{a.e.~in }Q_\mathcal{T} .
\end{align*}
From the above identity and from \cite[Thr.~10.20]{FeirNov09}, 
we conclude that $T^{(\tau)}$ is (up to subsequences) a.e.~convergent in $Q_\mathcal{T}$. The bounds \eqref{est.Ttau}, \eqref{est.Ttau.2} allow us
to deduce by dominated convergence that 
\begin{align*}
T^{(\tau)}\to T\quad\mbox{strongly in }L^{q_1}(Q_\mathcal{T})\quad\forall q_1<\beta + \frac{2}{3}K_2,\\
\frac{1}{T^{(\tau)}}\to\frac{1}{T}\quad\mbox{strongly in }L^{q_2}(Q_\mathcal{T})\quad\forall q_2<K_3+1,\\
T>0\quad\mbox{a.e.~in }Q_\mathcal{T} .
\end{align*}
Since we have already estimated the terms appearing in \eqref{app.1b}--\eqref{app.3b}, we can now take the limit $\tau\to 0$
in \eqref{app.1b}--\eqref{app.3b} and obtain
\begin{align}
\label{app.1c}
&\int_0^\mathcal{T}\langle \pa_t(S\rho_i) , 
\Phi\varphi\rangle dt
-\int_0^\mathcal{T}\int_\Omega( \rho_i\veb v + \veb J_i )\cdot\nabla\varphi dx dt\\ \nonumber
&\qquad
-\int_0^\mathcal{T}\int_\Omega r_{i,\eps}\varphi dx dt
+\int_0^\mathcal{T}\int_{\pa\Omega}\sum_{\ell=1}^N b_{i\ell}\left( \frac{\mu_\ell}{T}
- \frac{\mu_{0,\ell}}{T_0} \right)\varphi ds dt\\ \nonumber
&\quad = -\eps\int_0^\mathcal{T}\int_\Omega\left(
\nabla\frac{\mu_i}{T}\cdot\nabla\varphi + 
\frac{\mu_i}{T}\varphi\right)dx dt,
\end{align}
\begin{align}
\label{app.2c}
&\int_0^\mathcal{T}\langle \pa_t( \Phi E_{f,\eps} + 
(1-\Phi)E_s )\psi\rangle dt
-\int_0^\mathcal{T}\int_\Omega( (\rho e + p)\veb v + \veb q )\cdot\nabla\psi dx dt\\
\nonumber
&\quad +\alpha\int_0^\mathcal{T}\int_{\pa\Omega}(T - T_0)\psi ds dt\\ \nonumber
& = -\eps\int_0^\mathcal{T}\int_\Omega((1+T)\nabla\log T
\cdot\nabla\psi + (1+T^{-K_3})(\log T)\psi)dx dt\\
\nonumber
&\quad -\eps\int_0^\mathcal{T}\int_\Omega T^{-K_3} 
|\nabla\log T|^{K_3-1}\nabla\log T\cdot\nabla\psi dxdt,\\ 
\label{app.3c}
&\pa_t f(S) + P_{c,\eps}(S) + p = 0.
\end{align}
Let $\zeta\in C^1(\overline{Q_\mathcal{T}})$ arbitrary. Let us choose 
$\phi = \frac{\mu_i}{T}\zeta$, $\psi = -\frac{1}{T}\zeta$ in \eqref{app.1c},
\eqref{app.2c}, respectively. We point out that
thanks to \eqref{est.Ttau.2}, \eqref{est.JE.tau}, the choice
$\psi = -\frac{1}{T}\zeta$ as test function in \eqref{app.2c} is admissible. 

By proceeding like in the proof of Lemma \ref{lem.EBE}
we obtain the approximate entropy balance equation:
\begin{align}\label{app.EBE}
&\int_0^\mathcal{T}\langle\partial_t\left[\Phi S(\rho\eta) + (1-\Phi)(\rho \eta)_s \right],\zeta\rangle dt \\ \nonumber
&\qquad -\int_0^\mathcal{T}\int_\Omega\Big( (\rho \eta) \veb{v} - \sum_{i=1}^N \frac{\mu_i}{T} \veb{J}_i + \frac{\veb{q}}{T}
-\eps\sum_{i=1}^N\frac{\mu_i}{T}\nabla\frac{\mu_i}{T}\\
&\qquad\qquad -\eps\frac{1+T}{T}\nabla\log T
-\eps T^{-K_3-1}|\nabla\log T|^{K_3-1}\nabla\log T
\Big)\cdot\nabla\zeta dx dt  	\nonumber\\
& =\int_0^\mathcal{T}\int_\Omega\left( 
\sum_{i,j=1}^N L_{ij} \nabla \big( \frac{\mu_i}{T} \big) \cdot \nabla \big( \frac{\mu_j}{T} \big) 
+ \kappa(T)|\nabla\log T|^2\right)\zeta dx dt \nonumber\\ 
&\qquad + \int_0^\mathcal{T}\left( 
\frac{K}{T}\lambda(S,T)|\nabla p|^2
-\Phi \frac{1}{T} f'(S) \big( \partial_t S \big)^2
-\sum_{i=1}^N r_{i,\eps}\frac{\mu_i}{T}
\right)  \zeta dx dt\nonumber	\\
&\qquad +\int_0^{\mathcal{T}}\int_{\pa\Omega}\left(\alpha\frac{T_0 - T}{T}
+\sum_{i,j=1}^N b_{ij}\frac{\mu_i}{T}\left(\frac{\mu_j}{T} - \frac{\mu_{0,j}}{T_0}\right)
\right)\zeta d\sigma dt  \nonumber\\ \nonumber
&\qquad + \eps\int_0^\mathcal{T}\int_\Omega\Bigg( 
\left|\nabla\frac{\vec\mu}{T}\right|^2 +
\left|\frac{\vec\mu}{T}\right|^2 +
\frac{1+T}{T}|\nabla\log T|^2 + 
\left|\nabla\frac{1}{T}\right|^{K_3+1} \\ \nonumber
& \qquad\qquad +
(1+T^{-K_3})\frac{1}{T}\log\frac{1}{T}
\Bigg)\zeta  dx dt,\qquad
\forall\zeta\in C^1(\overline{Q_\mathcal{T}}).
\nonumber
\end{align}
Choosing $\psi=\psi(t')=\chf{[0,t]}(t')$ in \eqref{app.2c} yields the approximate integrated total energy balance:
\begin{align}
\label{app.totE}
&\int_\Omega (\Phi E_{f,\eps}(t) + (1-\Phi)E_s(t)) dx + 
\alpha\int_0^t\int_\Omega (T-T_0)dx dt' \\ \nonumber
&\quad = \int_\Omega (\Phi E_{f,\eps}^{in} + (1-\Phi)E_s^{in}) dx
-\eps\int_0^t\int_\Omega (1+T^{-K_3})\log T dx dt' ,\\
&\qquad\qquad t\in [0,\mathcal{T}].\nonumber
\end{align}

\subsection{Limit $\eps\to 0$.}
This limit is essentially the content of the weak stability analysis. The final equations \eqref{PC.Sw}, \eqref{weak.1}--\eqref{weak.EBE} will result from taking the limit $\eps\to 0$ in \eqref{app.1c}, \eqref{app.3c}, \eqref{app.EBE}, \eqref{app.totE}.
We only need to check the following two things: (i) that the regularizing terms in the equations vanish as $\eps\to 0$, and
(ii) that the regularizing terms do not destroy the argument proving the strong convergence of the approximate solution.
\medskip\\
{\em Step 1: error terms vanish.}
Let us begin by showing that the regularizing terms in \eqref{app.1c}, \eqref{app.3c}, \eqref{app.EBE}, \eqref{app.totE} tend to zero as $\eps\to 0$ (with the exception of the fourth integral on the right-hand side of \eqref{app.EBE}, which will converge towards a nonnegative Radon measure that will be included in the term $\xi$ appearing in \eqref{weak.EBE}).

Clearly \eqref{est.nablamutau} implies
\begin{align*}
-\eps\int_0^\mathcal{T}\int_\Omega\left(
\nabla\frac{\mu_i}{T}\cdot\nabla\varphi + 
\frac{\mu_i}{T}\varphi\right)dx dt\to 0\qquad\forall\varphi\in L^2(0,\mathcal{T}; H^1(\Omega)).
\end{align*}
Next, let us consider the regularizing term in \eqref{app.3c}.
This term is only present if $k_p=0$ in \eqref{hp.Pc.bound} (see the definition \eqref{def.Pceps} of $P_{c,\eps}$). 
Let $s_0$ like in \eqref{hp.Pcf2}. Since we are assuming $k_p=0$, \eqref{hp.Pc.bound} implies that $s_1\in (0,s_0)$
exists such that $\inf_{(0,s_1)}P_c>p_{at}$.
Let $\zeta : [0,1]\to [0,1]$ be a cutoff such that
$\zeta = 1$ in $[0,s_1/2]$, $\zeta = 0$ in $[s_1,1]$,
$\zeta$ is nonincreasing.
Let us test \eqref{app.3c} against $-\zeta(S^{(\eps)})\log S^{(\eps)}$ 
and recall that we consider a truncated initial datum for the saturation (given by \eqref{Sin.eps}) in place of $S^{in}$:
\begin{align*}
\int_\Omega\tilde f(S^{(\eps)}(\mathcal{T}))dx -
\int_0^\mathcal{T}\int_\Omega(P_c(S^{(\eps)})-\eps\log S^{(\eps)} )
\zeta(S^{(\eps)})\log S^{(\eps)}dx dt'
\\ = \int_\Omega\tilde f(S^{in,\eps})dx + 
\int_0^\mathcal{T}\int_\Omega p\zeta(S^{(\eps)})\log S^{(\eps)}dx dt'
\\ \leq \int_\Omega\tilde f(S^{in,\eps})dx - 
p_{at}\int_0^\mathcal{T}\int_\Omega 
\zeta(S^{(\eps)})\log S^{(\eps)}dx dt',
\end{align*}
where $\tilde f(s) = \int_{s_1}^s f'(u)\zeta(u)(-\log u)du$
and $s_0$ is like in \eqref{hp.Pcf2}.
Since $\tilde f\geq 0$, $\inf_{(0,s_1)}P_c>p_{at}$
and \eqref{lb.dfeps} holds, it follows
\begin{align*}
\eps\int_0^\mathcal{T}\int_\Omega (\log S^{(\eps)})^2 \zeta(S^{(\eps)})dx dt\leq
 \int_\Omega\tilde f(S^{in,\eps})dx\leq C.
\end{align*}
We deduce
\begin{align}
	\label{est.reg.S}
	\sqrt{\eps}\|\log S^{(\eps)}\|_{L^2(Q_\mathcal{T})}\leq C.
\end{align}
This implies that 
\begin{equation}
\label{reg.S.lim}
P_{c,\eps}(S^{(\eps)})-P_{c}(S^{(\eps)})\to 0
\quad\mbox{strongly in }L^2(Q_\mathcal{T}),
\end{equation}
as well as
\begin{equation}
\label{reg.Ef.lim}
E_{f,\eps} - E_f = \eps\int_{1/2}^{S^{(\eps)}} \log(\xi) d\xi\to 0\quad\mbox{strongly in }L^\infty(Q_\mathcal{T}).
\end{equation}
Let us now consider the error term in the entropy flux:
\begin{align}\label{R.eps}
R^{(\eps)} = -\eps\sum_{i=1}^N\frac{\mu_i^{(\eps)}}{T^{(\eps)}}\nabla\frac{\mu_i^{(\eps)}}{T^{(\eps)}} -\eps\frac{1+T^{(\eps)}}{T^{(\eps)}}\nabla\log T^{(\eps)}\\
\nonumber
-\eps (T^{(\eps)})^{-K_3-1}|\nabla\log T^{(\eps)}|^{K_3-1}\nabla\log T^{(\eps)} .
\end{align}
From \eqref{est.nablamutau} it follows
\begin{align*}
\eps\sum_{i=1}^N\frac{\mu_i^{(\eps)}}{T^{(\eps)}}\nabla\frac{\mu_i^{(\eps)}}{T^{(\eps)}}\to 0 \qquad\mbox{strongly in }
L^{ \frac{2a}{a+2} }(Q_\mathcal{T}).
\end{align*}
From \eqref{Est.A}, \eqref{est.Ttau.2} it follows immediately 
that 
\begin{align*}
\eps\frac{1+T^{(\eps)}}{T^{(\eps)}}\nabla\log T^{(\eps)}\to 0
\qquad\mbox{strongly in }L^2(Q_\mathcal{T}).
\end{align*}
H\" older inequality yields
\begin{align*}
&\eps\int_{Q_\mathcal{T}} (T^{(\eps)})^{-K_3-1}|\nabla\log T^{(\eps)}|^{K_3}dx dt =
\eps\int_{Q_\mathcal{T}} (T^{(\eps)})^{-1}|\nabla (T^{(\eps)})^{-1}|^{K_3}dx dt\\
&\leq C\eps
\left(\int_{Q_\mathcal{T}}(T^{(\eps)})^{-K_3-1}dx dt \right)^{\frac{1}{K_3+1}}
\left( \int_{Q_\mathcal{T}} |\nabla (T^{(\eps)})^{-1}|^{K_3+1}
dx dt \right)^{ \frac{K_3}{K_3+1} } .
\end{align*}
From \eqref{est.Ttau.2} we deduce
\begin{align*}
\eps\int_{Q_\mathcal{T}} (T^{(\eps)})^{-K_3-1}|\nabla\log T^{(\eps)}|^{K_3}dx dt\leq
C\left(\eps\int_{Q_\mathcal{T}}(T^{(\eps)})^{-K_3-1}dx dt \right)^{\frac{1}{K_3+1}} .
\end{align*}
On the other hand,
\begin{align*}
&\eps\int_{Q_\mathcal{T}}(T^{(\eps)})^{-K_3-1}dx dt\\
& = 
\eps\int_{Q_\mathcal{T}\cap\{ T^{(\tau)} < \eps^{\frac{1}{2(K_3+1)}} \}  }(T^{(\eps)})^{-K_3-1}dx dt
+\eps\int_{Q_\mathcal{T}\cap\{ T^{(\tau)} \geq \eps^{\frac{1}{2(K_3+1)}} \}  }(T^{(\eps)})^{-K_3-1}dx dt .
\end{align*}
From \eqref{est.Ttau.2} it follows
\begin{align*}
\eps & \int_{Q_\mathcal{T}\cap\{ T^{(\tau)} < \eps^{\frac{1}{2(K_3+1)}} \}  }(T^{(\eps)})^{-K_3-1}dx dt\\
&\leq
\frac{2(K_3+1)\eps}{ -\log\eps }
\int_{Q_\mathcal{T}\cap\{ T^{(\tau)} < \eps^{\frac{1}{2(K_3+1)}} \}  }(T^{(\eps)})^{-K_3-1}\log\frac{1}{T^{(\tau)}}  dx dt \leq \frac{C}{ -\log\eps } ,
\end{align*}
while
\begin{align*}
\eps\int_{Q_\mathcal{T}\cap\{ T^{(\tau)} \geq \eps^{\frac{1}{2(K_3+1)}} \}  }(T^{(\eps)})^{-K_3-1}dx dt\leq C\sqrt{\eps},
\end{align*}
so we deduce
\begin{align}\label{err.T.eps}
\eps\int_{Q_\mathcal{T}}(T^{(\eps)})^{-K_3-1}dx dt\leq 
C\left( \sqrt{\eps} + \frac{1}{ -\log\eps } \right)
\to 0\quad\mbox{as }\eps\to 0.
\end{align}
Therefore the error term $R^{(\eps)}$ in the entropy flux tends to 0 strongly in $L^1(Q_\mathcal{T})$ as $\eps\to 0$. The error term on the right-hand side of \eqref{app.EBE}, on the other hand, thanks to the estimates of the error terms previously derived, converges weakly-* in $\mathcal{M}(\overline{Q_\mathcal{T}})$ towards an element 
$\tilde\xi\in\mathcal{M}(\overline{Q_\mathcal{T}})$ satisfying $\langle\tilde\xi,\varphi\rangle\geq 0$ for every $\varphi\in C(\overline{Q_\mathcal{T}})$ such that $\varphi\geq 0$ in $\overline{Q_\mathcal{T}}$.

In the integrated energy balance \eqref{app.totE} the error term on the right hand side easily tends to zero as $\eps\to 0$
thanks to \eqref{err.T.eps}:
\begin{align*}
\eps\int_0^t\int_\Omega (1+T^{-K_3})\log T dx dt'\to 0
\qquad\mbox{as }\eps\to 0.
\end{align*}
Let us now turn our attention to the error terms in the thermodynamic quantities defined in \eqref{mu.eps}--\eqref{Es.eps}. Such error terms will tend to zero strongly in $L^1(Q_\mathcal{T})$ as $\eps\to 0$ provided that
\begin{align*}
\eps(\rho^{(\eps)})^{K_1} + \eps (T^{(\eps)})^{K_2}\to 0\quad\mbox{strongly in }L^1(Q_\mathcal{T}).
\end{align*}
From \eqref{Est.p}, \eqref{hp.K123}, \eqref{est.rhotau} it follows by interpolation
\begin{align*}
\|\rho^{(\eps)}\|_{L^{K_1}(Q_\mathcal{T})}\leq
\|\rho^{(\eps)}\|_{L^{\gamma}(Q_\mathcal{T})}^{ \frac{\gamma}{2K_1-\gamma}}
\|\rho^{(\eps)}\|_{L^{2K_1}(Q_\mathcal{T})}^{\frac{2K_1 - 2\gamma}{2K_1-\gamma}}
\leq C\eps^{-\frac{3\alpha_r}{4K_1}\frac{2K_1 - 2\gamma}{2K_1-\gamma}},
\end{align*}
which leads to
\begin{align*}
\int_{Q_\mathcal{T}}\eps(\rho^{(\eps)})^{K_1}dx dt\leq
C \eps^{1-\frac{3\alpha_r}{4}\frac{2K_1 - 2\gamma}{2K_1-\gamma}} .
\end{align*}
From \eqref{hp.kr}, \eqref{hp.K123} it follows that
\begin{align*}
\int_{Q_\mathcal{T}}\eps(\rho^{(\eps)})^{K_1}dx dt\to 0\quad\mbox{as }\eps\to 0.
\end{align*}
On the other hand, \eqref{Est.T} yields immediately that 
$\eps (T^{(\eps)})^{K_2}\to 0$ strongly in $L^1(Q_\mathcal{T})$
in the case that $K_2\leq \beta + \frac{2}{3}$. Therefore
let us consider the case $K_2 > \beta + \frac{2}{3}$,
$K_2<3\beta$ (recall \eqref{hp.K123}).
From \eqref{Est.T}, \eqref{est.Ttau}
we deduce by interpolation
\begin{align*}
&\|T^{(\eps)}\|_{L^{K_2}(Q_\mathcal{T})}\leq
\|T^{(\eps)}\|_{L^{ \beta + \frac{2}{3} }(Q_\mathcal{T})}^{1-\theta}
\|T^{(\eps)}\|_{L^{\beta + \frac{2}{3}K_2}(Q_\mathcal{T})}^{\theta}
\leq C\eps^{ -\frac{(2/3)\theta}{\beta + (2/3)K_2}}, \\
&\qquad\theta = \frac{(3\beta + 2K_2)(K_2-\beta-2/3)}{2K_2(K_2-1)},
\end{align*}
which implies
\begin{align*}
\eps\int_{Q_\mathcal{T}}(T^{(\eps)})^{K_2}dx dt\leq
C\eps^{1 - \frac{(2/3)K_2\theta}{\beta + (2/3)K_2}}.
\end{align*}
It is easy to see that $1 - \frac{(2/3)K_2\theta}{\beta + (2/3)K_2}>0$ if and only if $\beta>1/3$, which is true by assumption \eqref{hp.beta}. This means
\begin{align}\label{eps.TK2}
	\int_{Q_\mathcal{T}}\eps(T^{(\eps)})^{K_2}dx dt\to 0\quad\mbox{as }\eps\to 0.
\end{align}
It follows that the error terms in \eqref{mu.eps}--\eqref{Es.eps} tend to zero strongly in $L^1(Q_\mathcal{T})$ as $\eps\to 0$. Step 1 is therefore complete.\medskip\\
{\em Step 2: strong convergence holds.} The strong convergence of the saturation $S^{(\eps)}$ as well as the positive lower bound for $S^{(\eps)}$ work exactly like in the weak stability in Section \ref{sec.weakstab}.

In the proof of convergence of the total densities, it holds that the only additional error term appearing in $U^{(\eps)}$ can be found in the definition of the pressure $p^{(\eps)}$, and only contributes to the velocity field $\veb v^{(\eps)}$ via Darcy's Law; however for the regularized velocity field the same bounds hold, that hold for the not-regularized velocity, therefore the additional error term appearing in $U^{(\eps)}$ does not influence its integrability, which means
that the vector field $U^{(\eps)}$ is still bounded in $L^r(Q_\mathcal{T})$ for some $r>1$. Furthermore
\begin{align*}
\dv_{(t,x)}U^{(\eps)} = \eps\Delta\sum_{i=1}^N\frac{\mu_i^{(\eps)}}{T^{(\eps)}}
-\eps\sum_{i=1}^N\frac{\mu_i^{(\eps)}}{T^{(\eps)}},
\end{align*} 
and therefore, thanks to \eqref{est.nablamutau},
$\dv_{(t,x)}U^{(\eps)}\to 0$ strongly in $L^2(0,\mathcal{T}; H^1(\Omega)')$ (and a fortiori strongly in $W^{-1,z}(Q_\mathcal{T})$ for some $z>1$). This means that the Div-Curl Lemma can still be applied to the couple $U^{(\eps)}$, $V^{(\eps)}$ and argue like in the second step of the weak stability argument to deduce the strong convergence of the total density $\rho^{(\eps)}$. 

In a similar way the proof of the strong convergence of the partial mass densities $\rho_i^{(\eps)}$, $i=1,\ldots,N$, can be adapted to the approximate system.

Concerning the strong convergence of the temperature, the vector field $W^{(\eps)}$ is still bounded in $L^r(Q_\mathcal{T})$
for some $r>1$ since the only error term in $W^{(\eps)}$ (aside from the one in the velocity field, which can be ignored as in the case of $U^{(\eps)}$) appears in the skeleton entropy $(\rho\eta)_{s,\eps}$ and is proportional to $\eps (T^{(\eps)})^{K_2-1}$, which tends to zero (a fortiori, is bounded) strongly in $L^{K_2/(K_2-1)}(Q_\mathcal{T})$ thanks to
\eqref{eps.TK2}. On the other hand,
\begin{align*}
\dv_{(t,x)}W^{(\eps)} = 
-\dv_x R^{(\eps)} + \Xi^{(\eps)},
\end{align*}
where $R^{(\eps)}$ is given by \eqref{R.eps} and $\Xi$ is the right-hand side of \eqref{app.EBE}.
We have already proved that $\Xi^{(\eps)}$ is bounded in $L^1(Q_\mathcal{T})$, while $R^{(\eps)}\to 0$ strongly in $L^1(Q_\mathcal{T})$.
We can therefore apply the Div-Curl Lemma in the version of \cite{CoDoMu11} to the couple $(W^{(\eps)} , Y^{(\eps)}[G])$
and, by arguing like in the weak stability analysis, prove the strong convergence of the temperature $T^{(\eps)}$.

This finishes the proof of the existence theorem.








\section{Appendix}
To improve the readability of the paper,
we present here some proofs and results which are rather technical, but are neverthless needed for completeness. 
\subsection*{Proof of Proposition \ref{prop.5}}
Since \eqref{vow}, \eqref{mobility} hold and $\mu$ is uniformly positive, we deduce
\begin{align*}
\rho |\veb v| &\leq C \rho \sqrt{\lambda(S,T)T}\, \sqrt{\frac{\lambda(S,T)}{T}}|\nabla p| 
\leq C (k_r(S))^{1/2}\rho \, T^{1/2}\, \sqrt{\frac{\lambda(S,T)}{T}}|\nabla p|.
\end{align*}
From \eqref{hp.kr} and H\" older's inequality
(notice that $ \frac{\gamma + 6\beta + 3\beta\gamma}{6\beta\gamma}
= \frac{1}{\gamma} + \frac{1}{6\beta} + \frac{1}{2}$) we get
\begin{align*}
	\|\rho(t)& \veb v(t) \|_{ L^{\frac{6\beta\gamma}{\gamma + 6\beta + 3\beta\gamma}}(\Omega)}^{\frac{2\beta}{\beta+1}} \\
	&\leq
	C\|S(t)^{1/\gamma} \rho(t) \|_{L^\gamma(\Omega)}^{\frac{2\beta}{\beta+1}}
	\|\sqrt T(t)\|_{L^{6\beta}(\Omega)}^{\frac{2\beta}{\beta+1}}
	\left\|
	\sqrt{\frac{\lambda(S(t),T(t))}{T(t)}}\nabla p(t)
	\right\|_{L^2(\Omega)}^{\frac{2\beta}{\beta+1}}\\
	&\leq
	C\left(\esup_{t'\in [0,\mathcal{T}]}\|S(t')^{1/\gamma} \rho(t') \|_{L^\gamma(\Omega)}^{\frac{2\beta}{\beta+1}}\right)
	\|\sqrt T(t)\|_{L^{6\beta}(\Omega)}^{\frac{2\beta}{\beta+1}}
	\left\|
	\sqrt{\frac{\lambda(S(t),T(t))}{T(t)}}\nabla p(t)
	\right\|_{L^2(\Omega)}^{\frac{2\beta}{\beta+1}},
\end{align*}
which, by integration and using H\" older's inequality, leads to
\begin{align*}
\|\rho\veb v & \|_{L^{\frac{2\beta}{\beta+1}}(0,\mathcal{T};  L^{\frac{6\beta\gamma}{\gamma + 6\beta + 3\beta\gamma}}(\Omega))}^{\frac{2\beta}{\beta+1}} \\
&\leq
C\|S^{1/\gamma} \rho \|_{L^\infty(0,\mathcal{T}; L^\gamma(\Omega))}^{\frac{2\beta}{\beta+1}}
\|\sqrt T\|_{L^{2\beta}(0,\mathcal{T}; L^{6\beta}(\Omega))}^{\frac{2\beta}{\beta+1}}
\left\|
\sqrt{\frac{\lambda}{T}}\nabla p
\right\|_{L^2(Q_\mathcal{T})}^{\frac{2\beta}{\beta+1}}\\
&=C\| S\rho^\gamma \|_{L^\infty(0,\mathcal{T}; L^1(\Omega))}^{\frac{2\beta}{\gamma(\beta+1)}}
\|T^{\frac{\beta}{2}}\|_{L^{2}(0,\mathcal{T}; L^{6}(\Omega))}^{\frac{2}{\beta+1}}
\left\|\sqrt{\frac{\lambda}{T}}\nabla p\right\|_{L^2(Q_\mathcal{T})}^{\frac{2\beta}{\beta+1}}.
\end{align*}
We point out that $\frac{2\beta}{\beta+1}>1$, $\frac{6\beta\gamma}{\gamma + 6\beta + 3\beta\gamma}>1$ thanks to the assumptions \eqref{gamma}, \eqref{hp.Pcf}, \eqref{hp.beta} and the parameters $\gamma, q, \beta$.
From the Sobolev embedding $H^1(\R^d)\hookrightarrow L^6(\R^d)$ (valid for $d\leq 3$) and \eqref{Est.A} it follows
$$\|T^{\frac{\beta}{2}}\|_{L^{2}(0,\mathcal{T}; L^{6}(\Omega))}\leq C. $$
From the previous estimates as well as \eqref{Est.E.LinfL1}, 
\eqref{AJPP.1} we obtain \eqref{Est.rhov}.

Let us now estimate $\veb J_i$. From \eqref{flux.J} it follows
\begin{align*}
|\veb J_i|\leq C(|\nabla\log T| + |\nabla\Pi^N(\vec {\mu}/T)|),\quad i=1,\ldots,N,
\end{align*}
so from \eqref{Est.A}, \eqref{Est.B.2} we deduce \eqref{Est.J}.

From \eqref{cons_w}, \eqref{Est.rhov}, \eqref{Est.J} we easily conclude
that \eqref{Est.flux.1} and \eqref{Est.Srhot} hold.

An estimate \eqref{Est.St} for $\pa_t f(S)$ is immediately found from \eqref{PC.Sw} and
\eqref{Est.p}.

On the other hand, from the fact that $p\leq C(\rho e)$ and the bound \eqref{Est.E.LinfL1} we deduce that $S p$ is bounded in $L^\infty(0,\mathcal{T}; L^1(\Omega))$. As a consequence, multiplying \eqref{PC.Sw} by $S$ and taking the $L^\infty(0,\mathcal{T}; L^1(\Omega))$
norm leads to \eqref{Est.FSt}. 
\subsection*{Proof of Proposition \ref{prop.6}}
From \eqref{entropy} it follows
\begin{align*}
\|S(\rho\eta)\|_{L^{\frac{2\gamma}{\gamma+2}}(Q_\mathcal{T})} \leq 
\sum_{i=1}^N\|S\rho_i\log\rho_i\|_{L^{\frac{2\gamma}{\gamma+2}}(Q_\mathcal{T})}
+ c_w\|S\rho (\log T + 1)\|_{L^{\frac{2\gamma}{\gamma+2}}(Q_\mathcal{T})}.
\end{align*} 
However, for $\delta\in (0,\gamma/2]$, since $S\leq S^{(1+\delta)/\gamma}$ and
$\rho_i\log\rho_i\leq C(1+\rho^{1+\delta})$, it holds
\begin{align*}
\sum_{i=1}^N\|S\rho_i\log\rho_i\|_{L^{\frac{2\gamma}{\gamma+2}}(Q_\mathcal{T})}
&\leq C + C\|S^\frac{1+\delta}{\gamma}\rho^{1+\delta}\|_{L^{\frac{2\gamma}{\gamma+2}}(Q_\mathcal{T})}\\
& = C + C \|S^{1/\gamma}\rho\|_{L^{\frac{2\gamma(1+\delta)}{\gamma+2}}(Q_\mathcal{T})}^{
	1+\delta}.
\end{align*}
Since $\gamma>2$ and $0<\delta\leq \gamma/2$, from \eqref{Est.E.LinfL1} it follows
\begin{align*}
	\sum_{i=1}^N\|S\rho_i\log\rho_i\|_{L^{\frac{2\gamma}{\gamma+2}}(Q_\mathcal{T})}
	&\leq C .
\end{align*}
Moreover, H\" older's inequality yields
\begin{align*}
\|S\rho (\log T + 1)\|_{L^{\frac{2\gamma}{\gamma+2}}(Q_\mathcal{T})}
\leq \|S\rho\|_{L^\gamma(Q_\mathcal{T})}\|\log T + 1\|_{L^2(Q_\mathcal{T})}\leq C
\end{align*}
thanks to \eqref{Est.E.LinfL1}, \eqref{Est.A}. We conclude
\begin{equation}
\label{Est.rhoeta}
\|S(\rho\eta)\|_{L^{\frac{2\gamma}{\gamma+2}}(Q_\mathcal{T})} \leq C.
\end{equation}
From \eqref{Skel.Entr}, \eqref{Est.A} it follows that $(\rho\eta)_s$ is bounded in $L^2(Q_\mathcal{T})$, so we get \eqref{Est.entr}.

We will now find an estimate for the entropy flux.
We begin by considering 
\begin{align*}
(\rho\eta)\veb v = 
K\left(\sum_{i=1}^N\rho_i\log\rho_i - c_w\rho (\log T + 1)\right)
\lambda(S,T) \nabla p ,
\end{align*}
where the above equality holds thanks to \eqref{entropy}, \eqref{vow}.
Let $s\in\R$ be such that
\begin{equation}
\label{def.s}
s>1,\qquad
\frac{1}{\gamma} + \frac{1}{2\beta} + \frac{1}{2} < \frac{1}{s} .
\end{equation}
The above definition makes sense since \eqref{hp.beta} holds.
It follows via H\" older's inequality
\begin{align*}
\|(\rho\eta)\veb v\|_{L^s(Q_\mathcal{T})} \leq 
C\sum_{i=1}^N\|\rho_i (\log\rho_i) \sqrt{\lambda(S,T)T}\|_{L^{\frac{2s}{2-s}}(Q_\mathcal{T})}
\left\|\sqrt\frac{\lambda(S,T)}{T}\nabla p\right\|_{L^2(Q_\mathcal{T})}\\
+C\|\rho (\log T + 1)\sqrt{\lambda(S,T)T}\|_{L^{\frac{2s}{2-s}}(Q_\mathcal{T})}
\left\|\sqrt\frac{\lambda(S,T)}{T}\nabla p\right\|_{L^2(Q_\mathcal{T})} .
\end{align*}
However, since $\mu$ is uniformly positive and \eqref{AJPP.1} holds,
we get
\begin{align*}
\|(\rho\eta)\veb v\|_{L^s(Q_\mathcal{T})} \leq 
	C\sum_{i=1}^N\|\rho_i (\log\rho_i) \sqrt{k_r(S)T}\|_{L^{\frac{2s}{2-s}}(Q_\mathcal{T})}\\
	+C\|\rho (\log T + 1)\sqrt{k_r(S)T}\|_{L^{\frac{2s}{2-s}}(Q_\mathcal{T})} .
\end{align*}
Now, since for every $\delta>0$ there exists $C_\delta>0$ such that
$x|\log x| \leq C_\delta(1+x^{1+\delta})$ for $x>0$, H\" older's inequality and \eqref{def.s} allow us to state
\begin{align*}
\|(\rho\eta)\veb v\|_{L^s(Q_\mathcal{T})} \leq 
C_\delta\sum_{i=1}^N\|(1+\rho_i^{1+\delta}) \sqrt{k_r(S)} \, T^{1/2}\|_{L^{\frac{2s}{2-s}}(Q_\mathcal{T})}\\
+C_\delta\|\rho\sqrt{k_r(S)} (1+T^{(1+\delta)/2}) \|_{L^{\frac{2s}{2-s}}(Q_\mathcal{T})}\\
\leq C_\delta\sum_{i=1}^N\| (1+\rho_i^{1+\delta}) \sqrt{k_r(S)} \|_{L^{\gamma/(1+\delta)}(Q_\mathcal{T})}
\|T^{1/2}\|_{L^{2\beta}(Q_\mathcal{T})} \\
+C_\delta\| \rho\sqrt{k_r(S)}\|_{L^\gamma(Q_\mathcal{T})}
\| (1+T^{(1+\delta)/2}) \|_{L^{2\beta/(1+\delta)}(Q_\mathcal{T})},
\end{align*}
for some $\delta>0$ small enough. Assumption \eqref{hp.kr} and bounds \eqref{Est.E.LinfL1}, \eqref{Est.A} allow us to conclude
\begin{equation}
\label{Est.rhoetav}
\exists s>1:\quad
\|(\rho\eta)\veb v\|_{L^s(Q_\mathcal{T})} \leq C.
\end{equation}
Let us then consider
\begin{align*}
-\sum_{i=1}^N\frac{\mu_i}{T}\veb J_i = -\sum_{i=1}^N\frac{\mu_i}{T}L_{i0}\nabla\frac{1}{T}
+ \sum_{i,j=1}^N\frac{\mu_i}{T}L_{ij}\nabla\frac{\mu_j}{T},
\end{align*}
where the above equality comes from \eqref{flux.J}. Since \eqref{hp.Ltilde} holds and
$L_{ij}$ is symmetric and positive semidefinite, we obtain via Cauchy-Schwartz and \eqref{Ass.L} that
\begin{align*}
\left|-\sum_{i=1}^N\frac{\mu_i}{T}\veb J_i \right| &\leq
C\left|\Pi^N\frac{\vec\mu}{T}\right| |\nabla\log T|
+ C\left( \sum_{i,j=1}^N L_{ij}\frac{\mu_i\mu_j}{T^2} \right)^{\frac{1}{2}}
\left( \sum_{i,j=1}^N L_{ij}\frac{\nabla\mu_i\cdot\nabla\mu_j}{T^2} \right)^{\frac{1}{2}}\\
&\leq C\left|\Pi^N\frac{\vec\mu}{T}\right|\left(
|\nabla\log T| + \left|\nabla\Pi^N\frac{\vec\mu}{T}\right|
\right).
\end{align*}
Given that \eqref{Est.A}, \eqref{Est.B.2} hold, we obtain
\begin{equation}
\label{Est.muJ}
\left\|-\sum_{i=1}^N\frac{\mu_i}{T}\veb J_i \right\|_{L^{\frac{2a}{a+2}}(Q_\mathcal{T})}\leq C.
\end{equation}
Finally, let us consider 
\begin{align*}
\frac{\veb q}{T} = -\kappa(T)\nabla\log T + 
\sum_{j=1}^N \frac{L_{0j}}{T}\nabla\frac{\mu_j}{T} ,
\end{align*}
where the above equality comes from \eqref{flux.q}, \eqref{heat_cond}.
It follows from \eqref{kappa},  that
\begin{align*}
\left|\kappa(T)\nabla\log T\right| \leq C|\nabla\log T| + C T^{\beta/2}|\nabla T^{\beta/2}|,
\end{align*}
while \eqref{hp.Ltilde} implies
\begin{align*}
\left|\sum_{j=1}^N \frac{L_{0j}}{T}\nabla\frac{\mu_j}{T}\right|\leq
C \left|\nabla\Pi^N\frac{\vec\mu}{T}\right|.
\end{align*}
From \eqref{Est.A}--\eqref{Est.B.2} and H\" older's inequality we conclude
\begin{align}\label{Est.kappaT}
\left\|\kappa(T)\nabla\log T\right\|_{L^{\frac{2+3\beta}{1+3\beta}}(Q_\mathcal{T})}
+
\left\|\sum_{j=1}^N \frac{L_{0j}}{T}\nabla\frac{\mu_j}{T}\right\|_{L^2(Q_\mathcal{T})}\leq C,
\end{align}
which leads to
\begin{equation}
\label{Est.qT}
\left\|\frac{\veb q}{T}\right\|_{L^{\frac{2+3\beta}{1+3\beta}}(Q_\mathcal{T})}
\leq C .
\end{equation}
Putting \eqref{Est.rhoetav}, \eqref{Est.muJ}, \eqref{Est.qT} together
allows us to obtain the estimate for the entropy flux \eqref{Est.entrflux}.
\subsection*{Proof of Proposition \ref{prop.7}}
Let us now find an estimate for $\nabla\rho$. 
Eq.~\eqref{thp.1} yields:
\begin{align*}
	\nabla p 
	= (T + \gamma(\gamma-1)\rho^{\gamma-1})\nabla\rho + \rho\nabla T ,
\end{align*}
which implies
\begin{align*}
	|\nabla\rho^\gamma|\leq |(T + \gamma(\gamma-1)\rho^{\gamma-1})\nabla\rho|\leq
	\rho |\nabla T| + |\nabla p| .
\end{align*}
From \eqref{hp.kr} and the uniform boundedness of $\mu$ it follows
\begin{align*}
	\sqrt{k_r(S)}|\nabla\rho^\gamma|\leq C\rho S^{1/\gamma}|\nabla T|
	+ T^{1/2}\sqrt\frac{\lambda(S,T)}{T} |\nabla p|.
\end{align*}
From \eqref{Est.E.LinfL1}, \eqref{Est.A}, \eqref{AJPP.1}, \eqref{Est.T}
(as well as from H\" older's inequality) it follows
\begin{align*}
\|\rho S^{1/\gamma}\nabla T\|_{L^{\frac{2\gamma}{2+\gamma}}(Q_\mathcal{T})} \leq 
\|\rho S^{1/\gamma}\|_{L^\gamma(Q_\mathcal{T})}\|\nabla T\|_{L^2(Q_\mathcal{T})}\leq C,\\
\left\|
 T^{1/2}\sqrt\frac{\lambda(S,T)}{T} \nabla p
\right\|_{L^{\frac{2\beta + 4/3}{\beta + 5/3}}(Q_\mathcal{T})}\leq
\left\|T^{1/2}\right\|_{L^{2\beta + 4/3}(Q_\mathcal{T})}\left\|
\sqrt\frac{\lambda(S,T)}{T} \nabla p
\right\|_{L^{2}(Q_\mathcal{T})}\leq C,
\end{align*}
so we deduce
\begin{align*}
\|\sqrt{k_r(S)}\nabla\rho^\gamma\|_{L^{a_2}(Q_\mathcal{T})}\leq C,\quad
a_2\equiv\min\left\{\frac{2\gamma}{2+\gamma},
\frac{2\beta + 4/3}{\beta + 5/3}  \right\} > 1.
\end{align*}
In particular we deduce
\begin{align}\label{est.narho.1}
\|\sqrt{k_r(S)}\nabla G(\rho^\gamma)]\|_{L^{a_2}(Q_\mathcal{T})}\leq C[G],
\end{align}
where $G\in W^{1,\infty}(\R_+)$ is arbitrary, and the constant $C[G]>0$ depends on $G$.

On the other hand, since $1/\gamma + 1/q\leq 1$ (see \eqref{hp.Pcf}),
from \eqref{hp.krf}, \eqref{Est.E.LinfL1}, \eqref{Est.S} it follows
\begin{align*}
\|G(\rho^\gamma)\nabla \sqrt{k_r(S)}\|_{L^{\infty}(0,\mathcal{T}; L^q(\Omega))} &\leq \| G(\rho^\gamma) \|_{L^\infty(Q_\mathcal{T})}\|\nabla f(S)\|_{L^\infty(0,\mathcal{T}; L^q(\Omega))}\leq C[G],
\end{align*}
which, together with \eqref{est.narho.1}, leads to \eqref{Est.narho}.

Let us now find another gradient estimate related to $\rho_i/\rho$, $i=1,\ldots,N$. 
From \eqref{relchp} and \eqref{Est.B.2} we get
\begin{align*}
\int_0^\mathcal{T}\int_\Omega\left( \log\frac{\rho_i}{\rho} - \frac{1}{N}\sum_{j=1}^N\log\frac{\rho_j}{\rho} \right)^2 dx dt =
\int_0^\mathcal{T}\int_\Omega\left( \log\rho_i - \frac{1}{N}\sum_{j=1}^N\log\rho_j \right)^2 dx dt\leq C,
\end{align*}
for $i=1,\ldots,N$. Lemma \ref{lem.aux.1} allows us to deduce
\begin{align*}
\left\| \sum_{j=1}^N\log\frac{\rho_j}{\rho} \right\|_{L^2(Q_\mathcal{T})}\leq C,
\end{align*}
which, together with the previous estimate, leads to \eqref{est.relchp}.
Bound \eqref{Est.r} on the reaction terms $r_1,\ldots,r_N$ comes straightforwardly from \eqref{hp.r.2}, \eqref{Est.B.2}.

\subsection*{Proof of Lemma \ref{lem.estSk}.}


We show \eqref{Sceps}, i.e.~that $S^k\geq \eps>0$ a.e.~in $\Omega$, for $k\geq 0$, provided that $0<\eps<f^{-1}(p_{at}/\lambda_0)$.
Without loss of generality we can consider \eqref{app.3} on the set $\{S^k\leq s_0\}$ where $s_0$ is defined in \eqref{hp.Pcf2}.
From the definition of $P_{c,\eps}$, the positivity of $f$ in $(0,s_0)$ and \eqref{hp.Pcf2}, it follows
$$ \inf_{s\in (0,s_0)}\frac{P_{c,\eps}(s)}{f(s)}
\geq\inf_{s\in (0,s_0)}\frac{P_{c}(s)}{f(s)}\equiv \lambda_0 
> 0. $$
The positivity of $f$ in $(0,s_0)$, the fact that $p^k\geq -p_{at}$ and \eqref{app.3} imply
\begin{align*}
\frac{f(S^k) - f(S^{k-1})}{\tau} + \lambda_0 f(S^k) - p_{at} \leq 0,
\end{align*}
and so
\begin{align*}
f(S^k) \leq (1+\tau\lambda_0)^{-1}f(S^{k-1}) + p_{at}\tau(1+\tau\lambda_0)^{-1}.
\end{align*}
We prove by induction on $k\geq 0$ that
$f(S^k)\leq f(\eps)$ for every $k\geq 0$. 
Since $S^0 = S^{in,\eps}\geq \eps$, the statement is true for $k=0$.
Assume now the statement holds true for $k-1$: $f(S^{k-1})\leq f(\eps)$. It follows
\begin{align*}
f(S^k) \leq (1+\tau\lambda_0)^{-1}f(\eps) + p_{at}\tau(1+\tau\lambda_0)^{-1}.
\end{align*}
The right-hand side of the above inequality is $\leq f(\eps)$ if $f(\eps)\geq p_{at}/\lambda_0$,
which holds true since $0<\eps<f^{-1}(p_{at}/\lambda_0)$
and $f$ is decreasing.
Therefore the statement holds true for every $k\geq 0$.
Once again, being $f$ decreasing, it follows that 
\begin{align}\label{Sceps.appendix}
S^k\geq f^{-1}(f(\eps)) = \eps > 0,
\quad k\geq 0,\quad\mbox{a.e.~in }\Omega ,
\end{align}
provided that $\eps<f^{-1}(p_{at}/\lambda_0)$.

Bound \eqref{Sceps} allows us to immediately find a bound for the spatial integral of $f(S^k)$. Indeed, integrating \eqref{app.3} in $\Omega$ yields
\begin{align*}
\int_\Omega f(S^k)dx = \int_\Omega f(S^{k-1})dx
-\tau\int_\Omega( P_{c,\eps}(S^k) + p^k )dx .
\end{align*}
However, since \eqref{Sceps} holds, then $P_{c,\eps(S^k)}\leq C(\eps)$ for every $k\geq 0$, so
\eqref{est.f.k} holds.

Finally, we will show that  \eqref{est.naf.k} holds.
Let us now take the gradient of \eqref{app.3}, test the resulting equation times $\tau\frac{4}{3}(1+|\nabla f(S^{k})|^2)^{-1/3} \nabla f(S^{k})$
and use the elementary inequality (a consequence of the convexity of $x\mapsto (1+x^2)^{2/3}$)
$$\frac{4}{3}\frac{x}{(1+x^2)^{1/3}}(x-y)\geq (1+x^2)^{2/3} - (1+y^2)^{2/3},\qquad x,y\geq 0. $$
It holds (thanks to Young's inequality)
\begin{align*}
\int_\Omega &(1+|\nabla f(S^k)|^2)^{2/3} dx + \frac{4}{3}\tau \int_\Omega \frac{f'(S^k)P'_{c,\eps}(S^k)}{
	(1+|\nabla f(S^{k})|^2)^{1/3}}|\nabla S^k|^2 dx\\
&\leq \int_\Omega (1+|\nabla f(S^{k-1})|^2)^{2/3} dx - \frac{4}{3}\tau \int_\Omega 
\frac{\nabla p^k\cdot\nabla f(S^k)}{
	(1+|\nabla f(S^{k})|^2)^{1/3}}	dx\\
&\leq \int_\Omega (1+|\nabla f(S^{k-1})|^2)^{2/3} dx
+ \tau\int_\Omega |\nabla p^k|^{4/3} dx 
+ C\tau\int_\Omega \frac{|\nabla f(S^k)|^4}{
	(1+|\nabla f(S^{k})|^2)^{4/3}}dx .
\end{align*}
Since $f' P'_{c,\eps}\geq 0$ in $(0,1)$
and $x^4(1+x^2)^{-4/3}\leq (1+x^2)^{2/3}$ for $x\geq 0$,
we conclude that \eqref{est.naf.k} holds.



\subsection{Auxiliary results}\label{sec.aux}
We begin by stating a simple but useful algebraic property.
\begin{lemma}\label{lem.aux.1}
For every $\eps>0$ there exists $C_\eps>0$ such that
\begin{align*}
\sum_{i=1}^N\left(
\log u_i - \frac{1}{N}\sum_{j=1}^N\log u_j\right)^2
\geq C_\eps\left(\sum_{i=1}^N\log u_i \right)^2 - \eps	
\end{align*}
for every $\vec{u}\in (0,\infty)^N$, such that $\sum_{i=1}^N u_i = 1$.	
\end{lemma}
\begin{proof}
By contradiction. Assume there exists $\eps_0>0$ such that, for every
$n\in\N$, there is $u^{(n)}\in (0,\infty)^N$ such that
\begin{align*}
\sum_{i=1}^N\left(
\log u_i^{(n)} - \frac{1}{N}\sum_{j=1}^N\log u_j^{(n)}\right)^2
< \frac{1}{n}\left(\sum_{i=1}^N\log u_i^{(n)} \right)^2 - \eps_0,\quad
\sum_{i=1}^N u_i^{(n)} = 1.
\end{align*}
As a consequence $\left(\sum_{i=1}^N\log u_i^{(n)} \right)^2>0$, so
we can define 
$$
v_i^{(n)} = \frac{\log u_i^{(n)}}{ \sum_{k=1}^N\log u_k^{(n)} },\quad
i=1,\ldots,N,
$$
and it follows
\begin{align}\label{aux.1}
\sum_{i=1}^N\left(v_i^{(n)} - \frac{1}{N}\right)^2
< \frac{1}{n} - \frac{\eps_0}{\left(\sum_{i=1}^N\log u_i^{(n)} \right)^2} .
\end{align}
We point out that clearly $u_i^{(n)}\leq 1$ for $i=1,\ldots,N$, which implies that $v_i^{(n)}\geq 0$ for $i=1,\ldots,N$. Furthermore $\sum_{i=1}^N v_i^{(n)} = 1$ by construction, so the sequence $v^{(n)}$ is bounded. Therefore there exists a subsequence (not relabeled) of $v^{(n)}$ that is convergent: $v^{(n)}\to v$ as $n\to\infty$.
Taking the limit $n\to\infty$ in \eqref{aux.1} yields
\begin{align*}
\lim_{n\to\infty}v_i^{(n)} = \frac{1}{N}\quad i=1,\ldots,N,
\quad
\lim_{n\to\infty}\left|\sum_{k=1}^N\log u_k^{(n)} \right| = \infty .
\end{align*}
As a consequence,
\begin{align*}
-\log u_i^{(n)} = v_i^{(n)}\left|\sum_{k=1}^N\log u_k^{(n)} \right|
\to \infty\quad\mbox{as }n\to\infty,~~ i=1,\ldots,N.
\end{align*}
This means that $u_i^{(n)}\to 0$ as $n\to\infty$ for $i=1,\ldots,N$,
which is in contradiction with the fact that $\sum_{i=1}^N u_i^{(n)}=1$.
This finishes the proof of the Lemma.	
\end{proof}


The following result is a generalized version of Fatou's Lemma.
\begin{lemma}[Generalized Fatou's Lemma]\label{lem.Fatou.gen}
Let $Q\subset\R^d$ be a bounded domain, $(u_n)_{n\in\N}$, $(v_n)_{n\in\N}$ sequences of Lebesgue-measurable functions  $Q\to\R^m$, $Q\to\R^N$ (respectively) such that 
$u_n\to u$ a.e.~in $Q$ and $v_n\rightharpoonup v$ weakly in $L^p(Q)$, for some $p\geq 1$. Let $g : \R^m\times\R^N\to\R$
be a continuous nonnegative function such that 
\begin{align*}
\mbox{for every $u\in\R^m$,
	the mapping $v\in\R^N\mapsto g(u,v)\in\R$ is convex,}\\
\exists\omega\in C^0(\R_+),~~\omega(0)=0:\quad\
|g(u,v)-g(u',v)|\leq\omega(|u-u'|)|v|^p\quad\forall u,u'\in\R^m,~~ v\in\R^N,\\
\exists C>0:\quad\int_Q g(u_n,v_n) \leq C\quad\forall n\in\N.
\end{align*}
 Then
\begin{align}\label{goal.lem.Fatou.gen}
\liminf_{n\to\infty}\int_Q g(u_n,v_n)\geq \int_Q g(u,v) .
\end{align}
\end{lemma}
\begin{proof}
Since $u_n\to u$ a.e.~in $Q$, Egorov-Severini's Theorem implies that, for every $\eps>0$, there exists $E_\eps\subset Q$ such that $|Q\backslash E_\eps|<\eps$ and $u_n\to u$ in $L^\infty(E_\eps)$. Being $g$ nonnegative  it holds
\begin{align*}
\int_Q g(u_n,v_n) \geq \int_{E_\eps}g(u_n,v_n) = 
\int_{E_\eps}(g(u_n,v_n) - g(u,v_n))
+\int_{E_\eps}g(u,v_n) .
\end{align*}
Thanks to the assumptions on $g$, $u_n$, $v_n$
\begin{align*}
\left|\int_{E_\eps}(g(u_n,v_n) - g(u,v_n))\right|\leq
\int_{E_\eps}\omega(|u_n-u|)|v_n|^p\leq
\|\omega(|u_n-u|)\|_{L^\infty(E_\eps)}\|v_n\|_{L^p(Q)}^p
\to 0
\end{align*}
as $n\to\infty$, therefore
\begin{align*}
\liminf_{n\to\infty}\int_Q g(u_n,v_n) \geq
\liminf_{n\to\infty}\int_{E_\eps} g(u,v_n).
\end{align*}
Being $v\mapsto g(u,v)$ convex, the mapping 
$v\in L^p(E_\eps)\mapsto \int_{E_\eps}g(u,v)\in\R$ is
weakly lower semicontinuous and so
\begin{align*}
	\liminf_{n\to\infty}\int_Q g(u_n,v_n) \geq
	\int_{E_\eps} g(u,v).
\end{align*}
Since $\eps>0$ is arbitrary and $|Q\backslash E_\eps|<\eps$
we conclude that \eqref{goal.lem.Fatou.gen} holds. This finishes the proof of the Lemma.
\end{proof}



