# Modeling and control of in-situ decontamination of large water resources

Pedro Gajardo, Jérome Harmand, Hector Ramirez Cabrera, Alain Rapaport, Victor Riquelme, Antoine Rousseau

► To cite this version:

Pedro Gajardo, Jérome Harmand, Hector Ramirez Cabrera, Alain Rapaport, Victor Riquelme, et al. Modeling and control of in-situ decontamination of large water resources. ESAIM: Proceedings and Surveys, 2017, ETAMM 2016 - Emerging Trends in Applied Mathematics and Mechanics, 57, pp.70-85. 10.1051/proc/201657070 . hal-01459969v2

HAL Id: hal-01459969

https://inria.hal.science/hal-01459969v2

Submitted on 28 Feb 2017

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>
MODELING AND CONTROL OF IN-SITU DECONTAMINATION OF LARGE
WATER RESOURCES

PEDRO GAJARDO<fnref n="1" />, JÉRÔME HARMAND<fnref n="2" />, HÉCTOR RAMIREZ<fnref n="3" />, ALAIN RAPAPORT<fnref n="4" />,
VICTOR RIQUELME<fnref n="5" /> AND ANTOINE ROUSSEAU<fnref n="6" />

**Abstract.** We address the problem of the optimal control of *in situ* decontamination of water resources. We review several modeling, simulation and optimization techniques for this problem and their results. We show the benefit of combining tools from finite dimensional optimal control theory and numerical simulations of hydrodynamics equations, for providing simple and efficient feedback strategies.

**Résumé.** Nous considérons le problème de commande optimale d'une décontamination *in situ* de ressources hydriques. Nous passons en revue différentes techniques de modélisation, simulation et optimisation pour ce problème, ainsi que leurs résultats. Nous montrons le bénéfice à combiner les outils de la théorie de la commande optimale en dimension finie avec la simulation numérique des équations de l'hydrodynamique, afin de fournir des stratégies par retour d'état à la fois simples et efficaces.

INTRODUCTION

Water decontamination is one of the world largest issues due to scarcity of fresh water and access to drinking water. Industry of wastewater treatment has developed several technologies based on filtration, biological conversion, and the so-called tertiary treatments for the remediation of contaminated waters, depending on the pollution nature (accidental or regular) and former nature of the water (domestic, agricultural or industrial). In almost all cases, water decontamination is achieved in dedicated plants called wastewater treatment plants or WWTP: wastewaters are collected and sent to such plants while treated waters are released into the environment. For (possibly large) natural water resources such as lakes, lagoons, ponds, water tables, it is highly not desirable nor often even possible to empty reservoirs for treating contaminated waters in waste-water treatment plants (WWTP). Refilling with fresh water while the reservoir is emptied for the treatment of contaminated water could be an alternative, but not realistic or too expensive when dealing with very large volumes. Up to now, very few life-size depollution experiments of large water resources have been reported: very often, the experimenters

<footnote n="1">¹ Departamento de Matemática, Universidad Técnica Federico Santa María, Valparaíso, Chile; pedro.gajardo@usm.cl</footnote>

<footnote n="2">² LBE, INRA, 11100, Narbonne, France; jerome.harmand@inra.fr</footnote>

<footnote n="3">³ Departamento de Ingeniería Matemática and Centro de Modelamiento Matemático (UMI 2807, CNRS), Universidad de Chile.
Casilla 170-3 Santiago 3, Chile; hramirez@dim.uchile.cl</footnote>

<footnote n="4">⁴ UMR INRA/SupAgro 729 'MISTEA', Montpellier, France; alain.rapaport@inra.fr</footnote>

<footnote n="5">⁵ Departamento de Ingeniería Matemática and Centro de Modelamiento Matemático (UMI 2807, CNRS), Universidad de Chile.
Casilla 170-3 Santiago 3, Chile; vriquelme@dim.uchile.cl</footnote>

<footnote n="6">⁶ Inria team LEMON, Inria Chile and UMR 5149 France; antoine.rousseau@inria.fr</footnote>
are left to deal with pollution at the source, take measures to limit its effects (limit the use of chemicals, recycle waters in industries, save and minimize the water needs, multiply the number of WWTPs in the corresponding watershed...) and, in any case, treat systematically waters before being discharged into the environment, leaving "to the natural lifecycle the task of dealing with the pollution". Such a "solution" may take (possibly very long) time and the result is far from being guaranteed.

While there are a number of issues related to the use of living organisms to remove certain chemicals from waters, for instance some microorganisms are used to fight against oil spills, real life experiments are rare because we usually have little knowledge of possible consequences of introducing a microorganism in a medium from which it is naturally absent. This has led to the development of sometimes very detailed models such as for the Chinese Taihu lake for better understanding the ecology of the lake and trying to think of actions to restore its quality, [11]. Moreover, to the best of authors knowledge, such a solution has been rarely addressed in the literature in terms of simulation and optimization that very often rather used to identify and quantify pollution sources but not often to address water depollution at the resource level. Related to lake management, one of the only cases which is well documented is related to the use of brick powder (which plays the role of a nitrate-storing compound) to reduce phosphorus release from sediments and thus decrease the eutrophication of lakes (cf. instance [3] or [1] for recent reviews on the subject). Another field where the use of a recirculation system makes it very similar to that considered in the present paper is aquaculture. Indeed, fisheries make use of recirculation systems to optimize feedstock. After treatment, water with resting food is recirculated. However, in such systems, most of efforts are put on the potential recirculation of pathogens and pollutant removal is not the main consideration (cf. for instance [8]).

The idea developed in this paper is completely different: facing a contaminated water resource, an alternative could be to pump water from the reservoir, to treat it in dedicated and appropriate WWTPs (depending on the chemicals to be removed from water), and to release decontaminated water to the reservoir. This implies the use of a continuous water treatment system, typically in a (set of) process tank(s) called bioreactor(s).

Coupling a continuous decontamination process with the hydrodynamics of the water resource brings some issues, and especially the following ones.

(1) There is a trade-off in the choice of the flow rate: a large flow rate allows to quickly treat all the water of the reservoir but the treatment is less efficient (as the residence time in the process tank is short), while a small flow rate allows to refill a water of better quality but slowly.

(2) For large water reservoirs, the pollutant concentration is not spatially homogeneous and the locations of pumping and injection points may impact the hydrodynamics and thus the performance of decontamination.

