# A numerical scheme for the optimal control of groundwater pollution

Catherine Choquet, Moussa Mory Diédhiou, Houssein Nasser El Dine, Christophe Saint-Jean

► To cite this version:

Catherine Choquet, Moussa Mory Diédhiou, Houssein Nasser El Dine, Christophe Saint-Jean. A numerical scheme for the optimal control of groundwater pollution. Numerical Methods for Partial Differential Equations, 2023, 39 (4), pp.3037-3063. 10.1002/num.22999. hal-04202875

HAL Id: hal-04202875

https://hal.science/hal-04202875v1

Submitted on 30 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
## RESEARCH ARTICLE

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

# A numerical scheme for the optimal control of groundwater pollution

Catherine Choquet<fnref n="1" /> ID |

Moussa Mory Diédhiou<fnref n="2" /> |

Houssein Nasser El Dine<fnref n="3" /> |

Christophe Saint-Jean<fnref n="1" />

**Correspondence**

Catherine Choquet, La Rochelle Université, MIA, La Rochelle, France.

Email: cchoquet@univ-lr.fr

**Funding information**

New Aquitaine Region ECONAT

**Abstract**

A scheme is proposed for the numerical study of optimal control problems related to groundwater quality management in the case of nonpoint-source pollution. The corresponding state equation is a system of coupled non-linear partial differential equations. The approximation is given by a finite volume scheme based on a two-point flux approximation with upwind mobilities embedded in an iterative fixed point approximation. The analysis of the convergence of the scheme allows to establish existence and uniqueness results under reasonable assumptions. Numerical illustrations of the performance of the algorithm are given in realistic situations.

**KEYWORDS**

convergence analysis, finite volume scheme, nonlinearly coupled problems, numerical simulations, optimal control, parabolic and elliptic PDEs

## 1 | INTRODUCTION

The water resource is today threatened by different kinds of pollution linked to human activity, most of them coming from economically necessary production activities. A trade-off between the benefits of production and the decontamination costs must therefore be sought. Such a question enters in the framework of the optimal control theory.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2023 The Authors. *Numerical Methods for Partial Differential Equations* published by Wiley Periodicals LLC

<footnote n="1">¹La Rochelle Université, MIA, La Rochelle, France</footnote>

<footnote n="2">²Université Cheikh Anta Diop, LMA, Dakar, Fann, Senegal</footnote>

<footnote n="3">³Al Maaref University, Mathematics, Beirut, Lebanon</footnote>
Modeling the spread of a miscible pollutant in groundwater is equivalent to consider the flow of incompressible miscible fluids in porous media. Both convection, diffusion and reaction phenomena should be taken into account in order to be more realistic, leading to a system of coupled nonlinear partial differential equations. The corresponding model will play the role of the state equation constraining the optimization problem. Although there is a large literature considering the question of optimal policy for managing pollution, most works neglect important components of the problem. In particular many studies reduce the hydrogeological model to ordinary differential equations without taking into account any space dependency. Thus the delay between the application of any policy and its effects, induced by the small flow speed in aquifers is either neglected or modeled by a given fixed delay (e.g., [2, 6, 8, 20]). The importance of taking a spatial dimension into account was however underlined in [7]. The space dependency is fully included in the present article.

We treat the optimal control problem from the point of view of numerical analysis. This allows us both to show the relevance of a scheme and to prove a result of theoretical analysis. In summary, we propose a finite volume architecture embedded in an iterative fixed point approximation to handle the full 3D optimal control problem. The choice of the finite volume approximation is relevant because, by construction, finite volume schemes are locally conservative. This property ensures that the overall iterative process quickly converges to an accurate description of the optimal control. Also the scheme has been built to deal with the strong parameter and scale contrasts induced by concrete applications.

There is of course a substantial literature on the numerical approximation of equation-constrained optimal control problems. A major issue in this context is that the control approximations may be highly oscillating: see [18] (resp. [14]) and the references therein for finite differences schemes (resp. for finite element schemes), both dealing with linear elliptic problems. As emphasized in [18], it seems that the oscillations mainly occur when using an approach in the form “discretize then optimize”. Here we choose the opposite approach, optimize (that is derive the adjoint problem) then discretize. Unlike the work [14, 18], we do not compute precise error estimates but simply establish the convergence of the scheme. Indeed, we have chosen to keep reasonable assumptions and we have poor regularity results for the weak solutions we treat (compared to the $H^2$ regularity in [14] or even the $C^4$ or $C^2$ assumptions in [18]). Another important remark explains why we do not attempt to improve our error estimates. In general, the performance of our algorithm will be doubly constrained: by those of the fixed point algorithm and by the structure of the two coupled advection-diffusion equations, in general with dominant advection, the advection in the adjoint problem being equal to the negative advection of the direct problem. The latter structure is induced by a change in the time variable that we make for handling the anti-diffusive form of the adjoint equation. In our nonlinear parabolic setting, we thus find the difficulty already highlighted in [14] for the elliptic linear setting. The behavior of the numerical method applied to the coupled PDEs is worse than the one of the method applied to the governing PDE only. This structural difficulty is unlikely to be controllable, especially since we want an approach general enough to simulate multiple nonlinear phenomena.

The article is organized as follows. The optimal control problem is introduced in Section 2. Section 3 is devoted to the description of our numerical scheme and to a sketch of the main results of the article. These results are then properly stated and proved in Sections 4–6 through the numerical analysis of the scheme. Notice that we have chosen to detail the proofs to emphasize the ease with which the scheme can be adapted to other optimal control problems. Finally, numerical experiments are performed in Section 7.
## 2 | STATEMENT OF THE PROBLEM

### 2.1 | The optimal control problem

