M Pierre

R Texier-Picard

Global Existence for Degenerate Quadratic Reaction-Diffusion Systems

Keywords: reaction-diffusion systems, degenerate parabolic system, global existence, atmosphere modelling, quadratic nonlinearities 1991 MSC: MSC 35K45, 35K50, 35K57, 35K65, 86A10

We consider a class of degenerate reaction-diffusion systems with quadratic nonlinearity and diffusion only in the vertical direction. Such systems can appear in the modeling of photochemical generation and atmospheric dispersion of pollutants. The diffusion coefficients are different for all equations. We study global existence of solutions.

Nous considérons une classe de systèmes dégénérés de réaction-diffusion avec une non-linarité quadratique et avec diffusion uniquement dans la direction verticale. De tels systèmes peuvent apparaître dans la modélisation de la synthèse de polluants par réactions photochimiques et de leur dispersion atmosphérique. Les coefficients de diffusion sont différents pour chaque équation. Nous étudions l'existence globale de solutions.

Introduction

The goal of this paper is to analyze global existence in time of solutions to reaction-diffusion systems of the type considered in (4) and combined from various models in [START_REF] Seinfeld | Atmospheric Chemistry and Physics[END_REF], [START_REF] Kim | Computation accuracy and efficiency of the time splitting method in solving atmospheric transport equations[END_REF], [START_REF] Verwer | Numerical time integration for air pollution models[END_REF] and which describe, in particular, the atmospheric dispersion of ozone and other photochemically generated pollutants. Three main coupled difficulties appear in these systems, set in a three-dimensional spatial domain:

• First, the diffusion occurs only in the vertical direction: consequently, the system is 'degenerate'. • However, transport of species hold in all directions.

• Even if diffusion occurred in all directions (that is, even if the linear part of the system was strictly parabolic), global existence of solutions would not be obvious because of the structure of the nonlinear reaction terms. It is proved in (4) that, if the diffusion coefficients are all the same, then global existence of classical solutions does hold. Our goal here is to tackle the more difficult situation where these coefficients are different from each other and to give some global existence result.

Besides their interest with respect to the mentioned applications, these systems contain several theoretical questions of interest for lots of other reactiondiffusion systems. Let us describe more precisely an explicit family of these systems, and we will come back to more details on the corresponding difficulties.

Although we will consider more general systems, we focus on the atmospheric diffusion model of (4) to describe and comment the equations. If the geographical area to be considered is represented by a bounded regular open subset Ω of R 2 , the equations for atmospheric dispersion can be written in a cylindrical domain Q of the form Q = {(x, y, z, t) ∈ R 4 , (x, y) ∈ Ω, z ∈ (0, 1), t ∈ (0, T )}, where z = 0 represents the surface of the earth and z = 1 the limit of the troposphere.

Let the functions φ i , i = 1 to n, represent the molecular densities of the different species involved in the photochemical reaction. Then the reactionadvection-diffusion equations can be written in the form

∂φ i ∂t = d i ∂ 2 φ i ∂z 2 + ∇ • (ωφ i ) + f i (φ) + g i , (1.1)
where the velocity field ω = (ω 1 , ω 2 , ω 3 ) models the atmospheric current, d i is the diffusion coefficient for species i, f i (φ) are nonlinear reaction terms, representing the chemistry of the process, and the g i are source terms. Finally, the notation ∇ • u denotes the divergence of u in the spatial variables (x, y, z).

Throughout the paper we assume the following

                           ∀i = 1, ..., n ω i : Q → R is continuously differentiable f i : R n → R is continuously differentiable d i ∈ (0, +∞), g i ∈ L ∞ (Q).
(1.2)

In most references mentioned above, an incompressibility condition ∇ • ω = 0 is assumed. In this case, if the functions φ i are regular enough, (1.1) writes

∂φ i ∂t = d i ∂ 2 φ i ∂z 2 + ω 1 ∂φ i ∂x + ω 2 ∂φ i ∂y + ω 3 ∂φ i ∂z + f i (φ) + g i , (1.3)
We will not need this incompressibility assumption in our analysis. We only need to note that even if it is not satisfied, (1.1) can formally be written in the form (1.3) where the nonlinear term f i (φ) is replaced by

f i (φ) + (∇ • ω)φ i .
It can be checked easily that this operation does not affect the assumptions on the nonlinearity presented below. Hence in the following we will always consider equation (1.3) instead of (1.1).

Boundary and initial conditions have to be specified. Since the diffusion takes place in the vertical direction, boundary conditions are needed for z = 0 and z = 1. There is no horizontal diffusion, but conditions on the boundary of Ω are needed where the advection field -ω is inward (and only there). We denote by ∂ -Q the corresponding part of the boundary of Q. It is defined as follows: let ν = (ν 1 , ν 2 ) ∈ R 2 be the normal outward unitary vector at a point of the boundary ∂Ω; then

∂ -Q = {(x, y, z, t) ∈ ∂Ω × (0, 1) × (0, T ); ω 1 ν 1 + ω 2 ν 2 > 0} . (1.4)
We will choose the same boundary and initial conditions as in the references already mentioned, namely

φ i (x, y, z, 0) = φ 0 i (x, y, z) (t = 0), ∂φ i ∂z (x, y, 1, t) = 0 (z = 1), - ∂ ∂z [d i φ i ] + µ i φ i (x, y, 0, t) = e i (x, y, 0, t) (z = 0), φ i (x, y, z, t) = θ i (x, y, z, t), (x, y, z, t) ∈ ∂ -Q.                      (1.5)
The condition at the top of the cylinder is a simple noflux condition, while the condition at the surface of the earth takes into account a smooth nonnegative prescribed flux e i due to anthropogenic and biogenic emission, and a deposition flux µ i φ i . The coefficient µ i represents a positive deposition velocity.

Note that our approach could handle as well any kind of "reasonable" boundary conditions. In particular, the following could be chosen instead of (1.5) and may actually be more significant since it takes into account the transport in the flux conditions:

∂ ∂z [d i φ i ] + ω 3 φ i (x, y, 1, t) = 0 (at z = 1), - ∂ ∂z [d i φ i ] + ω 3 φ i + µ i φ i (x, y, 0, t) = e i (x, y, 0, t) (at z = 0).
The following conditions will be throughout assumed on the data:

                   ∀i = 1, ..., n, φ 0 i ∈ L ∞ (Ω × (0, 1)), e i ∈ L ∞ (Ω × (0, T )), θ i ∈ L ∞ (∂ -Q), φ 0 i ≥ 0, e i ≥ 0, µ i ∈ [0, +∞], θ i ≥ 0, ∃ ζ i ∈ [H 2 ∩ L ∞ ](Q)