Therefore, a central question is to decide which flow rate is the best for a decontamination objective. For large water reservoirs, which might require weeks or months to be decontaminated, an adequate objective is to minimize the time duration of the treatment process to obtain a pollutant concentration below a threshold. More specifically, one may wonder if a time-varying pumping strategy could have a significant impact on the total decontamination time, compared to a constant pumping strategy, and then how to adjust dynamically the flow rate. Such kind of questioning falls precisely in the field of automatic control, but with certain characteristics:

- If one considers the water resource as an “input-output” system (inputs are flow rate and pollutant concentration of the treated water and outputs are flow rate and pollutant concentration at the pumping location), the response of the system is rarely very well known and relies on complex infinite dimensional modeling of the hydrodynamics of the resource.

- On-line information is relatively poor: sensors could measure pollutant concentrations at some punctual locations but the accurate knowledge of the spatial distribution of the pollutant in the whole resource at any time is not accessible.

- The decontamination process imposes a (non-linear) constraint between the input flow rate and the pollutant concentration of the treated water.

Therefore, one may look for a feedback or adaptive strategy that could adjust the flow rate depending on on-line measurements of pollutant concentrations in order to guarantee short duration of the decontamination process.
Moreover, one may also want to investigate the benefit of having several measurement locations and/or several pumping or rejection locations.

The aim of the present paper is to review

- the challenges and the expected benefits of decision support tools based an mathematical modeling, simulations and optimization,

- the existing results and future possible extensions.

The paper is organized in six sections:

(1) Formulation and study of an optimal control problem for the choice of the best pumping strategies under a simplifying hypothesis of a perfectly homogeneous water resource

(2) Extensions of the former study to some simple spatial patterns of heterogeneity in the water resource

(3) Accurate simulations of spatial hydrodynamics under variable input flows and the decontamination constraint

(4) Use of feedback strategies obtained on simplified model with accurate simulations, and study of their performance

(5) Consideration of several pumps

(6) Future expected developments and extensions

# 1. HYPOTHESES, NOTATIONS AND FORMULATION OF THE PROBLEM

The modeling problem consists in coupling two dynamics:

(1) of the treatment unit, which receives a flow $Q_{in}$ of contaminated water from the resource to be decontaminated, and produces as output a flow $Q_{out}$ of treated water that returns to the resource,

(2) of the water resource that receives decontaminated water with flow $Q_{out}$ from the treatment unit, and from which a flow $Q_{in}$ is pumped to be decontaminated by the treatment unit.

A specificity of the problem is that one wants to preserve a constant volume $V$ of the resource at any time. Therefore this imposes to have $Q_{in} = Q_{out} := Q$. Then, the treatment unit consists in a continuous flow bioreactor of constant volume $V_T$ followed by a settler, that separates the decontaminated water from the biomass that exits from the bioreactor.

We shall assume the settler to be perfect, in the sense that it does not reduce the output flow $Q_{out}$ and that no biomass at all goes to the resource. This last point is quite important from the ecological and modelling viewpoints. A presence of biomass in the resource would participate to the decontamination of the resource but also to its eutrophication, consuming oxygen no longer available for other sources of life (fishes, algae...) and a proliferation of undesired biomass.

The interconnection scheme is summarized in Figure 1.

If one assumes that the pollutant concentration in the resource is homogeneous, the model that couples the time evolution of the pollutant concentrations ($S_T$) in the treatment unit and ($S$) in the resource can be written as follows:

$$ \dot{B} = \mu(S_T)B - \frac{Q}{V_T}B \quad (1a) $$

$$ \dot{S}_T = -\frac{\mu(S_T)B}{Y} + \frac{Q}{V_T}(S - S_T) \quad (1b) $$

$$ \dot{S} = \frac{Q}{V}(S_T - S) \quad (1c) $$

where $B$ represents the biomass concentration inside the bioreactor. The first two equations come from the classical (perfectly mixed) chemostat model (see a.e. [9]), where $\mu(\cdot)$ and $Y$ are respectively the growth rate function and the yield conversion factor, both specific to the bacterial strain that grows in the bioreactor.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Scheme of the interconnection between the resource and the decontamination unit (bioreactor + settler).</figcaption></figure>

We shall assume that the growth rate $\mu(\cdot)$ is monotonic, as it is often the case for a reasonable range of substrate concentrations, and define the usual break-even concentration as follows:

$$ \lambda(D) = \begin{cases} +\infty & \text{if } \{S > 0 \text{ s.t. } \mu(S) > D\} = \emptyset \\ \inf \{S > 0 \text{ s.t. } \mu(S) > D\} & \text{otherwise} \end{cases} $$

A typical instance of monotonic growth rate function is given by the Monod expression

$$ \mu(S) = \frac{\mu_{\max}S}{K_s + S}. \qquad (2) $$

When the volume of the bioreactor $V_T$ is very small compared to the volume $V$ of the resource, the first two equations can be approximated by the quasi-stationary asymptotically stable steady state $(B^{eq}, S_T^{eq})$ that is defined by the equations

$$ S_T^{eq}(Q) = \min(S, \lambda(Q/V_T)), \quad B^{eq}(Q) = Y(S - S_T^{eq}(Q)). $$

The interested reader is referred to [9] for the mathematical analysis of the dynamics and equilibria of the chemostat model. The reduced model is thus

$$ \dot{S} = \frac{Q}{V}(S_T^{eq}(Q) - S). \qquad (3) $$

When the flow rate $Q$ is larger than $\bar{Q}(S) = V\mu(S)$, one has $S_T^{eq} = S$ which amounts to state that the bioreactor is not converting any substrate at all for those values of $Q$ (washout regime). Such situations are clearly not of interest. So the flow rate $Q$ will be chosen in the interval $[0, \bar{Q}(S)]$ and the optimization problem can be stated as follows.

**Problem 1.** For a given threshold $\underline{S} > 0$ and an initial concentration $S^0$ in the resource, we look for a time varying flow rate $Q(\cdot)$ that realizes the minimization

$$ T_{\underline{S}}^{\text{opt}}(S^0) := \inf_{Q(\cdot) \in \mathcal{Q}(S^0)} T_{\underline{S}}(S^0, Q(\cdot)) $$
<figure><img src="image_3.png" /><figcaption>(A) Consideration of a simple gradient between input and (B) Consideration of active and dead zones (connected) in output locations in the resource. The resource.</figcaption></figure>

<caption>FIGURE 2</caption>

