\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2013 (2013), No. 268, pp. 1--24.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2013 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2013/268\hfil Renormalized and entropy solutions] {Renormalized and entropy solutions of nonlinear parabolic systems} \author[L. Shangerganesh, K. Balachandran \hfil EJDE-2013/268\hfilneg] {Lingeshwaran Shangerganesh, Krishnan Balachandran} % in alphabetical order \address{Lingeshwaran Shangerganesh \newline Department of Mathematics, Bharathiar University, Coimbatore - 641 046, India} \email{shangerganesh@gmail.com} \address{Krishnan Balachandran \newline Department of Mathematics, Bharathiar University, Coimbatore - 641 046, India} \email{kb.maths.bu@gmail.com} \dedicatory{Dedicated to Professor Sivaguru S. Sritharan} \thanks{Submitted February 1, 2013. Published December 5, 2013.} \subjclass[2000]{35K57, 35K65, 92D30} \keywords{Renormalized solutions; entropy solutions; cross-diffusion} \begin{abstract} In this article, we study the existence of renormalized and entropy solutions of SIR epidemic disease cross-diffusion model with Dirichlet boundary conditions. Under the assumptions of no growth conditions and integrable data, we establish that the renormalized solution is also an entropy solution. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{definition}[theorem]{Definition} \newtheorem{remark}[theorem]{Remark} \allowdisplaybreaks \section{Introduction} The primary concern of this paper is the study of a certain class of nonlinear parabolic systems — more precisely, the so-called cross-diffusion systems (see \cite{benjee}), which arise in various fields of sciences and have long been the subject of extensive mathematical studies. Many time and space dependent chemical, physical or biological processes can be described with partial differential equations, respectively, by reaction-diffusion equations. A special kind of reaction-diffusion equations is given by systems of parabolic partial differential equations with cross-diffusion effects. In recent years, considerable amount of work has been done on the existence of solutions of biological parabolic systems such as epidemiological model \cite{bendaaam,benjee,bencpaa,fitana,fitsiam,fitfields,lsgm2as} and its related chemotaxis, predator-prey model \cite{bencross,benjmaa,benade,del,horst,lsgaam,lsgaa,yxwang,zeng}, see the references therein. Here one of the topics of interest is the problem of cross-diffusion, the phenomenon in which a gradient in the concentration of one species induces a flux of another chemical species, has generally been neglected in the study of reaction-diffusion systems. It is worth mentioning that nowadays there is still no general theory available that covers all possible cross-diffusion models, even in the simplest case of only two coupled partial differential equations. Furthermore, let us point out that cross-diffusion effects are not that well studied in literature. In this particular context, there are very few contributions dealing, some extent, with cross-diffusion terms (see \cite{benjmpa,chen,del,gal,horst,ko,ruiz,wang,wu}) and also the references therein. Being inspired by the mathematical reaction-diffusion system with cross-diffusion terms \cite{benjee}, we show the existence of equivalence of renormalized and entropy solutions of the following parabolic system \begin{equation} \label{1} \begin{aligned} &\frac{\partial u_1}{\partial t} - \operatorname{div}(A_1(u_1,\nabla u_1))\\ &- \operatorname{div}\Big((k_1u_1+ \delta_{1,2}u_2+\delta_{1,3}u_3)\nabla u_1 + \delta'_{1,2}u_1\nabla u_2+ \delta'_{1,3}u_1\nabla u_3\Big) \\ &= -\sigma(u_1,u_2,u_3) - \mu_1u_1 +f \quad\text{in }Q_T, \\ &\frac{\partial u_2}{\partial t}- \operatorname{div}(A_2(u_2,\nabla u_2))\\ &- \operatorname{div}\Big((\delta_{2,1}u_1+ k_2u_2+\delta_{2,3}u_3)\nabla u_2 + \delta'_{2,1}u_2\nabla u_1+ \delta'_{2,3}u_2\nabla u_3\Big) \\ &= \sigma(u_1,u_2,u_3) - \varrho u_2 - \mu_2u_2 +g \quad\text{in }Q_T, \\ &\frac{\partial u_3}{\partial t}- \operatorname{div}(A_3(u_3,\nabla u_3))\\ &- \operatorname{div}\Big(\big(\delta_{3,1}u_1+ \delta_{3,2}u_2+k_3u_3\big)\nabla u_3 + \delta'_{3,1}u_3\nabla u_1+ \delta'_{3,2}u_3\nabla u_2\Big) \\ &= \varrho u_2 - \mu_3u_3 +h \quad\text{in }Q_T, \\ \end{aligned} \end{equation} with the initial and boundary conditions \begin{gather*} u_i(x,0)=u_{i,0}(x) \quad\text{in }\Omega,\; i=1,2,3, \\ u_i(x,t)= 0 \quad\text{on }\Sigma_T,\; i=1,2,3, \end{gather*} where $Q_T = \Omega \times (0,T)$, $\Sigma_T = \partial \Omega \times (0,T)$, $\Omega $ is a bounded domain in $\mathbb{R}^N$ with boundary $\partial \Omega$ (no smoothness assumed on the boundary $\partial \Omega$) and $T>0$. We consider the propagation of an epidemic disease in a spatially distributed population wherein $u_1(x,t)$ represents the density of susceptible individuals, $u_2(x,t)$ represents the density of infected individuals, and $u_3(x,t)$ represents the density of recovered individuals at location $x$ and time $t$. Here the homogeneous Dirichlet boundary condition means that the model \eqref{1} is self-contained and has no population on the boundary $\partial \Omega$. In the above system \eqref{1}, the diffusion operator is non-linear, $k_i>0$, $i=1,2,3$, denote the self-diffusion rates and $\delta {i,j}$, $\delta'_{i,j}$, $i\neq j$ are cross-diffusion rates. The duration of infectious stage is given by $1/r>0$ and the mortality rate $\mu_i>0$, $i=1,2,3$. There is no fertility function so that the system \eqref{1} models the evolution of a given group of individuals. The incidence function $\sigma(u_1,u_2,u_3)$, yielding the recruitment of newly infected individuals from the susceptible class takes a proportionate mixing or a frequency dependence form $\sigma(u_1,u_2,u_3)= \sigma_1 \frac{u_1u_2}{u_1+u_2+u_3}$ for some $\sigma_1>0$, and $u_1,u_2,u_3\geq 0$. According to the definition of the incidence function $\sigma$, it is easy to realize that $\sigma(u_1,u_2,u_3) \leq \alpha \min(u_1,u_2)$ for $u_1,u_2,u_3\geq 0$. For more details regarding the incidence term and their properties, one can see \cite{bendaaam,benjee,fitana,fitjam,fitsiam,fitfields}. For sake of simplicity in this work, without loss of generality, we assume that the coefficients $k_i \geq 1$ and $ \delta_{i,j}= \delta'_{i,j} =1$, $i\neq j$, $i,j=1,2,3$. It is worth mentioning that to prove the existence of renormalized and entropy solutions to system \eqref{1}, we are in need of the following hypotheses throughout this article. The divergence form operator $A_i(s,\zeta) :\mathbb{R}\times \mathbb{R}^N \to \mathbb{R}^N$ is a Caratheodory function (that is, it is continuous with respect to $s$ and $\zeta$) such that \begin{itemize} \item [(H1)] $A_i(s,\zeta) \zeta \geq \alpha_i |\zeta|^2$, for every $\zeta \in \mathbb{R}^N$, where $\alpha_i >0$ and $i=1,2,3$; \item [(H2)] For any $k>0$, there exists $\beta_k>0$ and a function $C_k(x,t) \in L^2(Q_T)$ such that $|A_i(s,\zeta)| \leq C_k(x,t) + \beta_k |\zeta|$, $i=1,2,3$; \item [(H3)] $[A_i(s,\zeta)- A_i(s,\zeta')][\zeta - \zeta'] \geq 0$, $i=1,2,3$. \item [(H4)] $u_{i,0}(x) \in L^1(\Omega)$, $i=1,2,3$. \item [(H5)] $f(x,t),g(x,t),h(x,t) \in L^1(Q_T)$. \end{itemize} for almost every $(x,t)$ in $Q_T$, for every $s$ in $\mathbb{R}$ and for every $\zeta, \zeta'$ in $\mathbb{R}^N, (\zeta \neq \zeta')$. Under these assumptions, establishing the weak solutions to the system \eqref{1} in the sense of distribution is really a tough task. To overcome this difficulty, we use the framework of renormalized solutions was introduced by Diperna and Lions \cite{lion} for Boltzmann equation and also the different notion of entropy solutions introduced and studied by B\'enilan et al.~\cite{benailan}. As mentioned earlier one of our mail goal of this paper is establishing the existence of renormalized solutions and further we want show the equivalence of renormalized solutions and entropy solutions of the cross-diffusion system \eqref{1}. As far as the notion of renormalized solution is concerned, the literature is very vast. After the pioneering work of Diperna and Lions \cite{lion}, the notion is strongly followed by Blanchard et al.\ for nonlinear parabolic equations with integrable data \cite{blanjde} and, for nonlinear elliptic and parabolic equations by Boccardo et al.~\cite{bocgjfa,bocjfa,bocaa}, see the references therein for more details. Apart from the literature mentioned above for the existence of solutions of partial differential equations, to the best of our knowledge, as far as the existence of solutions of parabolic system is concerned, only few articles have appeared; for example, Bendahmane and Karlsen established the existence of renormalized solutions of the reaction-diffusion system with $L^1$ data \cite{bencpaa,bensiam}, Bendahmane and Saad studied the existence of solutions of predator-prey system with $L^1$ data in \cite{benjmaa} and Bendahmane et al. proved existence of solutions of reaction-diffusion epidemic disease system with $L^1$ data in \cite{benade}. On the other hand regarding the equivalence between the renormalized solutions and entropy solutions only two papers are available in the literature for parabolic equations, see \cite{dro,zhang} and also for the system of parabolic equations concerned there is no paper available in the literature. Accordingly in our paper, in contrast to above results, we study the existence and equivalence of the renormalized and entropy solutions with no growth conditions and integrable data. It should be remarked that in the entire paper we use the generic constant $c$ instead of different constants. Now first we define the truncation function that we have used throughout in this paper and after that we give the definition of the renormalized solutions and the entropy solutions of the reaction-diffusion system with cross-diffusion terms \eqref{1} \begin{gather*} T_k(z) =\begin{cases} k & \text{if }z\geq k, \\ z &\text{if }|z| \leq k,\\ -k &\text{if }z \leq -k, \end{cases}\\ \tilde{T}_k(z) = \int ^z _0 T_k(s) ds =\begin{cases} \frac{z^2}{2} &\text{if }|z| \leq k, \\ k|z| - \frac{k^2}{2} &\text{if }|z| \geq k. \end{cases} \end{gather*} \begin{definition}\label{d1} \rm A renormalized solution of \eqref{1} is a triple function $(u_1,u_2,u_3)$ satisfying the following conditions: $u_1,u_2,u_3 \geq 0$ for a.e in $(x,t) \in Q_T$. For $i=1,2,3$, \begin{gather*} u_i \in L^\infty(0,T;L^1(\Omega)) \cap C([0,T], L^1(\Omega)), \\ T_k(u_i) \in L^2(0,T; H_0^1(\Omega)), \quad \text{for any } k \geq 0, \\ \int_{\{ n \leq |u_i| \leq n+1\}} A_i(u_i,\nabla u_i) \,dx\,dt \to 0 \quad \text{as } n \to \infty. \end{gather*} For all $S(u_i) \in C^\infty(\mathbb{R})$ with $\operatorname{supp}S'$ compact, \begin{equation} \label{2} \begin{aligned} &\frac{\partial S(u_1)}{\partial t} - \operatorname{div}(S' (u_1)A_1(u_1,\nabla u_1)) + S''(u_1)A_1(u_1,\nabla u_1) \nabla u_1 \\ &- \operatorname{div}\Big(S'(u_1) ((k_1u_1+u_2+u_3)\nabla u_1 +u_1\nabla u_2 +u_1 \nabla u_3)\Big) \\ & + S''(u_1)\Big((k_1u_1+u_2+u_3)\nabla u_1 +u_1\nabla u_2+u_1 \nabla u_3\Big) \nabla u_1 \\ &= (-\sigma(u_1,u_2,u_3) - \mu_1u_1 +f) S' (u_1) \quad\text{in }D'(Q_T), \\ &\frac{\partial S(u_2)}{\partial t} - \operatorname{div}(S' (u_2)A_2(u_2,\nabla u_2)) + S''(u_2)A_2(u_2,\nabla u_2) \nabla u_2\\ &- \operatorname{div}\Big(S'(u_2) ((u_1+k_2u_2+u_3)\nabla u_2 +u_2\nabla u_1 +u_2 \nabla u_3)\Big)\\ & + S''(u_2)((u_1+k_2u_2+u_3)\nabla u_2+ u_2\nabla u_1+u_2 \nabla u_3)\nabla u_2\\ & = (\sigma(u_1,u_2,u_3) - \varrho u_2 - \mu_2u_2 +g) S' (u_2) \quad\text{in }D'(Q_T), \\ &\frac{\partial S(u_3)}{\partial t} - \operatorname{div}(S' (u_3)A_3(u_3,\nabla u_3)) + S''(u_3)A_3(u_3,\nabla u_3) \nabla u_3 \\ &- \operatorname{div}\Big(S'(u_3) ((u_1+u_2+k_3u_3)\nabla u_3 +u_3\nabla u_1 +u_3 \nabla u_2)\Big) \\ & + S''(u_3)((u_1+u_2+k_3u_3)\nabla u_3+u_3\nabla u_1+u_3 \nabla u_2)\nabla u_3 \\ & = (\varrho u_2 - \mu_3u_3 +h) S' (u_3) \quad\text{in }D'(Q_T), \end{aligned} \end{equation} and the initial conditions $ S(u_i(x,o))= S(u_{i,0}(x))$, $i=1,2,3$, in $\Omega$ hold. \end{definition} \begin{definition}\label{d2} \rm An entropy solution of \eqref{1} is a triple function $(u_1,u_2,u_3)$ satisfying the following conditions, that is, for $i=1,2,3$, $$ u_i\in L^\infty (0,T;L^1(\Omega)) \cap C([0,T],L^1(\Omega)). $$ For any $k>0$ and for all $\phi_i \in C^1(Q_T)$ with $\phi_i = 0$ in $\Sigma_T$, % \label{3} \begin{align*} &\int_\Omega \tilde{T}_k(u_1-\phi_1) (T) dx - \int_\Omega \tilde{T}_k(u_1-\phi_1)(0) dx + \int_0^T \langle \phi_{1t}, T_k(u_1-\phi_1)\rangle dt \\ &+ \int _{Q_T} A_1(u_1,\nabla u_1) \nabla T_k(u_1-\phi_1) \,dx\,dt + \int_{Q_T}((k_1u_1+u_2+u_3)\nabla u_1 + u_1\nabla u_2 \\ &+ u_1 \nabla u_3)\nabla T_k(u_1-\phi_1) \,dx\,dt \\ &= \int_{Q_T} (-\sigma(u_1,u_2,u_3) - \mu_1u_1 +f) T_k(u_1-\phi_1) \,dx\,dt, \\ &\int_\Omega \tilde{T}_k(u_2-\phi_2) (T) dx - \int_\Omega \tilde{T}_k(u_2-\phi_2)(0) dx + \int_0^T \langle \phi_{2t}, T_k(u_2-\phi_2)\rangle dt \\ &+ \int _{Q_T} A_2(u_2,\nabla u_2) \nabla T_k(u_2-\phi_2) \,dx\,dt +\int_{Q_T}((u_1+k_2u_2+u_3)\nabla u_2 + u_2\nabla u_1 \\ & + u_2 \nabla u_3)\nabla T_k(u_2-\phi_2) \,dx\,dt\\ &= \int_{Q_T} (\sigma(u_1,u_2,u_3)- \varrho u_2 - \mu_2u_2 +g) T_k(u_2-\phi_2) \,dx\,dt, \\ &\int_\Omega \tilde{T}_k(u_3-\phi_3) (T) dx - \int_\Omega \tilde{T}_k(u_3-\phi_3)(0) dx + \int_0^T \langle \phi_{3t}, T_k(u_3-\phi_3)\rangle dt \\ &+ \int _{Q_T} A_3(u_3,\nabla u_3) \nabla T_k(u_3-\phi_3) \,dx\,dt + \int_{Q_T}((u_1+u_2+k_3u_3)\nabla u_3 + u_3\nabla u_1 \\ &+ u_3 \nabla u_2)\nabla T_k(u_3-\phi_3) \,dx\,dt \\ &= \int_{Q_T} (\varrho u_2 - \mu_3u_3 +h) T_k(u_1-\phi_1) \,dx\,dt, \end{align*} hold. \end{definition} It should be remarked that as far as the equivalence between the renormalized and the entropy solutions of reaction-diffusion system with cross-diffusion terms is concerned, there is no paper available in the literature. In this connection first we introduce the main results of the paper that is, the theorems which concern the existence of renormalized solutions, the existence of entropy solutions and finally the equivalence of the renormalized and entropy solutions. \begin{theorem} \label{t1} Under the hypotheses {\rm (H1)--(H5)}, there exists a renormalized solution of \eqref{1} in the sense of Definition \ref{d1}. \end{theorem} \begin{theorem} \label{t2} Under the hypotheses {\rm (H1)--(H5)}, the renormalized solution of \eqref{1} is also an entropy solution of the same system in the sense of Definition \ref{d2}. \end{theorem} \begin{remark}\label{r1}\rm The renormalized solution of \eqref{1} is equivalent to the entropy solution of given cross-diffusion epidemic system \eqref{1}. \end{remark} \begin{remark}\rm We also stress that the same method can be applied for density-dependent diffusion or $p$-Laplacian type diffusion under standard variational growth constraint. For parabolic equations with $p$-Laplacian diffusion, for example, see \cite{zhang} and for a system of nonlinear partial differential equations with anisotropic diffusivities, for example, see \cite{bencpaa}. \end{remark} \begin{remark} \rm The uniqueness of the renormalized and entropy solutions of the given cross-diffusion system, by following the classical methods is tough task due to the given hypothesis and cross-diffusion terms of the problem. However, there is an another method, which was introduced by Blanchard et al.~\cite{blanjde} for the uniqueness of renormalized solution of parabolic equation, but it cannot be applied here to prove the uniqueness of the given parabolic system. \end{remark} The article is organized as follows. In Section 2, we introduce the regularized problem of the system \eqref{1} and we prove the existence of solutions of regularized problem using Galerkin's approximation method. Also we establish some basic lemmas that we have used to prove the main results of the paper. In Section 3, we prove Theorem \ref{t1}, that is, the existence of renormalized solutions and finally, in Section 4, we prove Theorem \ref{t2}, that is, the existence of entropy solutions from the notion of renormalized solutions. \section{Approximation problem} In this section, first we introduce an approximation problem for the given r-action-diffusion system \eqref{1} and then we prove the existence of solutions of the approximation problem. For $\varepsilon>0$, let us introduce the following approximations on the data: \begin{itemize} \item[(H6)] $f^\varepsilon, g^\varepsilon, h^\varepsilon \in L^2(Q_T)$ and $f^\varepsilon \to f, g^\varepsilon\to g$ and $ h^\varepsilon \to h$ a.e in $Q_T$ and strongly in $L^1(Q_T)$ as $\varepsilon$ tends to zero; \item[(H7)] $u_{i,0}^\varepsilon \in L^2(\Omega), i=1,2,3$, and $u_{i,0}^\varepsilon \to u_{i,0},i=1,2,3$, a.e in $\Omega$ and strongly in $L^1(\Omega)$ as $\varepsilon$ tends to zero. \end{itemize} \begin{equation} \label{4} \begin{aligned} &\frac{\partial u_1^\varepsilon}{\partial t} - \operatorname{div}(A_1(u_1^\varepsilon,\nabla u_1^\varepsilon)) - \operatorname{div}\Big((k_1F^+_\varepsilon(u_1^\varepsilon)+ F^+_\varepsilon(u_2^\varepsilon)+F^+_\varepsilon(u_3^\varepsilon)) \nabla u_1^\varepsilon \\ &+ F^+_\varepsilon(u_1^\varepsilon)\nabla u_2^\varepsilon+ F^+_\varepsilon(u_1^\varepsilon)\nabla u_3^\varepsilon\Big) \\ &= -\sigma(u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon) - \mu_1u_1^\varepsilon +f^\varepsilon \quad\text{in }Q_T, \\ &\frac{\partial u_2^\varepsilon}{\partial t} - \operatorname{div}(A_2(u_2^\varepsilon,\nabla u_2^\varepsilon)) - \operatorname{div}\Big((F^+_\varepsilon(u_1^\varepsilon)+ k_2F^+_\varepsilon(u_2^\varepsilon) +F^+_\varepsilon(u_3^\varepsilon))\nabla u_2^\varepsilon \\ &+ F^+_\varepsilon(u_2^\varepsilon)\nabla u_1^\varepsilon+F^+_\varepsilon(u_2^\varepsilon)\nabla u_3^\varepsilon\Big)\\ & = \sigma(u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon) - \varrho u_2^\varepsilon - \mu_2u_2^\varepsilon +g^\varepsilon \quad\text{in }Q_T, \\ &\frac{\partial u_3^\varepsilon}{\partial t} - \operatorname{div}(A_3(u_3^\varepsilon,\nabla u_3^\varepsilon)) - \operatorname{div}\Big(\big(F^+_\varepsilon(u_1^\varepsilon)+ F^+_\varepsilon(u_2^\varepsilon)+k_3F^+_\varepsilon(u_3^\varepsilon)\big) \nabla u_3^\varepsilon \\ &+ F^+_\varepsilon(u_3^\varepsilon)\nabla u_1^\varepsilon + F^+_\varepsilon(u_3^\varepsilon)\nabla u_2^\varepsilon\Big)\\ &= \varrho u_2^\varepsilon - \mu_3u_3^\varepsilon +h ^\varepsilon\quad\text{in }Q_T, \end{aligned} \end{equation} with the initial and boundary conditions \begin{gather*} u_i^\varepsilon(x,0) =u_{i,0}^\varepsilon(x) \quad\text{in }\Omega, i=1,2,3, \\ u_i^\varepsilon(x,t) = 0 \quad\text{on }\Sigma_T, i=1,2,3, \end{gather*} where $ F^+_\varepsilon(a) = \max (0, \frac{a}{1+\varepsilon |a|})$. \begin{remark}\label{r2} \rm It is easy to understand that, from \cite{benjee}, the diffusion matrix of \eqref{4} can be replaced by \begin{align*} &\frac{\partial u_1^\varepsilon}{\partial t} - \operatorname{div}(A_1(u_1^\varepsilon,\nabla u_1^\varepsilon)) - \operatorname{div}(\alpha_{1,1}\nabla u_1^\varepsilon + \alpha_{1,2}\nabla u_2^\varepsilon+ \alpha_{1,3}\nabla u_3^\varepsilon)\\ & = -\sigma(u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon) - \mu_1u_1^\varepsilon +f^\varepsilon, \\ &\frac{\partial u_2^\varepsilon}{\partial t} - \operatorname{div}(A_2(u_2^\varepsilon,\nabla u_2^\varepsilon)) - \operatorname{div}(\alpha_{2,1}\nabla u_1^\varepsilon+ \alpha_{2,2}\nabla u_2^\varepsilon +\alpha_{2,3}\nabla u_3^\varepsilon)\\ & = \sigma(u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon) - \varrho u_2^\varepsilon - \mu_2u_2^\varepsilon +g^\varepsilon , \\ &\frac{\partial u_3^\varepsilon}{\partial t}- \operatorname{div}(A_3(u_3^\varepsilon,\nabla u_3^\varepsilon)) - \operatorname{div}(\alpha_{3,1}\nabla u_1^\varepsilon + \alpha_{3,2}\nabla u_2^\varepsilon+\alpha_{3,3}\nabla u_3^\varepsilon) \\ &= \varrho u_2^\varepsilon - \mu_3u_3^\varepsilon +h ^\varepsilon , \end{align*} where $\alpha_{i,i} = k_iF^+_\varepsilon(u^\varepsilon_i) + \sum _{j=1, j\neq i} ^3 F^+_\varepsilon (u^\varepsilon_j)$ for $i=1,2,3$. Now the diffusion matrix \[ \mathcal{M} = \begin{pmatrix} \alpha_{1,1} & F^+_\varepsilon(u_1^\varepsilon) & F^+_\varepsilon(u_1^\varepsilon) \\ F^+_\varepsilon(u_2^\varepsilon) & \alpha_{2,2} & F^+_\varepsilon(u_2^\varepsilon) \\ F^+_\varepsilon(u_3^\varepsilon) & F^+_\varepsilon(u_3^\varepsilon) & \alpha_{3,3} \end{pmatrix} \] is uniformly non-negative from the definition of the approximation problem. Using the assumptions on $k_i$, $ i=1,2,3$, and from the inequality $ab \geq -\frac{a^2}{2} - \frac{b^2}{2}$ for all $a,b \in \mathbb{R}$, we obtain \begin{equation} \label{5} \begin{aligned} \zeta^T \mathcal{M}\zeta &= \sum_{i=1}^3 \big(k_iF^+_\varepsilon(u_i^\varepsilon) + \sum_{j=1,j\neq i}^3 F^+_\varepsilon(u^\varepsilon_j)\big)\zeta_i^2 + F^+_\varepsilon (u_1^\varepsilon)(\zeta_1 \zeta_2 + \zeta_1\zeta_3) \\ &\quad + F^+_\varepsilon(u_2^\varepsilon) (\zeta_1\zeta_2+\zeta_3\zeta_2) + F^+_\varepsilon(u_3^\varepsilon) (\zeta_1 \zeta_3+\zeta_2\zeta_3) \\ &\geq \sum_{i=1}^3 (k_i-1) F^+_\varepsilon(u_i^\varepsilon) \zeta_i^2 + \sum_{i=1}^3\sum_{j=1,j\neq i}^3 \frac{F^+_\varepsilon(u_j)}{2} \zeta_i^2 \geq 0. \end{aligned} \end{equation} \end{remark} \begin{lemma}\label{l1} Assume that $u_{i,0}^\varepsilon \in L^2(\Omega)$, $i=1,2,3$, and $f^\varepsilon, g^\varepsilon, h^\varepsilon \in L^2(Q_T)$. Then the approximation problem \eqref{4} admits a unique weak solution $$ u_i^\varepsilon \in L^\infty(0,T;L^2(\Omega)) \cap L^2(0,T;H_0^1(\Omega)) \cap C([0,T],L^2(\Omega)), \quad i=1,2,3, $$ with $u^\varepsilon_{it} \in L^2(0,T; H^{-1}(\Omega))$ such that, for any $\phi_i \in L^2(0,T;H^1_0(\Omega))$, $ i=1,2,3$, \begin{align*} %\label{6} &\int_0^T \langle \partial_t u_1^\varepsilon, \phi_1\rangle dt + \int_{Q_T}A_1(u^\varepsilon_1,\nabla u_1^\varepsilon) \nabla \phi_1 \,dx\,dt \\ &+ \int_{Q_T} (\alpha_{1,1} \nabla u_1^\varepsilon +\alpha_{1,2} \nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon) \nabla \phi_1 \,dx\,dt \\ & = \int_{Q_T} (-\sigma(u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon) -\mu_1 u_1^\varepsilon +f^\varepsilon) \phi_1 \,dx\,dt, \\ &\int_0^T \langle \partial_t u_2^\varepsilon, \phi_2\rangle dt + \int_{Q_T}A_2(u^\varepsilon_2,\nabla u_2^\varepsilon) \nabla \phi_2 \,dx\,dt \\ &+ \int_{Q_T}(\alpha_{2,1}\nabla u_1^\varepsilon +\alpha_{2,2}\nabla u_2^\varepsilon + \alpha_{2,3}\nabla u_3^\varepsilon) \nabla \phi_2 \,dx\,dt \\ &= \int_{Q_T}(\sigma (u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon)-\varrho u_2^\varepsilon -\mu_2 u_2^\varepsilon +g^\varepsilon) \phi_2 \,dx\,dt, \\ &\int_0^T \langle \partial_t u_3^\varepsilon, \phi_3\rangle dt + \int_{Q_T}A_3(u^\varepsilon_3,\nabla u_3^\varepsilon) \nabla \phi_3 \,dx\,dt \\ &+ \int_{Q_T} (\alpha_{3,1}\nabla u_1^\varepsilon + \alpha_{3,2}\nabla u_2^\varepsilon + \alpha_{3,3} \nabla u_3^\varepsilon) \nabla \phi_3 \,dx\,dt \\ &= \int_{Q_T}( \varrho u_2^\varepsilon - \mu_3 u_3^\varepsilon + h^\varepsilon) \phi_3 \,dx\,dt, \end{align*} hold. \end{lemma} \begin{proof} The proof of this lemma is not new, we present here a self-contained sketch of the proof for the sake of simplicity and readability. For more details regarding this proof, one can refer the article \cite{benjee} and some details related to the following method see \cite{evans}. Here the proof is relies on using the Faedo-Galerkin approximation method. To use this method, a specific basis is required. For that we consider the approximate spectral problem, see \cite{benjee,evans} and for that problem, the corresponding eigen functions $w_l(x)$ form an orthogonal basis in $H_0^1(\Omega)$ and an orthonormal basis in $L^2(\Omega)$. Now we look for the finite-dimensional approximate solutions to problem \eqref{4} in the form of sequences $\{ u_{i,n}^\varepsilon\}$, $i=1,2,3$, defined for $t\geq 0$ and $ x \in \bar{\Omega}$ by $$ u^\varepsilon_{i,n} = \sum _{l=1}^n C_{i,n,l} (t) w_l(x), \quad i=1,2,3. $$ Our aim is to determine the set of coefficients $\{C_{i,n,l}\}_{l=1}^n$, $i=1,2,3$, such that for $m=1,2,\dots,n$, \begin{equation} \label{7} \begin{aligned} & \int_\Omega \partial_t u_{1,n}^\varepsilon w_m\, dx + \int_\Omega A_1(u^\varepsilon_{1,n},\nabla u_{1,n}^\varepsilon) \nabla w_m\, dx \\ &+ \int_\Omega (\alpha_{1,1,n} \nabla u_{1,n}^\varepsilon +\alpha_{1,2,n} \nabla u_{2,n}^\varepsilon + \alpha_{1,3,n}\nabla u_{3,n}^\varepsilon) \nabla w_m\, dx\\ &= \int_\Omega (-\sigma(u_{1,n}^\varepsilon,u_{2,n}^\varepsilon,u_{3,n}^\varepsilon) -\mu_1 u_{1,n}^\varepsilon +f^\varepsilon) w_m\, dx, \\ &\int_\Omega \partial_t u_{2,n}^\varepsilon w_mdx + \int_\Omega A_2(u^\varepsilon_{2,n},\nabla u_{2,n}^\varepsilon) \nabla w_m\, dx \\ &+ \int_\Omega (\alpha_{2,1,n}\nabla u_{1,n}^\varepsilon +\alpha_{2,2,n}\nabla u_{2,n}^\varepsilon + \alpha_{2,3,n}\nabla u_{3,n}^\varepsilon) \nabla w_m\, dx \\ &= \int_\Omega(\sigma (u_{1,n}^\varepsilon,u_{2,n}^\varepsilon,u_{3,n}^\varepsilon)-\varrho u_{2,n}^\varepsilon -\mu_2 u_{2,n}^\varepsilon +g^\varepsilon) w_m\, dx, \\ &\int_\Omega \partial_t u_{3,n}^\varepsilon w_m\, dx + \int_\Omega A_3(u^\varepsilon_{3,n},\nabla u_{3,n}^\varepsilon) \nabla w_m\, dx \\ &+ \int_\Omega (\alpha_{3,1,n}\nabla u_{1,n}^\varepsilon + \alpha_{3,2,n}\nabla u_{2,n}^\varepsilon + \alpha_{3,3,n} \nabla u_{3,n}^\varepsilon) \nabla w_m\, dx \\ & = \int_\Omega( \varrho u_{2,n}^\varepsilon - \mu_3 u_{3,n}^\varepsilon + h^\varepsilon) w_m\, dx, \end{aligned} \end{equation} where $\alpha_{i,i,n} = k_iF^+_\varepsilon(u^\varepsilon_{i,n}) + \sum _{j=1, j\neq i} ^3 F^+_\varepsilon (u^\varepsilon_{j,n})$, for $i=1,2,3$, and with initial conditions $u^\varepsilon_{i,n}(x,0) = u^\varepsilon_{i,0,n}(x) := \sum _{l=1}^n C_{i,n,l}(0)w_l(x)$, $i=1,2,3$. Further it should be remarked that the above form of the basis satisfies the required boundary conditions of the approximation problem \eqref{4}. Now, \eqref{7} can be rewritten in the form \begin{align*} %\label{8} C' _{1,n,m}(t) &= - \int_\Omega A_1(u^\varepsilon_{1,n},\nabla u_{1,n}^\varepsilon) \nabla w_m \,dx\\ &- \int_\Omega (\alpha_{1,1,n} \nabla u_{1,n}^\varepsilon +\alpha_{1,2,n} \nabla u_{2,n}^\varepsilon +\alpha_{1,3,n}\nabla u_{3,n}^\varepsilon)\nabla w_m \,dx \\ &\quad - \int_\Omega (-\sigma(u_{1,n}^\varepsilon,u_{2,n}^\varepsilon,u_{3,n}^\varepsilon) -\mu_1 u_{1,n}^\varepsilon +f^\varepsilon) w_m\, dx \\ &=: G_1^m\left(t, \{C_{1,n,l}\}_{l=1}^n,\{C_{2,n,l}\}_{l=1}^n, \{C_{3,n,l}\}_{l=1}^n\right),\\ \end{align*} \begin{align*} C' _{2,n,m}(t) &= - \int_\Omega A_2(u^\varepsilon_{2,n},\nabla u_{2,n}^\varepsilon) \nabla w_m \,dx\\ &- \int_\Omega (\alpha_{2,1,n}\nabla u_{1,n}^\varepsilon +\alpha_{2,2,n}\nabla u_{2,n}^\varepsilon +\alpha_{2,3,n}\nabla u_{3,n}^\varepsilon) \nabla w_m \,dx \\ &\quad - \int_\Omega(\sigma (u_{1,n}^\varepsilon,u_{2,n}^\varepsilon,u_{3,n}^\varepsilon) -\varrho u_{2,n}^\varepsilon -\mu_2 u_{2,n}^\varepsilon +g^\varepsilon) w_m\, dx \\ &=: G_2^m\left(t, \{C_{1,n,l}\}_{l=1}^n,\{C_{2,n,l}\}_{l=1}^n, \{C_{3,n,l}\}_{l=1}^n\right), \end{align*} \begin{align*} &C' _{3,n,m}(t)\\ &= - \int_\Omega A_3(u^\varepsilon_{3,n},\nabla u_{3,n}^\varepsilon) \nabla w_m\, dx\\ & - \int_\Omega (\alpha_{3,1,n}\nabla u_{1,n}^\varepsilon + \alpha_{3,2,n}\nabla u_{2,n}^\varepsilon + \alpha_{3,3,n} \nabla u_{3,n}^\varepsilon) \nabla w_m\, dx \\ & - \int_\Omega( \varrho u_{2,n}^\varepsilon - \mu_3 u_{3,n}^\varepsilon + h^\varepsilon) w_m\, dx\\ &=: G_3^m\left(t, \{C_{1,n,l}\}_{l=1}^n,\{C_{2,n,l}\}_{l=1}^n,\{C_{3,n,l} \}_{l=1}^n\right). \end{align*} Let $\rho \in (0,T)$ and set $U=[0,\rho]$. Choose $r>0$ large enough so that the ball $B_r \subset \mathbb{R}^n$ contains $\{C_{i,n,l}(0)\},i=1,2,3;$ then set $V= \bar{B_r}$. The components of $G_1$ can be bounded on $U \times V$, from $(H_1-H_7)$, we obtain \begin{align*} &|G_1^m\big(t, \{C_{1,n,l}\}_{l=1}^n,\{C_{2,n,l}\}_{l=1}^n,\{C_{3,n,l}\}_{l=1}^n \big)|\\ &\leq \Big( \int_\Omega |A_1\big(\sum_{l=1}^n C_{1,n,l} w_l, \sum_{l=1}^n C_{1,n,l} \nabla w_l\big)|^2 dx\Big)^{1/2} \Big( \int _\Omega |\nabla w_m|^2 dx\Big)^{1/2}\\ &\quad + \frac{(k_1+2)}{\varepsilon} \Big(\int_\Omega |\sum_{l=1}^n C_{1,n,l} \nabla w_l|^2dx\Big)^{1/2} \Big( \int _\Omega |\nabla w_m|^2 dx\Big)^{1/2} \\ &\quad + \frac{1}{\varepsilon}\Big(\Big(\int_\Omega |\sum_{l=1}^n C_{2,n,l} \nabla w_l|^2dx\Big)^{1/2}\\ &\quad +\Big(\int_\Omega |\sum_{l=1}^n C_{3,n,l} \nabla w_l|^2dx\Big)^{1/2}\Big) \Big( \int _\Omega |\nabla w_m|^2 dx\Big)^{1/2}\\ &\quad + \big(\sigma_1+\mu_1\big) \Big(\int_\Omega |\sum_{l=1}^n C_{1,n,l} w_l|^2dx)^{1/2}\Big) \Big( \int _\Omega | w_m|^2 dx\Big)^{1/2}\\ &\quad + \Big(\int_\Omega |f^\varepsilon|^2 dx \Big)^{1/2} \Big( \int _\Omega | w_m|^2 dx\Big)^{1/2} \\ &\leq c(r,n). \end{align*} where the constant depends only on $r$ and $n$. Similarly one can easily obtain \begin{gather*} |G_2^m\left(t, \{C_{1,n,l}\}_{l=1}^n,\{C_{2,n,l}\}_{l=1}^n, \{C_{3,n,l}\}_{l=1}^n\right)| \leq c(r,n),\\ |G_3^m\left(t, \{C_{1,n,l}\}_{l=1}^n,\{C_{2,n,l}\}_{l=1}^n, \{C_{3,n,l}\}_{l=1}^n\right)| \leq c(r,n). \end{gather*} Then the standard ODE theory shows that $\{C_{i,n,l}\}_{l=1}^n$, $i=1,2,3$, respectively satisfies \eqref{7} for a.e $t \in [0,\rho')$. Moreover we have \begin{align*} &C_{1,n,l}(t) \\ &= C_{1,n,l}(0) + \int_0^t G_1^l(\tau,\{C_{1,n,m}(\tau)\}_{m=1}^n,\{C_{2,n,m}(\tau)\}_{m=1}^n, \{C_{3,n,m}(\tau)\}_{m=1}^n) d \tau,\\ &C_{2,n,l}(t) \\ &= C_{2,n,l}(0) + \int_0^t G_2^l(\tau,\{C_{1,n,m}(\tau)\}_{m=1}^n,\{C_{2,n,m}(\tau)\}_{m=1}^n, \{C_{3,n,m}(\tau)\}_{m=1}^n) d \tau,\\ &C_{3,n,l}(t) \\ &= C_{3,n,l}(0) + \int_0^t G_3^l(\tau,\{C_{1,n,m}(\tau)\}_{m=1}^n,\{C_{2,n,m}(\tau)\}_{m=1}^n, \{C_{3,n,m}(\tau)\}_{m=1}^n) d \tau. \end{align*} This proves that the functions $(u_{1,n}^\varepsilon, u_{2,n}^\varepsilon, u_{3,n}^\varepsilon)$ are well-defined and approximate solutions to the problem \eqref{4} on $[0,\rho')$. Set $\phi_{i,n}(x,t) = \sum _{l=1}^n b_{i,n,l}(t) w_l(x),i=1,2,3$ where the coefficients $b_{i,n,l}, i=1,2,3$, are absolutely continuous functions. Then, from \eqref{7}, the approximate solutions satisfy the weak formulation \begin{equation} \label{9} \begin{aligned} &\int_\Omega \partial_t u_{1,n}^\varepsilon \phi_{1,n} dx + \int_\Omega A_1(u^\varepsilon_{1,n},\nabla u_{1,n}^\varepsilon) \nabla \phi_{1,n} dx \\ &+ \int_\Omega (\alpha_{1,1,n} \nabla u_{1,n}^\varepsilon +\alpha_{1,2,n} \nabla u_{2,n}^\varepsilon + \alpha_{1,3,n}\nabla u_{3,n}^\varepsilon) \nabla \phi_{1,n} dx \\ & = \int_\Omega (-\sigma(u_{1,n}^\varepsilon,u_{2,n}^\varepsilon,u_{3,n}^\varepsilon) -\mu_1 u_{1,n}^\varepsilon +f^\varepsilon) \phi_{1,n} dx, \\ &\int_\Omega \partial_t u_{2,n}^\varepsilon \phi_{2,n}dx + \int_\Omega A_2(u^\varepsilon_{2,n},\nabla u_{2,n}^\varepsilon) \nabla \phi_{2,n} dx \\ &+ \int_\Omega (\alpha_{2,1,n}\nabla u_{1,n}^\varepsilon +\alpha_{2,2,n}\nabla u_{2,n}^\varepsilon + \alpha_{2,3,n}\nabla u_{3,n}^\varepsilon) \nabla \phi_{2,n} dx \\ &= \int_\Omega(\sigma (u_{1,n}^\varepsilon,u_{2,n}^\varepsilon,u_{3,n}^\varepsilon) -\varrho u_{2,n}^\varepsilon -\mu_2 u_{2,n}^\varepsilon +g^\varepsilon) \phi_{2,n} dx, \\ &\int_\Omega \partial_t u_{3,n}^\varepsilon \phi_{3,n} dx + \int_\Omega A_3(u^\varepsilon_{3,n},\nabla u_{3,n}^\varepsilon) \nabla \phi_{3,n} dx \\ &+ \int_\Omega (\alpha_{3,1,n}\nabla u_{1,n}^\varepsilon + \alpha_{3,2,n}\nabla u_{2,n}^\varepsilon + \alpha_{3,3,n} \nabla u_{3,n}^\varepsilon) \nabla \phi_{3,n} dx \\ & = \int_\Omega( \varrho u_{2,n}^\varepsilon - \mu_3 u_{3,n}^\varepsilon + h^\varepsilon) \phi_{3,n} dx. \end{aligned} \end{equation} From now on, $\tilde{T}$ is an arbitrary time in the existence interval $[0,\rho')$ of Faedo-Galerkin solutions. Take $ \phi_{i,n}= u^\varepsilon_{i,n}, i=1,2,3$, respectively in \eqref{9} and using, Gronwall's lemma and Young's inequality \eqref{5}, we obtain \begin{gather*} \|u^\varepsilon_{i,n}\|_{L^\infty(0,\tilde{T};L^2(\Omega))} + \|u^\varepsilon_{i,n}\|_{L^2(0,\tilde{T};H^1_0(\Omega))} \leq c, \\ \| \partial _t u^\varepsilon_{i,n}\|_{L^2(0,\tilde{T};H^{-1}(\Omega))} + \|\sigma(u^\varepsilon_{1,n}, u^\varepsilon_{2,n}, u^\varepsilon_{3,n})\|_{L^2(Q_T)} \leq c, \end{gather*} for some constant $c>0$ and $i=1,2,3$. Moreover one can easily show, using similar approach of \cite{benjee} with the above estimates, the global existence of approximate solutions of the problem \eqref{4}. Hence, from the above arguments as $n \to \infty$, for $i=1,2,3$, we obtain \begin{gather*} u^\varepsilon_{i,n} \rightharpoonup u^\varepsilon_i \quad \text{weakly-$*$ in } L^\infty(Q_T), \\ u^\varepsilon_{i,n} \to u^\varepsilon_i \quad \text{a.e in $Q_T$ and strongly in } L^2(Q_T), \\ u^\varepsilon_{i,n} \rightharpoonup u^\varepsilon_i \quad \text{weakly in } L^2(0,T;H^1_0(\Omega)), \\ A_i(u^\varepsilon_{i,n},\nabla u^\varepsilon_{i,n}) \rightharpoonup \eta_i \quad \text{weakly in } L^2(Q_T), \\ \partial_t u^\varepsilon_{i,n} \rightharpoonup \partial_t u^\varepsilon_i \quad \text{weakly in } L^2(0,T;H^{-1}(\Omega)). \end{gather*} Using similar type of arguments as in \cite{bennhm} with the monotonicity assumption on $A_i$, we can show that $A_i(u_i^\varepsilon, \nabla u^\varepsilon_i) = \eta_i$, $i=1,2,3$. Since the solutions $u^\varepsilon_i \in L^2(0,T;H_0^1(\Omega)) \cap L^\infty(0,T;L^2(\Omega))$, $i=1,2,3$, using the approximation problem, we conclude that $u^\varepsilon_i \in C([0,T],L^2(\Omega)),i=1,2,3$. This establishes the existence of weak solutions $(u_1^\varepsilon, u_2^\varepsilon,u_3^\varepsilon)$ of the regularized problem \eqref{4}. \end{proof} \begin{lemma}\label{l2} Under hypotheses {\rm (H6)} and {\rm (H7)}, the functions $T_k(u^\varepsilon_i)$ and $\frac{\partial S(u^\varepsilon_i)}{\partial t}$, $i=1,2,3$, are bounded in $L^2(0,T;H^1_0(\Omega))$ and $L^1(Q_T) \cap L^2(0,T;H^{-1}(\Omega))$ respectively. \end{lemma} \begin{proof} Taking $T_k(u^\varepsilon_1)$ as a test function in the first equation of \eqref{4} and integrating over $Q_t = \Omega \times (0,t)$, we obtain \begin{equation} \label{10} \begin{aligned} &\int_{Q_t} u_{1s}^\varepsilon T_k(u_1^\varepsilon) dx ds + \int_{Q_t} A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla T_k(u^\varepsilon_1)\,dx\,ds\\ &+ \int_{Q_t} (\alpha_{1,1} \nabla u_1^\varepsilon+ \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon)\nabla T_k(u^\varepsilon_1) \,dx\,ds \\ &= \int_{Q_t} (-\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)-\mu_1u_1^\varepsilon +f^\varepsilon) T_k(u^\varepsilon_1)\,dx\,ds, \\ &\int_\Omega \tilde{T}_k(u^\varepsilon_1)(t) dx + \alpha_1 \int_{Q_t} |\nabla T_k(u_1^\varepsilon)|^2 \,dx\,ds \\ &+ \int_{Q_t}(\alpha_{1,1} \nabla u_1^\varepsilon+ \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon) \nabla T_k(u^\varepsilon_1) \,dx\,ds \\ &\leq \int_\Omega \tilde{T}_k(u^\varepsilon_{1,0})(x) dx + \int_{Q_t} (-\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3) -\mu_1u_1^\varepsilon+f^\varepsilon) T_k(u^\varepsilon_1)\,dx\,ds. \end{aligned} \end{equation} Similarly by considering the second and third equations of \eqref{4}, we obtain \begin{equation} \label{10a} \begin{aligned} &\int_\Omega \tilde{T}_k(u^\varepsilon_2)(t) dx + \alpha_2 \int_{Q_t} |\nabla T_k(u_2^\varepsilon)|^2 \,dx\,ds\\ & + \int_{Q_t}(\alpha_{2,1} \nabla u_1^\varepsilon+ \alpha_{2,2}\nabla u_2^\varepsilon + \alpha_{2,3}\nabla u_3^\varepsilon)\nabla T_k(u^\varepsilon_2) \,dx\,ds \\ &\leq \int_\Omega \tilde{T}_k(u^\varepsilon_{2,0})(x) dx + \int_{Q_t}(\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)-\varrho u_2^\varepsilon - \mu_2 u_2^\varepsilon + g^\varepsilon) T_k(u_2^\varepsilon) \,dx\,ds, \\ &\int_\Omega \tilde{T}_k(u^\varepsilon_3)(t) dx + \alpha_3 \int_{Q_t} |\nabla T_k(u_3^\varepsilon)|^2 \,dx\,ds \\ &+ \int_{Q_t}(\alpha_{3,1} \nabla u_1^\varepsilon+ \alpha_{3,2}\nabla u_2^\varepsilon + \alpha_{3,3}\nabla u_3^\varepsilon) \nabla T_k(u^\varepsilon_3) \,dx\,ds \\ &\leq \int_\Omega \tilde{T}_k(u^\varepsilon_{3,0})(x) dx + \int_{Q_t} (\varrho u_2^\varepsilon-\mu_3u^\varepsilon_3 +h^\varepsilon)T_k(u^\varepsilon_3) \,dx\,ds. \end{aligned} \end{equation} Summing the above three inequalities and using the Young's inequality, the properties of the functions $f^\varepsilon,g^\varepsilon,h^\varepsilon, u^\varepsilon_{i,0}(x)$, $i=1,2,3$, $\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)$, \eqref{5}, the boundedness of approximate solutions with the definition of the functions $\tilde{T}_k(u_i^\varepsilon),i=1,2,3$, we have \begin{equation} \label{11} \int_{Q_T} |\nabla T_k(u^\varepsilon_i)|^2 \,dx\,ds \leq c, \end{equation} for $i=1,2,3$, and for any constant $c>0$. This proves that $T_k(u^\varepsilon_i)$, $i=1,2,3$, are bounded in $L^2(0,T;H_0^1(\Omega))$. Now, multiplying the first equation of \eqref{4} by $S'(u^\varepsilon_1)$, we obtain \begin{equation} \label{12} \begin{aligned} \frac{\partial S(u^\varepsilon_1)}{\partial t} & = \operatorname{div}\big(S'(u_1^\varepsilon) A_1(u^\varepsilon_1, \nabla u_1^\varepsilon)\big) - S''(u^\varepsilon_1)A_1(u^\varepsilon_1, \nabla u_1^\varepsilon)\nabla u^\varepsilon_1 \\ &\quad + \operatorname{div}\big( S'(u^\varepsilon_1)(\alpha_{1,1} \nabla u_1^\varepsilon + \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon)\big) \\ &\quad - S''(u^\varepsilon_1) (\alpha_{1,1} \nabla u_1^\varepsilon + \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon)\nabla u_1^\varepsilon \\ &\quad + S'(u_1^\varepsilon)(-\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)-\mu_1u_1^\varepsilon+f^\varepsilon). \end{aligned} \end{equation} This may be rewritten in the following way by using the definition of $T_k(u_i^\varepsilon)$, $i=1,2,3$, \begin{equation}\label{13} \begin{aligned} &\frac{\partial S(u_1^\varepsilon)}{\partial t}\\ & = \operatorname{div}(S' (u_1^\varepsilon)A_1(T_k(u_1^\varepsilon),\nabla T_k(u_1^\varepsilon))) - S''(u_1^\varepsilon)A_1(T_k(u_1^\varepsilon),\nabla T_k(u_1^\varepsilon)) \nabla T_k(u_1^\varepsilon) \\ &\quad + \operatorname{div}(S'(u_1^\varepsilon) ((k_1F^+_\varepsilon(T_k(u_1^\varepsilon)) +F^+_\varepsilon(T_k(u_2^\varepsilon))+F^+_\varepsilon(T_k(u_3^\varepsilon)))\nabla T_k(u_1^\varepsilon) \\ &\quad +F^+_\varepsilon(T_k(u_1^\varepsilon))\nabla T_k(u_2^\varepsilon) +F^+_\varepsilon(T_k(u_1^\varepsilon)) \nabla T_k(u_3^\varepsilon))) - S''(u_1^\varepsilon)((k_1F^+_\varepsilon(T_k(u_1^\varepsilon))\\ &\quad +F^+_\varepsilon(T_k(u_2^\varepsilon))+F^+_\varepsilon(T_k(u_3^\varepsilon)))\nabla T_k(u_1^\varepsilon) +F^+_\varepsilon(T_k(u_1^\varepsilon))\nabla T_k(u_2^\varepsilon)\\ &\quad +F^+_\varepsilon(T_k(u_1^\varepsilon)) \nabla T_k(u_3^\varepsilon))\nabla T_k(u_1^\varepsilon) \\ &\quad +(-\sigma(T_k(u_1^\varepsilon),T_k(u_2^\varepsilon),T_k(u_3^\varepsilon)) - \mu_1T_k(u_1^\varepsilon) +f^\varepsilon) S' (u_1^\varepsilon). \end{aligned} \end{equation} For any $S\in C^\infty(\mathbb{R})$ with $\operatorname{supp}S'$ compact, \eqref{13} shows that $ \frac{\partial S(u^\varepsilon_1)}{\partial t}$ is bounded in $L^1(Q_T) \cap L^2(0,T;H^{-1}(\Omega))$, from the result \eqref{11}. Similar arguments on $ \frac{\partial S(u^\varepsilon_i)}{\partial t}$, $i=2,3$, proves the desired result. This completes the proof. \end{proof} \begin{lemma}\label{l3} The solutions $(u^\varepsilon_1,u^\varepsilon_2, u^\varepsilon_3)$ of the regularized problem \eqref{4} are non-negative. \end{lemma} \begin{proof} To prove the non-negativity of the solutions, we consider $u^{-\varepsilon}_i = \sup(-u^\varepsilon_i,0)$, $i=1,2,3$. Now, multiplying the first equation of \eqref{4} by $-T_k(u^{-\varepsilon}_1)$ and integrating over $\Omega$, we obtain \begin{align*} &\frac{d}{dt} \int_\Omega \tilde{T}_k(u^{-\varepsilon}_1)(t) dx + \int_\Omega A_1(u^{-\varepsilon}_1, \nabla u_1^{-\varepsilon}) \nabla T_k(u_1^{-\varepsilon}) dx + \int_\Omega( (k_1F^+_\varepsilon(u_1^{-\varepsilon})+ F^+_\varepsilon(u_2^{-\varepsilon})\\ &+ F^+_\varepsilon(u_3^{-\varepsilon}))\nabla u_1^{-\varepsilon} + F^+_\varepsilon(u_1^{-\varepsilon})\nabla u_2^{-\varepsilon} + F^+_\varepsilon(u_1^{-\varepsilon})\nabla u_3^{-\varepsilon})\nabla T_k(u_1^{-\varepsilon})dx\\ &=\int_\Omega(-\sigma(u_1^{-\varepsilon},u^{-\varepsilon}_2,u^{-\varepsilon}_3)-\mu_1u_1^{-\varepsilon} +f^\varepsilon)T_k(u_1^{-\varepsilon})dx, \\ &\frac{d}{dt} \int_\Omega \tilde{T}_k(u^{-\varepsilon}_1)(t) dx + \alpha_1 \int_\Omega |\nabla T_k(u_1^{-\varepsilon})|^2 dx + \int _\Omega ( (k_1F^+_\varepsilon(u_1^{-\varepsilon})+ F^+_\varepsilon(u_2^{-\varepsilon})\\ & + F^+_\varepsilon(u_3^{-\varepsilon}))\nabla u_1^{-\varepsilon} + F^+_\varepsilon(u_1^{-\varepsilon}) \nabla u_2^{-\varepsilon}+ F^+_\varepsilon(u_1^{-\varepsilon})\nabla u_3^{-\varepsilon})\nabla T_k(u_1^{-\varepsilon})dx \\ &\leq \int_\Omega(\sigma(u_1^{-\varepsilon},u^{-\varepsilon}_2,u^{-\varepsilon}_3) +\mu_1u_1^{-\varepsilon}+f^\varepsilon)T_k(u_1^{-\varepsilon})dx . \end{align*} By considering $-T_k(u_2^{-\varepsilon}), -T_k(u^{-\varepsilon}_3)$ respectively as the test functions of the other two equations of \eqref{4}, summing all the results, using the boundedness of solutions $u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon$, \eqref{5}, the definition of the functions $\sigma(u^\varepsilon_1,u^\varepsilon_2,u^\varepsilon_3), f^\varepsilon,g^\varepsilon,h^\varepsilon$ and finally from the non-negativity of the terms on the LHS in resulting inequalities, we obtain $$ \frac{d}{dt} \int_\Omega \tilde{T}_k(u^{-\varepsilon}_i)(t) dx \leq 0, \quad \text{for } i=1,2,3. $$ The non-negativity of the initial conditions $u^\varepsilon_{i,0},i=1,2,3$, with the above inequalities prove the non-negativity of the solutions $(u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon)$. \end{proof} \begin{definition}\label{d3} \rm Let us define the Lipschitz continuous function \[ \Theta _n (z) = T_{n+1}(z) - T_n(z) = \begin{cases} 0 &\text{if }|z| \leq n, \\ (|z|-n)\operatorname{sgn}(z) &\text{if }n\leq |z| \leq n+1, \\ \operatorname{sgn}(z) &\text{if }|z| \geq n+1. \end{cases} \] \end{definition} \begin{remark}\label{r3} \rm From the above definition, one can easily understand that the function $\Theta(z)$ satisfies $\|\Theta(z)\|_{L^\infty(\mathbb{R})} \leq 1$, for any $n\geq1$ and $\Theta(z) \to 0$, for any $z$ when $n\to \infty$. \end{remark} \begin{lemma}\label{l4} The Lipschitz continuous function $\Theta_n(u_i)$, $i=1,2,3$, for some $n>0$ and $\varepsilon>0$ satisfies \[ \lim _{n\to \infty} \limsup_{\varepsilon\to 0} \int_0^t \int_{(n\leq |u^\varepsilon_i|\leq n+1)} A_i(u^\varepsilon_i,\nabla u^\varepsilon_i) \nabla u^\varepsilon_i \,dx\,ds =0 \] and $\Theta_n(u_i) \to 0$, for $i=1,2,3$, strongly in $L^2(0,T;H^1_0(\Omega))$ as $n\to \infty$. \end{lemma} \begin{proof} Using $\Theta_n(u^\varepsilon_1)$ as a test function in the first equation of \eqref{4} and integrating over $\Omega$ and then over $(0,t)$, we have \begin{equation} \label{14} \begin{aligned} &\int_\Omega \tilde{\Theta}_n(u^\varepsilon_1)(t) dx + \int_{Q_t} A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla \Theta_n(u^\varepsilon_1)\,dx\,ds \\ &+ \int_{Q_t} (\alpha_{1,1} \nabla u_1^\varepsilon+ \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon)\nabla \Theta_n(u^\varepsilon_1) \,dx\,ds \\ &= \int_{Q_t} (-\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)-\mu_1u_1^\varepsilon+f^\varepsilon) \Theta_n(u^\varepsilon_1)\,dx\,ds + \int_\Omega \tilde{\Theta}_n(u^\varepsilon_{1,0}(x))dx, \end{aligned} \end{equation} for almost all $t$ in $(0,T)$ and $\varepsilon < \frac{1}{n+1}$. Since $\tilde{\Theta}_n(u^\varepsilon_i) \geq 0, i=1,2,3$, for all $x\in \Omega$, by taking in account of \eqref{5}, we obtain the inequality \begin{equation} \label{15} \begin{aligned} &\int_{Q_t} A_1(u^\varepsilon_1, \nabla u^\varepsilon_1) \nabla \Theta_n(u^\varepsilon_1) dx ds \\ &\leq \int_{Q_t}(\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)+\mu_1u_1^\varepsilon +f^\varepsilon)\Theta_n(u^\varepsilon_1)\,dx\,ds + \int_\Omega \tilde{\Theta}_n(u^\varepsilon_{1,0}(x))dx, \end{aligned} \end{equation} for all $t$ in $(0,T)$ and $ \varepsilon < \frac{1}{n+1}$. For any subsequences $u^\varepsilon_i$ (still denoted by $u^\varepsilon_i$, $i=1,2,3)$, Lemma \ref{l2} shows that \begin{equation} \label{16} \begin{gathered} (u^\varepsilon_1,u^\varepsilon_2,u^\varepsilon_3) \to (u_1,u_2,u_3) \quad \text{a.e in } Q_T, \\ T_k(u^\varepsilon_i) \rightharpoonup T_k(u_i) \quad \text{weakly in } L^2(0,T;H^1_0(\Omega)),\; i=1,2,3, \\ \Theta_n(u^\varepsilon_i)\rightharpoonup \Theta_n (u_i) \quad \text{weakly in } L^2(0,T;H_0^1(\Omega)),\; i=1,2,3, \end{gathered} \end{equation} as $\varepsilon \to 0$ for any $k>0$ and $n\geq 1$. From (H2), \[ |A_i(T_k(u^\varepsilon_i),\nabla T_k(u^\varepsilon_i))| \leq C_k(x,t)+\beta_k |\nabla T_k(u^\varepsilon_i)|,\quad i=1,2,3, \] where $C_k \in L^2(Q_T)$. It shows that $A_i(T_k(u^\varepsilon_i), \nabla T_k(u^\varepsilon_i))$, $i=1,2,3$ is bounded in $L^2(Q_T)$. Therefore, \begin{equation} \label{17} A_i(T_k(u_i^\varepsilon),\nabla T_k(u^\varepsilon_i)) \rightharpoonup \eta_{i,k} \quad \text{weakly in } L^2(Q_T),\; i=1,2,3, \end{equation} as $\varepsilon \to 0$ where $\eta_{i,k} \in L^2(Q_T)$, $i=1,2,3$. From \eqref{10}, \eqref{10a} and \eqref{5} with Lemma \ref{l2} we obtain \begin{align*} \sum _{i=1}^3 \int_\Omega \tilde{T}_k(u^\varepsilon_i)(t) dx &\leq \int_{Q_t}(\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3) +\mu_1u_1^\varepsilon+f^\varepsilon)T_k(u^\varepsilon_1)\,dx\,ds\\ &\quad + \int_{Q_T}(\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)+(\mu_2+\varrho) u_2^\varepsilon+g^\varepsilon) T_k(u^\varepsilon_2) \,dx\,ds \\ &\quad + \int_{Q_t} (\varrho u_2^\varepsilon+ \mu_3u_3^\varepsilon+ h^\varepsilon)T_k(u^\varepsilon_3)\,dx\,ds + \sum_{i=1}^3 \int_\Omega \tilde{T}_k(u^\varepsilon_{i,0}(x))dx. \end{align*} Using the boundedness of the solutions $(u^\varepsilon_1,u^\varepsilon_2,u^\varepsilon_3)$, Lemma \ref{l2}, the properties of $\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)$, $ f^\varepsilon,g^\varepsilon, h^\varepsilon$ and Young's inequality, we obtain \begin{align*} &\sum _{i=1}^3 \int_\Omega \tilde{T}_k(u^\varepsilon_i)(t) dx\\ &\leq C_k + k(\|f^\varepsilon\|_{L^2(Q_T)} + \|g^\varepsilon\|_{L^2(Q_T)} + \|h^\varepsilon\|_{L^2(Q_T)} + \sum _{i=1}^3 \|u^\varepsilon_{i,0}(x)\|_{L^2(\Omega)}), \end{align*} where $C_k$ is a constant independent of $\varepsilon$. Taking $\lim \inf$ as $\varepsilon$ tends to zero in the above estimate and using \eqref{16} and Lemma \ref{l2}, we obtain \begin{align*} &\sum _{i=1}^3 \int_\Omega \tilde{T}_k(u_i)(t) dx\\ &\leq C_k + k(\|f\|_{L^1(Q_T)} + \|g\|_{L^1(Q_T)} + \|h\|_{L^1(Q_T)} + \sum _{i=1}^3 \|u_{i,0}(x)\|_{L^1(\Omega)}). \end{align*} By using the definition of $\tilde{T}_k(u_i)$, we deduce that \begin{equation} \label{18} \begin{aligned} \sum _{i=1}^3 k \int_\Omega |u_i(x,t)|dx &\leq C_k + \frac{k^2}{2} \operatorname{meas}(\Omega) +k (\|f\|_{L^1(Q_T)}+ \|g\|_{L^1(Q_T)} \\ &\quad + \|h\|_{L^1(Q_T)} + \sum _{i=1}^3 \|u_{i,0}(x)\|_{L^1(\Omega)}), \end{aligned} \end{equation} for almost all $t$ in $(0,T)$. This proves that $u_i \in L^\infty(0,T;L^1(\Omega))$, $i=1,2,3$. Now equation \eqref{15} with \eqref{16} proves that \begin{equation} \label{19} \begin{aligned} &\int_{Q_t} A_1(u^\varepsilon_1, \nabla u^\varepsilon_1) \nabla \Theta_n(u^\varepsilon_1) dx ds \\ &\leq \int_{Q_t}(\sigma(u_1,u_2,u_3)+\mu_1u_1+f)\Theta_n(u_1)\,dx\,ds + \int_\Omega \tilde{\Theta}_n(u_{1,0}(x))dx. \end{aligned} \end{equation} Using (H1), $\nabla \Theta_n(u^\varepsilon_1) = \chi_{\{n \leq |u^\varepsilon_1| \leq n+1\}}\nabla u_1^\varepsilon $ and the weak convergence in \eqref{16}, we obtain \begin{equation} \label{20} \begin{aligned} &\alpha _1 \int_{Q_t} |\nabla \Theta_n(u_1)|^2dx ds \\ & \leq \int_{Q_t}(\sigma(u_1,u_2,u_3)+\mu_1u_1+f)\Theta_n(u_1)\,dx\,ds + \int_\Omega \tilde{\Theta}_n(u_{1,0}(x))dx. \end{aligned} \end{equation} Since $\Theta_n(u_1) \to 0, $ as $n \to \infty$ shows that $\Theta_n(u_1) \to 0$, weakly in $L^2(0,T;H^1_0(\Omega))$. This leads to the right-hand side of each term of \eqref{20}, that is, \begin{gather*} \int_{Q_T}\sigma(u_1,u_2,u_3)\Theta_n(u_1)\,dx\,ds \to 0,\quad \int_{Q_T}\mu_1u_1\Theta_n(u_1)\,dx\,ds \to 0, \\ \int_{Q_T}f\Theta_n(u_1)\,dx\,ds \to 0 \end{gather*} as $n\to \infty$ and $\|\Theta_n(u_{1,0})\|_{L^1(\Omega)} \leq \|u_{1,0}\|_{L^1(\Omega)} $ implying that $ \int_\Omega \tilde{\Theta}_n(u_{1,0}) dx \to 0$ as $n\to \infty$ which follows from the Lebesgue convergence theorem. Hence, passing to the $\lim \inf$ in \eqref{19} and \eqref{20}, we obtain \[ \lim _{n\to \infty} \limsup _{\varepsilon\to 0} \int_0^t \int_{\{n \leq |u^\varepsilon_1| \leq n+1\}}A_1(u^\varepsilon_1, \nabla u^\varepsilon_1) \nabla \Theta _n(u^\varepsilon_1) \,dx\,ds \to 0, \] $\Theta_n(u_1) \to 0$, strongly in $L^2(0,T;H^1_0(\Omega))$ as $n\to \infty$. Similarly, by considering $\Theta_n(u^\varepsilon_2),\Theta_n(u^\varepsilon_3)$ respectively test functions in the second and third equations of \eqref{4} and by adopting the above type of arguments to prove that for $i=1,2,3$, \begin{equation} \label{21} \begin{gathered} \lim _{n\to \infty} \limsup _{\varepsilon\to 0} \int_0^t \int_{\{n \leq |u^\varepsilon_i| \leq n+1\}}A_i(u^\varepsilon_i, \nabla u^\varepsilon_i) \nabla \Theta _n(u^\varepsilon_i) \,dx\,ds \to 0, \\ \Theta_n(u_i) \to 0, \quad \text{strongly in $L^2(0,T;H^1_0(\Omega))$ as $n\to \infty$}. \end{gathered} \end{equation} This proves the desired result of the lemma. \end{proof} \begin{definition} \rm We define the time regularization of the function, for $i=1,2,3$, by \[ (T_k(u_i))_\gamma = \gamma \int _{-\infty} ^t e^{ \gamma (s-t) } T_k(\overline{u_i(x,s)}) ds, \] where $\overline {u_i(x,s)} = \begin{cases} u_i(x,s) & \text{if } s>0, \\ u_{i,0}(x)& \text{if } s<0. \end{cases}$ \end{definition} Let us consider the unique solution $ (T_k(u_i))_{\gamma} \in L^\infty (Q_T) \cap L^2(0,T; H_0^1(\Omega))$ of the monotone problem \begin{equation} \label{22} \begin{gathered} \frac{\partial}{\partial t} \left( T_k(u_i)\right)_{\gamma} + \gamma ((T_k(u_i))_\gamma - T_k(u_i))=0 \quad\text{in } Q_T ,\\ (T_k(u_i(x,0)))_\gamma = T_k(u_{i,0}(x)) \quad \text{in } \Omega, \end{gathered} \end{equation} for $\gamma >0$ and $k>0$. From (\ref{22}) and Lemma \ref{l2}, we have $ \frac{\partial}{\partial t} ( T_k(u_i))_{\gamma} \in L^2(0,T;H_0^1(\Omega))$. \begin{remark}\label{r4} \rm For $i=1,2,3$, we have $ (T_k(u_i))_\gamma \to T_k(u_i)$ a.e in $Q_T$, weak - $^*$ in $L^\infty(Q_T)$ and strongly in $L^2(0,T;H_0^1(\Omega))$ as $ \gamma \to \infty$ and also $$ \|(T_k(u_i))_\gamma\|_{L^\infty(Q_T)} \leq \max ( \|T_k(u_i)\|_{L^\infty(Q_T)}, \|T_k(u_{i,0})\|_{L^\infty(\Omega)}) \leq k $$ for any $\gamma >0$ and $k\geq 0$. \end{remark} \begin{lemma}\label{l5} Let $k \geq 0$ be fixed and $S$ be an increasing $ C^\infty (\mathbb{R})$ function such that $ S(z)=z$ for $ |z|\leq k$ with $\operatorname{supp}S'$ compact. Then, for $i=1,2,3$, \[ \lim _{n \to \infty} \limsup _ {\varepsilon \to 0} \int _0^T \int _{Q_t} \frac{\partial S(u_i^\varepsilon) }{\partial t} (T_k(u_i^\varepsilon) - (T_k(u_i))_\gamma) \,dx\,ds\,dt \geq 0. \] \end{lemma} The proof of the above lemma is similar to that of the lemma in \cite{blanjde}, and is omitted here. \begin{lemma}\label{l6} For $i=1,2,3$, $\eta_{i,k}$, as defined in \eqref{17}, the subsequences $u_i^\varepsilon$ (still denoted by $u^\varepsilon_i$) satisfy \begin{equation} \label{23} \begin{aligned} &\limsup _{\varepsilon\to 0} \int_0^T \int_{Q_t} A_i(T_k(u^\varepsilon_i), \nabla T_k(u^\varepsilon_i)) \nabla T_k(u^\varepsilon_i) \,dx\,ds\,dt\\ &\leq \int_0^T \int_{Q_t} \eta_{i,k} \nabla T_k(u_i) \,dx\,ds\,dt. \end{aligned} \end{equation} \end{lemma} \begin{proof} Let $S_n$ be a sequence of increasing $C^\infty(\mathbb R)$ functions such that \begin{equation} \label{24} \begin{gathered} S_n(z) = z, \quad \text{for } |z| \leq n, \\ \operatorname{supp} S_n' \subset [ -(n+1), (n+1)] ,\quad \|S_n'' \|_{L^\infty(\mathbb{R})} \leq 1. \end{gathered} \end{equation} Multiply the first equation of \eqref{4} by $ S_n'(u_1^\varepsilon)$ to obtain \begin{equation} \label{25} \begin{aligned} \frac{\partial S_n(u^\varepsilon_1)}{\partial t} &= \operatorname{div}\big(S_n'(u_1^\varepsilon) A_1(u^\varepsilon_1, \nabla u_1^\varepsilon)\big) - S_n''(u^\varepsilon_1)A_1(u^\varepsilon_1, \nabla u_1^\varepsilon)\nabla u^\varepsilon_1 \\ &\quad + \operatorname{div}\big( S_n'(u^\varepsilon_1)(\alpha_{1,1} \nabla u_1^\varepsilon + \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon)\big)\\ &\quad - S_n''(u^\varepsilon_1) (\alpha_{1,1} \nabla u_1^\varepsilon + \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon)\nabla u_1^\varepsilon \\ &\quad + S_n'(u_1^\varepsilon)(-\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)-\mu_1u_1^\varepsilon+f^\varepsilon). \end{aligned} \end{equation} From \eqref{25}, we understand that $\frac{\partial S_n(u^\varepsilon_1)}{\partial t} \in L^1(Q_T) + L^2(0,T;H^{-1}(\Omega))$. Similarly it holds for $i=2,3$, $ \frac{\partial S_n(u^\varepsilon_i)}{\partial t} \in L^1(Q_T) + L^2(0,T;H^{-1}(\Omega))$. For fixed $k>0, \gamma >0$, and $\varepsilon>0$, we set \begin{eqnarray}\label{26} W^\varepsilon_{i,\gamma} = T_k(u^\varepsilon_i)- (T_k(u_i))_\gamma, \quad i=1,2,3. \end{eqnarray} Multiplying \eqref{25} by $W^\varepsilon_{1,\gamma}$ and integrating over $Q_t \times (0,T)$, we obtain \begin{equation}\label{27} \int_Q \frac{\partial S_n(u^\varepsilon_1)}{\partial t}W^\varepsilon_{1,\gamma}\,dx\,ds\,dt = I_1 +I_2+I_3+I_4 + I_5, \end{equation} where \begin{gather*} I_1 = - \int_Q S_n'(u_1^\varepsilon) A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla W^\varepsilon_{1,\gamma}\,dx\,ds\,dt, \\ I_2 = - \int_Q S_n''(u^\varepsilon_1)A_1(u^\varepsilon_1, \nabla u_1^\varepsilon)\nabla u^\varepsilon_1 W^\varepsilon_{1,\gamma}\,dx\,ds\, dt, \\ I_3 = - \int_Q S_n'(u^\varepsilon_1)(\alpha_{1,1} \nabla u_1^\varepsilon + \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon)\nabla W^\varepsilon_{1,\gamma}\,dx\,ds\,dt , \\ I_4 = -\int_Q S_n''(u^\varepsilon_1) (\alpha_{1,1} \nabla u_1^\varepsilon + \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon)\nabla u_1^\varepsilon W^\varepsilon_{1,\gamma}\,dx\,ds\,dt, \\ I_5 = \int_Q S_n'(u_1^\varepsilon)(-\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)-\mu_1u_1^\varepsilon +f^\varepsilon)W^\varepsilon_{1,\gamma}\,dx\,ds\,dt , \end{gather*} where $Q=Q_t \times (0,T)$. By \eqref{26}, for fixed $\gamma >0, W^\varepsilon_{i,\gamma} \rightharpoonup T_k(u_i) - (T_k(u_i))_\gamma$, $i=1,2,3$, weakly in $L^2(0,T;H_0^1(\Omega))$ as $\varepsilon \to 0$. By Remark \ref{r4}, we conclude that $\|W^\varepsilon_{i,\gamma}\|_{L^\infty(Q_T)} \leq 2k$, for any $\varepsilon>0$ and $\gamma >0$. This boundedness of $W^\varepsilon_{1,\gamma}$ shows that for fixed $\gamma >0, W^\varepsilon_{i,\gamma} \rightharpoonup T_k(u_i) - (T_k(u_i))_\gamma$, $i=1,2,3$, a.e in $Q_T$ and $L^\infty (Q_T)$ weak$^*$ as $\varepsilon \to 0$. By the definition of $S_n$, we have $\operatorname{supp}S_n''\subset[-(n+1), (n+1)] \cup [n,(n+1)]$, for any $n\geq 1$. As a consequence \[ |I_2| \leq T \|S_n''\|_{L^\infty(\mathbb{R})} \|W^\varepsilon_{1,\gamma}\|_{L^\infty(Q_T)} \int _{\{ (x,t); n \leq |u^\varepsilon_1| \leq n+1\}}A_1(u^\varepsilon_1, \nabla u^\varepsilon_1) \nabla u^\varepsilon_1 \,dx\,ds, \] for any $n\geq1,\varepsilon \leq \frac{1}{n+1}$ and $\gamma >0$. From $\|W^\varepsilon_{1,\gamma}\|_{L^\infty(Q_T)} \leq 2k$ and \eqref{24}, we easily obtain \[ \limsup _{\gamma \to \infty}\limsup _{\varepsilon \to 0} | I_2| \leq c \limsup _{\varepsilon \to 0} \int_{\{ (x,t); n \leq |u^\varepsilon_1| \leq n+1\}}A_1(u^\varepsilon_1, \nabla u^\varepsilon_1) \nabla u^\varepsilon_1 \,dx\,ds, \] for any $n \geq 1$, where the constant $c$ depends on $T$ and $k$. Hence, by Lemma \ref{l4}, we achieve that \begin{eqnarray}\label{28} \lim _{n \to \infty} \limsup _{\gamma \to \infty}\limsup _{\varepsilon \to 0} \int_Q S_n''(u^\varepsilon_1)A_1(u^\varepsilon_1, \nabla u_1^\varepsilon)\nabla u^\varepsilon_1 W^\varepsilon_{1,\gamma}\,dx\,ds\,dt=0. \end{eqnarray} For some $n\geq 1$, $I_3$ can be rewritten as \begin{align*} I_3&= \int_Q S_n' (u^\varepsilon_1) \Big((k_1 F^+_\varepsilon(T_{n+1}(u^\varepsilon_1)) +F^+_\varepsilon(T_{n+1}(u^\varepsilon_2))+F^+_\varepsilon(T_{n+1}(u^\varepsilon_3)) \nabla T_{n+1}(u^\varepsilon_1) \\ &\quad + F^+_\varepsilon(T_{n+1}(u^\varepsilon_1)) \nabla T_{n+1}(u^\varepsilon_2) + F^+_\varepsilon(T_{n+1}(u^\varepsilon_1)) \nabla T_{n+1}(u^\varepsilon_3)\Big) \nabla W^\varepsilon_{1,\gamma}\,dx\,ds\,dt, \end{align*} a.e in $Q_T$. Since $\operatorname{supp}S_n' \subset[-(n+1), (n+1)]$, for $i=1,2,3$, the definition of $F^+_\varepsilon(u^\varepsilon_i)$ and the results \eqref{16} lead to $S' _n(u^\varepsilon_1) F^+_\varepsilon(T_{n+1}(u^\varepsilon_i)) \to S_n' (u_1) T_{n+1}(u_i)$ a.e in $Q_T$ and in $L^\infty(Q_T)$ weak$^*$ as $\varepsilon \to 0$. This proves \begin{align*} &\lim _{\varepsilon \to 0} I_3 \\ &= \int_Q S_n'(u_1)\Big((k_1 T_{n+1}(u_1) +T_{n+1}(u_2) +T_{n+1}(u_3) )\nabla T_{n+1}(u_1) \\ &\quad + T_{n+1}(u_1) \nabla T_{n+1}(u_2) + T_{n+1}(u_1) \nabla T_{n+1}(u_3)\Big) (\nabla T_k(u_1) - \nabla (T_k(u_1))_\gamma ) \,dx\, ds\,dt, \end{align*} for any $\gamma >0$. By using the Remark \ref{r4}, this leads to \begin{equation} \label{29} \lim_{\gamma \to \infty} \lim_{\varepsilon \to 0} \int_Q S_n'(u^\varepsilon_1)(\alpha_{1,1} \nabla u_1^\varepsilon + \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon) \nabla W^\varepsilon_{1,\gamma}\,dx\,ds\,dt =0. \end{equation} Similarly, we can show that \begin{equation} \label{30} \begin{gathered} \lim_{\gamma \to \infty} \lim_{\varepsilon \to 0}\int_Q S_n''(u^\varepsilon_1) (\alpha_{1,1} \nabla u_1^\varepsilon + \alpha_{1,2}\nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon)\nabla u_1^\varepsilon W^\varepsilon_{1,\gamma}\,dx\,ds\,dt =0, \\ \lim_{\gamma \to \infty} \lim_{\varepsilon \to 0}\int_Q S_n'(u_1^\varepsilon) (\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)+\mu_1u_1^\varepsilon)W^\varepsilon_{1,\gamma}\,dx\,ds\,dt=0. \end{gathered} \end{equation} Since $f S_n' (u_1) \in L^1(Q_T)$, \eqref{16} and Remark \ref{r4} lead to \begin{equation} \label{31} \lim_{\gamma \to \infty} \lim_{\varepsilon \to 0} \int_Q f^\varepsilon S_n'(u_1^\varepsilon) W^\varepsilon_{1,\gamma} \,dx\,ds\,dt =0. \end{equation} Consequently, from Lemma \ref{l5} and the definition of $W^\varepsilon_{1,\gamma}$, we have \begin{equation} \label{32} \lim_{\gamma \to \infty} \lim_{\varepsilon \to 0} \int_Q \frac{\partial S_n(u^\varepsilon_1)}{\partial t}W^\varepsilon_{1,\gamma}\,dx\,ds\,dt \geq 0 \quad \text{for any } n\geq k. \end{equation} It is easy to understand that, from \eqref{28}-\eqref{31} along with \eqref{27} and \eqref{32}, we obtain \begin{equation} \label{33} \lim_{\gamma \to \infty} \lim_{\varepsilon \to 0} I_1 = \lim_{\gamma \to \infty} \lim_{\varepsilon \to 0} \int_Q S_n'(u_1^\varepsilon) A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla W^\varepsilon_{1,\gamma}\,dx\,ds\,dt \leq 0. \end{equation} Since \[ S_n'(u_1^\varepsilon) A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla T_k(u_1^\varepsilon) = A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla T_k(u^\varepsilon_1) \] for $k\leq \frac{1}{\varepsilon}$ and $k\leq n$ because of the definition of $S_n$ and from \eqref{33}, we obtain \begin{align*} &\limsup_{\varepsilon \to 0} \int_Q A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla T_k(u^\varepsilon_1) \,dx\,ds\,dt \\ &\leq \lim_{\gamma \to \infty} \lim_{\varepsilon \to 0} \int_Q S_n'(u_1^\varepsilon) A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla (T_k(u_1^\varepsilon))_\gamma \,dx\,ds\,dt \quad \text{for } k \leq n. \end{align*} For $\varepsilon \leq 1/(n+1)$, we know that \[ S'_n(u^\varepsilon_1) A_1(u_1^\varepsilon, \nabla u_1^\varepsilon) = S' _n(u_1^\varepsilon) A_1(T_{n+1}(u^\varepsilon_1), \nabla T_{n+1}(u^\varepsilon_1)) \] a.e in $Q$. Due to \eqref{17}, we have $S' _n(u_1^\varepsilon) A_1(T_{n+1}(u^\varepsilon_1), \nabla T_{n+1}(u^\varepsilon_1)) \rightharpoonup S_n'(u_1) \eta_{1,n+1}$ weakly in $L^2(Q_T)$ as $\varepsilon \to 0$. This help us to prove that, for any $n \leq 1$, \begin{equation} \label{34} \begin{aligned} &\lim _{\gamma \to \infty} \lim _{\varepsilon \to 0} \int_Q S_n'(u_1^\varepsilon) A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla (T_k(u_1^\varepsilon))_\gamma \,dx\,ds\,dt \\ & =\int _Q S'_n(u_1) \eta_{1,n+1} \nabla T_k(u_1) \,dx\,ds\,dt \\ &= \int _Q \eta_{1,n+1} \nabla T_k(u_1) \,dx\,ds\,dt \,. \end{aligned} \end{equation} For any $k \leq n$, we have $$ A_1(T_{n+1}(u^\varepsilon_1), \nabla T_{n+1}(u^\varepsilon_1)) _{\chi(|u_1^\varepsilon| \leq k)} = A_1(T_k(u^\varepsilon_1), \nabla T_k(u^\varepsilon_1))_{\chi(|u_1^\varepsilon| \leq k)} \quad \text{a.e in } Q_T. $$ The above equation with \eqref{16} and \eqref{17} implies that $ \eta_{{1, n+1}_{\chi\{|u_1^\varepsilon| \leq k\}}} = \eta_{{1,k}_{\chi\{|u_1^\varepsilon| \leq k\}}}$ a.e in $ Q_T - \{ |u^\varepsilon_1|=k\}$ for $k\leq n$ as $\varepsilon \to 0$. Therefore, \eqref{34} becomes \begin{align*} \limsup _{\varepsilon \to 0} \int_Q A_1(T_k(u_1^\varepsilon), \nabla T_k(u^\varepsilon_1)) \nabla T_k(u^\varepsilon_1) \,dx\,dsdt = \int _Q \eta_{1,k} \nabla T_k(u_1)\,dx\,ds\,dt. \end{align*} Similar arguments as we used to obtain the previous equation lead to, for $i=1,2,3$, \begin{align*} \limsup _{\varepsilon \to 0} \int_Q A_i(T_k(u_i^\varepsilon), \nabla T_k(u^\varepsilon_i)) \nabla T_k(u^\varepsilon_i) \,dx\,ds\,dt = \int _Q \eta_{i,k} \nabla T_k(u_i)\,dx\,ds\,dt. \end{align*} This completes the proof of the lemma. \end{proof} \begin{lemma}\label{l7} For any $k>0$ and $i=1,2,3$, we have \begin{align*} &\limsup_{\varepsilon \to 0} \int _Q \big[A_i(T_k(u^\varepsilon_i),\nabla T_k(u^\varepsilon_i)) - A_i(T_k(u^\varepsilon_i), \nabla T_k(u_i))\big] \\ &\times \big[\nabla T_k(u^\varepsilon_i) - \nabla T_k(u_i)\big] \,dx\,ds\,dt=0. \end{align*} \end{lemma} \begin{proof} The monotone assumption (H3) shows that, for $i=1,2,3$, \begin{equation} \label{35} \begin{aligned} &\limsup_{\varepsilon \to 0} \int _Q \big[A_i(T_k(u^\varepsilon_i),\nabla T_k(u^\varepsilon_i)) - A_i(T_k(u^\varepsilon_i), \nabla T_k(u_i))\big]\\ &\times \big[\nabla T_k(u^\varepsilon_i)- \nabla T_k(u_i)\big] \,dx\,ds\,dt \geq 0 \end{aligned} \end{equation} for any $k\geq0$. We remark that (H2) and first result of \eqref{16} implies that $$ A_i(T_k(u_i^\varepsilon), \nabla T_k(u_i)) \to A_i(T_k(u_i), \nabla T_k(u_i)), $$ a.e in $Q_T$ as $\varepsilon \to 0$ and that $A_i(T_k(u_i^\varepsilon), \nabla T_k(u_i)) \leq C_k(x,t) + \beta_k |\nabla T_k(u_i)|$ a.e. in $Q_T$, uniformly with respect to $\varepsilon$. Then, by Lemma \ref{l6}, \eqref{16} and \eqref{17}, we have, for $i=1,2,3$, \begin{align*} &\limsup_{\varepsilon \to 0} \int _Q \big[A_i(T_k(u^\varepsilon_i),\nabla T_k(u^\varepsilon_i)) - A_i(T_k(u^\varepsilon_i), \nabla T_k(u_i))\big]\\ &\times \big[\nabla T_k(u^\varepsilon_i)- \nabla T_k(u_i)\big] \,dx\,dsdt=0. \end{align*} This completes the proof. \end{proof} \begin{lemma}\label{l8} For fixed $k\geq0$ and $i=1,2,3$, we have \begin{gather*} \eta_{i,k} = A_i(T_k(u_i), \nabla T_k(u_i)) \quad \text{a.e in } Q_T, \\ A_i(T_k(u^\varepsilon_i), \nabla T_k(u^\varepsilon_i)) \nabla T_k(u^\varepsilon_i) \rightharpoonup A_i(T_k(u_i), \nabla T_k(u_i)) \nabla T_k(u_i)\quad \text{weakly in } L^1(Q_T). \end{gather*} \end{lemma} \begin{proof} For any $k>0$ and $0 < \varepsilon < \frac{1}{k}$, from Lemma \ref{l6}, convergence \eqref{16} implies that, for $i=1,2,3$, $$ \limsup _{\varepsilon \to 0} \int_Q A_i(T_k(u^\varepsilon_i), \nabla T_k(u^\varepsilon_i)) \nabla T_k(u^\varepsilon_i) \,dx\, ds\, dt = \int_Q \eta_{i,k} \nabla T_k(u_i) \,dx\, ds\, dt. $$ Using Minty's type arguments and \eqref{16}, \eqref{17}, from the above equation we obtain $A_k(T_k(u_i), \nabla T_k(u_i))= \eta_{i,k}$, for $i=1,2,3$, and any $k \geq 0$. This proves the first result of the present lemma. For any $k \geq 0, T' 0, \phi_i \in C'(\bar{Q}_T)$ with $\phi _i =0, i=1,2,3 $ in $\Sigma_T$. Now, integrating the equation \eqref{2} over $Q_T$, we obtain \begin{equation} \label{37} \begin{aligned} & \int_0^T \langle u^\varepsilon_{1t}, T_k(u^\varepsilon_1- \phi_1)\rangle dt +\int _{Q_T} A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla T_k(u^\varepsilon_1- \phi_1)\,dx\,dt \\ & + \int_{Q_T} (\alpha_{1,1}\nabla u_1^\varepsilon + \alpha_{1,2} \nabla u_2^\varepsilon + \alpha_{1,3}\nabla u_3^\varepsilon) \nabla T_k(u^\varepsilon_1-\phi_1)\,dx\,dt \\ &= \int_{Q_T}(-\sigma(u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon) - \mu_1u_1^\varepsilon +f^\varepsilon)T_k(u^\varepsilon_1-\phi_1)\,dx\,dt, \\ & \int_0^T \langle u^\varepsilon_{2t}, T_k(u^\varepsilon_2- \phi_2)\rangle dt +\int _{Q_T} A_2(u^\varepsilon_2, \nabla u_2^\varepsilon) \nabla T_k(u^\varepsilon_2- \phi_2)\,dx\,dt \\ &+ \int_{Q_T} (\alpha_{2,1}\nabla u_1^\varepsilon + \alpha_{2,2} \nabla u_2^\varepsilon + \alpha_{2,3}\nabla u_3^\varepsilon) \nabla T_k(u^\varepsilon_2-\phi_2)\,dx\,dt \\ &= \int_{Q_T}(\sigma(u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon) - \varrho u_2^\varepsilon - \mu_2u_2^\varepsilon +g^\varepsilon)T_k(u^\varepsilon_2-\phi_2)\,dx\,dt, \\ & \int_0^T \langle u^\varepsilon_{3t}, T_k(u^\varepsilon_3- \phi_3)\rangle dt +\int _{Q_T} A_3(u^\varepsilon_3, \nabla u_3^\varepsilon) \nabla T_k(u^\varepsilon_3- \phi_3)\,dx\,dt \\ &+ \int_{Q_T} (\alpha_{3,1}\nabla u_1^\varepsilon + \alpha_{3,2} \nabla u_2^\varepsilon + \alpha_{3,3}\nabla u_3^\varepsilon) \nabla T_k(u^\varepsilon_3-\phi_3)\,dx\,dt \\ &= \int_{Q_T}(\varrho u_2^\varepsilon - \mu_3 u^\varepsilon_3 +h^\varepsilon) T_k(u^\varepsilon_3-\phi_3)\,dx\,dt. \end{aligned} \end{equation} Note that, if $L= k+ \|\phi_1\|_{L^\infty(Q_T)} $ and $ u^\varepsilon_{1t} = (u^\varepsilon_1-\phi_1)_t + \phi_{1t}$, we have \begin{equation} \label{38} \begin{aligned} & \int _{Q_T} A_1(u^\varepsilon_1, \nabla u_1^\varepsilon) \nabla T_k(u^\varepsilon_1- \phi_1)\,dx\,dt \\ &= \int _{Q_T} A_1(T_L(u^\varepsilon_1), \nabla T_L(u_1^\varepsilon)) \nabla T_k(u^\varepsilon_1- \phi_1)\,dx\,dt. \\ & \int _0^T \langle u^\varepsilon_{1t}, T_k(u^\varepsilon_1- \phi_1)\rangle dt \\ &= \int _\Omega \tilde{T}_k(u^\varepsilon_1-\phi_1)(T)dx - \int_\Omega \tilde{T}_k(u^\varepsilon_1-\phi_1)(0)dx \\ &\quad + \int_0^T \langle \phi_{1t}, T_k(T_L(u^\varepsilon_1)-\phi_1)\rangle dt. \end{aligned} \end{equation} From \eqref{38}, the first equation of \eqref{37} can be re-written as \begin{align*} &\int _\Omega \tilde{T}_k(u^\varepsilon_1-\phi_1)(T)dx - \int_\Omega \tilde{T}_k(u^\varepsilon_1-\phi_1)(0)dx + \int_0^T \langle \phi_{1t}, T_k(T_L(u^\varepsilon_1)-\phi_1)\rangle dt \\ & +\int _{Q_T} A_1(T_L(u^\varepsilon_1), \nabla T_L(u_1^\varepsilon)) \nabla T_k(T_L(u^\varepsilon_1) - \phi_1)\,dx\,dt \\ & + \int_{Q_T} (\alpha_{1,1}\nabla T_L(u_1^\varepsilon) + \alpha_{1,2} \nabla T_L(u_2^\varepsilon ) + \alpha_{1,3}\nabla T_L(u_3^\varepsilon)) \nabla T_k(T_L(u^\varepsilon_1)-\phi_1)\,dx\,dt \\ &= \int_{Q_T}(-\sigma(u_1^\varepsilon,u_2^\varepsilon,u_3^\varepsilon) - \mu_1u_1^\varepsilon +f^\varepsilon)T_k(u^\varepsilon_1-\phi_1)\,dx\,dt. \end{align*} Using the fact that $\tilde{T}_k$ is Lipschitz continuous, (H7) and the results \eqref{16}, we obtain \begin{gather*} \int_\Omega \tilde{T}_k(u^\varepsilon_1 - \phi_1)(T)dx \to \int_\Omega \tilde{T}_k(u_1-\phi_1)(T) dx , \\ \int_\Omega \tilde{T}_k(u^\varepsilon_1 -\phi_1)(0)dx \to \int_\Omega \tilde{T}_k(u_1 -\phi_1)(0) dx, \quad \text{as }\varepsilon \to 0. \end{gather*} Finally, by considering the strong convergence of $f^\varepsilon$, \eqref{16}, \eqref{17}, and the definition of $\sigma(u_1^\varepsilon,u^\varepsilon_2,u^\varepsilon_3)$ and the Lemma \ref{l8}, we obtain the following result by extracting the limit as $\varepsilon \to 0$ \begin{equation} \label{39} \begin{aligned} &\int _\Omega \tilde{T}_k(u_1-\phi_1)(T)dx - \int_\Omega \tilde{T}_k(u_1-\phi_1)(0)dx + \int_0^T \langle \phi_{1t}, T_k(u_1-\phi_1)\rangle dt \\ &+\int _{Q_T} A_1(u_1, \nabla u_1) \nabla T_k(u_1- \phi_1)\,dx\,dt \\ &+ \int_{Q_T} ((k_1u_1+u_2+u_3)\nabla u_1 + u_1 \nabla u_2 + u_1\nabla u_3) \nabla T_k(u_1-\phi_1)\,dx\,dt \\ &= \int_{Q_T}(-\sigma(u_1,u_2,u_3) - \mu_1u_1 +f)T_k(u_1-\phi_1)\,dx\,dt. \end{aligned} \end{equation} Similar arguments on the other two equations of the \eqref{37} lead to %\label{40} \begin{align*} &\int_\Omega \tilde{T}_k(u_2-\phi_2) (T) dx - \int_\Omega \tilde{T}_k(u_2-\phi_2)(0) dx + \int_0^T \langle \phi_{2t}, T_k(u_2-\phi_2)\rangle dt \\ &+ \int _{Q_T} A_2(u_2,\nabla u_2) \nabla T_k(u_2-\phi_2) \,dx\,dt + \int_{Q_T}((u_1+k_2u_2+u_3)\nabla u_2 + u_2\nabla u_1 \\ &+ u_2 \nabla u_3)\nabla T_k(u_2-\phi_2) \,dx\,dt \\ &= \int_{Q_T} (\sigma(u_1,u_2,u_3)- \varrho u_2 - \mu_2u_2 +g) T_k(u_2-\phi_2) \,dx\,dt, \\ & \int_\Omega \tilde{T}_k(u_3-\phi_3) (T) dx - \int_\Omega \tilde{T}_k(u_3-\phi_3)(0) dx + \int_0^T \langle \phi_{3t}, T_k(u_3-\phi_3)\rangle dt \\ &+ \int _{Q_T} A_3(u_3,\nabla u_3) \nabla T_k(u_3-\phi_3) \,dx\,dt + \int_{Q_T}((u_1+u_2+k_3u_3)\nabla u_3 + u_3\nabla u_1 \\ &+ u_3 \nabla u_2)\nabla T_k(u_3-\phi_3) \,dx\,dt \\ &= \int_{Q_T} (\varrho u_2 - \mu_3u_3 +h) T_k(u_3-\phi_3) \,dx\,dt, \end{align*} for all $k>0$ and for $i=1,2,3, \phi_i \in C' (\bar{Q}_T)$ with $\phi_i =0$ in $\Sigma_T$ . This completes the existence of entropy solutions of the system \eqref{1}. \end{proof} \subsection*{Acknowledgments} The authors want to thank the anonymous referees for their useful comments and suggestion, which led to the improvement in the quality of this article. The first author is thankful to the the University Grants Commission (UGC), New Delhi, India for awarding the Basic Scientific Research Fellowship. \begin{thebibliography}{100} \bibitem{benailan} P. B\'enilan, L. Boccardo, T. Gallout, R. Gariepy, M. Pierre, J. L. Vazquez; \emph{An $L^1$-theory of existence and uniqueness of solutions of nonlinear elliptic equations}, Ann. Scuola Norm. Sup. Pisa Cl. Sci., 22(1995), 241-273. \bibitem{bencross} B. Andreianov, M. Bendahmane, R. Ruiz-Baier; \emph{Analysis of finite volume method for a cross-diffusion model in population dynamics}, Math. Models Meth. Appl. Sci., 21(2011), 307-344. \bibitem{bendaaam} M. Bendahmane, M. Saad; \emph{Mathematical analysis and pattern formation for a partial immune system modeling the spread of an epidemic disease}, Acta Appl. Math., 115(2011), 17-42. \bibitem{benjee} M. Bendahmane and M. Langlais; \emph{A reaction-diffusion system with cross-diffusion modeling the spread of an epidemic disease}, J. Evol. Equ., 10(2010), 883-904. \bibitem{benjmpa} M. Bendahmane, T. Lepoutre, A. Marrocco, B. Perthame; \emph{Conservative cross diffusions and pattern formation through relaxation}, J. Math. Pures Appl., 92(2009), 651-667. \bibitem{bencpaa} M. Bendahmane, K. H. Karlsen; \emph{Renormalized solutions of anisotropic reaction-diffusion-advection systems with $L^1$ data modelling the propagation of an epidemic disease}, Commun. Pure Appl. Anal., 5(2006), 733-762. \bibitem{bennhm} M. Bendahmane and K.H. Karlsen; \emph{ Analysis of a class of degenerate reaction-diffusion systems and the bidomain model of cardiac tissue}, Netw. Heterog. Media, 1(2006), 185-218. \bibitem{bensiam} M. Bendahmane, K. H. Karlsen; \emph{Renormalized entropy solutions for quasi-linear anisotropic degenerate parabolic equations}, SIAM J. Math. Anal., 36(2004), 405-422. \bibitem{benjmaa} M. Bendahmane, M. Saad; \emph{A predator-prey system with $L^1$ data}, J. Math. Anal. Appl., 277(2003), 272-292. \bibitem{benade} M. Bendahmane, M. Langlais, M. Saad; \emph{Existence of solutions for reaction-diffusion system with $L^1$ data}, Adv. Differential Equations, 7(2002), 743-768. \bibitem{blanjde} D. Blanchard, F. Murat, H. Redwane; \emph{Existence and uniqueness of a renormalized solution for a fairly general class of nonlinear parabolic problems}, J. Differential Equations, 177(2001), 331-374. \bibitem{bocgjfa} L. Boccardo, T. Gallouet; \emph{On some nonlinear elliptic and parabolic equations involving measure data}, J. Funct. Anal., 87(1989), 149-169. \bibitem{bocjfa} L. Boccardo, A. Dall\'Aglio, T. Gallouet,, L. Orsina; \emph{Nonlinear parabolic equations with measure data}, J. Funct. Anal., 147(1997), 237-258. \bibitem{bocaa} L. Boccardo, L. Orsina; \emph{Existence results for Dirichlet problems in $L^1$ via Minty's lemma}, Appl. Anal., 76(2000), 309-317. \bibitem{chen} L. Chen, A. Jungel; \emph{Analysis of a parabolic cross-diffusion population model without self-diffusion}, J. Differential Equations, 224(2006), 39-59. \bibitem{del} M. Delgado, M. Montenegro, A. Suarez; \emph{A Lotka-Volterra symbiotic model with cross-diffusion}, J. Differential Equations, 246(2009), 2131-2149. \bibitem{lion} R. J. DiPerna, P. L. Lions; \emph{On the Cauchy problem for Boltzmann equations: Global existence and weak stability}, Ann. Math., 130(1989), 321-366. \bibitem{dro} J. Droniou, A. Prignet; \emph{Equivalence between entropy and renormalized solutions for parabolic equations with smooth measure data}, Nonlinear Differ. Equ. Appl., 14(2007), 181-205. \bibitem{evans} L. C. Evans; \emph{Partial Differential Equations}, AMS, Providence, 2002. \bibitem{fitana} W. E. Fitzgibbon, M. Langlais; \emph{Diffusive SEIR models with logistic population control}, Commun. Appl. Nonlinear Anal., 4(1997), 1-16. \bibitem{fitjam} W. E. Fitzgibbon, M. Langlais, J. J. Morgan; \emph{Eventually uniform bounds for a class of quasi positive balanced reaction diffusion systems}, Japan J. Indust. Appl. Math., 16(1999), 225-241. \bibitem{fitsiam} W. E. Fitzgibbon, M. Langlais, J. J. Morgan; \emph{A mathematical model of the spread of Feline Leukemia virus through a highly heterogeneous domain}, SIAM J. Math. Analysis, 33(2001), 570-588. \bibitem{fitfields} W. E. Fitzgibbon, M. Langlais, J.J. Morgan; \emph{Modeling the spread of Feline Leukemia in heterogeneous habitats}, Fields Inst. Commun., 29(2001), 133-146. \bibitem{gal} G. Galiano; \emph{On a cross-diffusion population model deduced from mutation and splitting of a single species}, Comput. Math. Appl., 64(2012), 1927-1936. \bibitem{horst} D. Horstmann; \emph{Remarks on some Lotka-Volterra type cross-diffusion models}, Nonlinear Anal. RWA, 8(2007), 90-117. \bibitem{ko} W. Ko, K. Ryu; \emph{On a predator-prey system with cross diffusion representing the tendency of predators in the presence of prey species}, J. Math. Anal. Appl., 341(2008), 1133-1142. \bibitem{ruiz} R. Ruiz-Baier, C. Tian; \emph{Mathematical analysis and numerical simulation of pattern formation under cross-diffusion}, Nonlinear Anal. RWA, 14(2013), 601-612. \bibitem{lsgaam} L. Shangerganesh, K. Balachandran; \emph{Existence and uniqueness of solutions of predator-prey type model with mixed boundary conditions}, Acta Appl. Math., 116(2011), 71-86. \bibitem{lsgaa} L. Shangerganesh, N. Barani Balan, K. Balachandran; \emph{Weak-renormalized solutions for predator-prey system}, Appl. Anal., 92(2013), 441-459. \bibitem{lsgm2as} L. Shangerganesh, K. Balachandran; \emph{Solvability of reaction-diffusion model with variable exponents}, Math. Methods Appl. Sci., DOI: 10.1002/mma.2905, in press. \bibitem{wang} Y. Wang, J. Wang, L. Zhang; \emph{Cross diffusion-induced pattern in an SI model}, Appl. Math. Comput., 217(2010), 1965-1970. \bibitem{yxwang} Y. X. Wang, W. T. Li; \emph{Effects of cross-diffusion and heterogeneous environment on positive steady states of a prey-predator system}, Nonlinear Anal. RWA, 14(2013), 1235-1246. \bibitem{wu} Y. Wu; \emph{The instability of spiky steady states for a competing species model with cross diffusion}, J. Differential Equations, 213(2005), 289-340. \bibitem{zeng} X. Zeng, Z. Liu; \emph{Nonconstant positive steady states for a ratio-dependent predator-prey system with cross-diffusion}, Nonlinear Anal. RWA, 11(2010), 372-390. \bibitem{zhang} C. Zhang, S. Zhou; \emph{Renormalized and entropy solutions for nonlinear parabolic equations with variable exponents and $L^1$ data}, J. Differential Equations, 248(2010), 1376-1400. \end{thebibliography} \end{document}