%Final submitted to ARMA on September 26, 2002
%Submitted to ARMA on Feb. 12, 2002
%%February 12, 2002
%%Last update on November 20, 2003
\documentclass[10pt,twoside]{amsart}
\usepackage{amsmath, amsfonts, amsthm, amssymb, psfig, 
varioref,
enumerate,latexsym}
%\usepackage[mathcal]{euscript}
\newcommand{\tab}{\,\,\,\,}
\newcommand{\re}{\mathbb{R}}
\newcommand{\ze}{\mathbb{Z}}
\newcommand{\IM}{\mathbb{M}}
\newcommand{\ink}{\rule{0.40\baselineskip}{.40\baselineskip}}
\newcommand{\rhup}{\rightharpoonup} 
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea} {\begin{eqnarray}}
\newcommand{\nbea} {\begin{eqnarray*}}
\newcommand{\eea} {\end{eqnarray}}
\newcommand{\neea} {\end{eqnarray*}}
\newcommand{\intom} {\int_{\Omega}}
\newcommand{\M} {{\mathcal M}}
\newcommand{\linf} {\liminf_{n \to \infty}}  
\def\dist{\operatorname{dist}}
\newtheorem{definition}{Definition}[section]
   \newtheorem{lemma}[definition]{Lemma}