where $T_S(\cdot)$ is the time to reach the target:

$$T_S(S^0, Q(\cdot)) = \inf \{t \ge 0 \mid S^{S^0, Q(\cdot)}(t) = S\}.$$

Here, $S^{S^0, Q(\cdot)}(\cdot)$ denotes the solution of (3) with initial condition $S^0$ and time varying control $Q(\cdot)$. Optimal controls $Q(\cdot)$ are sought in the set

$$Q(S^0) = \{ Q : \text{measurable s.t. } Q(t) \in [0, \bar{Q}(S^{S^0, Q(\cdot)}(t))], \forall t \in [0, T_S(S^0, Q(\cdot))] \}.$$

As $\mu(\cdot)$ is monotonic, there is a one to one correspondence between the flow rate $Q$ and the steady state value $S_T^{eq}$, as long as $Q$ is chosen in $[0, \bar{Q}(S)]$ or equivalently $S_T^{eq}$ is chosen in the interval $[0, S]$. Consequently, we can consider that choosing the variable $Q$ in the interval $[0, \bar{Q}(S)]$ is equivalent to choosing the steady-state value $S_T^{eq}$ of the bioreactor in the interval $[0, S]$. This property will allow a simple characterization of an optimal control strategy as a state feedback given in Section 2.

The following Sections will consider extensions of this first model, with two different representations of the heterogeneity of the pollutant concentration in the water resource.

We first consider a concentration gradient between the input and output location points of the resource, defining two zones, as depicted on Figure 2a.

The reduced model that couples the time evolution of the pollutant concentrations in the bioreactor ($S_T$) and in the two zones ($S_1$ for the zone where water is extracted and $S_2$ for the one where water is injected) can be written as follows:

$$\dot{S}_1 = \frac{Q}{rV}(S_2 - S_1) \qquad (4a)$$

$$\dot{S}_2 = \frac{Q}{(1-r)V}(S_T^{eq}(Q) - S_2) \qquad (4b)$$
where $r \in (0, 1)$ represents the ratio of the volume if the first zone over the total volume $V$ of the water resource.

**Problem 2.** For a given threshold $\underline{S} > 0$ and initial concentrations $(S_1^0, S_2^0)$, we look for a time varying flow rate $Q(\cdot)$ that realizes the minimization

$$T_{\underline{S}}^{\text{opt}}(S_1^0, S_2^0) := \inf_{Q(\cdot) \in \mathcal{Q}(S_1^0, S_2^0)} T_{\underline{S}}(S_1^0, S_2^0, Q(\cdot))$$

where $T_{\underline{S}}(\cdot)$ is the time to reach the target:

$$T_{\underline{S}}(S_1^0, S_2^0, Q(\cdot)) = \inf \left\{ t \ge 0 \mid S_1^{(S_1^0, S_2^0, Q(\cdot))}(t) \le \underline{S} \text{ and } S_2^{(S_1^0, S_2^0, Q(\cdot))}(t) \le \underline{S} \right\}.$$

and $S_1^{(S_1^0, S_2^0, Q(\cdot))}(\cdot)$, $S_2^{(S_1^0, S_2^0, Q(\cdot))}(\cdot)$ denote the solution of (4) with initial condition $(S_1^0, S_2^0)$ and time varying control $Q(\cdot)$. Optimal controls $Q(\cdot)$ are sought in the set

$$Q(S_1^0, S_2^0) = \{ Q : \text{ measurable s.t. } Q(t) \in [0, \bar{Q}(S_2^{(S_1^0, S_2^0, Q(\cdot))}(t))], \forall t \in [0, T_{\underline{S}}(S_1^0, S_2^0, Q(\cdot))] \}.$$

We then consider instead a two patches representation: an “active” patch, where water is pumped out for purification and is re-injected back in the same patch after abatement; a “dead” patch that is supposed to be less directly influenced by advection but merely by passive diffusion with the active one, as depicted on Figure 2b.

This simple two patches model mimics the observation that one may roughly distinguish two zones that are advection or diffusion dominated. Typically, the active zone contains the shortest current lines that join input from output locations, while the dead one is expected to be away from those lines.

The reduced model that couples the time evolution of the pollutant concentrations in the bioreactor ($S_T$) and in the two zones ($S_A$ for the active one and $S_D$ for the dead one) can be written as follows:

$$\dot{S}_A = \frac{Q}{rV}(S_T^{eq}(Q) - S_A) + \frac{d}{rV}(S_D - S_A) \quad (5a)$$

$$\dot{S}_D = \frac{d}{(1-r)V}(S_A - S_D) \quad (5b)$$

where $r \in (0, 1)$ represents the ratio of the volume of the active patch over the total volume $V$ of the water resource, and $d$ is the diffusion rate of the pollutant between the two zones.

**Problem 3.** For a given threshold $\underline{S} > 0$ and initial concentrations $(S_A^0, S_D^0)$, we look for a time varying flow rate $Q(\cdot)$ that realizes the minimization

$$T_{\underline{S}}^{\text{opt}}(S_A^0, S_D^0) := \inf_{Q(\cdot) \in \mathcal{Q}(S_A^0, S_D^0)} T_{\underline{S}}(S_A^0, S_D^0, Q(\cdot))$$

where $T_{\underline{S}}(\cdot)$ is the time to reach the target:

$$T_{\underline{S}}(S_A^0, S_D^0, Q(\cdot)) = \inf \left\{ t \ge 0 \mid S_A^{(S_A^0, S_D^0, Q(\cdot))}(t) \le \underline{S} \text{ and } S_D^{(S_A^0, S_D^0, Q(\cdot))}(t) \le \underline{S} \right\}.$$

and $S_A^{(S_A^0, S_D^0, Q(\cdot))}(\cdot)$, $S_D^{(S_A^0, S_D^0, Q(\cdot))}(\cdot)$ denote the solution of (5) with initial condition $(S_A^0, S_D^0)$ and time varying control $Q(\cdot)$. Optimal controls $Q(\cdot)$ are sought in the set

$$Q(S_A^0, S_D^0) = \{ Q : \text{ measurable s.t. } Q(t) \in [0, \bar{Q}(S_A^{(S_A^0, S_D^0, Q(\cdot))}(t))], \forall t \in [0, T_{\underline{S}}(S_A^0, S_D^0, Q(\cdot))] \}.$$
## 2. OPTIMAL CONTROL UNDER HOMOGENEITY ASSUMPTION

