Thomas Carraro

Vladislav Olkhovskiy

Successive Approximation of Nonlinear Confidence Regions (SANCR)

Keywords: Parameter estimation, Nonlinear Confidence Region, Covariance Matrix, Differential Equations

In parameter estimation problems an important issue is the approximation of the confidence region of the estimated parameters. Especially for models based on differential equations, the needed computational costs require particular attention. For this reason, in many cases only linearized confidence regions are used. However, despite the low computational cost of the linearized confidence regions, their accuracy is often limited. To combine high accuracy and low computational costs, we have developed a method that uses only successive linearizations in the vicinity of an estimator. To accelerate the process, a principal axis decomposition of the covariance matrix of the parameters is employed. A numerical example illustrates the method.

Introduction

To simplify the notation, we consider a nonlinear model f (t, θ), with θ ∈ R n and t ∈ R, which does not depend on an additional (dynamical) system. We assume that f is differentiable with respect to θ and continuous with respect to t.

We consider the approximation of a confidence region about parameter values estimated by nonlinear least squares. The parameters are estimated by using experimental data y i in some given points t i with i = 1, . . . , m. The observed values contain unknown errors e i that we assume additive, so the response variable can be modeled by

y i = f (t i , θ true ) + e i , (1)
where θ true is the unknown true value of the parameters. Therefore, the least squares estimator θ is the value that solves the following problem

θ = argmin 1 2 S(θ), (2)
thomas.carraro@iwr.uni-heidelberg.de

where S(θ) is the residual sum of squares

S(θ) = m i=1 (y i -f (t i , θ)) 2 . ( 3
)
We assume that the model is correct and that the errors are normal, independent and identically distributed (iid) random variables with zero mean and variance σ 2 , i.e. e i ∼ N (0, σ 2 ). The confidence regions are here interpreted (from the frequentistic perspective [START_REF] Meeker | Teaching about approximate confidence regions based on maximum likelihood estimation[END_REF]) as the regions in the parameter space covering the true value of the parameters θ true , in large samples, with probability approximately 1 -α.

The use of linearized confidence regions with nonlinear algebraic models has been extensively treated in literature, see for example [START_REF] Bard | Nonlinear parameter estimation[END_REF][START_REF] Bates | Relative curvature measures of nonlinearity[END_REF][START_REF] Cook | A note on parameter-effects curvature[END_REF][START_REF] Jennrich | Asymptotic properties of non-linear least squares estimators[END_REF][START_REF] Gallant | Confidence regions for the parameters of a nonlinear regression model[END_REF][START_REF] Pázman | Nonlinear statistical models[END_REF]. In particular, it has been shown that confidence regions derived for the linear case can be used in linearized form also for nonlinear models, but in many cases with limited accuracy [START_REF] Rooney | Design for model parameter uncertainty using nonlinear confidence regions[END_REF]. Furthermore, there are approximation techniques for nonlinear models that are not based on linearizations [START_REF] Hamilton | Accounting for intrinsic nonlinearity in nonlinear regression parameter inference regions[END_REF][START_REF] Bates | Nonlinear regression analysis and its applications[END_REF][START_REF] Seber | Nonlinear regression[END_REF][START_REF] Potocký | Confidence regions in nonlinear regression models[END_REF].

To simplify the exposition, in this work we consider an algebraic model, but the method can be used for more complex models. In fact, the problem to approximate nonlinear confidence regions for implicit models, i.e. models based on a system of (differential) equations has been considered from different points of view and for different kind of applications by several authors. To cite only few of them, see the work [START_REF] Rooney | Design for model parameter uncertainty using nonlinear confidence regions[END_REF] and the references therein for the design under uncertainty, [START_REF] Vugrin | Confidence region estimation techniques for nonlinear regression in groundwater flow: Three case studies[END_REF] for an application to ground water flow, [START_REF] Marsili-Libelli | Confidence regions of estimated parameters for ecological systems[END_REF] for ecological systems, and [START_REF] Oberkampf | Measures of agreement between computation and experiment: Validation metrics[END_REF] for additional examples. Newly, it has been presented a method based on second-order sensitivity for the approximation of nonlinear confidence regions applied to ODE based models [START_REF] Kostina | Second-order sensitivity analysis of parameter estimation problems[END_REF]. It has been shown that higher order sensitivities give a higher accurate approximation of the confidence regions than methods using only the first order sensitivities.

With this work we show that the approximation using only linearized confidence regions can be substantially improved by a systematic successive application of linearizations, in the following called Successive Approximation of Nonlinear Confidence Regions (SANCR) method. We show results for the case with only two model parameters. An extension to more than two parameters is technically straightforward and could be partially parallelized, but the effect of successive linearizations in more than two (parameter space) dimensions has yet to be studied in this framework.

This paper is organized as follows i) In Section 2 we report the two methods on which our approach is based; ii) In Section 3 we describe the new method; iii) In Section 4 we show a numerical realization of the SANCR method.

Linearized confidence region and likelihood ratio test

As explained above, there are several methods to approximate (nonlinear) confidence regions. Our method is based on the following two approaches [START_REF] Seber | Nonlinear regression[END_REF].

For a given estimator θ of the parameter θ, we consider:

(i) The method derived from the likelihood ratio test (LR)

-2 log L(θ)/L( θ) ≤ γ 2 , (4)
from which it follows

S(θ) -S( θ) ≤ γ 2 . ( 5
)
where L is the likelihood function and γ 2 is the confidence level. (ii) The method based on the Wald test that leads to the linearized confidence regions (CL):

(θ -θ) T Cov -1 (θ -θ) ≤ γ 2 , ( 6
)
where Cov is the estimated covariance matrix of the parameters. There are several approximations of Cov [START_REF] Rooney | Design for model parameter uncertainty using nonlinear confidence regions[END_REF], we use the one based on the Jacobian J of f :

Cov = s 2 (J T J) -1 , (7)
where

J i,j = ∂f (t i , θ) ∂θ j . (8)
The level γ 2 = χ 2 1-α,n is given by the 1 -α percentile of the chi-square distribution with n degrees of freedom in case σ 2 is known, and it is

γ 2 = s 2 nF (1-α,n,m-n) in case σ 2 is unknown, but approximated by s 2 = S( θ)/(m -n).
It has been proved [START_REF] Cox | Theoretical statistics[END_REF] that these two confidence regions are asymptotically equivalent, but far from the asymptotic behavior, i.e. in case of a small number of data, they perform differently as presented in [START_REF] Rooney | Design for model parameter uncertainty using nonlinear confidence regions[END_REF]. Additionally, our method show the limitation of linearized confidence regions based only on [START_REF] Conrad | Accelerating asymptotically exact MCMC for computationally intensive models via local approximations[END_REF].

One of the major goals in defining the confidence regions is the reduction of the costs associated to their computation. From the perspective of the computational costs, the method CL is cheap since it needs only one evaluation of the covariance matrix at the parameter value θ, while the method LR is much more expensive because it is based on the evaluation of the functional S in an adequately high number of points θ in the vicinity of θ to produce a contour. In addition, the extension of the confidence region is not known a priori. In practice, the number of function evaluations needed for the method LR is in the order of several thousands, for example in our case with two parameters we use a grid of 10 4 points for the method LR.

On the contrary, as indicated in the expression [START_REF] Cook | A note on parameter-effects curvature[END_REF], the covariance matrix can be evaluated at the cost of building the Jacobian J. Therefore, the major computational costs for the method CL are given by the computation of the derivatives of the model f with respect to the parameters. Thus, we have few computations of a linearized model for the method CL while many thousand computations of a nonlinear model are needed for the method LR. Unfortunately, the accuracy of these two methods is inversely related to their computational costs, with the CL method being much more inaccurate if the model is highly nonlinear. We remind that both methods are only asymptotically exact for linear models and their quality decreases far from the asymptotic behavior.

Therefore, a compromise between computational costs and precision is highly required for many practical applications especially in case the model is based on differential equations. To this aim we established a new method combining low computational costs and high accuracy.

Successive linearizations of nonlinear confidence regions

The SANCR method is based on the use of successive linearizations of the confidence region, starting from the estimated parameter value θ (see expression ( 2)) combined with the likelihood ratio test (5) as explained below examplarily for a model with two parameters.

The likelihood ratio test is used to check whether a point belongs or not to the approximate nonlinear confidence region. Instead of testing all points in the vicinity of θ we use an educated guess, i.e. the likelihood ratio test is performed only on few points lying on the contour of the linearized confidence regions. In fact, linearized confidence regions are ellipsoids in the parameter space and the directions of the semi-axis are defined by the eigenvectors of the covariance matrix as can be deduced by the quadratic form [START_REF] Conrad | Accelerating asymptotically exact MCMC for computationally intensive models via local approximations[END_REF]. Note that the covariance matrix has dimension n × n, where n is the number of parameters to estimate. Therefore, starting from θ we determine the directions of the principal axes and their length which is given by i = γ λ i , where λ i is the eigenvalue corresponding to the i th eigenvector. We perform the likelihood ratio test for the extreme points of the semi-axes, see points θ A , θ B , θ C , θ D in Figure 1.

Let be θ A the first point to be processed. If this point passes the test, i.e. if the following condition is fulfilled

S(θ A ) -S( θ) ≤ γ 2 ,
it is considered for the construction of the confidence region and the procedure continues along the second axis. On the contrary, if the point θ A does not pass the test, it is discarded and a new candidate in the same direction θθ A is chosen.

A new point θ A along the selected semi-axis is taken by scaling 1 by a factor α < 1 as shown in Figure 2(a). This procedure is repeated with a new likelihood ratio test and possibly a rescaling (reducing α) until a point that satisfies the test

S(θ A ) -S( θ) ≤ γ 2 θ 1 θ 1 θ C θ B θ A 2 θ D θ 2
Fig. 1: Definition of the points to perform the likelihood ratio test for two parameters. is found. Once this point, say θ new , has been found, we linearize the confidence region around this new point. To this aim we calculate the Jacobian J(θ new ) (see [START_REF] Cox | Theoretical statistics[END_REF]) and the covariance Cov(θ new ) (see [START_REF] Cook | A note on parameter-effects curvature[END_REF]).