satisfying the conditions of (1.5).

(1.6)

The last regularity condition of (1.6) above could be weakened. It is only needed to reduce more easily the situation to homogeneous boundary conditions, keeping in particular the L 2 -regularity of the various derivatives involved.

Now, we come to the assumptions on the structure of the nonlinearity f = (f i ) 1≤i≤n . We recall in Section 4 the specific 20×20 nonlinearity mentioned in (4) and introduced in (9), [START_REF] Verwer | Numerical time integration for air pollution models[END_REF]. We will actually consider more general nonlinearities with the three following properties.

First, we assume that f preserves positivity of the solutions, namely quasipositivity ∀i = 1, ..., n, ∀φ ∈ [0, +∞) n with φ i = 0, f i (φ) ≥ 0.

(1.7)

This assumption is quite natural in this context.

Next, as one knows, one cannot expect global existence in time without any structure assumption on the nonlinearities. In (4), a so-called "intermediately quasi-conservative" triangular structure of the system of Section 4 is exploited to prove global existence (assuming also that all the d i 's are equal). Namely, five ordered relations are satisfied by the f i 's which allow to progressively make L ∞ estimates on the φ i 's when the d i 's are all equal; these five relations are recalled in (4.1).

Here, we will only assume ONE global dissipative property on the nonlinearity f , namely

     ∃ a 1 , • • • a n ∈ R, ∃b 1 , • • • b n ∈ R, ∀i = 1, ..., n, a i > 0, ∀φ ∈ R + n , n i=1 a i f i (φ) ≤ n i=1 b i φ i . (1.8)
Obviously, this assumption is satisfied by the nonlinearity in Section 4: add up for instance the five relations of (4.1). We may notice that, when all the coefficients d i 's are equal, say d i = d, and when, for instance ∀i, b i = 0 (to simplify), then the function

W = i a i φ i satisfies ∂W ∂t -d ∂ 2 W ∂z 2 + ω • ∇W ≤ i a i g i .
From this and the corresponding boundary conditions on W , we directly deduce an L ∞ -estimate on W and, by nonnegativity of all the φ i 's, an L ∞estimate on all of them. Global existence of solutions then easily follows.

The situation is quite more delicate when the diffusion coefficients are different from each other. It is known that, even for 2 × 2 systems with good nondegenerate diffusions and one (or even two!) conservative properties of type (1.8), blow up in finite time may occur (see e.g. ( 19)). Then, two kinds of results may be obtained for non degenerate diffusions:

-either, we assume more structure on the nonlinearity: for a 2 × 2 system for instance, we may moreover assume that an L ∞ -estimate is available for the first component; then, an L ∞ -estimate may also be deduced for the second component (see [START_REF] Martin | Influence of mixed boundary conditions in some reaction-diffusion systems[END_REF]) and global existence follows. This idea may be generalized to n × n systems with a so-called "triangular structure" which allow to successively obtain L ∞ estimates on the components, see (16; 5). In these situations, global classical solutions are obtained.

-another point of view is to look for global "weak solutions" which may not be in L ∞ for all time. This is the choice made in [START_REF] Pierre | Weak solutions and supersolutions in L 1 for reaction-diffusion systems[END_REF], [START_REF] Desvillettes | About Global Existence for Quadratic Systems of Reaction-Diffusion[END_REF] and it leads to global existence results for a quite larger class of nonlinearities. It may be proved in particular that global existence of weak solutions holds for systems with the (only) structure (1.8) and for which the nonlinear growth is at most quadratic.

A recent work [START_REF] Th | Regularity Analysis for Systems of Reaction-Diffusion Equations[END_REF] analyzes the size of the possible set of singularities of these weak solutions in a specific quadratic system.

Here, we will restrict our analysis to the (at most) quadratic situation and assume (like for the system of Section 4 and like in many more reactiondiffusion models arising in chemistry, biochemistry, biology, etc) that

∃k ≥ 0, ∀φ ∈ R n , |f (φ)| ≤ k [1 + |φ| 2 ]. (1.9)
This 'quadratic assumption' is consistent with our approach here which is mainly based on an a priori L 2 -estimate. Indeed, as noticed in [START_REF] Pierre | Blowup in reaction-diffusion systems with dissipation of mass[END_REF] 21) page 1108 and ( 6)), when one-dimensional lagrangian models are considered. We do not know how much it is necessary in our context, but it appears quite naturally in the analysis (see the discussion in Section 5).

Next section is devoted to a (classical) change of variable which reduces the first order part of the differential operator to a single ∂Φ/∂t operator. The main point is that, when (2.8) holds, then the diffusion operator in the vertical variable z remains invariant in the new set of variables and the new system reduces to a family of one-dimensional parameter-dependent nondegenerate reaction-diffusion systems to which we can apply more classical approaches. Global existence of solutions is then proved in Section 3. We recall in Section 4 the specific nonlinearity mentioned in (4). Then the last Section 5 is devoted to some remarks and to some open problems that we find interesting, together with a self-contained presentation of preliminary lemmas.

Change of variables

The system (1.3) can be viewed as a parabolic evolution problem, with degeneracies in the horizontal variables x, y, where only advection takes place. Following (4), we can also gather the horizontal advection terms and the time derivative. Our strategy here is not to view (1.3) as a degenerate elliptic system, as in ( 4), but to use a change of variables to get rid of the degeneracies. In this section we introduce new coordinates and characterize the domain in the new variables, then we write the system of equations in these new variables.

The domain in the new coordinates

The change of variables introduced in this section is equivalent to a "method of characteristics", where only the horizontal wind field is taken into account.

Let us define new coordinates (ξ, η) by

     x = β(ξ, η, z, t) y = γ(ξ, η, z, t)
where β, γ are the solutions of the Cauchy problem

∂β ∂t = -ω 1 (β(ξ, η, z, t), γ(ξ, η, z, t), z, t) (2.1) ∂γ ∂t = -ω 2 (β(ξ, η, z, t), γ(ξ, η, z, t), z, t) (2.2) β(ξ, η, z, 0) = ξ, (2.3) γ(ξ, η, z, 0) = η.
(2.4)