In this part, we analyze Problem 1, assuming the water resource is perfectly homogeneous.

Since $\mu(\cdot)$ is increasing, choosing $Q$ is equivalent to choosing $S_T^{eq}$ as a control variable. Then (3) can be written as

$$ \dot{S} = \alpha\mu(S_T^{eq})(S_T^{eq} - S), \quad S_T^{eq} \in (0, S) \qquad (6) $$

where we denote $\alpha = V_T/V$ and $S_T^{eq}$ is now the control variable.

Thus we analyze two cases: constant and state feedback controls.

### 2.1. Best constant controls

Notice that for a constant control $S_T^{eq}$ (i.e., $Q = V_T\mu(S_T^{eq})$), solutions of (6) can be made explicit:

$$ S(t) = S_T^{eq} + (S(0) - S_T^{eq})e^{-\alpha\mu(S_T^{eq})t}. $$

Then, the time to reach $S = \underline{S}$ is given by the expression

$$ T_f(S_T^{eq}) = \frac{1}{\alpha\mu(S_T^{eq})} \ln \left( \frac{S(0) - S_T^{eq}}{\underline{S} - S_T^{eq}} \right). \qquad (7) $$

Therefore, we obtain a characterization for the best constant control.

**Proposition 2.1.** The best constant control $Q^*$ is defined as $Q^* = V_T\mu(S_T^*)$, where $S_T^*$ is the unique minimum on the interval $(0, \underline{S})$ of the function $T_f(\cdot)$ given by (7).

The proof of this proposition can be found in [4], where the existence and uniqueness of the best constant control is shown.

### 2.2. State feedback optimal control

Regarding now the best state feedback control, we prove in [4] the following result.

**Proposition 2.2.** The optimal feedback fulfills $Q^{opt}(S) = V_T\mu(S_T^{opt}(S))$ with

$$ S_T^{opt}(S) \in \underset{S_T^{eq} \in (0, S)}{\operatorname{argmax}} \mu(S_T^{eq})(S - S_T^{eq}). \qquad (8) $$

Moreover, $t \mapsto Q^{opt}(S(t))$ is decreasing along any optimal trajectory.

The principle of the above proposition and its proof, is to establish that the optimal feedback $S \to S_T^{opt}(S)$ is the one that makes the time derivative of $S$ the most negative at any time in (6).

Observe that a necessary condition for a feedback $S_T^{opt}(S)$ to be optimal is to verify

$$ \mu(S_T^{opt}(S)) = \mu'(S_T^{opt}(S))(S - S_T^{opt}(S)). \qquad (9) $$

The condition (9) gives explicit expressions of the optimal policy for some well known increasing growth functions:

$$ \begin{array}{c|c} \text{linear: } \mu(s) = \mu s, & \text{Monod: } \mu(s) = \frac{\mu_{\max}s}{K + s}, \\ \hline S_T^{opt}(S) = S/2 & S_T^{opt}(S) = \sqrt{K^2 + KS - K} \end{array} $$

<caption>TABLE 1. Optimal feedback control for linear and Monod growth functions.</caption>
<table><caption>(A) Values of the parameters for the numerical simulations.</caption><thead><tr><th>μ<sub>max</sub></th><th>K</th><th>V</th><th>V<sub>T</sub></th><th>S(0)</th></tr></thead><tbody><tr><td>1.0</td><td>1.0</td><td>1000</td><td>1</td><td>1</td></tr></tbody></table>

<table><caption>(B) Comparison with several targets S.</caption><thead><tr><th>S&#773;</th><th>T<sub>f</sub><sup>*</sup></th><th>T<sup>opt</sup></th><th>gain</th></tr></thead><tbody><tr><td>0.5</td><td>5420</td><td>5293</td><td>1.9%</td></tr><tr><td>0.4</td><td>8064</td><td>7719</td><td>3.7%</td></tr><tr><td>0.3</td><td>12499</td><td>11606</td><td>7.2%</td></tr><tr><td>0.2</td><td>21663</td><td>19062</td><td>12%</td></tr><tr><td>0.1</td><td>51626</td><td>40425</td><td>22%</td></tr></tbody></table>

<caption>TABLE 2. Comparison of constant and feedback strategies</caption>

We show the benefit of considering non-constant flow rate $Q$ for optimizing the treatment time. Consider the Monod law (See Table 1) and the values of the problem's parameters given in Table 2a. The numerical comparisons of the time for optimal constant control $T_f^*$ to reach the target and the minimal time $T^{opt}$ are reported on Table 2b. We can see that the gain can be substantial when dealing with small set-points $S$.

### 3. OPTIMAL CONTROL FOR SIMPLE SPATIAL PATTERNS

In this section we study problems 2 and 3 that break the homogeneity assumption in the water resource $V$. In the first subsection we simply consider a gradient of concentration from the input to the output zone (Problem 2) while, in the second subsection, we consider a dead zone which cannot be reached by the depolluting unit (Problem 3). So, the depollution in that zone is only made indirectly via its connection with the treated zone.

#### 3.1. A simple gradient between input and output

Denoting $\alpha_i = V_T/V_i$ ($i = 1, 2$), (4) results in the following dynamics:

$$ 
\begin{aligned}
\dot{S}_1 &= \alpha_1 \mu(S_T^{eq})(S_2 - S_1) \\
\dot{S}_2 &= \alpha_2 \mu(S_T^{eq})(S_T^{eq} - S_2).
\end{aligned}
\qquad (10) $$

We can easily check that the domain $\mathcal{D} = \{(S_1, S_2) \in \mathbb{R}_+^2 | S_1 \ge S_2\}$ is invariant for any control $S_T^{eq}(\cdot)$ such that $S_T^{eq}(t) \in (0, S_2(t)]$ for any $t > 0$. We shall consider initial conditions in $\mathcal{D}$, and we define the target $\mathcal{T} = \{(S_1, S_2) \in \mathcal{D} | S_1 \le \underline{S}_t\}$. Also, for $p \in [0, 1]$ and $\tau \ge 0$, we define the following function:

$$ A(p, \tau) = \begin{cases} e^{-\alpha \tau} & \text{if } p = 0 \\ \frac{(1-p)e^{-\frac{\alpha \tau}{1-p}} - p e^{-\frac{\alpha \tau}{p}}}{1-2p} & \text{if } p \in (0, \frac{1}{2}) \\ \left(1+\frac{\alpha}{2}\tau\right) e^{-\frac{\alpha \tau}{2}} & \text{if } p = \frac{1}{2} \\ A(1-p, \tau) & \text{if } p \in (\frac{1}{2}, 1]. \end{cases} $$