Let $\Omega \subset \mathbb{R}^d$, $d \in \mathbb{N}$, $d \le 3$, be a bounded domain, representing the area affected by the pollution and containing the groundwater collection wells. We assume that the boundary $\partial\Omega$ of $\Omega$ is such that $\partial\Omega \in \mathcal{C}^1$. In the case of the 3D example treated in the last section, the boundary of $\Omega$ is divided into six subsets for taking into account different physical boundary conditions, namely $\partial\Omega = \cup_{i=1}^6 \Gamma_i$ where $\Gamma_1, \Gamma_2, \Gamma_3, \Gamma_4$ are the laterals faces and respectively $\Gamma_5, \Gamma_6$ are the top, bottom face. For the theoretical analysis however, it is sufficient to consider a division of $\partial\Omega$ into two components, $\Gamma_0$ and $\Gamma$. Moreover $\Gamma_d$ is a (possibly very small) part of $\Gamma_0$ where we use a Dirichlet boundary condition for the hydraulic head. We assume that $\Gamma_d$ is of nonzero measure for ensuring the validity of a Poincaré type inequality which facilitates the treatment of the elliptic equation of the model. The time horizon is denoted by $T$, $0 < T < \infty$. We set $\Omega_T = \Omega \times ]0, T[$. We consider the following control problem with a single control $p$: given the problem

$$ \mathrm{div}(v) = \theta, \quad v = -\kappa \nabla \phi \quad \text{in } \Omega_T, \qquad (1) $$

$$ R\psi \partial_t c + \mathrm{div}(v_c) = -r(c) + \chi_{\mathcal{C}} p, \quad v_c = -S(v)\psi \nabla c + cv \quad \text{in } \Omega_T, \qquad (2) $$

$$ c|_{t=0} = c_0 \quad \text{in } \Omega, \qquad (3) $$

$$ \begin{cases} v_c \cdot n = 0 & \text{on } \Gamma_0 = \Gamma_2 \cup \Gamma_4 \cup \Gamma_5 \cup \Gamma_6; \quad v_c \cdot n = f_1 & \text{on } \Gamma = \Gamma_1 \cup \Gamma_3 \\ v \cdot n = 0 & \text{on } \Gamma_0 \setminus \Gamma_d; \quad v \cdot n = v_1 & \text{on } \Gamma = \Gamma_1 \cup \Gamma_3; \quad \phi = 0 & \text{on } \Gamma_d, \end{cases} \qquad (4) $$

find the scalar control function $p$ which maximizes the functional

$$ J(p) = \int_{0}^{T} \left( \int_{\mathcal{L}} f(t, x, p(t, x))dx - \int_{\Omega} D(x, c(t, x))dx \right) e^{-\rho t} dt \\ - \zeta \int_{\Omega} D(x, c(t, x))dx \ e^{-\rho T}, \qquad (5) $$

subject to (2)–(4) and to the constraint

$$ 0 \le p \le \bar{p} \quad \text{in } \Omega_T. \qquad (6) $$

The state equations correspond to a hydrogeological model for the displacement of the pollution in the underground. They are given in the system (1)–(2) modeling the displacement of a miscible pollutant of concentration $c$ in the porous medium $\Omega$ under the coupled effect of convection, diffusion and reactions (see e.g., [4]). The velocity of the flow $v$ is related to the hydraulic head $\phi$ by the Darcy law for an incompressible fluid (1) where $\kappa$ is the mobility of the fluid in the soil and $\theta = q^l - q^p$ is a source term taking into account the recharge of the aquifer ($q^l$) and the water withdrawal in production wells ($q^p$). In (2), the dispersion effects due to the heterogeneity of the velocity field at the microscopic scale are given by the dispersion tensor $S(v) = S_m \mathbb{I}_{\mathbb{R}^d} + S_p(v)$, $\mathbb{I}_{\mathbb{R}^d}$ denoting the $d \times d$ identity matrix. For concision in the numerical analysis below, we retain only the diffusive component, thus setting $S_p \equiv 0$. Nevertheless the work extents to a general dispersive model under classical assumptions (see e.g., [1]). The porosity of the soil is given by the function $\psi$. The kinetic behavior of the pollutant is described by two terms, namely the retardation factor $R$ and the function $r(c)$, for taking into account reactions that may take place on different time scales. For the Neumann boundary conditions we denote by $n$ the outward normal to the boundary.
The pollution load is characterized by the rate of application $p$ and its location in a subset $\mathcal{L} \subset \Omega$, $\chi_{\mathcal{L}}$ being the characteristic function of $\mathcal{L}$. The function $p$ is actually the control of the problem. Note that this function depends on both space and time. Thus, in particular, it contains two pieces of information: the amount of pollution entering the soil and its location within $\mathcal{L}$, which can vary over time. The control is subject to the constraint (6) meaning from a mathematical point of view that it is bounded. From an application point of view, $\bar{p} > 0$ is the maximal pollution load (for instance the maximal fertilizer load that can be applied on a field). The non-negativity means that no remediation occurs.

The optimal control problem is in the Bolza form. We consider a central planner objective given by (5), following for instance the standard setting [20]. The function $D$ correspond to the cleaning costs depending on the position of the production wells and on the pollutant concentration in the groundwater $c$. In particular, function $D$ may contain the characteristic function of a subset of $\Omega$ for evaluating the costs only in a precise part of the domain of interest, for instance the water production wells. This allows to modulate the importance we give to environmental concerns in the model. The pollution is assumed to be linked with some economical activity. The corresponding benefits are described by the function $f$ and have to balance somehow the cleaning costs. The real numbers $\rho, 0 < \rho < 1$, and $\zeta \ge 0$ are respectively the social discount rate and a weight for the terminal costs.

## 2.2 | Mathematical framework

The state equations being a partial differential equations problem, we are naturally going to use the usual functional spaces associated to the mathematical analysis of PDEs. The set of admissible controls thus reads

$$E = \{p \in L^2(\Omega_T); 0 \le p(t, x) \le \bar{p} \text{ a.e. in } \Omega_T\}$$

and our problem is to find $p^*$ such that

$$p^* = \arg\max_{p \in E} J(p) \quad \text{subject to state constraints (1) - (4).} \qquad (7)$$

As already mentioned, we plan to use the optimality conditions of the first order instead of a direct approximation of (7). It is worth to mention that the characterization via any equivalent approach linked with the Pontryagin principle gives only necessary conditions for the optimal trajectory and optimal open-loop control. However, this approach makes sense because we can rely on an existence result for (7). Assume first:

(H0) For the sake of the simplicity in the notations, we treat the case where

$$\kappa = k \mathbb{I}_{\mathbb{R}^d}, \quad S(v) = S_m \mathbb{I}_{\mathbb{R}^d},$$

where $k$ and $S_m$ are given positive real numbers. Likewise the porosity $\psi$ and the factor $R$ are assumed to be positive real numbers. The initial data $c_0$ belongs to $L^2(\Omega)$.

(H1) The source term $\theta$ is a bounded function, $|\theta(x)| \le \theta^*$ a.e. in $\Omega$, $\theta^* \in \mathbb{R}^2$.

(H2) The function $f_1$ belongs to $L^2(\Gamma_1 \cup \Gamma_3)$ and $v_1$ is a positive real number.

(H3) The reaction kinetics is modelized by a continuous function $r$.

(H4) The benefit function $f$ is a continuous and concave w.r.t. $p$ function.

(H5) The function $D$ is a bounded, continuous and concave w.r.t. $c$ function.

Then the following result holds true (see lemma 5.4 in [1]).

**Proposition 2.1** (Assume (H1)–(H5)). There exists a solution $p^*$ of Problem (7).
The characterization of $p^*$ using the first order necessary optimality conditions of the first order now makes sense. Assume further:

(H6) The function $r$ (resp. $D$) is differentiable (resp. differentiable w.r.t. $c$).

(H7) The benefit function $f$ belongs to $\mathcal{C}^1$ and the codomain of the inverse function $(\partial_{\rho}f)^{-1}$ is a subset of $[0, \bar{p}]$.

The first order necessary conditions are derived using the Lagrangian associated with the problem (7) (see e.g., [3, 5, 9]). Canceling the variations of the Lagrangian associated to the problem with respect to the control $p$, to the state variable $c$ and to the final state $c(T, x)$, respectively provides the additional initial condition in (12), the optimality condition (10), the adjoint equation (11) and the boundary conditions in (12) below. With the state equations, we get the characterization of the optimal control $p^*$ by the following system in $\Omega_T$

$$ \mathrm{div}(v) = \theta, \quad v = -\kappa \nabla \phi, \qquad (8) $$

$$ R\psi \partial_t c + \mathrm{div}(v_c) = -r(c) + \chi_{\mathcal{S}} p^*, \quad v_c = -S_m \psi \nabla c + cv, \qquad (9) $$

$$ \partial_{\rho} f(x, p^*) = \mathcal{I}_T \mu, \qquad (10) $$

$$ \begin{aligned} R\psi \partial_t \mu + \mathrm{div}(v_\mu) &= -(r'(\mathcal{I}_{TC}) + R\psi \rho + \mathcal{I}_T \theta) \mu + \partial_c D(x, \mathcal{I}_{TC}), \\ v_\mu &= -S_m \psi \nabla \mu - \mu \mathcal{I}_T v, \end{aligned} \qquad (11) $$

completed by the following initial and boundary conditions

$$ \left\{ \begin{array}{ll} v_c \cdot n = 0 & \text{on } \Gamma_0, \quad v_c \cdot n = f_1 \quad \text{on } \Gamma, \quad v \cdot n = 0 & \text{on } \Gamma_0 \setminus \Gamma_d, \\ v \cdot n = v_1 & \text{on } \Gamma, \quad \phi = 0 & \text{on } \Gamma_d, \quad v_\mu \cdot n = f_2 \quad \text{on } \Gamma, \\ c|_{t=0} = c_0, \quad R\psi \mu|_{t=0} = \zeta \partial_c D(\cdot, c|_{t=T}). \end{array} \right. \qquad (12) $$

Notice that we have used a time variable change in the adjoint unknown to avoid handling an anti-diffusive equation for the adjoint variable. The time inversion operator, denoted by $\mathcal{I}_T$, is defined by

$$ \mathcal{I}_T f(t) = f(T - t), \quad \text{a.e. } t \in (0, T) $$

for any time-dependent function $f$.

We end this section with a proper definition of a weak solution of (8)–(12).

**Definition 2.2 (weak solution).** A function $(\phi, c, \mu)$ is a weak solution of (8)–(12) if $\phi \in L^\infty(0, T; H_1^1(\Omega))$, $c \in L^2(0, T; H^1(\Omega))$ and $\mu \in L^2(0, T; H^1(\Omega))$ satisfy: for all $(\varphi, \Phi) \in H^1(\Omega) \times L^2(0, T; H^1(\Omega))$ with $\varphi|_{\Gamma_d} = 0$,

$$ \int_{\Omega} \kappa \nabla \phi \cdot \nabla \varphi dx = \int_{\Omega} \theta \varphi dx + \int_{\Gamma} v_1 \varphi d\sigma, \qquad (13) $$

$$ \begin{aligned} & R\psi \int_0^T \int_{\Omega} \partial_t c \Phi + \int_0^T \int_{\Omega} S_m \psi \nabla c \cdot \nabla \Phi dxdt - \int_0^T \int_{\Omega} cv \cdot \nabla \Phi dxdt \\ & = - \int_0^T \int_{\Omega} (r(c) + \chi_{\mathcal{S}} (\partial_{\rho} f)^{-1} (\mathcal{I}_T \mu)) \Phi + \int_0^T \int_{\Gamma} f_1 \Phi d\sigma dt, \end{aligned} \qquad (14) $$
$$ 
\begin{aligned}
& R\psi \int_0^T \int_{\Omega} \partial_t \mu \Phi + \int_0^T \int_{\Omega} S_m \psi \nabla \mu \cdot \nabla \Phi dxdt + \int_0^T \int_{\Omega} \mu \mathcal{J}_{TV} \cdot \nabla \Phi dxdt \\
& = - \int_0^T \int_{\Omega} \left( (\rho' (\mathcal{J}_{TC}) + R\psi \rho + \theta) \mu - \partial_c D(x, \mathcal{J}_{TC}) \right) \Phi + \int_0^T \int_{\Gamma} f_2 \Phi \, d\sigma \, dt. \quad (15)
\end{aligned}
$$

### 3 | DEFINITION OF THE NUMERICAL SCHEME AND SKETCH OF THE MAIN RESULTS

We propose a discretization of the Euler-Lagrange's version (8)-(12) of the problem. The difficulty lies in the coupling of the diffusive Equation (9) with the antidiffusive Equation (11). Of course, there are not enough starting conditions available at either end point to produce a unique approximated solution. Such problems are known as two-point boundary value problems and are notably more difficult to solve than initial value problems. Due to the structure of the Pontryagin's maximum principle, the numerical solution of forward-backward problems is a classical issue in optimal control and a consequent literature exists on the subject. There are three main classes of methods, namely direct and indirect (multiple) shooting, relaxation, iteration. From a computational performance point of view, it is generally accepted that the first two methods are the most efficient. Nevertheless, all of these methods require an intelligent guess for the missing initial conditions (even if very effective refinements exist, see for instance [12]). Moreover the convergence speed calculations for the three categories of methods are more or less based on the same type of hypothesis, roughly speaking the function describing the solution trajectories of the differential system must be Lipschitz in time. This leads to problems that can become very troublesome when handling with PDEs.

As we want to prove Theorem 1 through the convergence of our iterative scheme without too many restrictive assumptions, we choose here to use another iterative method. Basically, the idea is to construct a finite volume architecture embedded in an iterative fixed point approximation for handling the nonlinear coupling due to (10) within the problem. The finite volume scheme is based on a two-point flux approximation with upwind mobilities. The time dynamics is treated with a Euler scheme. The mathematical analysis of the discrete problem, in particular its convergences properties, will allow us to ensure the existence of a solution to the optimal control problem in the sense of Definition 2.2. Then, the performance of the scheme will be illustrated at the end of the article by numerical illustrations.

We first introduce some notations. Let $\mathcal{T}$ be a regular and admissible mesh of the domain $\Omega$, constituted of open convex polygons called control volumes with maximum size (diameter) denoted by $h$. For all $K \in \mathcal{T}$, we denote by $|K|$ the measure of $K$. Let $N(K)$ be the set of the neighbors of $K$ that is, the set of cells of $\mathcal{T}$ which have a common interface with $K$. More precisely, $N_{int}(K)$ is the set of the neighbors of $K$ located in the interior of $\mathcal{T}$ while $N_{ext}(K)$ is the set of edges of $K$ on the boundary $\partial\Omega$. Recall (see [11]) that the admissibility of $\mathcal{T}$ implies that $\bar{\Omega} = \cup_{K \in \mathcal{T}} \bar{K}$, $K \cap L = \emptyset$ if $K$ and $L$ are two distinct elements of $\mathcal{T}$, that there exists a finite sequence of points $(x_K)_{K \in \mathcal{T}}$ such that $x_K \in K$ and such that the straight line $\overline{x_K x_L}$ is orthogonal to the edge $\sigma_{K,L}$. For all $L \in N_{int}(K)$ denote by $d_{K,L}$ the distance between $x_K$ and $x_L$, by $\sigma_{K,L}$ the interface between $K$ and $L$, by $\eta_{K,L}$ the unit normal vector to $\sigma_{K,L}$ outward to $K$. And for all $\sigma \in N_{ext}(K)$, denote by $d_{K,\sigma}$ the distance from $x_K$ to $\sigma$. Further, we assume the classical shape regularity property for the mesh:

$$ \min_{K \in \mathcal{T}, L \in N(K)} d_{K,L}/\mathrm{diam}(K) \ge \vartheta, \quad \text{for some } \vartheta \in \mathbb{R}_+^*. $$
Let $N_T$ be the number of time steps. We set $\Delta t = T/N_T$, $t^n = n\Delta t$, $0 \le n \le N_T$. The discrete unknowns are denoted by (for $w = c$ or $w = \mu$)

$$w_K^n \approx w(t^n, x_K), \quad w_h^n(t, x) = w(t^n, x_K), \quad \forall t \in [t^n, t^{n+1}], \forall x \in K.$$

To lighten the notations a bit, we will write the scheme in the case where the flow regime is independent of time, that is, the case where the source term $\theta$ does not depend on time. For the initialization of the fixed point algorithm, we use an initial control

$$p_0 \in L^2((0, T) \times \mathcal{L}), \quad 0 \le p_0(t, x) \le \bar{p} \quad \text{a.e. in } (0, T) \times \mathcal{L}.$$

The discrete problem reads as follows: find $P = (p_K^{m,n})_{K \in \mathcal{F}}$, $\Phi = (\phi_K)_{K \in \mathcal{F}}$, $C = (c_K^{m,n})_{K \in \mathcal{F}}$, and $\Upsilon = (\mu_K^{m,n})_{K \in \mathcal{F}}$ such that for any $n \in [0, N_T] \cap \mathbb{N}$, for any $K \in \mathcal{F}$ and for any $m \in \mathbb{N}$

$$p_K^{0,n} = \frac{1}{|K|} \int_K p_0(x) x, \quad \mu_K^{m,0} = \frac{1}{|K|} \int_K c_0(x) x, \quad (16)$$

$$- \sum_{L \in N(K)} |\sigma_{K,L}| \frac{\phi_L - \phi_K}{d_{K,L}} = |K|\theta_K + |\partial K \cap \Gamma|v_1, \quad (17)$$

$$R\psi|K| \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} - \sum_{L \in N(K)} |\sigma_{K,L}| \left( \psi S_m \frac{c_L^{m,n+1} - c_K^{m,n+1}}{d_{K,L}} + \frac{\phi_L - \phi_K}{d_{K,L}} c_K^{m,n+1} \right) \\ = -|K| \left( r(c_K^{m,n}) - p_K^{m,n} \right) + |\partial K \cap \Gamma|f_1, \quad (18)$$

$$\mu_K^{m,0} = \frac{\zeta}{R\psi|K|} \int_K \partial_c D(x, c_K^{m,N_T}) x, \\ R\psi|K| \frac{\mu_K^{m,n+1} - \mu_K^{m,n}}{\Delta t} - \sum_{L \in N(K)} |\sigma_{K,L}| \left( \psi S_m \frac{\mu_L^{m,n+1} - \mu_K^{m,n+1}}{d_{K,L}} - \frac{\phi_L - \phi_K}{d_{K,L}} \mu_K^{m,n+1} \right) \\ = -|K| \left( r'(c_K^{m,N_T-n}) + R\psi\rho + \theta_K \right) \mu_K^{m,n} + |K|\partial_c D(c_K^{m,N_T-n}) + |\partial K \cap \Gamma|f_2, \quad (19)$$

$$p_K^{m+1,n} = (\partial_{\rho} f)^{-1} (\mu_K^{m,N_T-n}). \quad (20)$$

**Sketch of the results.** The main results of this article may be summarized as follows:

* There exists at least one solution of the discrete problem (16)–(20). This result, along with estimates and stability properties, are proved in Section 4.

* It can be proven that the sequence $(p_h^{m,n})$ built from (16)–(20) converges in $L^2(\Omega_T)$ and almost everywhere in $\Omega_T$ as $m \to \infty$, $N_T \to \infty$, $h \to 0$, to a solution of the optimal control problem (7) (see Sections 5–6). Appropriate assumptions ensuring the uniqueness of the optimal control are also provided.

With these results, we prove on the one hand the existence of a solution to the problem in the sense of Definition 2.2 and we give on the other hand a theoretical basis to numerical simulations.

**Remark 3.1** The latter results improve the existence result stated in [1]: the assumptions are relaxed and a numerical scheme devoted to the approximation of the optimal control problem is proposed and analyzed.
4 | ANALYSIS OF THE DISCRETE PROBLEM (16)–(20)

Once again, we introduce some notations. We denote by $L_h^2(\Omega) \subset L^2(\Omega)$ the space of functions which are piecewise constant on each control volume $K \in \mathcal{T}$. For any function $P_h \in L_h^2(\Omega)$ and for all $K \in \mathcal{T}$, we denote by $P_K$ the constant value of $p_h$ in $K$. Conversely, we associate $\phi_h$, $c_h^{m,n}$, $\mu_h^{m,n}$ and $p_h^{m,n}$ with the set of values defined by the scheme (16)–(20). In the same way, we define $(\bar{c}_h, \bar{\mu}_h)$ associated with

$$ \bar{w}(t, x) := \frac{w(t_{n+1}, x) - w(t_n, x)}{\Delta t} (t - t_n) + w(t_n, x), \quad \text{for all } t \in [t_n, t_{n+1}], \quad w = c, \mu. $$

We define some approximate gradients which are piecewise constant on dual diamond cells. To this aim, if $K \in \mathcal{T}$ and $L \in N(K)$ with common (open) vertex (or face in the three-dimensional case) $\sigma_{K,L}$, let $T_{K,L}$ be the convex polygon defined by

$$ T_{K,L} = \{tx_K + (1-t)y, \ t \in [0, 1], y \in \sigma_{K,L}\} \cup \{tx_L + (1-t)y, \ t \in [0, 1], y \in \sigma_{K,L}\}. $$

If an edge $\sigma$ of $K$ belongs to $N_{ext}(K)$, we define $T_{K,\sigma}^{ext}$ by

$$ T_{K,\sigma}^{ext} = \{tx_K + (1-t)y, \ t \in [0, 1], y \in \sigma\}. $$

Notice that the $d$-dimensional measure $|T_{K,L}|$ equals to $|\sigma_{K,L}|d_{K,L}/d$. The discrete gradient $\nabla_h p_h$ of a constant per control volume function $p_h$ is defined as the following constant per diamond $T_{K,L}$ $\mathbb{R}^d$-valued function

$$ \nabla_h p_h(x) = \begin{cases} d \frac{P_L - P_K}{d_{K,L}} \eta_{K,L} & \text{if } x \in T_{K,L}, \\ d \frac{P_\sigma - P_K}{d_{K,\sigma}} \eta_{K,\sigma} & \text{if } x \in T_{K,\sigma}^{ext}. \end{cases} $$

The semi-norm $\|p_h\|_{H_h^1}$ coincides with the $L^2(\Omega)$ norm of $\nabla_h p_h$, actually:

$$ \begin{aligned} \|\nabla_h p_h\|_{L^2(\Omega)}^2 &= \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \int_{T_{K,L}} |\nabla_h p_h|^2 dx = d^2 \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} |T_{K,L}| \frac{|P_L - P_K|^2}{|d_{K,L}|^2} \\ &= d \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \frac{|\sigma_{K,L}|}{d_{K,L}} |P_L - P_K|^2 = \|p_h\|_{H_h^1(\Omega)}^2. \end{aligned} $$

We first derive some a priori estimates.

**Proposition 4.1** (discrete a priori estimates). There exists a constant real number $C_\Omega > 0$, depending on $\Omega$, $T$, $c_0$, $\mu_0$, $\theta$, $p_0$ and $\rho$, such that any solution of the finite volume scheme (16)–(20) satisfies for all $m \in \mathbb{N}$, for all $n \in [0, N_T]$ and for all $h > 0$:

$$ \|\phi_h\|_{H_h^1}^2 \le C_\Omega, \qquad (21) $$

and, if $\Delta t < R\psi/4\theta^*$,

$$ \|c_h^{m,n}\|_{L^2(\Omega)}^2 \le C_\Omega, \quad \|c_h^{m,n}\|_{H_h^1(\Omega)}^2 \le C_\Omega, \quad \|\partial_t \bar{c}_h^m\|_{L^2(0,T;L^2(\Omega))}^2 \le C_\Omega, \qquad (22) $$

$$ \|\mu_h^{m,n}\|_{L^2(\Omega)}^2 \le C_\Omega, \quad \|\mu_h^{m,n}\|_{H_h^1(\Omega)}^2 \le C_\Omega, \quad \|\partial_t \bar{\mu}_h^m\|_{L^2(0,T;L^2(\Omega))}^2 \le C_\Omega. \qquad (23) $$

Furthermore, the approximate control is such that

$$ 0 \le p_h^{m,n} \le \bar{p} \quad \text{for all } m \in \mathbb{N}, n \in [0, N_T], h > 0. \qquad (24) $$
**Remark 4.2** A careful look at the proof of (22)–(23) below shows that it is sufficient to ensure that $\Delta t < R\psi/4 \|\min\{0, \theta\}\|_{\infty}$, meaning that the constraint is necessary only when water withdrawals from the aquifer are greater than natural recharge.

*Proof.* The bound (24) is of course a direct consequence of the assumption (H7). To prove the estimate (21), we multiply the hydraulic charge equation (17) by $\phi_K$. Summing the resulting equation over $K \in \mathcal{T}$ and setting $\tau_{K,L} = \frac{|\sigma_{K,L}|}{d_{K,L}}$ for the description of the transmissivity through the interface $\sigma_{K,L}$, we get

$$ - \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \tau_{K,L} (\phi_L - \phi_K) \phi_K = \sum_{K \in \mathcal{T}} |K| \theta_K \phi_K + \sum_{K \in \mathcal{T}} |\partial K \cap \Gamma| v_1 \phi_K. $$

Then, by using a discrete integration by parts, the definition of $\theta_K$ and the Cauchy-Schwartz inequality and Young inequalities, we obtain

$$ \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \tau_{K,L} (\phi_L - \phi_K)^2 \le \frac{1}{\epsilon} \|\theta\|_{L^2(\Omega)}^2 + \epsilon \|\phi_h\|_{L^2(\Omega)}^2 + \frac{1}{\epsilon} \|v_1\|_{L^2(\Gamma)}^2 + \epsilon \|\phi_h\|_{H_h^1}^2, $$

for any fixed $\epsilon > 0$ (we will choose its value later). Then, according to the discrete Poincaré inequality, there exist a positive constant $C_1$ such that

$$ (1/d - \epsilon C_1) \|\phi_h\|_{H_h^1}^2 \le (\|\theta\|_{L^2(\Omega)}^2 + \|v_1\|_{L^2(\Gamma)}^2) / \epsilon. $$

Thus, choosing $\epsilon$ such that $1/d - \epsilon C_1 > 0$, we obtain (21).

For the pollutant estimate, we inject (17) in the convective terms of (18) and we multiply the resulting equation by $c_K^{m,n+1}$. After summation over $K \in \mathcal{T}$, we obtain

$$ \begin{aligned} R\psi \sum_{K \in \mathcal{T}} |K| \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} c_K^{m,n+1} - \psi S_m \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} |\sigma_{K,L}| \frac{c_L^{m,n+1} - c_K^{m,n+1}}{d_{K,L}} c_K^{m,n+1} & \\ + \sum_{K \in \mathcal{T}} |K| \theta_K (c_K^{m,n+1})^2 + \sum_{K \in \mathcal{T}} |\partial K \cap \Gamma| v_1 (c_K^{m,n+1})^2 & \\ = - \sum_{K \in \mathcal{T}} |K| (r(c_K^{m,n}) - p_K^{m,n}) c_K^{m,n+1} \sum_{K \in \mathcal{T}} |\partial K \cap \Gamma| f_1 c_K^{m,n+1}. & \end{aligned} $$

After a discrete integration by parts, the latter equation reads

$$ \begin{aligned} R\psi \sum_{K \in \mathcal{T}} |K| \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} c_K^{m,n+1} + \frac{\psi S_m}{2} \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \tau_{K,L} (c_L^{m,n+1} - c_K^{m,n+1})^2 & \\ + \sum_{K \in \mathcal{T}} |K| \theta_K (c_K^{m,n+1})^2 + \sum_{K \in \mathcal{T}} |\partial K \cap \Gamma| v_1 (c_K^{m,n+1})^2 & \\ = - \sum_{K \in \mathcal{T}} |K| (r(c_K^{m,n}) - p_K^{m,n}) c_K^{m,n+1} + \int_{\Gamma} f_1 c_h^{m,n+1} dx. & \end{aligned} $$

Using the Cauchy-Schwartz and Young inequalities together with the assumptions on $r$, $\bar{p}$, $\theta$ and $f_1$, we infer from the latter result that, for any $\epsilon, \gamma > 0$, there exist two positive real numbers $C_3$ and $C_4$ such that

$$ \begin{aligned} & \frac{R\psi}{2\Delta t} (\|c_h^{m,n+1}\|_{L^2(\Omega)}^2 - \|c_h^{m,n}\|_{L^2(\Omega)}^2) + \frac{\psi S_m}{2d} \|c_h^{m,n+1}\|_{H_h^1}^2 + v_1 \|c_h^{m,n+1}\|_{L^2}^2 \\ & \le \frac{C_3}{4\epsilon} + \epsilon \|c_h^{m,n+1}\|_{L^2(\Omega)}^2 + \frac{1}{4\gamma} \|f_1\|_{L^2(\partial\Omega)}^2 + (\gamma C_4 + \theta^*) \|c_h^{m,n+1}\|_{L^2(\Omega)}^2 + \gamma C_4 \|c_h^{m,n+1}\|_{H_h^1}^2, \end{aligned} $$
Notice that the last two terms of the right-hand side of the latter relation come from the continuity of the discrete trace function. The result also reads

$$ \frac{R\psi}{2\Delta t} \left( \|c_h^{m,n+1}\|_{L^2(\Omega)}^2 - \|c_h^{m,n}\|_{L^2(\Omega)}^2 \right) + \left( \frac{\psi S_m}{2d} - \gamma C_4 \right) \|c_h^{m,n+1}\|_{H_h^1}^2 + v_1 \|c_h^{m,n+1}\|_{L^2}^2 \\ \le \frac{C_3}{4\epsilon} + \frac{1}{4\gamma} \|f_1\|_{L^2(\partial\Omega)}^2 + (\epsilon + \gamma C_4 + \theta^*) \|c_h^{m,n+1}\|_{L^2(\Omega)}^2. \quad (25) $$

Now, choosing $\gamma$ small enough for ensuring $\psi S_m/2d - \gamma C_4 > 0$, setting $K_1 = (\frac{C_3}{\epsilon} + \frac{1}{\gamma}\|f_1\|_{L^2(\partial\Omega)}^2)/2R\psi$, $K_2 = 2(\epsilon + \gamma C_4 + \theta^*)/R\psi$, we claim that, according to the discrete Gronwall lemma,

$$ \|c_h^{m,n+1}\|_{L^2(\Omega)}^2 \le (c_h^{m,0} + K_1/K_2) \exp(K_2 T / (1 - K_2 \Delta t)) \quad (26) $$

as soon as $0 < \Delta t \le \Delta t_0$, $\Delta t_0 < 1/K_2$. Some computations show that the latter condition may be ensured if $\Delta t < R\psi/4\theta^*$. Next, with (26) at hands, (25) implies the existence of a positive real number $C_5$ such that

$$ \|c_h^{m,n+1}\|_{H_h^1(\Omega)}^2 \le C_5 (c_h^{m,0} + K_1/K_2) \exp(K_2 T / (1 - K_2 \Delta t)). \quad (27) $$

It remains to justify the estimate of the time derivative announced in (22). To this aim, multiply (18) by $(c_K^{m,n+1} - c_K^{m,n})/\Delta t$ and sum over $K$ in order to obtain

$$ \begin{aligned} & R\psi \sum_{K \in \mathcal{T}} |K| \left( \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} \right)^2 + \sum_{K \in \mathcal{T}} \theta_K c_K^{m,n+1} \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} \\ & + \sum_{K \in \mathcal{T}} |\partial K \cap \Gamma| v_1 c_K^{m,n+1} \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} \\ & - \psi S_m \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} |\sigma_{K,L}| \frac{c_L^{m,n+1} - c_K^{m,n+1}}{d_{K,L}} \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} \\ & = - \sum_{K \in \mathcal{T}} |K| \left( r(c_K^{m,n}) - p_K^{m,n} \right) \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} + \sum_{K \in \mathcal{T}} |\partial K \cap \Gamma| f_1 \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t}. \end{aligned} $$