These solutions are defined for all t ∈ (0, T ) and of class C 1 in all variables, according to the assumptions on ω 1 , ω 2 (see (1.2)). In the new variables, the relevant domain is no longer cylindrical. But, if we define

D z,t = {(ξ, η) ∈ R 2 , (β(ξ, η, z, t), γ(ξ, η, z, t)) ∈ Ω},
we can assert that, for fixed (z, t) ∈ (0, 1) × (0, T ), the transformation

Λ z,t : (ξ, η) → (x = β(ξ, η, z, t), y = γ(ξ, η, z, t))
defines a diffeomorphism from D z,t to the section Ω of the cylinder Q.

Indeed, for fixed (z, t), and for (x, y)

∈ Ω given, then (ξ, η) = [Λ z,t ] -1 (x, y) is uniquely determined by ξ = α(t), η = δ(t)
, where (α, δ) is the unique solution to the backward Cauchy problem

α s = ω 1 (α(t -s), δ(t -s), z, t -s), δ s = ω 2 (α(t -s), δ(t -s), z, t -s), α(0) = x, δ(0) = y.
We consider now the whole 4-dimensional domain

D = {(ξ, η, z, t) ∈ R 2 × (0, 1) × (0, T ), s.t.(ξ, η) ∈ D z,t }
and characterize its boundary. Obviously the faces z = 0, z = 1, t = 0 and t = T of Q are simply transported by the diffeomorphism to similar faces of D. To identify the part ∂ -Q of the boundary of Q, we introduce the distance function d(•) to ∂Q. It is regular near ∂Ω × (0, 1) × (0, T ) and satisfies ∇d = (ν 1 , ν 2 , 0, 0) on this part of the boundary. If we denote d(ξ, η, z, t) = d(β(ξ, η, z, t), γ(ξ, η, z, t), z, t), then, we have the characterization:

[(ξ, η, z, t) ∈ ∂D] ⇔ [ d(ξ, η, z, t) = 0].
Moreover, using the chain rule and equations (2.1)-(2.2), we obtain

∂ d ∂t = ∂d ∂t -ω 1 ∂d ∂x -ω 2 ∂d ∂x ,
so that, on ∂Ω × (0, 1) × (0, T ), we have

∂ d ∂t = -ω 1 ν 1 -ω 2 ν 2 . (2.5)
According to the definition (1.4) of ∂ -Q, we deduce:

Lemma 1 Define the function τ by ∀(ξ, η, z, t) ∈ D, τ (ξ, η, z, t) = inf{σ ≥ 0, ∀s ∈ [σ, t], (ξ, η, z, s) ∈ D}.
Then the following equivalence holds:

[(Λ z,t (ξ, η), z, t) ∈ ∂ -Q] ⇐⇒ [(ξ, η, z, t) ∈ ∂ -D],
where

∂ -D = {(ξ, η, z, τ (ξ, η, z, t)), (z, t) ∈ (0, 1) × (0, T ), (ξ, η) ∈ ∂D z,t }, (2.6)
Indeed, we know by (2.5) that

∂ -D = (ξ, η, z, t) ∈ ∂D; (z, t) ∈ (0, 1) × (0, T ), (ξ, η) ∈ ∂D z,t , ∂ d ∂t < 0 .
Since d is constant on ∂D, ∂ d/∂t is the fourth component of the normal outward vector to ∂D. It is negative at a point (ξ, η, z, t) of ∂D if and only if D is "above" this point, that is if t = τ (ξ, η, z, t). Whence the lemma.

The new system of equations

Now we write the system of equations (1.3-1.5) in the new variables ξ, η, z, t. Denoting

φ i (x, y, z, t) = φ i (β(ξ, η, z, t), γ(ξ, η, z, t), z, t) = Φ i (ξ, η, z, t), (2.7)
we differentiate Φ i with respect to t. Using the chain rule and equations (2.1)-(2.2), we obtain

∂Φ i ∂t (ξ, η, z, t) = ∂φ i ∂t -ω 1 ∂φ i ∂x -ω 2 ∂φ i ∂y (β(ξ, η, z, t), γ(ξ, η, z, t), z, t).
We now make the following crucial assumption T he f unctions ω 1 , ω 2 are independent of the vertical variable z. (2.8)

Then the functions β and γ do not depend either on z and we have the equalities:

∂Φ i ∂z (ξ, η, z, t) = ∂φ i ∂z (β(ξ, η, t), γ(ξ, η, t), z, t), ∂ 2 Φ i ∂z 2 (ξ, η, z, t) = ∂ 2 φ i ∂z 2 (β(ξ, η, t), γ(ξ, η, t), z, t).
With these notations we can rewrite equation (1.3) in the form

∂Φ i ∂t = d i ∂ 2 Φ i ∂z 2 + ω 3 ∂Φ i ∂z + F i (Φ) + G i , (2.9)
where the new nonlinear term is defined by F i (Φ) = f i (φ) and the function G i (ξ, η, z, t) is defined from g i as in (2.7) .

For simplicity we denote by Λ t instead of Λ z,t the diffeomorphism introduced in section 2.1 (which is now independent of z). Then the boundary conditions (1.5) can be written in the new coordinates in the form

Φ i (ξ, η, z, 0) = φ 0 i (Λ 0 (ξ, η), z) (2.10) ∂Φ i ∂z (ξ, η, 1, t) = 0 (2.11) - ∂ ∂z [d i Φ i ] + µ i Φ i (ξ, η, 0, t) = e i (Λ t (ξ, η), 0, t) (2.12) ∀(ξ, η, z, t) ∈ ∂ -D, Φ i (ξ, η, z, t) = θ i (Λ t (ξ, η), z, t), (2.13)
Note that for the problem (2.9)-(2.13), the variables ξ and η play the role of parameters. We have to deal with a usual nondegenerate parabolic problem in the variables (t, z), where the boundary conditions (2.11) and (2.12) have a usual form, namely

               ∂Φ i ∂t -d i ∂ 2 Φ i ∂z 2 -ω 3 ∂Φ i ∂z = F i (Φ) + G i , ∂Φ i ∂z = 0 for z = 1, - ∂ ∂z [d i Φ i ] + µ i Φ i = e i for z = 0. (2.14)
A particularity is that the domain D is not a cylinder in the variable t. For each (ξ, η), the evolution is to be solved in the open (possibly empty) subset