Finally, for $S_0 > \underline{S}$ we introduce mapping $B(S_T^{eq}) := \frac{S-S_T^{eq}}{S_0-S_T^{eq}}$ for $S_T^{eq} \in (0, \underline{S})$.

**Proposition 3.1.** Let $p = V_1/V$. For initial conditions such that $S_1(0) = S_2(0) = S_0 > \underline{S}$, the best constant control $Q^*$ and time $T_f^*$ to reach the target $\mathcal{T}$ are defined by $Q^* = V_T\mu(S_T^*)$, where $S_T^*$ is such that the graph of $B(\cdot)$ tangentially intersects the graph of $S_T^{eq} \mapsto A(p, \mu(S_T^{eq})T_f^*)$ at $S_T^{eq} = S_T^*$.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Graphical determination of $S_T^*$ and $T_f^*$.</figcaption></figure>

<table><caption>TABLE 3. Comparison of constant and feedback strategies</caption><thead><tr><th rowspan="2">S</th><th colspan="1">p = 0</th><th colspan="3">p = 0.2</th><th colspan="3">p = 0.5</th></tr><tr><th>T<sub>f</sub><sup>*</sup></th><th>T<sub>f</sub><sup>*</sup></th><th>T<sub>opt</sub></th><th>%</th><th>T<sub>f</sub><sup>*</sup></th><th>T<sub>opt</sub></th><th>%</th></tr></thead><tbody><tr><td>0.7</td><td>2387</td><td>2708</td><td>2661</td><td><b>1.7</b></td><td>2931</td><td>2881</td><td><b>1.7</b></td></tr><tr><td>0.6</td><td>3660</td><td>3816</td><td>3756</td><td><b>1.6</b></td><td>3993</td><td>3951</td><td><b>1.1</b></td></tr><tr><td>0.5</td><td>5420</td><td>5325</td><td>5194</td><td><b>2.5</b></td><td>5385</td><td>5317</td><td><b>1.2</b></td></tr><tr><td>0.4</td><td>8064</td><td>7567</td><td>7329</td><td><b>3.2</b></td><td>7395</td><td>7250</td><td><b>2.0</b></td></tr><tr><td>0.3</td><td>12499</td><td>11338</td><td>10754</td><td><b>5.4</b></td><td>10696</td><td>10320</td><td><b>3.6</b></td></tr><tr><td>0.2</td><td>21663</td><td>19137</td><td>17297</td><td><b>11</b></td><td>17389</td><td>16 180</td><td><b>7.4</b></td></tr></tbody></table>

For the sake of notation, for $S_1^0 > \underline{S} > S_2^0$, we define when $\alpha_1 \neq \alpha_2$:

$$f_0(S_1^0, S_2^0) = \frac{\alpha_2 \left( 1 - \left( \frac{S_1^0 - S_2^0}{\bar{S}_1 - S_2^0} \right) \frac{\alpha_1 - \alpha_2}{\alpha_1} \right)}{(\alpha_2 - \alpha_1)(S_1^0 - S_2^0)}, \quad \beta = \left( \frac{\alpha_1}{\alpha_2} \right)^{\frac{\alpha_1}{\alpha_1 - \alpha_2}}$$

and when $\alpha_1 = \alpha_2$: $f_0(S_1^0, S_2^0) = \frac{\ln\left(\frac{S_1^0 - S_2^0}{\bar{S}_1 - S_2^0}\right)}{(S_1^0 - S_2^0)}$, $\beta = e$.

The main result of this section, stated here below, concerns the optimality of feedback-type strategies. It can be found in [4].

**Proposition 3.2.** Under Assumption A1, the constant control $Q^{opt} = V_T\mu(\underline{S}_2)$ is optimal in the set