A discrete integration by parts yields

$$ \begin{aligned} & R\psi \sum_{K \in \mathcal{T}} |K| \left( \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} \right)^2 + \sum_{K \in \mathcal{T}} \theta_K c_K^{m,n+1} \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} \\ & + \sum_{K \in \mathcal{T}} |\partial K \cap \Gamma| v_1 c_K^{m,n+1} \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} \\ & + \frac{1}{2} \psi S_m \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} |\sigma_{K,L}| \left( \frac{c_L^{m,n+1} - c_K^{m,n+1}}{d_{K,L}} \right) \left( \frac{c_L^{m,n+1} - c_L^{m,n}}{\Delta t} - \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} \right) \\ & = - \sum_{K \in \mathcal{T}} |K| \left( r(c_K^{m,n}) - p_K^{m,n} \right) \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} + \sum_{K \in \mathcal{T}} |\partial K \cap \partial \Omega| f_1 \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t}. \end{aligned} $$

We rewrite $f_1$ in the form $f_1 = \sum_{i \in \mathcal{T}} f_{1,i} \chi_{\Gamma_i}$ where $\Gamma_1 \cup \Gamma_3 = \cup_{i \in \mathcal{T}} \Gamma_i$ and $f_{1,i}$ are given real numbers. Using the Cauchy-Schwartz and Young inequalities, we get
$$ (R\psi - 3\gamma)\|\partial_t \bar{c}_h^{m,n+1}\|_{L^2(\Omega)}^2 \\ + (\psi S_m/4\Delta t) \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \tau_{K,L} \left( (c_L^{m,n+1} - c_K^{m,n+1})^2 - (c_L^{m,n} - c_K^{m,n})^2 \right) \\ \le \left( C \left( \theta^*, v_1, \max_n \|c_h^{m,n}\|_{L^2(\Omega)}, \mathrm{Lip}(r), \bar{p} \right) + C \sum_{K \in \mathcal{T}} |\partial K \cap \partial \Omega| \sum_{i \in \mathcal{I}} f_{1,i}^2 \right) / \gamma $$