After performing the eigendecomposition of the new covariance matrix, the principal axes might have changed direction due to the nonlinearity of the model, see Figure 2(b). Following the new principal directions, we can analogously find the next candidate points belonging to the confidence region, i.e. the points θ new,A , θ new,C and θ new,D , see Figure 2(b). The point θ new,B is not considered because it is the opposite extremal point of the same principal axis. In fact, instead of taking θ new,B , we perform the same procedure starting from θ B to approximate the confidence region in the direction θθ B . Therefore, this procedure is repeated along all principal axes considering both directions.

Stopping criterion

The search along one principal axis is stopped if the distance of the next accepted point, let's say θ new,A , to the previous one is less than a given tolerance

|θ new,A -θ new,A | < T OL, (9)
then the point θ new,A is retained to define the nonlinear confidence region, see Figure 3.

Contour approximation

The countour of the nonlinear confidence region is approximated by connecting all retained points, in our case θ new,A , θ new,C , θ new,D , θ C , θ D and θ B . These points are linearly connected as shown in Figure 3(b).

θ θ 2 θ 1 α 1 1 θ A θ A (a) θ 2 θ 1 θ θ new,A θ new,B θ new,D θ new,C θ new (b)

Numerical results

As an example the following model is considered

y = θ 1 t θ2 ,
where the parameter θ 1 and θ 2 are estimated by the nonlinear least squares method. To simulate the parameter estimation process we have applied perturbed data generated using the "true" values of the parameters according to the following model response:

y i = f (t i , θ true ) + e i , (11)
where e i is a random variable distributed as N (0, σ 2 ). The Table 1 indicates the values θ true and σ 2 used in the calculations, and the least squares estimated values θ found by minimizing S(θ) for a realization of the observations y i . Additionally, the Table 2 includes the measurement positions t i . One stopping criterion of the SANCR method is that the distance of two successive candidates is smaller than a given tolerance T OL, see [START_REF] Gallant | Confidence regions for the parameters of a nonlinear regression model[END_REF]. We have used T OL = 0.15.

To evaluate the results of our approach we compare it with a Markov Chain Monte Carlo (MCMC) method described in [START_REF] Haario | Dram: Efficient adaptive mcmc[END_REF] using the associated MCMC toolbox for Matlab. In fact, an alternative way to perform a statistical analysis of nonlinear models is the use of the Bayes's theorem [START_REF] Box | Bayesian Inference in Statistical Analysis[END_REF]. Bayesian inference is not the focus of our work, therefore we refer for example to [START_REF] Conrad | Accelerating asymptotically exact MCMC for computationally intensive models via local approximations[END_REF] for a presentation of the Bayesian approach. Since the MCMC method does not allow to easily define a stopping criterion to assure convergence, we have set to 5 • 10 6 the number of model evaluations in the MCMC code.

In Figure 4 the approximations of the confidence region using the four methods can be qualitatively compared. The blue dots (for the colors see the electronic version) are the points of the MCMC method. The cyan ellipse is the linearized confidence region of the method CL. The green curve is the confidence region approximated by the method LR and the red curve is the confidence region approximated by the SANCR method.

One can observe that the linearized confidence region CL is much smaller than the MCMC approximation and that it is not centered in it. The SANCR method is an approximation of the confidence region defined by the method LR obtained at a much lower computational cost than the method LR itself. The computational costs are reported in tables 3 and 4. The method CL is very cheap with only one evaluation of the nonlinear model and the evaluations of the sensitivities with respect to the two parameters, but its quality is not satisfactory. The SANCR method uses 59 function evaluations and 42 ellipses. The latter correspond to 84 sensitivity evaluations according to the number of two parameters. The LR and the MCMC methods have been used here with 10 4 , respectively 5 • 10 6 , model evaluations.

Fig. 2 :

 2 Fig. 2: Scaling the semi-axes (a) and linearize at the new point (b).

Fig. 3 :

 3 Fig. 3: Stopping criterion (a) and interpolation (b).

Table 1 :

 1 Parameters and variance

	σ 2		θ true	θ
	0.1 (0.725, 4) (0.7279, 3.9974)
	Table 2: Position of measurement points
	x 1	x 2	x 3	x 4	x 5	x 6
	1.309 1.471 1.490 1.565 1.611 1.680

Table 3 :

 3 Model evaluations of the four methods SANCR CL LR MCMC 59 1 10 4 5 • 10 6

Table 4 :

 4 Derivatives computations of the four methods SANCR CL LR MCMC

	84	2 0	0

Fig. 4: Confidence region approximated by the four methods.


Igor Bock
email: igor.bock@stuba.sk

Regularized optimal design problem for a viscoelastic plate vibrating against a rigid obstacle

Keywords: Viscoelastic anisotropic plate, variable thickness, rigid foundation, regularization, optimal control, optimality conditions

We deal with a regularized optimal control problem governed by a nonlinear hyperbolic initial-boundary value problem describing behaviour of a viscoelastic plate vibrating against a rigid obstacle. A variable thickness of a plate plays the role of a control variable. The original problem for the deflection is regularized in order to have the uniqueness of a solution to the state problem and only the existence of an optimal thickness but also necessary optimality conditions.

Introduction

Shape design optimization problems belong to frequently solved problems with many engineering applications. We deal here with a regularized optimal design problem for a viscoelastic anisotropic plate vibrating against a rigid foundation. A variable thickness of a plate plays the role of a control variable. The corresponding state initial-boundary value contact problem represents one of the most natural problem of mechanics not frequently solved because of the hyperbolic character of the presented evolutional variational inequality. We deal here with a plate made of short memory viscoelastic material. It characterizes constructions made of concrete for example ( [START_REF] Christensen | Theory of Viscoelasticity[END_REF]). The dynamic contact for a viscoelastic bridge in a contact with a fixed road has been solved in [START_REF] Bock | Dynamic contact problem for a bridge modeled by a viscoelastic von Kármán system[END_REF]. The similar optimal control problems for the beams in a boundary contact are investigated in [START_REF] Bock | An optimization of a Mindlin-Timoshenko beam with a dynamic boundary contact on the boundary[END_REF] and [START_REF] Bock | An optimal design with respect to a variable thickness of a viscoelastic beam in a dynamic boundary contact[END_REF] respectively.

Due to the variable thickness e and the contact between a bottom of the plate and the obstacle represented by a function Φ the equation for the movement u of the middle surface and the complementarity conditions have the form 1 2 ρe(x)u tt - 1 12 div [e 3 (x)grad u tt ] + [e 3 (x)(A ijk u t,xixj + B ijk u xixj )]

x k x = F + G, 0 ≤ G ⊥ u -1 2 e -Φ ≥ 0 in (0, T ] × Ω,
where F and G express a perpendicular force acting on the plate and an unknown contact force respectively. In order to derive not only the existence of optimal variable thickness e but also the necessary optimality conditions we regularize the contact condition using the function

ω → g δ (ω), g δ (ω) =      0 for w ≤ 0 6 δ 3 ω 3 -8 δ 4 ω 4 + 3 δ 5 ω 5 for 0 < ω < δ 1 δ ω for ω ≥ δ.
in an analogous way as in [START_REF] Bock | Regularized optimal control problem for a beam vibrating against an elastic foundation[END_REF], where the control problem for an elastic beam vibrating against an elastic foundation of Winkler's type was considered. We remark that instead of the function g δ we can use any not negative nondecreasing function g ∈ C 2 (R) of the variable ω vanishing for ω ≤ 0 and equaled to 1 δ ω for ω ≥ δ.

Solving the state problem we apply the Galerkin method in the same way as in [START_REF] Bock | An optimization of a Mindlin-Timoshenko beam with a dynamic boundary contact on the boundary[END_REF], where the rigid obstacle acting against a beam is considered or in [START_REF] Bock | Unilateral dynamic contact of viscoelastic von Kármán plates[END_REF] where the problem for a viscoelastic von Kármán plate vibrating against a rigid obstacle has been solved. The compactness method will be used in solving the minimum problem for a cost functional. We apply the approach from [START_REF] Bock | Regularized optimal control problem for a beam vibrating against an elastic foundation[END_REF] in deriving the optimality conditions.

2 Solving the state problem

Setting of the state problem

We consider an anisotropic plate short memory viscoelastic plate with the middle surface Ω ⊂ R 2 . The variable thickness of the plate is expressed by a positive function x → e(x), x ∈ Ω, the positive constant ρ is the density of the material, A ijk , B ijk are the symmetric and positively definite tensor expressing the viscoelastic and elastic properties of the material. The plate is clamped on its boundary. Let F : (0, T ] × Ω → R be a perpendicular load per a square unit acting on the plate. Let u 0 , v 0 : Ω → R be the initial displacement and velocity, and

a = 1 6ρ , a ijk = 2 ρ A ijk , b ijk = 2 ρ B ijk , f = 2F ρ
be the new mechanical and material characteristics. Then the vertical displacement u : (0, T ] × Ω → R is a solution of the following regularized hyperbolic initial-boundary value problem

e(x)u tt -a 1 12 div [e 3 grad u tt ] + [e 3 (x)(a ijk u t,xixj + b ijk u xixj )] x k x -g δ ( 1 2 e(x) + Φ(x) -u) = f (t, x), in (0, T ] × Ω.
(1)

u(t, ξ) = ∂u ∂n (t, ξ) = 0, t ∈ (0, T ], ξ ∈ ∂Ω (2) u(0, x) = u 0 (x), u t (0, x) = v 0 (x), x ∈ Ω. (3)
We introduce the Hilbert spaces

H ≡ L 2 (Ω), H k (Ω) = {y ∈ H : D α y ∈ H, |α| ≤ k}, k ∈ N
((y, z)) = Ω y xixj (x)z xixj (x) dx, y = ((y, y)) 1/2 , y, z ∈ V.
We denote by V * the dual space of linear bounded functionals over V with duality pairing F, y * = F (y), F ∈ V * , y ∈ V. It is a Banach space with a norm • * . The spaces V, H, V * form the Gelfand triple meaning the dense and compact embedings