$$\mathcal{I} := \{(S_1^0, S_2^0) \in (\underline{S}, +\infty) \times (0, \underline{S}) : S_2^0 \le \bar{S}_2 \text{ or } \mu(S_2^0) f_0(S_1^0, S_2^0) \le \mu'(S_2^0)\}$$

where $\underline{S}_2$ is the value of $S_2$ when $\mathcal{I}$ is reached and $\bar{S}_2$ is the unique solution in $(0, \underline{S})$ of

$$\mu(\bar{S}_2) = \beta\mu'(\bar{S}_2)(\underline{S} - \bar{S}_2). \qquad (11)$$

Moreover, for any initial condition in $D \setminus \mathcal{T}$, the optimal control $Q^{opt}(\cdot)$ consists in reaching $\mathcal{I}$ is optimal until $S_1(\cdot)$ reaches $\underline{S}$. Moreover, $t \mapsto Q^{opt}(S_1(t))$ increases as it approaches the set $\mathcal{I}$.

## 3.2. Consideration of dead-zone

As before, considering $S_T^{eq}$ as the control in equations (5) we obtain the following system:
$$ \begin{aligned} \dot{S}_A &= \alpha\mu(S_T^{eq})(S_T^{eq} - S_A) + \beta_A(S_D - S_A) \\ \dot{S}_D &= \beta_D(S_A - S_D), \end{aligned} \qquad (12) $$

where $\alpha = \frac{V_T}{V_A}$, $\beta_A = \frac{d}{V_A}$ and $\beta_D = \frac{d}{V_D}$.

In this section we characterize the optimal strategy which solves our lake depollution problem under the situation described above. We thus obtain the somehow surprising result that this optimal strategy coincides with the optimal one obtained in the homogeneous case (see (9)). Indeed, by means of Pontryagin's Maximum Principle it was proven in [5] the following result.

**Proposition 3.3.** The optimal flow rate $Q^{opt}$ is given by $Q^{opt}(S_A) = V_T\mu(S_T^{opt}(S_A))$, where $S_T^{opt}$ is the feedback control characterized by

$$ S_T^{opt}(S_A) = \underset{S_T^{eq} \in (0, S_A)}{\operatorname{argmax}} \mu(S_T^{eq})(S_A - S_T^{eq}), \qquad (13) $$

Notice that, in particular, the optimal feedback policy $S_T^{opt}$ does not depend on the substrate concentration $S_D$ in the dead zone. Also, it belongs to the interval $[0, S_A(t)]$, for all $t \in [0, T]$.

# 4. ACCURATE MODELING AND SIMULATION OF THE HYDRODYNAMICS

In order to compare and validate the simplified models of Sections 2 and 3 above, we now consider a model with space dependency in the resource thanks to partial derivative equations (PDE). In the sequel we assume that $S(t, \mathbf{x})$ is a function of time and space<fnref n="1" />.

## 4.1. The system of PDE

The evolution of the pollutant concentration $S(t, \mathbf{x})$ in the resource is modeled by the following advection-diffusion equation:

$$ \frac{\partial S}{\partial t} + \mathbf{U} \cdot \nabla S - \nu_S \Delta S = 0, \qquad (14) $$

where $\nu_S$ is the pollutant diffusivity, and $\mathbf{U}(t, \mathbf{x}) = (u(t, \mathbf{x}), v(t, \mathbf{x}))$ is the fluid velocity in the resource, computed thanks to the Navier-Stokes equations (the reference model for incompressible viscous fluids, see [10]):

$$ \frac{\partial \mathbf{U}}{\partial t} + \mathbf{U} \cdot \nabla \mathbf{U} + \nabla p - \nu_u \Delta \mathbf{U} = 0, \qquad (15a) $$

$$ \nabla \cdot \mathbf{U} = 0. \qquad (15b) $$

Here, $p$ stands for the pressure exerted on the fluid, and $\nu_u$ denotes the water viscosity. Equations (14) and (15) are supplemented with appropriate initial conditions (at $t = 0$) and boundary conditions on $\partial\Omega$ (see [2] for details). In particular, for the velocity field the conditions on the input part of the domain boundary read<fnref n="2" />:

$$ \mathbf{U}(t, \mathbf{x}) = Q(t) \mathbf{U}_0(\mathbf{x}), \quad \forall \mathbf{x} \in \partial\Omega, $$

where $\mathbf{U}_0$ is a given function with a Poiseuille profile. These boundary conditions explicitly depend on the discharge $Q(t)$, inducing a coupling between equations (14) and (15) through the equations (1a) and (1b) that are kept in the bioreactor.

<footnote n="1">¹For the sake of simplicity, we only consider 2D (horizontal) models so that $\mathbf{x} \in \Omega \subset \mathbb{R}^2$.</footnote>

<footnote n="2">²The same kind of boundary condition holds for output boundary conditions.</footnote>
<figure><img src="image_5.png" /><figcaption>(A) Evolution of the (renormalized) pollutant concentration in the resource thanks to the coupled PDE-ODE system. The water is extracted from the left side of the domain, treated in the bioreactor and injected back through the right side. Snapshot of the pollutant concentration at $t_1 > 0$.</figcaption></figure>

<caption>(B) Time-evolution of the pollutant concentration in the resource with uniform initial value $S_0 = 15$. Reference PDE model (upper red), single zone model (lower black) and two zones model (upper blue). Parameters: $\nu_S = 5 \cdot 10^{-3}$, $(r, d) = (0.346, 0.228)$. The blue and red curves are almost indistinguishable</caption>

FIGURE 4

The numerical results provided by the coupled Navier-Stokes/Advection-Diffusion system (see Figure 4a) provide two main informations: i) it is not correct to assume that the pollutant concentration is homogeneous in the resource and ii) as a first approximation, the two zones model (dead for red and active for blue) can be considered.

## 4.2. Calibration of the O.D.E. models

In order to calibrate the two-zones model described in Section 3 above, one can optimize the values of the volume ratio $r = V_A/V$ and the diffusivity $d > 0$ and thus better approximate the reference numerical solution obtained with the complete Navier-Stokes equations. Figure 4b illustrates the fact that the two-zones model (with appropriate values of $r$ and $d$) is a much more accurate approximation of the reference model than the single zone one (corresponding to $r = 1, d = 0$). Thanks to previous works on optimal control, a minimal time strategy can be derived from the two-zones model that can then be applied to the reference model. As shown formerly, the optimal strategy that provides the online adjustment of the flow rate for minimizing the total decontamination time requires the single on-line measurement of the pollutant concentration at the output of the lake, that is provided by the reference model. Finally, this two zones approximation provides

* a simple and easily applicable feedback control of the flow rate,

* simple approximations of the input-output dynamics of the pollutant concentration in the resource and the expected minimal time of decontamination (with a system of two ODEs that offers quite light computation times), thanks to the combination of numerical simulations of a PDE based model of the hydrodynamics in the lake, and a calibrated system of ODEs for given geometrical data of the lake and diffusivity properties of the pollutant.

Even better approximations with higher dimensions of the ODE models (i.e. considering more than two zones) could be obtained with the same methodology. This approach with feedback control on simple models provide bases for future efficient decision making tools for the bioremediation of natural water reservoirs.
## 5. MULTI-PUMPS CONSIDERATION

In this Section we consider the scheme depicted on Figure 2b but considering that the output from the treatment process can be split directly in the two patches via two different pumps, whose flows are denoted $Q_1$ and $Q_2$. To preserve the volumes of each patch constant, we assume that the same flows $Q_1$ and $Q_2$ are pumped out, to be mixed and treated by the bioreactor as an unique input flow $Q = Q_1 + Q_2$ (see Figure 5). Therefore the two

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Consideration of two pumps dedicated to two patches in the resource.</figcaption></figure>

zones could be “active” simultaneously. Notice that in this framework the particular case $d = 0$ is relevant as it corresponds to the treatment of two independent resources with an unique treatment unit.

### 5.1. The optimal control problem

The reduced model that couples the time evolution of the pollutant concentrations in the bioreactor ($S_T$) and in the two patches (denoted $S_1$ and $S_2$) can be written as follows:

$$ \dot{S}_1 = \frac{Q_1}{rV}(S_T^{eq}(Q) - S_1) + \frac{d}{rV}(S_2 - S_1) \quad (16a) $$

$$ \dot{S}_2 = \frac{Q_2}{(1-r)V}(S_T^{eq}(Q) - S_2) + \frac{d}{(1-r)V}(S_1 - S_2) \quad (16b) $$

where $r \in (0, 1)$ represents the ratio of the volume of patch 1 over the total volume $V$ of the water resource and $d$ is the diffusion rate of the pollutant between the two zones. Notice that here the pollutant concentration that enters the bioreactor is $(Q_1S_1 + Q_2S_2)/Q$ and that one has now two controls: $Q_1$ and $Q_2$.

