Emmanuelle Augeraud-Véron
email: emmanuelle.augeraud@u-bordeaux.fr

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

Éloïse Comte
email: eloise.comte@inria.fr

Optimal buffer zone for the control of groundwater pollution from agricultural activities

Keywords: Pollution control, Water supply, Buffer zone, Nonlinear parabolic and elliptic equations, Optimal control, Upscaling process Mathematics Subject Classification (2010) MSC 49A22, 49B22, 91B76, 91B15, 76N25, 76S05

We consider an optimal control model of groundwater pollution due to agricultural activities, the objective of the optimal manager being the optimization of the trade-off between the fertilizer use and the cleaning costs. The size of a buffer zone defined around the water production sites for limiting pollution effects may also be chosen by the manager. The spread of the pollution from crops to captation wells is modeled using a convection-diffusion-reaction equation. The main hydrogeological features of the dynamics are taken into account: the process is convection-dominated, dispersion effects are included, generic nonlinear reaction terms are considered. The existence and the uniqueness of the optimal solution is proven. Using asymptotic analysis, we then rigorously prove that a one-dimensional static model can be substituted to the full dynamic 3D model for the long time study of the optimal solution. We prove that this new optimal control problem is also well-posed. All these theoretical results are used for characterizing the optimal buffer zone.

Introduction

Groundwater pollution is a worldwide major concern. In France for instance, between years 1994 and 2003, 7716 water captation wells had to be abandoned ( [START_REF]L'eau et les milieux aquatiques -Chiffres clés -Edition[END_REF]), mainly because of unacceptable nitrates concentrations. Yet 66% of nitrate found in French groundwater is due to agriculture ( [START_REF]dossiers scientifiques. L'eau ; Dégradations ; La pollution par les nitrates[END_REF]). This proportion reaches up to 80% in some European regions (Lankoski and Ollikainen [3]). Despite the application of the European Directive (91/676/CEE) and Water Framework Directive (2000/60/EC), a drastic number of captation wells still do not provide drinkable water (the nitrates concentration is above 50mg/L). Denitrification processes are expensive (0.5 euros/m 3 , according to Miquel and de Marsily [START_REF] Miquel | Rapport de l'Office parlementaire d'évaluation des choix scientifiques et technologiques[END_REF]). Such a critical situation holds in many part of the world. Therefore the trade-off between benefit from fertilizer use and cleaning cost is still a question that has to be tackled. The design of protected areas near boreholes, where polluting activities are prohibited, is one of the main effective management tool against water pollution (Correll [5]). It is worthwhile to study the impact of these buffer zones even when they are not active zones (Leandri [6]).

Since the seminal paper of Keeler et al. [START_REF] Keeler | The optimal control of pollution[END_REF], optimal management of nonpoint source pollutions has been widely studied. As it is stressed in Kossioris et al. [START_REF] Kossioris | Feedback Nash equilibria for non-linear differential games in pollution control[END_REF], advances in environmental economics emphasize the need for a realistic representation of the ecological system. Yet, the sensitivity of the problem to geographical elements, that is the development of models depending on space and time -and not solely on time-is seldom addressed. However, the space dependance is an obvious component of the problem since crops and wells are not located as the same place. The space dimension is all the more important that the flow velocity in aquifers is very low (Kim et al. [START_REF] Kim | The regulation of groundwater quality with delayed responses[END_REF]; in most of the aquifers, the maximal velocity is typically in the range from 10 meters per year to 100 meters per year, Lalbat et al. [START_REF] Lalbat | Description of a simple hydrochemical indicator to estimate groundwater residence time in carbonate aquifers[END_REF]) thus leading to an important time delay between the emission and the captation of the pollution in water production wells.

Nevertheless the first attempts for including the space dimension to optimal control and management models were unfortunately reduced to time ODE's models: Winkler ([11], [START_REF] Winkler | A note on the optimal control of stocks accumulating with a delay[END_REF]) and Augeraud-Véron and Leandri ([13]) consider an exogenous time-delayed equation for the pollution concentration.

Spacial dimension has been considered in optimal pollution control in shallow lakes, in a one-dimensional setting, by Brock and Xepapadeas [START_REF] Brock | Diffusion induced instability and pattern formation in infinite horizon recursive optimal control[END_REF]. Patterns formation is studied using a (linear) reaction diffusion system. Unlike shallow lakes models, groundwater models have to be ruled primary by advective phenomena. De Frutos and Martin-Herran [START_REF] De Frutos | Spatial effects and strategic behavior in a multiregional transboundary pollution dynamic game[END_REF] consider a groundwater pollution in a multiregional decentralized framework. They assume that the fluid velocity and the diffusion coefficient are given real numbers. They illustrate the behavior of the system through several numerical examples. However, due to microscopic soil heterogeneity, the diffusion operator should be replaced by a dispersion operator, nonlinearly depending on the velocity of the fluid which may be itself an unknown of the problem.

In the present paper, our starting point is a full time and space dependent model. After its mathematical analysis, we rigorously prove that for large time studies, corresponding in particular to the political time scale, a static model also makes sense, unlike most of the models developed in the literature which are rather only time dependent ODEs. We give some precisions in the next.

A high heterogeneity of the space and time scales is inherent to this kind of problem: compare for instance the large time scale of the political effects with the seasonality of agricultural activities, the small thickness or the small fluid displacements with regards to the length of most of the aquifers. In the present work, we rigorously derive an effective model adapted to this multiscale setting. The upscaling process enables us to reduce the problem to a static one with moreover the economy of one space-dimension. As emphasized by Lions in [START_REF] Lions | Asymptotic problems in distributed systems, Riviere Memorial Lecture[END_REF], the challenge of asymptotic questions arising in connection with the optimal control of PDEs systems is important. Nevertheless, settings leading to some kind of degeneracy in the structure of the problem were seldom addressed. For thin structures, we mention first Fonseca and Francfort who present in [START_REF] Fonseca | 3D-2D asymptotic analysis of an optimal design problem for thin films[END_REF], for an unconstrained problem, a setting where 3D-2D asymptotics and optimal design commute from a variational standpoint, using Γ -convergence theory. For constrained optimization problem in thin domains, we cite Kogut and Leugering [START_REF] Kogut | Optimal Control Problems for Partial Differential Equations on Reticulated Domains, Approximation and Asymptotic Analysis[END_REF] (with a linear parabolic state equation) and Casado-Díaz et al. [START_REF] Casado-Díaz | A rigorous asymptotic expansion for a control problem in a thin domain[END_REF] (with a linear elliptic state equation) who worked on the system formed by the direct and adjoint state equations. This is not the approach adopted in the present paper. Furthermore, since we also take into account a large time asymptotic, the singular perturbation problem induced by the thin structure has to be coupled with the time degeneracy of the problem. Basically, our parabolic state equation degenerates into an elliptic type equation. The present work thus meets the questions raised in [START_REF] Porretta | Long time versus steady state optimal control[END_REF] by Porretta and Zuazua. Up to our knowledge, the coupled large times and thin structure asymptotic analysis performed here for an optimal control problem was never addressed.

After the literature review let us go back to the present article's substance. It is devoted to an optimal control problem of groundwater pollution. The problem is constrained by the displacement model for contaminants in the aquifer. On the one hand, theoretical results are provided, with a real contribution from a mathematical point of view. On the other hand, qualitative results are developed for the analysis of the impact of a buffer zone (a fallow for instance).

The main theoretical result is the rigorous asymptotic analysis of the optimal control problem and of its solution. One of the main qualitative results is the analysis of the impact of the pollution at the captation well regarding the length of the buffer zone. We also compute the optimal size of the buffer zone in a social welfare approach.

The paper is organized as follows. In Section 2, the objective, taking into account the trade-off between fertilizer use benefits and cleaning costs of polluted water, is defined. It is associated with a hydrogeological state equations system describing the transport of the pollutant in groundwater. General boundary conditions are considered, allowing for distinct boundaries constraints for the two state variables, the fluid velocity and the pollutant concentration. The existence and the uniqueness of the solution of the corresponding optimal control problem is stated in Section 3. In Section 4, we explain how a one-dimensional static optimal control problem can be substituted to the full 3D dynamic problem for a wide range of large times studies (more precisely for times larger than the longitudinal characteristic time and for order one Peclet numbers). We prove that this new optimal control problem is also well-posed. Finally, optimality conditions are derived and closed form solutions are computed. Then the optimal buffer zone is characterized. The results are illustrated by numerical simulations. All the proofs have been postponed to the last section of the paper for favoring a more straightforward reading of the results.

The model

The domain We consider a three-dimensional area containing croplands, a groundwater reservoir and water collection wells. It is represented by a bounded domain Ω ⊂ R 3 . The cropland includes the subsoil's first centimeters and it is denoted S ⊂ Ω ⊂ R 3 . Indeed, due to the run-off process, the characteristic time of the displacement in this area of the domain is very small, either way negligible with regard to the characteristic time of the horizontal displacement in the whole aquifer (see Abbott et al. [START_REF] Abbott | An introduction to the European hydrological system "she" 2: Structure of a physically-based, distributed modelling system[END_REF]). We assume that the boundary ∂Ω of Ω is smooth, say C 2 for instance. Time horizon is denoted by T , with 0 < T ≤ ∞. We set Ω T = (0, T ) × Ω. If (t, x) ∈ Ω T , x stands for the location and t for the time.

The unknowns The pollutant concentration and the fertilizer load at time t and location x, are respectively denoted by c(t, x) and p(t, x), while the fluid velocity is v(t, x) ∈ R 3 .

The set of controls Let p > 0 be the maximal fertilizer load that can be applied on the field. Let p be a given real number. Because of the partial differential equations used for the hydrogeological model, the natural admissible set of control is