I ξ,η := {t ∈ (0, T ); (ξ, η, z, t) ∈ D}, (2.15)
which is a countable union of intervals, independent of z, each of them being of the form (τ -, τ + ) where τ -= τ (ξ, η, z, t) for some t (again τ -, τ + are independent of z). By virtue of (2.6), we can assert that the "initial" value of Φ i at the point (ξ, η, z, τ -) is prescribed as follows

     Φ i (ξ, η, z, τ -) = φ 0 i (Λ 0 (ξ, η), z) if τ -= 0 Φ i (ξ, η, z, τ -) = θ i (Λ τ -(ξ, η), z, τ -) if τ -> 0.
(2.16)

A technical remark: We will often have to make estimates of Q a(φ i ) for regular functions a : R → [0, +∞). Then, using the above change of variables, for some constant C depending only on the data (essentially on the sup norm of 10 the jacobian of Λ z,t and of its inverse), we may bound I = Q a(φ i )dx dy dz dt by

I ≤ C D a(Φ i )dξ dη dz dt = C Ω dξ dη (0,1)×I ξ,η a(Φ i )dz dt,
where

Ω = ∪{Λ z,t (Ω); (z, t) ∈ (0, 1) × (0, T )} = {(ξ, η) ∈ R 2 ; I ξ,η = ∅}.
If one has an estimate of the form

τ + τ - 1 0 a(Φ i )dt dz ≤ τ + τ - 1 0 h(ξ, η, z, t) dt dz, (2.17)
for some integrable function h and for all (ξ, η) and all intervals (τ -, τ + ) ∈ I ξ,η , then, summing up over all intervals (τ -, τ + ) and integrating in (ξ, η) leads to We also denote, for i = 1...n,

