# Approximation of the Equations of the Humid Atmosphere with Saturation

Roger Temam, Xiaoyan Wang

► To cite this version:

Roger Temam, Xiaoyan Wang. Approximation of the Equations of the Humid Atmosphere with Saturation. 27th IFIP Conference on System Modeling and Optimization (CSMO), Jun 2015, Sophia Antipolis, France. pp.21-42, 10.1007/978-3-319-55795-3_2 . hal-01626922

HAL Id: hal-01626922

https://inria.hal.science/hal-01626922

Submitted on 31 Oct 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Approximation of the equations of the humid atmosphere with saturation

Roger Temam and Xiaoyan Wang

Indiana University, Institute for Scientific Computing and Applied Mathematics, Rawles Hall, Bloomington, IN 47405, USA

temam@indiana.edu (R. Temam), wang264@indiana.edu (X. Wang)

**Abstract.** We investigate the numerical approximation of solutions to some variational inequalities modeling the humid atmosphere when the saturation of water vapor in the air is taken into account. Here we describe part of our work [31] and extend our former results to the case where the saturation $q_s$ evolves with time.

**Keywords:** Atmosphere equation, Variational inequality, Penalization, Regularization, Uniform estimates, Fractional step method

## 1 Introduction

The rigorous mathematical theory of the equations of humid atmosphere has been initiated in [21, 22] and has attracted the attention of a large number of researchers, see e.g., [1, 3–8, 13–16, 24, 30] and the references therein. These cited research works solved a large class of practical problems by investigating the system of partial differential equations based on different accuracies of the mathematical modelings [17, 18, 23, 26]. However, in the modelings in [17, 18, 23, 26], the saturation of vapor is not taken into account. As shown in [29, 32], the resulted systems of partial differential equations are not physically correct in the extreme cases where the atmosphere is totally dry, $q = 0$, or when the atmosphere is totally humid, $q = 1$. To remedy this drawback, we have proposed in [32] a new formulation of the problem in the context of the variational inequalities [2, 9, 19, 20, 25]. This new variational inequality formulation also involves discontinuities due to phase changes. In this work, we describe the numerical approximation of the solutions to the variational inequalities derived from the humidity equations when the saturation of water vapor in the air is taken into account. As explained above, a striking feature of our work here is that the problems we study contain discontinuities and involve inequalities which come from the changes of phases and the extreme cases for the vapor concentration respectively (see e.g., [10–12]). In addition, we manage to extend our recent study [31] to the case that the saturation concentration $q_s$ evolves according to the thermodynamical laws.

In [31], we proposed an implicit Euler scheme to approach the solutions to the system which involves a variational inequality. However, we can not simply proceed directly as usual due to the difficulties induced from the discontinuities
and physical requirement of vapor concentration (i.e. $0 \le q \le 1$). To overcome the difficulty caused by the discontinuities in our current modeling, we use a regularization method. This regularization enables us to first study a system of partial differential equations and then discuss the approximation of the solutions to the original variational inequality. The constraint requirement $q \in \mathcal{K}$ for almost every $t \in [0, t_1]$ for the vapor concentration $q = q(t, \mathbf{x})$ brings us great technical challenges. Here $t_1 > 0$ is an arbitrary but fixed constant. See Section 2.1 for more details about this physical range requirement. The source of challenges in our study is that this range requirement can not be preserved in the discretization procedure in the implicit Euler scheme. To deal with these challenges, we devised a penalization technique in the regularized Euler scheme. Together with delicate energy estimates, the penalization technique can elegantly help us achieve the physical requirement on $q$. We point out here that the forms and the signs of the penalization terms encode very elegant structural propositions and are crucial for us to obtain the desired energy estimates. Finally, when we extend the study to the case where the saturation vapor concentration $q_s$ evolves according to the thermodynamical laws, we emphasize that the discretization of the $q_s$-equation is of different nature from the discretization of the temperature equation on $T$ and vapor concentration equations on $q$. See Remark 2 for more detailed comments.

The rest of the article is organized as follows. In Section 2, we give the formulation of the problem. In Section 3, we introduce the Euler scheme and derive various uniform estimates for the functions associated with the penalized and regularized scheme. In Section 4, we investigate the convergence of the Euler scheme. We devote Section 5 to the study of the implicit Euler scheme in the case where $q_s$ depends on time.

# 2 The problem

## 2.1 Formulation of the system

Let $\mathcal{M} = \mathcal{M}' \times (p_0, p_1)$ where $\mathcal{M}' \subset \mathbb{R}^2$ is a bounded domain with smooth boundary and $p_0, p_1$ are two real numbers with $0 < p_0 < p_1$. We will use $\mathbf{x} = (x, y, p)$ to denote a typical point in $\mathcal{M}$, and use $n$ to denote the outward normal vector to $\partial \mathcal{M}$, the boundary of $\mathcal{M}$. Let $\mathcal{K}$ be the non-empty closed convex set in $H^1(\mathcal{M})$ defined as $\mathcal{K} = \{q \in H^1(\mathcal{M}); 0 \le q \le 1, \text{ a.e.}\}$. Given a fixed $t_1 > 0$, we consider the following problem:

To find $T: (0, t_1) \to H^1(\mathcal{M})$, $q: (0, t_1) \to \mathcal{K}$ and $h_q \in \mathcal{H}(q - q_s)$ such that for $q^b \in \mathcal{K}$, there hold

$$ \partial_t T + A_T T + \mathbf{v} \cdot \nabla T + \omega \partial_p T - \frac{R\omega}{c_p p} T = \frac{1}{p} \omega^{-} h_q \varphi(T), \quad (1) $$

$$ \langle \partial_t q, q^b - q \rangle + (A_q q + \mathbf{v} \cdot \nabla q + \omega \partial_p q, q^b - q) \geq \left( -\frac{1}{p} \omega^{-} h_q F(T), q^b - q \right), \quad (2) $$
with initial and boundary conditions to be specified. Here $\mathcal{H}$ is the multivalued Heaviside function such that $\mathcal{H}(\tau) = 0$ for $\tau < 0$, $\mathcal{H}(0) = [0, 1]$, $\mathcal{H}(\tau) = 1$ for $\tau > 0$.

For the sake of simplicity, the velocity field of the fluid $\mathbf{u} := (\mathbf{v}(\mathbf{x}, t), \omega(\mathbf{x}, t)) \in \mathbb{R}^3$ is considered as a given data in this article. Throughout the presentation, we assume that the time-dependent velocity field $\mathbf{u}$ satisfies $\mathbf{u} \in L^r(0, t_1; V) \cap L^\infty(0, t_1; H)$ for some given $r \in (4, +\infty]$; $\nabla = (\partial_x, \partial_y)$ and $\Delta = \partial_x^2 + \partial_y^2$ are the horizontal gradient and horizontal Laplace operators respectively. In this way, the operators $\mathcal{A}_T$ and $\mathcal{A}_q$ are defined as

$$ \mathcal{A}_T = -\mu_1 \Delta - \nu_1 \partial_p ((\frac{gp}{RT})^2 \partial_p), \quad \mathcal{A}_q = -\mu_2 \Delta - \nu_2 \partial_p ((\frac{gp}{RT})^2 \partial_p), \quad (3) $$

where $\mu_i, \nu_i, g, R, c_p$ are positive constants and $\bar{T} = \bar{T}(p)$ is the average temperature over the isobar with pressure $p$. We assume that $\bar{T}$ satisfies:

$$ \bar{T}_* \le \bar{T}(p) \le \bar{T}^*, \quad |\partial_p \bar{T}(p)| \le M, \quad \text{for some positive constants } \bar{T}_*, \bar{T}^*, M \text{ and } p \in [p_0, p_1]. \quad (4) $$

Concerning the right hand sides of equations (1)-(2), the functions $F$ and $\varphi$ both from $\mathbb{R}^1$ to $\mathbb{R}^1$ are defined as

$$ F(\zeta) = q_s \zeta \frac{RL(\zeta) - c_p R_v \zeta}{c_p R_v \zeta^2 + q_s L(\zeta)^2}, \quad \text{with } L(\zeta) = c_1 - c_2 \zeta; \quad \varphi(\zeta) = \frac{1}{c_p} L(\zeta) F(\zeta). \quad (5) $$

Above, $c_1, c_2, R_v, R_q$ are all strictly positive constants. It is easy to see that $F$ is bounded and that both functions $F$ and $\varphi$ are globally Lipschitz; $\omega^+ := \max\{\omega, 0\}$ refers to the positive part of $\omega$.

We partition the boundary of $\mathcal{M}$ as $\partial\mathcal{M} = \Gamma_i \cup \Gamma_u \cup \Gamma_l$ with $\Gamma_i, \Gamma_u$ and $\Gamma_l$ defined by