E T = {q ∈ L 2 (Ω T ), p ≤ q ≤ p a.e. in Ω T } if T < ∞, {q ∈ L 2 (Ω × (0, ∞)), p ≤ q ≤ p a.e. in Ω × (0, ∞)} if T = ∞. (1)
At first glance, the reader could find difficult to interpret cases where p < 0, especially since the control variable is a fertlizer "load". We thus give the interpretation of a negative load. In the cultivated area, the case q < 0 may correspond to the intercropping practice (for absorbing the surplus of nutrient, see Willey [START_REF] Willey | Resource use in intercropping systems[END_REF]) especially with service plants with the capacity to capture soil nitrogen and restore it to the next crop, see GEVES [START_REF] Geves | French Variety and Seed Study and Control Group, Service plants[END_REF] or Louison [START_REF] Louison | Modeling plant nutrient uptake: mathematical analysis and optimal control[END_REF]. For the non-cultivated area, the case q < 0 corresponds to an active buffer zone.

The objective function We consider a classical social optimization problem. The welfare is represented by a discounted utilitarian function which depends on the one hand on the private benefit (for the farmers) and on the other hand, on the cleaning cost of the environmental damage due to fertilizer pollution.

The objective of the central manager is denoted by J. Namely,

J(p, c) =    T 0 Ω f (x, p)χ S -D(x, c) e -ρt dx dt -νe -ρT Ω D(x, c(T, x) dx if T < ∞ ∞ 0 Ω f (x, p)χ S -D(x, c) e -ρt dx dt otherwise.
(

) 2
The instantaneous private benefit function Ω f (x, p)χ S dx is computed on the cropland part of the domain, thus the use of the characteristic function of the set S, χ S . Crops variety can be considered since f depends on the spatial variable x. For ensuring the integrability of the terms involving f in (2) it appears that it is sufficient to assume f such that x → f (x, q(x)) ∈ L1 (Ω) for any function q valued in (p, p). Function f dependency on the load fertilizer p is standard:

Assumption 1 Function p → f (•, p) belongs to C 1 (R)
, is increasing and strictly concave 1 .

The cleaning cost function D depends of course on the pollutant concentration but also on the location of the captation wells and on the pumping rates. We thus assume that D depends on x and c:

Assumption 2 Function D : Ω × (d -, ∞) → R + , where d -∈ R, is such that for a.e. x ∈ Ω, c ∈ (d -, ∞) → D(x, c) is a bounded, continuous function ∂ c D ≥ 0 and ∂ 2 cc D ≥ 0; ∃D + ∈ R + s.t. D(x, c(t, x)) L 1 (Ω T ) ≤ D + c L 2 (Ω T ) .
When the horizon is finite, T < ∞, the remaining cost of the pollution is taken into account in the scrapping value νe -ρT Ω D(x, c(T, x)) dx. It depends on parameter ν ≥ 0.

State dynamics

The pollutant concentration c is driven by a velocity field v, which is created by the space variations of the time and space dependent hydraulic head φ. The hydraulic head is a measure of both the fluid pressure and the gravity effects. The dynamics of the state variables c and φ, together with their dependency on the control p are described by the following set of equations in Ω T (see Augeraud-Véron et al. [START_REF] Augeraud-Véron | Optimal control for a groundwater pollution ruled by a convection-diffusion-reaction problem[END_REF] for the link between these PDEs and the conservation principles)

Rψ∂ t c + v • ∇c -div(ψS(v)∇c) = -r(c) -gc + pχ S + γ, (3)
div v = g, v = -κ∇φ. (4)
These equations describe a convection reaction diffusion process in the porous medium corresponding to soil and groundwater, taking into account the dispersion, a kind of diffusion depending on the velocity. Notice that v and φ are autonomous variables, depending exogenously on g. The soil porosity is characterized by the function ψ, which solely depends on the location x. Assumptions for ψ are the following. Assumption 3 Function ψ ∈ L ∞ (Ω) is such that there exists two real numbers (ψ -, ψ + ) satisfying 0 < ψ -≤ ψ(x) ≤ ψ + a.e. in Ω.

Displacements in porous media, hence the dynamics of (3), are dominated by the convective term v • ∇c. Moreover the diffusive term in (3), div(ψS(v)∇c), nonlinearly depends on the velocity v. Indeed, the classical diffusion operator is perturbed by the heterogeneity of the fluid velocity at the microscopic scale, where the soil appears like a highly heterogeneous medium made of voids and grains. We thus have to deal with a dispersion factor S(v) in the form

S(v) = S m Id + S p (v)
where the purely diffusive part contains the identity matrix Id and a diffusion coefficient S m , while the dispersive component S p (v), which is physically the leading order term, depends on the velocity v. More precisely, its longitudinal and transverse components are described as follows 2 (see Scheidegger [START_REF] Scheidegger | The Physics of Flow through Porous Media[END_REF]).

Assumption 4 The dispersion tensor reads

S(v) = S m Id + S p (v) with S p (v) = |v| α L |v| -2 v⊗v+α T (Id-|v| -2 v⊗v) , (S m , α L , α T ) ∈ R 3 such that S m > 0, α L ≥ α T ≥ 0. In particular, S(v)ξ • ξ ≥ (S m + α T |v|) |ξ| 2 , |S(v)ξ| ≤ (S m + α L |v|) |ξ| ∀ξ ∈ R 3 .
The time evolution in (3) is scaled by a parameter, the so-called retardation factor, R > 0. This scaling is induced by the chemical interactions (sorption) between the pollutants mixture and the soil. The possible occurrence of a reaction in the mixture is also considered and is modeled by the function r.

We assume general assumptions for r which enable to encompass most of the usual reaction kinetics, such as e.g. the nonlinear Langmuir and Freundlich isotherms (Bear and Verruijt [START_REF] Bear | Modeling Groundwater Flow and Pollution, Theory and Applications of Transport in Porous Media[END_REF]). Assumptions on r are the following.

Assumption 5 Function r belongs to C 1 ((d -, ∞)), is concave and there exists

r + ∈ R + such that |r(x)x| ≤ r + |x| 2 for any x ∈ (d -, ∞).
The first equation in (4) is the incompressibility equation for the mixture. The second one is the so-called Darcy law, where κ is the mobility of the fluid. It depends on the permeability of the underground and on the viscosity of the fluid as follows.

Assumption 6 There exists (κ -, κ

+ ) ∈ R 2 , 0 < κ -≤ κ + , such that κξ • ξ ≥ κ -|ξ| 2 and |κξ| ≤ κ + |ξ| , ∀ξ ∈ R 3 .
Furthermore the function κ satisfies one of the two following assumptions:

κ ∈ (C 1 ( Ω)) 3×3 or κ = κ * Id with κ * : Ω → R and κ * ∈ C 1 (Ω).
2 Here u ⊗ v denotes the tensor product, (u

⊗ v) ij = u i v j .
The second hypothesis could appear rather technical. Notice however that it is only required for ensuring that the mean velocity of the flow remains bounded, in L ∞ (Ω T ), which is of course physically realistic. Quantity γ ≥ 0 is a pollutant source term accounting for the contributions from soil itself and from other inputs. Nonnegative function g describes all the inputs containing other chemical elements than the studied pollutant. It is for instance related to natural degradation of organic elements in soils, to rain falls or to agricultural spreading. Assumptions on g is the following:

Assumption 7 Function g : (t, x) → g(t, x, p(t, x)) belongs to L ∞ (Ω T ) for any p ∈ L ∞ (Ω T ).
Equations ( 3)-( 4) are completed by initial and boundary conditions. The initial condition is specified for c,

c |t=0 = c 0 in Ω, (5)
where we assume that c 0 ∈ H 1 (Ω). In order to disentangle boundary assumptions on the two state variables, we consider two non-overlapping decompositions of the boundary

∂Ω of Ω, ∂Ω = Γ 1 ∪ Γ 2 , Γ 1 ∩ Γ 2 = ∅, and, indepen- dently, ∂Ω = Γ 3 ∪ Γ 4 , Γ 3 ∩ Γ 4 = ∅.
We assume Neumann boundary conditions on Γ 1 × (0, T ) and on Γ 3 × (0, T ) for the state variables c and φ. We assume Dirichlet boundary conditions on Γ 2 × (0, T ) and on Γ 4 × (0, T ). Namely:

S(v)∇c • n = 0 on Γ 1 × (0, T ), c = c -on Γ 2 × (0, T ) (6) v • n = -κ∇φ • n = v 1 on Γ 3 × (0, T ), φ = 0 on Γ 4 × (0, T ), (7)
where n denotes the unit exterior normal to ∂Ω. Function v 1 is such that v 1 ∈ W 1,∞ (Ω) and satisfies the compatibility condition Ω g dx = ∂Ω v 1 dσ which is a necessary condition for the well-posedness of the elliptic problem defining the velocity v.

Finally, an assumption is specifically introduced for ensuring a maximum principle for the concentrations solving (3)-( 6) :

Assumption 8 The data of the problem are assumed such that

c -≥ d -, c 0 ≥ c -a.e. in Ω, (γ + p)(1 -c -) -r(x) -gc -≥ 0 ∀x ∈ [d -, c -].
3 The optimal control problem, well-posedness result

In the present section, we present a generic existence and uniqueness result for the optimal control problem we are interested in.

Notice that for avoiding cumbersome notations, we do not separate the cases T < ∞ and T = ∞. We thus make in the following lines a notations abuse: when the time interval of the study is unbounded, the regularity results hold true for any T < ∞. We begin by a proposition stating the well-posedness of the state equations system. We refer to Augeraud et al [START_REF] Augeraud-Véron | Optimal control for a groundwater pollution ruled by a convection-diffusion-reaction problem[END_REF] for a detailed proof 3 .

Proposition 1 Under the aforementioned assumptions, for any p ∈ E T , there exists a unique global weak solution (c, φ) of (3)-( 7) belonging to the space

L ∞ (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 1 (Ω)) × L ∞ (0, T ; W 1,∞ (Ω))
. Moreover, the concentration satisfies c(t, x) ≥ c -for a.e. (t, x) ∈ Ω T .

We now deal with the optimal control problem. First we define it explicitly. The problem, named hereafter P E T , is the following:

P E T Find (p * , c * ) such that J(p * , c * ) = max p∈E T J(p, c
) where c is the concentration associated with p given by Proposition 1.

For the reader's convenience, recall that c associated with p satisfies

Rψ∂ t c + v • ∇c -div(ψS(v)∇c) = -r(c) + pχ S + γ -gc in Ω T , div(v) = g, v = -κ∇φ in Ω T , c |t=0 = c 0 in Ω, S(v)∇c • n = 0 on Γ 1 × (0, T ), c = c -on Γ 2 × (0, T ), v • n = -κ∇φ • n = v 1 on Γ 3 × (0, T ), φ = 0 on Γ 4 × (0, T ).
Note that, since φ (and thus v) is an autonomous variable, it is given once and for all by Proposition 1.

We claim and prove (see Subsection 5.1) the following result.

Theorem 1 There exists a unique solution

(p * , c * ) ∈ E T × C([0, T ]; L 2 (Ω)) ∩ L 2 (0, T ; H 1 (Ω))
to the optimal control problem P E T .

Optimal buffer zone: a large times study

In this section, we question the possible improvement of the setting by a buffer zone, that is a protected area where polluting activities (here the use of fertilizer) are prohibited. We aim to show that the optimal size of the buffer zone can be characterized.

On the other hand we do not yet have exploited the heterogeneity of the space and time scales involved in the problem: typically the large time scale of the political effects, the smallness of the thickness with regards to the length of most of the aquifers, the smallness of the Darcy velocity. In the present section, we rigorously derive an effective model taking into account the latter scales. The upscaling process let us reduce the problem to a static setting with moreover the economy of one space-dimension. We thus recover rigorously the setting which was formally introduced by Mäler et al. in the seminal paper [START_REF] Mäler | The economics of shallow lakes[END_REF].

Assuming moreover some symmetries in the reservoir, our problem reduces to a one-dimensional static problem. We then explicitly compute the optimal solution as a function of the length of the buffer zone and we determine the optimal size of the latter.

Large times, thin reservoir and order one Péclet number: effective model

We first focus on the hydrodynamical part of the model. For the sake of the clarity, we chose a less generic setting than in the previous section. We emphasize however that the only assumptions that are necessary for the upscaling process are those mentioned in Theorem 2 below. First, as we aim at computing analytical solutions for the effective model, we force it to be one-dimensional by assuming sufficient symmetries, in the geometry of the reservoir, in the structure of the soil characteristics and in the boundary conditions so that we can reduce the domain to a two-dimensional subset Ω :

Ω = {(x , z ) ∈ R 2 ; x ∈ (0, L ), z ∈ (h , H )}.
We denote by x ∈ R (resp. z ∈ R) the horizontal (resp. vertical) variable. The topography of the reservoir is described by functions h and H . We choose Neumann boundary conditions for all the state variables and we assume that the bottom and the top of the reservoir are impermeable. It is an interesting setting because such impermeable boundary conditions mean that the pollution remains trapped in the domain. Finally, we assume that the reaction terms are linear:

r(x) = αx, α ∈ R. Our starting point is thus 4 R ψ ∂ t c + v • ∇ x ,z c -div x ,z (ψ S(v )∇ x ,z c ) = -α c -g c + p χ S + γ in Ω × (0, T ), div x ,z v = g , v = -κ ∇ x ,z φ in Ω × (0, T ), S(v )∇ x ,z c • n = 0 on ∂Ω × (0, T ), c t=0 = c 0 in Ω , v • n = v 1 on ∂Ω × (0, T )
where v 1 (x , z ) = 0 if z ∈ {h (x ), H (x )}. We now introduce appropriate scales. We use reference values denoted by the subscript R. Setting

c = c /c R , v = v /v R , x = x /L R , z = z /H R , t = t /T R , . . . (8)
we obtain the dimensionless equations

R R ψ R L 2 R S R T R Rψ∂ t c + v R L R S R v x ∂ x c + v R L R S R L R H R v z ∂ z c -ψ R ∂ x ψS xx (v R v)∂ x c - L 2 R H 2 R ψ R ∂ z ψS zz (v R v)∂ z c - L R H R ψ R ∂ x ψS xz (v R v)∂ z c + ∂ z ψS xz (v R v)∂ x c = - α R L 2 R S R αc - g R L 2 R S R gc + p R L 2 R S R c R pχ S + γ R L 2 R S R c R γ in Ω T , (9)
κ Rxx L R ∂ x (κ xx ∂ x φ) + κ Rzz L R H 2 R ∂ z (κ zz ∂ z φ) + κ Rxz H R ∂ x (κ xz ∂ z φ) + ∂ z (κ xz ∂ x φ) = - g R L R φ R g in Ω T , ( 10
) v = - φ R v R t κ Rxx L R ∂ x φ + κ Rxz H R ∂ z φ, κ Rzz H R ∂ z φ + κ Rxz L R ∂ x φ in Ω T , (11)
S(v R v)∇ x,z c • n = 0 on ∂Ω × (0, T ), c t=0 = c 0 in Ω, (12)
v • n = v 1 on ∂Ω × (0, T ), (13)
where

Ω = {(x, z) ∈ R 2 ; x ∈ (0, L), z ∈ (h, H)} such that meas(Ω T ) = O(1)
and S ⊂ Ω is the rescaled set corresponding to S . The equations involve the non-dimensional Péclet and Damköhler numbers

Pe = v R L R /S R , Da = L R α R /v R
and the time scales

T L = characteristic longitudinal time scale = L R /v R , T A = characteristic adsorption time scale = 1/R R .
Since an overwhelming majority of aquifers are very thin with regard to their length, we define a small parameter by

= H R /L R .
We fix a large reference time T R T L /T A by setting

T L /(T A T R ) = .
Next we assume order one Péclet and Damköhler numbers (roughly speaking, these assumptions mean respectively a slow fluid displacement in the aquifer and an intermediate rate of conversion of reactants, thus are very realistic). We keep a porosity of order one, ψ R = O(1), we considerate quite moderate pollutants inputs, (p

R /c R )T L = O(1) and (γ R /c R )T L = O(1)
, this latter assumption being hopefully reasonable for an aquifer used for water supply. The permeability is assumed such that

κ Rxx /L R = O( ) and O(κ Rzz ) = O(κ Rxx ).
The pressure, thus φ R , is assumed of order one. For a simple setting corresponding to the latter scaling choices, one may for instance assume

H R = 1, L R = -1 , v R = O( ), O(S R ) = O(K R ) = 1, g R = p R /c R = γ R /c R = 2 .
The uniform ellipticity properties of Assumptions 4 and 6 are preserved; for the sake of simplicity in the notations, we simply assume that S(v) and κ are diagonal. We thus set κ z = κ zz and S z (v) = S zz (v). The problem ( 9)-( 13) then reads

Rψ∂ t c + v x ∂ x c + -1 v z ∂ z c -∂ x ψS x ( v)∂ x c --2 ∂ z ψS z ( v)∂ z c = -αc -gc + pχ S + γ in Ω T , (14)
∂ x v x + -1 ∂ z v z = g, v = -t (κ x ∂ x φ, -1 κ z ∂ z φ) in Ω T , (15)
S( v)∇ x,z c • n = 0 on ∂Ω × (0, T ), c t=0 = c 0 in Ω, (16)
v • n = v 1 on ∂Ω × (0, T ) (17)
where

v 1 (x, z) = 0 if z ∈ {h(x), H(x)}.
Let us now turn back to the economical part of the problem. Once again, we simplify the model to get at the limit a tractable problem for explicit computations (the following assumptions are not necessary for the rigorous derivation of the effective model). We assume that the domain contains only one well for the water supply, located around x = L w , in the form of a cylinder C(L w , w) of radius w. We assume that the cleaning costs depends on the averaged concentration of pollutants in the well. The interest reader can check straightforward that this minor change does not impact the proof of the results stated in the latter section. The dimensional form of (2) thus reads

J (p ) = T 0 Ω f (x , p )χ S dx dz -D 1 meas(C(L w , w)) C(L w ,w) c dx dz e -ρ t dt -ν e -ρ T Ω D 1 meas(C(L w , w)) C(L w ,w) c |t =T dx dz = T R T 0 L R H R Ω f (L R x, p R p)χ S dxdz -D c R meas(C(L w /L R , w/L R )) C(L w /L R ,w/L R ) c dxdz e -ρρ R T R t dt -ν e -ρρ R T R T D c R meas(C(L w /L R , w/L R )) C(L w /L R ,w/L R ) c |t=T R T dxdz
thanks to an obvious change of variable in the integration. The radius of the well is of course smaller than the thickness of the reservoir:

meas(C(L w /L R , w/L R )) = O( ) and e -ρ T R T /T R is at most of order . Set L w = L w /L R , C = C(L w , w/L R ), f (x, y) = L R H R f (L R x, p R y), ρ = ρ T R (the latter means that ρ R T R = O(1)), D(x) = D (c R x) and ν e -ρ T R T T -1 R D meas(C ) -1 C c |t=T R T dxdz = d .
The functional to maximize in the dimensionless form of the problem reads

J (p, c) = T 0 Ω f (x, p)χ S dxdz -D meas(C ) -1 C
c dxdz e -ρt dt + d (18) where meas(C ) → 0 as → 0 and d is uniformly bounded in R.

We now aim at letting → 0 for deriving an effective model exploiting the scale heterogeneity in the singularly perturbed problem P E T R T (14)- [START_REF] Kogut | Optimal Control Problems for Partial Differential Equations on Reticulated Domains, Approximation and Asymptotic Analysis[END_REF]. We easily guess that this process induces some kind of degeneracy in the model. It appears that the price to pay consists in assuming the strict concavity of D with respect to c. Assumption 2 is thus replaced by the following one:

Assumption 9 D(c) ∈ L 1 (Ω) if c ∈ L 2 (Ω), ∂ c D ≥ 0 and ∂ 2 cc D > 0.
We now claim and prove the following results describing the effective behaviour of the optimal control problem for the aforementioned scaling choices.

Theorem 2 Assume the following scalings

H R /L R = 1, Pe = O(1), Da = O(1), T R = -1 T L /T A T L /T A , that ψ R , g R T L , (p R /c R )T L and (γ R /c R )T
L are of order one and that K Rx /L R and K Rz /L R are of order . Assume that the reaction term is linear, r(x) = αx, α ∈ R. Assume that the space average of the hydraulic head is uniformly bounded with respect to , that α + g > 0 and g ∈ L ∞ (0, T ; L 2 (Ω)), that H -h > 0 and that H, h and κ x admit bounded derivatives in (0, L). Then the following effective problem, P eff , can be substituted to the problem P E T R T :

Find (p * , c * ) ∈ E eff × L ∞ (0, T ; H 1 (0, L)) s.t. J eff (p * , c * ) = max p∈E eff J eff (p, c), (19)
E eff = {q ∈ L ∞ ((0, T ) × (0, L)) s.t. 0 ≤ q ≤ p a.e.} , (20)
J eff (p, c) = T 0 L 0 f (x, p(x))χ S dx -D(c(L w )) e -ρt dt, ( 21
)
where c is the concentration associated with p through

-∂ x S m ψ∂ x c + v∂ x c = -α(H -h)c + pχ S + γ -ḡc, (22)
∂ x v = ḡ, v = -κ∂ x φ, (23)
∂ x c = 0, -κ(x)∂ x φ = v 1 , x ∈ {0, L} (24)
with ψ(x) =

H(x) h(x) ψ dz, γ(x) = H(x)
h(x) γ dz and S is the subset of (0, L) such that S = {(x, z) ; x ∈ S, z ∈ (h(x), H(x))}.

Notice that a solution of P eff a priori depends on t because some of the coefficients in ( 22)-( 24), for instance ḡ, may be time dependent. Nevertheless time t now only appears as a parameter in the problem.

The scaling choices listed in the assumptions of Theorem 2 are those we have used to derive ( 14)- [START_REF] Kogut | Optimal Control Problems for Partial Differential Equations on Reticulated Domains, Approximation and Asymptotic Analysis[END_REF]. The latter result has to be understood in the following sense.

Proposition 2 If (p * , c * ) is the optimal solution of the problem J (p * , c * ) = max p∈E T J (p, c) where J is defined in [START_REF] Kogut | Optimal Control Problems for Partial Differential Equations on Reticulated Domains, Approximation and Asymptotic Analysis[END_REF] and c is the solution of ( 14)- [START_REF] Fonseca | 3D-2D asymptotic analysis of an optimal design problem for thin films[END_REF], then it converges to the solution of P eff defined by ( 19)-( 24) as → 0 in the following sense:

v → v in (L 2 (Ω T )) 2 and a.e. in Ω T , c * c * weakly in L 2 (0, T ; H 1 (Ω)), p * p * weakly in L 2 (Ω T ) and lim sup →0 J (p * , c * ) = J eff (p * , c * ).
The latter proposition is proved in Subsection 5.2. The proof also contains the following well-posedness result.

Proposition 3 There exists a unique optimal solution (p * , c * ) belonging to

E eff × L ∞ (0, T ; H 1 (Ω)) of Problem P eff .

Optimal conditions for the effective problem

Set Ω eff = (0, L), where L is also the location of the captation well, L = L w . Denote S = (0, L 1 ) (with 0 ≤ L 1 ≤ L) the crop zone, that is the zone where the fertilizer is used. Thus the buffer zone, that is the uncultivated zone between the field and the captation well is

[L 1 , L].
At first glance, the reader may think that this additional parameter is useless and that we could use any initialization Ω\S for the buffer zone. Indeed, if the optimal setting includes a buffer zone, it should appear in the definition of the optimal control p * . More precisely, if p * vanishes in an area around the captation well, it defines a buffer zone. This is actually the case when dealing with the time dependent multidimensional in space problem corresponding to the setting of Theorem 1 with p = 0 (see the numerical results showing how the optimization process corresponding to Theorem 1 imposes a buffer zone in Comte [START_REF] Comte | Pollution agricole des ressources en eau : approches couples hydrogologique et conomique[END_REF] chap. 7). Nevertheless it appears below that the one-dimensional static problem induces stronger regularity properties for the optimal fertilizer load p * . In particular, its smoothness prevents p * from jumping to zero at the beginning of the buffer zone (see Remark 1 below). We thus have to look for the optimal pair load-concentration (p * , c * ) for any fixed L 1 , then to maximize the objective function with regard to L 1 for characterizing the optimal buffer zone.

We simplify the notations by removing the bars. The optimal effective solution is a priori space-structured and may depend on L 1 . We make some additional assumptions for the computational tractability of the problem, so that the state system reads

-ψSc + vc = -αc + pχ S + γ -gc in (0, L), ( 25
) v = g, v = -κφ in (0, L), (26)
c (0) = c (L) = 0, (κφ )(0) = (κφ )(L) = -v 1 , (27)
where S and ψ are given positive real numbers. The prime denotes the space derivation (we omit the dependence on the time parameter). We further consider v 1 ≥ 0 to ensure that the pollution is driven to the captation well.

The elliptic problem ( 26)-( 27) is solvable provided the compatibility condition L 0 g dx = 0. Once again, for simplifying the analytic computations that follow, we assume that g is a given positive real number. Due to the compatibility condition above, it implies that g = 0 and the fluid velocity is

v = v 1 . (28)
We set ρ = 0 and we no longer assume that f depends on x in the social welfare so that, due moreover to the non dependence on time of the coefficients in ( 25)-( 27), the study of J eff defined in ( 21) reduces to the one of

J eff (p, c) = L1 0 f (p(x)) dx -D(c(L)). ( 29
)
For the computations in the dual problem, Assumption 9 is replaced by the following one:

Assumption 10 D(c) ∈ L 1 (0, L) if c ∈ L 2 (0, L), ∂ c D > 0 and ∂ 2 cc D > 0.
Since the existence and the uniqueness of the optimal solution have been proved in Proposition 3, the optimal solution may be characterised by the first order necessary optimality conditions. They are presented in the following lemma (proved in Subsection 5.3 using the Lagrangian of the problem as in Benosman et al. [START_REF] Benosman | Optimization of cytostatic leukemia therapy in an advection-reaction-diffusion model[END_REF]).

Lemma 1 Let (p * , c * ) ∈ E eff × H 2 (0, L) be the optimal solution of the effective problem5 (25)- [START_REF] Bear | Modeling Groundwater Flow and Pollution, Theory and Applications of Transport in Porous Media[END_REF]. There exists µ * ∈ H 2 (0, L) solving the following dual problem 6ψSµ * + vµ * -αµ * = 0 in (0, L),

ψSµ * (0) + vµ * (0) = 0, Sψµ * (L) + vµ * (L) = D (c * (L)), (30)
f (p * )χ S = µ * χ S in (0, L). ( (31)
) 32
Remark 1 According to [START_REF] Kennedy | Dynamic programming: applications to agriculture and natural resources[END_REF] where f is assumed to be smooth, the optimal load p * inherits of the smoothness of the solution µ * of the ODE problem (30)- [START_REF] Barbu | Optimal control of population dynamics[END_REF]. In particular p * is a continuous function.

We now have the tools for computing the analytic expression of the optimal solution (p * , c * ) exhibited in Proposition 3, using the PDEs [START_REF] Augeraud-Véron | Optimal control for a groundwater pollution ruled by a convection-diffusion-reaction problem[END_REF], the first relation in ( 27) and ( 30)- [START_REF] Barbu | Optimal control of population dynamics[END_REF] completed by v = v 1 and (32). The following result is obtained (see Subsection 5.4).

Lemma 2 Assume moreover that the derivative f is positively valued. Let S = Sψ. The effective solution (p * , c * ) is explicitly given by:

p * (x, A, L 1 ) = (f ) -1 (Aϕ (x)) if x ∈ [0, L 1 ] 0 else , (33)
S (r 2 -r 1 ) c * (x, A, L 1 ) = - e -r2x -(r 2 /r 1 )e -r1x e -r2L -e -r1L × L 0 (f ) -1 (Aϕ (u))χ S (u) + γ e r2(u-L) - r 1 r 2 e r1(u-L) du + x 0 (f ) -1 (Aϕ (u))χ S (u) + γ e r2(u-x) -e r1(u-x) du, (34)
where A > 0 is the unique solution of

D L1 0 (f ) -1 (Aϕ (u))λ (u) du -θ 1 = Aθ 2 (35)
and

r 1 = -v 1 - √ ∆ 2S , r 2 = -v 1 + √ ∆ 2S , ∆ = v 2 1 + 4Sα, K = v 1 + Sr 1 v + Sr 2 , θ 1 = 1 Sr 1 r 2 γ < 0, θ 2 = (v 1 + Sr 1 )(e r1L -e r2L ) > 0, ϕ (u) = e r1u -Ke r2u > 0, λ (u) = (-r 1 e r1u + r 2 e r2u )e -(r1+r2)L Sr 1 r 2 e -r2L -e -r1L
> 0.

Optimal size of the buffer zone

We now exploit the analytic solution for the study of the impact of the length of the buffer zone on the optimal fertilizer load. The first result is the following lemma.

Lemma 3 Assume that the inverse function (f ) -1 is positive. Any reduction of the zone devoted to spreading induces a higher optimal fertilizer load in the remaining cropland, that is ∂ L1 p * < 0.

Proof First, we use the implicit formulation of A with respect to L 1 in equation [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF] for computing its derivative with respect to L 1 . We get

dA dL 1 = -D L1 0 (f ) -1 (Aϕ(u))λ(u)du -θ 1 (f ) -1 (Aϕ(L 1 ))λ(L 1 ) D L1 0 (f ) -1 (Aϕ(u))λ(u)du -θ 1 L1 0 ϕ(u)λ(u) f ((f ) -1 (Aϕ(u))) du -θ 2 .
The latter quantity is positive since the functions D , ϕ, λ, (f ) -1 and -f are positively valued and θ 2 > 0. Then deriving [START_REF] Godart | Rapport Final -APR GICC 2002, Convention de recherche MEDD num 02-00019[END_REF] with respect to L 1 and using the latter result, we obtain

∂p * ∂L 1 (x, L 1 ) = dA dL 1 ϕ(x) f ((f ) -1 (Aϕ (x))) < 0.
The lemma is proved.

Lemma 3 means that the larger the buffer zone, the more intensive fertilizer loads are required in the optimal setting. Notice that the result does not hold when admitting the possibility of active buffer zones in the sense that fertilizer load p, thus here (f ) -1 , may be nonpositive. Thus the impact and the "better" determination of a buffer zone is not straightforward for the social planner who has to deal with an objective summing up the one of the farmer and the one of the freshwater operator.

For the water operator's contribution in the objective, the cost, which can be defined as a function J Wa of L 1 , is given by D(c (L)). In the optimal setting, it reads as follows, according to equation [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF]:

J Wa (L 1 ) = D (D ) -1 (Aθ 2 ) since c * (L) = (D ) -1 (Aθ 2 ). Using ∂ L1 c(L) = (dA/dL 1 )θ 2 /D ((D ) -1 (Aθ 2 )), we compute J Wa (L 1 ) = D (D ) -1 (Aθ 2 ) ∂c ∂L 1 (L) = Aθ 2 ∂c ∂L 1 (L) > 0.
We thus recover that when L 1 is increasing, the buffer zone shrinks and the pollution and cleaning cost are increasing.

For the farmer's contribution in the objective, the evaluation of the impact is more tricky as there is a size effect (the larger the buffer zone, the smaller the crop field) and a quantity effect (the larger the buffer zone, the higher the fertilizer loads). Denote J F (L 1 ) the benefit for the farmer. We have

J F (L 1 ) = L1 0 f (f ) -1 (Aϕ(x)) dx and J F (L 1 ) = f (f ) -1 (Aϕ(L 1 )) + dA dL 1 L1 0 Aϕ(x) 2 f ((f ) -1 (Aϕ(x)))
dx.

The first term depends on the agricultural technical function f . It may be positive if a Mitscherlich production function f (p) = α (1 -e -p ) (Kennedy [START_REF] Kennedy | Dynamic programming: applications to agriculture and natural resources[END_REF]) or a Godart function ( [START_REF] Godart | Rapport Final -APR GICC 2002, Convention de recherche MEDD num 02-00019[END_REF]) f (p) = β J -(β J -α J )e -p(x)τ J is used. The second term is negative. Thus, there is no generic monotony result for J F . Finally, we consider the question of the optimal size of the buffer zone by working now on the trade-off between the farmer's and the operator's contributions to the objective. Let Obj be the function describing the indirect value of the social welfare as a function of L 1 :

Obj(L 1 ) = J F (L 1 ) -J Wa (L 1 ) .
The optimal size of the buffer zone, if it exists, is characterized by the variations of the function Obj. Indeed, if Obj admits a global maximum at point

L * 1 ∈ [0, L] then the optimal buffer zone is [L * 1 , L]. Of course, it makes only sense if L * 1 < L.
Otherwise, the crop field extends up to the captation well. We now point out how the existence of such an optimal value is related to the parameters and data of the problem. To this aim, we characterize the critical points of Obj. According to the previous computation, we have the following result.

Lemma 4

The variations with respect to L 1 of the indirect value of the social welfare function are characterized by

∂Obj ∂L 1 = f (f ) -1 (Aϕ (L 1 )) + dA dL 1 L1 0 ϕ(x) 2 f ((f ) -1 (Aϕ(x))) dx - θ 2 2 D ((D ) -1 (Aθ 2 )) ≤0 if (f ) -1 is nonnegative .
For instance, according to the computations in the proof of Lemma 3, when (f ) -1 is nonnegative that is if we only consider positive fertilizer loads, a necessary, of course not sufficient, condition for the existence of the optimal value L * 1 is that L * 1 cancels the latter derivative. In the next subsection, a numerical experiment consisting in tuning the cost function D provides the two possible settings, existence and nonexistence of the optimal value L * 1 .

Numerical simulations

Our aim in this paragraph is to qualitatively illustrate the results obtained previously. We show that the existence of an optimal buffer zone is strongly related to the choice of the technological parameters (production function and cleaning technology). Moreover, for a given technology, the size of the controlled area matters.

The tests are realized with Maple. We use the following data:

α = 0.1, S = 0.1, ψ = 1, γ = 0.05, v 1 = 0.05, L = 10. ( 36
)
Notice that these are the data for the dimensionless problem. Going back to the dimensional form, considering for instance an aquifer of L = 10 km length and H = 1 m height, this means that L R = 10 3 , H R = 1 and = 10 -3 . Assuming moreover that S = 10 -5 m 2 s -1 = 0.86 m 2 day -1 then implies that T > O(10 5 ) days and v = 4 • 10 -4 m day -1 . This latter value for the velocity is small and is very advantageous for the farmer since the pollutant is pushed quite slowly toward the freshwater production well. Below we plot the (optimal) concentration c * and fertilizer load p * solving the dimensionless problem. For turning back to the dimensional scale, one has to know the characteristic value of the natural input flux γ

: indeed c R = O(5γ • 10 7 ) and p R = O(5γ • 10 2 ).
We introduce realistic crop and damage functions. Let

f (p) = L R H R P β J -(β J -α J )e -2•10 6 p R pτ J /(2 • 10 6 ), D(c) = 24Q × 0.2(c R c) 2 /50 2 .
The choice for the production function comes from Godart et al. [START_REF] Godart | Rapport Final -APR GICC 2002, Convention de recherche MEDD num 02-00019[END_REF] where the crop yield is depending on α J and β J which are respectively the yield value without nitrogen input and the asymptotic value when the input becomes important. The maximal yield is reached at a velocity τ J . Values of the parameters depend on the crop species. For example, for the wheat crop, β J = 11.7888, α J = 3.1824 and τ J = 5.0465 • 10 -3 . We express here the output in tons per m 2 and day. We assume here that the land is cultivated 200 days per year. The price of the crop species per ton is P . For the computations, we consider that for the wheat crop P = 200 e. Concerning the damage function, we have chosen data from Ledoux et al. [START_REF] Ledoux | Agriculture and groundwater nitrate contamination in the Seine basin. The STICS-MODCOU modelling chain[END_REF] where only the cost to reach the norm of drinking water (less than 50mg/L nitrogen in Europe) is taken into account (and evaluated at 0.2e per m 3 ). Quantity Q denotes the production flow rate of the well (m 3 per hour).

We are going to question the existence of an optimal buffer zone when tuning the production rate Q of the well and thus the cleaning cost. We begin by a choice greatly favouring the farmer, a very small value for Q, namely Q = 17 m 3 hour -1 . The cleaning cost is thus little weighted compared to the production benefit. We get the results represented in Figure 1 (computed with c R = O(10 -2 ) and p R = O(10 -6 )). The shape of the graphs of the optimal concentration c * does not really change with the size of the buffer zone. Indeed, we recover here the classical curve of a diffused and advected solute. However, the values of the optimal concentration c * are influenced by the size of the buffer zone: as stated in Lemma 3, for a positive fertilizer load, the longer the buffer zone, the higher the spreading (see Figure 1-bottom) thus the higher the optimal concentration c * in the fertilized area (and only in this area, see for instance the curve corresponding to L 1 = 2 in Figure 1-left). This phenomenon is less pronounced if L 1 is close to L. However, the result of Lemma 3 does not hold for L * 1 = 10 because then p * is negative. It explains the shape of the black curve for L * 1 = 10 in Figure 1-left and Figure 1-bottom.

The objective function is increasing on [0, L] (Figure 1-right), which means that the optimal policy consists here in having no buffer zone (that is L * 1 = 10). The latter results are of course completely induced by the very low production rate of the water well which deeply favours the profit of the farmer in the model. The picture is very different for more realistic (realistic for the profit of the water operator) production rates. Set for instance Q = 10 3 m 3 hour -1 . The corresponding optimal fertilizer load computed by the algorithm is represented in Figure 2. Of course, one notices that too important negative loads are meaningless because unreachable, thus the quotation marks around 'optimal' in the caption. But their interpretation is clear. For important water well production rates, there is an optimal length for the buffer zone (here larger than L -8 = 2), unless the passive buffer zone is replaced by an active one. If no buffer zone exists, that is here if L * 1 = 10, the result in Figure 2 outlines the need of intercropping for highly absorbing the soil nutrients (see Willey [START_REF] Willey | Resource use in intercropping systems[END_REF] or GEVES [START_REF] Geves | French Variety and Seed Study and Control Group, Service plants[END_REF]). This point is confirmed by the Figure 3 depicting the variations of the objective, Obj, with respect to L 1 for Q = 10 3 and Q = 10 5 m 3 hour -1 (of course for these two latter computations we corrected the algorithm so that negative fertiliser loads are no more allowed). 5 Proof of the theoretical results

Proof of Theorem 1

Theorem 1 was already proved in [START_REF] Augeraud-Véron | Optimal control for a groundwater pollution ruled by a convection-diffusion-reaction problem[END_REF]. Nevertheless we reproduce here its main steps for the completeness of the article and, mainly, because they are reused in the proof of Proposition 2 below. The proof is decomposed into several lemmas. It is written for the setting T < ∞. The case T = ∞ may be treated in the same with only some minor changes due to the definition of the objective J.

First, let us turn back to the existence result for the state equations system given in Proposition 1 and point out some regularity results. The reader may easily check that

∂ t c ∈ L 2 (0, T ; H -1 (Ω)). It follows, since moreover c ∈ L 2 (Ω× (0, T )), that c ∈ C([0, T ]; L 2 (Ω)).
In particular, we can write c -≤ c(t, x) for almost every x ∈ Ω, for all t ∈ (0, T ) and the existence of c(T, x) is ensured almost everywhere in Ω. Now, we begin by focusing on the uniqueness result in Theorem 1.

Lemma 5

If there exists a solution to the problem P E T , then it is unique.

Proof Notice first that the hydraulic head φ appearing in the state system is uniquely defined since its equation does not depend on c nor on p. Now assume that there exist two optimal solutions (p * 1 , c * 1 ) and (p * 2 , c * 2 ) to P E T . Denote

J opt := J(p * 1 , c * 1 ) = J(p * 2 , c * 2 ).
We compute J((p * 1 + p * 2 )/2, c ) where c is the concentration issued from the state system with (p * 1 + p * 2 )/2 as fertilizer load, i.e. c solves:

     Rψ∂ t c + v • ∇c -div(ψS(v)∇c ) = -r(c ) + (p * 1 + p * 2 )χ S /2 -gc + γ in Ω T , S(v)∇c • n = 0 on Γ 1 × (0, T ), c = c -on Γ 2 × (0, T ), (37)
the velocity v being given by ( 4)- [START_REF] Keeler | The optimal control of pollution[END_REF]. We obtain

J 1 2 (p * 1 + p * 2 ), c = T 0 S f (x, 1 2 p * 1 + 1 2 p * 2 ) dx - Ω D(x, c ) dx e -ρt dt -νe -ρT Ω D(x, c (T, x)) dx.
If p * 1 = p * 2 almost everywhere in Ω T , the strict concavity of J implies that

J 1 2 (p * 1 + p * 2 ), c > T 0 S 1 2 f (x, p * 1 ) + 1 2 f (x, p * 2 ) dx dt - T 0 Ω D(x, c ) dx e -ρt dt -νe -ρT Ω D(x, c (T, x)) dx. ( 38
)
Using the equation satisfied by c = (c * 1 + c * 2 )/2 and (37), it can be easily shown that max(c -c, 0) = 0 and thus -D(x, c) ≤ -D(x, c ) a.e. in Ω T . We thus have

J 1 2 (p * 1 + p * 2 ), c > T 0 S f (x, p * 1 ) 2 + f (x, p * 2 ) 2 dx - Ω D(x, c) dx e -ρt dt -νe -ρT Ω D(x, c(T, x)) dx ≥ T 0 S 1 2 f (x, p * 1 ) dx e -ρt dt + T 0 S 1 2 f (x, p * 2 ) dx e -ρt dt - T 0 Ω 1 2 D(x, c * 1 ) dx e -ρt dt - T 0 Ω 1 2 D(x, c * 2 ) dx e -ρt dt - ν 2 e -ρT Ω D(x, c * 1 (T, x)) dx - ν 2 e -ρT Ω D(x, c * 2 (T, x)) dx = 2 × 1 2 J opt , (39)
the last inequality being justified by the convexity property of D. To summarize, we have J( 1 2 (p * 1 + p * 2 ), c ) > J opt where c solves the state problem associated with (p * 1 + p * 2 )/2. It is a contradiction with the definition of the maximal value J opt . The uniqueness of the solution is proved.

We now deal with the existence part in the proof of Theorem 1.

Lemma 6 There exists a solution to the control problem P E T .

Proof For any given p ∈ E T , denote by c p the concentration solution given by Proposition 1 for the fertilizer load p. According to the assumptions imposed to f and D, the set {J(p, c p ); p ∈ E T } is bounded in R. Indeed, as proved in Augeraud et al. [START_REF] Augeraud-Véron | Optimal control for a groundwater pollution ruled by a convection-diffusion-reaction problem[END_REF], the solution c exhibited in Proposition 1 is bounded in L 2 (Ω T ) (and the bound only depends on the data of the problem). The supremum of the set {J(p, c p ); p ∈ E T } thus exists and we denote it by J * . We aim at proving that there exists (p * , c * = c p * ) such that

J * = J(p * , c * ),
i.e. that the supremum is a maximum. According to the definition of the supremum, there exists a sequence (p n ) n∈N ⊂ E T , associated with (c n = c pn ) n∈N , such that lim

n→∞ J(p n , c n ) = J * .
Since p n is uniformly bounded in L ∞ (Ω T ) by the real number p, there exists a subsequence7 weakly converging to some function p * in L 2 (Ω T ). With classical energy estimates, one also checks that c n is uniformly bounded in L 2 (0, T ; H 1 (Ω)) ∩ H 1 (0, T ; H -1 (Ω)). Thus there exists a subsequence of c n weakly converging in the space L 2 (0, T ; H 1 (Ω)) to some function c. The Aubin's compactness argument (Simon [35]) shows that c n actually converges strongly to c in L 2 (Ω T ). Passing to the limit n → ∞ in the state system defining c n , we check that c is a concentration solution of the state system associated with p * . The uniqueness result in Proposition 1 allows the use of the notation c = c p * = c * . Moreover, c -≤ c * (t, x) almost everywhere in Ω, for all t ∈ [0, T ]. We now prove that

lim n→∞ J(p n , c n ) = J(p * , c * ).
Since c n converges to c * in L 2 (Ω T ), hence almost everywhere in Ω for all t ∈ [0, T ], and since function c → D(•, c) is continuous, Lebesgue's dominated convergence theorem implies that

J * = lim n→∞ T 0 S f (x, p n ) dx e -ρt dt - T 0 Ω D(x, c * ) dx e -ρt dt -νe -ρT Ω D(x, c * (T, x)) dx = lim n→∞ T 0 S f (x, p n ) dx e -ρt dt - T 0 Ω D(x, c * ) dx e -ρt dt -νe -ρT Ω D(x, c * (T, x)) dx.
Since f is concave, a weak lower semicontinuity theorem (see Tartar [START_REF] Tartar | Compensated compactness and applications to partial differential equations[END_REF]) implies that

T 0 S f (x, p * ) dx e -ρt dt ≥ lim n→∞ T 0 S f (x, p n ) dx e -ρt dt.
We infer from the two latter expressions that

J * ≤ T 0 S f (x, p * ) dx - Ω D(x, c * ) dx e -ρt dt -νe -ρT Ω D(x, c(T, x)) dx = J(p * , c * ).
But p * ∈ E T and c * = c p * . Thus, according to the definition of the supremum J * , we also have J(p * , c * ) ≤ J * . The equality J(p * , c * ) = J * is proven and (p * , c * ) is the desired optimal solution of Problem P E T .

Proof of Proposition 2 (thus of Theorem 2)

The first step consists in showing that ( 22)-( 24) is the effective hydrological state system. Let p be given in E. Let (c , φ ) solve ( 14)- [START_REF] Fonseca | 3D-2D asymptotic analysis of an optimal design problem for thin films[END_REF]. We aim at proving that (c , φ ) converges to a solution of ( 22)-( 24) as → 0. Using respectively c and φ as test functions in ( 14)-( 15) and integrating by parts, one easily checks that the following uniform estimates hold true:

∂ x φ L ∞ (0,T ;L 2 (Ω)) + -1 ∂ z φ L ∞ (0,T ;L 2 (Ω)) ≤ C, ∂ x c L 2 (Ω T ) + -1 ∂ z c L 2 (Ω T ) + √ α + g c L 2 (Ω T ) ≤ C.
Notice that the second part of Assumption 6 allows also uniform estimates for v in (L ∞ (Ω T )) 2 . Since we assume that the space average of the hydraulic head8 is uniformly bounded with respect to , that is Ω φ dxdz ≤ C, we infer from the first estimate above and from the Poincaré-Wirtinger inequality that φ is uniformly bounded in L ∞ (0, T ; L 2 (Ω)). Then, there exist subsequences, not relabelled for convenience, and limit functions φ ∈ L ∞ (0, T ;

H 1 (Ω)), c ∈ L 2 (0, T ; H 1 (Ω)), φ ∈ L 2 (Ω T ) and c ∈ L 2 (Ω T ) such that φ φ and c c weakly in L 2 (0, T ; H 1 (Ω)), -1 ∂ z φ φ , -1 ∂ z c c weakly in L 2 (Ω T ).
In particular the limit functions do not depend on the vertical variable:

∂ z φ = ∂ z c = 0 a.e. in Ω T .
This point is now exploited for tackling the singular structure of the microscopic problem. Let ϕ ∈ H 1 ((0, T ) × (0, L)) with ϕ |t=T = 0. The variational formulation of ( 14)-( 17) with ϕ as a test function is

L 0 H(x) h(x) κ x ∂ x φ dz ∂ x ϕ dx - H(0) h(0) v 1 ϕ |x=0 dz + H(L) h(L) v 1 ϕ |x=L dz = L 0 H(x) h(x)
g dz ϕ dx a.e. in (0, T ),

-

Ω T Rψc ∂ t ϕ dxdzdt - Ω Rψc 0 ϕ |t=0 dxdz - T 0 L 0 H(x) h(x) (κ x ∂ x φ ∂ x c -κ z -1 ∂ z φ -1 ∂ z c ) dz ϕ dxdt + T 0 L 0 H(x) h(x) ψS x ( v )∂ x c dz ∂ x ϕ dxdt = - T 0 L 0 H(x) h(x) αc + gc -pχ S -γ dz ϕ dxdt. (41)
Passing to the limit → 0, we recover ( 22)- [START_REF] Louison | Modeling plant nutrient uptake: mathematical analysis and optimal control[END_REF]. Notice that for computing the limit of the nonlinear terms in the second integral equation, we have to use a kind of strong convergence result for the velocity v . Multiply the first equation in ( 15) by (φ -φ) and integrate by parts over Ω. We obtain

Ω κ x ∂ x φ ∂ x (φ -φ) dxdz + -2 Ω κ z ∂ z φ ∂ z (φ -φ) dxdz - H(0) h(0) v 1 (φ -φ) |x=0 dz + H(L) h(L) v 1 (φ -φ) |x=L dz = Ω g(φ -φ) dxdz.
We rewrite the latter relation in the following form

Ω κ x |∂ x (φ -φ)| 2 dxdz + -2 Ω κ z |∂ z (φ -φ)| 2 dxdz - H(0) h(0) v 1 (φ -φ) |x=0 dz + H(L) h(L) v 1 (φ -φ) |x=L dz = Ω g(φ -φ) dxdz - Ω κ x ∂ x φ ∂ x (φ -φ) dxdz
where

- Ω κ x ∂ x φ ∂ x (φ -φ) dxdz = - Ω κ∂ x φ ∂ x κ-1 κ x (φ -φ) dxdz + Ω κ∂ x φ (φ -φ)∂ x κ-1 κ x dxdz = Ω ḡ κ-1 κ x (φ -φ) dxdz - H(0) h(0) v 1 κ-1 κ x (φ -φ) |x=0 dz + H(L) h(L) v 1 κ-1 κ x (φ -φ) |x=L dz + Ω κ∂ x φ (φ -φ)∂ x κ-1 κ x dxdz.
We thus have

Ω κ x |∂ x (φ -φ)| 2 dxdz + -2 Ω κ z |∂ z (φ -φ)| 2 dxdz = Ω g(φ -φ) dxdz + Ω ḡ κ-1 κ x (φ -φ) dxdz - H(0) h(0) v 1 (Id -κ-1 κ x )(φ -φ) |x=0 dz + H(L) h(L) v 1 (Id -κ-1 κ x )(φ -φ) |x=L dz + Ω κ∂ x φ (φ -φ)∂ x κ-1 κ x dxdz
where the right-hand side terms tend to zero with since φ -φ 0 weakly in H 1 (Ω). It follows that we have the following (strong) convergence results:

∂ x φ → ∂ x φ and -1 ∂ z φ → 0 in L 2 (Ω).
The latter result allows to pass to the limit in the singular terms of the concentration equation.

Before studying the limit behavior of the whole optimization problem, let us check that the results in Proposition 3 hold true. Since the counterpart of Proposition 1 is established by the convergence results above, the proof of Proposition 3 follows the lines of the one of Theorem 1 given in Subsection 5.1. Notice that the proof of the necessary compactness results is here much simpler since the problem is now only space-dependent. The compact embedding of H 1 (0, L) into L 2 (0, L) is thus sufficient for our purpose, with no need of Aubin's type arguments. For the uniqueness result, the proof follows the one for Lemma 5. However there is a difference since the limit process above shows the existence of a solution (c, φ) of problem ( 22)-( 24) but does not ensure its uniqueness. Turning back to the notations in the proof of Lemma 5,it

means that (p * 1 , c * 1 ) = (p * 2 , c * 2 ) does not imply p * 1 = p * 2 .
Thus the strict inequality is no more ensured in [START_REF] Stefanelli | The Brezis-Ekeland principle for doubly nonlinear equations[END_REF]. Nevertheless, we recover the strict inequality in the fundamental result (39) since we additionally assume here that D is strictly convex and because (p

* 1 , c * 1 ) = (p * 2 , c * 2 ) ⇔ p 1 = p * 2 or c * 1 = c * 2 .
The proof of Proposition 3 is complete.

We now perform the asymptotic analysis of the optimal control problem as → 0. First we note that, thanks to the mean value theorem, lim

→0 J (p, c) = J eff (p, c) for any (p, c) ∈ E T × L 2 (0, T ; C( Ω)).
Let (p * , c * ) be the optimal solution such that J (p * , c * ) = max p∈E T J (p, c) where J is defined in [START_REF] Kogut | Optimal Control Problems for Partial Differential Equations on Reticulated Domains, Approximation and Asymptotic Analysis[END_REF] and c is the solution of ( 14)- [START_REF] Fonseca | 3D-2D asymptotic analysis of an optimal design problem for thin films[END_REF]. Since p * ∈ E T and c * is a solution of ( 14)-( 17) with p = p * , one easily checks that we have sufficient boundedness result for ensuring the existence of a subsequence, not relabeled for convenience, and of a limit function (p , c ) such that

p * p weakly in L 2 (Ω T ), c * c weakly in L 2 (0, T ; H 1 (Ω))
and such that c is a solution of ( 22)-( 24) with p = p (follow the lines above for the study of the limit behavior of the hydrodynamical system). If (p * , c * ) is the solution of P eff given by Proposition 3, we thus have

J eff (p , c ) ≤ J eff (p * , c * ). (42)
Since (p , c ) does not depend on the vertical variable z, one may also consider that c is the solution associated to p through the scaled state system ( 14)-( 17), so that we also have J (p , c ) ≤ J (p * , c * ) and, passing to the upper limit lim sup

→0 J (p , c ) = J eff (p , c ) ≤ lim sup →0 J (p * , c * ). ( 43
)
We now aim at proving that lim sup J (p * , c * ) ≤ J eff (p , c ). It remains to pass to the limit in

Of course we have lim

T 0 D meas(C ) -1 C c * dxdz e -ρt dt. Let b * (t, x) = C ∩(h(x),H(x)) c * (t, x, z) dz, b (t, x) = C ∩(h(x),H(x)) c (t, x, z) dz.
Recall that, since we deal here with a dimensionless model, the height C ∩ (h, H) of the cylinder C does not depend on . Hence, we infer from the weak convergence of c * to c in L 2 (0, T ;

H 1 (Ω)) that b * b weakly in L 2 (0, T ; H 1 (0, L)). Now define the functions G , G 0 : H 1 (0, L) → (-∞, ∞] by G (u) = D -1 C ∩(0,L) u(x) dx , G 0 (u) = D u(L w ) .
Due to the properties of D in R and to the properties of the integral, these are two convex, proper and lower semicontinuous functions. Moreover (t, u) ∈ (0, T ) × H 1 (0, L) → e -ρt G (u), e -ρt G 0 (u) are normal convex integrands (see Rockafellar [START_REF] Rockafellar | Integrals which are convex functionals[END_REF]). Using the compact embedding of H 1 (0, L) in C([0, L]), we check that they also satisfy

G 0 (u) = lim →0 G (u ) ∀u u weakly in H 1 (0, L)
for any u ∈ H 1 (0, L). The Corollary 4.4 in Stefanelli [START_REF] Stefanelli | The Brezis-Ekeland principle for doubly nonlinear equations[END_REF] is thus applicable. We conclude that In the same way, since (p * , c * ) does not depend on the vertical variable z, one may also consider that c * is the solution associated to p * through the scaled state system ( 14)-( 17), so that we also have J (p * , c * ) ≤ J (p * , c * ). Passing to the upper limit and using (47), we get

J eff (p * , c * ) ≤ lim sup →0 J (p * , c * ) ≤ J eff (p , c ).
Combining the latter result with (42), we conclude that

J eff (p * , c * ) = J eff (p , c ).
Thanks to the uniqueness result proved below for P eff , we necessary have

(p * , c * ) = (p , c )
and the whole sequence (p * , c * ) converges to (p * , c * ). This ends the proof of Proposition 2.

Remark 2 For proving the uniqueness result of Proposition 3, we used the assumption ∂ cc D > 0 because we did no more (compared with the proof of Theorem 1) show the uniqueness of the solution to the state system ( 22)- [START_REF] Louison | Modeling plant nutrient uptake: mathematical analysis and optimal control[END_REF]. Nevertheless, notice that additional assumptions may be specified for recovering the uniqueness result. It is sufficient to ensure that the continuous bilinear form a 0 defined in H 1 (0, L) × H 1 (0, L) by

a 0 (c, ϕ) = L 0 ψS x (v)∂ x c • ∂ x ϕ dx + α L 0 (H -h)cϕ dx + L 0 v∂ x c ϕ dx
is uniformly elliptic. Indeed, such a property would allow the use of the Lax-Milgram theorem for proving the well-posedness of ( 22), [START_REF] Louison | Modeling plant nutrient uptake: mathematical analysis and optimal control[END_REF]. We compute

a 0 (c, c) = L 0 ψS x (v)|∂ x c| 2 dx + α L 0 (H -h)|c| 2 dx + L 0 v∂ x c c dx ≥ S m ψ -inf (0,L) (H -h) L 0 |∂ x c| 2 dx + α L 0 (H -h)|c| 2 dx - L 0 |v∂ x c c| dx.
Thanks to the Cauchy-Schwarz and Young inequalities, we write

L 0 |v∂ x c c| dx ≤ εS m ψ -inf (0,L) (H -h) L 0 |∂ x c| 2 dx + v 2 L ∞ (0,L) 4εS m ψ -inf (0,L) (H -h) L 0 |c| 2 dx
for any ε > 0. Turning back to a 0 (c, c), we get for any ε > 0:

a 0 (c, c) ≥ S m ψ -inf (0,L) (H -h)(1 -ε) L 0 |∂ x c| 2 dx + α inf (0,L) (H -h) - v 2 L ∞ ((0,L)) 4εS m ψ -inf (0,L) (H -h) L 0 |c| 2 dx.
Thus there exists a -> 0 such that a

0 (c, c) ≥ a -c 2 H 1 (0,L) if the data satisfy 4S m ψ -α( inf (0,L) (H -h)) 2 > v 2 L ∞ (0,L) . (48)

Proof of lemma 1

Let µ be the adjoint variable related to state variable c. The Lagrangian function of the problem is

L(c, p, µ) = J(p, c) + L 0 vc -ψSc + αc -pχ [0,L1] (x) -γ µ dx. ( 49
)
Integrating by parts and using the Neumann boundary conditions yield to

L(c, p, µ) = L1 0 f (p(x)) dx -D(c(L)) + L 0 (-vcµ -ψSµ c + αcµ -pχ [0,L1] µ -γµ) dx + v(L)µ(L)c(L) -v(0)µ(0)c(0)
+ ψSc(L)µ (L) -ψSc(0)µ (0).

The following expression is obtained by using Taylor's first order formula: Cancelling the Lagrangian variations with respect to the control p and to the state variable c, provides us the announced necessary conditions.

Proof of lemma 2

Recall that our aim is the resolution of the following two equations in [0, L] (we drop the stars for the sake of the simplicity of the notations):

-Sc + vc + αc -pχ S -γ = 0, ψSµ + vµ -αµ = 0, (50

)
where v is given by v = v 1 and p is such that f (pχ S ) = µχ S , completed by the boundary conditions:

D (c(L)) = ψSµ (L) + vµ(L), c (0) = c (L) = 0, ψSµ (0) + vµ(0) = 0.

Any solution for µ is in the form µ(x, A, L 1 ) = Ae r1x + Be r2x where r 1 = (-v 1 -√ ∆)/2S, r 2 = (-v 1 + √ ∆)/2S and ∆ = v 2 1 + 4Sα, A and B being real numbers. Thanks to the boundary condition satisfied by µ at x = 0, the parameter B is computed:

B = -AK where K = v 1 - √ ∆ v 1 + √ ∆ = v 1 + Sr 1 v 1 + Sr 2 < 0.
Hence, µ(x, A, L 1 ) = Aϕ(x), ϕ(x) = e r1x -Ke r2x .

Notice that K < 0 and r 1 < 0 < r 2 , thus e r1x -Ke r2x > 0 for every

Conclusion

We provide a time and space dependent model of groundwater pollution control. We prove in a general setting existence and uniqueness of the optimal solution. We also rigorously derive and justify an effective optimal control problem that may be substituted to the original one for taking into account the scale heterogeneity of the problem. This process allows to use a simpler spatial and static model to study analytically the problem, still based on a realistic model of pollutant displacement in groundwater, that encompass the spatial concern of such problems. The question of implementing buffer zones is tackled in order to increase the social welfare. We show that the optimal buffer zone does not correspond necessarily to the one that would prevail. Numerical simulations for the full 3D time dependent model are postponed to a forthcoming paper.

Fig. 1

 1 Fig. 1 Optimal concentration c * (left), objective with respect to L 1 (right) and optimal fertilizer load p * (bottom) if the water production flow is Q = 17 m 3 hour -1 .

Fig. 2 "

 2 Fig. 2 "Optimal" fertilizer load p * computed by the algorithm if the water production flow is Q = 10 3 m 3 hour -1 .

Fig. 3

 3 Fig. 3 Variations of the objective, Obj, with respect to L 1 , if the water production flow is Q = 10 3 m 3 hour -1 (left) and Q = 10 5 m 3 hour -1 (right). Here only nonnegative fertilizer loads are allowed.

f

  (x, p * ) dxdz e -ρt dt ≤ T 0 Ω f (x, p ) dxdz e -ρt dt. (45)

T 0 G 0 0 G 0 D meas(C ) - 1 Cc

 00001 (b )e -ρt dt ≤ lim inf →0 T (b * )e -ρt dt, that is T 0 D(c (L w ))e -ρt dt ≤ lim inf →0 T * dxdz e -ρt dt. (46) We infer from (44)-(46) that lim sup →0 J (p * , c * ) ≤ J eff (p , c ). (47) Now it follows from (43) and (47) that lim sup J (p * , c * ) = J eff (p , c ).

  δL(c, p, µ) = L1 0 (f (p(c)) -µ) δp dx -L 0 (vµ + ψSµ -αµ)δc dx -v(0)µ(0) + ψSµ (0) δc(0) + v(L)µ(L) + ψSµ (L) -D (c(L)) δc(L).

2 with θ 2 = (v 1 + 0 (f ) - 1 0 (f ) - 1 r 2 .

 22101012 x ∈ [0, L]. Moreover, the boundary condition D (c(L)) = ψSµ (L) + vµ(L) implies that D (c(L, A, L 1 )) = A (Sr 1 + v 1 )e r1L -K(Sr 2 + v 1 )e r2L = Aθ Sr 1 )(e r1L -e r2L ) ≥ 0. Since D is strictly increasing, quantity D (c(L, A, L 1 )) is positive. With e r1L -e r2L < 0, v 1 + Sr 1 = -Sr 2 < 0, we getA > 0. (51)We now focus on c. Since A > 0, (f ) -1 Aϕ(x) is well-defined for any x ∈ [0, L 1 ] and the equation for c also reads-Sc + v 1 c + αc -(f ) -1 Aϕ(x) χ [0,L1] (x) + γ = 0 so that c(x, A, L 1 ) = K 1 e -r2x + K 2 e -r1x + 1 S(r 2 -r 1 ) x 0 (f ) -1 Aϕ(u) χ [0,L1] (u) + γ (e r2(u-x) -e r1(u-x) ) du,(52)where K 1 , K 2 are real numbers that we now compute using the Neumann boundary conditions. According to (52), the derivative of c is∂c ∂x (x, A, L 1 ) = -K 1 r 2 e -r2x -K 2 r 1 e -r1x + 1 S(r 2 -r 1 ) × x Aϕ(u) χ [0,L1] (u) + γ (-r 2 e r2(u-x) + r 1 e r1(u-x) ) du (53)so that the Neumann boundary conditions specified for c at x = 0 and x = L respectively lead to K 2 = -K 1 /r 1 andK 1 = -L 0 ((f ) -1 Aϕ(u) χ [0,L1] (u) + γ) e r2(u-L) -r1 r2 e r1(u-L) du (e -r2L -e -r1L )S(r 2 -r 1 ).Inserting the values of K 1 and K 2 in (52), we obtainc(x, A, L 1 ) = 1 S(r 2 -r 1 ) x 0 (f ) -1 Aϕ(u) χ [0,L1] (u) + γ × e r2(u-x) -e r1(u-x) du -e -r2x -(r 2 /r 1 )e -r1x e -r2L -e -r1L × L Aϕ(u) χ [0,L1] (u) + γ e r2(u-L) -r 1 r 2 e r1(u-L) du . (54)For closing the resolution, it remains to characterize A. We turn back to the second boundary condition satisfied by µ, at x = L. It means thatD c(L, A, L 1 ) = Aθ 2 . (55)Now prove that A > 0 satisfying (55) exists and is unique. We have already shown in (51) that if A exists, then A > 0. Define the function η L1 in R * + by η L1 (a) = -aθ 2 + D (c(L, a, L 1 )).Equation (55) reads η L1 (A) = 0. For the study of the monotony ofη L1 we compute η L1 (a) = -θ 2 + D (c(L, a, L 1 ))∂ a c(L, a, L 1 ). -1 aϕ(u) λ(u) du -θ 1 , λ(u) = (r 2 e r2u -r 1 e r1u ) e -(r1+r2)L Sr 1 r 2 (e -r2L -e -r1L ) , θ 1 = γ Sr 1 Since ∂ a θ(u, a, L 1 ) < 0, we claim that ∂c ∂a (L, a, L 1 ) = e -(r1+r2)L S r 1 r 2 <0 e -r2L -e -r1L <0 × L 0 ∂((f ) -1 aϕ(u) χ [0,L1] (u) + γ) ∂a (-r 1 e r1u + r 2 e r2u ) >0 du < 0,(57)because, since e r1x -Ke r2x > 0, x ∈ [0, L], and since f is strictly concave, we have∂ (f ) -1 aϕ(x) χ [0,L1] (x) + γ ∂a = ϕ(x)((f ) -1 )aϕ(x) χ [0,L1] (x) < 0. Using (57) with the strict concavity of D in (56), we get η L1 (a) = -θ 2 + D (c(L, a, L 1 ))∂ a c(L, a, L 1 ) < 0. Hence, the function η L1 is strictly decreasing and it admits an inverse function. This inverse function and (55) allow to define A by A = η -1 L1 (0) provided that 0 belongs to the image of the function η L1 . It only remains to check that this final point holds true. Yet, lim a→0 η L1 (a) = lim a→0 D (c(L, a, L 1 )) > 0 and lim a→+∞ η L1 (a) = -∞ because c(L, a, L 1 ) → 0 as a → +∞, then D (c(L, a, L 1 )) → D (0) due to the continuity of D . According to the intermediate value theorem, 0 belongs to the image of the continuous function η L1 and A is uniquely defined by A = η -1 L1 (0).

The continuity can be relaxed by assuming that f is a bounded upper semi-continuous function.

The work[START_REF] Augeraud-Véron | Optimal control for a groundwater pollution ruled by a convection-diffusion-reaction problem[END_REF] only deals with the case p = 0 in[START_REF]L'eau et les milieux aquatiques -Chiffres clés -Edition[END_REF]. Nevertheless, a slight modification of the proof of the maximum principle for the concentration c gives Proposition 1.

We use the superscript diamond for the dimensional form of the equations.

The better regularity of c * , compared to the one stated in Proposition 3, comes obviously from the additional assumptions of the present subsection.

The terminology 'dual problem' is the one used in Benosman et al.[START_REF] Benosman | Optimization of cytostatic leukemia therapy in an advection-reaction-diffusion model[END_REF] andBarbu and Iannelli [31].

All along the paper, we denote the subsequences the same for the sake of the simplicity.

Bear in mind that the elliptic equation completed by a Neumann boundary condition defining φ may be completed by a prescription of the mean value of φ for ensuring the uniqueness of its solution.