for any $\gamma > 0$. Next, choosing $\gamma$ small enough, multiplying the latter relation by $\Delta t$ and summing over $n$ yields

$$ \|\partial_t \bar{c}_h^m\|_{L^2(0,T;L^2(\Omega))}^2 + \psi S_m \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \tau_{K,L} \left( (c_L^{m,N_T} - c_K^{m,N_T})^2 - (c_L^{m,0} - c_K^{m,0})^2 \right) \\ \le T C \left( \theta^*, v_1, \max_n \|c_h^{m,n}\|_{L^2(\Omega)}, \mathrm{Lip}(r), \bar{p} \right) + C \|f_1\|_{L^2(0,T;L^2(\partial\Omega))}^2. \quad \blacksquare $$

Estimate (22) is proved. The proof of (23) follows the same lines.

We now aim at stating an existence result for the discrete problem (16)–(19). First we recall the following Lemma.

**Lemma 4.3** (see [11]). Let $\mathcal{A}$ a finite-dimensional Hilbert space. Denote by $[\cdot, \cdot]$ its inner product and by $\|\cdot\|$ the associated norm. Let $\mathcal{F}$ be a continuous map from $\mathcal{A}$ in itself and a real number $r > 0$ such that $[\mathcal{F}(\xi), \xi] > 0$ if $\|\xi\| = r$. Then there exists $\xi \in \mathcal{A}$ with $\|\xi\| \le r$ such that $\mathcal{F}(\xi) = 0$.

We have sufficient tools for proving the following existence result for the discrete problem (16)–(19).

**Theorem 4.4** (existence of a discrete solution). Assume (H1)–(H7). There exists a solution $(\phi_K, c_K^{m,n}, \mu_K^{m,n})_{m,n,K}$ of the discrete problem (16)–(20).

*Proof.* The scheme (16)–(19) is explicit with regard to $m$ and $n$. It is thus sufficient to consider the question of the solution of the finite volume approximation. Since the equations are not strongly coupled we show the existence result of one unknown, namely the discrete concentration. The other proofs are similar. Assume that $(c_K^{m,n})_{K \in \mathcal{T}}$ are given real numbers satisfying (22). We set $\mathcal{M} := \mathrm{Card}(\mathcal{T}_h)$, $C_{\mathcal{M}} := \{c_K^{m,n+1}\}_{K \in \mathcal{T}_h} \in \mathbb{R}^{\mathcal{M}}$ and we define a map $\mathcal{F}_h : \mathbb{R}^{\mathcal{M}} \to \mathbb{R}^{\mathcal{M}}$ by

$$ \mathcal{F}_h(C_{\mathcal{M}}) = \{\mathcal{F}_K\}_{K \in \mathcal{T}_h}, \\ \mathcal{F}_K = R\psi|K| \frac{c_K^{m,n+1} - c_K^{m,n}}{\Delta t} - \psi S_m \sum_{L \in N(K)} |\sigma_{K,L}| \frac{c_L^{m,n+1} - c_K^{m,n+1}}{d_{K,L}} \\ - c_K^{m,n+1} \sum_{L \in N(K)} |\sigma_{K,L}| \frac{\phi_L - \phi_K}{d_{K,L}} + |K|(r(c_K^{m,n}) + p_K^{m,n}) - |\partial K \cap \partial \Omega|f_1. $$

Notice that $\mathcal{F}_h$ is a well-defined and continuous map. Bearing in mind Lemma 4.3, we aim at proving that $[\mathcal{F}_h(C_{\mathcal{M}}), C_{\mathcal{M}}] > 0$ for $\|C_{\mathcal{M}}\| = r$ if $r$ is large enough. Using Proposition 4.1 and repeating the computations that gave (25) we write
$$[\mathcal{F}_h(C_{\mathcal{M}}), C_{\mathcal{M}}] \geq \left(\frac{R\psi}{2\Delta t} - \gamma C_4 - \theta^* - \varepsilon\right) \|c_h^{m,n+1}\|_{L^2(\Omega)}^2 - \frac{R\psi}{2\Delta t} \|c_h^{m,n}\|_{L^2(\Omega)}^2 \\ + \left(\frac{\psi S_m}{2d} - \gamma C_4\right) \|c_h^{m,n+1}\|_{H_h^1}^2 - \frac{C_3}{4\varepsilon} - \frac{1}{4\gamma} \|f_1\|_{L^2(\partial\Omega)}^2.$$

If $\gamma$, $\varepsilon$ and $\Delta t$ are small enough we thus have

$$[\mathcal{F}_h(C_{\mathcal{M}}), C_{\mathcal{M}}] \geq C(\gamma, \varepsilon, \Delta t) \|c_h^{m,n+1}\|_{L^2(\Omega)}^2 + C_1(\gamma, \Delta t) \|c_h^{m,n+1}\|_{H_h^1}^2 \\ - R\psi \|c_h^{m,n}\|_{L^2(\Omega)}^2 / 2\Delta t - C_3 / 4\varepsilon - \|f_1\|_{L^2(\partial\Omega)}^2 / 4\gamma$$

where $C(\gamma, \varepsilon, \Delta t)$, $C_1(\gamma, \Delta t)$ and $C_3$ are given positive real numbers. Due to the regularity of the mesh, $\|c_h^{m,n+1}\|_{L^2(\Omega)}^2 \geq C\theta\|C_{\mathcal{M}}\|^2$, $C > 0$. Also, $\|c_h^{m,n}\|_{L^2(\Omega)}$ is bounded as in (22). Hence, if $r$ is large enough, $\|C_{\mathcal{M}}\| = r$ ensures $[\mathcal{F}_h(C_{\mathcal{M}}), C_{\mathcal{M}}] > 0$. The existence result in Theorem 4.4 follows from Lemma 4.3. ■

We now aim at analyzing the properties of the scheme. We begin by considering the stability of the scheme with regard to the initialization of the control.

**Proposition 4.5 (stability).** Let $\varepsilon \in \mathbb{R}$. For $n \in [0, N_T] \cap \mathbb{N}$, set $p_{\varepsilon,h}^{0,n} = p_h^{0,n} + \varepsilon$, $c_{\varepsilon,h}^{0,n} = c_h^{0,n}$ and $\mu_{\varepsilon,h}^{0,n} = \mu_h^{0,n}$. There exists a solution $(\phi_K, c_{\varepsilon,K}^{m,n+1}, \mu_{\varepsilon,K}^{m,n+1})_{K \in \mathcal{F}}$ of the following problem for any $m \in \mathbb{N}$ and any $n \in [0, N_T] \cap \mathbb{N}$

$$- \sum_{L \in N(K)} |\sigma_{K,L}| \frac{\phi_L - \phi_K}{d_{K,L}} = |K|\theta_K + \sum_{K \in \mathcal{F}} |\partial K \cap \Gamma| v_1 \phi_K, \quad (28)$$

$$R\psi|K| \frac{c_{\varepsilon,K}^{m,n+1} - c_{\varepsilon,K}^{m,n}}{\Delta t} - \sum_{L \in N(K)} |\sigma_{K,L}| \left( \psi S_m \frac{c_{\varepsilon,L}^{m,n+1} - c_{\varepsilon,K}^{m,n+1}}{d_{K,L}} + \frac{\phi_L - \phi_K}{d_{K,L}} c_{\varepsilon,K}^{m,n+1} \right) \\ = -|K| (r(c_{\varepsilon,K}^{m,n}) - p_{\varepsilon,K}^{m,n}) + |\partial K \cap \partial \Omega| f_1, \quad (29)$$

$$R\psi|K| \frac{\mu_{\varepsilon,K}^{m,n+1} - \mu_{\varepsilon,K}^{m,n}}{\Delta t} - \sum_{L \in N(K)} |\sigma_{K,L}| \left( \psi S_m \frac{\mu_{\varepsilon,L}^{m,n+1} - \mu_{\varepsilon,K}^{m,n+1}}{d_{K,L}} - \frac{\phi_L - \phi_K}{d_{K,L}} \mu_{\varepsilon,K}^{m,n+1} \right) \\ = -|K| \left( r' \left(c_{\varepsilon,K}^{m,N_T-n}\right) + R\psi\rho + \theta_K \right) \mu_{\varepsilon,K}^{m,n} + |K| \partial_c D \left(c_{\varepsilon,K}^{m,N_T-n}\right) + |\partial K \cap \partial \Omega| f_2, \quad (30)$$

$$p_{\varepsilon,K}^{m+1,n} = (\partial_{\mu} f)^{-1} (\mu_{\varepsilon,K}^{m,N_T-n}), \quad (31)$$

and it satisfies

$$\begin{align*} \|c_{\varepsilon,K}^{m,n+1} - c_K^{m,n+1}\|_{L^2(\Omega)} &\leq C(h, \Delta t)\varepsilon, & \|c_{\varepsilon,K}^{m,n+1} - c_K^{m,n+1}\|_{H_h^1(\Omega)} &\leq C_1(h, \Delta t)\varepsilon, \\ \|\mu_{\varepsilon,K}^{m,n+1} - \mu_K^{m,n+1}\|_{L^2(\Omega)} &\leq C_2(h, \Delta t)\varepsilon, & \|\mu_{\varepsilon,K}^{m,n+1} - \mu_K^{m,n+1}\|_{H_h^1(\Omega)} &\leq C_3(h, \Delta t)\varepsilon, \\ \|p_{\varepsilon,K}^{m,n+1} - p_K^{m,n+1}\|_{L^2(\Omega)} &\leq C_4(h, \Delta t)\varepsilon, \end{align*}$$

where $C(h, \Delta t)$ and $C_i(h, \Delta t)$ are positive real numbers and $(\phi_K, c_K^{m,n+1}, \mu_K^{m,n+1})$ is a solution of the discrete problem (16)–(19). Moreover, if the sequences $(c_{\varepsilon,K}^{m,n})$ and $(c_K^{m,n})$ converge a.e. in $\Omega_T$ as $(h, \Delta t) \to (0, 0)$, then

$$\lim_{(h, \Delta t) \to (0,0)} C(h, \Delta t) = \lim_{(h, \Delta t) \to (0,0)} C_i(h, \Delta t) = 0, \quad 1 \leq i \leq 4.$$
**Proof.** The existence of a solution to (28)–(30) is ensured by an obvious analogous of Proposition 4.1. Set $P_{\epsilon,K}^{m,n} = p_{\epsilon,K}^{m,n} - p_{K}^{m,n}$, $\Psi_{\epsilon,K}^{m,n} = c_{\epsilon,K}^{m,n} - c_{K}^{m,n}$ and $\Upsilon_{\epsilon,K}^{m,n} = \mu_{\epsilon,K}^{m,n} - \mu_{K}^{m,n}$. These quantities solve the following problem:

$$ R\psi|K| \frac{\Psi_{\epsilon,K}^{m,n+1} - \Psi_{\epsilon,K}^{m,n}}{\Delta t} - \psi S_m \sum_{L \in N(K)} |\sigma_{K,L}| \frac{\Psi_{\epsilon,L}^{m,n+1} - \Psi_{\epsilon,K}^{m,n+1}}{d_{K,L}} + |K|\theta_K \Psi_{\epsilon,K}^{m,n+1} \\ + \sum_{K \in \mathcal{T}} |\partial K \cap \partial \Omega| v_1 \Psi_{\epsilon,K}^{m,n+1} = -|K| \left( (r(c_{\epsilon,K}^{m,n}) - r(c_{K}^{m,n})) - P_{\epsilon,K}^{m,n} \right), \quad (32) $$

$$ R\psi|K| \frac{\Upsilon_{\epsilon,K}^{m,n+1} - \Upsilon_{\epsilon,K}^{m,n}}{\Delta t} - \psi S_m \sum_{L \in N(K)} |\sigma_{K,L}| \frac{\Upsilon_{\epsilon,L}^{m,n+1} - \Upsilon_{\epsilon,K}^{m,n+1}}{d_{K,L}} - |K|\theta_K \Upsilon_{\epsilon,K}^{m,n+1} \\ - \sum_{K \in \mathcal{T}} |\partial K \cap \partial \Omega| v_1 \Upsilon_{\epsilon,K}^{m,n+1} = -|K| \left( r' \left( c_{\epsilon,K}^{m,N_T-n} \right) - r' \left( c_{K}^{m,N_T-n} \right) \right) \Upsilon_{\epsilon,K}^{m,n} \\ - |K| (R\psi\rho + \theta_K) \Upsilon_{\epsilon,K}^{m,n} + |K| \left( \partial_c D \left( c_{\epsilon,K}^{m,N_T-n} \right) - \partial_c D \left( c_{K}^{m,N_T-n} \right) \right). \quad (33) $$

We have $P_{\epsilon,K}^{0,n} = \epsilon$. Thus multiplying (32) and (33) with $m=0$ by respectively $\Psi_{\epsilon,K}^{0,n+1}$ and $\Upsilon_{\epsilon,K}^{0,n+1}$ and summing for $K \in \mathcal{T}$ we obtain

$$ R\psi \sum_{K \in \mathcal{T}} |K| \frac{\Psi_{\epsilon,K}^{0,n+1} - \Psi_{\epsilon,K}^{0,n}}{\Delta t} \Psi_{\epsilon,K}^{0,n+1} - \psi S_m \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} |\sigma_{K,L}| \frac{\Psi_{\epsilon,L}^{0,n+1} - \Psi_{\epsilon,K}^{0,n+1}}{d_{K,L}} \Psi_{\epsilon,K}^{0,n+1} \\ + \sum_{K \in \mathcal{T}} |K| \theta_K (\Psi_{\epsilon,K}^{0,n+1})^2 + \sum_{K \in \mathcal{T}} |\partial K \cap \partial \Omega| v_1 (\Psi_{\epsilon,K}^{0,n+1})^2 \\ = - \sum_{K \in \mathcal{T}} |K| \left( (r(c_{\epsilon,K}^{0,n}) - r(c_{K}^{0,n})) - \epsilon \right) \Psi_{\epsilon,K}^{0,n+1}, \quad (34) $$

