Emmanuelle Augeraud-Véron
email: emmanuelle.augeraud@univ-lr.fr

Corresponding author Catherine Choquet
email: catherine.choquet@univ-lr.fr

Eloïse Comte
email: eloise.comte@univ-lr.fr

Éloïse Comte

Bernard Dacorogna

Optimal Control for a Groundwater Pollution Ruled by a Convection-Diffusion-Reaction Problem

Keywords: Optimal control problem, Hydrogeological state equations, Nonlinearly coupled problem, Parabolic and elliptic PDEs, Global existence, Fixed point theorem Mathematics Subject

come

Optimal control for a groundwater pollution ruled by a convection-diffusion-reaction problem

Introduction

According to the United States Environmental Protection Agency, agriculture is the leading source of water quality impacts on groundwater. In Europe, it contributes generally 50-80% of the total nitrogen and phosphorus loading to fresh water (Lankoski and Ollikainen [START_REF] Lankoski | Innovations in nonpoint source pollution policy-European perspectives, The magazine of food, farm, and resource[END_REF]). Xepapadeas [START_REF] Xepapadeas | The Economics of Non-Point-Source Pollution[END_REF] underlines that non point source pollution due to agriculture runoff is the major source of eutrophication and hypoxia. Despite increasing restrictions on fertilizer use (in Europe for example: Nitrate Directive (1991), Groundwater directive (2006), ...), nitrate and phosphorus remain present at high levels in intensive agricultural areas.

Optimal control of non point source pollution impacting groundwater has been studied for a long time ignoring the spatial dimension of the problem. Actually if, as in the seminal paper of Van der Ploeg and de Zeeuw [START_REF] Van Der Ploeg | A differential game of international pollution control[END_REF], pollution emission and water production wells are assumed at the same place, damages are instantaneous. However, empirical works have shown that, as groundwater flows are slow, a time delay between emission and polluted water intakes has to be taken into account (Bordave et al. [START_REF] Bordenave | Décalages temporels entre des modifications des pratiques agricoles et la diminution des nitrates dans les eaux superficielles[END_REF]). This time delay is mainly included through black-box models. Using delay equations, Winkler [START_REF] Winkler | Optimal control of pollutants with delayed stock accumulation[END_REF] proves that optimal pollution emission paths are monotonous if the objective is separable into an emission component and a stock component, but that they may display damping oscillations be considered on an infinite time length to take into account the long run consequences of the optimal policy. For the mathematical analysis, it means that the existence of a global in time solution has to be proved.

The paper is organized as follows: In Section 2, we present the hydrogeological state equations for a given load of fertilizer, thus modeling the transport and the displacement of the pollutant in groundwater. We consider a wide variety of boundary conditions. We next formulate our controlled problem, which presents the trade-off between benefits, that are due to fertilizer use to enhance agricultural production, and the cost due to polluted water treatment. Section 3 is devoted to the mathematical analysis of the state problem for a given source of fertilizer. An existence and uniqueness result of a global weak solution is proved, with moreover a physically realistic maximum principle. In Section 4, we present the main result of existence of an optimal control for both finite and infinite time horizons. Section 5 is devoted to the study of the asymptotic controllability of the state problem.

Description of the Model and Mathematical Assumptions

Transport and Displacement of an Exogenous Load of Fertilizer in the Underground