**Problem 4.** For a given threshold $\underline{S} > 0$ and initial concentrations $(S_1^0, S_2^0)$, we look for time varying flow rates $Q_1(\cdot)$ and $Q_2(\cdot)$ that realize the minimization

$$ T_{\underline{S}}^{opt}(S_1^0, S_2^0) := \inf_{(Q_1(\cdot), Q_2(\cdot)) \in \mathcal{Q}(S_1^0, S_2^0)} T_{\underline{S}}(S_1^0, S_2^0, Q_1(\cdot), Q_2(\cdot)) $$

where $T_{\underline{S}}(\cdot)$ is the time to reach the target:

$$ T_{\underline{S}}(S_1^0, S_2^0, Q_1(\cdot), Q_2(\cdot)) = \inf \left\{ t \ge 0 \mid S_1^{(S_1^0, S_2^0, Q_1(\cdot), Q_2(\cdot))}(t) \le \underline{S} \text{ and } S_2^{(S_1^0, S_2^0, Q_1(\cdot), Q_2(\cdot))}(t) \le \underline{S} \right\}. $$
and $S_i^{(S_1^0, S_2^0, Q_1(\cdot), Q_2(\cdot))}(\cdot)$ ($i = 1, 2$) denote the solution of (16) with initial condition $(S_1^0, S_2^0)$ and time varying controls $Q_1(\cdot), Q_2(\cdot)$. Optimal controls $Q_1(\cdot), Q_2(\cdot)$ are sought in the set

$$Q(S_1^0, S_2^0) = \left\{ \begin{array}{l} (Q_1, Q_2) \in \mathbb{R}_+^2 : \text{measurable s.t. } \forall t \in [0, T_S(S_1^0, S_2^0, Q_1(\cdot), Q_2(\cdot))], \\ Q_1(t) + Q_2(t) \le \bar{Q} \left( \sum_{i=1}^2 Q_i(t) S_i^{(S_1^0, S_2^0, Q_1(\cdot), Q_2(\cdot))}(t) / \sum_i Q_i(t) \right) \end{array} \right\}$$

As previously, it is convenient to replace the controls $(Q_1, Q_2)$ by the equivalent controls $(\alpha, S_T^{eq})$ such that $Q_1 = \alpha Q$, $Q = V_T \mu(S_T^{eq})$. Defining the function $\beta(\sigma, S_T^{eq}) = \mu(S_T^{eq})(\sigma - S_T^{eq})$ and changing the time $t$ by $tV/V_T$, the controlled dynamics can be re-written as follows

$$\frac{d}{dt} \begin{bmatrix} S_1 \\ S_2 \end{bmatrix} = \underbrace{\begin{bmatrix} -\frac{\alpha}{r}\beta(S_1, S_T^{eq}) \\ -\frac{1-\alpha}{1-r}\beta(S_2, S_T^{eq}) \end{bmatrix}}_{F(S,(\alpha,S_T^{eq}))} + \underbrace{\begin{bmatrix} \frac{S_2-S_1}{r} \\ \frac{S_1-S_2}{1-r} \end{bmatrix}}_{G(S)} \quad (17)$$

with $\bar{d} = d/V_T$ and where the control set is $U(S) = \{(\alpha, S_T^{eq}) | \alpha \in [0,1], S_T^{eq} \in [0, \alpha S_1 + (1-\alpha)S_2]\}$. Unfortunately, the velocity set $\cup_{(\alpha,S_T^{eq}) \in U(S)} F(S, (\alpha, S_T^{eq}))$ is non convex when $S_1 \neq S_2$, and consequently one cannot guarantee the existence of an optimal control. Instead we consider relaxed controls: $(\alpha^a, S_T^a, \alpha^b, S_T^b, p) \in \tilde{U}(S) = U(S)^2 \times [0,1]$ and the relaxed dynamics $\frac{dS}{dt} = \tilde{F}(S, \tilde{u}) + \bar{d}G(S)$ with $\tilde{F}(S, \tilde{u}) = pF(S, (\alpha^a, S_T^a)) + (1-p)F(S, (\alpha^b, S_T^b))$. In general, an optimal solution of the relaxed problems can be approached by sequences of “chattering” controls (but which are hardly practicable from an application viewpoint). Nevertheless, for this particular problem, we have been able to prove that there exists a state feedback with regular controls that is optimal for the relaxed problem. We denote the function $\gamma(\sigma) = \max_{S_T \in [0,\sigma]} \beta(\sigma, S_T)$ whose maximum is reached at an unique $\hat{S}_T(\sigma)$.

**Proposition 5.1.** [6] The feedback

$$(\alpha^{\text{opt}}, S_T^{\text{opt}}) = \begin{cases} (1, \hat{S}_T(S_1)) & \text{when } S_1 > S_2, \\ (r, \hat{S}_T(S_1 = S_2)) & \text{when } S_1 = S_2, \\ (0, \hat{S}_T(S_2)) & \text{when } S_1 < S_2. \end{cases}$$

is optimal.

The optimal strategy consists then in reaching $S_1 = S_2$ using only one input pump, and then to stay on the singular arc $S_1 = S_2$ using the two input pumps whatever the parameter $d$ is. This last property means that the knowledge of the diffusion parameter $d$ is not necessary to compute on-line the optimal flows, but the value of $d$ impacts the time to reach the target set $T := [0, \bar{S}]^2$. Consider the value function parameterized by $d$:

$$T_d^{\text{opt}}(S_0) := \inf_{u(\cdot)} \{t > 0 | S_{S_0, u, d}(t) \in T\}$$

where $S_{S_0, u, d}$ stands for the solution of (17) for initial vector $S_0$, control $u(\cdot) = (\alpha(\cdot), S_T^{eq}(\cdot))$ and parameter $d$. Define the functions

$$\bar{T}(s) := \int_{S}^{s} \frac{d\sigma}{\gamma(\sigma)}, \quad T(s) := \max(0, \bar{T}(s))$$

and the subset $\Delta := \{S \in \mathbb{R}_+^2 \text{ s.t. } S_1 = S_2\}$.

**Proposition 5.2.** [6] One has the following properties.
(1) $\mu(\cdot)$ concave implies that $T(\cdot)$ is concave on $[\underline{S}, +\infty)$

(2) $T_0^{opt}(S) = rT(S_1) + (1-r)T(S_2)$ for any $S \in \mathbb{R}_+^2$

(3) $T_0^{opt}(S) = r\bar{T}(S_1) + (1-r)\bar{T}(S_2)$ for $S \notin T$ with $S_i \neq \underline{S}$