$$ \Gamma_i = \{\mathbf{x} \in \overline{\mathcal{M}}; p = p_1\}, \quad \Gamma_u = \{\mathbf{x} \in \overline{\mathcal{M}}; p = p_0\}, \quad \Gamma_l = \{\mathbf{x} \in \overline{\mathcal{M}}; p_0 \le p \le p_1, (x, y) \in \partial\mathcal{M}'\}. $$

We supplement the system (1)-(2) with the boundary conditions

$$ \left\{ \begin{array}{l} \partial_p T = \alpha(T_* - T), \partial_p q = \beta(q_* - q) \text{ on } \Gamma_i, \\ \partial_p T = 0, \partial_p q = 0 \text{ on } \Gamma_u, \\ \partial_n T = 0, \partial_n q = 0 \text{ on } \Gamma_l, \end{array} \right. \quad (6) $$

and initial conditions

$$ T(\mathbf{x}, 0) = T_0(\mathbf{x}), \quad q(\mathbf{x}, 0) = q_0(\mathbf{x}). \quad (7) $$

If we allow $q_s$ to evolve, the dependence of the nonlinear functions $F$ and $\varphi$ on $q_s$ should be made explicit, i.e., $F = F(T, q_s)$, $\varphi = \varphi(T, q_s)$. In this case, we augment (1)-(2) with the following governing equation for the evolution of $q_s$ (see [17], [18], [26]):

$$ \frac{dq_s}{dt} = -\frac{\delta\omega^-}{p} h_q F(T, q_s). \quad (8) $$
We will also impose a further initial condition

$$q_s(\mathbf{x}, 0) = q_{s,0}(\mathbf{x}). \qquad (9)$$

We shall always assume $q_0 \in L^2(\mathcal{M})$, $0 \le q_0 \le 1$ for a.e. $\mathbf{x} \in \mathcal{M}$, and $q_{s,0} \in L^2(\mathcal{M}) \cap L^\infty(\mathcal{M})$, $0 < q_{s,0} < 1$ and $0 \le q_* \le 1$ for a.e. $\mathbf{x} \in \mathcal{M}$, and assume the boundary datum $T_*$ and $q_*$ to satisfy $T_*, q_* \in L^2(0, t_1; L^2(\Gamma_i))$. For the convenience of fixing the ideas and of presentation, we shall first assume $q_s$ is stationary during our study. In the last section, we will explain the case where $q_s$ evolves according to the governing equation (8).

## 2.2 Functional analytic framework

We denote as usual $H = L^2(\mathcal{M})$, $V = H^1(\mathcal{M})$. We use $(\cdot, \cdot)_{L^2}$ (regarded the same as $(\cdot, \cdot)_H$) and $|\cdot|_{L^2}$ to denote the usual scalar product and induced norm in $H$. In the space $V$, we will use $((\cdot, \cdot))$ and $\|\cdot\|$ to denote the scalar product adapted to the problem under investigation

$$((\varphi, \phi)) := (\nabla\varphi, \nabla\phi) + (\partial_p\varphi, \partial_p\phi) + \int_{\Gamma_i} \varphi\phi \, d\Gamma_i,$$

and the induced norm. The symbol $\langle \cdot, \cdot \rangle$ will denote the duality pair between a Banach space $E$ and its dual space $E^*$. We use the following standard function spaces for the vector field $\mathbf{u}$:

$$\mathbf{H} = \{\mathbf{u} \in H \times H \times H \mid \mathrm{div} \, \mathbf{u} = 0 \text{ and } \mathbf{u} \cdot \mathbf{n} = 0 \text{ on } \partial \mathcal{M}\},$$

$$\mathbf{V} = \{\mathbf{u} \in V \times V \times V \mid \mathrm{div} \, \mathbf{u} = 0 \text{ and } \mathbf{u} \cdot \mathbf{n} = 0 \text{ on } \partial \mathcal{M}\}.$$

For $T, T^b, q, q^b \in V$, we have the following specific forms for the duality pairs through integration by parts and in view of the Neumann boundary conditions:

$$\langle A_T T, T^b \rangle = \mu_1(\nabla T, \nabla T^b)_H + \nu_1 \int_{\mathcal{M}} \left(\frac{gp}{RT}\right)^2 \partial_p T \partial_p T^b \, d\mathcal{M} + \nu_1 \int_{\Gamma_i} \left(\frac{gp_1}{RT}\right)^2 \alpha(T-T_*) T^b \, d\Gamma_i, \qquad (10)$$

$$\langle A_q q, q^b \rangle = \mu_2(\nabla q, \nabla q^b)_H + \nu_2 \int_{\mathcal{M}} \left(\frac{gp}{RT}\right)^2 \partial_p q \partial_p q^b \, d\mathcal{M} + \nu_2 \int_{\Gamma_i} \left(\frac{gp_1}{RT}\right)^2 \beta(q-q_*) q^b \, d\Gamma_i. \qquad (11)$$

Consequently, we define the following bilinear forms

$$a_T(T, T^b) = \mu_1(\nabla T, \nabla T^b)_H + \nu_1 \int_{\mathcal{M}} \left(\frac{gp}{RT}\right)^2 \partial_p T \partial_p T^b \, d\mathcal{M} + \nu_1 \alpha \int_{\Gamma_i} \left(\frac{gp_1}{RT}\right)^2 TT^b \, d\Gamma_i, \qquad (12)$$

$$a_q(q, q^b) = \mu_2(\nabla q, \nabla q^b)_H + \nu_2 \int_{\mathcal{M}} \left(\frac{gp}{RT}\right)^2 \partial_p q \partial_p q^b \, d\mathcal{M} + \nu_2 \beta \int_{\Gamma_i} \left(\frac{gp_1}{RT}\right)^2 q q^b \, d\Gamma_i. \qquad (13)$$

Meanwhile, we set $U := (T, q)$, $U^b := (T^b, q^b)$ and introduce the bilinear form

$$a(U, U^b) := a_T(T, T^b) + a_q(q, q^b). \qquad (14)$$
As for the Navier-Stokes equation, we define

$$b(\mathbf{u}, U, U^b) := \int_{\mathcal{M}} (\mathbf{u} \cdot \nabla_{x,y,p} U) \cdot U^b d\mathcal{M} = b_T(\mathbf{u}, T, T^b) + b_q(\mathbf{u}, q, q^b), \quad (15)$$

where $b_T(\mathbf{u}, T, T^b)$ and $b_q(\mathbf{u}, q, q^b)$ are given by

$$b_T(\mathbf{u}, T, T^b) = \int_{\mathcal{M}} (\mathbf{v} \cdot \nabla T + \omega \partial_p T) T^b d\mathcal{M}, \quad b_q(\mathbf{u}, q, q^b) = \int_{\mathcal{M}} (\mathbf{v} \cdot \nabla q + \omega \partial_p q) q^b d\mathcal{M}. \quad (16)$$

In view of the last term in the left hand side of the (1), we introduce the following bilinear form

$$d(\omega, T, T^b) = \int_{\mathcal{M}} \frac{R\omega TT^b}{c_p p} d\mathcal{M}. \quad (17)$$

Similarly, in view of the last terms in (10) and (11), we define the linear functional

$$l(U^b) := l_T(T^b) + l_q(q^b) = \nu_1 \alpha \int_{\Gamma_i} \left(\frac{gp_1}{RT}\right)^2 T_* T^b d\Gamma_i + \nu_2 \beta \int_{\Gamma_i} \left(\frac{gp_1}{RT}\right)^2 q_* q^b d\Gamma_i. \quad (18)$$

Next, we consider the mapping relations related to the operators $A_T$, $A_q$ and the above defined functionals.

It is well-known that the linear operators $A_T$, $A_q: V \to V^*$ defined through the relations

$$\langle A_T u, v \rangle := a_T(u, v), \quad \langle A_q u, v \rangle := a_q(u, v), \forall u, v \in V, \quad (19)$$

are both bounded linear operators.

Similarly, the operators $B(\mathbf{u}, U) = (B_T(\mathbf{u}, U), B_q(\mathbf{u}, q)) : \mathbf{V} \times V^2 \to (V^*)^2$ and $D(\mathbf{u}, u) : \mathbf{H} \times V \to V'$ defined by

$$\langle B(\mathbf{u}, U), U^b \rangle := (b_T(\mathbf{u}, T, T^b), b_q(\mathbf{u}, q, q^b)), \quad \forall \mathbf{u} \in \mathbf{V}, U, U^b \in V^2, \quad (20)$$

and

$$\langle D(\mathbf{u}, u), v \rangle := d(\omega, u, v), \forall \mathbf{u} \in \mathbf{H}, u, v \in V, \quad (21)$$

are also bounded.

Due to the divergence free condition of $\mathbf{u}$, we easily see that for any $T, q \in V$,

$$b_T(\mathbf{u}, T, T) = 0, \quad b_q(\mathbf{u}, q, q) = 0. \quad (22)$$

Concerning the boundedness of the above functionals, we have the following lemma.

**Lemma 1 (Boundedness of the functionals).** Assume $U, U^b \in V^2$ and $\mathbf{u} \in \mathbf{V}$. There exist universal positive constants $\lambda$ and $K_i$, $1 \le i \le 6$ such that

$$|a_T(T, T^b)| \le K_1 \|T\| \|T^b\|, \quad a_T(T, T) \ge \lambda \|T\|^2; \quad (23)$$
$$|a_q(q, q^b)| \le K_2 \|q\| \|q^b\|, \quad a_q(q, q) \ge \lambda \|q\|^2; \qquad (24)$$

$$|b(\mathbf{u}, U, U^b)| \le K_3 \|\mathbf{u}\|_{V} |U|^{\frac{1}{2}} \|U\|^{\frac{1}{2}} \|U^b\|; \qquad (25)$$

$$|d(\omega, T, T^b)| \le K_4 |\omega|_{L^2} |T|^{\frac{1}{4}} \|T\|^{\frac{3}{4}} |T^b|^{\frac{1}{4}} \|T^b\|^{\frac{3}{4}}; \qquad (26)$$

$$|l_T(T^b)| \le K_5 \|T^b\|, \quad |l_q(q^b)| \le K_6 \|q^b\|. \qquad (27)$$

**Definition 1.** Let $(T_0, q_0) \in H \times H$ be such that $0 \le q_0 \le 1$ a.e. in $M$ and let $t_1 > 0$ be fixed. A vector $U = (T, q) \in L^2(0, t_1; V \times V) \cap C([0, t_1]; H \times H)$ with $(\partial_t T, \partial_t q) \in L^2(0, t_1; V^* \times V^*)$ is a solution to the initial and boundary value problem described by (1), (2), (6) and (7), if for almost every $t \in [0, t_1]$ and for every $(T^b, q^b) \in V \times K$, we have

$$\langle \partial_t T, T^b \rangle + a_T(T, T^b) + b_T(\mathbf{u}, T, T^b) - d(\omega, T, T^b) - l_T(T^b) = \left(\frac{1}{p}\omega^{-}(t)h_q\varphi(T), T^b\right), \qquad (28)$$

$$\langle \partial_t q, q^b - q \rangle + a_q(q, q^b - q) + b_q(\mathbf{u}, q, q^b - q) - l_q(q^b - q) \geq \left(-\frac{1}{p}\omega^{-}(t)h_q F(T), q^b - q\right), \qquad (29)$$

for some $h_q \in \mathcal{H}(q - q_s)$ and

$$U_0 = (T_0, q_0). \qquad (30)$$

# 3 Time discretization-The Euler scheme

## 3.1 Time-discretization

We assume that the velocity field $\mathbf{u}$ is given, time-dependent and satisfies $\mathbf{u} \in L^r(0, t_1; V) \cap L^\infty(0, t_1; H)$ for some given $r \in (4, +\infty]$.

Let $N$ be an integer which will eventually go to $+\infty$ and set $\Delta t := k = t_1/N$. We will define recursively a family of elements of $V \times K$, say $(T^0, q^0)$, $(T^1, q^1)$, ..., $(T^N, q^N)$, where $(T^m, q^m)$ will be in some sense an approximation of the functions $(T, q)$ we are looking for, on the interval $[(m-1)k, mk)$.

First, we define $\mathbf{u}^m = \frac{1}{k} \int_{(m-1)k}^{mk} \mathbf{u}(t) dt$, $m = 1, 2, \dots, N$. Our discretization is as follows:

We begin with $(T^0, q^0) := (T_0, q_0)$, i.e., the given initial datum. When $(T^0, q^0)$, $(T^1, q^1)$, ..., $(T^{m-1}, q^{m-1})$ are known, $T^m \in V$ and $q^m \in K$ are determined by:

$$\langle \frac{T^m - T^{m-1}}{k}, T^b \rangle + a_T(T^m, T^b) + b_T(\mathbf{u}^m, T^m, T^b) - d(\omega^m, T^{m-1}, T^b) - l_T(T^b) = \left(\frac{1}{p}[\omega^m]^{-} h_{Q^m} \varphi(T^{m-1}), T^b\right), \qquad (31)$$

$$\langle \frac{q^m - q^{m-1}}{k}, q^b - q^m \rangle + a_q(q^m, q^b - q^m) + b_q(\mathbf{u}^m, q^m, q^b - q^m) - l_q(q^b - q^m) \\ \geq \left(-\frac{1}{p}[\omega^m]^{-} h_{Q^m} F(T^{m-1}), q^b - q^m\right), \qquad (32)$$
for some $h_{Q^m} \in \mathcal{H}(Q^m - q_s)$ where $Q^m$ is either $q^{m-1}$ or $q^m$.

In the above construction of the discretization scheme (31)-(32), one shall pay special attention to the indices in the terms $d$ and $\varphi$. Notice that we have $d(\omega^m, T^{m-1}, T^b)$ and $\varphi(T^{m-1})$. Obviously, this choice of indices will have influence on our search for $T^m$ and $q^m$ recursively. More importantly for us, it is crucial for us to obtain energy estimates later: the required estimates would not be true if we changed the indices $m-1$ to be $m$ in $d$ and $\varphi$. However, the choice of the indices in $F$ and $h_{Q^m}$ is not so sensitive.

**Remark 1.** In the above discretization, we have to deal with variational inequalities due to the $q$-equation (32). Meanwhile, we shall keep in mind that the physical constraint on the function $q$ in our problem, $q^m \in \mathcal{K}$ is not preserved during the discretization (31)-(32). Finally, the problem we meet is nonlinear. The above three aspects form the main sources of difficulties for our study.

## 3.2 Regularization and penalization

In view of Remark 1, we proceed our investigation by way of regularization and penalization. Let $\varepsilon = (\varepsilon_1, \varepsilon_2)$ and $\varepsilon_i > 0$ be small for $i = 1, 2$. For $\varepsilon_2 > 0$, we define as follows the regularization $H_{\varepsilon_2}$ of $\mathcal{H}(\cdot) : \mathbb{R} \to [0, 1]$: equal to 0 for $\eta \ge 0$, to 1 for $\eta \ge \varepsilon_2$, and linear continuous between 0 and $\varepsilon_2$. and consider the associated regularized and penalized problem:

To find $T_{\varepsilon}^m, q_{\varepsilon}^m \in V$ such that

$$ \begin{aligned} & \left\langle \frac{T_{\varepsilon}^{m} - T_{\varepsilon}^{m-1}}{k}, T^{b} \right\rangle + a_{T}(T_{\varepsilon}^{m}, T^{b}) + b_{T}(\mathbf{u}^{m}, T_{\varepsilon}^{m}, T^{b}) - d(\omega^{m}, T_{\varepsilon}^{m-1}, T^{b}) - l_{T}(T^{b}) \\ & = \left(\frac{1}{p}[\omega^{m}]^{-} H_{\varepsilon_{2}}(Q_{\varepsilon}^{m} - q_{s})\varphi(T_{\varepsilon}^{m-1}), T^{b}\right), \end{aligned} \quad (33) $$

$$ \begin{aligned} & \left\langle \frac{q_{\varepsilon}^{m} - q_{\varepsilon}^{m-1}}{k}, q^{b} \right\rangle + a_{q}(q_{\varepsilon}^{m}, q^{b}) + b_{q}(\mathbf{u}^{m}, q_{\varepsilon}^{m}, q^{b}) - l_{q}(q^{b}) \\ & = \left( \frac{1}{\varepsilon_{1}}[q_{\varepsilon}^{m}]^{-}, q^{b} \right) - \left( \frac{1}{\varepsilon_{1}}[q_{\varepsilon}^{m} - 1]^{+}, q^{b} \right) - \left( \frac{1}{p}[\omega^{m}]^{-} H_{\varepsilon_{2}}(Q_{\varepsilon}^{m} - q_{s}) F(T_{\varepsilon}^{m-1}), q^{b} \right), \end{aligned} \quad (34) $$

for all $T^b, q^b \in V$.

Notice that we have two choices for $Q_{\varepsilon}^m$ either $Q_{\varepsilon}^m = q_{\varepsilon}^{m-1}$ or $Q_{\varepsilon}^m = q_{\varepsilon}^m$. The introduction of penalization in the scheme (33)-(34) is designed to remedy the difficulty brought by the physical range requirement for the humidity $q$. The regularization process will overcome the difficulty caused by the variational inequality and the requirement on $h_q$. Here one may suspect that the two penalization terms in (34) may be potentially dangerous due to the blowing up factor $\frac{1}{\varepsilon_1}$. However, we point out that we could still obtain elegant estimates which do not depend on $\varepsilon_1$ (and $\varepsilon_2, k$) though we have a blowing up factor $\frac{1}{\varepsilon_1}$ when we pass to the limit $\varepsilon \to (0+, 0+)$. These estimates will yield that the limit functions $q^m$ of $q_{\varepsilon}^m$ satisfy the range requirement, i.e., $0 \le q^m \le 1$ for $m = 1, 2, \dots, N$ and a.e. $\mathbf{x} \in \mathcal{M}$.
### 3.3 Validity of iteration

The scheme (33)-(34) yields elliptic system on $(T_{\varepsilon}^m, q_{\varepsilon}^m)$ when $(T_{\varepsilon}^{m-1}, q_{\varepsilon}^{m-1})$ is known. To carry out our program, the step of finding $(T_{\varepsilon}^m, q_{\varepsilon}^m)$ given $(T_{\varepsilon}^{m-1}, q_{\varepsilon}^{m-1})$ is indispensable. To realize this iteration step, we need some surjective or existence theorems. Typically, we can use the Minty-Browder surjective theorem, Lax-Milgram theorem or Galerkin method. Here we can realize the iteration step by different methods depending on the choices of $Q_{\varepsilon}^m$ in our scheme (33)-(34). When $Q_{\varepsilon}^m = q_{\varepsilon}^{m-1}$, the factor $H_{\varepsilon_2}(\cdot)$ is known when we proceed to obtain $T_{\varepsilon}^m$ and $q_{\varepsilon}^m$ once $T_{\varepsilon}^{m-1}$ and $q_{\varepsilon}^{m-1}$ are known. We can apply the Minty-Browder surjective theorem (see e.g., [20]) or Galerkin method (see e.g. [28]) to derive the existence of $(T_{\varepsilon}^m, q_{\varepsilon}^m)$. On the other hand, if $Q_{\varepsilon}^m = q_{\varepsilon}^m$, we only can proceed by Galerkin method (see e.g., [27, 28]) to derive the existence of $(T_{\varepsilon}^m, q_{\varepsilon}^m)$, since the factor $H_{\varepsilon_2}(\cdot)$ is not known when we proceed to obtain $T_{\varepsilon}^m$ and $q_{\varepsilon}^m$ even though $T_{\varepsilon}^{m-1}$ and $q_{\varepsilon}^{m-1}$ are known.

### 3.4 A priori estimate for $(T_{\varepsilon}^m, q_{\varepsilon}^m)$

The a priori estimates on $(T_{\varepsilon}^m, q_{\varepsilon}^m)$ independent of $k$ and $\varepsilon$ for the regularized and penalized problem (33)-(34) will be crucial for the processes of passing to the limits $\varepsilon \to (0+, 0+)$ and $k \to 0+$.

**Lemma 2.** We have the estimates

$$ 
\begin{aligned}
& |U_{\varepsilon}^j|_{L^2}^2 \le C(\mathbf{u}, U_0, t_1), \quad \forall 1 \le j \le N, \\
& \sum_{m=1}^N |U_{\varepsilon}^m - U_{\varepsilon}^{m-1}|_{L^2}^2 \le C(\mathbf{u}, U_0, t_1), \\
& k \sum_{m=1}^N \|U_{\varepsilon}^m\|^2 \le C(\mathbf{u}, U_0, t_1),
\end{aligned} 
\quad (35) 
$$

where $C(\mathbf{u}, U_0, t_1)$ is a finite constant depending on the given datum $\mathbf{u}, U_0$ and $t_1$, but independent of $\varepsilon$ and $k$.

In Lemma 2, the process to obtain the estimates on the $T^m$'s is more involved than that for the $q^m$'s. The reason lies in the fact that the function $F$ is bounded while $\varphi$ is not. Due to this reason, we need the following version of the so-called discrete Gronwall lemma [33]:

**Lemma 3 (Discrete Gronwall Lemma).** Let $\theta$ be any positive constant and $N_0 > 1$ be an integer. Suppose the three nonnegative number sequences $(X_m), (Y_m)$ and $(Z_m)$ for $m = 0, 1, 2, \dots, N_0$ satisfy the following relation

$$ X_m \le X_{m-1}(1 + \theta Y_m) + \theta Z_m. \quad (36) $$

Then for $m = 1, 2, 3, \dots, N_0$, the following estimates hold

$$ X_m \le X_0 \exp\left(\sum_{i=0}^{m-1} \theta Y_{i+1}\right) + \sum_{i=1}^{m-1} \theta Z_i \exp\left(\sum_{j=i}^{m-1} \theta Y_{j+1}\right) + \theta Z_n. \quad (37) $$
The iteration relation (36) in the Discrete Gronwall Lemma explains our choice of index in $\varphi(T^{m-1})$ and $d(\omega^m, T^{m-1}, T^b)$ in the initial discretization.

We have the following a priori bound for the norm $k \sum_{m=1}^{N} \| \frac{U_{\varepsilon}^{m} - U_{\varepsilon}^{m-1}}{k} \|_{V^{*}}^{2}$, which will be used in our compactness argument.

**Lemma 4.** For any $\varepsilon_1 > 0$ and any $\varepsilon_2 > 0$, the inequality

$$k \sum_{m=1}^{N} \left\| \frac{U_{\varepsilon}^{m} - U_{\varepsilon}^{m-1}}{k} \right\|_{V^{*}}^{2} \leq C(\mathbf{u}, U_{0}, t_{1}) < +\infty, \quad (38)$$

holds for some constant $C(\mathbf{u}, U_0, t_1)$ depending on $U_0, \mathbf{u}, t_1$, but not on $\varepsilon$ and $k$.

The main point of Lemma 4 is that the bound is independent of $\varepsilon = (\varepsilon_1, \varepsilon_2)$ and any $k$. As $\varepsilon_2$ comes into play through the regularization function $H_{\varepsilon_2}$ and $H_{\varepsilon_2}$ is bounded say by 1, it is easy to obtain the bound independent of $\varepsilon_2$. Therefore, the main issue here is to control the penalization terms which contain a blowing up factor $\frac{1}{\varepsilon_1}$ in the limit process $\varepsilon \to (0+, 0+)$. We have the following bounds for the penalization terms.

**Lemma 5.** The following bounds hold:

$$k \sum_{m=1}^{N} \left| \frac{[q_{\varepsilon}^{m}]^{-}}{\varepsilon_{1}} \right|_{L^{2}}^{2} \leq C |\omega|_{L^{2}(0,t_{1};H)}^{2}, \quad k \sum_{m=1}^{N} \left| \frac{[q_{\varepsilon}^{m}-1]^{+}}{\varepsilon_{1}} \right|_{L^{2}}^{2} \leq C |\omega|_{L^{2}(0,t_{1};H)}^{2}. \quad (39)$$

The proof of the estimates on the two penalization terms is subtle. Let us briefly illustrate this point. To prove the two estimates in Lemma 5, we choose the test function $q^b = [q^m]^{-}$ and $[q^m - 1]^+$ in (34) respectively. Then the terms $\langle \frac{q^m - q^{m-1}}{k}, [q^m]^{-} \rangle$ and $\langle \frac{q^m - q^{m-1}}{k}, [q^m - 1]^+ \rangle$ will appear. However, neither term has a favorable sign. Actually, we think that the two kinds of terms may not have the same sign for different $k = 1, 2, \dots, N$. Here we need more quantitative estimates. Interestingly, though for each fixed $k$, the above two kinds of terms may not have a definite sign, we have the following definite signs for their sums

$$\sum_{m=1}^{N} \langle \frac{q^m - q^{m-1}}{k}, [q^m]^{-} \rangle \leq 0, \quad -\sum_{m=1}^{N} \langle \frac{q^m - q^{m-1}}{k}, [q^m - 1]^+ \rangle \leq 0, \quad (40)$$

by writing $q^m = [q^m]^+ - [q^m]^{-}$ and the same for $q^{m-1}$ directly. Due to the form of the estimates in Lemma 5, the two relations in (40) are sufficient to derive the proof of Lemma 5.

## 3.5 Passage to the limit $\varepsilon \to (0+, 0+)$

Assume the time step $k > 0$ is fixed. Our goal, in this part, is to pass to the limit $\varepsilon \to (0+, 0+)$ in the scheme (33)-(34). The limit functions $(T^m, q^m)$ of $(T^m_\varepsilon, q^m_\varepsilon)$ will be solutions to the time discretized scheme (31)-(32). These solutions will serve as building blocks for us to construct approximate solutions to our original problem.
After extracting a finite number of subsequences, $\epsilon \to 0$, we infer from Lemma 2 that, for $m = 1, 2, \dots, N$ there exist functions $U^m \in V$ such that, as $\epsilon \to 0+$

$$ U_{\epsilon}^{m} \to U^{m} \text{ weakly in } V. \qquad (41) $$

We still use $\epsilon$ as the index for the subsequence.

Since the inclusion $V \subset H$ is compact and $U_{\epsilon}^m$ is weakly convergent in $V$, it is strongly convergent in $H$, i.e., we also have

$$ U_{\epsilon}^{m} \to U^{m} \text{ strongly in } H. \qquad (42) $$

By an additional extraction of subsequences we see that:

$$ U_{\epsilon}^{m}(x) \to U^{m}(x) \text{ a.e., } m = 1, 2, \dots, N. \qquad (43) $$

Meanwhile, we have $H_{\epsilon_2}(Q_{\epsilon}^m - q_s) \to h_{Q^m}$ weakly in $L^\infty(\mathcal{M})$.

Concerning the limit functions $q^m$, the second component of $U^m$ for $m = 1, 2, \dots, N$, we know from Lemma 5 that

$$ k \sum_{m=1}^{N} (|[q_{\epsilon}^{m}]^{-}|_{L^2}^{2} + |[q_{\epsilon}^{m} - 1]^{+}|_{L^2}^{2}) \leq C \epsilon_1^2 |\omega|_{L^2(0,t_1;H)}^2. \qquad (44) $$

As the real functions $g_{\pm}(\theta) = \theta^{\pm}$ are both Lipschitz functions with Lipschitz constant 1 on $\mathbb{R}$, i.e., $|g_{\pm}(\theta_1) - g_{\pm}(\theta_2)| \leq |\theta_1 - \theta_2|$, we have

$$ |[q_{\epsilon}^{m}]^{-} - [q^{m}]^{-}|_{L^2} \leq |q_{\epsilon}^{m} - q^{m}|_{L^2}, \quad |[q_{\epsilon}^{m} - 1]^{+} - [q^{m} - 1]^{+}|_{L^2} \leq |q_{\epsilon}^{m} - q^{m}|_{L^2}. $$

Consequently, with (42) we have $[q_{\epsilon}^m]^{-} \to [q^m]^{-}$ and $[q_{\epsilon}^m - 1]^+ \to [q^m - 1]^+$ in $H$. As $k > 0$ is a fixed number, we can pass to the limit on $\epsilon$ in (44) to obtain that

$$ \sum_{m=1}^{N} (|[q^m]^{-}|_{L^2}^2 + |[q^m - 1]^+|_{L^2}^2) = 0, $$

which implies

$$ 0 \leq q^m \leq 1, \text{ a.e. in } \mathbf{x} \in \mathcal{M}, \text{ i.e., } q^m \in \mathcal{K}. \qquad (45) $$

With the above preparations, we could pass to the limit from the scheme (33)-(34) to the scheme (31)-(32) term by term. We only point out the following three subtle points. First, here we obtain the strong convergence of the functions $U_{\epsilon}^m$ to their limits $U^m$ in $H$ by the compact embedding $V \subset H$. The strong convergence will imply the a.e. convergence of $U_{\epsilon}^m$ in $\mathcal{M}$ up to subsequences. The strong convergence and the a.e. convergence of $U_{\epsilon}^m$ are crucial when we pass to the limit in the nonlinear terms in the scheme (33)-(34). The reason is well-known: nonlinear mappings do not preserve weak convergences. Second, for the two penalization terms, we have for $q^b \in \mathcal{K}$ that

$$ \left(\frac{1}{\epsilon_1}[q_{\epsilon}^{m}]^{-}, q^{b} - q_{\epsilon}^{m}\right) \geq 0, \quad -\left(\frac{1}{\epsilon_1}[q_{\epsilon}^{m} - 1]^{+}, q^{b} - q_{\epsilon}^{m}\right) \geq 0. \qquad (46) $$
Due to weak convergence of $q^m$ in $V$ and the weak lower semi-continuity property of the norm, we also have

$$ 
\begin{aligned}
\limsup a_q(q_{\varepsilon}^m, q^b - q_{\varepsilon}^m) &= \lim a_q(q_{\varepsilon}^m, q^b) - \liminf a_q(q_{\varepsilon}^m, q_{\varepsilon}^m) \\
&\le a_q(q^m, q^b) - a_q(q^m, q^m) \qquad (47) \\
&= a_q(q^m, q^b - q^m).
\end{aligned}
$$

All the above three terms produce correct direction of inequalities during the limit process. Therefore, we obtain the desired variational inequalities in (32). Third, in order to show that $h_{Q^m} \in \mathcal{H}(Q^m - q_s)$, we shall use the idea of subdifferential for convex functions.

Summarizing the above arguments, we obtain from (33)-(34) via passing to the limit on $\varepsilon$ the existence of a solution $(T^m, q^m)$ to (31)-(32).

# 4 Convergence of the Euler Scheme

In this section, we want to prove the convergence of the solutions of the Euler scheme (31)-(32) to the solutions of the system (28)-(30). We shall use the same conventions on subsequences and indices as in the last section, that is, the limit process in this part is $N \to +\infty$ or equivalently $k \to 0+$ and up to subsequences.

Due to the weak lower semi-continuity property of the norms, we know that for the limit functions $U^m$ which now have no dependence on $\varepsilon$, the bounds in Lemmas 2 and 4 are now valid with $U_{\varepsilon}^m$ replaced by the limit functions $U^m$.

## 4.1 Construction of approximations

For each fixed $k$ (or $N$), we associate to the elements $U^0, U^1, U^2, \dots, U^N$ the following approximate functions $U_k = (T_k, q_k), \tilde{U}_k = (\tilde{T}_k, \tilde{q}_k)$ and $W_k = (\mathcal{T}_k, \mathcal{Q}_k)$ which are defined piecewise on $[0, t_1]$ and take values in the space $V^2$:

$$ U_k(t) = U^m, \quad \tilde{U}_k(t) = U^{m-1}, \quad \text{for } t \in [(m-1)k, mk), \quad m = 1, 2, \dots, N. \quad (48) $$

$$ W_k(t) = \frac{U^m - U^{m-1}}{k} (t - (m-1)k) + U^{m-1}, \quad \text{for } t \in [(m-1)k, mk), \quad m = 1, 2, \dots, N. \quad (49) $$

## 4.2 Reinterpretation of a priori estimates

First, we give a lemma measuring the distance in $L^2(0, t_1; H)$ of the functions $U_k, \tilde{U}_k$ and $W_k$ in the limit process $k \to 0+$.

**Lemma 6.** For the functions $U_k, W_k$ and $\tilde{U}_k$ defined above, there hold

$$ |U_k - W_k|_{L^2(0,t_1;H)} \le C(\mathbf{u}, U_0, t_1)\sqrt{k}, \quad |U_k - \tilde{U}_k|_{L^2(0,t_1;H)} \le C(\mathbf{u}, U_0, t_1)\sqrt{k}. $$

Now, we state a result concerning the boundedness of the functions $U_k, \tilde{U}_k$ and $W_k$.
**Lemma 7.** The functions $U_k, \tilde{U}_k$ and $W_k$ remain in a bounded set of $L^2(0, t_1; V) \cap L^\infty(0, t_1; H)$ as $k \to 0^+$. The functions $\partial_t W_k$ form a bounded set in $L^2(0, t_1; V^*)$ and $U_k - W_k \to 0$ in $L^2(0, t_1; H)$ strongly as $k \to 0^+$.

Lemmas 6 and 7 can be regarded as reinterpretations of Lemmas 2 and 4 in terms of the functions $U_k, W_k, \tilde{U}_k$.

Define $\mathbf{u}_k : [0, t_1] \to \mathbf{V}$ as follows:

$$ \mathbf{u}_k(t) = \mathbf{u}^m, \text{ for } t \in [(m-1)k, mk), m = 1, 2, \dots, N. \quad (50) $$

We have the following classical lemma:

**Lemma 8 (Convergence of $\mathbf{u}_k$).** For the functions $\mathbf{u}_k$ defined above, there holds

$$ \mathbf{u}_k \to \mathbf{u}, \text{ in } L^r(0, t_1, \mathbf{V}) \text{ as } k \to 0^+. \quad (51) $$

For later use, we also define the linear averaging map for the test functions $U^b = (T^b, q^b) \in L^2(0, t_1; V)$ that we will use below, that is, we define $U_k^b : [0, t_1] \to V^2$ piecewise by

$$ U_k^b(t) = \frac{1}{k} \int_{(m-1)k}^{mk} U^b(t) \, dt \text{ on } [(m-1)k, mk). $$

Similarly as in Lemma 8, we conclude that $U_k^b \to U^b$ strongly in $L^2(0, t_1; V^2)$ as $k \to 0$. Moreover, if $q^b \in \mathcal{K}$ for a.e. $t \in [0, t_1]$, we have $q_k^b \in \mathcal{K}$ for all $t \in [0, t_1]$.

### 4.3 Passage to the limit: $k \to 0^+$

We first reinterpret as follows the scheme (31)-(32) in terms of the functions $U_k = (T_k, q_k), \tilde{U}_k = (\tilde{T}_k, \tilde{q}_k), W_k = (\mathcal{T}_k, \mathcal{Q}_k)$ and $U_k^b = (T_k^b, q_k^b)$:

$$ \begin{aligned} & \langle \partial_t \mathcal{T}_k, T_k^b \rangle + a_T(T_k, T_k^b) + b_T(\mathbf{u}_k, T_k, T_k^b) - d(\omega_k, \tilde{T}_k, T_k^b) - l_T(T_k^b) \\ & \qquad = \left(\frac{1}{p}[\omega_k]^{-} h_{Q_k} \varphi(\tilde{T}_k), T_k^b \right), \end{aligned} \quad (52) $$

$$ \begin{aligned} & \langle \partial_t \mathcal{Q}_k, q_k^b - q_k \rangle + a_q(q_k, q_k^b - q_k) + b_q(\mathbf{u}_k, q_k, q_k^b - q_k) - l_q(q_k^b - q_k) \\ & \qquad \geq \left( -\frac{1}{p}[\omega_k]^{-} h_{Q_k} F(\tilde{T}_k), q_k^b - q_k \right), \end{aligned} \quad (53) $$

where $Q_k$ is either $\tilde{q}_k$ or $q_k$. Furthermore, $h_{Q_k}$ is defined by $h_{Q_k}(t) = h_{Q^m}$ when $t \in [(m-1)k, mk)$. Here we emphasize that we require $q_k^b \in L^2(0, t_1; \mathcal{K})$.

Due to Lemma 7, we have, up to subsequences, in the limit $k \to 0^+$, that

$$ U_k \to U = (T, q), \text{ weakly in } L^2(0, t_1; V) \text{ and weak-* in } L^\infty(0, t_1; H), \quad (54) $$

$$ W_k \to W = (\mathcal{T}, \mathcal{Q}), \text{ weakly in } L^2(0, t_1; V) \text{ and weak-* in } L^\infty(0, t_1; H), \quad (55) $$

and

$$ \partial_t W_k \to \partial_t W = (\partial_t \mathcal{T}, \partial_t \mathcal{Q}), \text{ weakly in } L^2(0, t_1; V^*). \quad (56) $$
Obviously, $\tilde{U}_k = U(\cdot - k)$ converges also to $U$ in $L^2(0, t_1; V)$ weakly and in $L^\infty(0, t_1; H)$ weak-*.

In view of Lemma 6, we know that

$$U = W. \qquad (57)$$

Now, we consider the inclusions $V \subset H \subset V^*$ where the first inclusion is compact and the second inclusion is continuous. In view of (55) and (56), we conclude, by applying the Aubin-Lions compactness theorem, that

$$W_k \to W, \text{ strongly in } L^2(0, t_1; H). \qquad (58)$$

By Lemma 6 again, we conclude that

$$U_k, \tilde{U}_k, W_k \to U, \text{ strongly in } L^2(0, t_1; H). \qquad (59)$$

With the above preparations, we can now pass to the limit $k \to 0+$ from (52)-(53) to (28)-(30) term by term. Here we also point several subtle points. First, we obtain the strong convergence, i.e. (59), in this step by the Aubin-Lions compactness argument. For further details, see [28]. Second, the limit function $q$ satisfies the required range condition. Indeed, regarded as a convex subset of $L^2(0, t_1; V)$, $L^2(0, t_1; K)$ is closed with respect to the strong topology induced by the $L^2(0, t_1; V)$-norm. Therefore, it is also closed with respect to the weak topology. Furthermore, in view of the fact that $q_k \in L^2(0, t_1; K)$ which is obvious from the definition and that $q_k$ converge to $q$ weakly in $L^2(0, t_1; V)$, we conclude that $q \in L^2(0, t_1; K)$. Third, the subtle point in this passage to the limit is to deal with the term $\int_0^{t_1} \langle \partial_t Q_k, q_k^b - q_k \rangle dt$ which is the sum of $\int_0^{t_1} \langle \partial_t Q_k, q_k^b - Q_k \rangle dt$ and $\int_0^{t_1} \langle \partial_t Q_k, Q_k - q_k \rangle dt$. Using integration by parts, (56) and the lower semi-continuity of the norm, we write

$$
\begin{align*}
\limsup \int_0^{t_1} \langle \partial_t Q_k, q_k^b - Q_k \rangle dt &= -\liminf \int_0^{t_1} \langle \partial_t Q_k, Q_k \rangle dt + \lim \int_0^{t_1} \langle \partial_t Q_k, q_k^b \rangle dt \\
&= -\liminf \frac{1}{2} |Q_k(t_1)|_{L^2}^2 + \frac{1}{2} |q_0|_{L^2}^2 + \int_0^{t_1} \langle \partial_t q, q^b \rangle dt \\
&\leq -\frac{1}{2} |q(t_1)|_{L^2}^2 + \frac{1}{2} |q_0|_{L^2}^2 + \int_0^{t_1} \langle \partial_t q, q^b \rangle dt \\
&= - \int_0^{t_1} \langle \partial_t q, q \rangle dt + \int_0^{t_1} \langle \partial_t q, q^b \rangle dt \\
&= \int_0^{t_1} \langle \partial_t q, q^b - q \rangle dt. \tag{60}
\end{align*}
$$

where we have used, in the second equality of (60), the observation

$$\lim \int_0^{t_1} \langle \partial_t Q_k, q_k^b \rangle dt \to \int_0^{t_1} \langle \partial_t q, q^b \rangle dt, \qquad (61)$$
which is a simple consequence of (56) and the strong convergence of $q_k^b$ to $q^b$ in $L^2(0, t_1; V)$.

A subtle point is the treatment of $\int_0^{t_1} \langle \partial_t Q_k, Q_k - q_k \rangle dt$. Though we have (56) (which implies in particular that $\partial_t Q_k$ is bounded in $L^2(0, t_1; V^*)$ and $Q_k - q_k \to 0$ weakly in $L^2(0, t_1; V)$, we can not conclude that the limit of $\int_0^{t_1} \langle \partial_t Q_k, Q_k - q_k \rangle dt$ is 0. Rather, we show, by the specific forms of $Q_k$ and $q_k$, that

$$ \limsup_{t_0} \int_{t_0}^{t_1} \langle \partial_t Q_k, Q_k - q_k \rangle dt \le 0. \quad (62) $$

Indeed, noticing that $\partial_t Q_k = \frac{q^m - q^{m-1}}{k}$ and $Q_k - q_k = \frac{q^m - q^{m-1}}{k}(t - mk)$ on the subinterval $[(m-1)k, mk)$ of $[0, t_1]$, we have:

$$ \begin{aligned} \int_0^{t_1} \langle \partial_t Q_k, Q_k - q_k \rangle dt &= \sum_{m=1}^N \int_{(m-1)t}^{mt} \langle \partial_t Q_k, Q_k - q_k \rangle dt \\ &= \sum_{m=1}^N \int_{(m-1)t}^{mk} \langle \frac{q^m - q^{m-1}}{k}, \frac{q^m - q^{m-1}}{k} (t - mk) \rangle dt \\ &= \sum_{m=1}^N \int_{(m-1)t}^{mk} \frac{|q^m - q^{m-1}|_{L^2}^2}{k^2} (t - mk) dt \\ &\le 0, \end{aligned} $$

which implies (62). From (60) and (62), we can conclude that

$$ \limsup_{t_0} \int_{t_0}^{t_1} \langle \partial_t Q_k, q_k^b - q_k \rangle dt \le \int_{t_0}^{t_1} \langle \partial_t q, q^b - q \rangle dt. \quad (63) $$

To sum up, we have proved the following theorem when $q_s$ is constant:

**Theorem 1.** Given $T_0, q_0 \in H$ with $0 \le q_0 \le 1$ a.e. in $M$, the Euler scheme (31)-(32) contains a subsequence which converges to a solution of the system (1)-(7).

# 5 The case where $q_s$ depends on time

In this part, we extend our former results to the case where $q_s$ is not constant.

## 5.1 The nonlinearities $\varphi$ and $F$

When $q_s$ evolves according to (8), the nonlinearities $F$ and $\varphi$ will also depend on the function $q_s$. The nonlinearities (see e.g., [17, 18]) $\varphi$ and $F: \mathbb{R} \times \mathbb{R} \to \mathbb{R}$ are defined as follows:

$$ F(T, q_s) = q_s G(T, q_s) = q_s T \frac{RL(T) - c_p R_v T}{c_p R_v T^2 + q_s L(T)^2}, \quad (64) $$
$$ \varphi(T, q_s) = \frac{L(T)}{c_p} F(T, q_s), \qquad (65) $$

where

$$ L(T) = c_1 - c_2 T, \quad G(T, q_s) = T \frac{RL(T) - c_p R_v T}{c_p R_v T^2 + q_s L(T)^2}. \qquad (66) $$

In the above, $c_1, c_2, R, c_p, R_v$ are all strictly positive constants. The additional dependence of $F$ and $\varphi$ on $q_s$ will bring us technical complexities during the passages to the limits.

Notice that the functions $F, G$, and $\varphi$ have a singularity at $(0,0)$; $F(T, q_s)$ is bounded but discontinuous at $(0,0)$ and $G(T, q_s)$ may blow up at $(0,0)$. To overcome this difficulty, we introduce the following regularized version $\varphi_r, F_r$ and $G_r$ for $\varphi, F$ and $G$.

$$ F_r(T, q_s) = q_s G_r(T, q_s) = q_s T \frac{RL(T) - c_p R_v T}{c_p R_v \max(T, \gamma)^2 + q_s L(T)^2}, \qquad (67) $$

$$ \varphi_r(T, q_s) = \frac{L(T)}{c_p} F_r(T, q_s), \qquad (68) $$

$$ G_r(T, q_s) = T \frac{RL(T) - c_p R_v T}{c_p R_v \max(T, \gamma)^2 + q_s L(T)^2}, \qquad (69) $$

where $\gamma > 0$ is smaller than any temperature on Earth. Once arriving at $q_s \ge 0$, we can derive $G_r$ is nonnegative by observing that $T \le \frac{LR}{c_p R_v}$ for any temperature on Earth. It is easy to see the rational function $F_r$ is bounded and globally Lipschitz on $\mathbb{R} \times \mathbb{R}$, i.e.,

$$ |F_r(\zeta_1, \xi_1) - F_r(\zeta_2, \xi_2)| \le C(|\zeta_1 - \zeta_2| + |\xi_1 - \xi_2|), \quad \forall \zeta_1, \zeta_2 \in \mathbb{R}, \xi_1, \xi_2 \in [0, \infty), \qquad (70) $$

and

$$ |F_r(\zeta, \xi)| \le C, \quad \forall \zeta \in \mathbb{R}, \xi \in [0, \infty). \qquad (71) $$

The function $\varphi_r$ is also globally Lipschitz,

$$ |\varphi_r(\zeta_1, \xi_1) - \varphi_r(\zeta_2, \xi_2)| \le C(|\zeta_1 - \zeta_2| + |\xi_1 - \xi_2|), \quad \forall \zeta_1, \zeta_2 \in \mathbb{R}, \xi_1, \xi_2 \in [0, \infty). \qquad (72) $$

In addition, as $F_r(0,0) = 0$, we have $\varphi_r(0,0) = 0$. Hence the Lipschitz function $\varphi_r$ also satisfies $|\varphi_r(\zeta, \xi)| \le C(|\zeta| + |\xi|)$.

**Definition 2.** Let $(T_0, q_0, q_{s,0}) \in H \times H \times H$ be such that $0 \le q_0 \le 1$, $0 < q_{s,0} < 1$ a.e. in $M$ and let $t_1 > 0$ be fixed. A vector $(T, q, q_s) \in L^2(0, t_1; V \times V) \cap C([0, t_1]; H \times H) \times L^\infty(M \times [0, t_1])$ with $0 < q_s < 1$, $(\partial_t T, \partial_t q, \partial_t q_s) \in L^2(0, t_1; V^* \times V^*) \times L^\infty(M \times [0, t_1])$ is a solution to the initial boundary value problem described by (1), (2), (8), (6), (7) and (9), if for a.e. $t \in [0, t_1]$ and for every $(T^b, q^b) \in V \times K$, we have

$$ \langle \partial_t T, T^b \rangle + a_T(T, T^b) + b_T(u, T, T^b) - d(\omega, T, T^b) - l_T(T^b) = \left( \frac{1}{p} \omega^{-}(t) h_q \varphi(T, q_s), T^b \right), \qquad (73) $$
$$ \langle \partial_t q, q^b - q \rangle + a_q(q, q^b - q) + b_q(\mathbf{u}, q, q^b - q) - l_q(q^b - q) \geq \left(-\frac{1}{p}\omega^-(t)h_q F(T, q_s), q^b - q\right), \quad (74) $$

$$ \frac{dq_s}{dt} = -\frac{1}{p}\omega^{-}h_q F(T, q_s) \quad \text{for a.e. } (t, \mathbf{x}) \in [0, t_1] \times \mathcal{M}, \quad (75) $$

for some $h_q \in \mathcal{H}(q - q_s)$ and

$$ U_0 = (T_0, q_0, q_{s,0}). \quad (76) $$

## 5.2 The discretization scheme

We begin with

$$ (T^0, q^0, q_s^0) := (T_0, q_0, q_{s,0}), \quad \text{i.e., the given initial datum.} \quad (77) $$

When $(T^0, q^0, q_s^0)$, $(T^1, q^1, q_s^1)$, ..., $(T^{m-1}, q^{m-1}, q_s^{m-1})$ are known, $T^m \in V$ and $q^m \in \mathcal{K}$ are determined by:

$$ \begin{aligned} & \left\langle \frac{T^m - T^{m-1}}{k}, T^b \right\rangle + a_T(T^m, T^b) + b_T(\mathbf{u}^m, T^m, T^b) - d(\omega^m, T^{m-1}, T^b) - l_T(T^b) \\ & = \left( \frac{1}{p} [\omega^m]^T - h_{q^{m-1}} \varphi(T^{m-1}, q_s^{m-1}), T^b \right), \end{aligned} \quad (78) $$

$$ \begin{aligned} & \left\langle \frac{q^m - q^{m-1}}{k}, q^b - q^m \right\rangle + a_q(q^m, q^b - q^m) + b_q(\mathbf{u}^m, q^m, q^b - q^m) - l_q(q^b - q^m) \\ & \geq \left( -\frac{1}{p} [\omega^m]^T - h_{q^{m-1}} F(T^{m-1}, q_s^{m-1}), q^b - q^m \right), \end{aligned} \quad (79) $$

and

$$ q_s^m := Z^m(mk), \quad (80) $$

where $Z^m(t)$ is the solution to the following initial value problem

$$ \begin{cases} \frac{dZ^m(t)}{dt} = -\frac{1}{p}[\omega^m]^T - h_{q^{m-1}} F(T^{m-1}, Z^m(t)), \\ Z^m((m-1)k) = q_s^{m-1}, \end{cases} \quad (81) $$

where $h_{q^{m-1}} \in H(q^{m-1} - q_s^{m-1})$.

**Remark 2.** We shall point out the distinct feature of the discretization for the $q_s$-equation. First, the $q_s$-equation shall be discretized once the $(T, q)$-equation is discretized as the $(T, q)$-equation depends on $q_s$ through the nonlinear functions $F(T, q_s)$ and $\varphi(T, q_s)$ and we do not allow time dependence in the discretized equation on $(T, q)$. Here we did not use the standard Euler algorithm for ordinary differential equations (ODEs). Rather, we define $q_s^m$ by successively solving the ordinary differential equation (81) and make evaluations at specific time points.
The latter will make the $q_s^m$'s satisfy the required range condition inherited from that of $q_s$. The Euler algorithm for ODEs cannot guarantee this range condition. Finally, we emphasize that we need to use $h_{q^{m-1}}$ in the scheme (78), (79) and (81). This choice makes the $(T^m, q^m, Z^m)$-equation and the $q_s^m$-equation decouple.

To show the existence of a solution $(T^m, q^m, Z^m)$ to (78)-(81), we consider the associated regularized and penalized problem:

To find $T_\epsilon^m, q_\epsilon^m \in V$ and $Z_\epsilon^m \in H$ such that

$$ \begin{aligned} & \left\langle \frac{T_{\epsilon}^{m} - T_{\epsilon}^{m-1}}{k}, T^{b} \right\rangle + a_{T}(T_{\epsilon}^{m}, T^{b}) + b_{T}(\mathbf{u}^{m}, T_{\epsilon}^{m}, T^{b}) - d(\omega^{m}, T_{\epsilon}^{m-1}, T^{b}) - l_{T}(T^{b}) \\ & = \left(\frac{1}{p}[\omega^{m}]^{-} H_{\epsilon 2}(q_{\epsilon}^{m-1} - q_{s,\epsilon}^{m-1})\varphi_{r}(T_{\epsilon}^{m-1}, q_{s,\epsilon}^{m-1}), T^{b}\right), \end{aligned} \quad (82) $$

$$ \begin{aligned} & \left\langle \frac{q_{\epsilon}^{m} - q_{\epsilon}^{m-1}}{k}, q^{b} \right\rangle + a_{q}(q_{\epsilon}^{m}, q^{b}) + b_{q}(\mathbf{u}^{m}, q_{\epsilon}^{m}, q^{b}) - l_{q}(q^{b}) \\ & = \left(\frac{1}{\epsilon_1}[q_{\epsilon}^{m}]^{-}, q^{b}\right) - \left(\frac{1}{\epsilon_1}[q_{\epsilon}^{m} - 1]^{+}, q^{b}\right) - \left(\frac{1}{p}[\omega^{m}]^{-} H_{\epsilon 2}(q_{\epsilon}^{m-1} - q_{s,\epsilon}^{m-1})F_{r}(T_{\epsilon}^{m-1}, q_{s,\epsilon}^{m-1}), q^{b}\right), \end{aligned} \quad (83) $$

and

$$ q_{s,\epsilon}^{m} := Z_{\epsilon}^{m}(mk), \quad (84) $$

where $Z_\epsilon^m$ is the solution to the following initial value problem

$$ \left\{ \begin{array}{l} \frac{dZ_{\epsilon}^{m}(t)}{dt} = -\frac{1}{p}[\omega^{m}]^{-} H_{\epsilon 2}(q_{\epsilon}^{m-1} - q_{s,\epsilon}^{m-1})F_{r}(T_{\epsilon}^{m-1}, Z_{\epsilon}^{m}(t)), \\ Z_{\epsilon}^{m}((m-1)k) = q_{s,\epsilon}^{m-1}. \end{array} \right. \quad (85) $$

## 5.3 Validity of the iteration

In order to have a valid scheme, we should be able to obtain $(T_\epsilon^m, q_\epsilon^m, q_{s,\epsilon}^m)$ when $(T_\epsilon^{m-1}, q_\epsilon^{m-1}, q_{s,\epsilon}^{m-1})$ are known. This is true. First, observe that the equation on $(T_\epsilon^m, q_\epsilon^m)$ and $q_{s,\epsilon}^m$ are decoupled. Noticing that the arguments of $F_r$ and $\varphi_r$ have index $m-1$ in (82) and (83), we can obtain $(T_\epsilon^m, q_\epsilon^m)$ either by the Minty-Browder Theorem or by the Galerkin method when $(T_\epsilon^{m-1}, q_\epsilon^{m-1}, q_{s,\epsilon}^{m-1})$ are known. The determination of $q_{s,\epsilon}^m$ is easy when $(T_\epsilon^{m-1}, q_\epsilon^{m-1}, q_{s,\epsilon}^{m-1})$ are known: we just need to use the initial value problem (85). The initial value problem (85) admits a unique solution $Z_\epsilon^m(t)$ on $[(m-1)k, T^*)$. The maximum time of existence $T^* = +\infty$ due to the form of $F_r$. Actually, we have

$$ F_r(T_{\epsilon}^{m-1}, Z_{\epsilon}^{m}(t)) = Z_{\epsilon}^{m}(t)G_r(T_{\epsilon}^{m-1}, Z_{\epsilon}^{m}(t)). $$

Therefore, we have the following integral form of (85):

$$ q_{s,\epsilon}^{m}(t) = q_{\epsilon}^{m-1} \exp \left\{ - \int_{(m-1)k}^{t} \frac{1}{p} [\omega^{m}]^{-} H_{\epsilon 2}(q_{\epsilon}^{m-1} - q_{s,\epsilon}^{m-1}) G_{r}(T_{\epsilon}^{m-1}, Z_{\epsilon}^{m}(t)) d\tau \right\}. \quad (86) $$
Due to the above expression of $q_{s,\epsilon}^m(t)$, no blow-up can happen as we always have

$$|Z_{\epsilon}^{m}(t)|_{L^{2}} \leq |q_{s,\epsilon}^{m-1}|_{L^{2}}, \quad |Z_{\epsilon}^{m}(t)|_{L^{\infty}} \leq |q_{s,\epsilon}^{m-1}|_{L^{\infty}} \text{ for } t \geq (m-1)k. \quad (87)$$

For our purpose, we just need $T^* \geq mk$ so that $q_{s,\epsilon}^m$ is well-defined.

## 5.4 A priori estimates for $(T_{\epsilon}^m, q_{\epsilon}^m, q_{s,\epsilon}^m)$

Integrating the first equation in (85) from $(m-1)k$ to $mk$, we find

$$q_{s,\epsilon}^{m} = q_{s,\epsilon}^{m-1} \exp\left\{-\int_{(m-1)k}^{mk} \frac{1}{p} [\omega^m]^{-H_{\epsilon_2}(q_{\epsilon}^{m-1} - q_{s,\epsilon}^{m-1})} G_r(T_{\epsilon}^{m-1}, Z_{\epsilon}^{m}(t)) d\tau\right\}. \quad (88)$$

As $0 < q_s^0 = q_{s,0} < 1$, we easily conclude that $0 < q_{s,\epsilon}^m < 1$ for all $1 \leq m \leq N$ by repeatedly using (88) for $m = 1, 2, \dots, N$. Actually, we have the following pointwise monotone relations

$$0 < q_{s,\epsilon}^{m}(\mathbf{x}) \leq q_{s,\epsilon}^{m-1}(\mathbf{x}) < 1, \quad m = 1, 2, \dots, N. \quad (89)$$

Now we aim to obtain a priori estimates on $(T_{\epsilon}^m, q_{\epsilon}^m)$ independent of $k$ and $\epsilon$ for the regularized and penalized problem (82)-(83). Due to the form of $\varphi_r$, $F_r$ and our estimate on $q_{s,\epsilon}^m$, we know that

$$|\varphi_r(T_{\epsilon}^{m-1}, q_{s,\epsilon}^{m-1})| \leq C(|T_{\epsilon}^{m-1}| + 1) \text{ and } |F_r(T_{\epsilon}^{m-1}, q_{s,\epsilon}^{m-1})| \leq C, \quad (90)$$

where $C$ is constant independent of $\epsilon$ and $k$. Therefore, Lemma 2 is still valid.

Now, we explain that Lemma 4 is still valid even when $q_s$ is not constant. Since $|F_r(T_{\epsilon}^{m-1}, q_{s,\epsilon}^{m-1})|$ can be bounded by a universal constant, we still have, for the penalization terms, that Lemma 5 holds. Then estimating the duality pair $\langle \frac{U_{\epsilon}^m - U_{\epsilon}^{m-1}}{k}, U^b \rangle$ where $U^b \in V^2$, we can derive the conclusion of Lemma 4.

## 5.5 Passage to the limit $\epsilon \to 0+$ and $k \to 0+$

The passage to the limit can be proceeded essentially as before. By the compact Sobolev embedding theorem as before, we can derive that $(T_{\epsilon}^m, q_{\epsilon}^m) \to (T^m, q^m)$ strongly in $H$. In addition, by the above estimates on $q_{s,\epsilon}^m$, we know that there exist $q_s^m \in H \cap L^\infty(\mathcal{M})$ such that, up to subsequences,

$$q_{s,\epsilon}^{m} \to q_{s}^{m} \text{ weakly in } H \text{ and weak-}^{*} \text{ in } L^{\infty}(\mathcal{M}). \quad (91)$$

During the passage to the limit: $k \to 0+$, besides defining $U_k, \tilde{U}_k, W_k, h_{\tilde{q}_k}$ as before, we should also define $q_{s,k} : [0, t_1] \to L^\infty(\mathcal{M}) \cap L^2(\mathcal{M})$ as follows:

$$q_{s,k}(t) = \begin{cases} Z^1(t), & \text{when } t \in [0, k), \\ Z^2(t), & \text{when } t \in [k, 2k), \\ \dots, & \\ Z^m(t), & \text{when } t \in [(m-1)k, mk), \\ \dots, & \\ Z^N(t), & \text{when } t \in [(N-1)k, Nk). \end{cases} \quad (92)$$
By our definition of $q_{s,k}$, we know that $q_{s,k}(0) = q_{s,0}$ and $q_{s,k}$ is piecewise differentiable and satisfies (84) on the whole interval $[0, t_1]$. Then, we can rein-terpret our scheme in terms of the above functions $U_k = (T_k, q_k)$, $\tilde{U}_k = (\tilde{T}_k, \tilde{q}_k)$, $W_k = (\mathcal{T}_k, \mathcal{Q}_k)$ and $U_k^b = (T_k^b, q_k^b)$ as follows

$$ \begin{aligned} & \langle \partial_t \mathcal{T}_k, T_k^b \rangle + a_T(T_k, T_k^b) + b_T(\mathbf{u}_k, T_k, T_k^b) - d(\omega_k, \tilde{T}_k, T_k^b) - l_T(T_k^b) \\ & = \left(\frac{1}{p}[\omega_k]\right)^{-} h_{\tilde{q}_k} \varphi_r(\tilde{T}_k, \tilde{q}_{s,k}), T_k^b, \end{aligned} \quad (93) $$

$$ \begin{aligned} & \langle \partial_t \mathcal{Q}_k, q_k^b - q_k \rangle + a_q(q_k, q_k^b - q_k) + b_q(\mathbf{u}_k, q_k, q_k^b - q_k) - l_q(q_k^b - q_k) \\ & \geq \left(-\frac{1}{p}[\omega_k]\right)^{-} h_{\tilde{q}_k} F_r(\tilde{T}_k, \tilde{q}_{s,k}), q_k^b - q_k, \end{aligned} \quad (94) $$

$$ \frac{dq_{s,k}}{dt} = -\frac{1}{p}[\omega_k]^{-} h_{\tilde{q}_k} F_r(\tilde{T}_k, q_{s,k}). \quad (95) $$

The two passages to the limit processes $\varepsilon \to 0+$ and $k \to 0+$ are paral-lel to those in the case that $q_s$ is constant. However, we shall pay attention to the convergences of the terms involving $\varphi_r(T_{\varepsilon}^{m-1}, q_{s,\varepsilon}^{m-1})$, $F_r(T_{\varepsilon}^{m-1}, q_{s,\varepsilon}^{m-1})$ dur-ing the passage to the limit $\varepsilon \to (0+, +)$, and $\varphi_r(\tilde{T}_k, \tilde{q}_{s,k})$, $F_r(\tilde{T}_k, \tilde{q}_{s,k})$ during the passage to the limit $k \to 0+$ respectively. In both of the two limit pro-cesses, we could achieve strong convergences for the sequences $(T_{\varepsilon}^m)_{\varepsilon}$ and $(\tilde{T}_k)_k$ up to subsequences by compact Sobolev embedding and Aubin-Lions argument respectively. While for $(q_{s,\varepsilon}^m)_{\varepsilon}$ and $(\tilde{q}_{s,k})_k$, we could not obtain strong conver-gences. Fortunately, we could obtain a.e. convergences for them for $\mathbf{x} \in \mathcal{M}$ and $(t, \mathbf{x}) \in [0, t_1] \times \mathcal{M}$ respectively. These pointwise convergences, together with the Lebesgue Dominated Convergence Theorem, will enable us to pass to the limits for these terms. The a.e. convergences are guaranteed by the following lemma.

**Lemma 9.** Consider the following equation

$$ \left\{ \begin{array}{l} \frac{dq_s^j(t)}{dt} = F(T^j(t), q_s^j(t)), \\ q_s^j(0) = q_{s,0}, \end{array} \right. \quad (96) $$

where $F(\cdot, \cdot)$ is a real-valued bounded Lipschitz function. Suppose $T^j = T^j(x, t)$ converges to some $T = T(x, t)$ strongly in $L^2(0, l; H)$ as $j \to \infty$. Then up to subsequences, we have

$$ q_s^j(x, t) \to q_s(x, t) \quad \text{for any } t \in [0, l] \quad \text{and a.e. } x \in \Omega \setminus \Omega_0, \quad (97) $$

where $q_s(x, t)$ is the solution to the initial value problem

$$ \left\{ \begin{array}{l} \frac{dq_s(t)}{dt} = F(T(t), q_s(t)), \\ q_s(0) = q_{s,0}, \end{array} \right. \quad (98) $$

and $\Omega_0$ is a subset of $\mathcal{M}$ which has Lebesgue measure 0 and is independent of $t$.
*Proof.* By assumption, we have up to a subsequence that

$$T^j(x, t) \to T(x, t) \text{ a.e. } (x, t).$$

Writing the ODE on $q_s^j$ in the integral equation form, we have

$$q_s^j(t) = q_s^j(0) + \int_0^t F(T^j(x, \tau), q_s^j(x, \tau)) d\tau. \quad (99)$$

By assumption, we have

$$\int_0^l \int_M |T^j(x,t) - T(x,t)|^2 dxdt = \int_M \left( \int_0^l |T^j(x,t) - T(x,t)|^2 dt \right) dx \to 0, \text{ as } j \to \infty.$$

By Fubini's theorem, we have for a.e. $x$ that

$$\int_0^l |T^j(x,t) - T(x,t)|^2 dt \to 0.$$

Now take $x_*$ such that

$$\int_0^l |T^j(x_*, t) - T(x_*, t)|^2 dt \to 0. \quad (100)$$

Notice that

$$\left\{ \begin{array}{l} \frac{dq_s^j(x_*, t)}{dt} = F(T^j(x_*, t), q_s^j(x_*, t)), \\ q_s^j(x_*, 0) = q_s(x_*, 0). \end{array} \right. \quad (101)$$

As $|\frac{dq_s^j}{dt}|_{L^\infty(0,l)} \le C$, we have up to a subsequence and with $x_*$ fixed that $q_s^{ji} \to q_s$ uniformly on $[0, l]$.

Since $F(T^{ji}(x_*, t))$ and $q_s^{ji}(x_*, t)$ converge for a.e. $t \in [0, l]$, we can pass to the limit in the following equation

$$q_s^{ji}(x_*, t) = q_s^j(x_*, 0) + \int_0^t F(T^{ji}(x_*, \tau), q_s^{ji}(x_*, \tau)) d\tau,$$

and by the Lebesgue Dominated Convergence Theorem, we have

$$q_s(x_*, t) = q_s(x_*, 0) + \int_0^t F(T(x_*, \tau), q_s(x_*, \tau)) d\tau. \quad (102)$$

By uniqueness of the solution to (99) and the fact that the limit solution is independent of the subsequence, we conclude that

$$q_s^j(\cdot, x_*) \to q_s(\cdot, x_*), \text{ uniformly in } t,$$

for any $x_*$ such that (5.5) holds, i.e., for $x_* \in \Omega \setminus \Omega_0$ where $|\Omega_0| = 0$, Hence,

$$q_s^j(x, t) \to q_s(x, t) \text{ for any } t \text{ and for a.e. } x \in \Omega \setminus \Omega_0. \quad (103)$$
**Remark 3.** A particular case of the above lemma is that $T^j = T^j(x)$, i.e., $T^j$ has no dependence on the time variable $t \in [0, l]$. In this case we can regard $T^j$ as a constant function of $t \in [0, l]$ and consequently, we still have $T^j \to T$ in $L^2(0, l; H)$.

In summary, when $q_s$ depends on time, we have the following theorem.

**Theorem 2.** Given $T_0, q_0, q_{s,0} \in H$ with $0 \le q_0 \le 1$, and $0 < q_{s,0} < 1$ a.e. in $M$, the scheme (78)-(81) contains a subsequence which converges to a solution of the system (1)-(9).

**Acknowledgments** The research was partially supported by the National Science Foundation under the grants NSF-DMS-1206438, NSF-DMS-1510249 and by the Research Fund of Indiana University.

## References

1. A. Bousquet, M. Coti Zelati, R. Temam, *Phase transition models in atmospheric dynamics*, Milan J. Math. 82 (2014) 99-128.

2. H. Brézis, *Problèmes unilatéraux*, J. Math. Pures. Appl. 9 (1972) 1-168.

3. C. Cao, E.S. Titi, *Global well-posedness of the three dimensional viscous primitive equations of large scale ocean and atmosphher dynamics*, Annals of Math. (2) 166 (2007) 245-267.

4. M. Coti Zelati, M. Frémond, R. Temam, J. Tribbia, *The equations of the atmosphere with humidity and saturation: uniqueness and physical bounds*, Physica D 264 (2013) 49-65.

5. M. Coti Zelati, A. Huang, I. Kukavica, R. Temam, M. Ziane, *The primitive equations of the atmosphere in presence of vapor saturation*, Nonlinearity (2015), http://dx.doi.org/10.1088/0951-7715/28/3/625, in press.

6. J. I. Diaz, *Mathematical analysis of some diffusive energy balance models in climatology*, in: J.I. Díaz, J.L. Lions (Eds.), Mathematics, Climate and Environment, in: Research Notes in Applied Mathematics, vol. 27, Masson, Paris, 1993, 28-56.

7. J. I. Diaz, L. Tello, *On a nonlinear parabolic problem on a Riemannian manifold without boundary arising in climatology*, Collect. Math. L (Fascicle 1) (1999) 19-51.

8. G. Duvaut, J.-L. Lions, *Inequalities in Mechanics and Physics*, translated from the French by C.W. John Springer-Verlag, Berlin-New York, 1976, 397 pp.

9. I. Ekeland, R. Temam, *Convex Analysis and Variational Problems*, Classics in Applied Mathematics, 1987, +416 pages.

10. E. Feireisl, J. Norbury, *Some existence, uniqueness and nonuniqueness theorems for solutions of parabolic equations with discontinuous nonlinearities*, Proc. Roy. Soc. Edinburgh, Sect. A 119 (1991) 1-17.

11. M. Frémond, *Phase Change in Mechanics*, Springer-Verlag, Berlin, Heidelberg, 2012.

12. R. Gianni, J. Hulshof, *The semilinear heat equation with a Heaviside source term*, European J. Appl. Math. 3 (1992) 369-379.

13. A. E. Gill, *Atmosphere-Ocean Dynamics*, International Geophysics Series, vol. 30, Academic Press, San Diego, 1982.

14. W. W. Grabowski, P. K. Smolarkiewicz, *A multiscale anelastic model for meteorological research*, Mon. Weather Rev. 130 (2002) 939-956.
15. B. Guo, D. Huang, Existence of weak solutions and trajectory attractors for the moist atmospheric equations in geophysics, J. Math. Phys. 47 (2006) 083508.

16. B. Guo, D. Huang, Existence of the universal attractor for the 3-D viscous primitive equations of large-scale moist atmosphere, J. Differential Equations 251 (2011) 457-491.

17. G. J. Haltiner, Numerical Weather Prediction, John Wiley and Sons, New York, 1971.

18. G. J. Haltiner, R. T. Williams, Numerical Prediction and Dynamic Meteorology, John Wiley and Sons, New York, 1980.

19. D. Kinderlehrer, G. Stampacchia, An Introduction to Variational Inequalities and Their Applications, Academic Press Inc., Harcourt Brace Jovanovich Publishers, New York, 1980.

20. J. L. Lions, Quelques méthodes de résolution des problèmes aux limites non linéaires, Dunod, Paris, 1969, reprinted in 2002.

21. J. Lions, R. Temam, S. Wang, New formulations of the primitive equations of atmosphere and applications, Nonlinearity 5 (2) (1992) 237-288.

22. J. Lions, R. Temam, S. Wang, Models for the coupled atmosphere and ocean, (CAO1), Comput. Mech. Adv. 1 (1) (1993) 5-119.

23. J. Pedlosky, Geophysical Fluid Dynamics, Springer-Verlag, New York, 1987.

24. M. Petcu, R. Temam, M. Ziane, Some mathematical problems in geophysical fluid dynamics, in: Computational Methods for the Atmosphere and the Oceans, in: Handbook of Numerical Analysis, Special vol. XIV, Elsevier, Amsterdam, 2008.

25. R. T. Rockafellar, Convex Analysis, Princeton University Press, Princeton, NJ, 1997.

26. R. R. Rogers, M. K. Yau, A Short Course in Cloud Physics, Pergamon Press, Oxford, New York, 1989.

27. R. Temam, Navier-Stokes Equations and Nonlinear Functional Analysis, SIAM, Philadelphia, 1995.

28. R. Temam, Navier-Stokes Equations: Theory and Numerical Analysis, AMS Chelsea Publishing, 2000, +408 pages.

29. R. Temam, J. Tribbia, Uniqueness of solutions for moist advection problems, Quarterly Journal of the Royal Meteorological Society, 140 (2014), 1315-1318.

30. R. Temam, J. Tribbia, The equations of the moist advection: A unilateral problem, Quarterly Journal of the Royal Meteorological Society, 00 (2015), 1-4. DOI: 10.1002/qj.2638.

31. R. Temam, X. Wang, Numerical approximation of a variational inequality related to the humid atmosphere, Submitted.

32. R. Temam, K. Wu, Formulation of the equations of the humid atmosphere in the context of variational inequalities, Journal of Functional Analysis, 269 (2015), 2187-2221.

33. F. Tone, D. Wirosoetisno, On the long-time stability of the implicit Euler scheme for the two-dimensional Navier-Stokes equations, SIAM. J. Numer. Anal, 44 (2006) 29-40.