V → → H → → V * .
We set I = (0, T ), Q = I × Ω. For a Banach space X we denote by L p (I; X) the Banach space of all functions y : I → X such that y(•) X ∈ L p (0, T ), p ≥ 1, by L ∞ (I; X) the space of essentially bounded functions with values in X, by C( Ī; X) the space of continuous functions y : Ī → X, Ī = [0, T ]. For k ∈ N we denote by C k ( Ī; X) the spaces of k-times continuously differentiable functions defined on Ī with values in X. If X is a Hilbert space we set

H k (I; X) = {v ∈ C k-1 ( Ī; X) : d k v dt k ∈ L 2 (I; X)}
the Hilbert spaces with the inner products

(u, v) H k (I,X) = I [(u, v) X + k j=1 (u j , v j ) X ] dt, k ∈ N.
We denote by ẇ, ẅ and ... w the first, the second and the third time derivative of a function w : I → X. In order to derive necessary optimality conditions in the next chapter we assume stronger regularity of data:

e ∈ E ad := e ∈ H 2 (Ω) : 0 < e min ≤ e(x) ≤ e max ∀x ∈ Ω, e 2 ≤ ê ;

u 0 ∈ V ∩ H 4 (Ω), u 0 (x) ≥ 1 2 e max + Φ(x) ∀x ∈ Ω; Φ ∈ C( Ω), Φ(ξ) ≤ 0 ∀ξ ∈ ∂Ω; , v 0 ∈ V, f ∈ H 1 (I; H). (4)
The symmetric and positively definite fourth-order tensors a ijk , b ijk fulfil

a ijk = a k ij = a jik , b ijk = b k ij = b jik , α 0 > 0, α 0 ε ij ε ij ≤ a ijk ε ij ε k ≤ α 1 ε ij ε ij ∀ {ε ij } ∈ R 2×2 sym , β 0 > 0, β 0 ε ij ε ij ≤ b ijk ε ij ε k ≤ β 1 ε ij ε ij ∀ {ε ij } ∈ R 2×2 sym , (5)
where the Einstein summation convention is employed and R 2×2 sym is the set of all second-order symmetric tensors. For e, u, y ∈ H 2 (Ω) we define bilinear forms A(e), B(e) by

A(e)(u, y) = e 3 a ijk u ij y k , B(e)(u, y) = e 3 b ijk u ij y k . Definition 1. A function u is a weak solution of the problem (1)-(3) if ü ∈ L 2 (I; H1 (Ω)), u ∈ L 2 (I; V ), there hold the identity Q e(x)üy + ae 3 (x)∇ü • ∇y + A(e; u, y) + B(e; u, y) dx dt = Q g δ ( 1 2 e(x) + Φ(x) -u) + f (t, x) y dx dt ∀ y ∈ L 2 (I; V ) (6)
and the initial conditions

u(0) = u 0 , u(0) = v 0 . (7)

Existence and uniqueness of the state problem

We verify the existence and uniqueness of a weak solution.

Theorem 1. There exists a unique solution u of the problem ( 6), [START_REF] Christensen | Theory of Viscoelasticity[END_REF]

such that u ∈ C 1 ( Ī; V ), ü ∈ L 2 (I; V ) ∩ L ∞ (I; H1 (Ω)) ∩ C 1 ( Ī; H), ... u ∈ L 2 (I; H)
and there hold the estimates

u C( Ī,V ) + u L2(I,V ) + u C( Ī, H1 (Ω)) ≤ C 0 (α 0 , α 1 , β 0 , β 1 , e min , e max , ê, u 0 , v 0 , f ), ( 8
) u C( Ī,V ) + ü L2(I,V ) + ü L∞(I, H1 (Ω)) + ... u L2(I,H) ≤ C 1 (δ, α 0 , α 1 , β 0 , β 1 , e min , e max , ê, u 0 , v 0 , f ). ( 9
)
Proof. Let {w i ∈ V ∩ H 4 (Ω); i ∈ N} be a basis of V . We introduce the Galerkin approximation u m of a solution in a form

u m (t) = m i=1 α i (t)w i , α i (t) ∈ R, i = 1, ..., m, m ∈ N, Ω e(x)ü m w i + ae 3 (x)∇ü m • ∇w i + A(e)( um , w i ) + B(e)(u m , w i ) dx = Ω g δ ( 1 2 e(x) + Φ(x) -u m ) + f (t) w i dx, i = 1, ..., m; u m (0) = u 0m , um (0) = v 0m ; u 0m → u 0 in H 4 (Ω), v 0m → v 0 in V.
A solution originally existing only locally can be prolonged to the whole time interval I with the a priori estimates

u m C( Ī,V ) + um L2(I,V ) + um C( Ī, H1 (Ω)) ≤ C 2 (α 0 , α 1 , β 0 , β 1 , e min , e max , ê, u 0 , v 0 , f ). ( 10
)
Better estimates can be achieved after differentiating the Galerkin equation with respect to t:

um C( Ī,V ) + üm L2(I,V ) + üm C( Ī, H1 (Ω)) ≤ C 3 (δ, α 0 , α 1 , β 0 , β 1 , e min , e max , ê, u 0 , v 0 , f ). ( 11
)
We proceed with the convergence of the Galerkin approximation. Applying the estimates (10), ( 11), the Aubin-Lions compact imbedding theorem [START_REF] Lions | Quelques méthodes de résolutions sep problémes aux limites non linéaires[END_REF], Sobolev imbedding theorems and the interpolation theorems in Sobolev spaces [START_REF] Eck | Unilateral Contact Problems in Mechanics. Variational Methods and Existence Theorems[END_REF] we obtain for a subsequence of {u m } (denoted again by {u m }) a function u ∈ C( Ī, V ) with u ∈ L ∞ (I, V ), ü ∈ L ∞ (I, H1 (Ω)) and the convergences

üm * ü in L ∞ (I, H1 (Ω)), üm ü in L 2 (I, V ), um * u in L ∞ (I; V ), u m → u in C( Ī; V ), u m → u in C 1 ( Ī; H 2-ε (Ω)) ∀ ε > 0, u m → u in C 1 ( Ī; C( Ω)). (12)
The convergence process (12) implies that a function u fulfils for a.e. t ∈ I Ω eüw + ae 3 (x)∇ü • ∇w + A(e; u, w) + B(e; u, w) dx

= Ω g δ ( 1 2 e(x) + Φ(x) -u) + f w dx, ∀ w ∈ V. (13)
The identity (6) follows directly after setting w ≡ y(t, •), y ∈ L 2 (I; V ) in (13).

The estimate (10) together with the convergences (12) implies the estimate [START_REF] Eck | Unilateral Contact Problems in Mechanics. Variational Methods and Existence Theorems[END_REF].

Due to the differentiability of g δ , f we obtain the third time derivative ... u ∈ L 2 (I; H) fulfilling Q ... u e(x)y -div (e 3 (x)∇y) + A(e; ü, y) + B(e; u, y) dx dt

= Q -g δ ( 1 2 e(x) + Φ(x) -u) u + ḟ (t, x) y dx dt ∀ y ∈ L 2 (I; V ). ( 14
)
The estimate ( 9) is then the consequence of (11) together with the convergences (12) and the relation ( 14). The proof of the uniqueness can be performed in a standard way using the Gronwall lemma.

Remark 1. The constant C 0 (α, β, e min , e max , ê, u 0 , v 0 , f, q) in the estimate [START_REF] Eck | Unilateral Contact Problems in Mechanics. Variational Methods and Existence Theorems[END_REF] does not depend on δ for δ ∈ (0, δ 0 ). It is possible to derive the existence of a variational solution u of the original problem with the rigid obstacle in a similar way as in [START_REF] Bock | Unilateral dynamic contact of viscoelastic von Kármán plates[END_REF], where the method of penalization was applied. We can use the limit of a subsequence of solutions {u δn }, δ n → 0+ to the problem ( 6), ( 7) for δ ≡ δ n , n ∈ N instead of the sequence of penalized solutions.

Let K be a closed convex set in L 2 (I; V ) of the form

K := {y ∈ L 2 (I; V ); ẏ ∈ L 2 (I; H1 (Ω)), y ≥ 1 2 e + Φ}. (15)
A function u ∈ K such that u ∈ L 2 (I; V ) and u(0, •) = u 0 solves the initial value problem for a nonstationary variational inequality

Q A(e; u, y -u) + B(e; u, y -u) -ae 3 ∇ u • ∇( ẏ -u) -e u( ẏ -u) dx dt + Ω ae 3 ∇ u • ∇(y -u) + e u(y -u) (T, •) dx ≥ Ω a∇v 0 • (∇y(0, •) -∇u 0 ) + v 0 (y(0, •) -u 0 ) dx + Q f (y -u) dx dt ∀y ∈ K. ( 16
)
3 Optimal control problem

The existence of an optimal thickness

We consider a cost functional J : L 2 (I; V ) × H 2 (Ω) → R fulfilling the assumption

u n u in L 2 (I; V ), e n e in H 2 (Ω) ⇒ J(u, e) ≤ lim inf n→∞ J(u n , e n ) ( 17
)
and formulate

Optimal control problem P : To find a control e * ∈ E ad such that J(u(e * ), e * ) ≤ J(u(e), e) ∀e ∈ E ad ,

where u(e) is a (unique) weak solution of the Problem (1)-(3).

Theorem 2. There exists a solution of the Optimal control problem P.

Proof. We use the weak lower semicontinuity property of the functional J and the compactness of the admissible set E ad of thicknesses in the space C( Ω). Let {e n } ⊂ E ad be a minimizing sequence for (18). i.e.

lim n→∞ J(u(e n ); e n ) = inf e∈E ad J(u(e), e). (19)
The set E ad is convex and closed and hence a weakly closed in H 2 (Ω) as the closed convex set. Then there exists a subsequence of {e n } (denoted again by {e n }) and an element e * ∈ E ad such that

e n e * in H 2 (Ω), e n → e * in C( Ω). ( 20
)
The a priori estimates (8), Sobolev imbedding theorems and the Ascoli theorem on uniform convergence on Ī imply the existence of a function