(4) $T_0^{opt}(S) = \max(r\bar{T}(S_1), (1-r)\bar{T}(S_2), r\bar{T}(S_1) + (1-r)\bar{T}(S_2))$ for $S \notin T$ with $S_1 = \underline{S}$ or $S_2 = \bar{S}$

(5) $T_d^{opt}(S) = T(S_1) = T(S_2)$ for any $S \in \Delta$ and $d \ge 0$

(6) $d \mapsto T_d^{opt}(S)$ is increasing for any $S \in (\underline{S}, +\infty)^2 \setminus \Delta$

(7) $T_{\infty}^{opt}(S) = \lim_{d \to +\infty} T_d^{opt}(S) = T(rS_1 + (1-r)S_2)$ for any $S \in \mathbb{R}_+^2$

Therefore, one has $T_0^{opt}(S) \le T_d^{opt}(S) < T_{\infty}^{opt}(S)$ for any $S \in (\underline{S}, +\infty)^2 \setminus \Delta$, but this might be not true when $S_i < \underline{S}$, as one can see on Figure 6.

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Level sets of the minimal time function for $\mu(S) = S/(1+S)$ and $r = 0.3$ ($d = 0$ on the left, $d = +\infty$ on the right)</figcaption></figure>

## 5.2. Simulation of the feedback strategy on the P.D.E. model

A system of PDEs similar to Section 4 has been implemented with several pumps. The only difference is that the input boundary is now separated in two or more parts, where the boundary conditions for the velocity read:

$$U(t, \mathbf{x}) = Q_n(t) U_0(\mathbf{x}), \quad \forall \mathbf{x} \in \partial\Omega_n,$$

where $\cup_{1 \le n \le N} \partial\Omega_n$ corresponds to the whole input boundary (with $N$ pumps). In Figure 7, we compare three strategies: one pump, 3 pumps with uniform distribution of the pumping power $Q_n(t) = Q(t)/3$, and 3 pumps with $Q(t)$ chosen with a bang-bang strategy as in [7].

## 6. CONCLUSION AND PERSPECTIVES

The main message of this short survey is that in-situ decontamination takes clearly advantage of real-time strategies that adjust dynamically the flow rate from the on-line measurement of the pollution at the pumping location. It is quite remarkable that simple o.d.e. models with little knowledge and sparse representations of the water hydrodynamics could lead to efficient treatment strategies, that are moreover very simple to be implemented. The key aspect is to calibrate these models on realistic numerical simulations of the hydrodynamics of the resource, based on the Navier-Stokes equations and their relevant boundary conditions. This approach seems rather new in the field and has not been yet tested in real situations.
<figure><img src="image_8.png" /><figcaption>FIGURE 7. One pump, 3 pumps with $Q_n(t) = Q(t)/3$, and 3 pumps with $Q(t)$ chosen as in [7].</figcaption></figure>

The criterion under consideration has been the single time to reach a prescribed low level of pollutant, while the pumping energy is most of the time non negligible from the economical view point. Therefore, extensions of current works would consider more general criterions to study if some compromise between the treatment duration and the spent energy could lead to significant different strategies.

The consideration of more than one pump in the optimal control problem leads to some mathematical technicalities (due to lack of convexity) for which the extension to more than two pumps is currently an open problem. From a practical viewpoint, the consideration of several points could be accompanied to the study of the best location of the pumps and even the possibility to change (optimally) these locations with time. Such a perspective should lead in the future to the development of computer decision making tools, coupling numerical simulations of the hydrodynamics and optimal control tools, in a more sophisticated way but in the same spirit than the one presented here.

## ACKNOWLEDGEMENT

The authors are grateful for the support of DYMECOS INRIA Associated Teams, project BIONATURE of CIRIC INRIA CHILE, and CONICYT grants REDES 130067 and 150011. The first, third and fifth authors were also supported by FONDECYT regular projects 1160204 and 1160567, and BASAL project (Centro de Modelamiento Matemático, Universidad de Chile).

The authors are also grateful to T. Bayen for fruitful discussions and insightful ideas.

## REFERENCES

[1] E.S. Bakker, E. Van Donk, and A.K. Immers. Lake restoration by in-lake iron addition: a synopsis of iron impact on aquatic organisms and shallow lake ecosystems. *Aquatic Ecology*, 50:121-135, 2016.

[2] S. Barbier, A. Rapaport, and A. Rousseau. Modelling of biological decontamination of a water resource in natural environment and related feedback strategies. *Journal of Scientific Computing*, 68(3):1267-1280, 2016.

[3] D. Copetti, K. Finsterle, L. Marziali, and et al. Eutrophication management in surface waters using lanthanum modified bentonite: A review. *Water Research*, 97:162-174, 2016.

[4] P. Gajardo, J. Harmand, H. Ramírez C., and A. Rapaport. Minimal time bioremediation of natural water resources. *Automatica J. IFAC*, 47(8):1764-1769, 2011.

[5] P. Gajardo, H. Ramírez, V. Riquelme, and A. Rapaport. Bioremediation of natural water resources via optimal control techniques. In *BIOMAT 2011*, pages 178-190. World Sci. Publ., Hackensack, NJ, 2012.

[6] H. Ramírez, A. Rapaport, and V. Riquelme. Optimal Feedback Synthesis and Minimal Time Function for the Bioremediation of Water Resources with Two Patches. *SIAM J. Control Optim.*, 54(3):1697-1718, 2016.
[7] A. Rapaport, A. Rousseau, and J. Harmand. Procédé de traitement d'une ressource fluide, programme d'ordinateur et module de traitement associés, 2014. FA 78 4546 - FR 13 55129.

[8] E. Rurangwa and M.C.J. Verdegem. Microorganisms in recirculating aquaculture systems and their management. *Reviews in aquaculture*, 7(2):117-130, 2015.

[9] H. . Smith and P. Waltman. *The theory of the chemostat*, volume 13 of *Cambridge Studies in Mathematical Biology*. Cambridge University Press, Cambridge, 1995. Dynamics of microbial competition.

[10] R. Temam. *Navier-Stokes equations*. AMS Chelsea Publishing, Providence, RI, 2001. Theory and numerical analysis, Reprint of the 1984 edition.

[11] H. Weiping. A review of the models for lake taihu and their application in lake environmental management. *Ecological Modelling*, 319:9-20, 2016.