I ≤ C D h(ξ, η, z, t)dξ dη dz dt ≤ C Q h(Λ -1 z,t (x,
∀ψ ∈ L 2 (Q), L i (ψ) = Λψ -d i ∂ 2 ψ ∂z 2 and L = (L i ) 1≤i≤n .
Theorem 2 Assume that (1.2),(1.6),(1.7),(1.8),(1.9) as well as (2.8) hold.

Then, for all g ∈ L ∞ (Q) n with g ≥ 0, there exists a unique nonnegative solution to the following system

             ∀i = 1, ..., n, φ i , ∂φ i /∂z, ∂ 2 φ i /∂z 2 , Λφ i ∈ L 2 (Q), φ i ∈ L ∞ (Q) Lφ = f (φ) + g on Q ∀i = 1, .
.., n, boundary and initial conditions (1.5) hold.

(3.1)

A starting point of the analysis is the following lemma on the linear part of the system.

Lemma 3 Assume that (1.2),(1.6) and (2.8) hold. Then, for all g i ∈ L 2 (Q), there exists a unique solution to the following system

             φ i , ∂φ i /∂z, ∂ 2 φ i /∂z 2 , Λφ i ∈ L 2 (Q) L i φ i = g i on Q
boundary and initial conditions (1.5) hold.

(3.2)

Moroever, [g i ≥ 0] ⇒ [φ i ≥ 0] and [g i ∈ L ∞ (Q)] ⇒ [φ i ∈ L ∞ (Q)].
This lemma follows from classical results that may be found in ( 17), ( 1), ( 12), [START_REF] Langlais | A degenerating elliptic problem with unilateral constraints[END_REF]. However, we give in Section 5 the main steps in its proof, and we point out how assumption (2.8) has to be used even in this lemma.

In order to prepare existence results for the complete nonlinear system, we first truncate the nonlinearities. Let M > 0 and

ϕ M ∈ C ∞ (R) such that 0 ≤ ϕ M ≤ 1, ϕ M = 1 on - M 2 , M 2 , ϕ M = 0 outside [-M, M ]. We define ∀φ ∈ R n , f M (φ) = f (φ) ϕ M (|φ| 2 ).
We denote f M = (f M i ) i etc.

Lemma 4 Assume that (1.2),(1.6),(1.7) and (2.8) hold. Then, for all g ∈ L ∞ (Q) n with g ≥ 0, there exists a unique solution to the following truncated system

             ∀i = 1, ..., n, φ i , ∂φ i /∂z, ∂ 2 φ i /∂z 2 , Λφ i ∈ L 2 (Q), φ i ∈ L ∞ (Q) Lφ = f M (φ) + g on Q ∀i = 1, .
.., n, boundary and initial conditions (1.5) hold.

(3.3)

Moreover, for all i = 1, ..., n, φ i ≥ 0.

This lemma is obtained in a standard way by a fixed point argument: we indeed have a Lipschitz perturbation of a "good" linear operator. Again, we indicate the main arguments in the last section. Preservation of positivity is due to the quasi-positivity of f M . Uniqueness is standard for reaction-diffusion systems in the family of classical uniformly bounded solutions.

Remark. To prove existence of solutions on Q with the nonlinearity f itself, the game consists in proving L ∞ -estimates on the solution of the approximate system (3.3) which do not depend on M . Then choosing M large enough, a solution of (3.3) will also be a solution of (3.1).

We will now denote by C any constant depending only on the data (but not on M ).

For simplicity, from now on, we will often write ψ t , ψ z , ψ zz , ... for the derivatives of a function ψ.

Thanks to the assumption (2.8), we may use the change of coordinates of Section 2. The system (3.3) is equivalent to the system (2.14-2.16) with F (Φ) replaced by F M (Φ) = f M (φ). For each (ξ, η), we will make estimates on the set G = {(z, t) ∈ (0, 1) × (τ -, τ + )} for all time intervals (τ -, τ + ).

We first start by the key L 2 -estimate which is valid essentially only under assumption (1.8).

Proposition 5 Let φ be the solution of the truncated system (3.3) and Φ the solution of the corresponding system in the new variables. Then

∀(ξ, η), f or all interval (τ -, τ + ), Φ L 2 (G) ≤ C(τ + -τ -) 1/2 .
(3.4)

Remark. As we know, Φ is defined for t in the domain I ξ,η (see (2.15)) which is the union of the intervals (τ -, τ + ). Summing up the above estimate on all of the sub-intervals (τ -, τ + ) of I ξ,η leads to

Φ 2 L 2 [(0,1)×I ξ,η ] ≤ C T. (3.5)
We can then go back to the initial function φ and, through the computations (2.17,2.18), obtain the bound φ 2 L 2 (Q) ≤ CT .

Proof: We set W = i a i Φ i , Z = i d i a i Φ i , where the coefficients a i are defined in (1.8). Summing all the equations in the Φ i 's (see (2.14) with F M instead of F ), we have, using also Assumption (1.8):

W t -Z zz -ω 3 W z = i a i F M i (Φ) + a i G i ≤ i b i Φ i + C.
Since a i > 0 for all i, one has i b i Φ i ≤ A W for some A > 0 depending only on the a i , b i . Then, we deduce the following inequality from which the L 2 estimate will follow:

W t -Z zz -ω 3 W z ≤ A W + C. (3.6)
We set σ := Z/W so that we may write

W t -σW zz -ω 3 W z ≤ A W + C. (3.7)
A main point is that, thanks to the positivity of the Φ i , a i ,

0 < min i d i ≤ σ ≤ max i d i . (3.8)
As in [START_REF] Pierre | Blowup in reaction-diffusion systems with dissipation of mass[END_REF], (3), we will exploit this property to get an L 2 (G)-estimate on W which obviously implies the expected L 2 -estimate on Φ. We do it by duality. For this, let Θ ∈ C ∞ 0 (G), Θ ≥ 0. We introduce the nonnegative solution θ of the dual problem

             -θ t -σθ zz + (ω 3 θ) z -Aθ = Θ on G θ z = r(2 z -1)θ f or z = 0, 1 θ = 0 at t = τ + , (3.9)
where r ∈ (0, +∞) is such that

r σ -ω 3 ≥ 0 f or z = 1, -r σ -ω 3 ≤ 0 f or z = 0. (3.10)
This choice of r is possible since σ is bounded from below and |ω 3 | is bounded from above. We denote R(z) = r(z 2 -z) so that R (z) = r(2z-1). In particular,

[θ z = r(2z -1)θ] ⇔ [e -R θ] z = 0. ( 3
≥ G θ[W t -Z zz -ω 3 W z -AW -C].
This writes

0 ≥ G (W Θ -Cθ) - 1 0 [θW ](z, τ -)dz + a, (3.12)
where a denotes the integrated terms in z, namely

a = τ + τ - {[-θZ z -ω 3 θ W + θ z Z](1, t) -[-θZ z -ω 3 θ W + θ z Z](0, t)} dt,
which, according to the boundary conditions on the Φ i 's and θ, is equal to

a = τ + τ - W θ[rσ -ω 3 ](1, t) + θ[W (ω 3 + rσ) + Σ i a i (µ i Φ i -e i )](0, t)dt.
By the choice of r as in (3.10) and thanks to µ i ≥ 0, e i ∈ L ∞ , we deduce that a ≥ -C τ + τθ(0, t) dt. Together with (3.12) and the fact that

W (•, τ -) is a given L ∞ -function, this leads to G W Θ ≤ C 1 0 θ(z, τ -)dz + τ + τ - θ(0, t)dt + G θ . (3.13)
The next task is now to prove

1 0 θ(z, τ -)dz + τ + τ - θ(0, t)dt + G θ ≤ C(τ + -τ -) 1/2 Θ L 2 (G) . (3.14)
Then, the expected estimate of Proposition 5 will follow by duality. To obtain (3.14), we actually prove maximal L 2 -regularity for the parabolic operator involved in (3.9), namely

G θ 2 + (θ t ) 2 + (θ z ) 2 + (θ zz ) 2 ≤ C G Θ 2 . (3.15)
Indeed, assuming (3.15), we then obtain, by setting ρ := (τ

+ -τ -) 1/2 : G θ ≤ |G| 1/2 G θ 2 1/2 ≤ ρ C G Θ 2 1/2 . 1 0 θ(z, τ -)dz = - τ + τ - 1 0 θ t ≤ ρ G (θ t ) 2 1/2 ≤ ρ C G Θ 2 1/2 .
Then, starting from the embedding estimate

∀t ∈ (τ -, τ + ), θ(0, t) ≤ 1 0 [θ 2 + (θ z ) 2 ](z, t) dz 1/2
, we have similarly, using again (3.15)

τ + τ - θ(0, t)dt ≤ τ + τ - dt 1 0 [θ 2 + (θ z ) 2 ] dz 1/2 ≤ ρ C G Θ 2 1/2
.

Whence (3.14). It remains to prove (3.15). We set U := e -R θ where R is defined in (3.11). Obviously, it is sufficient to prove (3.15) with θ replaced by U . The function U satisfies an equation of the form

-U t -σU zz + a(z, t)U + b(z, t)U z = Θ, U (•, τ + ) = 0, U z = 0 at z = 0 and z = 1,              (3.16)
where a, b are uniformly bounded by C and Θ = e -R Θ.

Multiply the equation in U by -U zz and integrate by parts to obtain:

1 0 -∂ t 1 2 [U z ] 2 + σ[U zz ] 2 = 1 0 [-Θ + aU + bU z ]U zz . (3.17)
Since σ is bounded from below, we can use Young's inequality to control the right-hand side as follows:

1 0 [ Θ + aU + bU z ]U zz ≤ 1 0 σ 2 [U zz ] 2 + C[ Θ2 + U 2 + (U z ) 2 ]. (3.18)
Absorbing the term in (U zz ) 2 and exploiting the linear Gronwall's inequality in the term (U z ) 2 , we deduce from (3.17), (3.18) and U (τ + ) = 0, the following estimate on G t = (t, τ + ) × (0, 1):

Gt [U zz ] 2 , sup s∈(t,τ + ) 1 0 [U z ] 2 (s) ≤ C Gt Θ2 + U 2 . (3.19)
Going back to (3.16), we derive Gt (U t ) 2 ≤ C Gt Θ2 + U 2 which, together with

U (•, t) = -τ + t U t (•, s)ds, leads to 1 0 U 2 (z, t)dz ≤ C τ + t 1 0 (U t ) 2 ≤ C τ + t 1 0 Θ2 + U 2 .
Integration of this Gronwall's inequality in t

→ 1 0 U 2 (t, z)dz gives G U 2 + (U t ) 2 ≤ C G Θ2 .
This inequality coupled with (3.19) leads to (3.15) with θ replaced by U , and this ends the proof of Proposition 5. Now, we recall classical estimates on the heat operator in one space dimension (see e.g. ( 14)):

Lemma 6 Let ψ be a solution on G of                    ψ t -dψ zz -ω 3 ψ z = F, ψ z = 0 for z = 1, -dψ z + µψ = e ∈ L ∞ (0, T ) for z = 0, ψ(z, 0) = ψ 0 ∈ L ∞ (0 , 1).
(3.20)

Then, for C = C(T, d, ω 3 , µ, e, ψ 0 ∞ ), ∀r ∈ (1, 3), ψ L r (G) ≤ C[1 + F L 1 (G) ], ∀p ∈ (1, 3/2), q = 3p/(3 -2p) ψ L q (G) ≤ C[1 + F L p (G) ], ∀p > 3/2, ψ L ∞ (G) ≤ C[1 + F L p (G) ].
Remark. These estimates may be found for instance in [START_REF] Ladyzenskaya | Linear and Quasi-linear Equations of Parabolic Type[END_REF]. The only remark to be done is that the constants above do not indeed depend on G and in particular do not depend on τ + -τ -even if this difference becomes small. A simple way to see it is to notice that a solution on the interval I = (τ -, τ + ) is the restriction to I of a solution on the fixed length interval I 1 = (τ -, τ -+ T ) with a left hand-side F equal to F on I and to zero on I 1 \ I. Then, for instance, the first inequality is obtained as follows from the same inequality on the 'fixed' domain G 1 = (0, 1) × I 1 :

ψ L r (G) ≤ ψ L r (G 1 ) ≤ C[1 + F L 1 (G 1 ) ] = C[1 + F L 1 (G) ].
Lemma 7 Let φ be the solution of the truncated system (3.3) and Φ the solution of the corresponding system in the new variables. Then

∀(ξ, η), f or all interval (τ -, τ + ), Φ L ∞ (G) ≤ C. (3.21)
Proof: We start with the L 2 -estimate of Proposition 5 where we even forget the precise dependence in (τ + -τ -), that is to say, we start with:

Φ L 2 (G) ≤ C.
By the 'quadratic hypothesis' (1.9):

F M (Φ) L 1 (G) ≤ C.
By the first statement of Lemma 6: ∀r ∈ (2, 3), Φ L r (G) ≤ C.

By (1.9) again:

∀p ∈ (1, 3/2), F M (Φ) L p (G) ≤ C.
By the second statement of Lemma 6: ∀q < +∞, Φ L q (G) ≤ C.

By (1.9): ∀m < +∞, F (Φ) L m (G) ≤ C.
Finally, by the last statement of Lemma 6:

Φ L ∞ (G) ≤ C. Lemma 7 follows.
Proof of Theorem 2: As explained in the Remark following Lemma 4, we apply the uniform estimate of Lemma 7 to the solution of the truncated system (3.3). We choose M larger than the constant C of Lemma 7. Then, the solution of the truncated system (3.3) is also a solution of the true system (3.1).

Uniqueness is classical in the class of uniformly bounded solutions.

Remark. In the context of air pollution models, the coefficients d i are the so-called "eddy diffusivities" introduced to model turbulent mixing, see for example [START_REF] Seinfeld | Atmospheric Chemistry and Physics[END_REF] and the references therein. In realistic models, these coefficients may depend on the vertical variable z. In this case, the diffusion part must be written in the divergence form,

L i ψ = Λψ - ∂ ∂z d i ∂ψ ∂z .
Adapting the proof of Proposition 5, it is easy to show that the results of Theorem 2 remain valid if the coefficients d i are differentiable functions of z, as long as the following conditions are satisfied:

∃α > 0, ∃C > 0, ∀i, ∀z ∈ (0, 1), d i (z) > α, |d i (z)| ≤ C. (3.22)
However, this approach fails as such when the diffusion coefficients vanish at z = 0, which is the case in the models stemming from the Monin-Obukhov similarity theory, see [START_REF] Seinfeld | Atmospheric Chemistry and Physics[END_REF] p. 869. Extensions of this method may be made to handle coefficients vanishing inside the domain (see [START_REF] Desvillettes | About Global Existence for Quadratic Systems of Reaction-Diffusion[END_REF]). Degeneracy at z = 0 would require further improvements.

A specific nonlinearity f

Here, as a main example, we reproduce the 20 × 20 system mentioned in [START_REF] Fitzgibbon | A degenerate reactiondiffusion system modeling atmospheric dispersion of pollutants[END_REF] and introduced in the context of pollutants models in [START_REF] Kim | Computation accuracy and efficiency of the time splitting method in solving atmospheric transport equations[END_REF] and [START_REF] Verwer | Numerical time integration for air pollution models[END_REF].

The nonlinear part is given as follows.

f 1 (φ) = -k 1 φ 1 + k 22 φ 19 + k 25 φ 20 + k 11 φ 13 + k 9 φ 11 φ 2 + k 3 φ 5 φ 2 +k 2 φ 2 φ 4 -k 23 φ 1 φ 4 -k 14 φ 1 φ 6 + k 12 φ 10 φ 2 -k 10 φ 11 φ 1 -k 24 φ 19 φ 1 ,
and

f 2 (φ) = k 1 φ 1 + k 21 φ 19 -k 9 φ 11 φ 2 -k 3 φ 5 φ 2 -k 2 φ 2 φ 4 -k 12 10 φ 2 f 3 (φ) = k 1 φ 1 + k 17 φ 4 + k 19 φ 16 + k 22 φ 19 -k 15 φ 3 f 4 (φ) = -k 17 φ 4 + k 15 φ 3 -k 16 φ 4 -k 2 φ 2 φ 4 -k 23 φ 1 φ 4 f 5 (φ) = 2k 4 φ 7 + k 7 φ 9 + k 13 φ 14 + k 6 φ 7 φ 6 -k 3 φ 5 φ 2 + k 20 φ 17 φ 6 f 6 (φ) = 2k 18 φ 16 -k 8 φ 9 φ 6 -k 6 φ 7 φ 6 + k 3 φ 5 φ 2 -k 20 φ 17 φ 6 -k 14 φ 1 φ 6 f 7 (φ) = -k 4 φ 7 -k 5 φ 7 + k 13 φ 14 -k 6 φ 7 φ 6 f 8 (φ) = k 4 φ 7 + k 5 φ 7 + k 7 φ 9 + k 6 φ 7 φ 6 f 9 (φ) = -k 7 φ 9 -k 8 φ 9 φ 6 f 10 (φ) = k 7 φ 9 + k 9 φ 11 φ 2 -k 12 φ 10 φ 2 f 11 (φ) = k 11 φ 13 -k 9 φ 11 φ 2 + k 8 φ 9 φ 6 -k 10 φ 11 φ 1 f 12 (φ) = k 9 φ 11 φ 2 f 13 (φ) = -k 11 φ 13 + k 10 φ 11 φ 1 f 14 (φ) = -k 13 φ 14 + k 12 φ 10 φ 2 f 15 (φ) = k 14 φ 1 φ 6 f 16 (φ) = -k 19 φ 16 -k 18 φ 16 + k 16 φ 4 f 17 (φ) = -k 20 φ 17 φ 6 f 18 (φ) = k 20 φ 17 φ 6 f 19 (φ) = -k 21 φ 19 -k 22 φ 19 + k 25 φ 20 + k 23 φ 1 φ 4 -k 24 φ 19 φ 1 f 20 (φ) = -k 25 φ 20 + k 24 φ 19 φ 1 .
As noticed and strongly used in (4), this nonlinearity satisfies the five following relations.

∂Φ i ∂t -d i ∂ 2 Φ i ∂z 2 -ω 3 ∂Φ i ∂z = G i , ∂Φ i ∂z = 0 for z = 1, - ∂ ∂z [d i Φ i ] + µ i Φ i = 0 for z = 0, Φ i (ξ, η, z, τ -) = 0.                              (5.1)
This problem is nondegenerate and it is well-known (see e.g. ( 14)) that for G i given in L 2 (G), it has a (unique) solution with derivatives in L 2 and the corresponding estimates, namely

τ + τ - 1 0 |(Φ i ) t | 2 + |(Φ i ) zz | 2 + |(Φ i ) z | 2 dt dz ≤ C τ + τ - 1 0 |G i | 2 dt dz, (5.2)
where C depends only on the data. The function G i is measurable with respect to (ξ, η). The above estimate preserves this measurability for the solution Φ: indeed, if G i depends continuously on (ξ, η), then so does the solution by (5.2). Now, if G is in L 2 in the four variables, it may be approximated in L 2 (G) and a.e in (ξ, η) by regular functions. Again, (5.2) ensures the convergence and the measurability of the limit. Moreover, according to the computations (2.17, 2.18) and to those of Section 2, we also have

Q |Λφ i | 2 + |(φ i ) z | 2 + |(φ i ) zz | 2 ≤ C.
For the last statement of Lemma 3, we may again use the change of variable of Section 2. We are back to a system of type (5.1) with

0 ≤ G i ≤ G i ∞ < +∞
and nonnegative bounded boundary data instead of zero data. It is well-known that they provide nonnegative bounded solutions.

Proof of Lemma 4. As in the previous proof, we use the change of variable of Section 2 and we are lead to the system

∂Φ i ∂t -d i ∂ 2 Φ i ∂z 2 -ω 3 ∂Φ i ∂z = F M i (Φ) + G i , ∂Φ i ∂z = 0 for z = 1, - ∂ ∂z [d i Φ i ] + µ i Φ i = e i for z = 0, Φ i (ξ, η, z, τ -) = Φ 0i (ξ, η, z).                              (5.3)
We have here a globally Lipschitz perturbation of the previous linear system. Since all data are in L ∞ , by a classical fixed point Theorem, we obtain unique uniformly bounded solutions on (0, T ). Measurability with respect to (ξ, η) is preserved in this approach.

To prove that the positivity is preserved when the data are nonnegative, instead of solving directly with the right-hand side F M , we first do it with F M replaced by F M • Π where Π : R n → (R + ) n is the orthogonal projection onto the positive cone (R + ) n of R n . The assumption (1.7) implies that, for all i and all Φ

∈ R n , [F M i Π(Φ) ]Φ - i ≥ 0.
Then, multipliying the i-th equation of the modified system by Φ - i and using boundary conditions lead to

- ∂ ∂t 1 2 1 0 (Φ - i ) 2 -d i 1 0 (Φ - i ) z 2 - ω 3 2 (Φ - i ) 2 z ≥ 0.
Integrating this inequality in t proves that Φ - i = 0, whence the positivity of the solutions of the modified system. But, since F M = F M • Π on the positive cone, Φ is also solution of the initial system. About the L 2 -estimate: The main ingredient in the proof of Theorem 2 is the L 2 -estimate given in Proposition 5. We could ask whether such an estimate would be true without the assumption (2.8). Let us make the question more precise. If, as in the proof of Proposition 5, we set W = i a i φ i , Z = i d i a i φ i , then summing the n equations in (1.3), we obtain similarly to (3.12):

W t -[σW ] zz + ω • ∇W -A W ≤ C, (5.4)
where 0 < min

i d i ≤ σ = Z W ≤ max i d i < +∞.
We do not know whether the inequality (5.4) implies an L 2 -estimate on W as in Proposition 5. The dual problem is

-θ t -σθ zz -∇ • (ω θ) -Aθ = Θ. ( 5.5)
The goal would be to prove L 2 -estimates for this equation (see (3.15)). Technically, when we multiply it by -θ zz as in the proof of Proposition 5, then it is not clear what to do with the new terms of the form ω 1 θ x θ zz . After integration by parts, they become θ z [ω 1z θ x + ω 1 θ xz ]. The last term can be treated easily, but it is not so clear for the first one... except if ω 1z = 0, whence the hypothesis (2.8).

We would like to emphasize that it is not even obvious to obtain L 2 -estimates for the main linear operator involved in the system (1.3). More precisely, let

θ t -θ zz -ω • ∇θ = Θ, (5.6)
with 'good boundary conditions' and Θ ∈ L 2 (Q). Does this imply "maximal regularity", namely

θ t -ω • ∇θ, θ zz ∈ L 2 (Q)? (5.7)
The above approach (multiplying by θ zz ) fails in general in exactly the same way as for the more general operator with the σ-term (and as just explained).

Therefore, here is a first open problem: for which ω does (5.6) imply (5.7) when Θ ∈ L 2 (Q)? One can find partial results in the literature from the theory of so-called ultraparabolic operators (see e.g. ( 2)). Strangely enough, it follows from this theory that, for instance, for ω(x, y, z, t) = (z, 0, 0), the answer would be positive, although ω 1z ≡ 0... Indeed, for each fixed y, the operator is then hypoelliptic in the variables (x, z, t) as shown by Kolmogorov in [START_REF] Kolmogorov | Zufällige Bewegungen (zur Theorie der Brownschen Bewegung)[END_REF] and generalized by Hörmander in [START_REF] Hörmander | Hypoelliptic second order differential equations[END_REF]. As proved in [START_REF] Rothschild | Hypoelliptic differential operators on nilpotent groups[END_REF] (see also [START_REF] Lanconelli | Linear and nonlinear ultraparabolic equations of Kolmogorov type arising in diffusion theory and in finance[END_REF]), this implies the expected L 2 (and even L p ) regularity in (x, z, t) and, therefore also in (x, y, z, t) after integration in y (see the technical remark (2.17,2.18)). This approach may be generalized to operators of the form (5.6) when they are hypoelliptic (see [START_REF] Rothschild | Hypoelliptic differential operators on nilpotent groups[END_REF]) or, equivalently, when they satisfy Hörmander's dimension condition for the space generated by the commutators of ∂ t -ω • ∇ and ∂ z up to an arbitrary order. However, this approach would not apply for instance to ω = (z, x-tz, 0) in which case the operator is not hypoelliptic (see [START_REF] Hörmander | Hypoelliptic second order differential equations[END_REF], [START_REF] Lanconelli | Linear and nonlinear ultraparabolic equations of Kolmogorov type arising in diffusion theory and in finance[END_REF]). Note that it is interesting to understand the optimal assumptions needed on ω in order to better handle nonlinear models where ω = ω(θ). Many such nonlinear models with degenerate diffusion are mentioned in Section 3 of [START_REF] Lanconelli | Linear and nonlinear ultraparabolic equations of Kolmogorov type arising in diffusion theory and in finance[END_REF].

Once the previous problem is understood, a second open problem would be to understand for which ω the solutions of (5.5) satisfy

Θ ∈ L 2 (Q) ⇒ θ, θ zz , θ t + ∇ • (ωθ) ∈ L 2 (Q).
This would be a first step in solving our main problem without Assumption (2.8).

A third open problem is to understand what happens in dimensions greater than 3. A main point in our approach is that we reduced the problem to a one-dimensional one. Assume more generally that we have a problem in dimension d > 3 with no diffusion in a number d of the directions. Then, with an assumption like (2.8), we may reduce the problem to a (d -d )-dimensional problem. We still keep the L 2 -estimate of Proposition ( 5), but, if d -d > 1, we cannot any more reach L ∞ -estimates as for d -d = 1. It is very likely that we can then obtain global "weak"-solution in the spirit of ( 18), ( 3), but this needs to be done.

Another interesting fourth problem is when the nonlinearities are superquadratic. Several models are mentioned in [START_REF] Seinfeld | Atmospheric Chemistry and Physics[END_REF] (see pages 143, 167, 209) which come with this structure. But the question is interesting for itself. Even when the diffusions are not degenerate, the problem is not yet completely understood. Some partial results are however known. For instance, if the right hand-side presents some "triangular structure", L ∞ -estimates may successively be obtained on all the components (see [START_REF] Martin | Influence of mixed boundary conditions in some reaction-diffusion systems[END_REF], ( 16)). This approach is based on L p -a priori estimates and it would be necessary to first extend them to a degenerate situation. It is probably necessary to first understand it in the case p = 2 (as in the previous problems just mentioned).

Another approach would be to look for global weak solutions. If the nonlinearity allows an a priori L 1 -estimate on the right-hand side, we may expect global weak-solutions as in [START_REF] Pierre | Weak solutions and supersolutions in L 1 for reaction-diffusion systems[END_REF]. But this needs to be done.

One more open problem would be interesting to look at, with respect to applications. The models considered here are approximate 'simplified' versions of more elaborate ones: for instance, saying that no diffusion occurs in the horizontal directions is just an approximation of the fact that diffusion does occur, but is small. Therefore, a natural question would be to study what happens exactly for the system with positive diffusion coefficients tending to zero in the horizontal directions. Since very different diffusions coefficients generate sharp behaviours in these systems, this question is not easy, but relevant.

  y), z, t)dx dy dz dt. (2.18) 3 Existence of global classical solutions In this section, we state and prove our main existence result of global classical solutions. For ψ ∈ L 2 (Q), we denote Λψ = ∂ψ ∂t -∇ • (ψ ω) + ψ ∇ • ω, computed in the sense of distributions. Note that, we formally have Λψ = ∂ψ ∂t -ω • ∇ψ, this being true only when ψ is regular enough.

. 11 )

 11 Multiplying(3.6) by θ ≥ 0 and using Z = σW lead to 0

  (Q) if the diffusions are nondegenerate. Here, the situation is more delicate since diffusion occurs only in one direction. We are not able to prove this L 2 -estimate in general. The considerations made in Section 5 indicate that there may be none in general. Strangely enough, if we assume that ω 1 , ω 2 do not depend on the vertical component z (see(2.8)), then this L 2 -estimate does hold. It follows that the nonlinear terms are bounded in L

	, (3),
	assumption (1.8) together with nonnegativity ensure that the solutions are
	bounded in L 2