u * ∈ C( Ī; V ) such that u ∈ L ∞ (I; V ) ∩ C( Ī; H1 (Ω)), ü ∈ L ∞ (I; H1 (Ω)) and the convergences ü(e n ) * ü * in L ∞ (I; H1 (Ω)), u(e n ) * u * in L ∞ (I; V ), u(e n ) → u * in C( Ī; H1 (Ω)), u(e n ) * u * in L ∞ (I; V ), u(e n ) → u * in C( Ī; C 1 ( Ω)) (21)
for a chosen subsequence. Functions u n ≡ u(e n ) solve the initial value state problem ( 6), ( 7) for e ≡ e n . We verify that u * solves the problem ( 6), ( 7) with e ≡ e * . The previous convergences together with the Lipschitz continuity of g δ imply

e n ün e * ü * in L 2 (Q), e 3 n ün,i e 3 * ü * i in L 2 (Q), i = 1, 2; e 3 n un,ij e 3 * u * ij in L 2 (Q), e 3 n u n,ij e 3 * u * ij in L 2 (Q), i, j ∈ {1, 2}; g δ ( 1 2 (e n + Φ -u n ) g δ ( 1 2 (e * + Φ -u * ) in L 2 (Q).
Then and the proof is complete.

Necessary optimality conditions

We introduce the Banach space W = {w ∈ H 1 (I;

ω(e) = 1 2 e + Φ -u(e), y ∈ L 2 (I; V ).

In a similar way as in [START_REF] Bock | Regularized optimal control problem for a beam vibrating against an elastic foundation[END_REF] or [START_REF] Bock | Optimal control of a viscoelastic plate bending with respect to a thickness[END_REF] the following theorem about Fréchet differentiability of the mapping e → u(e) can be verified.

r(h) = o(h) i.e. lim h 2→0 r(h) h 2 = 0 ⇒ z(h) = u (e)h.
In order to derive necessary optimality conditions we assume that the cost functional

J(•, •) : L 2 (I; V ) × H 2 (Ω) → R is Fréchet differentiable.
The optimal control problem can be expressed in a form j(e * ) = min e∈E ad j(e), j(e) = J(u(e), e).

The functional j in (25) is Fréchet differentiable and its derivative in e * ∈ E ad has the form j (e * ), h = J u (u(e * ), e * ), u (e * )h + J e (u(e * ), e * ), h -2 , h ∈ H 2 (Ω) (26) with the duality pairings •, • , •, • -2 between L 2 (I; V ) * and L 2 (I; V ), (H 2 (Ω)) * and H 2 (Ω) respectively.

The optimal thickness e * ∈ E ad fulfils the variational inequality j (e * ), e -e * -2 ≥ 0 ∀e ∈ E ad .

which can be expressed in a form then it is possible after using the variational inequality (28) to obtain for sufficiently large N the uniqueness of the Optimal control e * .

Remark 3. We have mentioned in Remark 1 that there is a sequence {δ n }, δ n → 0+ such that a corresponding sequence of regularized solutions u δn of ( 6), [START_REF] Christensen | Theory of Viscoelasticity[END_REF]

  with the standard inner products (•, •), (•, •) k and the norms | • | 0 , • k , H1 (Ω) = {y ∈ H 1 (Ω) : y(ξ) = 0, ξ ∈ ∂Ω (in the sense of traces)} and V ≡ H2 (Ω) = {y ∈ H 2 (Ω) : y(ξ) = ∂y ∂n (ξ) = 0, ξ ∈ ∂Ω (in the sense of traces)} with the inner product and the norm

  V ) : ẅ ∈ L 2 (I; V * )} with a norm w W = w H 1 (I;V ) + ẅ L2(I;V * ) and operators A(e) : W → L 2 (I; V * ), B(e) : H 2 (Ω) → L 2 (I; V * ) by A(e)z, y = Q z[ey -a div (e 3 ∇y)] dx dt+ Q [A(e)( ż, y) + B(e)(z, y) + g δ (ω(e))zy] dx dt, (22) B(e)h, y = Q h ü(e)y -3a div (e 2 ∇y) dx dt+ Q h A (e)( ż, y) + B (e)(z, y) -1 2 g δ (ω(e))y dx dt,

Theorem 3 .

 3 The mapping u(•) : E ad → W is Fréchet differentiable and its derivative z ≡ z(h) = u (e)h ∈ W, h ∈ H 2 (Ω) fulfils for every e ∈ E ad uniquely the operator equationA(e)z = -B(e)h, z(0) = ż(0) = 0 (24)Proof. The existence of a solution z to the equation (24) can be verified using the standard Galerkin method and its uniqueness by the Gronwall lemma.We proceed with the differentiability of e → z(e) : Let h ∈ H 2 (Ω) with e + h ∈ E ad and r(h) = u(e + h) -u(e) -z(h) We have A(e)r(h) = A(e)[u(e + h) -u(e)] + B(e)h, and verify

JTheorem 4 .Remark 2 .

 42 u (u(e * ), e * ), u (e * )(e -e * ) + J e (u(e * ), e * ), e -e * -2 ≥ 0 ∀e ∈ E ad . (28) Applying Theorem 3 we obtain necessary optimality conditions in a form of a system with an adjoint state p: The optimal thickness e * , the corresponding state (deflection) u * ≡ u(e * ) and the adjoint state p * ≡ p(e * ) are solutions of the initial value problem Q e * u * tt y + ae 3 ∇u * tt • ∇y + A(e)(u * t , y) + B(e)(u * , y) dx dt= Q g δ ( 1 2 e + Φ -u * ) + f (t, x) y dx dt ∀ y ∈ L 2 (I; V ), u * (0) = u 0 , u * t (0) = v 0 , A(e * )p * = -J u (u * , e * ); p * (T ) = p * t (T ) = 0, B(e * )(e -e * ), p * + J e (u * , e * ), e -e * -2 ≥ 0 ∀e ∈ E ad . If the partial derivative e → J e (u(e), e) is strongly monotone i.e. J e (u(e 1 ), e 1 ) -J e (u(e 2 ), e 2 ), e 1 -e 2 2 ≥ N e 2 2 ∀e 1 , e 2 ∈ H 2 (Ω), N > 0,

  converges to a solution u of the original problem. If e * ≡ e * (δ n ) is a sequence of optimal thicknesses tending to some ẽ * ∈ E ad then it is an open question if ẽ * ∈ E ad is a solution of the corresponding Optimal control problem connected with the Problem (16). In this case there is no uniqueness of solutions and hence this Optimal control problem has the form Optimal control problem P : To find a couple {ũ * , ẽ * } ∈ U × E ad such that J(ũ * , ẽ * ) ≤ J(u, e) ∀{u, e} ∈ U × E ad , where U = {u ∈ K; u ∈ L 2 (I; V ), u(•, 0) = u 0 and(16) holds}.

This work was financially supported by the Ministry of Education of Slovak Republic under grant VEGA-1/0159/16 and by grant of Science and Technology Assistance Agency no. APVV-0246-12.


Roger Temam
email: temam@indiana.edu

Xiaoyan Wang

Approximation of the equations of the humid atmosphere with saturation

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

 and extend our former results to the case where the saturation qs evolves with time.

Introduction

The rigorous mathematical theory of the equations of humid atmosphere has been initiated in [START_REF] Lions | New formulations of the primitive equations of atmosphere and applications[END_REF][START_REF] Lions | Models for the coupled atmosphere and ocean, (CAO1)[END_REF] 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 [START_REF] Haltiner | Numerical Weather Prediction[END_REF][START_REF] Haltiner | Numerical Prediction and Dynamic Meteorology[END_REF][START_REF] Rogers | A Short Course in Cloud Physics[END_REF]. However, in the modelings in [START_REF] Haltiner | Numerical Weather Prediction[END_REF][START_REF] Haltiner | Numerical Prediction and Dynamic Meteorology[END_REF][START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF][START_REF] Rogers | A Short Course in Cloud Physics[END_REF], the saturation of vapor is not taken into account. As shown in [START_REF] Temam | Uniqueness of solutions for moist advection problems[END_REF][START_REF] Temam | Formulation of the equations of the humid atmosphere in the context of variational inequalities[END_REF], 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 [START_REF] Temam | Formulation of the equations of the humid atmosphere in the context of variational inequalities[END_REF] a new formulation of the problem in the context of the variational inequalities [START_REF] Brézis | Problèmes unilatéraux[END_REF][START_REF] Ekeland | Convex Analysis and Variational Problems[END_REF][START_REF] Kinderlehrer | An Introduction to Variational Inequalities and Their Applications[END_REF][START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF][START_REF] Rockafellar | Convex Analysis[END_REF]. 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 the changes of phases and the extreme cases for the vapor concentration respectively (see e.g., [START_REF] Feireisl | Some existence, uniqueness and nonuniqueness theorems for solutions of parabolic equations with discontinuous nonlinearities[END_REF][START_REF] Frémond | Phase Change in Mechanics[END_REF][START_REF] Gianni | The semilinear heat equation with a Heaviside source term[END_REF]). In addition, we manage to extend our recent study [START_REF] Temam | Numerical approximation of a variational inequality related to the humid atmosphere[END_REF] to the case that the saturation concentration q s evolves according to the thermodynamical laws.

In [START_REF] Temam | Numerical approximation of a variational inequality related to the humid atmosphere[END_REF], 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 ≤ q ≤ 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 ∈ K for almost every t ∈ [0, t 1 ] for the vapor concentration q = q(t, 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.

The problem

Formulation of the system

Let M = M × (p 0 , p 1 ) where M ⊂ 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 x = (x, y, p) to denote a typical point in M, and use n to denote the outward normal vector to ∂M, the boundary of M. Let K be the non-empty closed convex set in H 1 (M) defined as K = {q ∈ H 1 (M); 0 ≤ q ≤ 1, a.e.}. Given a fixed t 1 > 0, we consider the following problem:

To find T : (0, t 1 ) → H 1 (M), q : (0, t 1 ) → K and h q ∈ H(q -q s ) such that for q b ∈ K, there hold

∂ t T + A T T + v • ∇T + ω∂ p T - Rω c p p T = 1 p ω -h q ϕ(T ), (1)
∂ t q, q b -q + A q q + v • ∇q + ω∂ p q, q b -q ≥ - 1 p ω -h q F (T ), q b -q , (2)
with initial and boundary conditions to be specified. Here H is the multivalued Heaviside function such that H(τ ) = 0 for τ < 0,

H(0) = [0, 1], H(τ ) = 1 for τ > 0.
For the sake of simplicity, the velocity field of the fluid u := (v(x, t), ω(x, t)) ∈ R 3 is considered as a given data in this article. Throughout the presentation, we assume that the time-dependent velocity field u satisfies u ∈ L r (0, t 1 ; V ) ∩ L ∞ (0, t 1 ; H) for some given r ∈ (4, +∞]; ∇ = (∂ x , ∂ y ) and ∆ = ∂ 2

x + ∂ 2 y are the horizontal gradient and horizontal Laplace operators respectively. In this way, the operators A T and A q are defined as

A T = -µ 1 ∆ -ν 1 ∂ p ( gp R T ) 2 ∂ p , A q = -µ 2 ∆ -ν 2 ∂ p ( gp R T ) 2 ∂ p , (3)
where µ i , ν i , g, R, c p are positive constants and T = T (p) is the average temperature over the isobar with pressure p. We assume that T satisfies:

T * ≤ T (p) ≤ T * , |∂ p T (p)| ≤ M, for some postive constants T * , T * , M and p ∈ [p 0 , p 1 ]. ( 4
) Concerning the right hand sides of equations ( 1)-( 2), the functions F and ϕ both from R 1 to R 1 are defined as

F (ζ) = q s ζ RL(ζ) -c p R v ζ c p R v ζ 2 + q s L(ζ) 2 , with L(ζ) = c 1 -c 2 ζ; ϕ(ζ) = 1 c p L(ζ)F (ζ). (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 ϕ are globally Lipschitz; ω + := max{ω, 0} refers to the positive part of ω.

We partition the boundary of M as ∂M = Γ i ∪ Γ u ∪ Γ l with Γ i , Γ u and Γ l defined by

Γ i = {x ∈ M; p = p 1 }, Γ u = {x ∈ M; p = p 0 }, Γ l = {x ∈ M; p 0 ≤ p ≤ p 1 , (x, y) ∈ ∂M }.
We supplement the system (1)-( 2) with the boundary conditions

     ∂ p T = α(T * -T ), ∂ p q = β(q * -q) on Γ i , ∂ p T = 0, ∂ p q = 0 on Γ u , ∂ n T = 0, ∂ n q = 0 on Γ l , (6)
and initial conditions

T (x, 0) = T 0 (x), q(x, 0) = q 0 (x). ( 7
)
If we allow q s to evolve, the dependence of the nonlinear functions F and ϕ on q s should be made explicit, i.e., F = F (T, q s ), ϕ = ϕ(T, q s ). In this case, we augment (1)-( 2) with the following governing equation for the evolution of q s (see [START_REF] Haltiner | Numerical Weather Prediction[END_REF], [START_REF] Haltiner | Numerical Prediction and Dynamic Meteorology[END_REF], [START_REF] Rogers | A Short Course in Cloud Physics[END_REF]):

dq s dt = - δω - p h q F (T, q s ). ( 8
)
We will also impose a further initial condition q s (x, 0) = q s,0 (x).

We shall always assume q 0 ∈ L 2 (M), 0 ≤ q 0 ≤ 1 for a.e. x ∈ M, and q s,0 ∈ L 2 (M) ∩ L ∞ (M), 0 < q s,0 < 1 and 0 ≤ q * ≤ 1 for a.e. x ∈ M, and assume the boundary datum T * and q * to satisfy T * , q * ∈ L 2 (0, t 1 ; L 2 (Γ 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).

Functional analytic framework

We denote as usual H = L 2 (M), V = H 1 (M). We use (•, •) L 2 (regarded the same as (•, •) H ) and | • | L 2 to denote the usual scalar product and induced norm in H. In the space V , we will use ((•, •)) and • to denote the scalar product adapted to the problem under investigation

((ϕ, φ)) := (∇ϕ, ∇φ) + (∂ p ϕ, ∂ p φ) + Γi ϕφ dΓ i ,
and the induced norm. The symbol •, • 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 u:

H = {u ∈ H × H × H div u = 0 and u • n = 0 on ∂M}, V = {u ∈ V × V × V div u = 0 and u • n = 0 on ∂M}.
For T, T b , q, q b ∈ V , we have the following specific forms for the duality pairs through integration by parts and in view of the Neumann boundary conditions:

A T T, T b = µ 1 (∇T, ∇T b ) H +ν 1 M gp R T 2 ∂ p T ∂ p T b dM+ν 1 Γi gp 1 R T 2 α(T -T * )T b dΓ i , (10)
A q q, q b = µ 2 (∇q, ∇q b ) H +ν 2 M gp R T 2 ∂ p q∂ p q b dM+ν 2 Γi gp 1 R T 2 β(q-q * )q b dΓ i .
(11) Consequently, we define the following bilinear forms

a T (T, T b ) = µ 1 (∇T, ∇T b ) H +ν 1 M gp R T 2 ∂ p T ∂ p T b dM+ν 1 α Γi gp 1 R T 2 T T b dΓ i , (12)
a q (q, q b ) = µ 2 (∇q, ∇q b ) H + ν 2 M gp R T 2 ∂ p q∂ p q b dM + ν 2 β Γi gp 1 R T 2 qq b dΓ i .
(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 ).

(

) 14
As for the Navier-Stokes equation, we define

b(u, U, U b ) := M (u • ∇ x,y,p U ) • U b dM = b T (u, T, T b ) + b q (u, q, q b ), (15)
where b T (u, T, T b ) and b q (u, q, q b ) are given by

b T (u, T, T b ) = M (v•∇T +ω∂ p T )T b dM, b q (u, q, q b ) = M (v•∇q+ω∂ p q)q b dM.
(16) In view of the last term in the left hand side of the (1), we introduce the following bilinear form

d(ω, T, T b ) = M RωT T b c p p dM. (17)
Similarly, in view of the last terms in [START_REF] Feireisl | Some existence, uniqueness and nonuniqueness theorems for solutions of parabolic equations with discontinuous nonlinearities[END_REF] and [START_REF] Frémond | Phase Change in Mechanics[END_REF], we define the linear functional

l(U b ) := l T (T b ) + l q (q b ) = ν 1 α Γi gp 1 R T 2 T * T b dΓ i + ν 2 β Γi gp 1 R T 2 q * q b dΓ i . (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 → V * defined through the relations

A T u, v := a T (u, v), A q u, v := a q (u, v), ∀ u, v ∈ V, (19)
are both bounded linear operators.

Similarly, the operators B(u, U ) = B T (u, U ), B q (u, q) : V × V 2 → (V * ) 2 and D(u, u) :

H × V → V defined by B(u, U ), U b := b T (u, T, T b ), b q (u, q, q b ) , ∀ u ∈ V, U, U b ∈ V 2 , (20)
and

D(u, u), v := d(ω, u, v), ∀ u ∈ H, u, v ∈ V, (21)
are also bounded. Due to the divergence free condition of u, we easily see that for any

T, q ∈ V , b T (u, T, T ) = 0, b q (u, q, q) = 0. ( 22
)
Concerning the boundedness of the above functionals, we have the following lemma.

Lemma 1 (Boundedness of the functionals). Assume U, U b ∈ V 2 and u ∈ V. There exist universal positive constants λ and K i ,

1 ≤ i ≤ 6 such that |a T (T, T b )| ≤ K 1 T T b , a T (T, T ) ≥ λ T 2 ; ( 23
)
|a q (q, q b )| ≤ K 2 q q b , a q (q, q) ≥ λ q 2 ; ( 24)

|b(u, U, U b )| ≤ K 3 u V |U | 1 2 L 2 U 1 2 U b ; (25) |d(ω, T, T b )| ≤ K 4 |ω| L 2 |T | 1 4 L 2 T 3 4 |T b | 1 4 L 2 T b 3 4 ; ( 26
) |l T (T b )| ≤ K 5 T b , |l q (q b )| ≤ K 6 q b . ( 27
)
Definition 1. Let (T 0 , q 0 ) ∈ H × H be such that 0 ≤ q 0 ≤ 1 a.e. in M and let

t 1 > 0 be fixed. A vector U = (T, q) ∈ L 2 (0, t 1 ; V × V ) ∩ C([0, t 1 ]; H × H) with (∂ t T, ∂ t q) ∈ L 2 (0, t 1 ; V * × V *
) is a solution to the initial and boundary value problem described by (1), ( 2), ( 6) and ( 7), if for almost every t ∈ [0, t 1 ] and for every (T b , q b ) ∈ V × K, we have

∂ t T, T b +a T (T, T b )+b T (u, T, T b )-d(ω, T, T b )-l T (T b ) = ( 1 p ω -(t)h q ϕ(T ), T b ), (28)
∂ t q, q b -q +a q (q, q b -q)+b q (u, q, q b -q)-l q (q b -q) ≥ (-

1 p ω -(t)h q F (T ), q b -q), (29)
for some h q ∈ H(q -q s ) and

U 0 = (T 0 , q 0 ). ( 30
)
3 Time discretization-The Euler scheme

Time-discretization

We assume that the velocity field u is given, time-dependent and satisfies u ∈ L r (0, t 1 ; V ) ∩ L ∞ (0, t 1 ; H) for some given r ∈ (4, +∞].

Let N be an integer which will eventually go to +∞ and set ∆t := k = t 1 /N . We will define recursively a family of elements of V × 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

u m = 1 k mk (m-1)k u(t) dt, m = 1, 2, • • • , 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 ∈ V and q m ∈ K are determined by:

T m -T m-1 k , T b + a T (T m , T b ) + b T (u m , T m , T b ) -d(ω m , T m-1 , T b ) -l T (T b ) = ( 1 p [ω m ] -h Q m ϕ(T m-1 ), T b ), (31)
q m -q m-1 k , q b -q m + a q (q m , q b -q m ) + b q (u m , q m , q b -q m ) -l q (q b -q m ) ≥ (- 1 p [ω m ] -h Q m F (T m-1 ), q b -q m ), (32)
for some h Q m ∈ 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 ϕ. Notice that we have d(ω m , T m-1 , T b ) and ϕ(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 ϕ. 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 [START_REF] Temam | Formulation of the equations of the humid atmosphere in the context of variational inequalities[END_REF]. Meanwhile, we shall keep in mind that the physical constraint on the function q in our problem, q m ∈ K is not preserved during the discretization ( 31)- [START_REF] Temam | Formulation of the equations of the humid atmosphere in the context of variational inequalities[END_REF]. Finally, the problem we meet is nonlinear. The above three aspects form the main sources of difficulties for our study.

Regularization and penalization

In view of Remark 1, we proceed our investigation by way of regularization and penalization. Let ε = (ε 1 , ε 2 ) and ε i > 0 be small for i = 1, 2. For ε 2 > 0, we define as follows the regularization H ε2 of H(•) : R → [0, 1]: equal to 0 for η ≥ 0, to 1 for η ≥ ε 2 , and linear continous between 0 and ε 2 . and consider the associated regularized and penalized problem:

To find

T m ε , q m ε ∈ V such that T m ε -T m-1 ε k , T b + a T (T m ε , T b ) + b T (u m , T m ε , T b ) -d(ω m , T m-1 ε , T b ) -l T (T b ) = ( 1 p [ω m ] -H ε2 (Q m ε -q s )ϕ(T m-1 ε ), T b ), (33)
q m ε -q m-1 ε k , q b + a q (q m ε , q b ) + b q (u m , q m ε , q b ) -l q (q b ) = ( 1 ε 1 [q m ε ] -, q b ) -( 1 ε 1 [q m ε -1] + , q b ) -( 1 p [ω m ] -H ε2 (Q m ε -q s )F (T m-1 ε ), q b ), (34)
for all T b , q b ∈ V . Notice that we have two choices for

Q m ε either Q m ε = q m-1 ε or Q m ε = q 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 1 ε1 . However, we point out that we could still obtain elegant estimates which do not depend on ε 1 (and ε 2 , k) though we have a blowing up factor 1 ε1 when we pass to the limit ε → (0+, 0+). These estimates will yield that the limit functions q m of q m ε satisfy the range requirement, i.e., 0 ≤ q m ≤ 1 for m = 1, 2, • • • N and a.e. x ∈ M.

Validity of iteration

The scheme ( 33)-(34) yields elliptic system on (T m ε , q m ε ) when

(T m-1 ε , q m-1 ε
) is known. To carry out our program, the step of finding (T m ε , q m ε ) given (T m-1 ε , q 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 m ε in our scheme ( 33)-(34).

When Q m ε = q m-1 ε
, the factor H ε2 (•) is known when we proceed to obtain T m ε and q m ε once T m-1 ε and q m-1 ε are known. We can apply the Minty-Browder surjective theorem (see e.g., [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]) or Galerkin method (see e.g. [START_REF] Temam | Navier-Stokes Equations: Theory and Numerical Analysis[END_REF]) to derive the existence of (T m ε , q m ε ). On the other hand, if Q m ε = q m ε , we only can proceed by Galerkin method (see e.g., [START_REF] Temam | Navier-Stokes Equations and Nonlinear Functional Analysis[END_REF][START_REF] Temam | Navier-Stokes Equations: Theory and Numerical Analysis[END_REF]) to derive the existence of (T m ε , q m ε ), since the factor H ε2 (•) is not known when we proceed to obtain T m ε and q m ε even though T m-1 ε and q m-1 ε are known.

A priori estimate for (T m

ε , q m ε ) The a priori estimates on (T m ε , q m ε ) independent of k and ε for the regularized and penalized problem ( 33)-(34) will be crucial for the processes of passing to the limits ε → (0+, 0+) and k → 0+.

Lemma 2. We have the estimates

|U j ε | 2 L 2 ≤ C(u, U 0 , t 1 ), ∀ 1 ≤ j ≤ N, N m=1 |U m ε -U m-1 ε | 2 L 2 ≤ C(u, U 0 , t 1 ), k N m=1 U m ε 2 ≤ C(u, U 0 , t 1 ), (35)
where C(u, U 0 , t 1 ) is a finite constant depending on the given datum u, U 0 and t 1 , but independent of ε 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 ϕ is not. Due to this reason, we need the following version of the so-called discrete Gronwall lemma [START_REF] Tone | On the long-time stability of the implicit Euler scheme for the two-dimensional Navier-Stokes equations[END_REF]: Lemma 3 (Discrete Gronwall Lemma). Let θ 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, • • • , N 0 satisfy the following relation X m ≤ X m-1 (1 + θY m ) + θZ m . ( 36
)
Then for m = 1, 2, 3, • • • , N 0 , the following estimates hold

X m ≤ X 0 exp( m-1 i=0 θY i+1 ) + m-1 i=1 θZ i exp( m-1 j=i θY j+1 ) + θZ n . (37)
The iteration relation (36) in the Discrete Gronwall Lemma explains our choice of index in ϕ(T m-1 ) and d(ω m , T m-1 , T b ) in the initial discretization.

We have the following a priori bound for the norm k

N m=1 U m ε -U m-1 ε k 2
V * , which will be used in our compactness argument. Lemma 4. For any ε 1 > 0 and any ε 2 > 0, the inequality

k N m=1 U m ε -U m-1 ε k 2 V * ≤ C(u, U 0 , t 1 ) < +∞, ( 38
)
holds for some constant C(u, U 0 , t 1 ) depending on U 0 , u, t 1 , but not on ε and k.

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

Lemma 5. The following bounds hold:

k N m=1 [q m ε ] - ε 1 2 L 2 ≤ C|ω| 2 L 2 (0,t1;H) , k N m=1 [q m ε -1] + ε 1 2 L 2 ≤ C|ω| 2 L 2 (0,t1;H) . ( 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 q m -q m-1 k

, [q m ] -and q m -q m-1 k , [q m -1] + 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, • • • , 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

N m=1 q m -q m-1 k , [q m ] -≤ 0, - N m=1 q m -q m-1 k , [q m -1] + ≤ 0, (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 ε → (0+, 0+)

Assume the time step k > 0 is fixed. Our goal, in this part, is to pass to the limit ε → (0+, 0+) in the scheme ( 33)-(34). The limit functions (T m , q m ) of (T m ε , q m ε ) will be solutions to the time discretized scheme (31)- [START_REF] Temam | Formulation of the equations of the humid atmosphere in the context of variational inequalities[END_REF]. These solutions will serve as building blocks for us to construct approximate solutions to our original problem.

After extracting a finite number of subsequences, ε → 0, we infer from Lemma 2 that, for m = 1, 2, • • • , N there exist functions U m ∈ V such that, as ε → 0+

U m ε U m weakly in V. (41)
We still use ε as the index for the subsequence.

Since the inclusion V ⊂ H is compact and U m ε is weakly convergent in V , it is strongly convergent in H, i.e., we also have

U m ε → U m strongly in H. ( 42
)
By an additional extraction of subsequences we see that:

U m ε (x) → U m (x) a.e., m = 1, 2, • • • , N. (43)
Meanwhile, we have

H ε2 (Q m ε -q s ) h Q m weak-* in L ∞ (M). Concerning the limit functions q m , the second component of U m for m = 1, 2, • • • , N , we know from Lemma 5 that k N m=1 |[q m ε ] -| 2 L 2 + |[q m ε -1] + | 2 L 2 ≤ Cε 2 1 |ω| 2 L 2 (0,t1;H) . (44)
As the real functions g ± (θ) = θ ± are both Lipschitz functions with Lipschitz constant 1 on R, i.e., |g

± (θ 1 ) -g ± (θ 2 )| ≤ |θ 1 -θ 2 |, we have |[q m ε ] --[q m ] -| L 2 ≤ |q m ε -q m | L 2 , |[q m ε -1] + -[q m -1] + | L 2 ≤ |q m ε -q m | L 2 .
Consequently, with (42

) we have [q m ε ] -→ [q m ] -and [q m ε -1] + → [q m -1] + in H.
As k > 0 is a fixed number, we can pass to the limit on ε in (44) to obtain that

N m=1 |[q m ] -| 2 L 2 + |[q m -1] + | 2 L 2 = 0, which implies 0 ≤ q m ≤ 1, a.e. in x ∈ M, i.e., q m ∈ K. (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 m ε to their limits U m in H by the compact embedding V ⊂ H. The strong convergence will imply the a.e. convergence of U m ε in M up to subsequences. The strong convergence and the a.e. convergence of U 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 ∈ K that 1 ε 1 [q m ε ] -, q b -q m ε ≥ 0, - 1 ε 1 [q m ε -1] + , q b -q m ε ≥ 0. ( 46
)
Due to weak convergence of q m in V and the weak lower semi-continuity property of the norm, we also have lim sup a q (q m ε , q b -q m ε ) = lim a q (q m ε , q b ) -lim inf a q (q m ε , q m ε ) ≤ a q (q m , q b ) -a q (q m , q m ) = a q (q m , q b -q m ).

(47) All the above three terms produce correct direction of inequalities during the limit process. Therefore, we obtain the desired variational inequalities in [START_REF] Temam | Formulation of the equations of the humid atmosphere in the context of variational inequalities[END_REF]. Third, in order to show that h Q m ∈ 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 ε the existence of a solution (T m , q m ) to (31)- [START_REF] Temam | Formulation of the equations of the humid atmosphere in the context of variational inequalities[END_REF].

Convergence of the Euler Scheme

In this section, we want to prove the convergence of the solutions of the Euler scheme ( 31)- [START_REF] Temam | Formulation of the equations of the humid atmosphere in the context of variational inequalities[END_REF] to the solutions of the system ( 28)- [START_REF] Temam | The equations of the moist advection: A unilateral problem[END_REF]. 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 → +∞ or equivalently k → 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 ε, the bounds in Lemmas 2 and 4 are now valid with U m ε replaced by the limit functions U m .

Construction of approximations

For each fixed k (or N ), we associate to the elements U 0 , U 1 , U 2 , • • • , U N the following approximate functions U k = (T k , q k ), Ũk = ( Tk , qk ) and W k = (T k , Q k ) which are defined piecewise on [0, t 1 ] and take values in the space V 2 :

U k (t) = U m , Ũk (t) = U m-1 , for t ∈ [(m -1)k, mk), m = 1, 2, • • • , N. (48) W k (t) = U m -U m-1 k (t-(m-1)k)+U m-1 , for t ∈ [(m-1)k, mk), m = 1, 2, • • • , N. (49)

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 , Ũk and W k in the limit process k → 0+.

Lemma 6. For the functions U k , W k and Ũk defined above, there hold

|U k -W k | L 2 (0,t1;H) ≤ C(u, U 0 , t 1 ) √ k, |U k -Ũk | L 2 (0,t1;H) ≤ C(u, U 0 , t 1 ) √ k.
Now, we state a result concerning the boundedness of the functions U k , Ũk and W k .

Lemma 7. The functions U k , Ũk and W k remain in a bounded set of L 2 (0,

t 1 ; V )∩ L ∞ (0, t 1 ; H) as k → 0+. The functions ∂ t W k form a bounded set in L 2 (0, t 1 ; V * ) and U k -W k → 0 in L 2 (0, t 1 ; H) strongly as k → 0+.
Lemmas 6 and 7 can be regarded as reinterpretations of Lemmas 2 and 4 in terms of the functions U k , W k , Ũk .

Define u k : [0, t 1 ] → V as follows:

u k (t) = u m , for t ∈ [(m -1)k, mk), m = 1, 2, • • • , N. (50)
We have the following classical lemma:

Lemma 8 (Convergence of u k ). For the functions u k defined above, there holds

u k → u, in L r (0, t 1 , V) as k → 0 + . ( 51
)
For later use, we also define the linear averaging map for the test functions

U b = (T b , q b ) ∈ L 2 (0, t 1 ; V ) that we will use below, that is, we define U b k : [0, t 1 ] → V 2 piecewise by U b k (t) = 1 k mk (m-1)k U b (t) dt on [(m -1)k, mk).
Similarly as in Lemma 8, we conclude that

U b k → U b strongly in L 2 (0, t 1 ; V 2 ) as k → 0. Moreover, if q b ∈ K for a.e. t ∈ [0, t 1 ], we have q b k ∈ K for all t ∈ [0, t 1 ].

Passage to the limit: k → 0+

We first reinterpret as follows the scheme ( 31)- [START_REF] Temam | Formulation of the equations of the humid atmosphere in the context of variational inequalities[END_REF] in terms of the functions

U k = (T k , q k ), Ũk = ( Tk , qk ), W k = (T k , Q k ) and U b k = (T b k , q b k ): ∂ t T k , T b k + a T (T k , T b k ) + b T (u k , T k , T b k ) -d(ω k , Tk , T b k ) -l T (T b k ) = ( 1 p [ω k ] -h Q k ϕ( Tk ), T b k ), ( 52
)
∂ t Q k , q b k -q k + a q (q k , q b k -q k ) + b q (u k , q k , q b k -q k ) -l q (q b k -q k ) ≥ (- 1 p [ω k ] -h Q k F ( Tk ), q b k -q k ), (53)
where

Q k is either qk or q k . Furthermore, h Q k is defined by h Q k (t) = h Q m when t ∈ [(m -1)k, mk).
Here we emphasize that we require q b k ∈ L 2 (0, t 1 ; K). Due to Lemma 7, we have, up to subsequences, in the limit k → 0+, that

U k U = (T, q), weakly in L 2 (0, t 1 ; V ) and weak- * in L ∞ (0, t 1 ; H), (54)
W k W = (T , Q), weakly in L 2 (0, t 1 ; V ) and weak- * in L ∞ (0, t 1 ; H), (55)
and

∂ t W k ∂ t W = (∂ t T , ∂ t Q), weakly in L 2 (0, t 1 ; V * ). ( 56
)
Obviously, Ũk = U (• -k) converges also to U in L 2 (0, t 1 ; V ) weakly and in L ∞ (0, t 1 ; H) weak- * .

In view of Lemma 6, we know that

U = W. ( 57
)
Now, we consider the inclusions V ⊂ H ⊂ 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 → W, strongly in L 2 (0, t 1 ; H). ( 58
)
By Lemma 6 again, we conclude that

U k , Ũk , W k → U, strongly in L 2 (0, t 1 ; H). ( 59
)
With the above preparations, we can now pass to the limit k → 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 [START_REF] Temam | Navier-Stokes Equations: Theory and Numerical Analysis[END_REF]. 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 ∈ 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 ∈ L 2 (0, t 1 ; K). Third, the subtle point in this passage to the limit is to deal with the term

t1 0 ∂ t Q k , q b k -q k dt which is the sum of t1 0 ∂ t Q k , q b k -Q k dt and t1 0 ∂ t Q k , Q k -q k dt.
Using integration by parts, (56) and the lower semicontinuity of the norm, we write lim sup

t1 0 ∂ t Q k , q b k -Q k dt = -lim inf t1 0 ∂ t Q k , Q k dt + lim t1 0 ∂ t Q k , q b k dt = -lim inf 1 2 |Q k (t 1 )| 2 L 2 + 1 2 |q 0 | 2 L 2 + t1 0 ∂ t q, q b dt ≤ - 1 2 |q(t 1 )| 2 L 2 + 1 2 |q 0 | 2 L 2 + t1 0 ∂ t q, q b dt = - t1 0 ∂ t q, q dt + t1 0 ∂ t q, q b dt = t1 0 ∂ t q, q b -q dt. ( 60
)
where we have used, in the second equality of (60), the observation

lim t1 0 ∂ t Q k , q b k dt → t1 0 ∂ t q, q b dt, (61)
which is a simple consequence of (56) and the strong convergence of q b k to q b in L 2 (0, t 1 ; V ).

A subtle point is the treatment of

t1 0 ∂ t Q k , Q k -q k dt.
Though we have (56) (which implies in particular that ∂ t Q k is bounded in L 2 (0, t 1 ; V * ) and Q k -q k 0 weakly in L 2 (0, t 1 ; V ), we can not conclude that the limit of

t1 0 ∂ t Q k , Q k - q k dt
is 0. Rather, we show, by the specific forms of Q k and q k , that lim sup

t1 0 ∂ t Q k , Q k -q k dt ≤ 0. (62)
Indeed, noticing that

∂ t Q k = q m -q m-1 k and Q k -q k = q m -q m-1 k (t -mk) on the subinterval [(m -1)k, mk) of [0, t 1 ], we have: t1 0 ∂ t Q k , Q k -q k dt = N m=1 mt (m-1)t ∂ t Q k , Q k -q k dt = N m=1 mk (m-1)t q m -q m-1 k , q m -q m-1 k (t -mk) dt = N m=1 mk (m-1)t |q m -q m-1 | 2 L 2 k 2 (t -mk) dt ≤ 0,
which implies (62). From ( 60) and (62), we can conclude that lim sup

t1 0 ∂ t Q k , q b k -q k dt ≤ t1 0 ∂ t q, q b -q dt. (63)
To sum up, we have proved the following theorem when q s is constant:

Theorem 1. Given T 0 , q 0 ∈ H with 0 ≤ q 0 ≤ 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.

The nonlinearities ϕ and F

When q s evolves according to [START_REF] Duvaut | Inequalities in Mechanics and Physics, translated from the French[END_REF], the nonlinearities F and ϕ will also depend on the function q s . The nonlinearities (see e.g., [START_REF] Haltiner | Numerical Weather Prediction[END_REF][START_REF] Haltiner | Numerical Prediction and Dynamic Meteorology[END_REF]) ϕ and F : R × R → R are defined as follows:

F (T, q s ) = q s G(T, q s ) = q s T RL(T ) -c p R v T c p R v T 2 + q s L(T ) 2 , (64)
ϕ(T, q s ) = L(T ) c p F (T, q s ), (65)
where

L(T ) = c 1 -c 2 T, G(T, q s ) = T RL(T ) -c p R v T c p R v T 2 + q s L(T ) 2 .
(66)

In the above, c 1 , c 2 , R, c p , R v are all strictly positive constants. The additional dependence of F and ϕ on q s will bring us technical complexities during the passages to the limits. Notice that the functions F , G, and ϕ 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 ϕ r , F r and G r for ϕ, F and G.

F r (T, q s ) = q s G r (T, q s ) = q s T RL(T ) -c p R v T c p R v max (T, γ) 2 + q s L(T ) 2 , ( 67
) ϕ r (T, q s ) = L(T ) c p F r (T, q s ), (68)
G r (T, q s ) = T RL(T ) -c p R v T c p R v max (T, γ) 2 + q s L(T ) 2 , (69)
where γ > 0 is smaller than any temperature on Earth. Once arriving at q s ≥ 0, we can derive G r is nonegative by observing that T ≤ LR cpRv for any temperature on Earth. It is easy to see the rational function F r is bounded and globally Lipschitz on R × R, i.e.,

|F r (ζ 1 , ξ 1 ) -F r (ζ 2 , ξ 2 )| ≤ C(|ζ 1 -ζ 2 | + |ξ 1 -ξ 2 |), ∀ ζ 1 , ζ 2 ∈ R, ξ 1 , ξ 2 ∈ [0, ∞), ( 70
) and |F r (ζ, ξ)| ≤ C, ∀ ζ ∈ R, ξ ∈ [0, ∞). (71)
The function ϕ r is also globally Lipschitz,

|ϕ r (ζ 1 , ξ 1 ) -ϕ r (ζ 2 , ξ 2 )| ≤ C(|ζ 1 -ζ 2 | + |ξ 1 -ξ 2 |), ∀ ζ 1 , ζ 2 ∈ R, ξ 1 , ξ 2 ∈ [0, ∞). (72
) In addition, as F r (0, 0) = 0, we have ϕ r (0, 0) = 0. Hence the Lipschitz function ϕ r also satisfies |ϕ r (ζ, ξ)| ≤ C(|ζ| + |ξ|). Definition 2. Let (T 0 , q 0 , q s,0 ) ∈ H × H × H be such that 0 ≤ q 0 ≤ 1, 0 < q s,0 < 1 a.e. in M and let t 1 > 0 be fixed. A vector (T, q, q s ) ∈ L 2 (0,

t 1 ; V × V ) ∩ C([0, t 1 ]; H × H) × L ∞ (M × [0, t 1 ]) with 0 < q s < 1, (∂ t T, ∂ t q, ∂ t q s ) ∈ L 2 (0, t 1 ; V * × V * ) × L ∞ (M × [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 ∈ [0, t 1 ] and for every (T b , q b ) ∈ V × K, we have

∂ t T, T b +a T (T, T b )+b T (u, T, T b )-d(ω, T, T b )-l T (T b ) = ( 1 p ω -(t)h q ϕ(T, q s ), T b ), (73)
∂ t q, q b -q +a q (q, q b -q)+b q (u, q, q b -q)-l q (q b -q) ≥ (-1 p ω -(t)h q F (T, q s ), q b -q), (74)

dq s dt = - 1 p ω -h q F (T, q s ) for a.e. (t, x) ∈ [0, t 1 ] × M, (75)
for some h q ∈ H(q -q s ) and U 0 = (T 0 , q 0 , q s,0 ). (76)

The discretization scheme

We begin with (T 0 , q 0 , q 0 s ) := (T 0 , q 0 , q s,0 ), i.e., the given initial datum.

(77) When (T 0 , q 0 , q 0 s ), (T 1 , q 1 , q 1 s ), • • • , (T m-1 , q m-1 , q m-1 s

) are known, T m ∈ V and q m ∈ K are determined by:

T m -T m-1 k , T b + a T (T m , T b ) + b T (u m , T m , T b ) -d(ω m , T m-1 , T b ) -l T (T b ) = ( 1 p [ω m ] -h q m-1 ϕ(T m-1 , q m-1 s ), T b ), (78)
q m -q m-1 k , q b -q m + a q (q m , q b -q m ) + b q (u m , q m , q b -q m ) -l q (q b -q m ) ≥ (- 1 p [ω m ] -h q m-1 F (T m-1 , q m-1 s ), q b -q m ), (79)
and

q m s := Z m (mk), (80)
where Z m (t) is the solution to the following inital value problem

dZ m (t) dt = -1 p [ω m ] -h q m-1 F (T m-1 , Z m (t)), Z m ((m -1)k) = q m-1 s , (81)
where h q m-1 ∈ H(q m-1 -q m-1 s

).

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 ϕ(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 m s by successively solving the ordinary differential equation (81) and make evaluations at specific time points. The latter will make the q m s '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 )-equation and the q m s -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 m ε , q m ε ∈ V and Z m ε ∈ H such that

T m ε -T m-1 ε k , T b + a T (T m ε , T b ) + b T (u m , T m ε , T b ) -d(ω m , T m-1 ε , T b ) -l T (T b ) = ( 1 p [ω m ] -H ε2 (q m-1 ε -q m-1 s,ε )ϕ r (T m-1 ε , q m-1 s,ε ), T b ), (82)
q m ε -q m-1 ε k , q b + a q (q m ε , q b ) + b q (u m , q m ε , q b ) -l q (q b ) = ( 1 ε 1 [q m ε ] -, q b ) -( 1 ε 1 [q m ε -1] + , q b ) -( 1 p [ω m ] -H ε2 (q m-1 ε -q m-1 s,ε )F r (T m-1 ε , q m-1 s,ε ), q b ), (83) and q m s,ε := Z m ε (mk), (84)
where Z m ε is the solution to the following initial value problem

dZ m ε (t) dt = -1 p [ω m ] -H ε2 (q m-1 ε -q m-1 s,ε )F r (T m-1 ε , Z m ε (t)), Z m ε ((m -1)k) = q m-1 s,ε . (85)

Validity of the iteration

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

F r (T m-1 ε , Z m ε (t)) = Z m ε (t)G r (T m-1 ε , Z m ε (t)).
Therefore, we have the following integral form of (85):

q m s,ε (t) = q m-1 ε exp - t (m-1)k 1 p [ω m ] -H ε2 (q m-1 ε -q m-1 s,ε )G r (T m-1 ε , Z m ε (t)) dτ . (86)
Due to the above expression of q m s,ε (t), no blow-up can happen as we always have

|Z m ε (t)| L 2 ≤ |q m-1 s,ε | L 2 , |Z m ε (t)| L ∞ ≤ |q m-1 s,ε | L ∞ for t ≥ (m -1)k. ( 87
)
For our purpose, we just need T * ≥ mk so that q m s,ε is well-defined.

5.4 A priori estimates for (T m ε , q m ε , q m s,ε ) Integrating the first eqation in (85) from (m -1)k to mk, we find

q m s,ε = q m-1 s,ε exp - mk (m-1)k 1 p [ω m ] -H ε2 (q m-1 ε -q m-1 s,ε )G r (T m-1 ε , Z m ε (t)) dτ .
(88) As 0 < q 0 s = q s,0 < 1, we easily conclude that 0 < q m s,ε < 1 for all 1 ≤ m ≤ N by repeatedly using (88) for m = 1, 2, • • • , N . Actually, we have the following pointwise monotone relations

0 < q m s,ε (x) ≤ q m-1 s,ε (x) < 1, m = 1, 2, • • • , N. (89)
Now we aim to obtain a priori estimates on (T m ε , q m ε ) independent of k and ε for the regularized and penalized problem (82)-(83). Due to the form of ϕ r , F r and our estimate on q m s,ε , we know that

|ϕ r (T m-1 ε , q m-1 s,ε )| ≤ C(|T m-1 ε | + 1) and |F r (T m-1 ε , q m-1 s,ε )| ≤ C, (90)
where C is constant independent of ε 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 m-1 ε , q m-1 s,ε )| can be bounded by a universal constant, we still have, for the penalization terms, that Lemma 5 holds. Then estimating the duality pair

U m ε -U m-1 ε k
, U b where U b ∈ V 2 , we can derive the conclusion of Lemma 4.

5.5 Passage to the limit ε → 0+ and k → 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 m ε , q m ε ) → (T m , q m ) strongly in H. In addition, by the above estimates on q m s,ε , we know that there exist q m s ∈ H ∩ L ∞ (M) such that, up to subsequences, q m s,ε q m s weakly in H and weak-* in L ∞ (M).

During the passage to the limit: k → 0+, besides defining U k , Ũk , W k , h qk as before, we should also define q s,k :

[0, t 1 ] → L ∞ (M) ∩ L 2 (M) as follows: q s,k (t) =                    Z 1 (t), when t ∈ [0, k), Z 2 (t), when t ∈ [k, 2k), • • • , Z m (t), when t ∈ [(m -1)k, mk), • • • , Z N (t), when t ∈ [(N -1)k, N k).
(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 reinterpret our scheme in terms of the above functions U k = (T k , q k ), Ũk = ( Tk , qk ),

W k = (T k , Q k ) and U b k = (T b k , q b k ) as follows ∂ t T k , T b k + a T (T k , T b k ) + b T (u k , T k , T b k ) -d(ω k , Tk , T b k ) -l T (T b k ) = ( 1 p [ω k ] -h qk ϕ r ( Tk , qs,k ), T b k ), (93)
∂ t Q k , q b k -q k + a q (q k , q b k -q k ) + b q (u k , q k , q b k -q k ) -l q (q b k -q k ) ≥ (- 1 p [ω k ] -h qk F r ( Tk , qs,k )), q b k -q k ), (94)
dq s,k dt = - 1 p [ω k ] -h qk F r ( Tk , q s,k ). ( 95
)
The two passages to the limit processes ε → 0+ and k → 0+ are parallel to those in the case that q s is constant. However, we shall pay attention to the convergences of the terms involving ϕ r (T m-1 ε , q m-1 s,ε ), F r (T m-1 ε , q m-1 s,ε ) during the passage to the limit ε → (0+, +), and ϕ r ( Tk , qs,k ), F r ( Tk , qs,k ) during the passage to the limit k → 0+ respectively. In both of the two limit processes, we could achieve strong convergences for the sequences (T m ε ) ε and ( Tk ) k up to subsequences by compact Sobolev embedding and Aubin-Lions argument respectively. While for (q m s,ε ) ε and (q s,k ) k , we could not obtain strong convergences. Fortunately, we could obtain a.e. convergences for them for x ∈ M and (t, x) ∈ [0, t 1 ] × 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. = F (T j (t), q j s (t)), q j s (0) = q s,0 ,

where F (•, •) 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 → ∞. Then up to subsequences, we have q j s (x, t) → q s (x, t) for any t ∈ [0, l] and a.e. x ∈ Ω \ Ω 0 ,

where q s (x, t) is the solution to the initial value problem dqs(t) dt = F (T (t), q s (t)), q s (0) = q s,0 ,

and Ω 0 is a subset of 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) → T (x, t) a.e. (x, t).

Writing the ODE on q j s in the integral equation form, we have q j s (t) = q j s (0) + t 0 F (T j (x, τ ), q j s (x, τ )) dτ.

(99)

By assumption, we have

Notice that dq j s (x * ,t) dt = F (T j (x * , t), q j s (x * , t)), q j s (x * , 0) = q s (x * , 0).

As | dq j s dt | L ∞ (0,l) ≤ C, we have up to a subsequence and with x * fixed that q ji s → q s uniformly on [0, l].

Since F (T ji (x * , t)) and q ji s (x * , t) converge for a.e. t ∈ [0, l], we can pass to the limit in the following equation q ji s (x * , t) = q j s (x * , 0) + t 0 F (T ji (x * , τ ), q ji s (x * , τ )) dτ, and by the Lebesgue Dominated Convergence Theorem, we have q s (x * , t) = q s (x * , 0) + t 0 F (T (x * , τ ), q s (x * , τ )) dτ.

By uniqueness of the solution to (99) and the fact that the limit solution is independent of the subsequence, we conclude that q j s (•, x * ) → q s (•, x * ), uniformly in t, for any x * such that (5.5) holds, i.e., for x * ∈ Ω \ Ω 0 where |Ω 0 | = 0, Hence, q j s (x, t) → q s (x, t) for any t and for a.e. x ∈ Ω \ Ω 0 .

(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 ∈ [0, l]. In this case we can regard T j as a constant function of t ∈ [0, l] and consequently, we still have T j → 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 ∈ H with 0 ≤ q 0 ≤ 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).

Lemma 9 .

 9 Consider the following equation dq j s (t) dt

l 0 M

 0 |T j (x, t)-T (x, t)| 2 dxdt = M l 0 |T j (x, t)-T (x, t)| 2 dt dx → 0, as j → ∞.By Fubini's theorem, we have for a.e. x thatl 0 |T j (x, t) -T (x, t)| 2 dt → 0. Now take x * such that l 0 |T j (x * , t) -T (x * , t)| 2 dt → 0.

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.