$$ R\psi \sum_{K \in \mathcal{T}} |K| \frac{\Upsilon_{\epsilon,K}^{0,n+1} - \Upsilon_{\epsilon,K}^{0,n}}{\Delta t} \Upsilon_{\epsilon,K}^{0,n+1} - \psi S_m \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} |\sigma_{K,L}| \frac{\Upsilon_{\epsilon,L}^{0,n+1} - \Upsilon_{\epsilon,K}^{0,n+1}}{d_{K,L}} \Upsilon_{\epsilon,K}^{0,n+1} \\ - \sum_{K \in \mathcal{T}} |K| \theta_K (\Upsilon_{\epsilon,K}^{0,n+1})^2 - \sum_{K \in \mathcal{T}} |\partial K \cap \partial \Omega| v_1 (\Upsilon_{\epsilon,K}^{0,n+1})^2 \\ = - \sum_{K \in \mathcal{T}} |K| \left( (r' (c_{\epsilon,K}^{0,N_T-n}) - r' (c_{K}^{0,N_T-n}) - (R\psi\rho + \theta_K + r' (c_{\epsilon,K}^{0,N_T-n}))) \Upsilon_{\epsilon,K}^{0,n} \Upsilon_{\epsilon,K}^{0,n+1} \right. \\ \left. + (\partial_c D (c_{\epsilon,K}^{0,N_T-n}) - \partial_c D (c_{K}^{0,N_T-n})) \Upsilon_{\epsilon,K}^{0,n+1} \right). \quad (35) $$

Hence, after a discrete integration by parts in (34), using the Lipschitz character of $r$ and the Gronwall Lemma we get first for any $n \in [0, N_T - 1]$

$$ \| \Psi_{\epsilon,h}^{0,n+1} \|_{L^2(\Omega)} \le C(h, \Delta t)\epsilon, \quad \| \Psi_{\epsilon,h}^{0,n+1} \|_{H_h'(\Omega)} \le C_1(h, \Delta t)\epsilon. \quad (36) $$

Next, using the continuity of $r'$ and $\partial_c D$ and (23), we infer from (35)–(36) that

$$ \| \Upsilon_{\epsilon,K}^{0,n+1} \|_{L^2(\Omega)} \le C_2(h, \Delta t)\epsilon, \quad \| \Upsilon_{\epsilon,K}^{0,n+1} \|_{H_h'(\Omega)} \le C_3(h, \Delta t)\epsilon. \quad (37) $$

Then, using the continuity of the map $(\partial_p f)^{-1}$ defining the discrete controls $p$ with regard to the discrete adjoint unknowns $\mu$, we infer from the latter estimates that there exists a positive real number $C_4(h, \Delta t)$ such that

$$ \| P_{\epsilon,h}^{1,n} \|_{L^2(\Omega)} \le C_4(h, \Delta t)\epsilon. \quad (38) $$
The estimates in Proposition 4.5 follows by induction with similar computations. Checking that $\lim_{(h,\Delta t)\to(0,0)} C(h, \Delta t) = \lim_{(h,\Delta t)\to(0,0)} C_i(h, \Delta t) = 0$ is obvious if the sequences $(c_{\epsilon,K}^{m,n})$ and $(c_K^{m,n})$ converge a.e. in $\Omega_T$.

The next result is a first stone towards the construction of compactness properties.

**Proposition 4.6** (space and time translations). There exist positive real numbers $C, C_1, C_2$, and $C_3$ such that, for any $m \in \mathbb{N}$, for any $\zeta \in \mathbb{R}^d$ and any $\eta \in ]0, T[$, we have

$$ \iint_{\Omega' \times [0, T]} |\bar{c}_h^m(t, x + \zeta) - \bar{c}_h^m(t, x)|^2 dxdt \le C|\zeta|(|\zeta| + 2h), \quad (39) $$

$$ \iint_{\Omega' \times [0, T]} |\bar{\mu}_h^m(t, x + \zeta) - \bar{\mu}_h^m(t, x)|^2 dxdt \le C_1 |\zeta|(|\zeta| + 2h), \quad (40) $$

$$ \iint_{\Omega \times [0, T - \eta]} |\bar{c}_h^m(t + \eta, x) - \bar{c}_h^m(t, x)|^2 dxdt \le C_2(\eta + \Delta t), \quad (41) $$

$$ \iint_{\Omega \times [0, T - \eta]} |\bar{\mu}_h^m(t + \eta, x) - \bar{\mu}_h^m(t, x)|^2 dxdt \le C_3(\eta + \Delta t). \quad (42) $$

*Proof.* With the estimates in Proposition 4.1 at hands, follow the lines of the proofs of Lemmas 4.3 and 4.7 in the classical textbook [11].

## 5 | CONVERGENCE ANALYSIS

In this section, we claim and prove various convergence results. They are summarized in the following theorem.

**Theorem 5.1** (convergence results). Assume (H1)–(H7). Assume further the continuity of the derivatives $r'$ and $\partial_c D$, with $r'$ bounded and $\partial_c D$ with a at most affine growth with respect to $c$. There exists a subsequence of solutions $(\phi_h, c_h^{m,n}, \mu_h^{m,n}, p_h^{m,n})$ of (16)–(20) that converges in $L^2(\Omega_T)$ and almost everywhere in $\Omega_T$ as $m \to \infty, N_T \to \infty, h \to 0$, to a solution $(\phi, c, \mu, p)$ of the Euler–Lagrange’s problem (8)–(12).

*Proof.* For the proof of Theorem 5.1 we could mind a priori that the outcome depends on the way we pass to the limit $m \to \infty$, $(h, \Delta t) \to (0, 0)$ in (16)–(20). We thus detail all the possible combinations.

**Convergence 1: letting** $(h, \Delta t) \to (0, 0)$ **then** $m \to \infty$.

Let $c_h^m$ the function defined in $\Omega_T$ by $c_h^m = c_K^{m,n}$, for any $K \in \mathcal{T}$ and $n \in [1, N_T] \cap \mathbb{N}$. This function should rather be denoted by $c_{h,\Delta t}^m$, but, for the sake of the simplicity and since the limit behaviors of the function for vanishing $h$ and $\Delta t$ are not studied separately, we choose this slight notations abuse. The proofs for the $\phi, c$ and $\mu$ approximations follow the same lines. We thus focus on the limit behavior of $c_h^m$. Let $m > 0$. Due to the first estimate in (22), there exists some subsequence, not relabeled for convenience, and some $c_h^m \in L^2(\Omega_T)$ such that

$$ c_h^m \to c^m \text{ weakly in } L^2(\Omega_T) \text{ as } h, \Delta t \to 0. $$
Moreover, according to Proposition 4.1, there exists $\lambda^m \in (L^2(\Omega_T))^d$ such that, up to a subsequence, $\nabla_h c_h^m$ weakly converges to $\lambda^m$ in $(L^2(\Omega_T))^d$ as $h, \Delta t \to 0$. It remains to identify $\nabla c^m$ and $\lambda^m$ in the sense of distributions. To this aim, it is enough to show that, as $h \to 0$,

$$E_h := \iint_{\Omega_T} \nabla_h c_h^m \cdot \varphi \, dxdt + \iint_{\Omega_T} c_h^m \operatorname{div} \varphi \, dxdt \to 0, \quad \forall \varphi \in (D(\Omega_T))^d.$$

Let $h$ be small enough so that $\varphi$ vanishes in $T_{K,\sigma}^{ext}$ for all $K \in \mathcal{T}$. Since $\eta_{K,L} = -\eta_{L,K}$ we have on the one hand for all $t \in ]t^n, t^{n+1}[$

$$\begin{aligned} \int_{\Omega} c_h^m \operatorname{div} \varphi(t, x) dx &= \sum_{K \in \mathcal{T}} \int_K c_h^m \operatorname{div} \varphi(t, x) dx = \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} c_K^m \int_{\sigma_{K,L}} \varphi(t, s) \cdot \eta_{K,L} ds \\ &= \frac{1}{2} \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} (c_K^m - c_L^m) \int_{\sigma_{K,L}} \varphi(t, s) \cdot \eta_{K,L} ds. \end{aligned}$$

On the other hand, we infer from the definition of the discrete gradient that

$$\begin{aligned} \int_{\Omega} \nabla_h c_h^m \cdot \varphi(t, x) dx &= \frac{1}{2} \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \int_{T_{K,L}} \nabla_h c_h^m \cdot \varphi(t, x) dx \\ &= \frac{1}{2} \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \frac{d}{dK_L} (c_L^m - c_K^m) \int_{T_{K,L}} \varphi(t, s) \cdot \eta_{K,L} dx. \end{aligned}$$

Inserting the two latter relations in the definition of $E_h$, we get

$$E_h = \sum_{n=0}^{N_T-1} \frac{\Delta t}{2} \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \sigma_{K,L} (c_K^m - c_L^m) \left( \frac{1}{\sigma_{K,L}} \int_{\sigma_{K,L}} \varphi(t,s) \cdot \eta_{K,L} ds - \frac{1}{T_{K,L}} \int_{T_{K,L}} \varphi(t,s) \cdot \eta_{K,L} dx \right).$$

Due to the smoothness of $\varphi$, we have

$$\left| \frac{1}{\sigma_{K,L}} \int_{\sigma_{K,L}} \varphi(t,s) \cdot \eta_{K,L} ds - \frac{1}{T_{K,L}} \int_{T_{K,L}} \varphi(t,s) \cdot \eta_{K,L} dx \right| \le Ch.$$

Thus, with the Cauchy–Scharwz inequality and the first estimate in (22), we obtain

$$\begin{aligned} |E_h| &\le Ch \sum_{n=0}^{N_T-1} \Delta t \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \sigma_{K,L} (c_K^m - c_L^m) \\ &\le Ch \sum_{n=0}^{N_T-1} \Delta t \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \frac{\sigma_{K,L}}{d_{K,L}} (c_K^m - c_L^m)^2 + Ch \sum_{n=0}^{N_T-1} \Delta t \sum_{K \in \mathcal{T}} \sum_{L \in N(K)} \sigma_{K,L} d_{K,L} \\ &\le Ch|\Omega|T \end{aligned}$$

which vanishes as $h \to 0$. We thus claim that

$$\nabla c_h^m \to \nabla c^m \text{ weakly in } (L^2(\Omega_T))^d \text{ as } h, \Delta t \to 0. \quad (43)$$

Thanks to the last estimate in (22), similar computations lead to

$$\partial_t \bar{c}_h^m \to \partial_t c^m \text{ weakly in } L^2(\Omega_T) \text{ as } h, \Delta t \to 0. \quad (44)$$
Let us now look for a sequential compactness result for handling the nonlinearities in the problem. Writing for any $\Omega' \in \Omega$ and $\eta > 0$

$$ \iint_{\Omega' \times (0, T-\eta)} |\bar{c}_h^m(t+\eta, x+\zeta) - \bar{c}_h^m(t,x)|^2 \\ \le \iint_{\Omega' \times [0, T]} |\bar{c}_h^m(t, x+\zeta) - \bar{c}_h^m(t,x)|^2 + \iint_{\Omega \times (0, T-\eta)} |\bar{c}_h^m(t+\eta, x) - \bar{c}_h^m(t,x)|^2 $$

allows the use of the space and time translation results (39) and (41) to get

$$ \| \bar{c}_h^m(t + \eta, x + \zeta) - \bar{c}_h^m(t, x) \|_{L^2(\mathbb{R}^{d+1})} \to 0 \quad \text{as} \quad (\eta, \zeta, \Delta t) \to 0. $$

The relative compactness of $(c_h^m)$ in $L^2(\Omega_T)$ then follows by the Riesz-Fréchet-Kolmogorov compactness criterion. There exists a subsequence $(c_h^m)$ such that

$$ c_h^m \to c^m \text{ in } L^2(\Omega_T) \text{ and a.e in } \Omega_T \text{ as } h, \Delta t \to 0. \quad (45) $$

Similar arguments show the existence of the limit functions $\phi \in L^2(0, T; H^1(\Omega))$ and $\mu^m \in L^2(0, T; H^1(\Omega)) \cap H^1(0, T; L^2(\Omega))$ such that, up to a subsequence, the following convergence results hold true as $h, \Delta t \to 0$:

$$ \nabla_h \phi_h \to \nabla \phi \text{ weakly in } (L^2(\Omega_T))^d, \quad (46) $$

$$ \mu_h^m \to \mu^m \text{ in } L^2(\Omega_T) \text{ and a.e in } \Omega_T, \quad (47) $$

$$ \nabla_h \mu_h^m \to \nabla \mu^m \text{ and } \partial_t \bar{\mu}_h^m \to \partial_t \mu^m \text{ weakly in } (L^2(\Omega_T))^d. \quad (48) $$

Now, let $\varphi \in C_c^\infty([0, T[ \times \bar{\Omega})$ and $\varphi_K^n = \varphi(t^n, x_K)$, $K \in \mathcal{T}_h$, $n \in [0, N_T] \cap \mathbb{N}$. Notice the strong convergence of $\varphi_h$ and $(\nabla \varphi)_h$ to respectively $\varphi$ and $\nabla \varphi$ as $h, \Delta t \to 0$. Multiplying the Equations (16)–(19) by $\varphi_K^{n+1}$ and summing the results over $K \in \mathcal{T}_h$ and $n=0, \dots, N_T-1$, we can rewrite the scheme (16)–(20) as follows (see e.g., Section 3.6 in [19])

$$ \int_{\Omega_T} \kappa \nabla \varphi_h \cdot (\nabla \varphi)_h dx = \int_{\Omega} \theta \varphi_h dx + \int_{\Gamma} v_1 \varphi_h d\sigma, \quad (49) $$

$$ \begin{aligned} & R\psi \int_{\Omega_T} \partial_t \bar{c}_h^m \varphi_h dxdt + \int_{\Omega_T} S_m \psi \nabla c_h^m \cdot (\nabla \varphi)_h dxdt - \int_{\Omega_T} c_h^m \nabla \varphi_h \cdot (\nabla \varphi)_h dxdt \\ & = - \int_{\Omega_T} \left( r(c_h^m) + \chi_{\mathcal{T}} p_h^m \right) (t - \Delta t, x) \varphi_h(t, x) dxdt + \int_0^T \int_{\Gamma} f_1 \varphi_h d\sigma dt, \end{aligned} \quad (50) $$

$$ \begin{aligned} & R\psi \int_{\Omega_T} \partial_t \bar{\mu}_h^m \varphi_h dxdt + \int_{\Omega_T} S_m \psi \nabla \mu_h^m \cdot (\nabla \varphi)_h dxdt + \int_{\Omega_T} \mu_h^m \nabla \varphi_h \cdot (\nabla \varphi)_h dxdt \\ & = - \int_{\Omega_T} \left( (r'(\mathcal{T}_T c_h^m) + R\psi \rho + \theta) \mu_h^m (t - \Delta t, x) - \frac{\partial D}{\partial c}(x, \mathcal{T}_T c_h^m) \right) \varphi_h dxdt \\ & \quad + \int_0^T \int_{\Gamma} f_2 \varphi_h d\sigma dt, \end{aligned} \quad (51) $$

where $p_h^{m,n}$ is defined by (20). The convergence results (43)–(48) are sufficient for passing to the limit $(h, \Delta t) \to (0, 0)$ in the latter three equations. The bounded nonlinear terms
may be handled with the Lebesgue's dominated convergence theorem and the term containing $\partial_c D$ with the Severini-Egorov's theorem (since $\partial_c D$ has at most an affine growth with regard to $c$, if $A$ is the arbitrarily small part of the domain where the uniform convergence of $c_h^m$ is not ensured, the quantity $|\int_A \partial_c D(x, \mathcal{F}_T c_h^m) \varphi_h|$ is also arbitrarily small). We obtain an analogous of the variational problem given in Definition 2.2 with $(c^m, \mu^m)$ instead of $(c, \mu)$.

It remains to let $m \to \infty$. Since all the estimates in Proposition 4.1 are uniform with respect to $m$, they also hold for $(c^m, \mu^m)$ (property of the lim inf). Hence

$$ \|c^m\|_{L^2(\Omega)}^2 \le C_\Omega, \quad \|\nabla c^m\|_{L^2(\Omega)}^2 \le C_\Omega, \quad \|\partial_t c^m\|_{L^2(0,T;L^2(\Omega))}^2 \le C_\Omega, \qquad (52) $$

$$ \|\mu^m\|_{L^2(\Omega)}^2 \le C_\Omega, \quad \|\nabla \mu^m\|_{L^2(\Omega)}^2 \le C_\Omega, \quad \|\partial_t \mu^m\|_{L^2(0,T;L^2(\Omega))}^2 \le C_\Omega. \qquad (53) $$

The classical Aubin's compactness argument then allows to pass to the limit $m \to \infty$, once again with the Lebesgue's and the Severini-Egorov's theorems for handling with the nonlinearities.

**Convergence 2: letting simultaneously $(h, \Delta t, m) \to (0, 0, \infty)$**

The reader may check that we can follow the lines developed for (43)–(48) to obtain the analogous results as $(h, \Delta t, m) \to (0, 0, \infty)$, namely, up to a subsequence:

$$ \begin{aligned} & \nabla_h \phi_h \to \nabla \phi \text{ weakly in } (L^2(\Omega_T))^d, \\ & c_h^m \to c \text{ and } \mu_h^m \to \mu \text{ in } L^2(\Omega_T) \text{ and a.e in } \Omega_T, \\ & \nabla_h c_h^m \to \nabla c \text{ and } \nabla_h \mu_h^m \to \nabla \mu \text{ weakly in } (L^2(\Omega_T))^d, \\ & \partial_t \bar{c}_h^m \to \partial_t c \text{ and } \partial_t \bar{\mu}_h^m \to \partial_t \mu \text{ weakly in } L^2(\Omega_T). \end{aligned} $$

Hence, passing to the limit as $(h, \Delta t, m) \to (0, 0, \infty)$ in (49)–(51) is the same work than as $(h, \Delta t) \to (0, 0)$.

**Convergence 3: letting $m \to \infty$ then $(h, \Delta t) \to (0, 0)$**

With the Convergence 2 at hands, it is sufficient to prove that the limit of $(c_h^{m,n}, \mu_h^{m,n})$ as $m \to \infty$ exists for any given $(h, \Delta t)$. First, the sequence $(p_h^{m,n})_m$ is uniformly bounded since all its elements belong to $[0, \bar{p}]$. Next, turning back once again to the estimates in Proposition 4.1, we consider them in the case where $h$ and $\Delta t$ are fixed values. Due to the regularity assumption on the mesh, the $L^2$ estimates imply that the sequences $(c_h^{m,n})_m$ and $(\mu_h^{m,n})_m$ are also uniformly bounded in $\mathbb{R}$. Hence there exist $c_h^n$ and $\mu_h^n$ such that (once again we do not relabel the subsequences), for any fixed $h$ and $n$

$$ c_h^{m,n} \to c_h^n \quad \text{and} \quad \mu_h^{m,n} \to \mu_h^n \quad \text{as} \quad m \to \infty $$

where $c_h^n$ and $\mu_h^n$ are defined by the following strongly coupled scheme: for any $n \in [0, N_T] \cap \mathbb{N}$, for any $K \in \mathcal{T}$

$$ \begin{aligned} c_K^0 &= \frac{1}{|K|} \int_K c_0(x) x, \quad - \sum_{L \in N(K)} |\sigma_{K,L}| \frac{\phi_L - \phi_K}{d_{K,L}} = |K|\theta_K + |\partial K \cap \Gamma|v_1, \\ R\psi|K| \frac{c_K^{n+1} - c_K^n}{\Delta t} - \sum_{L \in N(K)} |\sigma_{K,L}| \left( \psi S_m \frac{c_L^{n+1} - c_K^{n+1}}{d_{K,L}} + \frac{\phi_L - \phi_K}{d_{K,L}} c_K^{n+1} \right) \\ &= -|K| \left( r(c_K^n) - (\partial_p f)^{-1}(\mu_K^{N_T-n}) \right) + |\partial K \cap \partial \Omega|f_1, \end{aligned} $$
$$R\psi|K|\frac{\mu_K^{n+1} - \mu_K^n}{\Delta t} - \sum_{L \in N(K)} |\sigma_{K,L}| \left( \psi S_m \frac{\mu_L^{n+1} - \mu_K^{n+1}}{d_{K,L}} - \frac{\phi_L - \phi_K}{d_{K,L}} \mu_K^{n+1} \right) = -|K| \left( r' (c_K^{N_T-n}) + R\psi\rho + \theta_K \right) \mu_K^n + |K| \partial_c D (c_K^{N_T-n}) + |\partial K \cap \partial \Omega|f_2.$$

Because of the strong coupling, this scheme has a strongly implicit form. But in any case, we can also pass to the limit $m \to \infty$ in the estimates of Proposition 4.1. We get

$$\begin{aligned} \|c_h^n\|_{L^2(\Omega)}^2 &\le C_\Omega, & \|c_h^n\|_{H_h^1(\Omega)}^2 &\le C_\Omega, & \|\partial_t \bar{c}_h\|_{L^2(0,T;L^2(\Omega))}^2 &\le C_\Omega, \\ \|\mu_h^n\|_{L^2(\Omega)}^2 &\le C_\Omega, & \|\mu_h^n\|_{H_h^1(\Omega)}^2 &\le C_\Omega, & \|\partial_t \bar{\mu}_h\|_{L^2(0,T;L^2(\Omega))}^2 &\le C_\Omega. \end{aligned}$$

Furthermore, $(\partial_t f)^{-1}(\mu_K^{N_T-n}) \in [0, \bar{p}]$. The analogous of Theorem 4.4 may be retrieved and the arguments developed for passing to the limit $(h, \Delta t) \to (0, 0)$ in the Convergence 1's setting may be reproduced. Due to the result of Convergence 2, there is a subsequence that satisfies

$$\begin{aligned} & \nabla_h \phi_h \to \nabla \phi \text{ weakly in } (L^2(\Omega_T))^d, \\ & c_h^n \to c \text{ and } \mu_h^n \to \mu \text{ in } L^2(\Omega_T) \text{ and a.e in } \Omega_T, \\ & \nabla_h c_h^n \to \nabla c \text{ and } \nabla_h \mu_h^n \to \nabla \mu \text{ weakly in } (L^2(\Omega_T))^d, \\ & \partial_t \bar{c}_h^n \to \partial_t c \text{ and } \partial_t \bar{\mu}_h^n \to \partial_t \mu \text{ weakly in } L^2(\Omega_T). \end{aligned}$$

Theorem 5.1 is proved. ■

## 6 | A MANGASARIAN'S TYPE RESULT FOR TURNING BACK TO THE CONTROL PROBLEM

The latter section ensures the convergence of the numerical scheme (16)–(20). Of course, it remains to find the link between the limit function $p = (\partial_t f)^{-1}(\mu)$ and the optimal control $p^*$ solution of our original problem (7). It is given in the following theorem.

**Theorem 6.1** (characterization of the limit). Assume (H1)–(H7). Assume the continuity of the derivatives $r'$ and $\partial_c D$, with $r'$ bounded and $\partial_c D$ with a at most affine growth with respect to $c$. We have the following results:

(i) Assume further that the function $r$ is concave and that the cost function $D$ is increasing with the pollutant concentration $c$. Then the limit $p$ of $(p_h^{m,n})$ built from (16)–(20) as $m \to \infty, N_T \to \infty, h \to 0$, solves the optimal control problem (7).

(ii) If one of the functions $-r, -f$ and $D$ is strictly convex, $p$ is the unique solution of (7) and all the sequence $(\phi_h, c_h^{m,n}, \mu_h^{m,n}, p_h^{m,n})$ converges in $L^2(\Omega_T)$ and almost everywhere in $\Omega_T$.

*Proof.* We turn back to the maximization problem. Recall $(c, p)$, defined in Theorem 5.1, satisfies the Euler–Lagrange's problem (8)–(12). Let $q \in E$. We denote by $c_q$ the concentration solution of the state equation associated with the control $q$, Using the convexity properties of $f$ and $-D$, we compute
$$
\begin{aligned}
J(p) - J(q) &= \int_0^T \left( \int_{\mathcal{S}} (f(x, p(t, x)) - f(x, q(t, x))) dx - \int_{\Omega} \left( D(x, c(t, x)) \right. \right. \\
&\qquad \left. \left. - D(x, c_q(t, x)) dx \right) e^{-\rho t} dt - \zeta e^{-\rho T} \int_{\Omega} \left( D(x, c(T, x)) - D(x, c_q(T, x)) \right) dx \right. \\
&\geq \int_{\Omega_T} \left( \chi_{\mathcal{S}} \partial_{\rho} f(x, p)(p-q) + \partial_c D(c) (c_q - c) \right) dx e^{-\rho t} dt \\
&\qquad + \zeta e^{-\rho T} \int_{\Omega} \partial_c D(c(T, x)) (c_q(T, x) - c(T, x)) dx.
\end{aligned}
$$

Bearing in mind that $\partial_{\rho}f(x, p) = \mathcal{J}_T\mu$ and $\partial_c D(c(T, x)) = R\psi\mu(0, x)$, that $\partial_c D(c)$ is given by the adjoint equation (11), integrating by parts, we obtain

$$
\begin{aligned}
J(p) - J(q) &\geq \int_{\Omega_T} \chi_{\mathcal{S}} \mathcal{J}_T \mu (p-q) e^{-\rho t} dxdt - \int_{\Omega_T} \mathcal{J}_T \mu (R\psi \partial_t (c_q - c) \\
&\quad + \operatorname{div} (v_{c_q} - v_c) - r'(c) (c_q - c)) e^{-\rho t} dxdt - \int_{\Omega} R\psi (\mathcal{J}_T \mu)(T,x) (c_q(T,x) \\
&\quad - c(T,x)) e^{-\rho T} dx + e^{-\rho T} \int_{\Omega} R\psi \mu(0,x) (c_q(T,x) - c(T,x)) dx \\
&= - \int_{\Omega_T} \mathcal{J}_T \mu (r(c) - r(c_q) - r'(c) (c_q - c)) e^{-\rho t} dxdt.
\end{aligned}
\tag{54}
$$

We thus study the sign of the last term in (54). Setting $\mu^- = \max\{0, -\mu\}$, we multiply the adjoint equation (11) by $\mu^-$ and we integrate by parts over $\Omega_T$. Since $R\psi\mu^-|_{t=0} = \zeta(\partial_c D(c|t=T))^- = 0$, we get

$$
\begin{aligned}
& \frac{1}{2} \int_{\Omega} R\psi |\mu^{-}|^2 dx + \int_{\Omega_T} S_m \psi |\nabla \mu^{-}|^2 dxdt + \int_{\Omega_T} \partial_c D(\mathcal{J}_{TC}) \mu^{-} dxdt \\
& = - \int_{\Omega_T} \mu \mathcal{J}_{T\nu} \cdot \nabla \mu^{-} dxdt - \int_{\Omega_T} (r'(\mathcal{J}_{TC}) + R\psi\rho + \mathcal{J}_{T\theta}) |\mu^{-}|^2 dxdt.
\end{aligned}
$$

Here, since $\operatorname{div}\nu \in L^\infty(\Omega_T)$ and $\operatorname{curl}\nu = 0$, the velocity belongs to $L^\infty(\Omega_T)$. Since moreover $r'(\mathcal{J}_{TC}) + R\psi\rho + \mathcal{J}_{T\theta} \in L^\infty(\Omega_T)$, using the Cauchy-Schwarz and the Young inequalities, we infer from the latter relation that, if $\partial_c D \ge 0$ (that is if the cost function $D$ is increasing with the pollutant concentration $c$, which seems like a reasonable choice),

$$ \frac{1}{2} \int_{\Omega} R\psi |\mu^{-}|^2 dx + \int_{\Omega_T} \frac{S_m}{2} \psi |\nabla \mu^{-}|^2 dxdt \le C \int_{\Omega_T} R\psi |\mu^{-}|^2 dxdt. $$

The Gronwall lemma then ensures that $\mu^- = 0$, that is $\mu \ge 0$ almost everywhere in $\Omega_T$. Turning back to (54), we conclude that if $r$ is a concave function then

$$ J(p) \ge J(q) \quad \text{for any } q \in E, $$

meaning that the limit value $p$ obtained by the scheme (16)–(20) necessary solves the optimal control problem (7). Furthermore, if one of the functions $-f$, $D$, $-r$ is strictly convex, we get $J(p) > J(q)$ for any $q \in E$ meaning that the solution of (7) is unique. It implies that there is no need to extract subsequences in the previous subsection. This ends the proof of Theorem 6.1. ■
7 | NUMERICAL EXPERIMENTS

In this part, some numerical results are given for illustrating the efficiency of the scheme (16)–(20). The scheme was implemented in the Python language using the sparse matrices defined in the Scipy module. This choice was guided by the impossibility of representing the different systems to be solved in a dense way in memory, especially in the three-dimensional case. For a reasonable discretization (up to $2 \cdot 10^4$ unknowns), the simulations was run on a laptop computer (2.3 GHz processor, 16 GB memory). For more important simulations ($\emptyset (10^5)$ unknowns), we deported the calculations on a modest computer (Xeon E5-2630 quadcore 2.60GHz, 32 Go). The codes were neither optimized nor parallelized. For the examples below, the calculations took about 12 hours in the 3D case, a few hours in the 2D case. We also checked the stability of the results with mesh refinement for comforting this discretization choice. By the way, we are currently studying the possibility of carrying out our simulations on GPUs (many libraries exist but are not always stable) in order to significantly improve the speed of execution and thus increase the discretization.

First notice that for all the simulations the time and space steps for the discretization are given by $\Delta t = h = 10^{-2}$. The number of loops (that is the number of increments $m$ in (20)) necessary to arrive at the stopping criterion is denoted by $M$. Namely $M = \min_{m \in \mathbb{N}} \{m + 1\}$, such that $\|p_h^{m+1} - p_h^m\|_{L^2(0,T;L^2(\Omega))} \le \varepsilon$ where $\varepsilon$ is a given threshold. We used $\varepsilon = 5.10^{-6}$.

The aquifer is simply represented by a parallelepiped, namely $(x, y) \in ]0,900[^2$, $z \in ]-9,0[$, the values are given here in meters. We consider such a shallow aquifer because most of the groundwater reservoirs are thin geological formations. The aquifer is rescaled into the unit cube for the readability. More precisely, we introduce reference quantities, hereafter denoted with the subscript $R$ for scaling the model. For instance, $c_R$ is the reference concentration such that $c = c_R \hat{c}$. The surface length scale is $x_R = 900$, the vertical length scale is $z_R = 9$, the time scale is $t_R$. We set $\hat{x} = x/x_R \in (0, 1)$ and $\hat{y} = y/x_R \in (0, 1)$, $\hat{z} = z/z_R + 1 \in (0, 1)$ and $\hat{t} = t/t_R \in (0, T/t_R)$. We skip the details. The rescaled version of equation (9) reads for instance

$$R_R \hat{R} \psi_R \hat{\psi} c_R \partial_t \hat{c} - \frac{t_R}{x_R^2} c_R \operatorname{div} \hat{v}_{\hat{x}, \hat{y}} \left( S_{m,R} \hat{S}_m \psi_R \hat{\psi} \nabla_{\hat{x}, \hat{y}} \hat{c} + v_{R,x,y} \hat{c} \hat{v}_{x,y} \right) \\ - \frac{t_R}{\zeta_R^2} c_R \partial_z \left( S_{m,z,R} \hat{S}_{m,z} \psi_R \hat{\psi} \partial_z \hat{c} + v_{R,z} \hat{c} \hat{v}_z \right) = -r_R \hat{r}(\hat{c}) + \chi \hat{\mathcal{F}} p_R \hat{p}^*$$

in $(0, T/t_R) \times (0, 1)^3$. The discretization with $h = 10^{-2}$ (resp. $\Delta t = 10^{-2}$) is then used in $[0, 1]^3$ (resp. in $(0, T/t_R)$).

A water production well is located on the parcel, at the point of coordinates (450; 450). It penetrates 2/5 of the thickness of the aquifer. Furthermore, we assume that the water is pumped all along the borehole and not only in deep areas that could be for a while unaffected by pollution. The water withdrawal in production wells is thus in the form $q^P = 100\beta\chi_{\text{wells}}$. We will use this expression in the definition of the cost function below. The parameter $\beta$ is introduced for handling different pumping rates. Our simulations correspond to a production rate in the range $10^3-10^4 \text{ m}^3 \text{ h}^{-1}$. For the sake of the simplicity we assume that the recharge occurs in all the domain $\Omega$. Globally the balance between outputs and inputs does not put the reservoir at risk: $\theta = q^I - q^P$ is assumed positive, in average of the order $10^{-3} \text{ m}^3 \text{ h}^{-1}$. Such values correspond in particular to French data (see [17]). The physical parameters have been chosen in agreement with the classical literature (see for instance Bear [4]). The soil characteristics are $\kappa = 39.04 \text{ m} \cdot \text{day}^{-1}$ and $\psi = 0.3$. In such a shallow aquifer the fluid displacement is essentially horizontal. Here moreover, for ensuring the interpretability of the results, the boundary conditions for the fluid velocity equation (1) are nonhomogeneous Neumann boundary conditions (scaled by parameters $v_\ell$ for the inlet side and $v_r$ for the outlet side) such that the quasi-horizontal
flow goes from left to right (see Figure 1) with a Darcy velocity of about 600 m year⁻¹. Once the inlet and outlet flow are known, the boundary conditions (4) for the concentration are of Robin's type. We choose $f_1$ such that they correspond to convective conditions, namely $f_1 = v_x c_{ext}$, $x = \ell$ or $r$, $c_{ext}$ being initialized at the value of the natural concentration, that is before pollution. This latter value is also used for $c_0$.

The water molecular diffusion is $S_m = 9.4 \cdot 10^{-8} \text{ m}^2 \text{ s}^{-1}$. It follows that the Péclet number is of order 600. For such a value, the dispersion cannot be neglected. Nevertheless, it is known to be well approximated by the Fickian hypothesis. In the present shallow aquifer with essentially horizontal water displacement, only the longitudinal dispersion is dominant. The last two points lead to consider a diffusion tensor rather than a diffusion coefficient, of the form:

$$S_m = \begin{pmatrix} S & 0 & 0 \\ 0 & S & 0 \\ 0 & 0 & S_m \end{pmatrix},$$

where $S = 10^3 S_m$ to account for dispersive effects (see [15] and references therein). One easily checks that using the matrix $S_m$ instead of a real number $S_m$ in the theoretical analysis do not change the results obtained in the previous sections.

We choose to illustrate the problem of agricultural pollution due to nitrates fertilization. Nitrates are almost not adsorbed by the soil and we thus set $R = 1$. For the reaction term $r$, we use the function $r(c) = \gamma_r c^2$ with $\gamma_r = 10^{-3}$. The corresponding Damköhler number is thus of order $5 \cdot 10^{-3}$. Small values of the Damköhler number, that is smaller than one, correspond to the fact that nitrate removal within the groundwater is controlled by water flux rather than reaction rate. A test using $\gamma_r = 10^{-1}$ is also performed. The initial condition is $c_0 = 5 \text{ mg L}^{-1}$. The benefit function is defined by

$$f(p) = \alpha \begin{cases} K_1 - K_2 K_3 10^{-4} e^{-K_3 \bar{p}} e^{-10^4 p} - K_2 + K_2 K_3 10^{-4} e^{-K_3 \bar{p}} & \text{if } p < 0, \\ K_1 - K_2 e^{-K_3 p} & \text{if } p \in [0, 2\bar{p}], \\ K_1 - K_2 K_3 10^{-4} e^{-2K_3 \bar{p}} (e^{10^4(2\bar{p}-p)} - 1) - K_2 e^{-2K_3 \bar{p}} & \text{if } p > 2\bar{p}, \end{cases}$$

with $K_1 = 11.7888$, $K_2 = 8.6 \cdot 10^{-3}$, $K_3 = 50.465$ and $\bar{p} = 1.5$. It comes from Godart et al. [13] where the crop yield is depending on the yield value without nitrogen input and the asymptotic value when the input becomes important. Values of the parameters depend on the crop species. Here we choose the example of wheat crops. Notice that the Godart function has been modified for ensuring that the

<figure><img src="image_3.png" /><figcaption>FIGURE 1 Representation of the Darcy velocity $v$ in the aquifer rescaled into the cube $[0, 1]^3$ for the readability</figcaption></figure>
fertilizer load $p$ remains in the interval $[0, \bar{p}]$ (using loads outside this range is highly penalized). The cost function, on the other hand, is given by

$$D(c) = 100\beta c^2 \chi_{\text{wells}}.$$

The real number $\alpha$ depends on the price of cultivated species per ton and $\beta$ accounts for the production flow rate of the well and cost of treatment per cubic meter. The real numbers $\alpha$ and $\beta$ in the definition of $f$ and $D$ may also be used for possibly tuning the weights given to the economic and environmental concerns in the model. For instance, choosing $\alpha \ll 1$ and $\beta \approx 1$ allows to privilege the quality of the water. In all cases we used reasonable parameters to seek the trade-off between the polluter's benefits and the water treatment cost. Thus, for example, in all the simulation results there is indeed fertilizer application but the immediate vicinity of the water production well is systematically protected. The other parameters in the functional $J$ are $\zeta$ and $\rho$. We fix $\rho = 0.05$. The value of $\zeta$ was chosen in the set $\{10^{-2}, 10, 100\}$ for giving more or less weight to the legacy of pollution.

The results of our first experiments are represented in Figures 2 and 3. They were especially designed to verify that the code gives intuitively logical results. But for decision support, a rigorous calculation is sometimes more convincing than an intuition. The performance of the iterative algorithm is illustrated by the second line of figures with in all cases the maximum number of tests $M$ such

<figure><img src="image_4.png" /><figcaption>FIGURE 2 Influence of the weights on the “polluter” profit $\alpha$, and on the groundwater protection $\beta$: optimal fertilizer load on the line $x = 450$ as a function of time (top). The blue lines figure the boarder of the well. The quality of the convergence is illustrated by the decrease of the ‘error’ $\|p_h^{m+1} - p_h^m\|_{L^2(0,T;L^2(\Omega))}$ represented in log scale (middle) and by the increase of the objective function $J$ (bottom). The simulation duration is $T = 100$ days. The weight on the pollution legacy is $\zeta = 1$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIGURE 3 Influence of the weight $\zeta$ on the pollution legacy: optimal fertilizer load on the line $x = 450$ as a function of time. The simulation duration is $T = 100$ days</figcaption></figure>

that $M \le 50$. The greater the ratio $\beta/\alpha$, the more a policy of water resource protection is favored. The simulations results given in Figure 2 respect this point. They also show that the owner of the left parcel is the most penalized. This can be explained by the direction of the groundwater velocity, from left to right in our figures. Thus, everything that is spread on the left part is transported to the water production well. This penalty becomes less severe over time because what is applied in the last few days does not have time to reach the well before the end of the study period. In the last few days there is even a patch of high concentration (around $y = 16$ decameters) indicating that the farmer can afford to compensate for the lack of application during the first few days. In Figure 3, we focus on the influence of the prize of the pollution legacy by varying the parameter $\zeta$. Once again, the results are qualitatively in line with what one would expect. However, we note that the legacy cost may be as important as the immediate cleanup costs (the reader may also compare the results in Figures 2 and 3 where the pollution legacy weight was set as $\zeta = 1$). Such results make sense of a taxation policy that is uncorrelated with the immediate cost of water treatment.

For all remaining simulations the parameter $\zeta$ was set to the value 1. The duration of the experiment in Figures 2–3 is $T = 100$ days which corresponds to the fertilizer application period for most cereal crops. In Figure 4 however, the simulations are performed for a longer duration, namely for $T = 3000$ days. First of all, there is a great stability in time of the quantity of fertilizer to be spread
$$\alpha = 1, \beta = 30$$

$$\alpha = 1, \beta = 100$$

$$\alpha = 0.1, \beta = 100$$

<figure><img src="image_6.png" /><figcaption>FIGURE 4 Example of a longer simulation: the duration is $T = 3000$ days. Optimal fertilizer load on the line $x = 450$ as a function of time: during the complete time interval $(0, T)$ (top line), between day 2980 and day 2990 (second line), during the last 10 days (third line). The error as a function of $m \in \{0, \dots, M\}$ is plotted at the bottom</figcaption></figure>

and the areas concerned. The only notable changes occur during the last 10 days of spreading, a wind-fall effect since the pollution will not have time to reach the water extraction well. Obviously we also notice that the spreading is totally forbidden on the left plot, except during the last 10 days. The optimal solution in long time is thus very different from the one in shorter time. From the point of view of the algorithm on the other hand, the convergence is not at all impacted by the change of scale.
So far we have taken advantage of the symmetries of the problem to make only 2D representations. In order to convince the reader of the validity of this approach, two figures finally give the 3D representation of the aquifer. Figure 5 provides an opportunity to verify that the essential mechanism controlling the nitrate removal is the water flux rather than the reaction rate (bear in mind that the flow is oriented from the left side of the cube to the right side; as $v_ℓ = v_r$, the smaller $v_ℓ$ is, the higher the water velocity is). This is even clearer when looking at the nitrate concentration of the input (note the colorbar in Figure 6). There is another way to analyze these results. Indeed, we have chosen $α = 1$, $β = 100$ and yet the owner of the plot on the right is not particularly advantaged (compared to the 2D

$$ \gamma_r = 0.001, v_ℓ = v_r = 0.1 \qquad \gamma_r = 0.1, v_ℓ = v_r = 0.1 \qquad \gamma_r = 0.1, v_ℓ = v_r = 0.01 $$

<figure><img src="image_7.png" /><figcaption>FIGURE 5 Three-dimensional representation of the optimal fertilizer load. The simulation duration is $T = 100$ days, the results are figured at the 1st, 20th, and 80th days. Changing the parameter $γ_r$ in the reaction term highlights the influence of the natural remediation of the soil. Changing the entry and exit velocities $v_ℓ$ and $v_r$ illustrates the influence of the water flux</figcaption></figure>

$$ \gamma_r = 0.001, v_ℓ = v_r = 0.1 \qquad \gamma_r = 0.1, v_ℓ = v_r = 0.1 \qquad \gamma_r = 0.1, v_ℓ = v_r = 0.01 $$

<figure><img src="image_8.png" /><figcaption>FIGURE 6 Same experiment as in Figure 5. Nitrate concentrations on day 20</figcaption></figure>
case presented in Figures 2–4). To arrive at this result we assumed for the 3D simulations that the geo-
logical nature of the soil to the right of the studied area is much less permeable to nitrate propagation.
For the mathematical model this is taken into account by weighting the Robin boundary condition (9)
so that it approaches a homogeneous Neumann boundary condition. For interpretation, this means that
the importance of the nature of the soils adjacent to the plots studied should not be forgotten, since it
can even balance out the effects of groundwater flow.

8 | CONCLUSION

We have constructed a numerical scheme for the optimal control of groundwater pollution. By proving its convergence, we have also proved the existence and uniqueness of a solution to the optimal control problem. Numerical simulations were provided. Looking at the log scale representations of the convergence curves, we can see the good performance of the algorithm with a convergence speed that approaches the geometric convergence of a fixed point algorithm.

We would like to point out that our work can be extended to many other applications. The non-
point sources of groundwater pollution linked to a profit are obviously numerous, agriculture, industry,
mining and so forth. On the other hand, environmental pollution has many facets, behind the obvious
cleaning costs that were considered in the present article. Health risk issues and thus cost of care could
be included in the objective functional (see [16]). More surprisingly and perhaps less philanthropi-
cally, work has already been developed to express the devaluation of property values as a function of
pollution levels (see [10]).

ACKNOWLEDGMENT

This work was partially supported by the New Aquitaine region and Europe (CPER ECONAT).

DATA AVAILABILITY STATEMENT

Data sharing not applicable to this article as no datasets were generated or analysed during the current
study

ORCID

Catherine Choquet https://orcid.org/0000-0003-0632-0358

REFERENCES

[1] E. Augeraud-Véron, C. Choquet, and É. Comte, *Existence, uniqueness and asymptotic analysis of optimal control problems for a model of groundwater pollution*, ESAIM: Control Opt. Calc. Var. 25 (2019), 1–20.

[2] E. Augeraud-Véron and M. Leandri, *Optimal pollution control with distributed delays*, J. Math. Econom. 55 (2014), 24–32.

[3] V. Barbu and M. Iannelli, *Optimal control of population dynamics*, J. Optim. Theory Appl. 102 (1999), 1–14.

[4] J. Bear and A. Verruijt, “Theory and applications of transport in porous media,” *Modeling of groundwater flow and pollution*, Reidel, Dordrecht, 1987.

[5] C. Benosman, B. E. Aïnseba, and A. Ducrot, *Optimization of cytostatic leukemia therapy in an advection-reaction-diffusion model*, J. Opt. Theory Appl. 167 (2015), 296–325.

[6] C. Bourgeois and P. A. Jayet, *Regulation of relationships between heterogeneous farmers and an aquifer accounting for lag effects*, Australian J. Agric. Res. Econom. 60 (2016), 39–59.

[7] W. Brock and A. Xepapadeas, *Pattern formation, spatial externalities and regulation in coupled economic-ecological systems*, J. Env. Econom. Manage. 59 (2010), 149–164.

[8] S. R. Carpenter, D. Ludwig, and W. A. Brock, *Management of eutrophication for lakes subject to potentially irreversible change*, Ecol. Appl. 9 (1999), 751–771.
[9] F. H. Clarke, *A new approach to Lagrange multipliers*, Math. Opt. Res. 1 (1976), 165-174.

[10] D. M. Dornbusch, S. M. Barrager: *Benefit of water pollution control on property values*, Socioeconomic Environmental Studies Series. Office of Research and Development U.S. Environmental Protection Agency Washington, D.C., EPA-600/5-73-005, 1973.

[11] R. Eymard, T. Gallouët, and R. Herbin, *Finite volume methods*, Handbook Numer. Anal. 7 (2000), 713-1018.

[12] M. Gander, F. Kwok, and J. Salomon, *PARAOPT: A parareal algorithm for optimality systems*, SIAM J. Sci. Comp. 42 (2020), A2773-A2802.

[13] C. Godart, P. A. Jayet, B. Niang, L. Bamiere, S. DeCara, E. Debove, E. Baranger: *Rapport Final, APR GICC 2002*, Convention de recherche MEDD num 02-00019 INRA-INAPG, 2002.

[14] M. Heinkenschloss and D. Leykekhman, *Local error estimates for SUPG solutions of advection-dominated elliptic linear-quadratic optimal control problems*, SIAM J. Numer. Anal. 47 (2010), 4607-4638.

[15] M. Icardi, G. Boccardo, D. L. Marchisio, T. Tosco, and R. Sethi, *Pore-scale simulation of fluid flow and solute dispersion in three-dimensional porous media*, Phys. Rev. E 90 (2014), 013032.

[16] D. T. Jamison, J. G. Breman, A. R. Measham, et al., Eds., *Disease control priorities in developing countries. 2nd edition. Washington (DC): The International Bank for Reconstruction and Development / the World Bank*, Oxford University Press, New York, 2006.

[17] E. Ledoux, E. Gomez, J. M. Monget, C. Viavattene, P. Viennot, A. Ducharne, M. Benoit, C. Mignolet, C. Schott, and B. Mary, *Agriculture and groundwater nitrate contamination in the seine basin, the STICS-MODCOU modelling chain*, Sci. Total Environ. 375 (2007), 33-47. https://doi.org/10.1016/j.scitotenv.2006.12.002.

[18] J. Liu and Z. Wang, *Non-commutative discretize-then-optimize algorithms for elliptic PDE-constrained optimal control problems*, J. Comput. Appl. Math. 362 (2018), 596-613.

[19] H. Nasser El Dine, *Étude mathématique et numérique pour le modèle Darcy-Brinkman pour les écoulements diphasiques en milieu poreux (Doctoral dissertation)*, École Centrale de Nantes, 2017.

[20] R. Winkler, *A note on the optimal control of stocks accumulating with a delay*, Macroecon. Dyn. 15 (2011), 565-578.

**How to cite this article:** C. Choquet, M. M. Diédhiou, H. N. El Dine, and C. Saint-Jean, *A numerical scheme for the optimal control of groundwater pollution*, Numer. Methods Partial Differ. Eq. 39 (2023), 3037-3063. https://doi.org/10.1002/num.22999