\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
International Conference on Applications of Mathematics to Nonlinear Sciences,\newline
\emph{Electronic Journal of Differential Equations},
Conference 24 (2017), pp. 11--22.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2017 Texas State University.}
\vspace{8mm}}

\begin{document} \setcounter{page}{11}
\title[\hfilneg EJDE-2017/conf/24\hfil Degenerate coupled flows in porous media]
{Global weak solutions to degenerate coupled diffusion-convection-dispersion
processes and heat transport in porous media}

\author[M. Bene\v{s}, Luk\'a\v{s} Krupi\v{c}ka \hfil EJDE-2017/conf/24\hfilneg]
{Michal Bene\v{s}, Luk\'a\v{s} Krupi\v{c}ka}

\address{Michal Bene\v{s} \newline
Department of Mathematics,
Faculty of Civil Engineering,
Czech Technical University in Prague,
 Th\'akurova 7, 166 29 Prague 6, Czech Republic}
\email{michal.benes@cvut.cz}

\address{Luk\'a\v{s} Krupi\v{c}ka \newline
Department of Mathematics,
Faculty of Civil Engineering,
Czech Technical University in Prague,
Th\'akurova 7, 166 29 Prague 6, Czech Republic}
\email{lukas.krupicka@fsv.cvut.cz}

\thanks{Published November 15, 2017.}
\subjclass[2010]{5A05, 35D05, 35B65, 35B45, 35B50, 35K15, 35K40}
\keywords{Initial-boundary value problems for second-order parabolic systems;
\hfill\break\indent global solution, smoothness and regularity of solutions;
coupled transport processes; 
\hfill\break\indent porous media}

\begin{abstract}
In this contribution we prove the existence of weak solutions
to degenerate parabolic systems
arising from the coupled moisture movement,
transport of dissolved species and heat transfer
through partially saturated porous materials.
Physically motivated mixed Dirichlet-Neumann boundary conditions
and initial conditions are considered.
Existence of a global weak solution
of the problem is proved by means of semidiscretization in time
and by passing to the limit from discrete approximations.
Degeneration occurs in the nonlinear transport coefficients
which are not assumed to be bounded below and above by positive
constants. Degeneracies in all transport coefficients
are overcome by proving suitable a priori $L^{\infty}$-estimates for
the approximations of primary unknowns of the system.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\allowdisplaybreaks

\section{Introduction}

Let $\Omega$ be a bounded domain in $\mathbb{R}^2$,
$\Omega \in C^{0,1}$ and let $\Gamma_D$ and $\Gamma_N$ be open disjoint
subsets of $\partial\Omega$ (not necessarily connected) such that 
$\Gamma_D\neq\emptyset$
and the $\partial\Omega \backslash (\Gamma_D\cup\Gamma_N)$ is a
finite set. Let $T\in(0,\infty)$ be fixed throughout the paper, $I=(0,T)$
and $Q_T = \Omega\times I$ denotes the space-time cylinder,
$\Gamma_{DT} = \Gamma_D\times I$ and
$\Gamma_{NT}=\Gamma_N\times I$.

We shall study the following initial boundary value problem in $Q_{T}$,
\begin{gather}
\partial_t b(u)
=\nabla\cdot[a(\theta)\nabla u], \label{strong:eq1a}\\
\partial_t [b(u)w]
=\nabla\cdot[b(u)D_w(u)\nabla w]
+\nabla\cdot[ w a(\theta)\nabla u], \label{strong:eq1b}\\
\partial_t [ b(u) \theta + {\varrho} \theta ]
=\nabla \cdot [ \lambda(\theta,u) \nabla \theta ]
+\nabla\cdot[\theta a(\theta)\nabla u ],\label{strong:eq1c}
\end{gather}
with the mixed-type boundary conditions
\begin{gather}
u = 0,\quad w = 0,\quad \theta = 0 \quad \text{on } \Gamma_{DT},\label{strong:eq1f}
\\
\nabla u \cdot \mathbf{n} = 0,\quad \nabla w \cdot \mathbf{n} = 0,
\quad \nabla \theta \cdot \mathbf{n} =0 \quad \text{on } \Gamma_{NT}
\label{strong:eq1i}
\end{gather}
and the initial conditions
\begin{equation}
u(\cdot,0) = u_0,\quad w(\cdot,0) = w_0,\quad \theta(\cdot,0) = \theta_0
\quad \text{in } \Omega.\label{strong:eq1l}
\end{equation}

System \ref{strong:eq1a}--\ref{strong:eq1l} arises from the coupled
 moisture movement,
transport of dissolved species and heat transfer
through the porous system \cite{Bear,PinderGray}.
For simplicity,
the gravity terms and external sources are not included since they do not
 affect the analysis.
For specific applications we refer the reader to e.g. \cite{Ozbolt}.
Here $u : Q_{T} \to \mathbb{R}$, $w : Q_{T} \to \mathbb{R}$
and $\theta : Q_{T} \to \mathbb{R}$ are the unknown functions.
In particular, $u$ corresponds to the Kirchhoff transformation
of the matric potential \cite{AltLuckhaus1983}, $w$ represents
concentration of dissolved species
and $\theta$ represents the temperature of the porous system.
Further, $a:\mathbb{R}\to \mathbb{R}$,
$D_w:\mathbb{R}\to \mathbb{R}$,
$b:\mathbb{R}\to \mathbb{R}$,
$\lambda:\mathbb{R}^2\to \mathbb{R}$,
$u_{0} : \Omega \to \mathbb{R}$,
$w_{0} : \Omega \to \mathbb{R}$,
and $\theta_{0} : \Omega\to \mathbb{R}$
are given functions,
$\varrho$ is a real positive constant
and $\mathbf{n}$ is the outward unit normal vector.
In this paper we study the existence
of the weak solution to \eqref{strong:eq1a}--\eqref{strong:eq1l}.

Nowadays, description of heat, moisture or soluble/non-soluble
contaminant transport in concrete, soil or rock porous matrix is
frequently based on time dependent models.
Coupled transport processes (diffusion processes, heat conduction, moister flow,
contaminant transport or coupled flows through porous media) are typically
 associated with systems of strongly nonlinear degenerate parabolic partial
differential equations of type (written in terms of operators
${A}$, ${\Psi}$, ${F}$)
\begin{equation}\label{eq:par1}
 \partial_t {\Psi}(\mathbf{u}) -\nabla\cdot {A}(\mathbf{u},\nabla \mathbf{u})
= {F}(\mathbf{u}),
\end{equation}
where $\mathbf{u}$ stands for the unknown vector of state variables.
There is no complete theory for
such general problems. However,
some particular results assuming special structure of
operators ${A}$ and ${\Psi}$ and growth
conditions on ${F}$ can be found in
the literature, see \cite{Roubicek2005}.

Most theoretical results on parabolic systems exclude the case of non-symmetrical
parabolic parts \cite{AltLuckhaus1983,FiloKacur1995,Kacur1990a}.

Giaquinta and Modica \cite{GiaquintaModica1987} proved the local-in-time
solvability of quasilinear
diagonal parabolic systems with nonlinear boundary conditions
 (without assuming any growth condition), see also \cite{Weidemaier1991}.

The existence of weak solutions to more general non-diagonal systems
like \eqref{eq:par1} subject to mixed
boundary conditions has been proven in \cite{AltLuckhaus1983}. The authors
proved an existence result assuming the operator ${\Psi}$
to be only (weak) monotone and subgradient. This
result has been extended in \cite{FiloKacur1995}, where the authors presented
the local existence of the weak
solutions for the system with nonlinear Neumann boundary conditions
and under more general growth conditions on nonlinearities in $\mathbf{u}$.
These results, however, are not applicable if ${\Psi}$ does not take the
subgradient structure, which is typical of coupled transport models in porous media.
Thus, the analysis needs to exploit the specific structure of such problems.

The existence of a local-in-time strong solution
for moisture and heat transfer in multi-layer porous structures
modelling by the doubly nonlinear parabolic system
is proven in \cite{BenesZeman}.
In \cite{Vala2002}, the author
proved the existence of the solution to the purely diffusive
hygro-thermal model allowing non-symmetrical operators $\Psi$, but requiring
non-realistic symmetry in the elliptic part.
In \cite{degond,jungel2000}, the authors studied the existence, uniqueness
and regularity of coupled quasilinear equations modeling evolution of
fluid species influenced by thermal, electrical
and diffusive forces.
In \cite{LiSun2010,LiSunWang2010,LiSun2012}, the authors studied a model
of specific structure
of a heat and mass transfer arising from textile industry and
proved the global existence for one-dimensional problems in
\cite{LiSun2010,LiSunWang2010}
and three-dimensional problems in \cite{LiSun2012}.

In the present paper we extend our previous existence result for coupled
 heat and mass flows in porous media \cite{BenesKrupicka2016}
to more general problem (including the convection-dispersion equation)
modeling coupled moisture, solute and heat transport in porous media.
This leads to a fully nonlinear degenerate parabolic system
with natural (critical) growths and degeneracies in all transport coefficients.

The rest of this paper is organized as follows. In Section \ref{sec:prelim},
we introduce basic notation and suitable function spaces
and specify our assumptions on data
and coefficient functions in the problem.
In Section~\ref{sec:main_result}, we formulate the problem in the
variational sense and state the main result, the global-in-time
existence of the weak solution.
The main result is proved by an approximation procedure in
Section \ref{sec:proof_main}.
First we formulate the semi-discrete scheme and prove
the existence of its solution.
The crucial a priori estimates and uniform boundness
of time interpolants are proved in part~\ref{sec:estimates}.
Finally, we conclude that
the solutions of semi-discrete scheme converge and the limit
is the solution of the original problem (Subsection~\ref{subsec:limit}).

\begin{remark} \label{rmk1.1} \rm
The present analysis can be straightforwardly extended to a setting with
nonhomogeneous boundary conditions (see \cite{BenesKrupicka2016} for details).
Here we work with homogeneous boundary conditions, ignoring the gravity terms
and excluding external sources to simplify
the presentation and avoid unnecessary technicalities in the existence result.
\end{remark}

\section{Preliminaries} \label{sec:prelim}

\subsection{Notation and some properties of Sobolev spaces}
\label{notation}

Vectors and vector functions are denoted by boldface letters.
Throughout the paper, we will always use positive constants $C$,
$c$, $c_1$, $c_2$, $\dots$, which are not specified and which may
differ from line to line.
Throughout this paper we suppose
$s,q,s'\in [1,\infty]$, $s'$ denotes the conjugate exponent to $s>1$,
${1}/{s} + {1}/{s'} = 1$.
$L^s(\Omega)$ denotes the usual
Lebesgue space equipped with the norm $\|\cdot\|_{L^s(\Omega)}$ and
$W^{k,s}(\Omega)$, $k\geq 0$ ($k$ need not to be an integer, see
\cite{KufFucJoh1977}), denotes the usual
Sobolev-Slobodecki space with the norm $\|\cdot\|_{W^{k,s}(\Omega)}$.
We define
\[
W^{1,2}_{\Gamma_D}(\Omega):=\big\{
v \in W^{1,2}(\Omega): v \big|_{\Gamma_D} = 0 \big\}.
\]
By $E^*$ we denote the space of all continuous, linear forms on Banach space $E$
and by $\langle \cdot,\cdot \rangle$ we denote the
duality between $E$ and $E^*$.
By $L^s(I;E)$ we denote the Bochner space (see \cite{AdamsFournier1992}).
Therefore, $L^s(I;E)^*=L^{s'}(I;E^*)$.


\subsection{Structure and data properties}
\label{Structure and data properties}

We start by introducing our assumptions on functions in
\eqref{strong:eq1a}--\eqref{strong:eq1l}.
\begin{itemize}

\item[(i)]
$b \in C^{1}(\mathbb{R})$, $ 0 < b'(\xi) < b_* $ and
\begin{equation*}
 0 < b(\xi) \leq b_{2} < +\infty \quad
\forall \xi \in \mathbb{R} \quad (b_{2},b_* = {\rm const}).
\end{equation*}

\item[(ii)]
$a$, $D_w$ $\in C(\mathbb{R})$ and $\lambda$ $\in C(\mathbb{R}^2)$ such that
\begin{gather*}
0 < a(\xi) , \quad 0 < D_w(\xi) \quad \forall \xi \in \mathbb{R},\\
0 < \lambda(\xi,\zeta) \quad \forall \xi,\zeta \in \mathbb{R}.
\end{gather*}

\item[(iii)] (Initial data)
Assume
$u_0, w_0, \theta_0 \in L^{\infty}(\Omega)$,
such that
\begin{equation}
 -\infty < u_{1} < u_0 < 0 \qquad \text{a.e. in } \Omega \quad
(u_{1}={\rm const}).
\end{equation}
\end{itemize}

\subsection{Auxiliary results}

\begin{remark}[{\cite[Section 1.1]{AltLuckhaus1983}}] \label{rmk2.1} \rm
Let us note that {\rm (i)} implies that
there is a (strictly) convex $C^1$-function
$\Phi:\mathbb{R}\to \mathbb{R}$, $\Phi(0)=0$,
$\Phi'(0)=0$, such that
$b(z) - b(0) = \Phi'(z)$ for all $z \in \mathbb{R}$.
Introduce the Legendre transform
$$
B(z):=\int_{0}^{1}(b(z)-b(sz)) z \, {\rm d}s
=\int_{0}^{z}(b(z)-b(s))\, {\rm d}s .
$$
\end{remark}

Let us present some properties of $B$ \cite{AltLuckhaus1983}:
\begin{gather*}
B(z):=\int_{0}^{1}(b(z)-b(sz)) z \, {\rm d}s \geq 0 \quad \forall z \in \mathbb{R},\\
B(s) - B(r) \geq (b(s)-b(r))r \quad \forall r,s \in \mathbb{R}, \\
b(z)z - \Phi(z) + \Phi(0) = B(z) \leq b(z)z \quad \forall z \in \mathbb{R}.
\end{gather*}

\section{Main result}\label{sec:main_result}

The aim of this paper is to prove the existence of a weak solution to problem
\eqref{strong:eq1a}--\eqref{strong:eq1l}.
First we formulate our problem in a variational sense.

\begin{definition} \label{def_weak_solution}\rm
A weak solution of \eqref{strong:eq1a}--\eqref{strong:eq1l}
is a triplet $[u,w,\theta]$ such that
\begin{gather*}
u \in L^2(I;W_{\Gamma_D}^{1,2}(\Omega)) ,\quad
w \in L^2(I;W_{\Gamma_D}^{1,2}(\Omega)) \cap L^{\infty}({Q_T}),\\
\theta \in L^2(I;W_{\Gamma_D}^{1,2}(\Omega)) \cap L^{\infty}({Q_T}),
\end{gather*}
which satisfies
\begin{equation} \label{weak_form_01}
-\int_{Q_T} b(u) \partial_t\phi\,{\rm d}x {\rm d}t
+\int_{Q_T} {a(\theta){\nabla u}} \cdot\nabla\phi
\,{\rm d}x {\rm d}t
=\int_{\Omega} b(u_{0}) \phi(\mathbf{x},0) \,{\rm d}x
\end{equation}
for any $\phi \in L^2(I;W^{1,2}_{\Gamma_D}(\Omega))
\cap W^{1,1}(I;L^{\infty}(\Omega))$ with $\phi(\cdot,T)=0$;
\begin{equation} \label{weak_form_02}
\begin{aligned}
&-\int_{Q_T} b(u)w \partial_t\eta\,{\rm d}x {\rm d}t
+\int_{Q_T}b(u)D_w(u)\nabla w \cdot\nabla \eta \,{\rm d}x {\rm d}t\\
&+ \int_{Q_T}w {a(\theta){\nabla u}} \cdot\nabla\eta\,{\rm d}x {\rm d}t\\
&=\int_{\Omega} b(u_0)w_0 \eta(\mathbf{x},0) \,{\rm d}x
\end{aligned}
\end{equation}
for any $\eta \in L^2(I;W^{1,2}_{\Gamma_D}(\Omega))\cap W^{1,1}(I;L^{\infty}(\Omega))$
with $\eta(\cdot,T)=0$;
\begin{equation}\label{weak_form_03}
\begin{aligned}
&-\int_{Q_T}[ b(u)\theta + {\varrho} \theta ] \partial_t\psi \,{\rm d}x {\rm d}t
+\int_{Q_T} \lambda(\theta,u) \nabla\theta \cdot \nabla\psi \,{\rm d}x {\rm d}t\\
&+\int_{Q_T}\theta {a(\theta){\nabla u}}\cdot \nabla\psi\,{\rm d}x {\rm d}t \\
&=\int_{\Omega}[ b(u_{0}) \theta_{0} + {\varrho} \theta_{0} ] \psi(\mathbf{x},0)
 \,{\rm d}x
\end{aligned}
\end{equation}
for any $\psi \in L^2(I;W^{1,2}_{\Gamma_D}(\Omega))
\cap W^{1,1}(I;L^{\infty}(\Omega))$
with $\psi(\cdot,T)=0$.
\end{definition}

The main result of this paper reads as follows.

\begin{theorem}\label{main_result}
Let assumptions {\rm (i)--(iii)} be satisfied.
Then there exists at least one weak solution of the
system \eqref{strong:eq1a}--\eqref{strong:eq1l}.
\end{theorem}

To prove the main result of the paper we use the method
of semidiscretization in time by constructing temporal approximations
and limiting procedure.
The proof can be divided into three steps.
In the first step we approximate our problem by means of a semi-implicit time
 discretization
scheme (which preserve the pseudo-monotone structure of the discrete problem)
and prove the existence and $W^{1,s}(\Omega)$-regularity (with some $s>2$)
of temporal approximations.
In the second step we construct piecewise constant time interpolants and
derive suitable a priori estimates.
The key point is to establish $L^{\infty}$-estimates to overcome degeneracies
in transport coefficients.
Finally, in the third step we pass to the limit from discrete approximations.

\section{Proof of the main result}\label{sec:proof_main}

\subsection{Approximations}\label{sec:approximations}

Let us fix $p \in \mathbb{N}$ and set $\tau:= T/p$ (a time step).
Further, let us consider
$u^0_{p} := u_{0}$, $w^0_{p} := w_{0}$,
$\theta^{0}_{p} := \theta_{0}$ a.e. on $\Omega$.
We approximate our evolution problem by a semi-implicit time discretization scheme.
Then we define, in each time step $n=1,\dots,p$, a triplet
$[u^n_{p},w^n_{p},\theta^n_{p}]$ as a solution of the
following recurrence steady problem.


For a given triplet
$[u^{n-1}_{p},w^{n-1}_{p},\theta^{n-1}_{p}]$,
$n=1,\dots,p$,
$u^{n-1}_{p} \in L^{\infty}(\Omega)$,
$w^{n-1}_{p} \in L^{\infty}(\Omega)$,
$\theta^{n-1}_{p} \in L^{\infty}(\Omega)$,
find $[u^n_{p},w^n_{p},\theta^n_{p}]$,
such that
$u^n_{p} \in W_{\Gamma_D}^{1,s}(\Omega)$,
$w^n_{p} \in W_{\Gamma_D}^{1,s}(\Omega)$,
$\theta^n_{p} \in W_{\Gamma_D}^{1,s}(\Omega)$ with some $s>2$
and
\begin{equation} \label{approximate_problem_01}
\int_{\Omega}\frac{b(u_{p}^n) - b(u_{p}^{n-1})}{\tau} \phi
\,{\rm d}x
+\int_{\Omega} a(\theta_{p}^{n-1}) \nabla u_{p}^n
\cdot\nabla\phi \,{\rm d}x
=0
\end{equation}
for any $\phi \in {W_{\Gamma_D}^{1,2}(\Omega)}$;
\begin{equation}\label{approximate_problem_02}
\begin{aligned}
&\int_{\Omega}\frac{ b(u_{p}^n)w^n_{p} - b(u_{p}^{n-1})w_{p}^{n-1} }{\tau} \eta
\,{\rm d}x \\
&+\int_{\Omega}b(u_{p}^{n-1})D_w(u_{p}^{n-1})
\nabla w_{p}^n \cdot \nabla\eta \,{\rm d}x
+\int_{\Omega} w_{p}^n a(\theta_{p}^{n-1}) \nabla u_{p}^n
\cdot \nabla \eta \,{\rm d}x
=0
\end{aligned}
\end{equation}
for any $\eta \in {W_{\Gamma_D}^{1,2}(\Omega)}$;
\begin{equation}\label{approximate_problem_03}
\begin{aligned}
&\int_{\Omega} \frac{ b(u_{p}^n)\theta^n_{p}
 - b(u_{p}^{n-1})\theta_{p}^{n-1} }{\tau} \psi\,{\rm d}x
+\int_{\Omega}{\varrho}\frac{ \theta_{p}^n - \theta_{p}^{n-1} }{\tau} \psi
\,{\rm d}x \\
&+\int_{\Omega} \lambda(\theta_{p}^{n-1},u_{p}^{n-1}) \nabla
\theta_{p}^n \cdot \nabla\psi{{\rm d}\Omega}
+\int_{\Omega}\theta_{p}^n a(\theta_{p}^{n-1}) \nabla u_{p}^n
\cdot \nabla\psi {{\rm d}\Omega}
=0
\end{aligned}
\end{equation}
for any $\psi \in {W_{\Gamma_D}^{1,2}(\Omega)}$.

Next we show the existence of the solution to
\eqref{approximate_problem_01}--\eqref{approximate_problem_03}.

\begin{theorem} \label{thm:aprox}
Let $u^{n-1}_{p} \in L^{\infty}(\Omega)$,
$w^{n-1}_{p} \in L^{\infty}(\Omega)$,
$\theta^{n-1}_{p} \in L^{\infty}(\Omega)$
be given and the assumptions {\rm (i)--(iii)} be satisfied.
Then there exists $[u^n_{p},w^n_{p},\theta^n_{p}]$,
such that
$u^n_{p} \in W_{\Gamma_D}^{1,s}(\Omega)$,
$w^n_{p} \in W_{\Gamma_D}^{1,s}(\Omega)$
and
$\theta^n_{p} \in W_{\Gamma_D}^{1,s}(\Omega)$ with some $s>2$
satisfying
\eqref{approximate_problem_01}--\eqref{approximate_problem_03}.
\end{theorem}

\begin{proof}
The proof rests on the
$W^{1,p}$-regularity of elliptic problems presented in \cite{gallouet,groger1989}
and the embedding
$W_{\Gamma_D}^{1,s}(\Omega) \subset L^{\infty}(\Omega)$ if $s>2$
(recall that $\Omega$ is a bounded domain in $\mathbb{R}^2$).

The existence of $u^n_{p} \in W_{\Gamma_D}^{1,s}(\Omega)$ with some $s>2$
and $\theta^n_{p} \in W_{\Gamma_D}^{1,2}(\Omega)$, solutions to problems
\eqref{approximate_problem_01} and \eqref{approximate_problem_03}, respectively,
is proven in \cite{BenesKrupicka2016}.
The existence of $w^n_{p} \in W_{\Gamma_D}^{1,2}(\Omega)$,
the solution to \eqref{approximate_problem_02},
can be handled in the same way.

Now, with $w^n_{p} \in W_{\Gamma_D}^{1,2}(\Omega)$ in hand,
rewrite the equation \eqref{approximate_problem_02} in the form
(transferring the lower-order terms to the right hand side)
\begin{align*}
&\int_{\Omega} b(u_{p}^{n-1})D_w(u_{p}^{n-1})
\nabla w_{p}^n \cdot \nabla\eta \,{\rm d}x\\
&=-\int_{\Omega}\frac{ b(u_{p}^n)w^n_{p} - b(u_{p}^{n-1})w_{p}^{n-1} }{\tau} \eta
\,{\rm d}x
-\int_{\Omega}w_{p}^n a(\theta_{p}^{n-1}) \nabla u_{p}^n\cdot \nabla \eta
\,{\rm d}x.
\end{align*}
Since $u^{n-1}_{p} \in L^{\infty}(\Omega)$,
$u^n_{p} \in W_{\Gamma_D}^{1,s}(\Omega)$ with some $s>2$,
$w^{n-1}_{p} \in L^{\infty}(\Omega)$,
$\theta^{n-1}_{p} \in L^{\infty}(\Omega)$,
both integrals on the right hand side
make sense for any $\eta \in W_{\Gamma_D}^{1,r'}(\Omega)$, $r'=r/(r-1)$
with some $r>2$.
Now we are able to apply \cite[Theorem~4]{gallouet} to obtain
$w^n_{p} \in W_{\Gamma_D}^{1,s}(\Omega)$ with some $s>2$.
Analysis similar to the above implies that
$\theta^n_{p} \in W_{\Gamma_D}^{1,s}(\Omega)$ with some $s>2$.
\end{proof}


\subsection{A priori estimates} \label{sec:estimates}
In this part we prove some uniform estimates (with respect to $p$) for
the time interpolants of the solution.
In the following estimates, many different constants
will appear. For simplicity of notation,
$C$ represents generic constants which may change their numerical value
from one formula to another
but do not depend on $p$ and
the functions under consideration.

\subsubsection{Construction of temporal interpolants}

With the sequences $u^n_{p},w^n_{p},\theta^n_{p}$ constructed
in Section~\ref{sec:approximations}, we define
the piecewise constant interpolants
$\bar{\phi}_{p}(t) = \phi_{p}^{n}$ for $t \in ((n - 1)\tau, n\tau]$
and, in addition, we extend $\bar{\phi}_{p}$ for $t\leq 0$ by
$\bar{\phi}_{p}(t) = \phi_{0}$ for $t \in (-\tau, 0]$.
For a function $\varphi$ we often use the simplified notation
$\varphi := \varphi(t)$, $\varphi_\tau(t) := \varphi(t-\tau)$,
$\partial_t^{-\tau}\varphi(t) := \frac{\varphi(t) - \varphi(t-\tau)}{\tau}$,
$\partial_t^{\tau}\varphi(t) := \frac{\varphi(t+\tau) - \varphi(t)}{\tau}$.
Then,
following \eqref{approximate_problem_01}--\eqref{approximate_problem_03},
the piecewise constant time interpolants
$\bar{u}_{p} \in L^{\infty}(I;W_{\Gamma_D}^{1,s}(\Omega))$,
$\bar{w}_{p} \in L^{\infty}(I;W_{\Gamma_D}^{1,s}(\Omega))$
and
$\bar{\theta}_{p} \in L^{\infty}(I;W_{\Gamma_D}^{1,s}(\Omega))$ (with some $s>2$)
satisfy the equations
\begin{equation}
\label{eq:101}
\int_{\Omega}
\partial_t^{-\tau} {b}(\bar{u}_{p}(t)) \phi \,{\rm d}x
+
\int_{\Omega}
a(\bar{\theta}_{p}(t-\tau))\nabla \bar{u}_{p}(t) \cdot\nabla \phi
\,{\rm d}x
=0
\end{equation}
for any $\phi \in {W_{\Gamma_D}^{1,2}(\Omega)}$,
\begin{equation}\label{eq:102}
\begin{aligned}
&\int_{\Omega} \partial_t^{-\tau} [ {b}(\bar{u}_{p}(t))\bar{w}_{p}(t) ] \eta
\,{\rm d}x
+\int_{\Omega}
{b}(\bar{u}_{p}(t-\tau))D_w(\bar{u}_{p}(t-\tau)) \nabla \bar{w}_{p}(t) \cdot \nabla\eta
\,{\rm d}x\\
&+\int_{\Omega}\bar{w}_{p}(t) a(\bar{\theta}_{p}(t-\tau))\nabla \bar{u}_{p}(t)
\cdot \nabla\eta\,{\rm d}x
=0
\end{aligned}
\end{equation}
for any $\eta \in {W_{\Gamma_D}^{1,2}(\Omega)}$
and
\begin{equation}\label{eq:103}
\begin{aligned}
&\int_{\Omega}\partial_t^{-\tau} \left[ {b}(\bar{u}_{p}(t))\bar{\theta}_{p}(t)
+{\varrho}\bar{\theta}_{p}(t)\right] \psi\,{\rm d}x \\
&+\int_{\Omega}
\lambda(\bar{\theta}_{p}(t-\tau),\bar{u}_{p}(t-\tau)) \nabla
\bar{\theta}_{p}(t) \cdot \nabla\psi\,{\rm d}x\\
&+\int_{\Omega}\bar{\theta}_{p}(t)
a(\bar{\theta}_{p}(t-\tau))\nabla \bar{u}_{p}(t)
\cdot \nabla\psi\,{\rm d}x
=0
\end{aligned}
\end{equation}
for any $\psi \in {W_{\Gamma_D}^{1,2}(\Omega)}$.

\subsubsection{$L^{\infty}$-bound for $\bar{u}_{p}$, $\bar{w}_{p}$
and $\bar{\theta}_{p}$}

First we prove the $L^{\infty}$-estimate for $\bar{u}_{p}$.
Let us set
\begin{equation}
\phi := [ b(\bar{u}_{p}) - b(u_{1}) ]_{-}
=\begin{cases}
b(\bar{u}_{p}) - b(u_{1}) , & \bar{u}_{p} < u_{1} ,\\
0 , & \bar{u}_{p} \geq u_{1},
\end{cases} 
\end{equation}
as a test function in \eqref{eq:101}. Note that $\phi$ vanishes on $\Gamma_D$.
It is a simple matter to derive
\[
\frac{1}{2}\int_{\Omega} [ b(\bar{u}_{p}(t))-b(u_{1}) ]_-^2{\rm d}{x}
+\int_{Q_t}a(\bar{\theta}_{p}(s-\tau)) b'(\bar{u}_{p}(s)) |\nabla \bar{u}_{p}(s)|^2
\chi_{\left\{ \bar{u}_{p} < u_{1} \right\}}
{\rm d}x{\rm d}s
\leq 0
\]
for almost every $t\in I$.
Hence we conclude that the set $\left\{x \in \Omega: \bar{u}_{p}(x,t) < u_{1} \right\}$
has a measure zero for almost every $t\in I$.

Now setting
\begin{equation}
\phi = [ b(\bar{u}_{p})-b(0) ]_{+}
=\begin{cases}
b(\bar{u}_{p}) - b(0) , & \bar{u}_{p} > 0 ,\\
0 , & \bar{u}_{p} \leq 0,
\end{cases}
\end{equation}
we obtain, using similar arguments,
$$
\frac{1}{2}\int_{\Omega} [ b(\bar{u}_{p}) - b(0) ]_+^2{\rm d}{x} = 0
\quad \text{for almost every }t\in I.
$$
Hence the set $\left\{x \in \Omega: \bar{u}_{p}(x,t) > 0 \right\}$
has a measure zero for almost every $t\in I$.
Finally, combining the previous arguments, we deduce
\begin{equation}\label{est:uniform_bound_u_02}
\|\bar{u}_{p}\|_{L^{\infty}({Q_T})} \leq {C},
\end{equation}
where ${C}$ does not depend on $p$.

Now we prove a similar estimate for $\bar{w}_{p}$.
Let $\ell$ be an odd integer.
Using $\phi = [\ell/(\ell+1)] (\bar{w}_{p})^{\ell+1}$
as a test function in \eqref{eq:101}
and
$\eta = (\bar{w}_{p})^{\ell}$ in \eqref{eq:102}
and combining both equations we obtain
\begin{equation} \label{eq:501}
\begin{aligned}
&\frac{1}{\tau}\frac{1}{\ell+1}\int_{\Omega}
 {b}(\bar{u}_{p}(s))[\bar{w}_{p}(s)]^{\ell+1}\,{\rm d}x\\
&-\frac{1}{\tau}\frac{1}{\ell+1}
\int_{\Omega} {b}(\bar{u}_{p}(s-\tau))[\bar{w}_{p}(s-\tau)]^{\ell+1}
\,{\rm d}x\\
&+\frac{1}{\tau}
\frac{1}{\ell+1}
\int_{\Omega} {b}(\bar{u}_{p}(s-\tau))[\bar{w}_{p}(s-\tau)]^{\ell+1}
\,{\rm d}x\\
&+\frac{1}{\tau}\frac{\ell}{\ell+1}\int_{\Omega}
 {b}(\bar{u}_{p}(s-\tau))[\bar{w}_{p}(s)]^{\ell+1}\,{\rm d}x\\
&-\frac{1}{\tau}\int_{\Omega}
 {b}(\bar{u}_{p}(s-\tau))\bar{w}_{p}(s-\tau)[\bar{w}_{p}(s)]^{\ell}\,{\rm d}x\\
&+\int_{\Omega}
\ell [\bar{w}_{p}(s)]^{\ell-1}
{b}(\bar{u}_{p}(s-\tau))D_w(\bar{u}_{p}(s-\tau)) \nabla \bar{w}_{p}(s)
\cdot \nabla\bar{w}_{p}(s)\,{\rm d}x
=0.
\end{aligned}
\end{equation}
Applying the Young's inequality we can write for the term in the third line
\begin{equation}\label{eq:502}
\begin{aligned}
&\frac{1}{\tau}\int_{\Omega}
 {b}(\bar{u}_{p}(s-\tau))\bar{w}_{p}(s-\tau)[\bar{w}_{p}(s)]^{\ell}
\,{\rm d}x \\
&\leq \frac{1}{\tau}\frac{1}{\ell+1}\int_{\Omega}
 {b}(\bar{u}_{p}(s-\tau))[\bar{w}_{p}(s-\tau)]^{\ell+1}
\,{\rm d}x\\
&+\frac{1}{\tau}\frac{\ell}{\ell+1}\int_{\Omega}
 {b}(\bar{u}_{p}(s-\tau))[\bar{w}_{p}(s)]^{\ell+1}
\,{\rm d}x.
\end{aligned}
\end{equation}
Combining \eqref{eq:501} and \eqref{eq:502} we deduce
\begin{equation} \label{eq:503}
\begin{aligned}
&\frac{1}{\tau}\frac{1}{\ell+1}
\int_{\Omega} {b}(\bar{u}_{p}(s))[\bar{w}_{p}(s)]^{\ell+1}\,{\rm d}x\\
&-\frac{1}{\tau}\frac{1}{\ell+1}\int_{\Omega}
 {b}(\bar{u}_{p}(s-\tau))[\bar{w}_{p}(s-\tau)]^{\ell+1}\,{\rm d}x\\
&+\int_{\Omega}\ell [\bar{w}_{p}(s)]^{\ell-1}
{b}(\bar{u}_{p}(s-\tau))D_w(\bar{u}_{p}(s-\tau)) \nabla \bar{w}_{p}(s)
 \cdot \nabla\bar{w}_{p}(s)\,{\rm d}x
\leq 0.
\end{aligned}
\end{equation}
Now, integrating \eqref{eq:503} over $s$ from $0$ to $t$ we obtain
\begin{equation} \label{app:eq_bound_w_00}
\begin{aligned}
&\int_{\Omega}(\bar{w}_{p}(t))^{\ell+1}b(\bar{u}_{p}(t)){\rm d}x\\
&+\int_{\Omega_t}(\ell+1)\ell [\bar{w}_{p}(s)]^{\ell-1}
{b}(\bar{u}_{p}(s-\tau))D_w(\bar{u}_{p}(s-\tau))
|\nabla \bar{w}_{p}(s)|^2{\rm d}x{\rm d}s\\
& \leq\int_{\Omega}(w_0)^{\ell+1}b\left(u_0\right){\rm d}{x}.
\end{aligned}
\end{equation}
Note that the second integral in
\eqref{app:eq_bound_w_00} is nonnegative ($\ell$ is supposed to be the odd integer).
Moreover, from \eqref{app:eq_bound_w_00} and \eqref{est:uniform_bound_u_02}
it follows that
\begin{equation}\label{est:unform_bound_w_01}
\|\bar{w}_{p}\|_{L^{\infty}(0,T;L^{\ell+1}(\Omega))} \leq C,
\end{equation}
where the constant $C$ is independent of $\ell$ and $p$.
Now, let $\ell \to +\infty$ in \eqref{est:unform_bound_w_01}, we obtain
\begin{equation}\label{est:unform_bound_w_05}
\|\bar{w}_{p}\|_{L^{\infty}({Q_T})} \leq C.
\end{equation}
In the same manner we arrive at
the estimate for $\bar{\theta}_{p}$, i.e.
\begin{equation}\label{est:unform_bound_theta_05}
\|\bar{\theta}_{p}\|_{L^{\infty}({Q_T})} \leq C.
\end{equation}

\subsubsection{Energy estimates for $\bar{u}_{p}$, $\bar{w}_{p}$ and 
$\bar{\theta}_{p}$}
We test \eqref{eq:101} with $\phi = \bar{u}_{p}(t)$ and integrate 
\eqref{eq:101} over $t$
from $0$ to $s$. For the parabolic term we can write
\begin{equation}\label{est:301}
\int_0^s \int_{\Omega}
\partial_t^{-\tau} {b}(\bar{u}_{p}(t)) \bar{u}_{p}(t)
\,{\rm d}x{\rm d}t
\geq \frac{1}{\tau}\int_{s-\tau}^s
\int_{\Omega}{B}(\bar{u}_{p}(t)) - {B}( u_0 )
\,{\rm d}x{\rm d}t.
\end{equation}
Further, using \eqref{est:uniform_bound_u_02} and \eqref{est:301},
applying the usual estimates for the elliptic
part (see also \cite{AltLuckhaus1983}), we obtain the a priori estimate
\begin{equation}\label{energy_estimate_u}
\sup_{0\leq t \leq T} \int_{\Omega}{B}(\bar{u}_{p}(t)) {\rm d}x
+\int_0^T\int_{\Omega} |\nabla \bar{u}_{p}(t)|^2 {\rm d}x{\rm d}t
\leq {C}.
\end{equation}
Now it follows that there exists a function
$u \in L^2(I;W^{1,2}_{\Gamma_D}(\Omega))$
such that, along a selected
subsequence (letting $p \to \infty$), we have
$\bar{u}_{p}(t) \rightharpoonup u$ weakly in 
$L^2(I;W^{1,2}_{\Gamma_D}(\Omega))$.

Now we prove similar result for $\bar{w}_{p}(t) $.
Using $\eta(t) = 2 \bar{w}_{p}(t) $
as a test function in \eqref{eq:102} we obtain
\begin{equation} \label{est:_w_01}
\begin{aligned}
&\int_{\Omega}\partial_t^{-\tau}{b}(\bar{u}_{p}(t)) 2 \bar{w}_{p}(t)^2
\,{\rm d}x
+\int_{\Omega}\partial_t^{-\tau} \bar{w}_{p}(t) 2\bar{w}_{p}(t)
b(\bar{u}_{p}(t-\tau))\,{\rm d}x\\
&+2\int_{\Omega} {b}(\bar{u}_{p}(t-\tau))D_w(\bar{u}_{p}(t-\tau))
 \nabla \bar{w}_{p}(t)\cdot\nabla\bar{w}_{p}(t) \,{\rm d}x\\
&+\int_{\Omega} a(\bar{\theta}_{p}(t-\tau)) \nabla \bar{u}_{p}(t)
\cdot 2\bar{w}_{p}(t) \nabla \bar{w}_{p}(t)\,{\rm d}x
=0.
\end{aligned}
\end{equation}
One is allowed to use $\phi(t) = \bar{w}_{p}(t)^2$
as a test function in \eqref{eq:101} to obtain
\begin{equation}\label{est:_w_02}
\int_{\Omega} [\partial_t^{-\tau} {b}(\bar{u}_{p}(t))]
\bar{w}_{p}(t)^2\,{\rm d}x
+\int_{\Omega}a(\bar{\theta}_{p}(t-\tau)) \nabla \bar{u}_{p}(t)
\cdot\nabla \bar{w}_{p}(t)^2 \,{\rm d}x
=0.
\end{equation}
Combining \eqref{est:_w_01} and \eqref{est:_w_02} we deduce
\begin{equation} \label{est:w_100}
\begin{aligned}
&\int_{\Omega}\partial_t^{-\tau}\left[
\bar{w}_{p}(t) ^2b(\bar{u}_{p}(t))\right]\,{\rm d}x
+\int_{\Omega}\frac{1}{\tau}\left[\bar{w}_{p}(t)-\bar{w}_{p}(t-\tau)\right]^2
b(\bar{u}_{p}(t-\tau))\,{\rm d}x\\
&+2\int_{\Omega}{b}(\bar{u}_{p}(t-\tau))D_w(\bar{u}_{p}(t-\tau))
\nabla \bar{w}_{p}(t) \cdot \nabla \bar{w}_{p}(t)\,{\rm d}x
=0.
\end{aligned}
\end{equation}
In view of \eqref{est:uniform_bound_u_02} we have
\begin{equation}
b(\bar{u}_{p}(t)), \;
{b}(\bar{u}_{p}(t-\tau)), \;
D_w(\bar{u}_{p}(t-\tau)) > C \qquad \text{in } \Omega \times (-\tau,T).
\end{equation}
Recall that $C$ does not depend on $p$.
Now, integrating \eqref{est:w_100} with respect to time $t$ we obtain
\begin{equation*}
\sup_{0\leq t \leq T}
\int_{\Omega} |\bar{w}_{p}(t)|^2 {\rm d}\Omega
+\int_0^T \|\bar{w}_{p}(t)\|^2_{W^{1,2}_{\Gamma_D}(\Omega)} {\rm d}\Omega
\leq {C}.
\end{equation*}
From this we can write
\begin{equation}\label{est20_w_a}
\| \bar{w}_p \|_{L^2(I;W^{1,2}_{\Gamma_D}(\Omega))} \leq {C}.
\end{equation}
Similarly, we use $\psi(t) = 2\bar{\theta}_{p}(t)$
as a test function in {\eqref{eq:103}} to obtain
\begin{equation} \label{approximate_estimate_02}
\begin{aligned}
&\int_{\Omega}\partial_t^{-\tau} b(\bar{u}_{p}(t))2 \bar{\theta}_{p}(t)^2
\,{\rm d}x
+\int_{\Omega}\partial_t^{-\tau} \bar{\theta}_{p}(t) 2\bar{\theta}_{p}(t)
b(\bar{u}_{p}(t-\tau))\,{\rm d}x\\
&+2\int_{\Omega}\lambda(\bar{\theta}_{p}(t-\tau),\bar{u}_{p}(t-\tau))
\nabla \bar{\theta}_{p}(t)
\cdot\nabla \bar{\theta}_{p}(t)\,{\rm d}x\\
&+\int_{\Omega}a(\bar{\theta}_{p}(t-\tau)) \nabla \bar{u}_{p}(t)
\cdot 2 \bar{\theta}_{p}(t) \nabla\bar{\theta}_{p}(t)\,{\rm d}x
\leq 0.
\end{aligned}
\end{equation}
Using $\phi(t) = \bar{\theta}_{p}(t)^2$
as a test function in \eqref{eq:101} we obtain
\begin{equation}\label{approximate_estimate_03}
\int_{\Omega}
\partial_t^{-\tau} {b}(\bar{u}_{p}(t)) \bar{\theta}_{p}(t)^2
\,{\rm d}x
+\int_{\Omega}
a(\bar{\theta}_{p}(t-\tau)) \nabla \bar{u}_{p}(t)
\cdot\nabla \bar{\theta}_{p}(t)^2
\,{\rm d}x
=0.
\end{equation}
Combining \eqref{approximate_estimate_02} and \eqref{approximate_estimate_03}
we deduce
\begin{equation}\label{approximate_estimate_05}
\begin{aligned}
&\int_{\Omega}\partial_t^{-\tau} \left[
\left( \bar{\theta}_{p}(t) \right)^2
b(\bar{u}_{p}(t))\right]\,{\rm d}x
+\int_{\Omega}\frac{1}{\tau}\left[\bar{\theta}_{p}(t)-\bar{\theta}_{p}(t-\tau)
\right]^2b(\bar{u}_{p}(t-\tau))\,{\rm d}x\\
&+2\int_{\Omega}\lambda(\bar{\theta}_{p}(t-\tau),\bar{u}_{p}(t-\tau))
\nabla \bar{\theta}_{p}(t)\cdot\nabla \bar{\theta}_{p}(t)\,{\rm d}x
\leq0.
\end{aligned}
\end{equation}

Integrating \eqref{approximate_estimate_05} with respect to time $t$ 
we obtain the a priori estimate
(using \eqref{est:uniform_bound_u_02} and \eqref{est:unform_bound_theta_05})
\begin{equation}\label{apriori_est_theta_01}
\sup_{0\leq t \leq T} \int_{\Omega} |\bar{\theta}_{p}(t)|^2 {\rm d}x
+
\int_0^T \|\bar{\theta}_{p}(t)\|^2_{W^{1,2}_{\Gamma_D}(\Omega)} {\rm d}t
\leq {C}.
\end{equation}
From this we have
\begin{equation}\label{est20_theta_a}
\| \bar{\theta}_p \|_{L^2(I;W^{1,2}_{\Gamma_D}(\Omega))} \leq {C}.
\end{equation}


\subsubsection{Further estimates}

To show that $\bar{u}_{p}$ converges to $u$
almost everywhere on $Q_T$ we follow \cite{AltLuckhaus1983}.
Let $k \in \mathbb{N}$ and use
$$
\phi(t) = \partial^{k\tau}_t \bar{u}_{p}(s)
$$
for $j\tau \leq t \leq (j+k)\tau$ with $(j-1)\tau \leq s \leq j\tau$ and 
$1\leq j\leq\frac{T}{\tau}-k$, as a test function in \eqref{eq:101}.
For the parabolic term, we can write
\begin{align*}
&\int_{j \tau}^{(j+k)\tau}\int_{\Omega}
\partial_t^{-\tau} {b}(\bar{u}_{p}(t)) \, \partial^{k\tau}_t \bar{u}_{p}(t)
\,{\rm d}x{\rm d}t \\
&=\frac{1}{k \tau^2} \int_{(j-1)\tau}^{j \tau}
\int_{\Omega}\left({b}(\bar{u}_{p}(t+k\tau)) - {b}(\bar{u}_{p}(t)) \right)
\left(\bar{u}_{p}(t+k\tau) - \bar{u}_{p}(t)\right)
\,{\rm d}x{\rm d}t.
\end{align*}
Hence, summing over $j=1, \dots, p-k$ we obtain the estimate
\begin{equation} \label{est:401}
\begin{aligned}
&\sum_{j=1}^{p-k}\int_{j \tau}^{(j+k)\tau}
\int_{\Omega}\partial_t^{-\tau} {b}(\bar{u}_{p}(t))
\partial^{k\tau}_t \bar{u}_{p}(t) \,{\rm d}x{\rm d}t\\
&\geq \frac{1}{k \tau^2}\int_{0}^{T - k\tau}
\int_{\Omega}\left({b}(\bar{u}_{p}(t+k\tau)) - {b}(\bar{u}_{p}(t)) \right)
\left(\bar{u}_{p}(t+k\tau) - \bar{u}_{p}(t)\right)
\,{\rm d}x{\rm d}t.
\end{aligned}
\end{equation}
Similarly, for the elliptic term, after a little lengthy but straightforward
computation we obtain
\begin{equation} \label{est:402}
\begin{aligned}
&\sum_{j=1}^{p-k}\int_{j \tau}^{(j+k)\tau}
\int_{\Omega}a(\bar{\theta}_{p}(t-\tau)) \nabla \bar{u}_{p}
\cdot\nabla \partial^{k\tau}_t \bar{u}_{p}\,{\rm d}x{\rm d}t\\
&=\sum_{\ell=1}^{k} \sum_{j=1}^{p-k}\int_{(j+\ell-1)\tau}^{(j+\ell)\tau}
\int_{\Omega}\left( a(\bar{\theta}_{p}(t-\tau))
\nabla\bar{u}_{p}\right)\cdot\nabla
\partial^{k\tau}_t \bar{u}_{p}{\rm d}x{\rm d}t\\
&=\sum_{\ell=1}^{k}\int_{\ell\tau}^{T-k\tau+\ell\tau}
\int_{\Omega}a(\bar{\theta}_{p}(t-\tau))\nabla \bar{u}_{p}(t)
\cdot\nabla\partial^{k\tau}_t \bar{u}_{p}(t-\ell\tau)
\,{\rm d}x{\rm d}t \\
&\leq \frac{c_1}{\tau}\int_{Q_T}
|a(\bar{\theta}_{p}(t-\tau))\nabla \bar{u}_{p}|^2
\,{\rm d}x{\rm d}t
+\frac{c_2}{\tau}\int_{Q_T}|\nabla \bar{u}_{p}|^2\,{\rm d}x{\rm d}t\\
&\leq\frac{C}{\tau}.
\end{aligned}
\end{equation}

Combining \eqref{est:401}--\eqref{est:402} and using \eqref{energy_estimate_u} 
we obtain
\begin{equation} \label{est:402b}
\int_0^{T-k\tau}
\left({b}(\bar{u}_{p}(s+k\tau)) - {b}(\bar{u}_{p}(s)) \right)
(\bar{u}_{p}(s+k\tau) - \bar{u}_{p}(s)){\rm d}s
\leq C k \tau.
\end{equation}
Using the compactness argument one can show in the same way
as in \cite[Lemma~1.9]{AltLuckhaus1983} and
\cite[Eqs. (2.10)--(2.12)]{FiloKacur1995}
\begin{equation}\label{eq555}
b(\bar{u}_{p}) \to b(u) \quad \text{in }L^1(Q_T)
\end{equation}
and almost everywhere on $Q_T$.
Since $b$ is strictly monotone,
it follows from \eqref{eq555}
that \cite[Proposition 3.35]{Kacur1990a}
\begin{equation}\label{conv:u00}
\bar{u}_{p} \to u
\quad\text{ almost everywhere on } Q_T.
\end{equation}
Further, in much the same way as
in \eqref{est:402b}, we arrive at
\begin{equation}\label{est21_w}
\int_0^{T-k\tau}|{b}(\bar{u}_{p}(s+k\tau))\bar{w}_{p}(s+k\tau)
- {b}(\bar{u}_{p}(s))\bar{w}_{p}(s) |^2{\rm d}s
\leq C k \tau.
\end{equation}
From this we conclude, using \eqref{est:unform_bound_w_05}, that
\begin{equation}\label{est22_w}
\int_0^{T-k\tau} |\bar{w}_{p}(s+k\tau) - \bar{w}_{p}(s) |^2
{\rm d}s \leq C k \tau.
\end{equation}
Finally, in a similar way, using \eqref{est:unform_bound_theta_05}, we arrive at
\begin{equation}\label{est21_theta}
\int_0^{T-k\tau}|\bar{\theta}_{p}(s+k\tau) - \bar{\theta}_{p}(s)|^2{\rm d}s
\leq {C} k \tau.
\end{equation}

\subsection{Passage to the limit} \label{subsec:limit}

The a priori estimates
\eqref{est:unform_bound_w_05},
\eqref{est:unform_bound_theta_05},
\eqref{energy_estimate_u},
\eqref{est20_w_a},
\eqref{est20_theta_a},
\eqref{est:402b},
\eqref{est22_w},
\eqref{est21_theta}
allow us to conclude that there exist
$u \in L^2(I;W^{1,2}_{\Gamma_D}(\Omega))$,
$w \in L^2(I;W^{1,2}_{\Gamma_D}(\Omega)) \cap L^{\infty}({Q_T})$
and
$\theta \in L^2(I;W^{1,2}_{\Gamma_D}(\Omega)) \cap L^{\infty}({Q_T})$ such that,
letting $p \to +\infty$ (along a selected subsequence),
\begin{gather*}
\bar{u}_{p}  \rightharpoonup u \quad \text{weakly in }
 L^2(I;W^{1,2}_{\Gamma_D}(\Omega)),\\
\bar{u}_{p}  \to u \quad \text{almost everywhere on } Q_T,\\
\bar{w}_{p}  \rightharpoonup w \quad \text{weakly in }
  L^2(I;W^{1,2}_{\Gamma_D}(\Omega)), \\
\bar{w}_{p}  \rightharpoonup w \quad \text{weakly star in } L^{\infty}({Q_T}),\\
\bar{w}_{p}  \to w \quad \text{almost everywhere on } Q_T, \\
\bar{\theta}_{p}  \rightharpoonup \theta \quad \text{weakly in }
 L^2(I;W^{1,2}_{\Gamma_D}(\Omega)), \\
\bar{\theta}_{p}  \rightharpoonup \theta \quad \text{weakly star in }
 L^{\infty}({Q_T}), \\
\bar{\theta}_{p}  \to \theta \quad \text{almost everywhere on } Q_T.
\end{gather*}
The above established convergences
are sufficient for taking the limit $p \to \infty$
in \eqref{eq:101}--\eqref{eq:103}
(along a selected subsequence) to get the weak solution
of the system \eqref{strong:eq1a}--\eqref{strong:eq1l}
in the sense of Definition~\ref{def_weak_solution}.
This completes the proof of the main result stated in 
Theorem~\ref{main_result}.

\subsection*{Acknowledgments}
This research was supported by the project GA\v{C}R~16-20008S
(Michal Bene\v{s}) and by the grant SGS17/001/OHK1/1T/11 provided
by the Grant Agency of the Czech Technical University in Prague 
(Luk\'a\v{s} Krupi\v{c}ka).

\begin{thebibliography}{99}

\bibitem{AdamsFournier1992}
A.~Adams, J.F.~Fournier;
\emph{Sobolev spaces},
Pure and Applied Mathematics 140, Academic Press, 2003.


\bibitem{AltLuckhaus1983}
H.W.~Alt, S.~Luckhaus;
\emph{Quasilinear elliptic-parabolic differential equations},
Mathematische Zeitschrift,
{183} (1983), pp. 311--341.


\bibitem{Aubin1963}
J.-P.~Aubin;
\emph{Un th\'{e}or\`{e}me de compacit\'{e}},
Comptes Rendus de l'Acad�mie des Sciences,
{256} (1963), pp. 5042�-5044.


\bibitem{Bear}
J.~Bear;
\emph{Dynamics of Fluids in Porous Media},
Courier Corporation, 1972.


\bibitem{BenesZeman}% article
M.~Bene\v{s}, J.~Zeman;
\emph{Some properties of strong solutions to nonlinear heat and moisture transport in
multi-layer porous structures},
Nonlinear Analysis: Real World Applications,
13 (2012), pp. 1562--1580.


\bibitem{BenesKrupicka2016}
M.~Bene\v{s}, L.~Krupi\v{c}ka;
\emph{Weak solutions of coupled dual porosity flows in fractured rock mass and structured porous media},
Journal of Mathematical Analysis and Applications,
433 (2016), pp. 543--565.


\bibitem{degond}
P.~Degond, S.~G\'{e}nieys, A.~J\"{u}ngel;
\emph{A system of parabolic equations in nonequilibrium thermodynamics
including thermal and electrical effects},
Journal de Math\'{e}matiques Pures et Appliqu\'{e}es,
76 (1997), pp. 991--1015.


\bibitem{FiloKacur1995}
J.~Filo, J.~Ka\v{c}ur;
\emph{Local existence of general nonlinear parabolic systems},
Nonlinear Analysis,
{24} (1995), pp. 1597--1618.


\bibitem{gallouet}
T.~Gallou\"{e}t, A.~Monier;
\emph{On the regularity of solutions to elliptic equations},
Rendiconti di Matematica,
19 (1999), pp. 471--488.


\bibitem{GiaquintaModica1987}
M.~Giaquinta, G.~Modica;
\emph{Local existence for quasilinear
parabolic systems under nonlinear boundary conditions},
Annali di Matematica Pura ed Applicata,
{149} (1987), pp. 41--59.


\bibitem{groger1989}
K.~Gr\"{o}ger;
\emph{A $W^{1,p}$-estimate for solutions to mixed boundary value
problems for second order elliptic differential equations},
Mathematische Annalen,
283 (1989), pp. 679--687.


\bibitem{jungel2000}
A.~J\"{u}ngel;
\emph{Regularity and uniqueness of solutions to a parabolic system in nonequilibrium thermodynamics},
Nonlinear Analysis,
41 (2000), pp. 669--688.


\bibitem{Kacur1990a}
J.~Ka\v{c}ur;
\emph{On a solution of degenerate elliptic-parabolic systems in Orlicz-Sobolev spaces. I},
Mathematische Zeitschrift,
203 (1990), pp. 153--171.


\bibitem{KufFucJoh1977}
A.~Kufner, O.~John, S.~Fu\v{c}\'{i}k;
\emph{Function Spaces},
Academia, 1977.


\bibitem{LiSun2010}
B.~Li, W.~Sun;
\emph{Global existence of weak solution for nonisothermal multicomponent flow
in porous textile media},
SIAM Journal on Mathematical Analysis,
42 (2010), pp. 3076--3102.

\bibitem{LiSunWang2010}
B.~Li, W.~Sun, Y.~Wang;
\emph{Global existence of weak solution to
the heat and moisture transport system in fibrous porous media},
Journal of Differential Equations,
249 (2010), pp. 2618--2642.


\bibitem{LiSun2012}
B.~Li, W.~Sun;
\emph{Global weak solution for a heat and sweat transport
system in three-dimensional fibrous porous media with condensation/evaporation
and absorption},
SIAM Journal on Mathematical Analysis,
{44} (2012), pp. 1448--1473.


\bibitem{necas1967}
J.~Ne\v{c}as;
\emph{Les methodes directes en theorie des equations elliptiques},
Academia, Prague 1967.


\bibitem{Ozbolt}
J.~O\v{z}bolt, G.~Balabani\'{c}, G.~Peri\v{s}ki\'{c}, M.~Ku\v{s}ter;
\emph{Modelling the effect of damage on transport processes in concrete},
Construction and Building Materials,
24 (2010), pp. 1638--48.


\bibitem{PinderGray}
G.F.~Pinder, W.G.~Gray;
\emph{Essentials of Multiphase Flow in Porous Media},
Wiley-Interscience, New Jersey, 2008.


\bibitem{Vala2002}
J.~Vala;
\emph{On a system of equations of evolution with a
non-symmetrical parabolic part occuring in the analysis of moisture
and heat transfer in porous media},
Applications of Mathematics,
{47} (2002), pp. 187--214.


\bibitem{Roubicek2005}
T.~Roub\'{i}\v{c}ek,
Nonlinear Partial Differential Equations with Applications,
Birkh\"{a}user, 2005.


\bibitem{Weidemaier1991}
P.~Weidemaier;
\emph{Local existence for parabolic problems with
fully nonlinear boundary condition; An $L_p$-approach},
Annali di Matematica Pura ed Applicata,
{160} (1991), pp. 207--222.


\end{thebibliography}


\end{document}