The space-time dynamics of the pollutant in groundwater is presented in this paragraph. To this aim, we use a miscible displacements model for flow in porous media. The area containing croplands, subsurface and groundwater collection wells is represented by a bounded domain Ω ⊂ R N (N ≤ 3 for our application). We assume that the boundary ∂Ω of Ω is such that ∂Ω ∈ C 2 . Time horizon is denoted by T , with 0 < T ≤ ∞. We set Ω T = Ω×]0, T [. Pollutant concentration in groundwater, c(t, x), is driven by the velocity v(t, x) of the mixture.

Here x ∈ Ω stands for position and t ∈ [0, T ] for time. Due to the mass conservation principle, concentration c satisfies ϕ∂ t (c + R(c)) + div(vc) -div(ϕS(v)∇c) = -r(c) + γ + p, [START_REF] Lankoski | Innovations in nonpoint source pollution policy-European perspectives, The magazine of food, farm, and resource[END_REF] where p is the fertilizer load. The other source term, γ, enables to take into account contribution from soil itself and other inputs. It is represented by a non-negative function such that

γ ∈ L ∞ (Ω T ).
The dispersion tensor is S(v). Following [START_REF] Scheideger | The Physics of Flow through Porous Media[END_REF], we consider a nonlinear dependence of the longitudinal and transverse components of the dispersion with regard to the velocity: tensor

S(v) is such that 1 S(v) = S m Id + S p (v), S p (v) = |v| α L v ⊗ v v.v + α T (Id - v ⊗ v v.v ) ( 2
)
where S m , α L and α T are respectively the diffusion coefficient, the longitudinal and transverse dispersion factors, with

S m > 0, 0 ≤ α T ≤ α L .
Soil porosity is ϕ ∈ L ∞ (Ω). We denote by ϕ -and ϕ + real numbers such that

0 < ϕ -≤ ϕ(x) ≤ ϕ + , a.e. x ∈ Ω. (3)
A chemical reaction (adsorption) of the pollutant with the soil is included in the model through the term R(c). Following the arguments in [START_REF] De Marsily | Quantitative Hydrogeology: Groundwater Hydrology for Engineers[END_REF] (page 251), we assume that R is a linear function, thus

1 + R := R ∈ R.
Coefficient R is the so-called retardation factor. Another chemical reaction, transforming the pollutant into other miscible species, is also taken into account. It is represented by the term 1 Here u ⊗ v denotes the tensor product defined by (u ⊗ v) ij = u i v j while u • v denotes the scalar product

u • v = N i=1 u i v i . The idendity matrix is Id.
r(c). Classical isotherms (see e.g. [START_REF] Williams | Combustion Theory: the fundamental Theory of chemically reacting Flow Systems[END_REF]) are characterized by linear models, in the form r(c) = kc, or by Freundlich models, in the form r(c) = kc k , or by Langmuir models, in the form r(c) = kc/(1 + k c), (k, k ) ∈ R 2 + (k commonly less than unity for Freundlich's isotherms). In what follows, we thus assume without restriction that r is a continuous function satisfying for some

r + ∈ R + |r(x) x| ≤ r + |x| 2 , ∀x ∈ R.
Such an assumption covers all the isotherms classically employed (see [START_REF] Bear | Netherlands, Theory and Applications of Transport in Porous Media[END_REF]).

We now have to derive an equation ruling the velocity v.

ϕR∂ t c + div(vc) -div(ϕS(v)∇c) = r(c) + g(t, x, p). (4)
This component is either produced as the result of the reaction of the first pollutant (thus the term r(c)), or brought in by natural inputs and/or in the fertilizer mixture (thus the source term g(t, x, p)). We assume that g is a non-negative function such that

(t, x) → g t, x, p(t, x) ∈ L ∞ (Ω T ) for any p ∈ L ∞ (Ω T ).
Summing up equations ( 1) and ( 4) and bearing in mind that c + c = 1 yield the following equation for the velocity v:

div(v) = γ + p + g. (5)
It is the expression of the total mass conservation principle. Its homogeneous form is known as the incompressibility equation.

Using div(vc) = v • ∇c + c div(v) and (5) in equation ( 1), one gets:

ϕR∂ t c + v • ∇c -div(ϕS(v)∇c) = -r(c) + (γ + p)(1 -c) -gc.
We complete the model with the Darcy law, which links, in porous media, the hydraulic head φ with the velocity v. We obtain the following set of coupled equations for the unknowns c and

φ ϕR∂ t c + v • ∇c -div(ϕS(v)∇c) = -r(c) + (γ + p)(1 -c) -gc, (6) div
(v) = γ + p + g, v = -κ∇φ. (7)
Here κ ∈ L ∞ (Ω; R N ) is the mobility of the fluid, taking into account the permeability of the underground and the viscosity of the fluid. Let us assume that there exists (κ -, κ

+ ) ∈ R 2 , 0 < κ -≤ κ + , such that κξ • ξ ≥ κ -|ξ| 2 and |κξ| ≤ κ + |ξ| , ∀ξ ∈ R N . ( 8
)
Notice that information on the heterogeneity of the soil (especially its microscopic geometry) is contained in functions ϕ and κ.

We complete the system ( 6)-( 7) by initial and boundary conditions. In the following lines, we mostly use Neumann conditions for c (no-flux) and Dirichlet conditions for φ, namely

S(v)∇c • n = 0 on ∂Ω×]0, T [, (9)
φ = φ 1 on ∂Ω×]0, T [, ( 10
)
where n is the unit exterior normal to ∂Ω and φ 1 ∈ L ∞ (0, T ; L 2 (∂Ω)). Choosing such mixed conditions for c and φ leads indeed to the most technical analysis. For the sake of simplicity, we still denote φ 1 an appropriate lifting such that

φ 1 ∈ L ∞ (0, T ; H 1 (Ω)).
Nevertheless, we also consider other choices for boundary conditions, namely:

c = c 1 on ∂Ω×]0, T [, (11)
-κ∇φ • n = v 1 on ∂Ω×]0, T [, (12)
where c 1 and v 1 may be associated to appropriate liftings, with the regularity

c 1 ∈ H 1 (0, T ; L 2 (Ω) ∩ L 2 (0, T ; W 1,∞ (Ω)) if α T = 0, c 1 ∈ H 1 (0, T ; L 2 (Ω) ∩ L 2 (0, T ; W 1,4 (Ω)) if α T > 0,
and v 1 ∈ L 2 (∂Ω). In the present paper the analysis is performed for the possible combinations of boundary conditions ( 9)- [START_REF] Scheideger | The Physics of Flow through Porous Media[END_REF]. We could of course use other conditions. First, by dividing the boundary ∂Ω into different parts associated with different boundary conditions. We could also consider Robin type boundary conditions, in the form

∂ n u + σ u u = u 1 , σ u ∈ R, u = c or φ. Their
analysis is not detailed in the present paper for the sake of brevity. It would result of the same arguments as those presented here. Initial condition is given in L 2 (Ω) by

c(0, x) = c 0 (x), x ∈ Ω. ( 13
)
The latter data satisfies the following conditions, which are physically realistic for any concentration:

0 ≤ c 0 (x) ≤ 1 a.e. x ∈ Ω, 0 ≤ c 1 (t, x) ≤ 1 a.e. (x, t) ∈ ∂Ω×]0, T [.

The Control Problem

Let p > 0 be the maximal fertilizer load that can be applied on the field. Such a quantity exists due to obvious practical constraints imposed to the farmer. The natural admissible set of control is defined as

E 0 = p ∈ L 2 (Ω T ); 0 ≤ p(t, x) ≤ p a.e. in Ω T .
It appears for technical reasons that the optimal control should rather be looked for in the subset

E ⊂ E 0 defined by E = p ∈ L 2 (Ω T ) ∩ W 1,1 (0, T ; X); 0 ≤ p(x, t) ≤ p a.e. in Ω T , α∂ t p L 1 (0,T ;X) ≤ C p ,
and

Ω (γ + p + g) dx = ∂Ω v 1 dσ if (12) is used
where C p ≥ 0 is a given real number, X is a functional space such that the embedding L 2 (Ω) ⊂ X is continuous, and α is defined by

α =          0 if g = g(t,
x, p) is such that p + g does not depend on p

1 else. ( 14
)
Notice that case α = 1 corresponds to settings where the fertilizer load affects the velocity of the underground fluid.

We consider a standard central planner objective, which takes into accounts benefits due to fertilizers use on the one hand and cleaning costs on the other hand. The instantaneous benefits is modeled by a function f (x, p), such that

∂ 2 f ∂p 2 < 0.
Function f depends on x because fertilizer loads are localized on some part of the soil surface.

We consider that cleaning costs are described by a function D depending on the position of the production wells and on the pollutant concentration at this position. Function D satisfies:

(t, x) → D(x, c(t, x)) ∈ L 1 (Ω T ) if c ∈ L 2 (Ω T ), ∂D ∂c > 0 and ∂ 2 D ∂c 2 > 0.
Such a specific localization of the fertilizer load and of the production wells has never been included in previous works on spatial economic models of pollution. The delay type approaches model this gap between pollution emission and pumping by a fixed time delay. In the present work, such an empirical choice is avoided by associating the control problem with the full model for the mixture displacement in the aquifer described in the previous section.

The intertemporal objective is thus written as

J(p) = T 0 Ω (f (x, p) -D(x, c)) dx e -ρt dt
where ρ is the social discount rate, 0 < ρ < 1, and, to sum up, we consider the following problem.

Definition 2.1 (Problem P) Find

max p∈E J(p)
subject to state constraints

         ϕR∂ t c + v.∇c -div(ϕS(v)∇c) = -r(c) + (γ + p)(1 -c) -gc, div(v) = γ + p + g, v = -κ∇φ, (15)
in Ω T , completed by the initial condition ( 13) and boundary conditions (( 9) or ( 11)) and ( [START_REF] De Frutos | Pollution control in a multiregional setting: A differential game with spatially distributed controls[END_REF] or ( 12)).

Existence of a Bounded Solution for a Given Control

In this section, for any given p in the admissible set E, we prove the existence of a solution to the state problem [START_REF] Bear | Netherlands, Theory and Applications of Transport in Porous Media[END_REF] with initial and boundary conditions chosen in ( 9)-( 13). Let us denote this problem P(p).

When handling with estimates, we denote by C a generic constant which only depends on the data of the problem.

Definition

First step should consist in giving an appropriate definition of a weak solution of problem P(p).

For the sake of concision, we only give below the definition when considering Neumann boundary conditions for c and Dirichlet boundary conditions for φ.

Definition 3.1 For any given p ∈ E, (c, φ) ∈ C(0, T ; L 2 (Ω))∩L 2 (0, T ; H 1 (Ω))×L ∞ (0, T ; H 1 (Ω))
is a weak solution of problem P(p) consisting in ( 15), ( 9)-( 10), [START_REF] De Marsily | Quantitative Hydrogeology: Groundwater Hydrology for Engineers[END_REF] if:

(i) the concentration c is such that, for any ψ ∈ H 1 (Ω T ) with ψ |t=T = 0,

- Ω T ϕRc∂ t ψ dxdt - Ω Rc 0 ψ |t=0 dx + Ω T (v • ∇c)ψ dxdt + Ω T ϕS(v)∇c • ∇ψ dxdt = Ω T -r(c) + (γ + p)(1 -c) -gc ψ dxdt; (16) (ii) the hydraulic head φ is such that φ-φ 1 ∈ L ∞ (0, T ; H 1 0 (Ω)) and, for any Ψ ∈ L ∞ (0, T ; H 1 0 (Ω)), Ω T κ∇φ • ∇Ψ dxdt = Ω T (γ + p + g)Ψ dxdt. (17)

Existence and Uniqueness of the Hydraulic Head φ

As, for a given p in E, the hydraulic head φ dynamics is not coupled with the concentration dynamics, we first study its behavior. We recall that it is ruled by the equation

div(-κ∇φ) = γ + p + g in Ω T , (18)
completed by [START_REF] De Frutos | Pollution control in a multiregional setting: A differential game with spatially distributed controls[END_REF] or [START_REF] Scheideger | The Physics of Flow through Porous Media[END_REF] on ∂Ω×]0, T [. We claim and prove the following result.

Proposition 3.1 There exits a unique function φ ∈ L ∞ (0, T ; H 1 (Ω)) solving problem [START_REF] Amann | Dynamic theory of quasilinear parabolic systems, III. Global existence[END_REF], whatever the boundary conditions [START_REF] De Frutos | Pollution control in a multiregional setting: A differential game with spatially distributed controls[END_REF] or [START_REF] Scheideger | The Physics of Flow through Porous Media[END_REF]. Moreover, for Dirichlet boundary conditions (10), we have: Proof The proof of existence and uniqueness relies on the Lax-Milgram theorem. It is standard (see for instance [START_REF] Ladyzhenskaya | Linear and Quasilinear Elliptic Equations[END_REF]) thanks to the first property in [START_REF] Brock | Diffusion-induced instability and pattern formation in infinite horizon recursive optimal control[END_REF] and to the regularity assumed for φ 1 or v 1 . Let us now focus on additional regularity results. Setting f = γ + p + g + div(-κ∇φ 1 ),

1. If φ 1 ∈ W 2,p ( Ω) with p > N and κ ∈ (C 1 ( Ω)) N ×N , then φ belongs to L ∞ (0, T ; W 2,p (Ω)). In particular, v = -κ∇φ belongs to L ∞ (Ω T ). 2. If φ 1 ∈ W 2,p ( Ω) with p > N/2, then φ belongs to L ∞ (Ω T ). If moreover κ = κ * Id with κ * : Ω → R, κ * ∈ C 1 ( Ω), then φ belongs to L ∞ (0, T ; W 1,q (Ω)) with q > N . In particular, v = -κ∇φ belongs to L ∞ (Ω T ).
u = φ -φ 1 solves the homogenous Dirichet problem div(-κ∇u) = f in Ω T , u = 0 on ∂Ω×]0, T [.
We begin by the first item listed in the proposition. If

φ 1 ∈ W 2,p ( Ω), then f ∈ L p (Ω).
According to Nirenberg's theorem, it follows from this regularity (see e.g. [START_REF] Ladyzhenskaya | Linear and Quasilinear Elliptic Equations[END_REF]) that u ∈ L ∞ (0, T ; W 2,p 0 (Ω))

and thus φ ∈ L ∞ (0, T ; W 2,p (Ω)). Since p > N , it follows from Sobolev's embeddings that

φ ∈ L ∞ (0, T ; W 1,∞ (Ω)) and thus v ∈ L ∞ (Ω T ).
In the second item, the first assertion is a classical Stampacchia's result. If moreover κ = κ * Id,

κ * ∈ C 1 ( Ω), we have as above φ ∈ L ∞ (0, T ; W 2,p (Ω)
). We know that the Darcy velocity satisfies

div(v) = -div(κ∇φ) = γ + p + g ∈ L ∞ (Ω T ). Due to the specific structure of κ, we have furthermore curl(v) = -curl(κ∇φ) = -∇κ ∧ ∇φ ∈ L ∞ (0, T ; W 1,p (Ω)) → L ∞ (0, T ; L q (Ω))
where q > N since p > N/2. In view of the regularity of the divergence and of the curl of v, we claim that v belongs to L ∞ (0, T ; W 1,q (Ω)) with q > N . The result follows from (8).

Existence of the Concentration c

In the present subsection, for a given p ∈ E and thus a given φ (see Proposition 3.1), we perform the mathematical analysis of the concentration part of problem (15) recalled below:

ϕR∂ t c + v.∇c -div(ϕS(v)∇c) = -r(c) + (γ + p)(1 -c) -gc, (19)
completed by one of the following set of conditions:

S(v)∇c • n = 0 on ∂Ω×]0, T [, c(0, x) = c 0 (x) in Ω, (20)
c = c 1 on ∂Ω×]0, T [, c(0, x) = c 0 (x) in Ω. (21)
For explaining the strategy of the analysis below, let us emphasize that the main difficulty is the control of the convective term in [START_REF] Amann | Erratum Dynamic theory of quasilinear parabolic systems, III. Global existence[END_REF]. This mathematical difficulty can not be pushed aside by physical arguments since the fluid displacement is very often convection dominated. We thus begin by the study of an auxiliary problem, which is actually the same but assuming that the velocity belongs to L ∞ (Ω T ) (which is not always ensured, see Proposition 3.1). Next we show that the existence of a weak solution for problem P(p) is ensured if we pay the price at some point: if N > 2, handling with the convective term requires that we ensure either the control in L ∞ of the velocity or a maximum principle for the concentration. As explained later things are easier if the boundary conditions for c and φ are of the same type.

We begin by the following auxiliary result. Proof Since the dispersion tensor S is assumed in the form [START_REF] Xepapadeas | The Economics of Non-Point-Source Pollution[END_REF] with

S m > 0, α L ≥ α T ≥ 0,
showing the existence of the solution c to [START_REF] Amann | Erratum Dynamic theory of quasilinear parabolic systems, III. Global existence[END_REF], (20) or ( 19), (21) with v ∈ L ∞ (Ω T ) is quite straightforward, for instance using the Faedo-Galerkin method. Let us concentrate on the possible maximum principles. We write the proof for problem [START_REF] Amann | Erratum Dynamic theory of quasilinear parabolic systems, III. Global existence[END_REF], (20) (the one for ( 19), (21) follows the same lines, using moreover 0 ≤ c 1 ≤ 1 a.e. for eliminating the boundary terms). We begin with item 1. Since we aim at proving that c ≥ 0 a.e. in Ω T , we can consider for the sake of simplicity and without loss of generality that the reactive term is replaced by

r(c) =        r(c) if c > 0, r(0) else.
Let c -= min{0, c}. We multiply ( 19) by c -and we integrate over Ω×]0, τ [, 0 < τ ≤ T . Since

cc -= (c -) 2 and Ω S(v)∇c • ∇c -dx = Ω S(v)∇c -• ∇c -dx ≥ S m Ω |∇c -| 2 dx, we obtain R 2 τ 0 d dt Ω ϕ c -2 dxdt + τ 0 Ω (v • ∇c)c -dxdt + ϕ -S m τ 0 Ω |∇c -| 2 dxdt + τ 0 Ω r(c) -γ c -dxdt - τ 0 Ω pc -dxdt + τ 0 Ω γ + p + g |c -| 2 dxdt ≤ 0.
Since γ, p, g and -c -are non-negative, since c -(x, 0) = min{c 0 (x), 0} = 0, and thanks to the definition of r, the latter relation reduces to

R 2 Ω ϕ c -(τ, x) 2 dx + τ 0 Ω (v • ∇c)c -dxdt + ϕ -S m τ 0 Ω |∇c -| 2 dxdt + τ 0 Ω r(0) -γ c -dxdt ≤ 0.
Using the Cauchy-Schwarz and the Young inequalities, we compute

τ 0 Ω (v • ∇c)c -dxdt ≤ ϕ -S m 2 τ 0 Ω |∇c -| 2 dxdt + v 2 L ∞ (Ω T ) 2ϕ -S m τ 0 Ω |c -| 2 dxdt.
Thus we get

R 2 Ω ϕ c -(τ, x) 2 dx + ϕ -S m 2 τ 0 Ω |∇c -| 2 dxdt + τ 0 Ω r(0) -γ c -dxdt ≤ C τ 0 Ω |c -| 2 dxdt and, if r(0) -γ ≤ 0, R 2 Ω ϕ c -(τ, x) 2 dx + ϕ -S m 2 τ 0 Ω |∇c -| 2 dxdt ≤ C τ 0 Ω |c -| 2 dxdt.
Thanks to (3) and to c -(x, 0) = 0, we apply the Gronwall lemma and we conclude that (H5) the boundary conditions for the hydraulic head are given by [START_REF] Scheideger | The Physics of Flow through Porous Media[END_REF] and: either v 1 may be associ-

ated to a lifting such that v 1 ∈ (L ∞ (Ω)) N , div(v 1 ) ∈ L ∞ (Ω); or v 1 is essentially non-negative in the sense that min{0, v 1 } L ∞ (∂Ω) < ϕ -S m
Then there exists a weak solution to problem P(p) in the sense of Definition 3.12 . Function φ is unique and we have the additional regularity results for function c: α

T |v| 1/2 ∇c ∈ (L 2 (Ω T )) N , c ∈ H 1 (0, T ; W -1,4 (Ω)).
Moreover, we can assert that the velocity v belongs to (L ∞ (Ω T )) N if (H1) or (H2) is satisfied;

-if r(0) -γ ≤ 0 then c(t, x) ≥ 0 a.e. (x, t) ∈ Ω T ; -if r(1) + g ≥ 0 then c(t, x) ≤ 1 a.e. (x, t) ∈ Ω T .
Remark 3.2 Some of these assumptions could be weakened but we have chosen to restrict the statement to cases where the uniqueness of the solution is also ensured (see Proposition 3.3 below).

Proof The existence of φ is already ensured by Proposition 3.1. Furthermore, if assumption (H1)

or (H2) is satisfied, we know by Proposition 3.1 that the velocity v belongs to L ∞ (Ω T ). The existence of c is then stated in Lemma 3.1. It remains to explore cases where v does not belong to L ∞ (Ω T ). To this aim, we introduce a regularization v of v, for instance by convolution, such that in particular

v ∈ (L ∞ (Ω T )) N for any > 0, v → v in L 2 (Ω T ) as → 0.
For any > 0, we associate with v the solution c of

ϕR∂ t c + v • ∇c -div(ϕS(v )∇c ) = -r(c ) + (γ + p)(1 -c ) -gc in Ω T , (22)
S(v )∇c • n = 0 on ∂Ω×]0, T [, c (0, x) = c 0 (x) in Ω. ( 23
)
We now look for sufficient uniform estimates to pass to the limit → 0 in ( 22)-( 23). Let τ ∈]0, T [.

We multiply (22) by c (t, x)χ [0,τ ] (t) and we integrate by parts over Ω T . Thanks to the definition of tensor S(v), we obtain

R 2 Ω ϕ |c (τ, x)| 2 dx + τ 0 Ω ϕ -(S m + α T |v |)|∇c | 2 dx dt + τ 0 Ω (v • ∇c )c dx dt + τ 0 Ω r(c )c dx dt - τ 0 Ω (p + γ)(1 -c )c dx dt + τ 0 Ω g(c ) 2 dx dt ≤ R 2 Ω ϕ |c 0 | 2 dx. ( 24
)
As p, γ, g ∈ L ∞ (Ω T ) and |r(x)x| ≤ r + x 2 for any x ∈ R, Cauchy-Schwarz and Young inequalities yield

τ 0 Ω r(c )c dx dt + τ 0 Ω (p + γ)(1 -c )c + g(c ) 2 dx dt ≤ C + C τ 0 Ω |c | 2 dx dt.
We thus infer from (24) that

Rϕ - 2 Ω |c (τ, x)| 2 dx + τ 0 Ω ϕ -(S m + α T |v |)|∇c | 2 dx dt ≤ τ 0 Ω (v • ∇c )c dx dt + C + C τ 0 Ω |c | 2 dx dt. ( 25
)
As already mentioned, the difficulty for the estimates lies in the convective term. We distinguish three cases, corresponding to assumptions (H3)-(H5):

Assume first (H3). Using the Cauchy-Schwarz and Young inequalities, we write

τ 0 Ω (v • ∇c )c dx dt ≤ 1 2ϕ -α T τ 0 Ω |v ||c | 2 dx dt + ϕ -α T 2 τ 0 Ω |v | |∇c | 2 dx dt ≤ C v L ∞ (0,τ ;L 2 (Ω)) τ 0 c (t, •) 2 L 4 (Ω) dt + ϕ -α T 2 τ 0 Ω |v | |∇c | 2 dx dt ≤ C τ 0 c (t, •) 2 L 4 (Ω) dt + ϕ -α T 2 τ 0 Ω |v | |∇c | 2 dx dt.
The Gagliardo-Nirenberg inequality reads, if

N = 2, c (t, •) 2 L 4 (Ω) ≤ C c (t, •) L 2 (Ω) c (t, •) H 1 (Ω)
a.e. in ]0, T [. It gives in the latter relation:

τ 0 Ω (v • ∇c )c dx dt ≤ C τ 0 c (t, •) 2 L 2 (Ω) dt + C τ 0 c (t, •) L 2 (Ω) ∇c (t, •) (L 2 (Ω)) N dt + ϕ -α T 2 τ 0 Ω |v | |∇c | 2 dx dt.
Finally, the Cauchy-Schwarz and Young inequalities give

τ 0 Ω (v • ∇c )c dx dt ≤ C τ 0 Ω |c | 2 dx dt + ϕ -S m 2 τ 0 Ω |∇c | 2 dx dt + ϕ -α T 2 τ 0 Ω |v | |∇c | 2 dx dt. ( 26
)
Assume now (H4). Since the boundary condition is now c = c 1 on ∂Ω×]0, T [, we rather multiply (22) by c -c 1 and we integrate by parts other Ω×]0, τ [. The convective term thus gives:

τ 0 Ω (v • ∇c )(c -c 1 ) dx dt = τ 0 Ω (v • ∇(c -c 1 ))(c -c 1 ) dx dt + τ 0 Ω (v • ∇c 1 )(c -c 1 ) dx dt ≤ τ 0 Ω v • ∇((c -c 1 ) 2 /2) dx dt + C + τ 0 Ω |c | 2 dx dt = - 1 2 τ 0 Ω div(v )(c -c 1 ) 2 dx dt + C + τ 0 Ω |c | 2 dx dt ≤ C + C τ 0 Ω |c | 2 dx dt. ( 27
)
The last estimate is based on the fact that div(v ), where v is a regularization by convolution of v, may be estimated as div(v) = γ + p + g ∈ L ∞ (Ω T ). Thus,

τ 0 Ω div(v )c 2 dx dt ≤ C τ 0 Ω |c | 2 dx dt. If (H5) is satisfied, we compute τ 0 Ω (v • ∇c )c dx dt = 1 2 τ 0 Ω v • ∇(c 2 ) dx dt = - 1 2 τ 0 Ω div(v )c 2 dx dt + 1 2 τ 0 ∂Ω (v • n)c 2 |∂Ω dσ dt.
Since, once again, div(v ) is uniformly bounded in L ∞ (Ω T ), the first term in the righthand side of the latter relation does not bring new difficulty. The second term is

τ 0 ∂Ω v 1 c 2 |∂Ω dσ dt/2. It is non-negative if v 1 ≥ 0 and then it does not modify (25). If min{0, v 1 } L ∞ (∂Ω) < ϕ -S m , we
have thanks to the trace theorem

1 2 τ 0 ∂Ω v 1 c 2 |∂Ω dσ dt < ϕ -S m c |∂Ω 2 L 2 (∂Ω×(0,T )) ≤ ϕ -S m c 2 L 2 (0,T ;H 1 (Ω)) .
The latter bound may be absorbed in the lefthand side of (25). Finally, if we assume that

v 1 = v1 • n with v1 ∈ (L ∞ (Ω)) N and div(v 1 ) ∈ L ∞ (Ω), we write 1 2 τ 0 ∂Ω v 1 c 2 |∂Ω dσ dt = τ 0 Ω (v 1 • ∇c )c dx dt + 1 2 τ 0 Ω div(v 1 )c 2 dx dt, thus - 1 2 τ 0 Ω div(v )c 2 dx dt + 1 2 τ 0 ∂Ω v 1 c 2 |∂Ω dσ dt ≤ C τ 0 Ω |c | 2 dx dt + ϕ -S m 2 τ 0 Ω |∇c | 2 dx dt. (28)
Now, inserting (26) or ( 27) or (28) in (25), we get, whatever the assumption (H3), (H4) or (H5):

Rϕ - 2 Ω |c (τ, x)| 2 dx + ϕ - 2 τ 0 Ω (S m + α T |v |)|∇c | 2 dx dt ≤ C + C τ 0 Ω |c | 2 dx dt.
Since ϕ -> 0, we then use the Gronwall lemma and we conclude that the following uniform estimates hold true

c L ∞ (0,T ;L 2 (Ω))∩L 2 (0,T ;H 1 (Ω)) ≤ C, (29)
α T |v | 1/2 ∇c (L 2 (Ω T )) N ≤ C. (30)
It follows that there exists a subsequence, still denoted c for convenience, and some function

c ∈ L ∞ (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 1 (Ω)) such that c c weakly in L 2 (0, T ; H 1 (Ω)).
This result is however not sufficient to pass to the limit in the nonlinear equation ( 22). For reaching some compactness result for c , we first prove that ϕ∂ t c is uniformly bounded in L 2 (0, T ; W -1,4 (Ω)). This point may be easily checked by multiplying ( 22) by some test function ψ in L 2 (0, T ; W 1,4 0 (Ω)) and by integrating by parts over Ω T . Indeed (for some lines symbol ∈ also means "uniformly bounded"

), |v | 1/2 ∈ L ∞ (0, T ; L 4 (Ω)) follows from v ∈ (L ∞ (0, T ; L 2 (Ω))) N .
Then, on the one hand, since

α T |v | 1/2 ∇c ∈ (L 2 (Ω T )) N , S(v )∇c ∈ (L 2 (0, T ; L 4/3 (Ω))) N . On the other hand, v •∇c ∈ L 2 (0, T ; L 1 (Ω)) and thanks to Sobolev embedding, ψ ∈ L 2 (0, T ; L ∞ (Ω)).
Thus ,4 (Ω), the first embedding being compact, we can use an advanced version of the Aubin compactness criterium (see for instance Lemma 3.1 in [START_REF] Galusinski | On a degenerate parabolic system for compressible, immiscible, two-phase flows in porous media[END_REF]) for proving that c is sequentially compact in L p (0, T ; L 2 (Ω)) for any 1 ≤ p < 2. In particular, there exists c ∈ C(0, T ;

ϕ∂ t c L 2 (0,T ;W -1,4 (Ω)) ≤ C. ( 31
) Since 0 < ϕ -≤ ϕ(x) ≤ ϕ + a.e. x ∈ Ω, since H 1 (Ω) ⊂ L 2 (Ω) ⊂ W -1
L 2 (Ω)) such that c → c in L 2 (Ω T ), c → c a.e. in Ω T .
Using in particular Lebesgue's convergence theorem, we pass to the limit → 0 in problem ( 22 Proof The uniqueness of φ has already been proved in Proposition 3.1. Assume that there exist two solutions c 1 and c 2 of ( 19)-( 20) such that

c i (Ω T ) ⊂ I, i = 1, 2. Let c = c 1 -c 2 .
This function satisfies the equation

Rϕ∂ t c + v • ∇c -div(S(v)∇c) + r(c 1 ) -r(c 2 ) + (p + γ)c + gc = 0 in Ω T , (32)
completed by one of the two following conditions

S(v)∇c • n = 0 on ∂Ω×]0, T [, c |t=0 = 0 in Ω, (33)
c = 0 on ∂Ω×]0, T [, c |t=0 = 0 in Ω. (34)
Let τ ∈]0, T [. We multiply (32) by χ [0,τ ] c and integrate over Ω T . We get, using the properties of S(v), φ and c |t=0 = 0:

Rϕ - 2 Ω |c(τ, x)| 2 dx + τ 0 Ω ϕ -(S m + α T |v|)|∇c| 2 dx dt + τ 0 Ω (v • ∇c)c dx dt + τ 0 Ω (r(c 1 ) -r(c 2 ))c dx dt + τ 0 Ω (p + γ + g)c 2 dx dt = 0.
Using the additional property of r, we make the following estimate

τ 0 Ω (r(c 1 ) -r(c 2 ))c dx dt ≤ r ∞ τ 0 Ω |c 1 -c 2 | |c| dx dt = C τ 0 Ω |c| 2 dx dt.
Once again, it remains to estimate the convective term. If v belongs to L ∞ (Ω T ) (assumption (H1) or (H2)), we write

τ 0 Ω (v • ∇c)c dx dt ≤ S m ϕ - 2 τ 0 Ω |∇c| 2 dx dt + C v 2 ∞ 2S m ϕ - τ 0 Ω |c| 2 dx dt.
If N = 2 and α T > 0 (assumption (H3)), we get the analogous to (26)

τ 0 Ω (v • ∇c)c dx dt ≤ C τ 0 Ω |c| 2 dx dt + ϕ -S m 2 τ 0 Ω |∇c| 2 dx dt + ϕ -α T 2 τ 0 Ω |v| |∇c| 2 dx dt.
If homogeneous Dirichlet conditions (34) are imposed on c (assumption (H4)), then

τ 0 Ω (v • ∇c)c dx dt = - 1 2 τ 0 Ω div(v)c 2 dx dt = - 1 2 τ 0 Ω (γ + p + g)c 2 dx dt and τ 0 Ω (v • ∇c)c dx dt ≤ γ + p + g L ∞ (Ω T ) τ 0 Ω |c| 2 dx dt.
Finally, if Neumann conditions are prescribed for φ (assumption (H5)), we have

τ 0 Ω (v • ∇c)c dx dt = - 1 2 τ 0 Ω div(v)c 2 dx dt + 1 2 τ 0 ∂Ω (v • n)c 2 |∂Ω dσ dt = - 1 2 τ 0 Ω (γ + p + g)c 2 dx dt + 1 2 τ 0 ∂Ω v 1 c 2 |∂Ω dσ dt.
We have already estimated the first term in the righthand side of the latter relation. For the second one, either we write

τ 0 ∂Ω v 1 c 2 |∂Ω dσ dt = τ 0 ∂Ω max{0, v 1 } + min{0, v 1 } c 2 |∂Ω dσ dt where τ 0 ∂Ω max{0, v 1 }c 2 |∂Ω dσ dt ≥ 0 and τ 0 ∂Ω min{0, v 1 }c 2 |∂Ω dσ dt < ϕ -S m c 2 L 2 (∂Ω×(0,T )) ≤ ϕ -S m τ 0 Ω |∇c| 2 dx dt + ϕ -S m τ 0 Ω |c| 2 dx dt,
or we write

1 2 τ 0 ∂Ω v 1 c 2 |∂Ω dσ dt = τ 0 Ω (v 1 • ∇c)c dx dt - 1 2 τ 0 Ω div(v 1 )c 2 dx dt ≤ S m ϕ - 2 τ 0 Ω |∇c| 2 dx dt + C v1 2 ∞ 2S m ϕ - T 0 Ω |c| 2 dx dt + div(v 1 ) ∞ 2 T 0 Ω |c| 2 dx dt.
We infer from all the latter estimates that there exists some non-negative real number C such that

Rϕ - 2 Ω |c(τ, x)| 2 dx ≤ C τ 0 Ω |c| 2 dx dt.
We conclude with the Gronwall lemma that c = 0 and thus c 1 = c 2 almost everywhere in Ω T .

The uniqueness result is proved.

Remark 3.3

The interested reader may check that all the classical reaction models (for instance those mentioned in Subsection 2.1) are such that r has a bounded derivative on any bounded subset of R.

Existence of an Optimal Control

We now have the tools for the study of the optimal control problem. We adopt a fixed point strategy for conciliating the maximization problem with the system of hydrogeological constraints.

In the latter section, we have studied the constraints system when it is decoupled of the optimization question. We begin this section by the following auxiliary result where the optimization problem is considered without constraints. Proof Since c is given, the problem is equivalent to: find p ∈ E such that p = arg min q∈E J(q)

where J is defined for q ∈ L 2 (Ω T ) by

J (q) = - T 0 Ω f (x, q) dx e -ρt dt.
Now the proof is standard. First, E is a closed convex subset of L 2 (Ω T ) (for the strong topology of L 2 (Ω T )). Next one checks easily (using Lebesgue's theorem of convergence) that J is a continuous function for the same topology. Finally one infers from the strict concavity of f the strict convexity of J.

Let us first assume that time horizon T is finite. The first main result of the paper is the following.

Theorem 4.1 Let T > 0. Assume that one of the assumptions listed in Proposition 3.2 is satisfied. Assume that there exists a convex set I ⊂ R such that the derivative of r exists and is bounded in I. Furthermore, assume that I = R, or that r(0) -γ ≤ 0 and R + ⊂ I, or that r(1) + g ≥ 0 and ] -∞, 1] ⊂ I, or that r(0) -γ ≤ 0 and r(1) + g ≥ 0 and [0, 1] ⊂ I. Then there exists (p * , c * , φ * ) solving problem P as described in Definition 2.1. The control p * belongs to E.

Properties of (c * , φ * ) are given in Propositions 3.

Proof Since E is a closed bounded subset of L 2 (Ω T ) , E is a compact set for the weak topology of L 2 (Ω T ) . Moreover E is a convex set. In view of the existence and uniqueness results proved in Lemma 4.1, and Propositions 3.2-3.3, we can define T : E → E by

T : p 0 ∈ E → p = arg max q∈E Ω T (f (x, q) -D(x, c
)) e -ρt dx dt where (c, φ) is the unique solution of

         Rϕ∂ t c + v • ∇c -div(ϕS(v)∇c) = -r(c) + (γ + p 0 )(1 -c) -gc in Ω T , div(v) = γ + p 0 + g, v = -κ∇φ in Ω T , (35)
with boundary conditions ( 9) or ( 11) and ( 10) or [START_REF] Scheideger | The Physics of Flow through Porous Media[END_REF], and initial condition [START_REF] De Marsily | Quantitative Hydrogeology: Groundwater Hydrology for Engineers[END_REF].

Now the important point is to prove the continuity of T , here for the weak topology of

L 2 (Ω T ). Let (p 0,n ) be a sequence in E such that p 0,n p 0 weakly in L 2 (Ω T ).
In order to prove that T is continuous, we need to prove that T (p 0,n ) T (p 0 ) weakly in L 2 (Ω T ). Let p n = T (p 0,n ) and p = T (p 0 ). It means that p n = arg max q∈E Ω T (f (x, q) -D(x, c n )) e -ρt dx dt where (c n , φ n ) is the weak solution of

Rϕ∂ t c n + v n • ∇c n -div(ϕS(v n )∇c n ) = -r(c n ) + (γ + p 0,n )(1 -c n ) -gc n in Ω T , (36) div(v n ) = γ + p 0,n + g, v n = -κ∇φ n in Ω T , (37)
with boundary conditions ( 9) or ( 11) and ( 10) or [START_REF] Scheideger | The Physics of Flow through Porous Media[END_REF], and initial condition [START_REF] De Marsily | Quantitative Hydrogeology: Groundwater Hydrology for Engineers[END_REF].

Multiplying (37) by φ n and integrating by parts over Ω, we immediately get the following uniform estimate:

φ n L ∞ (0,T ;H 1 (Ω)) ≤ C.
Thus there exist a subsequence of φ n , not relabeled for convenience, and some function φ ∈

L ∞ (0, T ; H 1 (Ω)) such that φ n φ weakly in L q (0, T ; H 1 (Ω)), v n v = -κ∇φ weakly in L q (0, T ; (L 2 (Ω)) N ), ∀q ≥ 1.
Let us prove that ∇φ n actually strongly converges to ∇φ. Let us first assume that α = 0 which means that the source term p 0,n + g in (37) does not depend on p 0,n . Then φ n = φ does not depend on n and the result is straightforward. Assume now α = 1. Then φ n -φ is defined by

-div(κ∇(φ n -φ)) = p 0,n -p 0 in Ω T , φ n -φ = 0 on ∂Ω.
We multiply the latter equation by φ n -φ and we integrate by parts over Ω T . Thanks to (8), we write

κ - T 0 Ω |∇(φ n -φ)| 2 dx dt ≤ T 0 Ω (p 0,n -p 0 )(φ n -φ) dx dt. (38)
Since we only know that p 0,n -p 0 0 weakly in L 2 (Ω T ), we can not pass directly to the limit in the righthand side of (38). Nevertheless we know, because of the definition of E, that ∂ t (p 0,n -p 0 ) is uniformly bounded in L 1 (0, T ; X). Yet function ∂ t (φ n -φ) satisfies

-div κ∇ ∂ t (φ n -φ) = ∂ t (p 0,n -p 0 ) in Ω T , ∂ t (φ n -φ) = 0 on ∂Ω×]0, T [.
Thus ∂ t (φ n -φ) is uniformly bounded in L 1 (0, T ; X). We already knew that φ n -φ is uniformly bounded in L ∞ (0, T ; H 1 (Ω)). Since H 1 (Ω) ⊂ L 2 (Ω) ⊂ X, the first embedding being compact, the classical Aubin's theorem implies that there exists a subsequence of φ n (still denoted the same below) such that

φ n -φ → 0 in L 2 (Ω T ).
According to (38), we thus have

lim n→∞ T 0 Ω |∇(φ n -φ)| 2 dx dt = 0.
This point proves that φ n → φ in L 2 (0, T ; H 1 (Ω)), thus v n → v in (L 2 (Ω T )) N and almost everywhere in Ω T . Now, following the lines of the study of problem ( 22)-( 23), we can show that there

exists c ∈ L 2 (0, T ; H 1 (Ω)) such that c n c weakly in L 2 (0, T ; H 1 (Ω)) and c n → c in L 2 (Ω T )
and almost everywhere in Ω T . The latter convergence results are sufficient to pass to the limit n → ∞ in system (36)-(37). We retrieve (35) completed by the initial and boundary conditions (thanks to the continuity of the trace function). We have chosen the setting of Proposition 3.3 which ensures that the solution of this problem is unique. Thus the whole sequence (c n , φ n , v n ) converges towards (c, φ, v).

Up to now, we have focused on the limit behavior of the constraints system. Let us go back to the whole optimization problem. Since p n = T (p 0,n ) belongs to E, it is uniformly bounded in L 2 (Ω T ). There exists p ∈ L 2 (Ω T ) such that (up to some subsequence)

p n p weakly in L 2 (Ω T ).
We aim at proving that p = p = T (p 0 ). For any q ∈ E, by definition of T ,

Ω T f (x, p n ) -D(x, c n ) e -ρt dx dt ≥ Ω T f (x, q) -D(x, c n ) e -ρt dx dt. ( 39
)
Function f being concave and p n p weakly in L 2 (Ω T ), a lower semi-continuity argument gives Ω T f (x, p)e -ρt dx dt ≥ lim n→∞ Ω T f (x, p n )e -ρt dx dt. Since D is continuous with regard to the concentration variable and since c n → c in L 2 (Ω T ), we thus have

Ω T f (x, p) -D(x, c) e -ρt dx dt ≥ lim n→∞ Ω T f (x, p n ) -D(x, c) e -ρt dx dt = lim n→∞ Ω T f (x, p n ) -D(x, c n ) e -ρt dx dt.
Using (39), we get, for any q ∈ E,

Ω T f (x, p) -D(x, c) e -ρt dx dt ≥ lim n→∞ Ω T f (x, q) -D(x, c n ) e -ρt dx dt = Ω T f (x, q) -D(x, c) e -ρt dx dt.
The latter relation shows that p is actually equal to the unique solution of the arg max problem p = T (p 0 ). Hence T (p 0,n ) T (p 0 ) and T is continuous (for the weak topology of L 2 (Ω T )).

We conclude the proof using the Schauder fixed point theorem. Since E is a convex and compact set for the weak topology of L 2 (Ω T ), since T is a continuous function for the same topology, T has a fixed point. This fixed point is the optimal solution exhibited in Theorem 4.1.

This ends the proof.

The previous result has been obtained for a finite T. The following result applies for an infinite planning horizon. Proof Because of assumptions r(0) -γ ≤ 0 and r(1) + g ≥ 0, any solution given in Theorem 4.1 satisfies 0 ≤ c * ≤ 1 a.e. in Ω×]0, T [ (see Lemma 3.1). The control in L ∞ of the norm of (p * , c * , φ * ) then ensures that any maximal weak solution is actually a global solution (see also [START_REF] Amann | Dynamic theory of quasilinear parabolic systems, III. Global existence[END_REF][START_REF] Amann | Erratum Dynamic theory of quasilinear parabolic systems, III. Global existence[END_REF]).

is controllable for the functional J by the set of admissible controls L 2 (Ω T ), for any T > 0. On the one hand, replacing the set of admissible controls E 0 by L 2 (Ω T ) does not allow anymore to prove that the concentration satisfies 0 ≤ c(t, x) ≤ 1 almost everywhere in Ω T , whatever the assumptions on r, γ and g. But on the other hand, we can prove the asymptotic controllability of (42). Recall that a system Lc = F (q) is asymptotically controllable in a Banach space B by a set of admissible controls E adm if, for any arbitrary c ∈ B, there exists a sequence of controls

p n ∈ E adm such that lim n→∞ c(p n ) -c B = 0
where we denote by c(q) the solution of Lc = F (q).

We claim and prove the following result for problem (42).

Theorem 5.1 Assume that one of the assumptions listed in Proposition 3.2 is satisfied. Problem

(42) is asymptotically controllable in B = {f ∈ L 2 (0, T ; H 1 (Ω)), S(v)∇f • n = 0 or f = c 1 on ∂Ω, f | t=0 = c 0 }
by the set E adm = L 2 (Ω T ), for any T > 0 and for any given c 0 ∈ L 2 (Ω).

Proof Let T > 0. Let v being given in (L 2 (Ω T )) N such that (41) holds. Let c ∈ B. By density arguments, there exist sequences c

1,n ∈ C 1 (Ω T ), c 0,n ∈ C 1 (Ω) and c n ∈ C 1 (Ω T ) with S(v)∇c n • n = 0 or c n = c 1,n on ∂Ω×]0, T [ and c n | t=0 = c 0,n such that lim n→∞ c 1,n -c 1 H 1 (0,T ;L 2 (Ω))∩L 2 (0,T ;W 1,∞ (Ω)) = 0, lim n→∞ c 0,n -c 0 L 2 (Ω) = 0, (43)
and

lim n→∞ c -c n L 2 (0,T ;H 1 (Ω)) = 0. ( 44
)
Let L be the operator defined in C 1 (Ω T ) by

Lc = ϕR∂ t c -div(ϕS(v)∇c) + v • ∇c + r(c) -γ(1 -c) + gc. One checks straightforward that Lc n ∈ L 2 (0, T ; H -1 (Ω)) for any n ∈ N. Since L 2 (Ω T ) is dense in L 2 (0, T ; H -1 (Ω)) there exists a sequence p n,k ∈ L 2 (Ω T ) such that Lc n -p n,k L 2 (0,T ;H -1 (Ω)) ≤ k for any n ∈ N , where lim k→∞ k = 0. ( 45
)
Let c(p n,k ) be the solution of (42) for q = p n,k . Our proof of Theorem 5.1 consists in showing that lim n,k→∞ c(p n,k ) -c L 2 (0,T ;H 1 (Ω)) = 0. To this aim, we intend to compute the limit of each term in the righthand side of the following relation

c -c(p n,k ) L 2 (0,T ;H 1 (Ω)) ≤ c -c n L 2 (0,T ;H 1 (Ω)) + c n -c(p n,k ) L 2 (0,T ;H 1 (Ω)) . (46)
Due to (44), we only have to study c n -c(p n,k ) L 2 (0,T ;H 1 (Ω)) . Function C n defined by

C n = c n -c(p n,k ) satisfies ϕR∂ t C n -div(ϕS(v)∇C n ) + v • ∇C n = -r(c n ) -r(c(p n,k )) -(γ + g)C n + (Lc n -p n,k ) in Ω T , S (v)
Rϕ - 2 Ω |(C n -C 1,n )(τ, x)| 2 dx + ϕ - τ 0 Ω S m + α T |v| |∇(C n -C 1,n )| 2 dx dt ≤ R 2 Ω |C 0,n (x) -C 1,n (0, x)| 2 dx - τ 0 Ω v • ∇(C n -C 1,n ) (C n -C 1,n ) dx dt + r + γ + g ∞ τ 0 Ω |C n -C 1,n | 2 dx dt + r + γ + g ∞ τ 0 Ω |C 1,n | |C n -C 1,n | dx dt + τ 0 Ω Lc n -p n,k -R∂ t C 1,n -v • ∇C 1,n + div(ϕS(v)∇C 1,n ) (C n -C 1,n ) dx dt. ( 47
)
Let us estimate the terms in the righthand side of the latter relation. We have Following the lines of the treatment of the convective term in the proof of Proposition 3.2, we prove that for some 0 < < 1,

τ 0 Ω v • ∇(C n -C 1,n ) (C n -C 1,n ) dx dt ≤ ϕ - τ 0 Ω S m + α T |v| |∇(C n -C 1,n )| 2 dx dt +C C 1,n 2
H 1 (0,T ;L 2 (Ω)) + C τ 0 Ω |C n -C 1,n | 2 dx dt.
Thanks to Cauchy-Schwarz and Young inequalities, we have

r + γ + g ∞ τ 0 Ω |C 1,n | |C n -C 1,n | dx dt ≤ C C 1,n 2
L 2 (Ω T ) + C τ 0 Ω |C n -C 1,n | 2 dx dt.
Next, using moreover (45), Since k → 0 as k → ∞, we can choose k large enough to ensure that 1 -k ≥ C > 0. We then infer from the Gronwall lemma that

C n -C 1,n 2
L ∞ (0,T ;L 2 (Ω)) + ∇(C n -C 1,n ) (L 2 (Ω T )) N ≤ CL n + C k e CT .
Inserting this result in (46), we get c -c(p n,k ) L 2 (0,T ;H 1 (Ω)) ≤ c -c n L 2 (0,T ;H 1 (Ω))

+ C T k + CL n ,
where, in view of (43), ( 44) and (45), lim n,k→∞ c -c n L 2 (0,T ;H 1 (Ω)) + C T k + L n = 0. The theorem is proved.

Conclusion

In this paper, the existence of a global optimal solution has been proved as well as an asymptotic controllability result. An outlook of this work is to extend this fine hydrogeological modeling to a game situation, with several farmers. Such a game situation is worth studying as far as the non-cooperative situation is the one that prevails in real life. It would be then interesting to compare the outcome of the Nash equilibrium strategy that occurs in a non-cooperative game to a cooperative situation as the one presented in this paper to determine whether any distortion appears. If it is the case, optimal fiscal policy (taxation, subvention, ...) should have to be computed.

  For the sake of simplicity, let us assume that the mixture contains only one other component. It means that we are only interested in one specie, the pollutant. Denote by c the concentration of the other component. It satisfies an equation similar to (1):

Remark 3 . 1

 31 Of course, there exist additional regularity results for Neumann boundary conditions (see e.g.[START_REF] Ladyzhenskaya | Linear and Quasilinear Elliptic Equations[END_REF] chapter 10). Nevertheless they are not necessary for our purpose.

Lemma 3 . 1

 31 Assume p ∈ E and v ∈ L ∞ (Ω T ). There exists a weak solution c to problem (19), (20) or (19), (21). Furthermore, we have 1. If r(0) -γ ≤ 0 and c 0 ≥ 0 a.e. in Ω, then c ≥ 0 a.e. in Ω T . 2. If r(1) + g ≥ 0 and c 0 ≤ 1 a.e. in Ω, then c ≤ 1 a.e. in Ω T .

τ 0 ΩProposition 3 . 2

 032 |c -| 2 dxdt = 0 a.e. τ ∈]0, T [ that is c -= 0 a.e. in Ω T . Item 1 is proved. The proof of item 2 is similar using the test function max{c -1, 0} instead of c -. Now we turn back to the generic setting for the concentration problem. We claim and prove the following existence result. Let p ∈ E. Assume that one of the following assumptions is satisfied, (H1) the boundary conditions are (9)-(10), φ 1 ∈ W 2,p ( Ω) with p > N and κ ∈ (C 1 ( Ω)) N ×N (H2) the boundary conditions are (9)-[START_REF] De Frutos | Pollution control in a multiregional setting: A differential game with spatially distributed controls[END_REF],φ 1 ∈ W 2,p ( Ω) with p > N/2, κ = κ * Id with κ * : Ω → R, κ * ∈ C 1 ( Ω)(H3) the boundary conditions are (9)-[START_REF] De Frutos | Pollution control in a multiregional setting: A differential game with spatially distributed controls[END_REF], N = 2 and α T > 0 (H4) the boundary conditions for the concentration are given by[START_REF] Benosman | Optimization of cytostatic leukemia therapy in an advection reaction diffusion model[END_REF], whatever it is associated with[START_REF] De Frutos | Pollution control in a multiregional setting: A differential game with spatially distributed controls[END_REF] or[START_REF] Scheideger | The Physics of Flow through Porous Media[END_REF]

Proposition 3 . 3

 33 )-(23) and we prove that c is a weak solution of problem (19)-(20). The additional regularity result comes from (30)-(31) and from the maximum principle statement in Lemma 3.1. Proposition 3.2 is proved. 3.4 Uniqueness of the Weak Solution of P(p) Assume that there exists a convex set I ⊂ R such that the derivative of r exists and is bounded in I. Then, under one of the assumptions listed in Proposition 3.2, any weak solution (c, φ) of problem P(p) satisfying c(Ω T ) ⊂ I is unique.

Lemma 4 . 1

 41 Let c ∈ L 2 (Ω T ) be given. The problem find p ∈ E such that p = arg max q∈E T 0 Ω (f (x, q) -D(x, c)) dx e -ρt dt admits a unique solution p ∈ E.

Proposition 4 . 1

 41 If assumptions of Theorem 4.1 are satisfied, if moreover r(0) -γ ≤ 0 and r(1)+g ≥ 0, then there exists a global solution (p * , c * , φ * ) to problem P as described in Definition 2.1 such that 0 ≤ c * (t, x) ≤ 1 a.e. in [0, ∞[×Ω.

  ∇C n • n = 0 or C n = C 1,n = c 1,n -c 1 on ∂Ω×]0, T [, C n|t=0 = C 0,n = c 0,n -c 0 in Ω.Let us now detail only the computations when choosing Dirichlet boundary condition, that is C n = C 1,n on ∂Ω×]0, T [. This is indeed the most technical case since it imposes to treat a boundary term. We multiply the first equation of the latter system by C n -C 1,n and we integrate by parts in Ω×]0, τ [, τ ∈]0, T [. Since in particular the derivative of r is assumed bounded, we write

R 2 Ω 2 L 2 ( 2 H 1

 22221 |C 0,n (x) -C 1,n (0, x)| 2 dx ≤ C C 0,n Ω) + C C 1,n (0,T ;L 2 (Ω)) .

( 0 Ω≤ C + C τ 0 ΩFinally, we write τ 0 Ω≤ C C 1,n 2 H 1 0 Ω 2 ΩC τ 0 Ω

 00021020 Lc n -p n,k ) (C n -C 1,n ) dx dt ≤ Lc n -p n,k L 2 (0,T ;H -1 (Ω)) C n -C 1,n L 2 (0,T ;H 1 (Ω)) ≤ k C n -C 1,n L 2 (0,T ;H 1 (Ω)) ≤ k 4ϕ -S m + k ϕ -S m C n -C 1,n 2 L 2 (0,T ;H 1 (Ω)) ≤ k 4ϕ -S m + k ϕ -S m τ |C n -C 1,n | 2 dx dt + k ϕ - τ 0 Ω S m |∇(C n -C 1,n )| 2 dx dt |C n -C 1,n | 2 dx dt + k ϕ - τ 0 Ω (S m + α T |v|) |∇(C n -C 1,n )| 2 dx dt. R∂ t C 1,n + v • ∇C 1,n -div(S(v)∇C 1,n ) (C n -C 1,n ) dx dt ≤ C C 1,n H 1 (0,T ;L 2 (Ω))∩L 2 (0,T ;W 1,∞ (Ω)) C n -C 1,n L 2 (Ω×(0,τ )) (0,T ;L 2 (Ω))∩L 2 (0,T ;W 1,∞ (Ω)) + τ |C n -C 1,n | 2 dx dt.Due to (43), all the latter estimates give in (47)Rϕ - |(C n -C 1,n )(τ, x)| 2 dx + (1 -k )ϕ - τ 0 Ω S m + α T |v| |∇(C n -C 1,n )| 2 dx dt ≤ |C n -C 1,n | 2 dx dt + C k + CL nwhere lim n→∞ L n = 0.

or comparable to Definition

3.1 if boundary conditions are chosen in (11)-(12).

Asymptotic Controllability

In the present section, we restrict the problem to the cases where equations ( 6) and ( 7) are decoupled, meaning that the fertilizer load does not impact the global dynamic of the underground flow. We thus consider the following problem

)) dx e -ρt dt, where c satisfies:

the velocity v being the solution of div v = g + γ, v = -κ∇φ, with boundary conditions [START_REF] De Frutos | Pollution control in a multiregional setting: A differential game with spatially distributed controls[END_REF] or [START_REF] Scheideger | The Physics of Flow through Porous Media[END_REF].

(41)

Quantity v thus may be considered as a data of the problem, given in (L 2 (Ω T )) N by Proposition 3.1. Notice that we have relaxed the constraints imposed on the set where we look for the control p which is now L 2 (Ω T ). Indeed, for problem (40) the parameter α defined in ( 14) is α = 0.

Assuming p ∈ L ∞ (Ω T ) is also useless. Following the lines of the proof of Theorem 4.1, we easily recover the following existence result. of the optimal control problem (40) such that c * ∈ L 2 (0, T ; H 1 (Ω)) ∩ H 1 (0, T ; W -1,4 (Ω)).

The latter result means that problem