1

. Since we are able to reduce the problem to a one-dimensional situation, it is then possible to bootstrap the estimates and reach L ∞ . Thus, global existence of classical solutions based on an a priori L 2 -estimate is the main contribution of this paper. Global existence in time is obtained no matter the various values of the d i 's in (0, +∞). Note that assumption (2.8) is discussed in some situations (see e.g. (

Summing up these five relations yields the assumption (1.8). Let us emphasize that we only use this last global dissipative relation, and not the five ones separately. On the other hand, we exploit the quadratic structure (see next Section for more comments on other kinds of systems).

Comments and open problems

In this section, we first give direct and elementary proofs of the preliminary Lemmas 3 and 4. Then we make some comments and indicate open problems.

Proof of Lemma 3. Preliminary remark: One way to prove this linear lemma is to use the results in (1), ( 12), [START_REF] Langlais | A degenerating elliptic problem with unilateral constraints[END_REF]. However, in these references, solutions are a priori less regular than stated in the lemma. Roughly speaking, Λφ i , ∂ 2 φ i /∂z 2 are found in an H -1 -type space rather than in L 2 . We do not know how to prove the L 2 regularity in general (see open problems below). On the other hand, under the assumption (2.8), we are able to do it. For simplicity, we take advantage of the change of variable introduced in Section 2 and we provide a more selfcontained proof. (1.6). Then, setting by translation ψ i := φ i -ζ i , in order to prove existence of φ i as in Lemma 3, we may assume all the data equal to zero except for the right-hand side g i ∈ L 2 (Q). Now by Assumption (2.8) and thanks to the change of variables described in Section 2, we are reduced to solving the new following system on each interval [τ -, τ + ):