\newtheorem{theorem}[definition]{Theorem}
\newtheorem{proposition}[definition]{Proposition}
\newtheorem{corollary}[definition]{Corollary}
\newtheorem{remark}[definition]{Remark}
\newtheorem{example}[definition]{Example}
\newtheorem{claim}{Claim}
\def\BbbZ{{\Bbb{z}}}
\def\email#1{E-mail: {\tt #1}}
\def\address#1{{\def\\{\cr}\halign{##\hfill\cr #1}}}
\def\setA{V}
\def\setB{{V'}}
\def\setC{{V''}}
\def\beg{\begin{equation}}
\oddsidemargin 	.5in
   \evensidemargin .5in
   \textwidth	5.85in
   \topmargin     -0.45in
   \textheight	9.00in
   \numberwithin{equation}{section}
\title[Navier-Stokes Equations for Compressible, Reacting Flow]
{\Large{Global Solutions to A Model For Exothermically Reacting, 
compressible Flows with large discontinuous initial data }}

\subjclass{35B40, 35D05, 76N10, 35B45}
\keywords{Global discontinuous solutions, large 
discontinuous initial data,
evolution of large jump discontinuities, asymptotic 
behavior, combustion, 
Navier-Stokes equations, difference approximations, energy 
estimates, 
total variation estimates, uniform bounds.}
\begin{document}
\maketitle
\author \centerline{Gui-Qiang \ Chen \qquad David \ Hoff
   \qquad Konstantina Trivisa}
\begin{abstract}
We prove the global existence of solutions of the 
Navier-Stokes
equations describing the dynamic combustion of a 
compressible,
exothermically reacting fluid, and we study the large-time 
behavior of
solutions, giving necessary and sufficient conditions for 
complete
combustion in certain cases.  The adiabatic constants and 
specific
heats of the burned (product) and unburned (reactant) 
fluids may
differ, and the initial data may be large and discontinuous.
\end{abstract}

\vskip20pt

\section {\bf Introduction.}

We prove the global existence of solutions of the 
Navier-Stokes
equations describing the dynamic combustion of a 
compressible,
exothermically reacting fluid, and we study the large-time 
behavior of
solutions, giving necessary and sufficient conditions for 
complete
combustion in certain cases.  The adiabatic constants and 
specific
heats of the burned (product) and unburned (reactant) 
fluids may
differ, and the initial data may be large and discontinuous.

In the present context, the Navier-Stokes equations express 
the
conservation of mass, the balance of momentum and energy, 
and
two-species chemical kinetics as follows:
%% 1.1
\begin{equation}\label{E:1.1.1}
\begin{cases}
v_t-u_x=0,\\
u_t+p(v,\theta,Z)_x=\displaystyle\Big(\frac {\varepsilon 
u_x}{v}\Big)_x,\\
(e+\frac{1}{2}u^2+qZ)_t+(up(v,\theta,Z))_x=
\displaystyle\Big(\frac{\lambda \theta_x
+\varepsilon uu_x}{v}\Big)_x,\\
Z_t+K\phi(\theta)Z=0.
\end{cases}
\end{equation}
\noindent
Here $v=1/\rho,u,\theta,$ and $Z$ represent specific 
volume, velocity,
temperature, and reactant mass fraction, respectively, and 
are the
basic unknown functions of $x$ and $t$; $e$ and $p$ are the 
specific
internal energy and pressure, and $\phi$ is the 
reaction-rate
function (the precise assumptions we make on the state 
functions
$p(v,\theta,Z), e(\theta,Z)$ and $\phi(\theta)$ will be given below);  
$\varepsilon,\lambda,K$, and $q$ are positive constants
describing viscosity, thermal conduction, reaction rate, 
and the difference in heats between  reactant and product.  
The above equations
are written in Lagrangian coordinates, so that $x=$ constant
corresponds to the trajectory of a fixed fluid particle.

First, in Theorem 1.1 we establish the global-in-time existence 
of solutions
with quite general initial data, obtaining a number of 
results
concerning the regularity of solutions, as well as bounds 
which in
certain cases are independent of time.  In Theorem 1.2 we 
give a result
concerning the propagation of singularities in solutions 
for which $v$
and $Z$ are initially piecewise $H^1$.  These results show 
that, while
there is a sort of ``parabolic regularization'' of $u$ and 
$\theta$,
singularities in $v,p,Z,e,u_x,$ and $\theta_x$ persist for 
all time,
convecting along particle trajectories, and evolving 
according to 
certain dynamics derivable from the corresponding 
Rankine-Hugoniot
conditions.  In Theorem 1.3 we discuss the large-time 
behavior of
solutions under the assumption that the regularity 
estimates of Theorem 1.1 are independent of time, 
showing that the solution converges as
$t\rightarrow \infty$ to an equilibrium determined by the 
initial mass and energy.  
Finally, in Theorem 1.4 we give both necessary conditions
and sufficient conditions that complete combustion occurs, 
that is,
that $Z\equiv 0$ in the equilibrium state, in terms of the 
total
initial energy.
\bigskip

We begin by describing some of the physical background and 
deriving the
specific forms of the state functions $p(v,\theta, Z)$ and
$e(\theta,Z)$.
In this discussion, $x$ will be the usual (Eulerian) 
spatial coordinate
in physical space, and $\rho,u,\theta,$ and $Z$ will be 
unknown functions
of $x$ and $t$.  The density $\rho$ should satisfy the 
standard continuity
equation
%1.2
\begin{equation}\label{E:1.1.2}
\rho_t+(\rho u)_x=0,
\end{equation}
%
and the reactant mass fraction of a fixed fluid particle 
will decrease
according to standard chemical kinetics
%1.3
\begin{equation}\label{E:1.1.3}
Z_t + uZ_x=-K\phi(\theta)Z.
\end{equation}
%
Here $K$ is a rate constant, and the reaction function 
$\phi$ is
assumed to be an increasing, Lipschitz function on 
$[0,\infty)$
satisfying
%1.4
\begin{equation}\label{E:1.1.4}
C^{-1}\phi (\theta) \le \theta^2\phi^\prime(\theta) \le 
C\phi
(\theta)
\end{equation}
%
for a constant $C$, and
%1.5
\begin{equation}
\begin{cases}\label{E:1.1.5}
\phi (\theta)=0, \quad  0\le \theta \le\theta_i,\\
\phi (\theta) >0, \quad  \theta > \theta_i,
\end{cases}
\end{equation}
%
where $\theta_i>0$ is the ignition temperature.  These 
conditions are
satisfied, for example, by the Arrhenius function $\phi
(\theta)=e^{-A/\theta},\theta \gg \theta_i.$
%
Thus, when $Z=1$, for a given fluid particle, that particle 
consists
entirely of reactant.  When its temperature rises above the 
ignition
temperature, combustion will occur, resulting in the 
conversion of some
or all of its mass to the product species, and decreasing 
$Z$
accordingly.  When $Z=0$, complete combustion has occurred, 
and the
particle consists entirely of the product species.

We may regard the reactant and product as separate fluids 
with
respective densities $\rho_1=Z\rho$ and $\rho_2=(1-Z)\rho,$ 
having the same
velocities and temperatures at each point, and 
characterized by
polytropic equations of state
%1.6
\begin{equation}\label{E:1.1.6}
\begin{cases}
P_j=c_j(\gamma_j-1)\rho \theta,\\
e_j=c_j\theta, \ \qquad \quad\quad j=1,2.
\end{cases}
\end{equation}
%
Here $P_j$ and $e_j$ are pressures and specific internal
energies, $\gamma_j>1$ are adiabatic constants, and $c_j$ 
are
(constant) specific heats.  The two fluids satisfy the 
momentum and
energy equations
%1.7
\begin{equation}\label{E:1.1.7}
(\rho_ju)_t + (\rho_ju^2+P_j)_x=\varepsilon_ju_{xx},
\end{equation}
%1.8
\begin{equation}\label{E:1.1.8}
(\rho_jE_j)_t +(\rho_jE_ju + P_j u)_x=\lambda_j\theta_{xx} +
\varepsilon_j(uu_x)_x+(-1)^j q_j \rho \phi (\theta)Z,
\end{equation}
%
where $E_j=e_j+u^2 /2,$ and $\varepsilon_j, \lambda_j,$ and 
$q_j$ are
positive constants.  The last term on the right side of 
(\ref{E:1.1.8})
represents the rate of energy lost to the reactant or 
gained by the
product as a result of the chemical reaction.

Now define macroscopic pressure and energy by
%1.9
\begin{equation}\label{E:1.1.9}
\begin{cases}
P=P_1 +P_2  \,\,\, \hbox{(Dalton's law)},\\
\rho E=\rho_1E_1+\rho_2E_2,\\
e=E-u^2 /2,
\end{cases}
\end{equation}
%
and system parameters
%no #
\begin{equation*}
\varepsilon=\varepsilon_1 +\varepsilon_2, \quad 
\lambda=\lambda_1 +\lambda_2, \quad q=q_2-q_1.
\end{equation*}
%
Summing in (\ref{E:1.1.7}) and (\ref{E:1.1.8}), we then 
obtain that 
%1.10
\begin{equation}\label{E:1.1.10}
(\rho u)_t + (\rho u^2 +P)_x=\varepsilon u_{xx},
\end{equation}
%
and
%1.11
\begin{equation}\label{E:1.1.11}
(\rho E)_t + (\rho Eu +Pu)_x=\lambda \theta_{xx} +
q\rho\phi(\theta)Z.
\end{equation}
%
\noindent We define a variable specific heat $c(Z)$ by
%1.12
\begin{equation}\label{E:1.1.12}
c(Z)=c_1Z +c_2(1-Z),
\end{equation}
%
so that, by (\ref{E:1.1.6}) and (\ref{E:1.1.9}),
%1.13
\begin{equation}\label{E:1.1.13}
e(\theta,Z)=Z e_1 + (1-Z)e_2=c(Z)\theta.
\end{equation}
%
We also define $\gamma(Z)$ by
%1.14
\begin{equation}\label{E:1.1.14}
\gamma(Z)= \frac{Zc_1 \gamma_1 + (1-Z)c_2 \gamma_2}{Zc_1 + 
(1-Z)c_2},
\end{equation}
%
so that, again by (\ref{E:1.1.6}) and (\ref{E:1.1.9}),
%1.15
\begin{equation}\label{E:1.1.15}
P(\rho,\theta,Z)= a(Z)\rho\theta,
\end{equation}
%
where $a(Z)=c(Z)(\gamma(Z)-1)$.
Thus the functions of state $P$ and $e$ for the composite
fluid are exactly those of a polytropic fluid, as in 
(\ref{E:1.1.6}),
but with adiabatic ``constant'' $\gamma$ and specific
heat $c$ depending on $Z$.  The equations (\ref{E:1.1.2}),
(\ref{E:1.1.3}), (\ref{E:1.1.10}), and (\ref{E:1.1.11}), 
together with
the state equations (\ref{E:1.1.13}) and (\ref{E:1.1.15}), 
then comprise
a closed system for the unknown functions $\rho,u, \theta,$ 
and $Z$.
Rewriting in terms of the Lagrangian coordinate, still 
denoted by $x$,
and taking
%1.16
\begin{equation}\label{E:1.1.16}
p(v,\theta, Z)=P(v^{-1},\theta,Z)=a(Z)\theta/v,
\end{equation}
%
we then obtain the system (\ref{E:1.1.1}).

In the present paper we consider a fluid of finite total 
mass confined
to a bounded region, so that, without loss of generality, 
$0<x<1$.  We impose boundary conditions
%
%1.17 (1.2.1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{equation}\label{E:1.2.1}
u(x,t)=0, \quad  \theta_x(x,t)=0, \quad  x=0,1.
\end{equation}
%
The solutions we obtain will be weak solutions in the 
usual
sense: writing the system (\ref{E:1.1.1}) in the form
%NO #
\begin{equation*}
w_t^j + f^j(w)_x=(b^{jk}(w)w^{k}_x)_x
+g^j(w),
\end{equation*}
%
we say that $w$ is a weak solution provided that $w$ is 
suitably
integrable, and
%1.18 (1.2.2)
\begin{equation}\label{E:1.2.2}
\int^{1}_{0} 
w^j(x,\cdot)\psi(x,\cdot)dx|^{t_2}_{t_1}=\int^{t_2}_{t_1}
\int^1_0(w^j\psi_t + (f^j-b^{jk}w^k_x)\psi_x + g^j \psi)dxdt
\end{equation}
%
for all $t_2 \ge t_1 \ge 0$ and all $\psi \in 
C^1([0,1]\times [0,\infty))$,
with the exception that, when $w^j=u,$ (\ref{E:1.2.2}) is 
required to hold
only for $\psi$ satisfying the boundary conditions $\psi
(0,t)=\psi(1,t)=0$.

We consider initial data
%1.19 (1.2.3)
\begin{equation}\label{E:1.2.3}
(v,u,\theta,Z)|_{t=0}=(v_0,u_0,Z_0,\theta_0)
\end{equation}
%
satisfying
%1.20 (1.2.4)
\begin{equation}\label{E:1.2.4}
\begin{cases}
C_0^{-1} \le v_0(x) \le C_0, \ \theta_0(x)\ge C_0^{-1}, \ 
0\le Z_0(x)\le 1,\\
||\theta_0 ||_{L^2} + || u_0 ||_{L^4} \le C_0
\end{cases}
\end{equation}
%
for a constant $C_0 >0$, which may be arbitrarily large.  
Finally, we
define the following functionals for a given weak solution
$(v,u,\theta,Z)$ of (\ref{E:1.1.1}), (\ref{E:1.2.1}),
and (\ref{E:1.2.3}):
%No #
%\begin{equation*}
\begin{align*}
\mathcal{D}(t)=&\sup_{0\le s \le t}\Big(||u(\cdot,s)^2||^2 
                 +||\theta(\cdot,s)||^2\Big)   
% + ||u(\cdot,s) \theta (\cdot,s)||^2
\quad  \qquad
 \qquad \ \ \notag\\
 & + \int^1_0 \int^t_0 \Big(\theta^{-1} u_x^2 + u^2 u_x^2 + 
\theta_x^2
 + \phi(\theta)Z^2\Big)dxdt,\notag 
\end{align*}
%\end{equation*}
%1.21 (1.2.5)
%\begin{equation}
\begin{align}\label{E:1.2.5} 
\mathcal{E}(t)=&\sup_{0\le s \le t}\left(\sigma(s)||u_x
(\cdot,s)||^2 + \sigma^2(s)||\theta_x(\cdot,s)||^2\right)\\ 
 %
& + \int^1_0\left(||(u_x,\theta_x)(\cdot,s)||^2 +
 \sigma(s)||u_t(\cdot,s)||^2
  + \sigma^2(s)||\theta_t(\cdot,s)||^2\right)ds,\notag\\
%
\mathcal{F}(t)=&\sup_{0\le s \le t}
\left(\sigma^2(s)||u_t(\cdot,s)||^2
      +\sigma^3(s)||\theta_t(\cdot,s)||^2\right)\notag \\
%
& + \int^t_0\left(\sigma^2(s)||u_{xt}(\cdot,s)||^2 
  + \sigma^3(s)||\theta_{xt}(\cdot,s)||^2\right)ds,\notag
\end{align}
%\end{equation}
%
where $\sigma(t)=\min\{1,t\}$, and $||\cdot||$ denotes the 
usual norm
in $L^2([0,1])$.

The following theorem then gives our main results 
concerning the
existence and regularity of solutions of the system (\ref{E:1.1.1}).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bigskip
\begin{theorem}\label{T:1.1}
Assume that the functions $p(v,\theta,Z),e(\theta,Z)$, and
$\phi(\theta)$ are as described above in 
{\rm (\ref{E:1.1.4})}, {\rm (\ref{E:1.1.5})}, and 
{\rm (\ref{E:1.1.12})}--{\rm (\ref{E:1.1.16})}, 
and let initial data $(v_0,u_0,\theta_0,Z_0)$ be given, 
as in {\rm (\ref{E:1.2.4})}.
Then the initial-boundary value problem 
{\rm (\ref{E:1.1.1})}, {\rm (\ref{E:1.2.1})}, 
and {\rm (\ref{E:1.2.3})} has a global weak solution 
$(v,u,\theta,Z)$ for which
%
%1.22 (1.2.55)
\begin{equation}\label{E:1.2.55}
\begin{cases}
v\in C([0,\infty);L^p([0,1])),& 1\le p < \infty,\\
u \in C ([0,\infty); L^p([0,1])),& 1\le p<4,\\
\theta \in C ((0,\infty);L^2([0,1])), & \theta(\cdot,t)
\rightharpoonup \theta_0  \ \text{weakly in} \ L^2 \  
\text{as}  \ t\to 0,\\
Z\in C([0,\infty);L^p([0,1])),& 1\le p \le \infty.
\end{cases}
\end{equation}
%
 Also, given $T>0$, there is a constant $M$ depending on 
$T,C_0,$ and
on the system parameters such that, for $t \in(0,T]$,
%1.23 (1.2.6)
\begin{equation}\label{E:1.2.6}
\begin{cases}
M^{-1} \le v(x,t)\le M, \quad a.e. \, x \in [0,1],\\
M^{-1} \le \theta (x,t) \le M \max \{1,t^{-1}\}, \quad  x \in [0,1],\\
0\le Z(x,t)\le 1, \quad a.e. \ x \in [0,1],
\end{cases}
\end{equation}
%
and
%1.24 (1.2.7)
\begin{equation}\label{E:1.2.7}
\mathcal{D}(T) + \mathcal{E}(T) + \mathcal{F}(T) \le M.
\end{equation}
%
In particular, $u$ and $\theta$ are H$ \ddot{o}$lder 
continuous on compact
sets in $[0,1]\times(0,\infty),$ and
%No #
\begin{equation*}
u(\cdot,t), \ \theta(\cdot,t), \ \big(\frac{\varepsilon 
u_x}{v}-p\big) (\cdot,t)
\in H^1([0,1])
\end{equation*}
%
for $t>0$.  Finally, in the case that $c_1=c_2$ and 
$\gamma_1=\gamma_2$ {\rm(}see {\rm (\ref{E:1.1.12})} and 
{\rm (\ref{E:1.1.14}))}, the
bounds in {\rm (\ref{E:1.2.6})} and {\rm (\ref{E:1.2.7})} 
hold with a constant $M$ which is independent of $T$.  
\end{theorem}
%%%%%%%%%%%%%%%%%%%%%%%% end therorem 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We prove Theorem \ref{T:1.1} in Sections 2 and 3, by first 
deriving
certain a priori bounds for suitable approximate solutions, 
then
applying these bounds to obtain the solutions of the 
theorem as strong
limits of the approximate solutions.  Both aspects of the 
analysis involve interesting features.  
The most important of the a priori estimates are the pointwise 
bounds for the temperature and the specific
volume, and are achieved by a modification of the analysis 
in Chen, Hoff, and Trivisa \cite{CHT} for nonreacting flows
and Chen \cite{C} for reacting flows, which was 
motivated by
an idea of Kazhikov and Shelukhin \cite{KS}.  
These
pointwise bounds depend crucially on the existence of an 
entropy, which
is a convex function of the dependent variables satisfying 
the
Clausius-Duhem inequality, the latter guaranteeing that the 
integral of the entropy over space is monotone in time.  
We show in the appendix
that such an entropy can indeed be constructed, but only 
for $v,\theta$
in a compact set in $(0,\infty) \times (0,\infty)$;  we 
also show that there is in fact no single entropy with the 
required properties for all
$(v,\theta)\in (0,\infty) \times (0,\infty)$.  In other 
words, the existence
of a suitable entropy, which is required for the derivation 
of pointwise
bounds for $v$ and $\theta$, depends on these very 
pointwise bounds.
This circularity leads to apparently unavoidable Gronwall 
inequalities,
resulting in the time-dependence of the constant $M$ in 
Theorem~
\ref{T:1.1}.  These difficulties do not occur if $c_1=c_2$ 
and $\gamma_1=\gamma_2,$ however, and in this case 
{\it time-independent} pointwise bounds can be obtained for
$v$ and $\theta$.
%\cite{CHT2}.  

Once these pointwise bounds have been 
established,
time independent or not, the higher-order regularity 
estimates in (\ref{E:1.2.5}) and (\ref{E:1.2.7}) can be 
derived by
straightforward but rather long and technical energy 
estimates.  All
these a priori bounds are given in Section 2.

The other noteworthy aspect of the proof of 
Theorem~\ref{T:1.1}
is the proof of strong compactness of approximate solutions 
in the
absence of uniform regularity for the variables $Z$ and 
$v$.  Now,
since $Z$ can be represented in terms of $\theta$ by 
integrating its
evolution equation along particle trajectories, compactness 
for
approximate $Z$'s can be obtained from the corresponding 
compactness of
approximate $\theta$'s.  For the specific volume $v$, the 
argument is
much more subtle, and strong compactness is shown to result
from certain structural properties of approximate solutions 
derived
from considerations relating to the cancellation, 
persistence, and
evolution of singularities in solutions.  These compactness 
arguments
are given in Section 3.

The basic facts concerning the propagation of singularities 
in
solutions of (\ref{E:1.1.1}) are best understood at the 
level of
piecewise smooth solutions.  These facts are given in 
Theorem
\ref{T:1.2} below and are similar to corresponding results 
for
nonreacting flows (see \cite{Ho2} example).  Still, there 
are
interesting new features which arise, due to the 
variability of the
mass fraction.  To describe these results, we first suppose 
that
$(v,u,\theta,Z)$ is a piecewise smooth solution, having 
one-sided limits
along a curve $(y(t),t)$, and we derive at the heuristic 
level the
various conditions that guarantee that $(v,u, \theta, Z)$ 
is a weak
solution across $(y(t),t)$. First, in light of the 
regularity
estimates in Theorem \ref{T:1.1} for $u$ and $\theta$,
%
%1.25 (1.3.0)
\begin{equation}\label{E:1.3.0}
[u(y,t)] = [\theta(y,t)]=0,
\end{equation}
%
where $[u(y,t)] = u(y+,t)-u(y-,t)$, etc. Next, applying the 
standard
Rankine-Hugoniot conditions to the mass equation in 
(\ref{E:1.1.1}), we
find that $\dot{y}(t)[v(y,t)]=[u(y,t)]=0$.  Combining this 
with the observation
that singularities in $Z$ must remain at fixed points of 
$x-$space, we
see that  $y(t)=y$ is constant.  We then conclude that 
singularities remain stationary in Lagrangian space, 
and therefore convect along
particle trajectories in physical space.  Applying the
Rankine-Hugoniot condition to the other equations in 
(\ref{E:1.1.1}),
we then find that, at $(y,t)$,
%
%1.26 
\begin{equation}\label{E:1.3.1}
\displaystyle \Big[\frac{\varepsilon u_x}{v}-p \Big]=0, \quad 
\Big[\frac{\theta_x}{v}\Big]=0, \quad  [Z]_t + K \phi(\theta) [Z]=0.
\end{equation}

\noindent We solve these equations formally as follows.  
First, 
%No #
%1.27 (1.3.2)\
\begin{equation} \label{E:1.3.2}
\frac{d}{dt}[\log v (y,t)] = \Big[ \frac{v_t}{v} \Big ]= 
\Big [\frac{u_x}{v}\Big ] = \varepsilon^{-1}[p]
= \alpha [\log v] + \beta[Z],
\end{equation}
%
where
%1.28 (1.3.3)
\begin{equation}\label{E:1.3.3}
\alpha (y,t) = \varepsilon^{-1}\theta (y,t) \ \displaystyle
\frac{c_+(\gamma_+ -1) + c_- (\gamma _- -1)}{2} \ \
\displaystyle \frac {[v^{-1}(y,t)]}{[\log v(y,t)]},
\end{equation}
%
and
%1.29 (1.3.4)
\begin{equation}\label{E:1.3.4}
\beta(y,t) 
=\varepsilon^{-1}
\displaystyle 
 \frac{v_+ + v_-} {2 v_+v_-} \big (c_1 (\gamma_1 -1)+c_2(\gamma_2-1)\big)
\theta (y,t) 
\end{equation}
%
(here $c_+ = c(Z(y+,t)$), etc.).  Solving the ode's in 
(\ref{E:1.3.1}) and
(\ref{E:1.3.2}), we then obtain
%
%1.30 (1.3.5)
%\begin{equation}
\begin{align}\label{E:1.3.5}
[\log v(y,t)]= & \  e^{-\int^t_0\alpha (y,s)ds} [\log 
v_0(y)]\notag \\
& \qquad + \int^t_0 e^{-\int^t_s\alpha(y,\tau)d\tau} 
\beta(y,s)[Z(y,s)]ds,
\end{align}
%\end{equation}
%
and
%1.31 (1.3.6)
\begin{equation}\label{E:1.3.6}
[Z(y,t)]= \ e^{-\int^t_0 K\phi(\theta(y,s))ds} [Z_0(y)].
\end{equation}
%
Explicit expressions for $[u_x]$ and $[\theta_x]$ can then 
be obtained
from \eqref{E:1.3.1}.  Observe that $\alpha$ is strictly 
negative, so
that the strengths of jump discontinuities in $v$ and $Z$, 
hence in
$u_x$ and $\theta_x$, decrease in time. Unlike the case of 
nonreacting
flows, however, these strengths may not decrease to zero as 
$t\to
\infty$, because $\theta(y,t)$ may be asymptotically less 
than
$\theta_i$.  We shall discuss this point further, following 
the
statement of Theorem 1.4.
Also see [3,7] for the asymptotic behavior of
solutions for inviscid flows.


These observations are made precise in the following 
theorem.
%%%%%%%%%%%%%%%%%%%%% theorem 1.2 %%%%%%%%%%%%%%%
\begin{theorem}\label{T:1.2}
Assume that the hypotheses of Theorem {\rm \ref{T:1.1}} are 
in force, and suppose in
addition that $v_0$ and $Z_0$ are piecewise in $H^1,$ 
having simple
jump discontinuities at points $y_1 < y_2 <\cdots <y_n$ in 
$(0,1)$ (thus $v_0,Z_0\in H^1((y_{j-1},y_j))$ for each $j$).  
Then, for the solution $(v,u,\theta,Z)$ of Theorem~{\rm \ref{T:1.1}}, 
the quantities $v(\cdot,t), p(\cdot,t), Z(\cdot, t), u_x (\cdot,t),$ 
and $\theta_x(\cdot,t)$ are piecewise in
$H^1$ for $t>0$ and therefore have one-sided limits at each 
$y_j$ for
$t>0$, and the jump conditions {\rm (\ref{E:1.3.0})} and
{\rm(\ref{E:1.3.1})} and the representations 
{\rm(\ref{E:1.3.5})} and
{\rm(\ref{E:1.3.6})} hold in a strict pointwise sense at 
each $y=y_j$ for
$t>0$.
\end{theorem}
%
Next we turn to the study of the large-time behavior of 
solutions.
This will require the assumption that the estimates of 
Theorem \ref{T:1.1} hold
independently of time, and as indicated in 
Theorem~\ref{T:1.1}, this
time-independence is known so far only for the special case 
that
$c_1=c_2$ and $\gamma_1=\gamma_2$.  We emphasize, however, 
that the
large-time results given below in Theorems \ref{T:1.3} and 
\ref{T:1.4}
do not depend directly upon the latter assumptions in any 
way, but
instead are consequences of the weak forms of the equations
(\ref{E:1.1.1}) and time-independent estimates
(\ref{E:1.2.6}) and (\ref{E:1.2.7}).  
In the following theorem, we apply these
estimates to show that the solution tends to a 
time-asymptotic
steady-state solution, and we characterize this 
steady-state in terms
of the asymptotic mass-fraction and the initial mass and 
energy, given
by
%1.32 (1.4.1)
\begin{equation}\label{E:1.4.1}
\bar{v}_0=\int^1_0 v_0(x)dx,
\end{equation}
%%1.33 (1.4.2)
\begin{equation}\label{E:1.4.2}
E_0=\int^1_0 \big(c(Z_0(x)) \theta_0(x) + {\textstyle 
\frac{1}{2}} u_0(x)^2 + qZ_0(x)\big)dx.
\end{equation}
%
%%Theorem 1.3%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bigskip
\begin{theorem}\label{T:1.3}
Let $(v,u,\theta,z)$ be a weak solution of  \ {\rm 
\eqref{E:1.1.1}} as in
Theorem {\rm \ref{T:1.1}}, satisfying
{\rm (\ref{E:1.2.6})} and {\rm (\ref{E:1.2.7})} with $M$ 
independent of $T$,
and let $\bar{v}_0$ and $E_0$ be as above in {\rm 
(\ref{E:1.4.1})} and {\rm (\ref{E:1.4.2})}.  
Then there is a function $Z_\infty$ such that, as $t\to\infty$,
%
%no #
\begin{equation*}
Z(x,t)\to Z_\infty(x), \quad  a.a. \ x;
\end{equation*}
%
and, if $v_\infty(x)$ and $\theta_\infty$ are defined by
% 
%1.34 (1.4.3)
\begin{equation}
\theta_\infty = \displaystyle \frac{E_0 - q \displaystyle 
\int_0^1
Z_\infty(x)dx}{\displaystyle \int_0^1 c(Z_\infty(x))dx},
\end{equation}
%
and
%%1.35 (1.4.4)
\begin{equation}
v_\infty(x) = \displaystyle \frac{\gamma(Z_\infty(x))-1}
{\displaystyle \int_0^1 (\gamma (Z_\infty(x))-1)dx} \ \bar{v}_0,
\end{equation}
%
then, as $t\to \infty$,
%1.36 (1.4.5)
\begin{equation}
||(u(\cdot,t), \theta(\cdot,t)-\theta_\infty)||_{H^1([0,1])}\to 0,
\end{equation}
%
%1.37 (1.4.6)
\begin{equation}
||v(\cdot,t)-v_\infty||_{L^p ([0,1])}\to 0, \quad 1\le p<\infty.
\end{equation}
%
\smallskip
	
\end{theorem}
Of particular interest are the conditions that guarantee that 
the steady-state 
$(v_\infty,0,\theta_\infty, Z_\infty)$ is constant in 
space and is characterized by complete combustion, 
$Z_\infty \equiv 0$.  In view
of the fourth equation in (\ref{E:1.1.1}), this should be 
true when
the asymptotic temperature $\theta_\infty$ is greater than 
the ignition 
temperature $\theta_i,$ and the latter should hold when the 
total
energy $E_0$ in the system is sufficiently large.  This is 
indeed the
case; the precise statements are given in the following 
theorem.
\bigskip
%%%%%%%%%%%% theorem 1.4 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{theorem}\label{T:1.4}
Let the hypotheses and notations of Theorem {\rm \ref{T:1.3}} be in
force and assume that the initial mass-fraction $Z_0$ is not
identically zero. 

{\rm (a)} \quad
If $Z_\infty \equiv 0,$  then $E_0 \ge c_2 \theta_i.$ 

{\rm (b)}\quad Conversely, if
\begin{equation} \label{E:1.4.7}
E_0> \max \lbrace c_2\theta_i,\Big(c_2+(c_1-c_2)\int^1_0 Z_0dx\Big)
\theta_i + q\int_0^1 Z_0 dx\rbrace,
\end{equation}
then $Z_\infty \equiv 0$ and $v_\infty \equiv \bar{v}_0.$
\end{theorem}
\smallskip

The gap between the necessary condition in (a) and 
the sufficient condition in (b) can be closed 
for the special class of
solutions for which $v=v_0$ is constant, $u=0,$ and 
$\theta,Z,$ and $e$
depend only on $t$.  In this case, the fluid has neither 
kinetic energy nor elastic potential energy, and all 
effects are
thermal:
% in article Remark 1.6
\bea
\frac{de}{dt} &\!\! = \!\! & K q \phi({\theta}) Z,\nonumber 
\\
\frac{dZ}{dt} &\!\! = \!\! & - K \phi({\theta}) Z, 
\nonumber\\ 
{\theta}\ \ &\!\! = \!\!  &  \frac{e}{(c_1-c_2)Z + c_2}. 
\nonumber 
\eea
The solution $(e,Z)$ thus remains on the line segment 
$$ e+ qZ = E_0 = \,\mbox{constant},$$
so that 
$$
e(t) = E_0 - q Z(t) \,\quad  \mbox{and} \,\quad 
\theta(e(t), Z(t)) = \frac{E_0 - qZ(t)}{(c_1 - c_2) Z(t) + c_2}. 
$$
The trajectory $\{(e(t), Z(t)): 0\le t < \infty \}$ is 
therefore disjoint
from the set of points $\{(e,Z): \theta(e,Z) \le {\theta}_i 
\}$ 
if and only if 
$$ 
E_0 - qZ > (c_2 + (c_1 - c_2) Z){\theta}_i 
$$
for $Z \in (Z_{\infty}, Z_0].$ It then follows easily that 
$Z_{\infty} =0$ if and only if 
$$
E_0 \ge c_2 {\theta}_i \,\quad \mbox{and} \,\quad 
E_0 > (c_2 + (c_1-c_2)Z_0) {\theta}_i + q Z_0.
$$
%
\smallskip

Of course, the results of the present paper apply to 
nonreacting flows
as well: if $Z_0\equiv 0$, then $Z(\cdot,t)\equiv 0$ for all
$t,c(Z)\equiv c_2,$ and $\gamma(Z)\equiv \gamma_2$.  The 
bounds in
(\ref{E:1.2.6}) and (\ref{E:1.2.7}) are then independent of 
time, and the
time-asymptotic state given in Theorem \ref{T:1.3} is simply
$(\bar{v}_0,0,E_0/c_2,0)$.  In particular, singularities 
disappear in the
time-asymptotic limit (see (\ref{E:1.3.5})).  This same 
``asymptotic
compactness'' is known even when an external force is 
added, at least
for the case of barotropic flows, and is the basis for the 
existence of
a global attractor (Hoff and Ziane \cite{HZ}).  For the 
more general
case of reacting flows considered here, the situation is 
quite a bit
more complicated:  when $\theta_\infty <\theta_i$ the 
steady-state will
not be constant in space, and, more importantly, 
singularities in
solutions will not disappear in the time-asymptotic limit 
(see
(\ref{E:1.3.5}) and (\ref{E:1.3.6})).

The existence and asymptotic behavior of global smooth 
solutions for
reacting, compressible Navier-Stokes flows, with 
constant adiabatic
exponent and additional species diffusion, was first 
established by Chen \cite{C}.  
The case of {\em nonreacting} flows in one space
dimension has been treated by many authors.  
See, for example, Chen, Hoff, and Trivisa \cite{CHT} on the 
well-posedness and asymptotic behavior of global solutions with 
large discontinuous initial data,
Matsumura and Yanagi \cite{MY} for the isothermal case, 
Amosov and Zlotnick \cite{AZ} for the existence of weak solutions 
for certain initial-boundary data, and  Hoff \cite{Ho1}--\cite{Ho4} 
for the global well-posedness and large-time behavior for small 
discontinuous initial data, including the case of several 
space variables. 

We also remark that the results of the present paper carry
over to the initial boundary value problem with Dirichlet 
boundary data for the temperature $\theta$ for (\ref{E:1.1.1}), 
as well as to the corresponding system written in Eulerian 
coordinates, (\ref{E:1.1.2}), (\ref{E:1.1.10}), and (\ref{E:1.1.11}).
Some related discussions for dynamic combustion of reacting fluids
can be found in [2,7,8,9,11,12,21,22,24].

\bigskip



%%%%%%%%%%%%%%%%%  section 2  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section {\bf A Priori Estimates for Approximate Solutions} 
 
%
In this section we begin the proof of the existence result,
Theorem~1.2, by constructing approximate solutions based 
upon a
suitable semidiscrete difference scheme.  We shall show in 
a sequence
of lemmas that these approximate solutions satisfy all the 
estimates of
Theorems~1.2, independently of the spatial discretization.  
These
estimates will then be applied in Section 3 to obtain the 
solution of
Theorems~1.2 as the strong limit of these approximate 
solutions, and
the more general solutions of Theorem~1.1 will be obtained 
as the
limits of smooth solutions corresponding to mollified 
initial data.

\medskip
Let $h$ be an increment in $x$ such that $Qh = 1$ for some
$Q=Q(h) \in {\ze}_+$, $x_k = kh$ for $k \in \{0,1,\cdots, 
Q-1\}$, and 
$x_j = jh$ for $j\in\{\frac{1}{2}, \frac{3}{2},\cdots, 
Q-\frac{1}{2}\}.$
Approximations $(v_j,u_k, {\theta}_j,Z_j)(t)$ to 
$(v(x_j,t),u(x_k,t), {\theta}(x_j,t),Z(x_j,t))$ are then 
constructed 
as follows:
%2.1
\bea
& &\dot{v}_j=\delta u_j,  \label{E:2.1}\\
%2.2
& &\dot{u}_k+\delta p_k
   =\varepsilon\delta\left(\frac{\delta u}{v}\right)_{\! 
k}, \label{E:2.2}\\
%2.3
& & c_v(Z_j) \dot{\theta}_j + p_j \delta u_j  = 
 \varepsilon\frac{(\delta u_j)^2}{v_j} + \lambda \delta
\left(\frac{\delta \theta}{v}\right)_{\!j} + 
K \phi({\theta}_j) Z_j(q 
+c_v'(Z_j){\theta}_j),\label{E:2.3}\\ 
%2.4
& & \dot{Z}_j + K \phi({\theta}_j) Z_j  =  0. \label{E:2.4}
\eea 
Here $p_j=p(v_j,{\theta}_j, Z_j)$,  $v_k$ is taken
to be the average $v_k = \frac{v_{k+\frac{1}{2}} + 
v_{k-\frac{1}{2}}}{2}$
with  $j \in \{\frac{1}{2}, \frac{3}{2}, \cdots, 
Q-\frac{1}{2}\}$ 
and $k\in\{0,1,\cdots,Q-1\}$, and $\delta$ is the operator 
defined by
$$
\delta w_l = \frac{w_{l+ \frac{1}{2}} - w_{l - 
\frac{1}{2}}}{h},
\qquad l=k,\,\, \hbox{or}\,\, j.
$$ 
For the time being we assume only that initial data 
$(v_j,u_k, {\theta}_j,Z_j)(0)$ for 
\eqref{E:2.1}--\eqref{E:2.4} have been 
specified and satisfy
%2.5 
\bea 
u_0 = u_Q =  0, \qquad
\delta {\theta}_0 = \delta {\theta}_Q = 0, 
\label{E:2.5} 
\eea
and
%2.6
\bea
\begin{cases}
C_0^{-1} \le v_j(0) \le C_0,  \tab {\theta}_j(0)\ge 
C_0^{-1}, \quad 0\le
Z_j(0) \le 1, \\
\sum_k u_k^4(0)h +\sum_j {\theta}_j^2(0) h \le C_0. 
\end{cases}
\label{E:2.6}
\eea 
We also assume that there are distinguished points 
$0<x_{k_1}<x_{k_2}<\cdots < x_{k_N}<1, N=N(h), N^4h\le 1$,
such that
%2.7
\bea 
\sum_{k=k_i}\left(|[v_{k_i}(0)]|+|[Z_{k_i}(0)]|\right)+
\sum_{k\ne k_i}\left(|\delta v_k(0)|^2+|\delta 
Z_k(0)|^2\right) h  \le C_0. 
\label{E:2.7}
\eea

Standard theory of ordinary differential equations applies 
to show 
that the initial value problem \eqref{E:2.1}--\eqref{E:2.7} 
has a unique 
solution $(v_j,u_k, {\theta}_j,Z_j)(t)$, defined at least 
for small time. 
The a-priori bounds to be derived in this section will show 
that these 
solutions exist {\it globally} in time, and will provide 
sufficient 
compactness both to extract limiting solutions as $h\to 0$ 
and
to determine their asymptotic behavior.

Here and in what follows, $C>0$ will denote a generic 
constant
independent of $h$ and $t$, $M=M(m_0)$ is a constant 
independent of $h$
but depending on $T$, for a given interval $[0,T]$, only 
through a 
parameter $m_0=m_0(T)$, which in turn is also independent 
of $h$ and $N$.  

Using the difference equations \eqref{E:2.1}--\eqref{E:2.7} 
and the boundary 
conditions, we derive some basic but important energy 
estimates described
in Lemmas \ref{L2.1}--\ref{L2.2}. 
\medskip
%%%%%%%%%%%%%%%%%  lemma 2.1
\begin{lemma}\label{L2.1}
%2.8
\bea 
& &  \sum_j v_j(t) h = 1, \label{E:2.8}\\
%2.9
& & 0 \le Z_j(t) \le 1,  \label{E:2.9}\\
%2.10
& & \sum_j Z_j(t)h+\int_0^t\sum_j 
K\phi({\theta_j}(s))Z_j(s)h ds 
      = \sum_j Z_j(0)h,  \label{E:2.10} \\ 
%2.11
& & \sum_j\left(c_v(Z_j){\theta}_j(t)+q Z_j(t)\right)h 
+\frac{1}{2}\sum_k u_k^2(t)h 
          \label{E:2.11} \\
& &\qquad \,\,\,\,  = 
\sum_j\left(c_v(Z_j(0)){\theta}_j(0)+qZ_j(0)\right)h 
   +\frac{1}{2} \sum_k u_k^2(0)h, 
           \nonumber\\
%2.12
&&  \sum_j \frac{1}{2} Z_j^2(t) h + \int_0^t \sum_j
 K \phi(Z_j(s)) Z_j^2(s)hds = \sum_j \frac{1}{2} Z_j^2(0)h. 
\label{E:2.12} 
\eea
\end{lemma}
\noindent
{\bf Proof:}
{}From equation \eqref{E:2.1}, we obtain
$\frac{d}{dt}\sum_jv_jh=\sum_j\delta u_jh=u_Q-u_0= 0,$
which yields \eqref{E:2.8}.
Regarding $Z_j$, the set $[0,1]$ is
invariant for equation \eqref{E:2.4}, for each $j$, since 
$\phi(\theta_j)\ge 0$;
and then \eqref{E:2.9} follows.
Summing \eqref{E:2.4} and integrating the resulting 
equation on
the interval $[0,t]$, we conclude \eqref{E:2.10}.
Summing \eqref{E:2.3} and integrating the resulting 
equation on $[0,t]$,
we obtain
%2.13
\bea 
&&\sum_j\left(c_v(Z_j(t)){\theta}_j(t) +q Z_j(t)\right)h
  +\int_0^t\sum_j (p_j \delta u_j) h ds 
\label{E:2.13}\nonumber\\
&&\,\,\,\, = \sum_j \left(c_v(Z_j(0)) {\theta}_j(0)+ qZ_j(0)\right)h 
   + \int_0^t \sum_j \frac{\varepsilon (\delta 
u_j)^2}{v_j}h ds.
\eea    
Now, summing  \eqref{E:2.2} and integrating the resulting 
equation on
$[0,t]$,  we  have
%2.14 
\bea
\frac{1}{2} \sum_k u_k^2(t) h - \int_0^t \sum_j (p_j \delta 
u_j)h ds =
 \frac{1}{2} \sum_k u_k^2(0) h    - \varepsilon \int_0^t 
\sum_j
\frac{(\delta u_j)^2}{v_j} h ds. \label{E:2.14}
\eea 
Adding \eqref{E:2.13} and \eqref{E:2.14} together yields 
\eqref{E:2.11}.
Finally, using \eqref{E:2.4}, we have
\bea
\frac{d}{dt} \left( \sum_j \frac{1}{2} Z_j^2(t) h \right) = 
\sum_jZ_j(t)\dot{Z}_j(t)h=-K\sum_j\phi({\theta}_j(t))Z_j^2(t)h. 
\nonumber
\eea 
Integrating this last equation on the interval $[0,t]$, 
we obtain \eqref{E:2.12}.
\hfill $\blacksquare$
%%%%%%%%%%%%%%%%%%%%%%%%%Lemma 2.2%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{lemma}\label{L2.2}
There exists $m_0=m_0(T)$, independent of $h$ and $N$,
such that, for each $T>0$ and any $t\in [0,T]$,
%2.15
\bea
E(t)+\int_0^t(V(s)+W(s))ds \le m_0< \infty,\label{E:2.15}
\eea
with
%No #
\bea
\begin{cases}
E(t) =\sum_j \left(c_v(Z_j) ({\theta}_j -1 - 
\log{{\theta}_j}) 
 + a(Z_j)(v_j-1-\log{v_j})\right) h + \frac{1}{2} \sum_k 
u_k^2 h, \\
V(t) = \sum_j \left( \frac{\varepsilon (\delta u_j)^2}{v_j 
e_j} \right)h 
+ \lambda \sum_k \left(\frac{(\delta {\theta}_k)^2}{v_k
{\theta}_{k+\frac{1}{2}}{\theta}_{k-\frac{1}{2}}}\right)h,\\ 
W(t) =-K\sum_j\left(\frac{{\theta}_j-1}{{\theta}_j}\right)
  \phi({\theta}_j)Z_j (q+c_v'(Z_j)\theta_j).
\end{cases} \nonumber
\eea
In the case that $\gamma$ and $c_v$ are independent of $Z$, 
that is,  $\gamma_1 = \gamma_2$ and $c_1 = c_2$, 
the parameter $m_0$ is independent of $T$.
% that yields the time independent bounds.  
\end{lemma}
\noindent
{\bf Proof:}
The choice of the energy $E=E(t)$ under consideration is
inspired by the analysis in the appendix which leads us to 
a 
precise form of a physical entropy for \eqref{E:1.1.1}.
Differentiating  $E=E(t)$ with respect to $t$ and taking 
into
consideration the relation \eqref{E:1.1.2}  and the 
momentum 
equation \eqref{E:2.3}, we obtain
%no #
\bea
\frac{d}{dt} E(t)=\!\! && \!\!\!\!\!\!
\sum_j \left( c_v{'}(Z_j)({\theta}_j-1-\log{{\theta}_j}) 
         + a'(Z_j)(v_j-1-\log{v_j})\right)
          \dot{Z}_j h \nonumber \\
& &\!\! - 
\sum_j \frac{1}{{\theta}_j}\left( \frac{\varepsilon (\delta 
u_j)^2}{v_j} 
  +\lambda\delta \left( \frac{\delta \theta}{v}\right)_j
  + K\phi({\theta}_j) Z_j(q+c_v'(Z_j)\theta_j) \right)h 
\nonumber\\
& & \!\! +
   \sum_j\left(c_v(Z_j) \dot{\theta}_j
         + a(Z_j) \dot{v}_j \right)h 
    + \sum_k u_k \dot{u}_k h.\nonumber
\eea 
Integrating the last relation on the interval $[0,t]$, we 
get
$$ 
E(t)-E(0)=\sum_{l=1}^4 A_l(t),
$$
where
%no #
\bea
& & \!\!\!\!\!\!\!\!\!\!\!\!  A_1(t)
= \sum_j \int_0^t a(Z_j)\dot{v}_j ds\, h,
\nonumber\\
& & \!\!\!\!\!\!\!\!\!\!\!\! 
A_2(t)=-\sum_j \int_0^t \frac{1}{{\theta}_j} 
\left(\frac{\varepsilon (\delta u_j)^2}{v_j} 
   +\lambda\delta\left(\frac{\delta \theta}{v}\right)_j 
\right)ds\, h, 
   \nonumber \\
& & \!\!\!\!\!\!\!\!\!\!\!\! 
A_3(t) =\sum_j \int_0^t \left(c_v(Z_j) \dot{\theta}_j 
-\frac{1}{{\theta}_j}K\phi({\theta}_j)Z_j(q+c_v'(Z_j)\theta_j)\right)ds\, h 
 + \int_0^t \sum_k u_k \dot{u}_k h ds, \nonumber\\ 
& & \!\!\!\!\!\!\!\!\!\!\!\!
A_4(t)=\sum_j\int_0^t\left(c_v{'}(Z_j)
   ({\theta}_j-1-\log{{\theta}_j}) 
   + a'(Z_j)(v_j -1-\log{v_j})\right)\dot{Z_j}ds\, h. 
\nonumber
\eea
Now, $A_2(\cdot)$ is equivalent to
\bea
A_2(t)
=-\int_0^t\left(\sum_j\frac{\varepsilon(\delta u_j)^2}{v_j 
{\theta}_j} h +\lambda\sum_k
 \frac{(\delta{\theta}_k)^2}{v_k{\theta}_{k+\frac{1}{2}}
   {\theta}_{k-\frac{1}{2}}} h \right)ds. \nonumber 
\eea
Using Lemma \ref{L2.1}, we have
\bea
\!\!\!\!\!\! \!\!\!\!\!\!A_3(t) && \!\!\!\!\!\! 
\!\!\!\!\!\!  
=\sum_j \left( 
c_v(Z_j(t))\theta_j(t)-c_v(Z_j(0))\theta_j(0)\right)h
 +\left(\frac{1}{2}\sum_ku_k^2(t)h-\frac{1}{2}\sum_k 
u_k^2(0) h \right)
   \nonumber\\
\!\!\!\!\!\! \!\!\!\!\!\!& & \!\!\!\!\!\!  \!\!\!\!\!\! 
 \quad +\sum_j\int_0^t \left(-c_v{'}(Z_j)\dot{Z}_j 
{\theta}_j 
    -\frac{1}{\theta_j}K\phi(\theta_j)Z_j
     \left(q+c_v'(Z_j)\theta_j\right) \right) ds h 
\nonumber\\
\!\!\!\!\!\! \!\!\!\!\!\!& & \!\!\!\!\!\!  \!\!\!\!\!\! 
=-\sum_kqZ_j(t)h +\sum_k qZ_j(0) h 
 +\sum_j\int_0^t 
K\phi(\theta_j)Z_j\left(c_v{'}(Z_j){\theta}_j 
     -\frac{1}{\theta_j}\left(q+c_v'(Z_j)\theta_j\right)    
        \right) ds h \nonumber\\
\!\!\!\!\!\! \!\!\!\!\!\!& & \!\!\!\!\!\!  \!\!\!\!\!\! 
=\sum_j\int_0^t 
KZ_j\phi(\theta_j)\frac{\theta_j-1}{\theta_j}
       \left(q+c_v'(Z_j)\theta_j\right)ds\, h. \nonumber
\eea
Therefore, by Lemma \ref{L2.1}, we obtain
%2.16
\bea
&& \label{E:2.16}\\
& &
\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!
\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\
E(t) + \int_0^t \left(V(s) + W(s)\right)ds \nonumber \\
& &  \!\!\!\!\!\!\!\!\!
=E(0) +\sum_j\int_0^t a(Z_j) \dot{v}_jds\,h \nonumber\\
&&\,\,\,\,   
-\int_0^t\sum_j\left(c_v{'}(Z_j)({\theta}_j-1-\log{{\theta}_j}) 
   + a'(Z_j)(v_j-1-\log{v_j})\right)K\phi(\theta_j)Z_j h 
ds\nonumber\\
& &  \!\!\!\!\!\!\!\!\!
\le 
E(0)+C_1\int_0^t\sum_j\left(c_v(Z_j)({\theta}_j-1-\log{{\theta}_j}) 
   + a(Z_j)(v_j-1-\log{v_j})\right)h ds \nonumber\\
&&\,\,\,\,\,  +\sum_ja(Z_j(t))v_j(t)h
   -\sum_ja(Z_j(0))v_j(0)h
   -\int_0^t\sum_j K \phi(\theta_j) a'(Z_j) Z_j v_j h\,ds. 
 \nonumber
\eea
Hence,
\bea
& & E(t) + \int_0^t \left(V(s) + W(s)\right)ds \nonumber \\
& & \tab \tab \le  E(0) 
+C_1\int_0^t\sum_j\left(c_v(Z_j)({\theta}_j-1-\log{{\theta}_j}) 
              + a(Z_j)(v_j-1-\log{v_j})\right)h ds +C_2T.
\nonumber
\eea
The result \eqref{E:2.15} now follows by using the Gronwall 
inequality.

Observe that, if one assumes that $\gamma$ and $c_v$ are 
independent of $Z$,
that is, ${\gamma}_1 = {\gamma}_2$ and $c_1 = c_2$, 
then \eqref{E:2.16} and Lemma \ref{L2.1} imply that $m_0$ 
is independent of $T$.
%which yields the time independent bounds. 
The detailed analysis 
of the qualitative behavior of the solutions in this 
case including  the stability results are presented in 
\cite{CHT2}. Some further discussion on the model and its
variants is presented in the appendix.

\hfill $\blacksquare$

\medskip
Next we obtain  pointwise bounds for the specific 
volume $v=v_j(t)$  with the aid of Lemmas 
\ref{L2.1}--\ref{L2.2}.
This estimate will be crucial  in the forthcoming analysis.
%%%%%Lemma 2.3%%%%%%%%%%%
\begin{lemma}\label{L2.3}
There exists $M=M(m_0)>0$, independent of $h$, such that, 
for each $T >0$, 
%2.17 
\bea
M^{-1}  \le v_j(t) \le M <\infty, \qquad 
\hbox{for all }\,  t \in [0,T].
\label{E:2.17}%3.8
\eea
\end{lemma}

\noindent
{\bf Proof:}
The proof is based on the  derivation of two different 
representations
for the specific volume $v$. The pointwise bounds are 
obtained by combining 
the analysis given in Chen, Hoff, and Trivisa \cite{CHT} and Chen \cite{C}
with  an idea in Kazhikov and Shelukhin \cite{KS} 
and using equation \eqref{E:2.2}, Lemmas 
\ref{L2.1}--\ref{L2.2}, and 
boundary condition  \eqref{E:2.5}
with special consideration of the features of the system.
In particular, the pointwise bounds are essential to guarantee 
the continuous dependence of the solutions on the initial data (cf. [16]).
Here we present only the main steps of the proof. 

\noindent
{\em Step 1.}
Lemma \ref{L2.2} yields that, for any $t > 0$, there exist 
$j_1(t)$ and 
positive constants $\alpha, \beta$ with $\alpha < \beta$ 
such that 
\bea
\alpha \le v_{j_1}, \quad \!\!\! {\theta}_{j_1} \le \beta. 
\nonumber
\eea
Following now  similar line of arguments to the one 
presented in \cite{CHT},
we get
%2.18  
\bea
v_j(t)=\frac{1+ \frac{1}{\varepsilon}\int_0^t 
a(Z_j(s))P(s)Q_j(s){\theta}_j(s)ds}
   {P(t) Q_j(t)}, \label{E:2.18}%3.9
\eea
where
%2.19
\bea
\begin{cases}
P(t) =  v_{j_1(t)} (0) {\mbox{exp}} \{\frac{1}{\varepsilon}
\int_0^t p_{j_1(t)}(s) ds \}, \\  
 Q_j(t) =  \frac{1}{v_j(0) v_{j_1}(t)} 
\mbox{exp}\{\frac{1}{\varepsilon}\sum_{(j_1,j)}(u_k(0)-u_k(t))h\}.
\end{cases}
\label{E:2.19}%3.10
\eea
Furthermore, using the momentum equation, we obtain that 
there 
exists $j_2 = j_2(t)$ such that 
%2.20
\bea
\sum_{k<j_2(t)}u_k(0)h+\int_0^t{\sigma}_{j_2(t)}(s)ds
=& & \!\!\!\!\! \sum_j\left(v_j(0)\sum_{k < 
j}u_k(0)h\right)h \nonumber \\
 & & \!\!\!\!\! -  \int_0^t \sum_j 
\left(\frac{u_{j+\frac{1}{2}}^2 +
u_{j-\frac{1}{2}}^2}{2} + a(Z_j) {\theta}_j \right)hds + 
O(h^{1/2}),
\label{E:2.20}%3.17
\eea 
where
${\sigma}_j=\left(\frac{\varepsilon\delta u}{v}\right)_j-p_j$,
and $O(h^{1/2})$ denotes the terms bounded by $Mh^{1/2}$.

{\em Step 2.}
Next, we derive (cf. \cite{CHT}) another representation for
$v_j=v_j(t)$. For  any $t \ge 0$,
%2.21 
\bea
\label{E:2.21}&&\\%v
&&v_j(t)=(1+O(h^{1/2}))D_j(t) \mbox{exp} 
\left\{- \frac{1}{\varepsilon} \int_0^t 
\sum_l\left(\frac{u_{l+\frac{1}{2}}^2+u_{l-\frac{1}{2}}^2}{2} 
  +a(Z_l) {\theta}_l\right)h ds \!\right\} \nonumber\\
&& \times 
\left\{1 + \frac{1 + O(h^{1/2})}{\varepsilon} \int_0^t
\frac{a(Z_j) 
{\theta}_j}{D_j(s)}\mbox{exp}\left\{\frac{1}{\varepsilon}
\int_0^s\sum_l\left(\frac{u_{l+\frac{1}{2}}^2+u_{l-\frac{1}{2}}^2}{2} 
   + a(Z_l){\theta}_l \right)h d{\tau}\right\} ds \!
\right\},\nonumber 
\eea 
where
%no #
\bea
D_j(t) = v_j(0) \mbox{exp}\left\{\frac{1}{\varepsilon} 
\left(
\sum_{(j_2(t), j)}u_k(t)h - \sum_{k<j}u_k(0) h 
+\sum_lv_l(0)\left(\sum_{k<l}u_k(0)h\right)h\right)\right\}, 
\nonumber
\eea
and $j_2(t)$ is determined in \eqref{E:2.20}.

This result is obtained by following similar arguments to 
the one given
 in  Lemma \ref{L2.2} in \cite{CHT}.

\medskip
\noindent
{\em Step 3.} Next we obtain pointwise bounds for 
$v_j(\cdot)$.
We first set
\bea
m_v(t) = \min_{x \in[0,1]} v(x,t), \tab\ 
m_{\theta}(t) = \min_{x \in[0,1]} {\theta}(x,t),\nonumber\\
M_v(t) = \max_{x \in[0,1]} v(x,t), \tab\   
M_{\theta}(t) = \max_{x \in[0,1]} {\theta}(x,t). \nonumber
\eea
Using Lemma \ref{L2.2}, we conclude that there exist 
$M=M(m_0) >0, 
\hat{\beta} = \hat{\beta}(m_0)$, and 
$\tilde{\beta}=\tilde{\beta}(m_0)$ 
such that 
$$ 
M^{-1} \le D_j(t), Q_j(t) \le M, 
\qquad \hat{\beta} \le \sum_j {\theta}_j(t)h 
\le\tilde{\beta}.
$$

The second representation \eqref{E:2.21} for $v_j(\cdot)$ 
implies
%2.22
\bea
M_v(t) \le M\mbox{exp}\{- a_1 t\} 
\left(1+\int_0^tM_{\theta}(s)
\mbox{exp}\{a_1 s\} ds \right),
\label{E:2.22}%3.24
\eea
where 
$a_1=\frac{1}{\varepsilon}\hat{\beta}\min_{0\le Z\le 
1}a(Z).$

Next, we obtain bounds for $M_\theta(\cdot)$ and
$m_{\theta}(\cdot),$ namely,
%2.23
\bea
M_{\theta}(t)
& & \!\!\!\!\!\!\!\!\!  \le
 M (1 + M_v(t) B(t)), \label{E:2.23}\\%3.25
%2.24
m_{\theta}(t) & & \!\!\!\!\!\!\!\!\!  \ge
 M^{-1} \left(1-M_v(t) B(t)\right),  \label{E:2.24}%3.26a
\eea
where
$$
B(t)= \sum_k \frac{(\delta {\theta}_k)^2}{ 
{\theta}_{k+\frac{1}{2}}
{\theta}_{k-\frac{1}{2}} {v_k}} h.
$$
Substituting \eqref{E:2.23} into \eqref{E:2.22}, we get
\bea 
M_v(t)\le M\exp\{-a_1 t\} 
\left\{1+\int_0^t\exp\{a_1 s\}
\left(1+M_v(s)B(s)\right) ds \right\}. \nonumber
\eea
Setting  now
$A(t) = M_v(t) \exp \{a_1 t\}$ and taking into 
consideration that 
$\int_0^t B(\tau) d \tau \le M, $
we get
%2.25 
\bea
M_v(t) = \exp \{-a_1 t\} A(t) \le M. \label{E:2.25}
\eea 
Next, by using \eqref{E:2.24}, \eqref{E:2.25}, and
the second representation \eqref{E:2.21} for $v_j(\cdot)$ 
yields
that there exists $t_0>0$ such that
%2.26
\bea
m_v(t)\ge \frac{M^{-1}}{2}, \qquad t\ge t_0>0.
\label{E:2.26}%3.28
\eea

On the other hand,
\bea
P(t) 
& & \!\!\!\!\!\! = \sum_j P(t) v_j(t) h 
=\sum_j\left(\frac{1+
  \frac{1}{\varepsilon}\int_0^t a(Z_j(s)){\theta}_j(s) P(s) 
Q_j(s) ds}{Q_j(t)}
\right) h \nonumber \\
& & \!\!\!\!\!\! \le M\left(1+\int_0^t P(s)(\sum_j 
{\theta}_j(s)h)ds\right) 
  \le M \left(1 + \int_0^t P(s)ds\right),
\nonumber 
\eea
which gives
$$ 
P(t) \le M \exp \{Mt\}.
$$
Therefore, from \eqref{E:2.18},
$$ 
v_j(t) \ge \frac{M^{-1}}{P(t)} \ge M^{-1} \exp \{-Mt\}
\ge M^{-1}\exp\{-Mt_0\}, \quad \hbox{when}\quad 0\le t\le 
t_0.
$$
Combining this last relation with \eqref{E:2.26},
we obtain the result. 
\hfill $\blacksquare$
\medskip

The next lemma establishes a lower bound for the 
temperature $\theta$. 
It will be clear in the sequel that, away from the initial 
line $\{t=0\}$,
the temperature $\theta$ has indeed uniform bounds.
\smallskip

\begin{lemma}\label{L2.4} 
There exists $M=M(m_0)>0$, with $m_0$ independent of $h$ 
and $N$, 
such that, for all $t$,
%2.27
\bea
{\theta}_j (t) \ge \frac{1}{M(t+1)}.  \label{2.27}%3.30
\eea
\end{lemma}

\noindent
{\bf Proof:}
Set $\omega_j = \frac{1}{{\theta}_j}$.
Multiplying (2.3) by $\{-{\theta}_j^{-2} \}$, we obtain
%no number
\bea
c_v(Z_j)\dot{\omega}_j = a(Z_j)\frac{\omega_j \delta u_j}{v_j}
-\frac{\varepsilon\, {\omega}_j^2 (\delta u_j)^2}{v_j} 
+ \lambda \delta \left(\frac{\delta \omega}{v}\right)_j 
- K q\phi({\theta}_j)Z_j \omega_j \left(\omega_j + 
\frac{c_v{'}(Z_j)}{q} \right), 
   \nonumber
\eea
which implies
%2.28  
\bea
c_v(Z_j)\dot{\omega}_j \le \frac{M}{v_j} 
+ \lambda \delta \left(\frac{\delta \omega}{v}\right)_j 
-K q \phi({\theta}_j)Z_j \omega_j \left(\omega_j + 
\frac{c_v{'}(Z_j)}{q}\right). 
 \label{E:2.28}%
\eea 
Multiplying the last inequality by 
$\{2p{\omega}_j^{2p-1}\}, p>1$, we obtain
\bea
\frac{d}{dt} (c_v(Z_j)\omega_j^{2p}) 
\le 2pM \frac{\omega_j^{2p-1}}{v_j}
 + 2p\lambda\delta \left(\omega^{2p-1} \frac{\delta \omega}{v}\right)_j
  -(2p+1)K \frac{\phi({\theta}_j)}{\theta_j} Z_j 
c_v'(Z_j)w\omega_j^{2p-1}. 
 \nonumber
\eea
Summing over all $j$ and using Lemma \ref{L2.3}, we obtain 
\bea
\sum_j \frac{d}{dt} (c_v(Z_j)\omega_j^{2p}) h 
\le M p \sum_j \omega_j^{2p-1} h.  \nonumber
%+ p M_2 \sum_j 
%\frac{\phi({\theta}_j)}{{\theta}_j}c_v(Z_j)w^{2p-1} 
%\left(w_j + \frac{2p-1}{2p} \frac{c_v{'}(Z_j)}{q} \right) h. 
\eea
Now,
\bea
\left(\sum_j {\omega_j}^{2p-1}h \right) 
\le\left(\sum_j \omega_j^{2p} h\right)^{\frac{2p-1}{2p}}
 \left(\sum_j h\right)^{\frac{1}{2p}}
\le\left(\sum_j \omega_j^{2p} h\right)^{\frac{2p-1}{2p}}. 
\nonumber
\eea
%By Lemma \ref{L3},
%\bea
%\left(\sum_j {w_j}^{2p-1}h \right) \le 
%\left(\sum_j w_j^{2p} h\right)^{\frac{2p-1}{2p}}
% \left(\sum_j (pM)^{2p}h\right)^{\frac{1}{2p}}. \nonumber
%\eea
Set
$$
g(t)=\left(\sum_j c_v(Z_j)\omega_j^{2p} h \right)^{\frac{1}{2p}}.
$$
Then 
%2.29
\bea
\frac{d}{dt} \left(g(t)^{2p}\right) \le Mp 
g(t)^{2p-1}.\label{E:2.29}
% \,\,\,\, \mbox{or} \,\,\,\, \frac{d}{dt} g(t) \le M. 
\label{te}
\eea
Therefore,
$$
\left(\sum_j{\omega}_j(t)^{2p}h\right)^{\!\! \frac{1}{2p}}
\le C_1^{\frac{1}{2p}}
  \left(\sum_j{\omega}_j(0)^{2p}h\right)^{\!\! 
\frac{1}{2p}}+Mt,
$$ 
where $C_1=\max_{0\le Z\le 1}c_v(Z)/\min_{0\le Z\le 
1}c_v(Z)$.
Letting now $p\to\infty$, we get
\bea
{\| \omega(t) \|}_{L^{\infty}} \le 
\|\omega(0)\|_{L^{\infty}}+Mt
\le M(1+t).  \nonumber
\eea
This completes the proof.
\hfill $\blacksquare$


\medskip
The next lemma will be useful in establishing estimates on 
the variation of $Z$ (Lemma \ref{L2.7}).
Note, Lemmas \ref{L2.1}--\ref{L2.4}, in combination with the 
bounds
%2.30
\bea
&& \label{E:2.30}\\%3.31
&& c_v(Z_j){\theta}_i(t)
\le\left(\sum_j c_v(Z_j){\theta}_j(t)h\right)h^{-1} \le 
Mh^{-1}, \,\, 
u_i^2(t)\le\left(\sum_j u_j^2(t)h\right)h^{-1}\le Mh^{-1}, 
\nonumber
\eea
show that the system of ordinary differential equations
\eqref{E:2.1}--\eqref{E:2.4} is solvable for all $t>0$ for 
fixed $h>0$.

%%%%%%%%%%%%%%%%% lemma 2.5
\begin{lemma}\label{L2.5} 
\noindent
For ${\hat u}_j$ determined by
$u_{j+ \frac{1}{2}}^3-u_{j-\frac{1}{2}}^3 
  =3{\hat u}_j^2(u_{j+\frac{1}{2}}-u_{j-\frac{1}{2}}),$
there exists $M=M(m_0)$ such that
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%**
%2.31
\bea 
\sum_j (c_v(Z_j)({\theta}_j - 1)^2  + ({\theta}_j -1)u^2_{j 
+ \frac{1}{2}}
& \!\!\!\!  + \!\!\!\! & L u_{j+ \frac{1}{2}}^4 + Z_j^2) h
 +   \int_0^t \sum_k (\delta {\theta}_k)^2hds
\nonumber \\
& \!\!\!\! + \!\!\!\! & 
\int_0^t\sum_j\left(K{\phi}_j Z_j^2+\hat{u}_j^2{\theta}_j^2 
 +\hat{u}_j^2(\delta u_j)^2\right)hds \le M, 
\label{E:2.31}%3.32
\eea
where $L > (4 \min_{0\le Z\le 1} c_v(Z))^{-1}$.
\end{lemma}

\noindent
{\bf Proof:}
The result \eqref{E:2.31} is obtained by combining 
the {\em energy bounds},  which are derived by using  the 
momentum equation 
\eqref{E:2.2},  with the results of Lemma \ref{L2.1}.

Multiplying the momentum difference equation \eqref{E:2.2} 
by $\{u_k^3\}$, 
summing over all $k$, and integrating on $[0,t]$, we obtain
%2.32
\bea
\sum_k u_k^4(t) h + \int_0^t \sum_j {\hat u_j}^2 (\delta 
u_j)^2 h ds
\le M + M \int_0^t \sum_j {\hat u}^2 e_j^2 h ds. 
\label{E:2.32}%*1 
\eea
Next, multiplying \eqref{E:2.2} by $\{u_k\}$, we get
%2.33
\bea
{\left(\frac{1}{2} u_k^2 \right)}_{\!\!t} + u_k \delta p_k 
= - \varepsilon
\frac{(\delta u_{k-\frac{1}{2}})^2}{v_{k- \frac{1}{2}}} + 
\varepsilon
\delta \left( {\tau}_{-}u \frac{\delta u}{v}\right)_{\!\!k},
\label{E:2.33}%*2
\eea 
where $({\tau}_{-} u)_j = u_{j - \frac{1}{2}}.$
Notice that \eqref{E:2.3} implies that
%2.34
\bea
c_v(Z_j)({\theta}_j - 1)_t + p_j \delta u_j 
= \frac{\varepsilon (\delta u_j)^2}{v_j} +
\lambda \delta {\left(\frac{\delta 
\theta}{v}\right)}_{\!\!j} 
+ K\phi(\theta_j) Z_j (q+c_v'(Z_j)\theta_j). 
\label{E:2.34}%*3
\eea
Multiplying \eqref{E:2.34} by $\{2({\theta}_j -1)\}$, we get
%2.35 
\bea
&& c_v(Z_j)(({\theta}_j -1)^2)_t +2 \lambda \frac{(\delta 
{\theta}_j)^2}{v_j}
-2\lambda\delta\left(({\theta}_j -1) \frac{\delta 
{\theta}}{v}\right)_j
 \label{E:2.35} \\%*4
&&= -2 ({\theta}_j -1) p_j \delta u_j 
+2 \varepsilon \frac{({\theta}_j -1) (\delta u_j)^2}{v_j}
+2K\phi({\theta}_j)({\theta}_j-1)Z_j(q+c_v{'}(Z_j){\theta}_j). 
\nonumber
\eea
Therefore,
\bea
& &\!\!\!\!\!\!\!\!\!\! \left(c_v(Z_j)({\theta}_j 
-1)^2\right)_t +
2 \lambda \frac{(\delta {\theta}_j)^2}{v_j} -
2 \lambda \delta \left(({\theta}_j -1) \frac{\delta 
{\theta}}{v} \right)_j
 \nonumber \\
& & \,\,\,\,\,\!\!\!\!\!
=-2 ({\theta}_j -1) \delta(u {\tau}_+ p)_j 
 + 2 ({\theta}_j -1) u_{j+\frac{1}{2}} \delta 
p_{j+\frac{1}{2}} 
 + 2\varepsilon \frac{({\theta}_j -1) (\delta u_j)^2}{v_j} 
\nonumber\\
&& \quad 
+2K\phi({\theta}_j)({\theta}_j-1)Z_j(q+c_v{'}(Z_j){\theta}_j)
     -K\phi(\theta_j)(\theta_j-1)^2c_v'(Z_j)Z_j.
   \nonumber
\eea
Using \eqref{E:2.33}, we now arrive at
\bea
& &\!\!\!\!\! (c_v(Z_j)({\theta}_j -1)^2)_t
 +2 \lambda \frac{(\delta {\theta}_j)^2}{v_j} -
2 \lambda \delta \left(({\theta}_j -1) \frac{\delta 
{\theta}}{v} \right)_j
\nonumber \\
& &\,\,\,\,\, 
= -2 ({\theta}_j -1) \delta(u {\tau}_+ p)_j 
 + 2 ({\theta}_j -1) \left(-\left(\frac{1}{2} 
u^2_{j+\frac{1}{2}}\right)_t -
\varepsilon \frac{(\delta u_j)^2}{v_j} + \varepsilon \delta 
\left({\tau}_- u 
\frac{\delta u}{v}\right)_{j+\frac{1}{2}}\right) \nonumber\\
& & \qquad
+ 2 \varepsilon \frac{({\theta}_j -1) (\delta u_j)^2}{v_j} 
+ 
K\phi({\theta}_j)Z_j({\theta}_j-1)\left(2q+c_v{'}(Z_j)({\theta}_j+1)\right),
 \nonumber
\eea
which, using 
$({\theta}_j -1) (u^2_{j+\frac{1}{2}})_t 
=\left(({\theta}_j-1) u^2_{j+\frac{1}{2}}\right)_t
  - u^2_{j+\frac{1}{2}} \dot{\theta}_j$, yields that 
%2.36
\bea
& & \left(c_v(Z_j)({\theta}_j - 1)^2 + ({\theta}_j - 1) 
u^2_{j+\frac{1}{2}}\right)_t 
 +2 \lambda \frac{(\delta {\theta}_j)^2}{v_j} -
2 \lambda \delta \left(({\theta}_j -1) \frac{\delta 
{\theta}}{v} \right)_j
+ 2 \delta(({\theta}_j -1) u {\tau}_+ p)_j  \nonumber \\
& & \tab \tab 
=2 u_k p_{k+\frac{1}{2}} \delta {\theta}_k - 2 \varepsilon 
u_{j-\frac{1}{2}} 
 \frac{\delta u_j}{v_j} \delta {\theta}_{j - \frac{1}{2}} 
 + K \phi({\theta}_j) Z_j
({\theta}_j -1)\left(2 q 
+c_v{'}(Z_j)(\theta_j+1)\right)\nonumber \\
& & \tab \tab \quad
+ \frac{u^2_{j+\frac{1}{2}}}{c_v(Z_j)} 
\left(- p_j \delta u_j + \varepsilon \frac{(\delta 
u_j)^2}{v_j} + 
\lambda \delta \left(\frac{\delta \theta}{v}\right)_j 
+K\phi({\theta}_j)Z_j(q + c_v{'}(Z_j){\theta}_j)\right). 
\label{E:2.36}%*5  
\eea
Summing \eqref{E:2.36} over all $j$ and integrating over 
$[0,t]$, we
get
 %2.37
\bea
& & \sum_j \left(c_v(Z_j)({\theta}_j - 1)^2 
+ ({\theta}_j -1) u^2_{j+\frac{1}{2}}\right)(t)h\nonumber\\
&& =
\sum_j \left(c_v(Z_j)({\theta}_j - 1)^2 
+ ({\theta}_j -1) u^2_{j+\frac{1}{2}}\right)(0)h 
-2 \varepsilon \int_0^t \sum_j \left(u_{j-\frac{1}{2}} 
 \frac{\delta u_j}{v_j} \delta {\theta}_{j - 
\frac{1}{2}}\right) h ds 
\nonumber\\
& & \quad
- 2 \lambda \int_0^t \sum_j \frac{(\delta 
{\theta}_j)^2}{v_j} h ds + 2
\int_0^t \sum_j (u_k p_{k+\frac{1}{2}} \delta {\theta}_k)h 
ds 
-\int_0^t \sum_j \frac{1}{c_v(Z_j)}
 \left(u^2_{j+\frac{1}{2}} p_j \delta u_j\right) h ds 
\nonumber\\
&&\quad +\int_0^t\sum_j
  \left(u^2_{j+\frac{1}{2}}\frac{(\delta u_j)^2}{c_v(Z_j)v_j}\right)h ds
 +\int_0^t \sum_j \left( 
\frac{u^2_{j+\frac{1}{2}}}{c_v(Z_j)} 
 \delta\left(\frac{\delta \theta}{v}\right)_j\right) h ds 
\nonumber \\
& &\quad + \int_0^t \sum_j K \phi({\theta}_j) Z_j 
\left((\theta_j - 1)\left(2q+c_v'(Z_j)(\theta_j+1)\right)
+\frac{u^2_{j+\frac{1}{2}}}{c_v(Z_j)}(q+c_v{'}(Z_j){\theta}_j)\right)h ds.
\label{E:2.37}%*6
\eea
Now, take $L$ such that $L> (4\min_{0\le Z\le 
1}{c_v(Z)})^{-1}$
and consider the quantity
$$
w_j(t)=\left(c_v(Z_j)({\theta}_j-1)^2+({\theta}_j-1)u^2_{j+\frac{1}{2}} 
       + L u^4_{j+\frac{1}{2}}\right)^{1/2}.
$$
Then, multiplying \eqref{E:2.32} by $L$ adding to 
\eqref{E:2.36},
following the line of arguments as in \cite{CHT}, 
\cite{Ho1}, \cite{Ho2}, 
and noting that
$u_j u_{j-\frac{1}{2}} = u_j^2 - \frac{h}{2} u_j \delta 
u_j$,
we arrive at the estimate
\begin{eqnarray}
& & \!\!\!\!\! \sum_j w_j^2(t) h + \int_0^t \sum_k
(\delta {\theta}_k)^2hds+\int_0^t\sum_j( u_j^2 + 
\hat{u}^2_j) (\delta u_j)^2 h
ds \nonumber \\
& & \le M +M\int_0^t \sum_k \left(|u_k p_{k+ \frac{1}{2}}| 
|\delta e_k| 
+|u_k{\bar u}_{k+\frac{1}{2}}p_{k+\frac{1}{2}}\delta 
u_{k+\frac{1}{2}}|\right)
 hds \nonumber\\
&& + M\int_0^t\sum_j\left(h|u_j||\delta u_j|^3 
   + |u_{j-\frac{1}{2}}| |\delta u_j||\delta 
{\theta}_{j-\frac{1}{2}}|
 + (\theta_j^2+1 
+u^2_{j+\frac{1}{2}})(\theta_j+1)\phi(\theta_j)Z_j\right)hds,
\nonumber
\end{eqnarray}
where $\bar{u}_k = \frac{ u_{k+1} + u_k}{2}.$
The last relation implies (cf. \cite{CHT}, \cite{CHT2}, 
\cite{Ho1}, \cite{Ho2})
that
%2.38
\bea
&&\sum_j\left(c_v(Z_j)({\theta}_j-1)^2+({\theta}_j-1)u^2_{j+\frac{1}{2}} 
   + L u^4_{j+\frac{1}{2}}\right) h  \nonumber\\
&&\tab\tab\tab 
 +\int_0^t\left(\sum_k(\delta {\theta}_k)^2 h +\sum_j 
(u^2_j + \hat{u}_j^2) 
(\delta u_j)^2 h \right) ds \le M. \label{E:2.38}%*7
\eea
Finally, from Lemma \ref{L2.1},
\bea
\frac{1}{2} \sum_j Z_j^2 h + \int_0^t \sum_j  K 
\phi(\theta_j)Z_j^2 h ds =
\frac{1}{2} \sum_j Z_j^2(0)h. 
\nonumber
\eea
Adding the last two relations, we obtain the result.
\hfill $\blacksquare$


\medskip
Next we derive  some a-priori estimates, 
which provide essential information on the evolution 
of large jumps of discontinuities.
More specifically, the next lemmas   show that the
magnitude of any jump of
discontinuity of $v=v_j(\cdot)$ and $Z=Z_j(\cdot)$: 
\bea
[Z_{k}] =  Z_{k + \frac{1}{2}} - Z_{k-\frac{1}{2}}, \qquad
[v_{k}] =  v_{k + \frac{1}{2}} - v_{k-\frac{1}{2}}, 
\nonumber
\eea
at time $t$ can be controlled by the magnitude of jumps of
discontinuity of $v=v_j(\cdot)$ and $Z=Z_j(\cdot)$  at time 
$t=0.$
Without ambiguity, we denote 
$f_k=(f)_k=\frac{f_{k+\frac{1}{2}}+f_{k-\frac{1}{2}}}{2}$
for all other  variables except $u_k=u(x_k,t)$.
%%%%%%%%%%%%% lemma 2.6  %%%%%%%%%%
\medskip
\begin{lemma}\label{L2.6} 
For any $k \in \{0,1, \cdots, Q \},$
%2.39
\bea
&&\label{E:2.39}  \\ %zz
&&[Z_{k}(t)]=\exp\left\{-\int_0^t  K {\phi}_k (\tau)\, 
d\tau\right\}[Z_{k}(0)] 
+\int_0^t \exp \left\{-\int_s^t K {\phi}_k (\tau)\, 
d\tau\right\}
    Q_{k}(s)[{\theta}_{k}]\, ds, \nonumber
\eea
where
%2.40
\begin{equation}
Q_k(t)=-KZ_{k}(t)\int_0^1\phi{'}(\alpha{\theta}_{k+\frac{1}{2}}(t)
       +(1-\alpha){\theta}_{k-\frac{1}{2}}(t)) d \alpha.    
\label{E:2.40}%zz2
\end{equation}
\end{lemma}
\noindent
{\bf Proof:}
Equation \eqref{E:2.4} implies 
\bea
[Z_{k}]_t = -K [(\phi Z)_{k}] = -K{\phi}_{k}[Z_{k}] - K
Z_{k} [{\phi}_{k}].
\nonumber
\eea
Now, 
\bea
[{\phi}_{k}]
& = & 
\phi(\theta_{k+\frac{1}{2}})-\phi(\theta_{k-\frac{1}{2}})
=\left\{\int_0^1\phi{'}(\alpha{\theta}_{k+\frac{1}{2}}
     +(1-\alpha){\theta}_{k-\frac{1}{2}})d\alpha\right\} 
  [{\theta}_k]. \nonumber
\eea
Therefore,
\bea
\frac{d}{ds} 
\left\{\exp\left\{\int_0^s K {\phi}_k (\tau)d\tau \right\}
[Z_{k}] \right\} = 
\exp \left\{\int_0^s K{\phi}_k(\tau)d\tau\right\}Q_{k}[{\theta}_{k}]. 
\nonumber
\eea
Integrating the last relation on the interval $[0,t]$, we 
obtain 
\eqref{E:2.39}.
\hfill $\blacksquare$

%%%%%%%%%%%%%% lemma 2.7%%%%%%%%%%%%%%%%%
\begin{lemma}\label{L2.7}
There exists $M=M(m_0)>0,$  $m_0$ independent of $h$ and 
$N$, such that, 
for the distinguished discontinuities 
$ 0 < x_{k_1} <  x_{k_2} < \dots < x_{k_N} < 1, $
%\bea
%|[Z_{k_i}(t)]| \le M |[Z_{k_i}(0)]|. \label{na}
%\eea
%Away from the distinguished discontinuities,
%2.41
\bea
\sup_t\left(\sum_{k=k_i}|\delta Z_k|h
+\sum_{k\ne k_i}|\delta Z_k|^2h \right) \le M.
\label{E:2.41}%nana
\eea
\end{lemma}

\noindent
{\bf Proof:}
{\em Step 1.} We first prove that $\quad 
\sup_t\sum_{k=k_i}|\delta Z_k|h \le M$.
Notice that
%2.42
\bea
%&=&\frac{Z_{k+ \frac{1}{2}} - Z_{k - \frac{1}{2}}}{h} 
%\nonumber\\
%&=&\frac{1}{h}\left(Z_{k+\frac{1}{2}}(0)
% \exp\left\{-\int_0^t K\phi({\theta}_{k+\frac{1}{2}})ds\right\} 
%  -Z_{k-\frac{1}{2}}(0)
%  \exp\left\{-\int_0^t K\phi({\theta}_{k-\frac{1}{2}})d\tau\right\}\right) 
%   \nonumber \\
\delta Z_k 
&=&\frac{1}{h} \left(Z_{k+\frac{1}{2}}(0) - 
Z_{k-\frac{1}{2}}(0)\right) 
  \left(\exp\left\{-\int_0^t 
K\phi(\theta)ds\right\}\right)_k\nonumber\\
&& + \frac{1}{h} Z_{k}(0)
   \left(\exp\left\{-\int_0^t 
K\phi({\theta}_{k+\frac{1}{2}})ds\right\} 
   -\exp\left\{-\int_0^t 
K\phi({\theta}_{k-\frac{1}{2}})ds\right\}\right)
 \nonumber \\
&=&\delta Z_k(0)
 \left(\exp\left\{-\int_0^t
K{\phi}(\theta)ds\right\}\right)_k\nonumber \\
&& + Z_k(0)\exp\left\{-\int_0^t 
K{\phi(\tilde{\theta}_k)}ds\right\}
\left\{-\int_0^t  K {\phi}{'}(\tilde{\theta}_k) \delta{\theta}_k 
d \tau \right\}.
\label{E:2.42}%2.9
\eea
Therefore,
\bea
|\delta Z_k|
& & \!\!\!\!\!  \le |\delta Z_k(0)|
 +   Z_k(0) 
\int_0^t \mbox{exp}
  \left\{-\int_0^tK{\phi}(\tilde{\theta}_k)ds\right\}
   K |{\phi}{'}(\tilde{\theta}_k)| |\delta {\theta}_k| d \tau 
\nonumber\\
& & \!\!\!\!\! \le
|\delta Z_k(0)| +  Z_k(0)
 \int_0^t \mbox{exp}
  \left\{-\int_0^t K{\phi}(\tilde{\theta}_k)ds\right\}
 K |{\phi}{'}(\tilde{\theta}_k)|  \frac{|\delta 
{\theta}_k|}{\sqrt{{\theta}_{k+\frac{1}{2}}{\theta}_{k-\frac{1}{2}}}}
\sqrt{{\theta}_{k+\frac{1}{2}}{\theta}_{k-\frac{1}{2}}} d 
\tau \nonumber\\
& & \!\!\!\!\! \le  
|\delta Z_k(0)| +  Z_k(0)
\int_0^t \mbox{exp}
  \left\{-\int_0^t K{\phi}(\tilde{\theta}_k)ds\right\}
K |{\phi}{'}(\tilde{\theta}_k)|  \frac{|\delta 
{\theta}_k|}{\sqrt{{\theta}_{k+\frac{1}{2}}{\theta}_{k-\frac{1}{2}}}}
(\tilde{\theta}_k +  ch |\delta {\theta}_k|) d \tau. 
\nonumber
\eea
Summing over all $k$, we obtain
\bea
\sum_{k} |\delta Z_k|h 
\le && \!\!\!\!\!\!\! \sum_k|\delta Z_k(0)|h  
  + \sum_k Z_k(0) (I_{1_k} +  I_{2_k})h. \nonumber
\eea
Furthermore,
\bea
&& \!\!\!\!\!\!\!\!\!\!\! \sum_k Z_k(0) I_{1_k} h  
\nonumber\\
&& \!\!\!\!\!\!\!\!\!\!\!=\sum_k Z_k(0)
\int_0^t \exp\left\{-\int_0^t 
K\phi({\tilde{\theta}}_k)ds\right\}
{\phi}{'}(\tilde{\theta}_k)\tilde{\theta}_k
\frac{|\delta{\theta}_k|}
{\sqrt{{\theta}_{k+\frac{1}{2}}{\theta}_{k-\frac{1}{2}}}} 
 d\tau h \nonumber\\
&&\!\!\!\!\!\!\!\!\!\!\!\le \sum_k Z_k(0) \int_0^t 
\frac{|\delta {\theta}_k|^2}{{{\theta}_{k+\frac{1}{2}}}
  {{\theta}_{k-\frac{1}{2}}}}d\tau h
  + \sum_k Z_k(0) \int_0^t
\exp\left\{-2\int_0^t K\phi({\tilde{\theta}}_k)ds\right\} 
|{\phi}{'}(\tilde{\theta}_k) \tilde{\theta}_k|^2 d\tau h, 
  \nonumber 
\eea
while
\bea
\sum_k  Z_k(0) I_{2_k}h 
&\!\!\!\!\! = \!\!\!\!\! & \sum_k Z_k(0)
\int_0^t \exp\left\{-\int_0^t 
K\phi({\tilde{\theta}}_k)ds\right\}
|{\phi}{'}(\tilde{\theta}_k)|
\frac{|\delta 
{\theta}_k|^2}{\sqrt{{\theta}_{k+\frac{1}{2}}{\theta}_{k-\frac{1}{2}}}} 
h d\tau h \nonumber \\
&\!\!\!\!\! = \!\!\!\!\! & \sum_k Z_k(0) 
\int_0^t \exp\left\{-\int_0^t 
K\phi({\tilde{\theta}}_k)ds\right\}
|{\phi}{'}(\tilde{\theta}_k)|
\frac{|\delta 
{\theta}_k|^2}{{\theta}_{k+\frac{1}{2}}{\theta}_{k-\frac{1}{2}}} 
\sqrt{{\theta}_{k+\frac{1}{2}}{\theta}_{k-\frac{1}{2}}} h d 
\tau h. \nonumber
\eea
The result now follows by using  Lemma \ref{L2.2}.

\medskip
{\em Step 2.} Similarly, we now prove 
$\quad \sup_t\sum_{k=k_i}|\delta Z_k|^2h \le M$.

\bea
(\delta Z_k)^2 
&\le& M (\delta Z_k(0))^2 \left( \exp \left\{-2 \int_0^t K
 \phi(\theta) ds \right\} \right)_k \nonumber \\ 
&& + \,  M Z_k^2(0)\exp\left\{-2\int_0^t 
K\phi(\tilde{\theta}_k)ds\right\}
 \left\{\int_0^t K {{\phi}'}(\tilde{\theta}_k)
 \delta\theta_k ds\right\}^2\nonumber \\
& \le& M(\delta Z_k(0))^2+CZ_k^2(0) \int_0^t\exp
  \left\{-2\int_0^t K\phi(\tilde{\theta}_k)ds\right\} 
|{\phi}{'} (\tilde{\theta}_k)|^2 |\delta {\theta}_k|^2 ds. 
\nonumber
\eea
Thus,
\bea
& & \!\!\!\!\!\!\!\!\!\!\!\! \sum_k |\delta Z_k|^2 h  
\nonumber \\
& & \!\!\!\!\!\!\!\!\!\!\!\! \le  M \sum_k |\delta 
Z_k(0)|^2 h 
+ M\sum_kZ_k^2(0)  \int_0^t \exp
  \left\{-2\int_0^t K\phi(\tilde{\theta}_k)ds\right\} 
|{\phi}{'}(\tilde{\theta}_k)|^2 |\delta {\theta}_k|^2 d 
\tau h.\nonumber
\eea
Using Lemma \ref{L2.5}, we obtain the result.
\hfill $\blacksquare$
\medskip

For the study of the large jumps of discontinuities, 
we also need estimates on  $[v], [u_x],$ and $[{\theta}_x].$
The following lemma addresses this issue.
\smallskip
%%%%%%%%%%%%%%%%%%%% lemma 2.8%%%%%%%%%%%%%%%%%%%%
\begin{lemma}\label{L2.8} 
There exists $M=M(m_0)>0,$ $m_0$ independent of $h$ and $N$,
such that  the following estimates hold:

\noindent
{\rm (a)} For the distinguished discontinuities,
$0<x_{k_1}<x_{k_2}<\cdots <x_{k_N}<1$, 
\bea
\begin{cases}\label{E:2.43}
[\log{v_k}(t)]
={\mu}^{-1}_k(t) [\log{v_k}(0)] 
+ {\mu}^{-1}_k(t)\int_0^t{\mu}_k(s)(R_{k,h}(s)+\beta_k(s)[Z_k(s)])ds,\\ 
 \\
|[v_k(t)]| \le {\mu_k}^{-1}(t) |[v_k(0)]| + {\mu_k}^{-1}(t) 
\int_0^t {\mu}_k(s) {\beta}_k(s) [Z_k(s)] ds + Mh^{1/2}, \\
  \\
|[u_x(t)]|+|[{\theta}_x(t)]|
 \le M\min\left(\exp \{-M^{-1}t^{1/2}\}, 
        \;\sigma(t)^{-3/2}\exp\{-M^{-1} t\} \right),
\end{cases} 
\eea
where $\sigma(t)=\min(1,t)$ and
%2.44 
\bea
&&{\mu}_k(t)=\exp\left\{-\int_0^t {\alpha}_k(s) 
{\theta}_k(s) ds \right\},
\qquad
R_{k,h}(t)=\frac{a(Z)_k}{\varepsilon}[{\theta}_k]\left(\frac{1}{v}\right)_k 
          + \frac{h}{\varepsilon} \dot{u}_k , \nonumber\\
&&{\alpha}_k(t)=\frac{a(Z)_k\left[\frac{1}{v}\right]_k}{\varepsilon[\log{v_k}]},
\qquad \beta_k(t) = a(\tilde{Z})_k 
\left(\frac{\theta}{v}\right)_k;
\label{**B}
\eea

\medskip
\noindent
{\rm (b)} Away from the distinguished discontinuities,
\bea
\qquad\sup_t\sum_{k\ne k_i} (\delta v_k)^2 h + \int_0^t 
\sum_{k \not= k_i}
(1+ {\theta}_k) (\delta v_k)^2 h ds + \int_0^t \sum_j 
(\delta u_j)^2 h ds \le
M (1 + N h^{1/2}t). 
\label{5.1}
\eea
\end{lemma}

\noindent
{\bf Proof:}
\noindent
{\rm (a)}
Fix a jump point $k=k_i$. Then, for $w_k = \log{v_k}$,
\bea
[\varepsilon w_k]_t = \frac{\varepsilon \dot{v}_k}{v_k} 
=\frac{\varepsilon\delta u_k}{v_k} = [p_k] + h \dot{u}_k.
\label{**1}
%2.45
\eea
The momentum equation \eqref{E:2.2} yields
\bea
[p_k]
=a(Z)_k {\theta}_k \left[\frac{1}{v}\right]_k 
+ \frac{a(Z)_k}{\varepsilon} [{\theta}_k] 
\left(\frac{1}{v}\right)_k + 
a{'}(\tilde{Z}_k) [Z_k] \left(\frac{\theta}{v}\right)_k. 
\label{**2}
%2.47
\eea
Therefore,
\iffalse
\bea
[\dot{w}_k] = {\alpha}_k(t) {\theta}_k [w_k] +
\frac{a(Z)_k}{\varepsilon} [{\theta}_k] 
\left(\frac{1}{v}\right)_k
+a{'}(\tilde{Z}_k)[Z_k] \left(\frac{\theta}{v}\right)_k + 
\frac{h}{\varepsilon} \dot{u}_k,
\eea
\fi
\bea
[\dot{w}_k]= \alpha_k(t) {\theta}_k(t) [w_k] + \beta_k(t) 
[Z_k] + R_{k,h}(t),
%2.48 
\eea
where  ${\alpha}_k = {\alpha}_k(t)$, 
${\beta}_k={\beta}_k(t)$, and
$R_{k,h} = R_{k,h}(t)$ are given in \eqref{**B}
 with
$$a(Z)_k = \frac{a(Z(y_k+,t)) + a(Z(y_k-,t))}{2}.$$ 


Now let $\mu_k(t) = \exp \{- \int_0^t 
\alpha_k(s){\theta}_k(s)ds\}.$
Then
\bea
\frac{d}{dt}(\mu_k(t)[w_k]) = \mu_k(t)[\dot{w}_k] 
+\mu_k(t)\left(- \alpha_k(t){\theta}_k(t)\right)[w_k]
 = \mu_k(t) (R_{k,h}(t) + \beta_k(t) [Z_k]). 
\nonumber
\eea
Integrating with respect to $t$, we obtain
$$
[w_k(t)]=\mu_k^{-1}(t) [w_k(0)]  
+\mu_k^{-1}(t) \int_0^t \mu_k(s) \left(R_{k,h}(s) 
+ {\beta}_k(s) [Z_k(s)]\right) ds.
$$
By the mean-value theorem, we have 
\bea
[(\log{v})_k] = \log v_{k+\frac{1}{2}}-\log 
v_{k-\frac{1}{2}}
=\frac{1}{\tilde{v}_k} [v_k], \nonumber
\eea
which implies that
%2.49 
\bea
|[v_k]| = |\tilde{v}_k [w_k]| \le M |[w_k]|. 
\label{E:2.49}%5.8
\eea
Next, we estimate the quantity
$$
\left\{\int_0^t {\mu}_k(s) R_{k,h}(s)ds\right\}
$$
using the H\"{o}lder inequality, 
the estimates $\int_0^t |\delta {\theta}_k|^2h$ 
or $|\delta {\theta}_k| \le C/\sqrt{h}$, and the line of 
arguments presented in \cite{CHT} and \cite{CHT2} to obtain
$$
\mu_k^{-1}(t) \int_0^t {\mu}_k(s) R_{k,h}(s) ds \le M h^{1/2}.
$$
Therefore,
\bea
|[v_k(t)]| \le {\mu_k}^{-1}(t) |[v_k(0)]| + {\mu_k}^{-1}(t) 
\int_0^t {\mu}_k(s) {\beta}_k(s) [Z_k(s)] ds + Mh^{1/2}.\nonumber
\eea
\medskip

\noindent
{\rm (b)}
Let $x_{k_1} < x_{k_2}< \dots < x_{k_N}$ be distinguished 
nodes at which
discontinuities in $v(\cdot, t)$ are modeled and  
$[w_{k_i}]$ 
 the jump $w_{k_i + \frac{1}{2}} - w_{k_i - \frac{1}{2}}$ 
in a 
sequence $\{w_j\}$.
Here and in what follows, we denote
$$ 
\sum_k{\!'} w_k \equiv \sum_{k \not= k_i}w_k.
$$ 

\noindent
{\em Step 1.} We obtain an estimate for the quantity 
$\sum(\delta v_k)^2 h$.
Set $w=\log{v}.$ Then \eqref{E:2.2} implies that
$$ 
\varepsilon \delta \dot{w}_k = \dot{u}_k + \delta p_k.
$$ 
Multiplying the above relation by $\delta w_k$ and using 
\eqref{E:1.1.2},
we obtain
%2.50
\bea
& & \!\!\!\!\!\!\!\!\!\!\! \frac{\varepsilon}{2} \sum{\!'} 
(\delta w_k)^2 h
\left|_0^t \right.=\sum{\!'} u_k
\delta w_k h \left|_0^t \right. - \int_0^t \sum{\!'} u_k 
\delta \dot{w}_k h ds
+ \int_0^t \sum{\!'} \delta p_k \delta w_k h 
ds\label{E:2.50}.%**8
\eea
Since
\bea
\delta p_k 
=-a(Z)_k{\theta}_k\frac{\delta 
v_k}{v_{k+\frac{1}{2}}v_{k-\frac{1}{2}}}
 +\frac{a(Z)_k}{\varepsilon} \delta 
{\theta}_k\left(\frac{1}{v}\right)_k
 + a{'}(\tilde{Z}_k) \delta Z_k 
\left(\frac{\theta}{v}\right)_k, \nonumber
\eea
we have
\bea 
& &\!\!\!\!\!\!\!\!\!\!\! \sum{\!'}(\delta w_k)^2 
h+\int_0^t\sum{\!'}
 \frac{a(Z)_k{\theta}_k \delta w_k \delta v_k}
 {v_{k+\frac{1}{2}}v_{k-\frac{1}{2}}}hds  \nonumber \\
& &\!\!\!\!\!\le M+\frac{1}{\eta} \sum u_k^2(t)h 
           + \eta \sum{'} (\delta w_k)^2(t) h \nonumber \\
& & \!\!\!\!\!\,\,\, + \int_0^t \sum_j|\delta u_j 
\dot{w}_j|h ds 
    - \sum_i \int_0^t u_{k_i}[\dot{w}_{k_i}]ds 
 +M\int_0^t\sum{\!'}|\delta w_k|\left(a(Z)_k|\delta 
{\theta}_k| 
 +a{'}(\tilde{Z})_k {\theta}_k |\delta Z_k|\right) h ds, 
  \nonumber 
\eea
where $\eta$ is sufficiently small.
Therefore,
\bea
& & \!\!\!\!\!\!\!\!\!\! \sum{\!'}(\delta v_k)^2 h 
  + M\int_0^t\sum{\!'} {\theta}_k (\delta v_k)^2 h 
ds\nonumber\\
&&\le M 
+ M \int_0^t \sum_j (\delta u_j)^2 h 
-M'\sum_i\int_0^tu_{k_i}[\dot{w}_{k_i}]hds 
 \nonumber \\ 
& &\quad  + M \int_0^t\sum{\!'}|\delta v_k|
   (a(Z)_k|\delta {\theta}_k|+ 
a{'}(\tilde{Z}_k){\theta}_k|\delta Z_k|) hds, 
  \nonumber
\eea
where $M'$ is a constant which is similar to $M$ but may be
different.

The third term in the above relation can be controlled by 
using 
Lemmas \ref{L2.1}, \ref{L2.2}, and \ref{L2.4}.
Taking into consideration 
$$
[\dot{w}_{k_i}] 
={\alpha}_{k_i}(t) {\theta}_{k_i}(t)[w_{k_i}]
+{\beta}_{k_i}(t) [Z_{k_i}] + R_{k_i,h}(t),
$$ 
one has
$$
\sum_i \int_0^t u_{k_i} [\dot{w}_{k_i}] h ds 
=\sum_i\int_0^t u_{k_i}(s)
\left( {\alpha}_{k_i}(s) {\theta}_{k_i}(s) [w_{k_i}]+
{\beta}_{k_i}(s) [Z_{k_i}]+ R_{k_i,h}(s)\right)hds. 
$$
The first and third terms  can be handled by using the line 
of argument 
in \cite{CHT} and \cite{CHT2}. Here we estimate only the 
reacting term
%2.51
\bea
& & \!\!\!\!\! \int_0^t \sum_i u_{k_i}(s) 
{\beta}_{k_i}(s)[Z_{k_i}]
   h ds
=\int_0^t \sum_i u_{k_i}(s) a{'}(\tilde{Z})_{k_i}(s) 
 \left(\frac{{\theta}_{k_i}}{v}\right)(s) [Z_{k_i}]h ds 
\nonumber\\
& & \quad\le
M \int_0^t \|u\|_{L^{\infty}} \left(\sum_i  |\delta 
Z_{k_i}|h \right)ds  
\le M \sup_t \left(\sum_i |\delta Z_{k_i}| h \right)
\int_0^t \left(\sum_j (\delta u_j)^2  h \right)^{1/2} ds.
\eea
The quantity $\int_0^t \sum_j (\delta u_j)^2 h ds$ can
be estimated by following the line of argument in 
\cite{Ho1} and \cite{Ho2}.
The result now follows by using Lemma \ref{L2.7}.
\hfill $\blacksquare$

\medskip
Since the  initial data is discontinuous, we need to 
introduce 
an auxiliary function $\sigma=\sigma(t)=\min(t,1), t>0,$
which will serve as a weight for the following
regularity estimates.
The estimates derived in this section will be also crucial 
in the 
study of the large-time behavior of the solutions 
to \eqref{E:1.1.1}--\eqref{E:1.1.7} in Section \ref{S3}.
\noindent 
%%%%%%%%%%%%%%  Lemma 2.9%%%%%%
\begin{lemma}\label{L2.9}
There exists $M=M(m_0)$, $m_0$ independent of $h$ and $N$,
such that the following estimates hold:
\bea
& & \!\!\!\!\!\!\!\!\! \!\!\!{\rm (a)}\tab\tab
\sup_{t}\left(\! \sigma(t)\sum_j(\delta
u_j)^2(t)h \! \right) 
  +\!\int_0^t \sigma(s)\sum_k\dot{u}_k^2(s)h \;ds 
   \le M(1+N h^{1/2}t), \nonumber \\
&& \!\!\!\!\!\!\!\!\!\!\!\!  {\rm (b)} \tab\tab
\sup_{t}\left({\sigma(t)}\sum_k(\delta 
{\theta}_k)^2(t)h\right)
   +\int_0^{t}\sigma^2(s)\sum_j\dot{\theta}_j^2(s)h\;ds
   \le M (1 + N h^{1/2} {t}), \nonumber \\
& & \!\!\!\!\!\!\!\!\!\!\!\! {\rm
(c)}\tab\tab\sup_{t}\left({\sigma^2(t)} \left( \sum_k 
\dot{u}_k^2(t)  h
+ \sum_j (\delta u_j)^4(t)h\right) \right)   
 +\int_0^{t} 
\!{\sigma^2(s)}\sum_j(\delta\dot{u}_j)^2(s)h\;ds 
  \le M (1 + N  h^{1/2}{t}), \nonumber \\
& & \!\!\!\!\!\!\!\!\!\!\!\!
 {\rm 
(d)}\tab\tab\sup_t\left({\sigma}^3(t)\sum_j\dot{\theta}_j^2(t) 
h\right) 
+ \int_0^t {\sigma}^3(s) \sum_j 
\left(\frac{{(\delta\dot{\theta}_j)^2}}{v_j}
   \right)(s) h \; ds \le M (1+Nh^{1/2}t), \nonumber\\
& & \!\!\!\!\!\!\!\!\!\!\!\! {\rm (e)} \tab\tab 
  \sup_t\left(\sum_j\dot{Z}_j^2(t)h\right) 
  +\int_0^t 
\sum_j\left(K\phi({\theta}_j)\dot{Z}_j^2\right)(s)hds 
  \le M(1+Nh^{1/2}t). \nonumber 
\eea
\end{lemma}
\noindent

\noindent
{\bf Proof:}
The proofs of statements {\rm (a)--(d)} are 
standard and  follow similar lines of arguments as in 
\cite{CHT, Ho1, Ho2} 
with the aid of Lemmas \ref{L2.1}--\ref{L2.8}.
Special attention has been taken to accommodate the {\em 
large}
discontinuous initial data and the fact that $\gamma$ and 
$c_v$  
depend on $Z$ and therefore vary in space and time.
Result {\rm (e)} is a direct corollary of equation 
\eqref{E:2.4}, 
Lemma \ref{L2.1}, and the boundedness of $\phi(\theta)$. 
\hfill $\blacksquare$

%\noindent 
\begin{lemma}\label{L2.10}
There exists $M=M(m_0)>0$, $m_0$ independent of $h$ and 
$N$, such that, 
for all $t \in (0,T]$ and for distinguished discontinuities 
$0<x_{k_1}<x_{k_2}<\cdots<x_{k_N}<1$,
\bea
|[Z_{k_i}(t)]|\le M(|[Z_{k_i}(0)]|+ T^{1/2}h^{1/2}).
\nonumber
\eea
\end{lemma}

\medskip
\noindent
{\bf Proof:}
Lemma \ref{L2.6} yields that,
for the distinguished discontinuities,
\bea
& & |[Z_{k_i}(t)]| \le M |[Z_{k_i}(0)]| \nonumber \\
& & \tab \tab \tab \tab - K \int_0^t \exp 
  \left\{-\int_s^t K {\phi}_{k_i}(\tau)\,d\tau\right\}  
  \left( Z_{k_i}\int_0^t {\phi}{'}(\alpha{\theta}_{k_i-}
  +(1-\alpha) {\theta}_{k_i+})d\alpha \right)
    |[\theta_{k_i}(s)]|\, ds,\nonumber  
\eea
which implies
\bea
|[Z_{k_i}(t)]| \le M|[Z_{k_i}(0)]|+ M\int_0^t 
|[{\theta}_{k_i}]| ds. \nonumber
\eea
Using Lemma \ref{L2.9}, we have  
$$
|[{\theta}_{k_i}]|\le M\sqrt{\frac{h}{t}}.
$$
Therefore,
$$
|[Z_{k_i}(t)]|\le M(|[Z_{k_i}(0)]|+ T^{1/2} h^{1/2}).
$$
\hfill $\blacksquare$
%%%%%%%%%%%%%%%%%%%%%%%   section 3   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% section 3
\section {\bf Existence of Solutions:  Proofs of Theorems 
\label{S3}
\ref{T:1.1} and \ref{T:1.2}}

In this section we apply the estimates of Section 2 to 
prove the
existence results, Theorems \ref{T:1.1} and \ref{T:1.2}.  
We shall
prove Theorem \ref{T:1.2} first, obtaining in particular 
solutions with
smooth initial data.  We then construct the more general 
solutions of
Theorem \ref{T:1.1} as limits of these smooth solutions.  
The proof of
Theorem \ref{T:1.2}, which is quite standard, is nearly 
identical to
that of related constructions in the literature, and 
therefore will be
just sketched.  The proof of Theorem \ref{T:1.1}, on the 
other hand,
involves an interesting argument establishing strong 
convergence of
densities in the absence of regularity, and will therefore 
be presented
in detail.  
\bigskip

{\noindent \bf Proof of Theorem \ref{T:1.2}:}  The 
construction of solutions here
is nearly identical to that carried out in \cite{Ho1,Ho2}, 
and in
earlier work, for non-reacting, viscous flows.  Briefly, 
initial data
$(v_j(0), u_k(0), \theta_j(0), Z_j(0))$ for the 
semidiscrete system
(\ref{E:2.1})--(\ref{E:2.4}) is generated from the initial 
data
$(v_0,u_0,\theta_0,Z_0)$ in such a way that the constant 
$C_0$ in
(\ref{E:2.6})--(\ref{E:2.7}) is bounded by the constant 
$C_0$ in the
hypothesis of Theorem~\ref{T:1.2}.  The results of Section 2 
then apply
to show that the system (\ref{E:2.1})--(\ref{E:2.4}) has a 
global solution
$(v_j(t), u_k(t), \theta_j(t),Z_j(t))$ defined for all $t$, 
and
satisfying all the conclusions of Lemmas 
\ref{L2.1}--\ref{L2.10}.
These mesh functions are then extended to approximate 
solutions
$(v^h(t),u^h(t),\theta^h(t),Z^h(t))$, which then satisfy 
all the bounds
in (\ref{E:1.2.6})--(\ref{E:1.3.6}), independently of $h$ 
(modulo terms
which are $0(h)$).  The solution $(v,u,\theta,Z)$ of 
Theorem \ref{T:1.2} is then obtained as the strong limit as 
$h_i\to 0$ for a
suitable sequence $\{h_i\}$.  The details, which are 
routine, are nearly
identical to those of \cite{Ho1}, Section 3, and therefore 
are omitted.
\hfill $\blacksquare$

\bigskip
{\noindent \bf Proof of Theorem \ref{T:1.1}:}  First we 
construct smooth
approximations 
$(v^\delta_0,u_0^\delta,\theta_0^\delta,Z_0^\delta)$ to
the initial data $(v_0,u_0,\theta_0,Z_0)$ of Theorem 
\ref{T:1.1} by
mollifying appropriate extensions of the data with a 
smoothing kernel of
width $\delta$.  We then apply the result of Theorem 
\ref{T:1.2} to
obtain global solutions 
$(v^\delta,u^\delta,\theta^\delta,Z^\delta)$ of
\eqref{E:1.1.1} with initial data 
$(v^\delta_0,u^\delta_0,\theta^\delta_0,
Z^\delta_0)$ which satisfy the bounds in 
\eqref{E:1.2.6}--\eqref{E:1.2.7}
independently of $\delta$.  We shall obtain the solution 
$(v,u,\theta,Z)$ of Theorem \ref{T:1.1} as
the limit, in  an appropriate sense, of these smooth 
solutions.

First, the uniform bounds in \eqref{E:1.2.7} for 
$\displaystyle \int_0^1
(u^\delta_x)^2(x,t)dx$ for $t>0$ and for $\displaystyle 
\int_\tau^T\int^1_0 
(u_t^\delta)^2dxdt$ for $\tau>0$ show that $u^\delta$ 
are
uniformly H${\rm \ddot{o}}$lder continuous on compact sets 
in
$[0,1]\times (0,\infty)$.  There is therefore a subsequence 
$\delta_k \to 0$ such that $u^{\delta_k}$ converges uniformly, say to 
$u$,
uniformly on these sets.  Then $u^\delta_x(\cdot,t)\to 
u_x(\cdot,t)$ in
$\mathcal{D}^\prime ((0,1))$, and since $\{u^\delta _x 
(\cdot,t)\}$ is
weakly compact in $L^2([0,1])$, again by the uniform bounds 
in
(\ref{E:1.2.6}), $u^\delta_x(\cdot,t)\rightharpoonup 
u_x(\cdot,t)$ weakly in
$L^2([0,1])$.  In particular, $u(\cdot,t) \in H^1 ([0,1])$ 
for
$t>0$.  Similar arguments apply to the sequence 
$\theta^\delta$.
We have therefore identified a sequence $\delta_k = \delta \to 0$ 
such that
%3.1
\begin{equation}\label{E:3.1}
\begin{cases}
u^\delta \to u, \ \theta^\delta \to \theta  \ 
\text{uniformly on compact
sets in $[0,1] \times (0,\infty)$};\\
%
u^\delta _x(\cdot,t)\rhup u_x(\cdot,t), \ \theta^\delta 
_x(\cdot,t)\rhup 
\theta_x(\cdot,t) \ \text{weakly in $L^2([0,1])$}  \ 
\text{for $t>0$};\\
%
u^\delta_x \rhup u_x, \ \theta ^\delta_x\rhup \theta_x  \ 
\text{weakly in
$L^2([0,1] \times [0,T])$} \ \text{for $T>0$}.
\end{cases}
\end{equation}
%
Next, we show that, for this same sequence $\{\delta_k\}$,
$Z^{\delta_k}(\cdot,t)$ converges, say to $Z(\cdot,t)$, 
strongly in
$L^p([0,1])$ for all $t\ge0$ and all $p\in[1,\infty)$.  To 
prove this,
we solve the fourth equation in (\ref{E:1.1.1}) for
$Z^{\delta_j}-Z^{\delta_k}$ to obtain that
%No Number
\begin{equation*}
\begin{split}
Z^{\delta_j}(x,t)&-Z^{\delta_k}(x,t)=e^{-K \int^t_0\phi
(\theta^{\delta _j}(x,s))ds}
\big(Z^{\delta_j}(x,0)-Z^{\delta_k}(x,0)\big) \\
&+\displaystyle \int^t_0 e^{-K \int^t_s \phi
(\theta^{\delta_j}(x,\tau))d\tau}
\big(\phi(\theta^{\delta_k}(x,s))-\phi(\theta^{\delta_j}(x,s))\big)
Z^{\delta_k}(x,s)ds,
\end{split}
\end{equation*}
%
so that
%3.2
\begin{equation}\label{E:3.2}
\begin{split}
\|Z^{\delta_j}(\cdot,t)-Z^{\delta_k}(\cdot,t)\|_{L^1} \le \ 
&\|Z^{\delta_j}(\cdot,0)-Z^{\delta_k}(\cdot,0)\|_{L^1}\\
&+C \int^t_0
\|\theta^{\delta_j}(\cdot,s)-\theta^{\delta_k}(\cdot,s)\|_{L^1}ds.
\end{split}
\end{equation}
%
The first term on the right here goes to zero as 
$\delta_j,\delta_k\to
0$ because $Z^\delta(\cdot,0) \to Z_0$ in $L^1$, and the 
second term is
bounded by
%3.3
\begin{equation}\label{E:3.3}
\begin{split}
\int^t_\tau \int^1_0 
|&\theta^{\delta_j}-\theta^{\delta_k}|dxds + C 
\int^\tau_0 \int^1_0 (\theta^{\delta_j} + \theta^{\delta_k}) dxds\\
&\le \int^t_\tau \int^1_0 
|\theta^{\delta_j}-\theta^{\delta_k}|dxds + C
\tau
\end{split}
\end{equation}
%
by the uniform bound in (\ref{E:1.2.7}) for $\mathcal{D}$. 
This together
with (\ref{E:3.1}) then shows that the second term on the 
right side of
(\ref{E:3.2}) approaches zero as $\delta_j,\delta_k\to 0$.  
This proves
that $Z^\delta (\cdot,t) \to Z(\cdot,t)$ in $L^1$; 
convergence in
$L^p$ for $p<\infty$ then follows from the fact that $0 \le 
Z^\delta\le
1$.

Finally, we claim that there is a further subsequence 
$\delta_j \to 0$
such that $v^{\delta_j}(\cdot,t)$ converges, say to 
$v(\cdot,t)$,
strongly in $L^p$ for all $t\ge 0$ and all 
$p\in[1,\infty)$.  To prove this, we first define 
%3.4
\begin{equation}\label{E:3.4}
F^\delta = \frac{\varepsilon
u^\delta_x}{v^\delta}-p(v^\delta,\theta^\delta,Z^\delta).
\end{equation}
%
Then $F^\delta_x = u^\delta_t$, so that, by (\ref{E:1.2.6}),
%NO Number
\begin{equation*}
\|F^\delta_x(\cdot,t)\|_{L^2}\le C(\tau,T), 
\qquad  \tau \le t \le T,
\end{equation*}
%
and
%
%NO Number
\begin{equation*}
|F^\delta_t | \le C(|u^\delta_{xt}| + |u^\delta_x |^2 + 
|u^\delta_x | +
|\theta^\delta_x | + |Z^\delta |),
\end{equation*}
%
so that, again by \eqref{E:1.2.6},
%3.5
\begin{equation}\label{E:3.5}
\int^T_\tau \int (F^\delta_\tau)^2 dxdt \le C(\tau, T), \quad 
0<\tau<T.
\end{equation}
%
Thus, $\{F^\delta\}$ is uniformly bounded in $H^1 ([0,1]
\times[\tau,T])$, and there is a further subsequence 
$\delta_j\to 0$
such that $F^{\delta_j}\to F$ strongly in $L^2([0,1] \times 
[1/k,k])$
for all $k\in \ze_+$.  Note also that
%no #
\begin{equation*}
\int^\tau_0 \int^1_0 |F^\delta|dxdt \le C\int^\tau_0
\int^1_0(|u^\delta_x| + \theta^\delta)dxdt \le C\tau^{1/2}.
\end{equation*}
%
It therefore follows that
%3.6
\begin{equation}\label{E:3.6}
F^{\delta_j} \to F \ \text{in}  \ L^1([0,1] \times [0,T]) \ 
\quad \text{for} \ T>0
\end{equation}
%
(the limit $F$ could be identified in terms of 
$v,\theta,Z$, but the
present argument does not require it).  Now fix $j$ and 
$k$, let
$v^{\delta_j} = v^j, v^{\delta_k} = v^k$, etc., and let 
$L=\log v$.
Then, from the first equation in \eqref{E:1.1.1} and from the
definition of $F$,
%3.7
\begin{equation}\label{E:3.7}
\varepsilon \frac{\partial}{\partial t}(L^j-L^k) = 
(F^j-F^k) + (p^j-p^k).
\end{equation}
%
Since $p=a(Z)\theta/v$, we have
%N0 #
\begin{equation*}
p^j-p^k = a^k\theta^k 
\Bigg(\frac{(v^j)^{-1}-(v^k)^{-1}}{L^j-L^k}\Bigg)
(L^j-L^k) + O (|\theta^j||Z^j-Z^k| + |\theta^j-\theta^k|).
\end{equation*}
%
Letting $-\alpha$ denote the coefficient of $L^j-L^k$ on 
the right, so
that $\alpha \ge 0$, and substituting into \eqref{E:3.7}, 
we then
obtain that
%
\begin{equation*}
\varepsilon \frac{\partial}{\partial t} (L^j - L^k)= 
-\alpha(L^j-L^k) + 
(F^j-F^k)+O(|\theta^j||Z^j-Z^k| + |\theta^j-\theta^k|),
\end{equation*}
%
and therefore that
%3.8
\begin{equation}\label{E:3.8}
\begin{split}
\|(L^j-L^k)(\cdot,t)\|_{L^1}\le & \ 
\|(L^j-L^k)(\cdot,0)\|_{L^1}\\
+&C\int^t_0 \big(\|(F^j-F^k)(\cdot,s)\|_{L^1} +
 \|\theta^j(\cdot,s)\|_{L^\infty}\|(Z^j-Z^k)(\cdot,s)\|_{L^1}\\
&\qquad \qquad \qquad \qquad \qquad \qquad \qquad \
+\|(\theta^j-\theta^k) (\cdot,s)\|_{L^1}\big)ds.
\end{split}
\end{equation}
%
The first term on the right here approaches zero as $j,k\to 
\infty$
because $v^\delta(\cdot,0)\to v(\cdot,0)$ in $L^1$ as 
$\delta\to 0,$
and the first and third terms in the time integral have 
already been
shown above in \eqref{E:3.2}, \eqref{E:3.3}, and 
\eqref{E:3.6} to approach zero.  
For the remaining term, we apply the bounds in
\eqref{E:1.2.6} again to obtain, for $t>0$ fixed,
%3.9
\begin{equation}\label{E:3.9}
\|\theta^\delta\|_{L^\infty} 
\le \Big(\int_0^1(\theta^\delta)^2 dx\Big)^{1/2} 
+ \Big(\int_0^1(\theta^\delta)^2 dx\Big)^{1/4}
\Big(\int_0^1(\theta^\delta_x)^2 dx\Big)^{1/4}
\le C(1+t^{-1/2}).
\end{equation}
%
The second term in the time integral in \eqref{E:3.8} 
therefore goes to
zero by the $L^1$ convergence $Z^j(\cdot,t)\to Z(\cdot,t)$ 
and by the
dominated convergence theorem.  It therefore follows from 
\eqref{E:3.8}
that $\{v^{\delta_j}(\cdot,t)\}$ is strongly Cauchy in 
$L^1([0,1])$,
hence in $L^p ([0,1])$ for all $p\in[1,\infty)$, and for 
each $t\ge
0$.

It is now routine to check that $(v,u,\theta, Z)$ is indeed 
a weak
solution of \eqref{E:1.1.1} with initial data 
$(v_0,u_0,\theta_0,Z_0)$,
and that $(v,u,\theta,Z)$ inherits all the bounds 
 \eqref{E:1.2.6} and \eqref{E:1.2.7} in Theorem 
\ref{T:1.1} from the
smooth solutions 
$(v^\delta,u^\delta,\theta^\delta,Z^\delta)$.  The
results  in
\eqref{E:1.2.55} concerning regularity in time are then 
easily derived
from the bounds in \eqref{E:1.2.6}, \eqref{E:1.2.7}, and the weak 
forms of the equations \eqref{E:1.1.1}.
\hfill $\blacksquare$
%%%%%%%%%%%%%%%%%%%%%%% section 4   
%%%%%%%%%%%%%%%%%%%%%%%%%
\section {\bf Large-Time Behavior}  \label{S4}

In this section we prove the large-time behavior result, 
Theorem
\ref{T:1.3}, for a weak solution of \eqref{E:1.1.1} which 
satisfies the
bounds in \eqref{E:1.2.6} and \eqref{E:1.2.7} 
with a constant $M$ which is now assumed
to be independent of time.  As indicated in the statement 
of Theorem
1.1, $M$ is indeed independent of time when $c_1=c_2$ and
$\gamma_1=\gamma_2$.  The results of Theorem \ref{T:1.3} do 
not
otherwise depend on these conditions, however.  We also 
prove
Theorem~1.4 giving necessary conditions and sufficient 
conditions
for complete combustion in the time-asymptotic solutions.
\bigskip

\noindent
{\bf Proof of Theorem 1.3:}  \  
We fix a weak solution $(v,u,\theta,Z)$ satisfying the
conclusions of Theorem~\ref{T:1.1}, and we assume that  
$M$ in 
\eqref{E:1.2.6} and \eqref{E:1.2.7} is independent of time.  
The bound
%No #
\begin{equation*}
\int^\infty_1 \int (u_x^2 + u_{xt}^2)dxdt \le M
\end{equation*}
%
then shows that the function 
$\displaystyle \int_0^1 (u_x)^2(x,t)dx$ is in
$(L^1\cap BV) ([1,\infty))$, and consequently approaches 
zero as $t\to
\infty$.  That is, $u(\cdot,t)\to 0$ in $H^1 ([0,1])$ as 
$t\to \infty$.

To describe the large-time behavior of $Z$, we again let
$(v^\delta,u^\delta,\theta^\delta,Z^\delta)$ be the smooth 
solution of
\eqref{E:1.1.1} corresponding to mollified initial data, so 
that, for a
suitable sequence $\delta\to 0, Z^\delta (\cdot,t)\to 
Z(\cdot, t)$ in
$L^p([0,1])$ for $p\in[1,\infty)$, and $\theta^\delta\to 
\theta$
uniformly on compact sets in $[0,1] \times (0,\infty)$, 
just as in
Section \ref{S3}.  Then, from \eqref{E:1.1.1}, 
%No #
\begin{equation*}
Z^\delta(x,t)=Z^\delta(x,0)e^{-K\int_0^t 
\phi(\theta^\delta(x,s))ds}.
\end{equation*}
%
The integral in the exponent here converges to $\int^t_0
\phi(\theta(x,s))ds$ as $\delta\to 0$ because $\phi$ is 
Lipschitz
and $\|\theta^\delta(\cdot,s)\|_{L^\infty} \le \ 
C(1+s^{-1/2})$ (see
\eqref{E:3.9}). 
Thus, if $E$ is the set of Lebesgue points of $Z_0$, then, 
for a particular representative $Z(\cdot,t)$,
%4.1
\begin{equation}\label{E:4.1}
Z(x,t) = Z_0(x)e^{-K\int^t_0 \phi(\theta (x,s))ds}, \quad x\in E.
\end{equation}
%
$Z(x,\cdot)$ is therefore a decreasing function of $t$ for 
$x\in E$, and so
converges as $t\to \infty$, say to $Z_\infty (x)$.

To describe the large-time behavior of $\theta$, we apply 
the bounds
%
\begin{equation*}
\int^\infty_1 \int^1_0 (\theta_x^2 + \theta^2_{xt})dxdt \le 
M
\end{equation*}
%
to conclude that $\theta(\cdot,t) - \bar{\theta}(t)\to0$ in 
$H^1([0,1])$
as $t\to\infty$, where
%
\begin{equation*}
\bar{\theta}(t) = \int^1_0 \theta(x,t)dx.
\end{equation*}
%
On the other hand, the conservation of energy,
%4.2
\begin{equation}\label{E:4.2}
\int^1_0\big(c(Z(x,\cdot))\theta(x,\cdot) + \textstyle 
\frac{1}{2} u(x,\cdot)^2 +
qZ(x,\cdot)\big)dx  \ \Big|^t_0=0
\end{equation}
%
holds for smooth solutions, and therefore for the weak 
solutions of
Theorem \ref{T:1.1} as well, by the strong convergence of 
smooth
solutions described in Section \ref{S3}.  Taking the limit 
as $t\to
\infty$ in \eqref{E:4.2}, we therefore obtain that
%No #
\begin{equation*}
\left(\int^1_0 c(Z_\infty(x))dx\right)\, \bar{\theta}(t) 
+ q\int^1_0 Z_\infty(x)dx \to E_0,
\end{equation*}
%
where $E_0$ is the total initial energy, so that
%No #
\begin{equation*}
\bar{\theta}(t)\to  \frac{E_0-q\displaystyle 
\int^1_0Z_\infty
(x)dx}{\displaystyle \int^1_0 c(Z_\infty (x))dx} \ ,
\end{equation*}
% 
and therefore 
%4.3
\begin{equation}\label{E:4.3}
\theta(\cdot,t)\to \theta_\infty \equiv \displaystyle
\frac{E_0-q\displaystyle\int^1_0 Z_\infty(x) 
dx}{\displaystyle\int_0^1 c(Z_\infty(x))dx}
\end{equation}
%
in $H^1([0,1])$ as $t\to\infty$.  It then follows that
%No #
\begin{equation*}
e(x,t) = c(Z(x,t))\theta(x,t) \to 
c(Z_\infty(x))\theta_\infty
\end{equation*}
%
as $t\to\infty$, for $x\in E$.

Before deriving the large-time behavior of $v$, we examine 
that of the
quantity $F=\displaystyle \frac{\varepsilon u_x}{v}-p$ 
introduced in Section \ref{S3}.  
First, since $F_x=u_t$, we have from \eqref{E:1.2.6} that
%No #
\begin{equation*}
\displaystyle \int^\infty_1 \displaystyle \int_0^1 F^2_x dxdt\le M.
\end{equation*}
%
Also,
%no #
\begin{equation*}
\begin{split}
\displaystyle \int^\infty_1 \big|\frac{d}{dt} 
\displaystyle &\int^1_0 F_x^2 dx \big| dt = 2\displaystyle 
\int^\infty_1
\big|\displaystyle \int_0^1 u_{xt}F_t dx \big| dt\\
&\le C\displaystyle \int^\infty_1 \displaystyle \int_0^1
|u_{xt}|\big(|u_{xt}|+|u_x|^2 
+|Z_t|+|u_x|+|\theta_t|\big)dxdt \le M
\end{split}
\end{equation*}
%
by routine estimates based on 
\eqref{E:1.2.6} and \eqref{E:1.2.7}.  Thus
$\displaystyle \int_0^1 F_x(x,t)^2 dx \in (L^1 \cap 
BV)([1,\infty))$  (actually, the
above estimates are carried out for the smooth solutions
$(v^\delta,u^\delta,\theta^\delta,Z^\delta)$ of Section 
\ref{S3}, but
the conclusion persists in the limit as $\delta \to 0$).  
Thus, for the
weak solution under consideration, $\displaystyle \int^1_0 
F_x(x,t)^2 dx \to 0$, and
%4.4
\begin{equation}\label{E:4.4}
F(\cdot,t)-\int^1_0 F(x,t)dx \to 0 \quad \text{in} \ H^1([0,1])
\end{equation}
%
as $t\to \infty$.

To derive the large-time behavior of $v,$ we define
%No #
\begin{equation*}
\bar{v}_0=\int_0^1 v_0(x)dx = \int_0^1 v(x,t)dx
\end{equation*}
%
and
%no #
\begin{equation*}
v_\infty (x)=  \frac{\gamma(Z_\infty(x))-1}{\displaystyle 
\int^1_0
(\gamma(Z_\infty(x))-1)dx} \bar{v}_0 \ ,
\end{equation*}
%
and we show that $v(\cdot,t)\to v_\infty$ as $t\to \infty$ 
in $L^p$ for
$p=2$, hence for all $p\in[1,\infty)$.  We denote by $o(1)$ 
any term
which approaches zero in $L^2$ as $t\to \infty$.  Then, from
\eqref{E:4.4}, the definition of $F$, and the fact that
$u_x(\cdot,t)=o(1)$,
%No #
\begin{equation*}
\begin{split}
o(1)&= F(\cdot,t) -\int^1_0 F(x,t) dx\\
%
&= o(1) + \int^1_0 \Bigg(
\frac{\gamma(Z(x,t))-1}{v(x,t)}\Bigg)\theta(x,t)dx - \Bigg(
\frac{\gamma(Z(\cdot,t))-1}{v(\cdot,t)}\Bigg)\theta(\cdot,t
)\\
%
	&= o(1) + \int^1_0 \Bigg(\Bigg(\frac {\gamma(Z)
	-\gamma(Z_\infty)}{v}\Bigg)\theta + 
\Bigg(\frac{\gamma(Z_\infty)-1}
	{v}\Bigg)(\theta-\theta_\infty) + \Bigg(
	\frac{\gamma(Z_\infty)-1}{v}\Bigg)\theta_\infty \Bigg)dx\\
%
	&\qquad \qquad \qquad-\Bigg(\Bigg(\frac {\gamma(Z)
	-\gamma(Z_\infty)}{v}\Bigg)\theta + 
\Bigg(\frac{\gamma(Z_\infty)-1}
	{v}\Bigg)(\theta-\theta_\infty) + \Bigg(
	\frac{\gamma(Z_\infty)-1}{v}\Bigg)\theta_\infty \Bigg)\\
%
	&=o(1) + \Bigg( \int^1_0\Bigg( 
	\frac{\gamma(Z_\infty(x))-1}{v(x,t)}\Bigg)dx - 
	\frac{\gamma(Z_\infty(\cdot))-1}{v(\cdot,t)}\Bigg)\theta_\infty \ ,\\
\end{split}
\end{equation*}
%
since $Z(\cdot,t)\to Z_\infty$ in $L^1$ and $\theta\to 
\theta_\infty$ in
$H^1$.  Thus
%4.5
\begin{equation}\label{E:4.5}
v(\cdot,t) - \displaystyle \frac{\gamma(Z_\infty 
(\cdot))-1}{\displaystyle \int^1_0
 \Bigg(\frac {\gamma(Z_\infty(x))-1}{v(x,t)}\Bigg)dx}  
\quad \ \to 0
\end{equation}
%
in $L^2$ as $t\to \infty$.  Integrating with respect to 
$x$, we find that
that
%No #
\begin{equation*}
\int^1_0 \Bigg( 
\frac{\gamma(Z_\infty(x))-1}{v(x,t)}\Bigg)dx \to
\frac{1}{\bar{v}_0} \int^1_0 (Z_\infty(x)-1\Big) dx,
\end{equation*}
%
and, then substituting back into \eqref{E:4.5}, that
$v(\cdot,t)\to v_\infty$ in $L^2$ as $t\to \infty$.
\hfill $\blacksquare$
%
\bigskip

{\bf \noindent Proof of Theorem \ref{T:1.4}:}  First assume 
that $Z_\infty \equiv
0$.  Then $\theta_\infty\ge \theta_i$, for otherwise, since
$\theta(\cdot,t)\to \theta_\infty$ uniformly in $x$, there 
would be a
time $T$ such that $\theta(x,t) \le \theta_i$ and
$\phi(\theta(x,t))=0$ for all $x$ and all $t\ge T$.  It 
would then
follow from \eqref{E:4.1} that $Z_\infty >0$ on a set of 
positive
measure.  Thus $\theta_\infty\ge \theta_i$, and so by the 
conservation of mass \eqref{E:4.2},
%No #
\begin{equation*}
E_0 = \int^1_0 \Big(c(Z)\theta + \textstyle\frac{1}{2}u^2 
+qZ \Big)dx
\to c(0)\theta_\infty=c_2\theta_\infty \ge c_2\theta_i,
\end{equation*}
%
so that $E_0\ge c_2\theta_i$.

To prove (b) of Theorem \ref{T:1.4}, we first 
observe from 
\eqref{E:4.3}
that
%4.6
\begin{equation}\label{E:4.6}
\theta_\infty> \theta_i \iff E_0> c_2\theta_i + 
[(c_1-c_2)\theta_i + q]\int Z_\infty dx.
\end{equation}
%
Now, $E_0> c_2\theta_i$ by hypothesis, so that in the case 
that the
term in the brackets above is nonpositive, the contingency 
on the right hand side
of \eqref{E:4.6} holds, and $\theta_\infty > \theta_i$.  In 
the
other case, that 
%4.7
\begin{equation}\label{E:4.7}
(c_1-c_2)\theta_i +q>0,
\end{equation}
%
we obtain from \eqref{E:1.4.7} that
%No #
\begin{equation*}
\begin{split}
E_0& > c_2\theta_i + ((c_1-c_2) \theta_i+q)\int_0^1 Z_0 dx\\
&\ge c_2\theta_i + ((c_1-c_2) \theta_i+q)\int_0^1 Z_\infty 
dx
\end{split}
\end{equation*}
%
by \eqref{E:4.7} and the fact that $\displaystyle \int^1_0 
Z(x,t)dx$ is decreasing in time. 
The fact \eqref{E:4.6} then applies to show that 
$\theta_\infty> \theta_i$
in this case as well. It then follows from \eqref{E:4.1} 
that
$Z_\infty=0, \ a.e.$
\hfill $\blacksquare$
%%*** section below copied from article.tex

\section*{Appendix: Existence of Entropies}

In this section we establish the existence of a physical 
entropy for our 
system \eqref{E:1.1.1}. Let us start with a system of the 
form
%5.1
\bea
\begin{cases}\label{E:5.1}
v_t - u_x   = 0, \\
u_t - {\sigma}_x = 0, \\
\left(e+\frac{u^2}{2} +qZ \right)_t - (u \sigma)_x = Q_x,  
\\ 
Z_t + K\phi Z = 0.
\end{cases}%sy
\eea
Here $v, u, \theta, e$ are described as before, 
while $\sigma$ and $Q$ denote the stress and the heat flux, 
respectively.
In the process of establishing the presence of a physical 
entropy for 
system \eqref{E:5.1}, we will also specify what some 
appropriate choices
are for the stress ${\sigma}$ and the heat flux $Q.$
Here the internal energy, stress, and heat flux are 
determined through the constitutive relations
%5.2
\bea
\begin{cases}\label{E:5.2}
e = \hat{e}(v, \theta, {\theta}_x, Z),  \\  
{\sigma}=-\hat{p}(v, \theta, {\theta}_x, Z, Z_x) + 
\frac{\varepsilon u_x}{v},  \\
Q = \hat{Q} (v, \theta, u_x, {\theta}_x),  \\ 
\end{cases}
\eea
%
while ${\phi}={\phi}({\theta}).$
The constitutive variables are subject to restrictions 
arising from
the second law of  thermodynamics.
We seek a physical entropy $\eta$ for which the 
Clausius-Duhem
inequality 
%5.3 
\bea \label{E:5.3}
{\eta}_t-\left(\frac{Q}{{\theta}}\right)_x- K\frac{q\phi 
Z}{\theta} \ge 0
\eea
%
is satisfied.
Using \eqref{E:5.3}, the momentum equation yields 
%No #
\bea
e_t + u u_t - u {\sigma}_x - u_x \sigma
=Q_x+ Kq\phi Z\le \theta {\eta}_t + \frac{Q 
{\theta}_x}{\theta},\nonumber
\eea
%
or
%5.4
\bea \label{E:5.4}
e_t - u_x \sigma  \le \theta {\eta}_t + \frac{Q 
{\theta}_x}{\theta}. 
\eea
In terms of the Helmholtz free energy $ {\Psi} = e - \theta 
\eta $,
the last relation \eqref{E:5.4} can be rewritten as
%5.5
\bea \label{E:5.5}
{\Psi}_t + \eta {\theta}_t + p u_x - \frac{\varepsilon 
u_x^2}{v}
 - \frac{Q {\theta}_x}{\theta} \le 0. 
\eea
%
Assuming that
%no #
$$ 
\Psi = \hat{\Psi}(v, \theta, u_x, {\theta}_x, Z),
$$
then \eqref{E:5.5} implies that
%5.6
\bea \label{E:5.6}
{\Psi}_v v_t + {\Psi}_{\theta} {\theta}_t + {\Psi}_Z Z_t +
{\Psi}_{{\theta}_x} {\theta}_{xt} + {\Psi}_{u_x} u_{xt} +
\eta {\theta}_t + p u_x - \frac{\varepsilon u_x^2}{v} 
- \frac{Q {\theta}_x}{\theta} \le 0.  
\eea
%
At this point we have  to require certain conditions 
to guarantee  the sign in \eqref{E:5.6}, namely,
%5.7
\bea \label{E:5.7}
\begin{cases}
\eta = -{\Psi}_{\theta}, \tab {\Psi}_Z > 0,\\
p = - {\Psi}_v, \tab {\Psi}_{u_x} = 0,\\
Q = \frac{{\lambda} {\theta}_x}{v}, \tab  
{\Psi}_{{\theta}_x} = 0. \\
\end{cases} 
\eea
%
Here $Q$ is a multiple of ${\theta}_x$ and hence satisfies 
the 
Fourier's law of heat flux.  The conditions in 
\eqref{E:5.7} yield
that the Helmholtz free energy $\Psi$ is independent of 
$u_x$ and ${\theta}_x,$
that is,
$$ %no # 
\Psi = \hat{\Psi}(v, \theta, Z).
$$
Relations  \eqref{E:5.2} and \eqref{E:5.7}   imply that one 
has to look for  
$(\eta, e, \theta)$ such that
%No #
\bea 
\begin{cases}
{\eta}_{\theta} = \frac{e_{\theta}}{\theta},\\
(\theta \eta)_v = p + e_v, \\
(\theta \eta)_z \le e_z,\\
\eta >0.
\end{cases}\nonumber
\eea
Now, by the Dalton law
%No #
\bea
e_{\theta} = c_v(Z)\equiv c_1 Z + c_2 (1-Z). \nonumber
\eea
Therefore,
%No #
\bea
{\eta}_{\theta} = \frac{1}{\theta} \left(c_1 Z + c_2 (1-Z) 
\right), \nonumber
\eea
which implies
%5.8
\bea \label{E:5.8}
\eta = c_v(Z) \log{\theta} + f(v,Z). 
\eea
%
Now, by the Dalton law,
%no #
\bea
\begin{cases}
(\theta \eta)_v = p = (\gamma(Z)-1) \frac{e}{v} 
= \frac{a(Z)}{c_v(Z)} \frac{e}{v},\\
%({\gamma}_1 -1) c_1 Z + ({\gamma}_2 -1) c_2 (1-Z)}{c_1 Z + c_2 (1-Z)}, \\
(\theta \eta)_Z \le c_v{'}(Z) \theta = (c_2 - c_1) \theta. 
\end{cases}\nonumber
\eea
Therefore,
%no #
\bea
(\theta f)_v 
=- \left(c_v(Z)\theta \log \theta \right)_v 
+ a(Z)\frac{\theta}{v},\nonumber
%\left( ({\gamma}_1 -1) c_1 Z + ({\gamma}_2 -1) c_2 (1-Z) \right), \nonumber
\eea
and then
%no #
\bea
f(v,Z) = - c_v(Z)\log{\theta} + a(Z)\log{v}
%(({\gamma}_1 -1)c_1 Z + ({\gamma}_2-1)c_2(1-Z)) \log{v}
+ \frac{\omega({\theta,Z})}{{\theta}}, \nonumber
\eea
or
%no #
\bea
f(v,Z) = a(Z)\log{v} + h(Z), \nonumber
%\left( ({\gamma}_1 -1) c_1 Z + ({\gamma}_2 -1) c_2 (1 -Z) \right) \log{v} 
\eea
with
%no #
\bea
h(Z) = \frac{\omega({\theta},Z)}{\theta} - c_v(Z) \theta 
\log{\theta}. 
\nonumber
\eea
%
Relation $({\theta} \eta)_Z \le c_v{'}(Z) {\theta}$ is 
equivalent to the 
condition
%no #
\bea
h^{'}(Z) & \!\! \le \!\! & (c_2 - c_1) \log {\theta} 
+(({\gamma}_1-1)c_1-({\gamma}_2 -1)c_2) \log{v} + (c_2 - 
c_1)\nonumber\\
& \!\! = \!\! &  (c_2 - c_1) (1 + \log {\theta}) 
 + (({\gamma}_1 -1) c_1 - ({\gamma}_2 -1)c_2) 
\log{v}.\nonumber
\eea
Therefore, the entropy we are seeking is of the form
%no #
\bea
\eta = c_v(Z) \log{\theta} + a(Z)\log{v} + h(Z), \nonumber 
\eea
for appropriate function $h=h(Z).$
\hfill $\blacksquare$

Clearly $\eta = \eta(\cdot)$ is not in general a convex 
function, 
a remark closely related to the fact that the estimates in 
our analysis 
are not in general time-independent.

\begin{remark}\label{R:5.1}
As is well-known, the constitutive equations of a real 
gas are fairly well 
approximated within moderate ranges of $\theta$ and $v$ by 
the model of 
an ideal gas, in which
%5.9
\bea \label{E:5.9}
e=c_v\theta, \tab \sigma=- p(v,\theta)+ \frac{{\varepsilon} 
u_x}{v}, 
\tab Q= \frac{\lambda {{\theta}_x}}{v} 
\eea
%
with suitable constants $c_v, \varepsilon, \lambda.$
However, as our asymptotic results also indicate, under 
very high 
temperatures and densities,  the equations in \eqref{E:5.9}
become inadequate, since in particular specific heat, 
conductivity,
and viscosity vary with temperature and density.
The model introduced in this work is certainly more 
realistic since  it 
takes into consideration the dependence of $c=c_v(Z), 
\gamma=\gamma(Z),
p=p(v,\theta, Z)$. A more realistic model than 
\eqref{E:1.1.1}
would be a linearly viscous gas (or Newtonian fluid)
%no #
\bea
\sigma(v,\theta, u_x, Z) = - p(v, \theta, Z) + 
\frac{\varepsilon(v,\theta)}{v}u_x,
\nonumber
\eea
%
satisfying the Fourier's law of heat flux
%5.10
\bea \label{E:5.10}
Q(v,\theta, {\theta}_x) = \frac{\lambda(v, 
\theta)}{v}{\theta}_x.
\eea 
It would be interesting to investigate such a model and 
compare the difference 
of solution behavior between this model and the previous 
models.
%The investigation of such a model will be the focus of a future research
%project.
\end{remark}
\bigskip

 
\section*{Acknowledgments}
Gui-Qiang Chen's research was supported in part by the 
National Science 
Foundation under Grants DMS-0204225, DMS-0204455, and INT-9987378.
David Hoff's research was supported in part by the National 
Science
Foundation under Grant DMS-9986658.
Trivisa was supported in part by the National Science 
Foundation under
Grants DMS-0072496 and INT-9987378 and under the Alfred P. 
Sloan Foundation 
Research Fellowship.





   

%%%%%%%%%%%%%%%%%%% Bibliography   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
\begin{thebibliography}{99}

\bibitem{AZ} A.~A. Amosov and A.~A. Zlotnick,
A semidiscrete method for solving equations of the 
one-dimensional 
motion of a non-homogeneous viscous heat-conducting gas 
with 
nonsmooth data, 
Izv. Vyssh. Uchebn. Zaved. Mat. 1997, 3--19 (Russian); 
transl. in Russian Math. (Iz. VUZ) 41 (1997), 1--17.

\bibitem{BMR} A. \ Bourlioux, A. J. \ Majda, and V. \ 
Roytburd,
Theoretical and numerical structure for unstable
one-dimensional detonations, 
{\em SIAM J. Appl. Math.} {\bf 51} (1991), 303--343.

\bibitem{CH} T.\ Chang and L.\ Hsiao,
The Riemann Problem and Interaction of Waves in Gas 
Dynamics,
Pitman Monographs and Surveys in Pure and Applied 
Mathematics, 41,
Longman Scientific \& Technical, Harlow; copublished in the 
United 
States with John Wiley \& Sons, Inc.: New York, 1989. 

\bibitem{C} G.-Q. \ Chen, Global solutions to the 
compressible 
Navier-Stokes equations for a reacting mixture,
{\em SIAM  J. Math. Anal.} {\bf 23} (1992), 609--634.

\bibitem{CHT}  G.-Q. \ Chen, D.\ Hoff, and K. \ Trivisa, 
Global
solutions of the compressible Navier-Stokes equations with 
large 
discontinuous initial data, {\em Comm. Partial Diff. Eqs.} 
{\bf 25} (2000), 2233--2257.

\bibitem{CHT2}  G.-Q. \ Chen, D.\ Hoff, and K. \ Trivisa,
On the Navier-Stokes equations for exothermically reacting 
compressible 
fluids, {\em Acta Math. Appl. Sinica},
{\bf 18} (2002), 15-36.
 
\bibitem{CW}  G.-Q. Chen and D. Wagner,
Global entropy solutions to exothermically reacting, 
compressible
Euler equations, {\em J. Diff. Eqs.} 2003  (to appear).


\bibitem{C-M-R} P. \ Collela, A.\ Majda  and  V. \ 
Roytburd, 
Theoretical and numerical structure for reacting shock 
waves, 
{\em SIAM J. Sci. Stat. Comput.} {\bf 7} (1986), 1059--1080.

\bibitem{C-F} R. \ Courant and K. O. \ Friedrichs,
Supersonic Flow and Shock Waves,
{\em Interscience: New York}, 1948.

\bibitem{Da} C.\ M.\ Dafermos,
Global smooth solutions to the initial--boundary value 
problem 
for the equations of one--dimensional nonlinear 
thermoviscoelasticity,
{\em SIAM J. Math. Anal.}  {\bf 13} (1982), 397--408.

\bibitem{G} J. \ Glimm,  
The continuous structure of discontinuities, 
{\em Lecture Notes in Physics} {\bf 344} (1989), 177--186.

\bibitem{GR} E. \ Godlewski and P. \ Raviart, 
Numerical Approximation of Hyperbolic Systems of 
Conservation Laws,
{\em Appl. Math. Sc.} {\bf 118}, Springer-Verlag: New York, 
1996. 

\bibitem{Ho1} D.\ Hoff, Global well--posedness of the 
Cauchy problem 
for nonisentropic gas dynamics with discontinuous initial 
data, {\em J. Diff. Eqs.}  {\bf 95} (1992), 33-74.

\bibitem{Ho2} D.\ Hoff, Discontinuous solutions of the 
Navier--Stokes 
equations for compressible flow, 
{\em Arch. Rational Mech. Anal.} {\bf 114} (1991), 15--46. 

\bibitem{Ho4} D.\ Hoff. Discontinuous solutions of the 
Navier-Stokes
equations for multidimensional heat-conducting fluids,
Arch. Rational Mech. Anal. {\bf 139} (1997), 303--354. 

\bibitem{HS} D. \ Hoff and D. \ Serre. The failure of 
continuous
dependence on the initial data for the Navier Stokes 
equations of
compressible flow, {\em SIAM J. Appl. Math.} {\bf 51} 
(1991), 887-898.

\bibitem{HZ} D. \ Hoff and M. \ Ziane, The global 
attraction and finite determining nodes for the Navier-Stokes 
equations of one
dimensional compressible flow with discontinuous initial 
data, Indiana
Univ. Math. J. {\bf 49(3)} (2000), 843-889.

\bibitem{Ka} Y. Kanel, On a model system of equations of 
one--dimensional  
gas motion, {\em J. Diff. Eq.} {\bf 4} (1968), 374--380.

\bibitem{K} A. V.\ Kazhikhov. On  the theory of 
initial--boundary--value 
problems for the equations of one--dimensional 
nonstationary motion 
of a viscous heat-conductive gas, 
{\em Din. Sploshnoi Sredy} {\bf 50} (1981), 37--62 
(Russian).

\bibitem{KS} A. V. \ Kazhikhov and V. V. \ Shelukhin,
Unique global solution with respect to time of 
initial-boundary-value 
problems for one--dimensional equations of a viscous gas, 
{\em J. Appl. Math. Mech.} {\bf 41} (1977), 273--282.

\bibitem{Maj} A. Majda. High Mach number combustion, 
{\em Lecture Notes in Appl. Math.} {\bf 24} (1986), 
109--184.

\bibitem{L} G. S. S. Ludford. Low Mach number combustion, 
{\em Lecture
Notes in Appl. Math.} {\bf 24} (1986), 3--74.

\bibitem{MY} A. Matsumura  and S. Yanagi,
Uniform boundedness of the solutions for a one-dimensional 
isentropic 
model system of compressible viscous gas, 
{\em Comm.  Math. Phys.} {\bf 175} (1996), 259--274. 

\bibitem{Wi} F.~A. Williams, {\em Combustion Theory},
     Addison-Wesley, Reading, MA, 1965.


\end{thebibliography}

\bigskip
\bigskip
\medskip
\address{\relax
\noindent
Gui-Qiang Chen\\
Department of Mathematics, Northwestern University\\ 
Evanston, IL 60208-2730, USA\\
\email{gqchen@math.northwestern.edu}\\
\\
David Hoff\\
Mathematics Department,  Indiana University  \\
Bloomington, IN 47405-4301, USA\\
\email{hoff@indiana.edu }\\
\\
Konstantina Trivisa\\
Department of Mathematics, University of Maryland\\
College Park, MD 20742-4015, USA\\
\email{trivisa@math.umd.edu}\\
}


\end{document}

