WATER ARTIFICIAL CIRCULATION FOR EUTROPHICATION CONTROL

. This work analyzes, from a mathematical point of view, the artiﬁcial mixing of water - by means of several pairs collector/injector that set up a circulation pattern in the waterbody - in order to prevent the undesired eﬀects of eutrophication. The environmental problem is formulated as a constrained optimal control problem of partial diﬀerential equations, where the state system is related to the velocity of water and to the concentrations of the diﬀerent species involved in the eutrophication processes, and the cost function to be minimized represents the volume of recirculated water. In the main part of the work, the wellposedness of the problem and the existence of an optimal control is demonstrated. Finally, a complete numerical algorithm for its computation is presented, and some numerical results for a realistic problem are also given.

1. Introduction. Artificial circulation and aeration are management techniques for oxygenating eutrophic waterbodies subject to quality problems, such as loss of oxygen, sediment accumulation and algal blooms. These techniques are typically applied in stratified shallow waters to mitigate problems associated with oxygen depletion. Artificial circulation disrupts stratification and minimizes the development of stagnant zones that may be subject to water quality problems. The movement of water and/or air is accomplished by use of bottom diffusers and water pumps to create a circulation pattern in shallow waters [11].
The general purpose of artificial mixing and oxygenation is to increase the dissolved oxygen content of the water. Various systems are available to help do this by mechanically mixing or agitating the water. Artificial circulation and aeration can increase fish and other aquatic animal habitat, prevent fishkills, reduce algal blooms, and improve the quality of domestic and industrial water supplies, decreasing treatment costs.
This work deals with artificial circulation as a shallow water aeration technique. Large waterbodies (for instance, lakes or reservoirs) get much of their oxygen from the atmosphere through diffusion processes. Artificial circulation increases water's oxygen by forcefully circulating the water to expose more of it to the atmosphere. Proper choice and design of an artificial circulation system depends on the management goals and the geophysical characteristics. Two techniques are the most common: air injection and mechanical mixing. The former has been analyzed, from an optimal control viewpoint in a few works (see, for instance, [1] and the references therein). However, in this paper we will focus our attention on the latter that, as far as we know, has remained unaddressed in the mathematical literature.
In air injection systems, a compressor delivers air through lines connected to perforated pipes or diffusers placed near the bottom. The rising air bubbles cause water in the bottom layer to also rise, pulling this water into the surface. This aeration technique is sometimes referred to as the air-lift method of circulation, since bottom waters are lifted to the lake surface through the action of the injected air [13]. On the contrary, mechanical mixing systems use a top-down approach to set up a circulation pattern. A flow pump takes water from the epilimnion (well aerated upper layer of water) by means of a collector, injecting it into the hypolimnion (poorly oxygenated bottom layer), setting up a circulation pattern that prevents stratification. In this way, oxygen-poor water from the bottom is circulated to the surface, where oxygenation from the atmosphere can occur. Essentially, this technique can be understood as an artificial promoter of water movement, recirculating bottom water to the surface [17].
Most of the effects of destratification and circulation are beneficial for the environment, and include two that turn out to be fundamental for our purposes of eutrophication remediation: (a) an increase in dissolved oxygen through the water column, and (b) a reduction of algal growth by reducing its exposure to light, by changing water chemistry (pH, CO 2 and so on), or by mixing of algae-grazing zooplankton.
So, in next section we will introduce a mathematical formulation of the environmental problem as a control/state constrained optimal control problem of partial differential equations, analyzing the corresponding state system, in order to set its wellposedness. Then, we will prove the existence of solution for the optimal control problem. Finally, we will deal with the numerical resolution of the problem, presenting a complete numerical algorithm and a realistic computational example.
2. Mathematical formulation of the control problem.
2.1. The physical domain. We consider a domain Ω ⊂ R 3 corresponding, for instance, to a lake. We suppose the existence of a set of N CT pairs collectorinjector, assuming that each water collector is attached to a unique injector by a pipe with a pumping/turbine group.
We assume a smooth enough boundary ∂Ω, such that it can be split into four disjoint subsets ∂Ω = Γ S ∪ Γ B ∪ Γ C ∪ Γ T , where: • Γ S corresponds to the surface of the lake, • Γ B corresponds to the bottom of the lake, • Γ C correponds to the part of the boundary where the water collectors C k , k = 1, . . . , N CT , can be located, • Γ T correponds to the part of the boundary where the water injectors T k , k = 1, . . . , N CT , can be located. In order to simulate the dynamics of water we consider the modified Navier-Stokes equations following the Smagorinsky model of turbulence: where v(x, t) is the velocity of water, F(x, t) if the source term, n represents the unit outward normal vector to the boundary ∂Ω, µ(S) denotes the usual measure of a generic set S, the water velocity on the surface v S satifies that v S · n = 0, v 0 is the initial velocity, and, for each k = 1, . . . , N CT , g k (t) ∈ H s (0, T ), with s ≥ 1 a suitable index, represents the volume of water displaced by pump k at each time t. This volume of water can be positive g k (t) > 0 (we say that we are turbinating: water enters by the collector and is turbinate by the pipeline to the injector), or negative g k (t) < 0 (in this case, we say that we are pumping: water enters by the injector and is pumped to the collector). Finally, the term Ξ(v) is given by: where D is a potential function (for instance, in the particular case of the classical Navier-Stokes equations, D( ) = ν [ : ] and, consequently, Ξ(v) = 2ν (v)). In our case, the Smagorinsky turbulence model, the potential function is [14]: so, where β( (v)) = 2ν + 2ν tur [ (v) : (v)] 1/2 . From a mathematical point of view, the advantage of considering the modified Navier-Stokes equations, besides being more appropriate for turbulent flows, lies in the fact that, if the potential function fulfills certain properties (see, for instance, [10,9] and the references therein), it is possible to demonstrate the uniqueness of solution, in addition to gain in regularity with respect to that obtained for the Navier-Stokes equations. The potential function (3) for the Smagorinsky turbulence model considered in this work, as we will see when analyzing the existence of solution for the state equations, fulfills the necessary conditions to guarantee the uniqueness and regularity of the solution.
2.2.2. The eutrophication model: Michaelis-Menten kinetics. We consider the following species concentrations: • u 1 for a generic nutrient (usually, nitrogen and/or phosphorus), • u 2 for phytoplankton • u 3 for zooplankton, • u 4 for organic detritus, and • u 5 for dissolved oxygen. Interactions between five above species can be modelled by the following system of coupled partial differential equations [4]: for i = 1, . . . , 5, where sign(y) denotes the sign function: is defined by the following nonlinear expression: where: • C oc is the oxygen-carbon stoichiometric relation, • C nc is the nitrogen-carbon stoichiometric relation, • C f z is the zooplankton grazing efficiency factor, • K rd is the detritus regeneration rate, • K r is the phytoplankton endogenous respiration rate, • K mf is the phytoplankton death rate, • K mz is the zooplankton death rate (including predation), • K z is the zooplankton predation (grazing), • K F is the phytoplankton half-saturation constant, • K N is the nitrogen half-saturation constant, • Θ is the detritus regeneration thermic constant, • µ i , i = 1, . . . , 5, are the diffusion coefficients of each species, • f i (x, t), i = 1, . . . , 5, are the source terms of each species, • θ(x, t) is the water temperature, • L is the luminosity function, given by: with I 0 the incident light intensity, I s the light saturation, C t the phytoplankton growth thermic constant, ϕ 1 the light attenuation due to depth, ϕ 2 the light attenuation due to phytoplankton, and µ the maximum phytoplankton growth rate. All above coefficients are supposed to be nonnegative, except for the half-saturation constants that are strictly positive, Moreover, for the sake of simplicity in the theoretical part, we assume without loss of generality that the attenuation constant ϕ 2 = 0. Remark 1. It is worthwhile noting here that the boundary conditions of problem (5) associated to collectors and injectors change depending on the pump operation mode (turbination, pumping or rest): • Turbination mode (g k (t) > 0): in this case, sign(g k (t)) = 1, thus the corresponding boundary conditions read, for each i = 1, . . . , 5: That is, concentration of species i in the injector is the mean concentration in the collector, where water enters the pipeline to be turbinated. • Pumping mode (g k (t) < 0): in this case, sign(g k (t)) = −1, thus: That is, concentration of species i in the collector is the mean concentration in the injector, where water enters the pipeline to be pumped. • Rest mode (g k (t) = 0): in this case, sign(g k (t)) = 0, thus: That is, the system is isolated. 3. The optimal control problem. Our main aim is related to mitigating the adverse effects of eutrophication. In order to do this, many strategies are possible, for example, optimizing the locations and/or sizes of the collectors and injectors, minimizing the pumped flow, etc. To fix ideas, we will focus in this last problem: the optimization of the volume of pumped water in order, for instance, to minimize the energy cost. The effects of eutrophication can be measured, for example, on the basis of the dissolved oxygen level in the deeper areas (which is usually low due to the decomposition of organic detritus). The idea is to get some circulation of water from the upper layers (rich in oxygen due to photosynthesis realized by phytoplankton and to the contact with the atmosphere) to lower areas (poor in oxygen). In addition, with this artificial circulation, we intend to break with the stratification of the species, which results in a better quality of water.
So, to address the mitigation of the harmful effects of eutrophication by controlling the flow pumped by the injectors/collectors, we assume that their geometry and position are fixed beforehand, and that we can only act on the pumped flow rate. Thus, the following optimal control problem needs to be solved: with U a suitable functional space (that will be detailed in next section) and the positive constants c 1 , c 2 related to technological constraints on the pumps, where J(g) is the cost functional given by: representing the energy saving related to pumps' operation, and where K C ⊂ R 5 is the rectangle: with λ m k and λ M k the minimal and maximal mean concentrations allowed for the species k, for k = 1, . . . , 5, in the control domain Ω C ⊂ Ω. Finally, (v, u) are the solutions of the state system (1) and (5) associated to control g.

3.
Mathematical analysis of the state system. In this section we proceed to analyze the existence and uniqueness of solution for the state system (1) and (5).
3.1. The hydrodynamic model: Modified Navier-Stokes equations. Next, we will analyze the existence of solution of the hydrodynamic model (1). For this, we begin analyzing the existence of solution in the homogeneous Dirichlet case, then we analyze the existence for the case in which the recirculated flow is of constant sign, to finally give a result of existence for the general case. where Definition 3.1. Within the framework introduced at the beginning of this section, we will say that an element v ∈ W 1,∞,2 (0, T ; V, [L 2 (Ω)] 3 ) is a solution of the problem (10) if it satisfies the following conditions: in Ω.
2. v verifies the variational formulation: where we denote (see the monograph of Roubíček [16] for more details on the Sobolev-Bochner spaces): Next, we will obtain a set of estimates that are not only valid for the solutions of (10), but also for Galerkin approximation of the problem (technique that we will use to prove the existence of solution).

Lemma 3.2.
Let v be a solution of (10) in the sense of definition 3.1, then there exist continuous functions ∂v Proof. First, we consider z = v(t) as a test function. Thus, a.e. t ∈ (0, T ): Now, from the definition (11) of β and Korn's inequality, we have: where C 1 and C 2 are positive constants. Finally, taking into account estimates (16), (17) and Gronwall's Lemma, we can derive the following estimate: where Φ 1 is a continuous function. (It should be noted here that the L 3 boundedness of the velocity gradient is achieved by the presence of the term associated with turbulent viscosity. We will see, later on, that this boundedness is essential when obtaining additional regularity.) Let us take now z = ∂v/∂t(t) as a test function (this is possible since ∇ · z = 0 and z |∂Ω = 0). Arguing as in [9], we obtain: Now, if we use the following inequalities: we obtain: Integrating in the time interval (0, t): from where, using the definition (3) of operator D, we have: Using now Korn's inequality: where C 1 , C 2 , C 3 and C 4 are positive constants. On the other hand, where: It is clear now that (and here turbulent viscosity has a crucial effect on increasing the regularity of the solution: without turbulent viscosity the estimate could not be done): and, using Gronwall's Lemma, we obtain that: where Φ 2 is also a continuous function. That is, we are able to obtain continuous estimates with respect to data in the functional space which ends the proof.
of the problem (10), in the sense of definition 3.1, continuous with respect to data, that verifies the estimates (14) and (15).
Proof. The solvability of system (10) is proved, using the estimates obtained in Lemma 3.2, with the classical Galerkin method. That is, a sequence of functions is constructed, which solve finite-dimensional approximations of (10). Then, it is shown that a subsequence of {v n } converges to an element v that is a solution of the system (10). Proof. Let us assume that Then, substracting the variational formulations corresponding to v 1 and v 2 , we have that: However, on the one hand, and, on the other hand, by classical integration results, where C 1 and C 2 are two positive constants. Therefore, if we choose as a test function z = v(t), take into account the previous estimates, and integrate in [0, t] the equality (23), we obtain that: we obtain, by Korn's inequality, that , and V ⊂⊂ L q (Ω) for any 1 < q < ∞ (cf. Corollary 1.22 in [16]), in particular, v 1 ∈ L p (0, T ; [L q (Ω)] 3 ) for any 1 < p, q < ∞. Thus, bearing in mind that, for domains Ω ⊂ R 3 , we have (cf. [16]) that where C 4 and C 5 are positive constants.
, there will be a positive constant C 5 such that: 3.1.2. Existence and uniqueness of solution for the non-homogeneous case. Let us assume that the flow rates g k (t) are, for all values of k = 1, . . . , N CT , continuous functions of constant sign over the entire interval [0, T ]. We then analyze the following problem: where g k (0) = 0, for k = 1, . . . , N CT , v 0 ∈ V, and the surface velocity v S is assumed null for the sake of simplicity. In order to establish the concept of solution, we need to introduce the following space: Remark 2. We must note that V ⊂ W. Thus, considering the initial condition in V makes sense. In this way, it is also evident the need for the compatibility condition between the initial and the boundary conditions. Now, to prove the existence of solution of the problem with non-homogeneous Dirichlet conditions, we will use the techniques developed in [7]. That is, we will reconstruct from the values associated with the Dirichlet condition a function -in a suitable space -whose trace coincides with the given values. Once the above reconstruction is recovered, it will be used to reduce the search for a solution of the problem (25) to a homogeneous problem of the type (10). Remark 3. In the following we will commit a notation abuse, denoting by g n ∈ U to the element verifying that (i.e., we will denote equal to the space associated to the vector that to the space of its normal component, since its tangential component is zero): . Then, the vector satisfying that its normal component is equal to the vector g n will be denoted by g, that is, g = g n n ∈ U.
If we take into account the previous notation, we can reformulate the problem (25) in the following terms: The idea, as we commented before, is to reconstruct a function so that its trace coincides with the previous function g, and such that the regularity of the reconstruction is sufficient to reach that of the homogeneous solution. This reconstruction can be done on the basis of Theorem 2.2 of [7] which, in our particular case, establishes the following result, that we state in the form of Lemma: Lemma 3.5. There exists a linear, continuous mapping The problem now is that establishing a normal component constant in space may be incompatible with previous reconstruction theorem for the trace. To overcome this problem, we will consider the following regularization of the normal component: Definition 3.6. For the function g n associated to the Dirichlet condition, introduced in Remark 3, and for any value s ≥ 3/2, we will define its regularization Moreover, we will assume that this element g n is consistent with the initial conditions of the problem, that is, g n (0) = 0 a.e. x ∈ ∂Ω.
From now on, we will assume that the normal component of the velocity associated with the Dirichlet condition satisfies the regularity established in Lemma 3.5. For the sake of simplicity in the notation, we will denote by g that regularized normal component, and, given a value of s, by: the space associated to the Dirichlet conditions (clearly depending on s, although this is not shown in order to simplify the notation).
Lemma 3.7. We have the following continuous injections for s ≥ 7/2: . Proof. For s = 7/2, we only need to consider the embedding: As remarked at the beginning of the section, we look for a solution of problem (25) in the way where is the reconstruction of the trace g, as given in Lemma 3.5, and z ∈ W 1,∞,2 (0, is a function to be determined, solution of a homogeneous problem similar to (10).
, the reconstruction of the trace given in Lemma 3.5 for values of s ≥ 7/2. 2. z(0) = v 0 − ζ g (0) a.e. in Ω. 3. z verifies the following variational formulation: Remark 4. It is worthwhile noting here that, thanks to the regularity required for the reconstruction of the trace, the second member associated to the variational formulation (30) lies in L 2 (0, T ; [L 2 (Ω)] 3 ), that is: Moreover, it is straightforward that: Let z be a solution of (30), then, for s ≥ 7/2, there exist continuous functions Φ 1 and Φ 2 such that: Proof. We divide the demonstration into two parts: one for each estimate. a) First, in order to obtain the energy estimate, we consider η = z(t) as a test function. So, where C 1 is a positive constant. Now, in one hand, Thus, fitting the values of in previous expression, and using triangle inequality for the norms and Korn's Lemma, we obtain the existence of a positive constant C 2 such that: and, consequently, from Gronwall's Lemma, we derive the existence of a continuous function Φ 1 verifying: b) Second, for deriving the estimate on the time derivative, we will consider η = ∂z/∂t(t) as a test function. Then, arguing as in the homogeneous case we have: , we have the following estimates: Integrating in the time interval (0, t): β(ζ g (s) + z(s)) (ζ g (s) + z(s)) : ∂ζ g ∂t (s) dx ds, from where: where C 1 and C 2 are positive constants. We only have to estimate the last two terms of the previous inequality. So, for the first one it is very similar to the homogeneous case, that is,  (ζ g (s) + z(s)) :

ds,
with C 5 a positive constant. Taking into account all above estimates, we can establish the following: where C 6 , C 7 and C 8 are positive constants. Finally, using triangular inequality, Korn's Lemma and Gronwall's Lemma, we can conclude the existence of a continuous function Φ 2 , verifying the desired estimate:  (25), continuous with respect to data, that verifies the following estimates: with Φ 3 and Φ 4 continuous functions.
Proof. Using the estimates of Lemma 3.9 combined with the Galerkin method, we can obtain the existence of z ∈ W 1,∞,2 (0,  We are now ready to analyze the existence of solution of the general problem (1). Logically, the existence of solution of the general problem can only be achieved for certain pumping/turbination actuation regimes. That is, we will start from the base on which only a finite number of changes in the regime of operation of the pumps takes place. In order to simplify the exposure of the results, and without losing generality, we will assume that all pumps change their operating mode in unison. It is always possible to consider such a situation by introducing an appropriate refinement of time intervals in which it is possible for all groups to keep unchanged its regime during each time period. Let us assume then that g = (g 1 , . . . , g N CT ) ∈ [C([0, T ])] N CT are functions verifying that there exist a finite number of time instants t 1 , . . . , t R ∈ (0, T ) with 0 = t 0 < t 1 < . . . < t R < t R+1 = T , such that: where c k ∈ {0, 1, −1} is a constant depending on the pumping operation mode of the group. The problem is reduced, therefore, to solve in a chained way the hydrodynamic model in each one of the intervals: Based on the results obtained in the previous subsection, we can guarantee the existence and uniqueness of the solution in each one of the time intervals, verifying the estimates (38)-(39). It should be noted that the time continuity of the solution in space W is a key element to demonstrate the existence of solution. We recall again that the only restriction we ask the set of pumps is that the operation mode changes occur at the same time: this restriction is met naturally if we require that the groups can only change state a finite number of times in the whole time [0, T ]. Then, we define the function v ∈ W 1,∞,2 (0, It is evident that the solution constructed in this way is the only solution of the hydrodynamic model (1). Moreover, it is continuous with respect to the data in the sense that it verifies estimates analogous to those obtained in (38)-(39). So, we have demonstrated the following existence/uniqueness result:  (1), continuous with respect to data, that verifies the estimates (38) and (39).

The eutrophication model: Michaelis-Menten kinetics.
In this section we will analyze the existence and uniqueness of the solution to the problem (5), for this, we will start by analyzing a certain type of parabolic equation that will later be used as the basis for the study of the existence of solution of the above problem.
3.2.1. Analysis of a particular type of parabolic equation. Let us start by analyzing the existence and uniqueness of solution of the following particular equation that we will use to analyze the existence of solution of the eutrophication model. Basically, this is a reaction-diffusion-convection type equation affected, through the boundary conditions associated to collectors and injectors, by the operation mode of the pumps. We will assume, without loss of generality, that all groups are in turbination mode (the other cases are completely analogous), that is, taking into account Remark 1, we are going to consider the following equation: where h k (t) is a function only depending on time (constant in space) representing the concentration in the injectors.
To look for solutions of the problem (42) in space W 1,2,2 (0, T ; H 1 (Ω), H 1 (Ω) ), we take as the space of test functions for the variational formulation associated with the problem (42). We have the following technical result: Lemma 3.13. The following injections are compact: • W 1,2,2 (0, T ; H 1 (Ω), H 1 (Ω) ) ⊂⊂ L Proof. The demonstration for the first inclusion is an immediate consequence of the compactness of injection of H 1 (Ω) into L 6− (Ω) and of the interpolation Lemma 7.8 of [16].
In order to prove the second and third compact inclusions, we can use the Aubin-Lions Lemma 7.7 of [16], along with the compact inclusions of H 1 (Ω) into L 6− (Ω) and L 4− (∂Ω).
Definition 3.14. We will say that an element u ∈ W 1,2,2 (0, T ; H 1 (Ω), H 1 (Ω) ) is a solution of problem (42) if it satisfies the following conditions: where φ 6 is a positive constant, possibly depending on upper bound k max .
Proof. The demonstration for the case of all groups in turbination mode can be performed using standard techniques for parabolic equations. An analogous result for the case of homogeneous Neumann boundary conditions has been proved -employing these techniques -by the authors in their recent work [4]. Interested readers can consult further details, for instance, in the monograph of Roubíček [16].
The demostration for the other cases (with groups also pumping or at rest) is completely similar to the previous one, since the only changes are related to the particular location of the nonhomogeneous Dirichlet and the homogeneous Neumann conditions, but the other issues are exactly the same.

Analysis of the model for the different species.
In this section we will demonstrate the existence and uniqueness of solution for the system of equations (5). As in the case of the hydrodynamic model, we will start by analyzing the case in which all groups are turbinating (g k (t) > 0 for all k = 1, . . . , N CT ), being analogous the case where there are groups pumping or at rest. (In case where all groups are at rest, we refer the reader to the reference [4]). Once analyzed the case where all the groups maintain their operation mode throughout the time period under study, we will analyze the general case.
Then, let us begin by analyzing the existence and uniqueness of solutions for the eutrophication model, assuming that all groups are turbinating, that is, for i = 1, . . . , 5: where is the solution of the hydrodynamic model (25), f i ∈ L 2 (0, T ; L 2 (Ω)) represents the source term for each species, and initial condition u i 0 lies in a regular enough subspace of L 2 (Ω), for i = 1, . . . , 5.
3. u i verifies the variational formulation: , with φ 7 a positive constant depending on the parameters of the problem.
• u is locally Lipschitz continuous with respect to data.
Proof. In order to demonstrate the existence of solution, we consider the mapping: where B is the set defined by:  ∈ (0, T ), ∀i = 1, . . . , 4, ∀k = 1, . . . , N CT , and the image of K is given by: 5 ) the solution of the following uncoupled system of equations (solved in the order where the reaction terms are given by: , and the source term is given by: The non-negativity of u Ω is a direct consequence of the Theorem 3.15, and for proving the non-negativity of u T we only have to regularize the equation (42) in order to obtain a continuous solution and then pass to the limit. Now we will show that the mapping K has a fixed point in a certain subset B of B. On the one hand, thanks to the Theorem 3.15 (for the equations corresponding to nutrients, phytoplankton, zooplankton and organic detritus), and to classic results for parabolic equations (for the equation associated with dissolved oxygen), we have that the application K is well defined.
In addition, following the order of resolution mentioned above, we can obtain the following estimates: • For phytoplankton: • For zooplankton: • For organic detritus: • For generic nutrient: • For dissolved oxygen: Above constants C j , for j = 1, . . . , 10 can possibly depend on the initial conditions, on the water velocity, and on the source terms. In addition, if we assume that they satisfy the compatibility conditions: then, by simple arithmetic operations, we can assure the existence of new positive constants C j , for j = 1, . . . , 10, such that, if we consider the subset B of B defined by: To prove this, if we consider an element (z Ω , z T ) ∈ B, we have that (u Ω , u T ) = K(z Ω , z T ) verifies the following estimates: Then, we can take the constants C j , for j = 1, . . . , 10, defined as the solution of the following 10 × 10 linear system: (1 + C 7 + C 8 + C 9 + C 10 + C 3 ) = C 5 , C 6 C 1 (1 + C 6 + C 7 + C 8 + C 9 + C 3 ) = C 6 , C 7 C 2 (1 + C 7 + C 3 ) = C 7 , C 8 C 3 (1 + C 8 ) = C 8 C 9 C 4 (1 + C 7 + C 8 + C 9 + C 3 ) = C 9 , C 10 C 5 (1 + C 7 + C 8 + C 9 + C 10 + C 3 ) = C 10 . Previous system can be written as a linear system Ax = b, where x = ( C 1 , . . . , C 10 ) T , By solving explicitely the system, we obtain the following solution: , , , .
We can observe that, since C j C 5+j − 1 < 0, ∀j = 1, . . . , 5, then C i > 0, ∀i = 1, . . . , 10, and then its clear that Moreover, the continuity of the application K can be demonstrated by techniques analogous to those described in [4], and its relative compactness is a straightforward consequence of the compact inclusions demonstrated in the Lemma 3.13. Thus, we can apply the Schauder fixed-point Theorem, which guarantees the existence of a fixed point for the application K that is, by its definition, a solution of the state system (46).
Finally, for the uniqueness of solution and its continuity with respect to data, we refer the reader to the reference [4], since their demonstration can performed using analogous techniques.

3.2.3.
Analysis of the model for the general case. The existence and uniqueness of solution for the general model (5) -that is, the case in which some pairs collectorinjector are pumping, other ones turbinating, and other ones at rest -can be proved following the same idea as that developed in the section 3.1.3, mainly due to the fact that W 1,2,2 (0, T ; H 1 (Ω), H 1 (Ω) ) ⊂ C([0, T ]; L 2 (Ω)), and that the initial conditions u i 0 are taken in L 2 (Ω), for i = 1, . . . , 5.
4. Mathematical analysis of the optimal control problem. Let us assume, without loss of generality, that all the groups are in turbination mode (the generalization to the rest of the cases does not present problems). We begin by establishing the precise formulation of the control problem (P) in the following terms: (P) min{J(g) : g ∈ U * ad }, where: • Control constraints: As we said at the beginning of these work, we will act on the pumped/turbinated flow for each of the pairs collector-injector. Since we are assuming that we are in a situation in which all groups are turbinating, we will consider the following set of admissible controls g (normal component of velocity on the boundary): where c 1 , c 2 , c 3 > 0 are technological constants, and the parameter s ≥ 7/2. Given a normal velocity on the boundary g ∈ U ad , we denote by g = gn + 0τ the associated velocity, and by ζ g ∈ W 1,2,2 (0, 3 ) its corresponding reconstruction, as given in Lemma 3.5.
• State equations: Given an admissible control g ∈ U ad , we denote by the solutions of, respectively, the state systems (25) and (46), in the sense of Lemma 3.8 and Definition 3.16. • Cost functional: We consider the following cost function: • State constraints: We define the mapping: that is, in order to evaluate G, it is necessary to consider the following chain of compositions: Then, we consider the following state constraints: where K C ⊂ R 5 is the rectangle given by , and assume that the set U * ad = {g ∈ U ad : G(g)(t) ∈ K C , ∀t ∈ [0, T ]} = ∅. So, we can demonstrate the following existence result for the control problem: and is the solution of the variational formulation: So, there exist subsequences -still denoted in the same way for simplicity -such that verify the following convergences: • g n g in H s−1/2 (∂Ω × (0, T )) ∩ L 2 (0, T ; H 1 (∂Ω)).
• u n → u in C w ([0, T ]; L 2 (Ω)). Now, on the one hand, based on the strong convergence of g n to g, on its weak convergence in U, and on the weak lower semicontinuity of the norm, we have that g ∈ U ad .
On the other hand, the pass to the limit in the terms associated with the variational formulations (65) and (66) can be performed based on above convergences, using standard techniques similar to which can be found in [10], [3] and [2]. So, we obtain that v = ζ g + z and u are the states associated to g.
Finally, passing to the limit in the cost functional J and in the state constraints given by G is straightforward, due to above convergences. Thus, we have demonstrated that the admissible control g ∈ U * ad is a (not necessarily unique) solution of the optimal control problem (P), with associated states v and u.

5.
Numerical solution of the optimal control problem. In this section we will see how to solve numerically the control problem (P). For this, we will consider a space-time discretization combining the method of the characteristics (for the time discretization) and the finite element method (for the space discretization). For the numerical resolution of the nonlinear coupling between the velocity components and those of the eutrophication model we will use a fixed-point algorithm. Finally, the numerical resolution of each of the Stokes-type problems derived from the fixedpoint algorithm associated with the hydrodynamic model will be carried out using a penalty method [18,5,6].
For the discretization of the problem, let us consider a partition 0 = t 0 < t 1 < . . . < t N = T of the time interval [0, T ] such that t n+1 − t n = ∆t = 1/α, ∀n = 0, . . . , N − 1, and a family of meshes for the domain Ω with characteristic size h. Associated to this family of meshes, we also consider three compatible finite element spaces W h , M h and Z h corresponding, respectively, to the velocity and the pressure of water, and to the concentrations of the species involved in the eutrophication model.
In order to simplify the presentation of the numerical problem, we will assume that all the groups collector-injector are in turbination mode, with the controls only depending on the time variable. Therefore, we try to solve the following constrained nonlinear minimization problem resulting from discretizing the control problem (P): (Q) min{J(g) : g ∈ U * ad }, where: • Discretized control constraints: As already commented, we will assume that the inlet/outlet velocities for the collectors and injectors is constant at each step of time. We will then consider the following set of admissible discrete controls: with c a bound related to pumps' technical characteristics. • Discretized state equations: Assumed known the initial velocity v 0 ∈ W h and the initial concentrations for the species u 0 ∈ Z h , and given an admissible control g = (g 1,1 , g 2,1 , . . . , g N CT ,1 g 1 , . . . , g 1,N , g 2,N , . . . , g N CT ,N we compute the associated state in the following way: (i) The hydrodynamic model: For each n = 0, 1, . . . , N − 1, the pair veloc- is the solution of: where the space , and a penalty parameter λ > 0 is considered in order to regularize the problem (see [8], for instance).
(ii) The eutrophication model: For each n = 0, . . . , N , the concentration u n+1 ∈ Z h , with is the solution of: where the space X h = {η ∈ Z h : η |∪ N CT k=1 T k = 0} and the discrete characteristic (u n • X n )(x) = u n (x − ∆t v n (x)).

Remark 6.
It should be noted that we employ a different number of time steps in both models, that is, in the eutrophication model one more step of time is solved than in the hydrodynamic model. This shift is motivated by the following dependency scheme: It is easy to see that, due to the time discretization used, the control associated to the first step of time, begins to have influence from the second step of time for the species of the eutrophication model. If we consider the same number of steps of time for the two models, the control associated with the last time step has no effect on the state constraints. To avoid this inconsistency, we propose to solve the eutrophication model an additional step of time, and to displace one time step the state constraints. Finally, as a direct consequence of above scheme, the Jacobian associated to state constraints is a block lower triangular matrix: From the computational viewpoint, for the generation of the mesh associated to the domain and for the numerical resolution of the state equations we have used FreeFem++ [12], and for the numerical resolution of the nonlinear minimization problem, the interior-point algorithm IPOPT [19]. In final section we give further details associated to the numerical test that we present there.
One of the necessary requirements to use the IPOPT algorithm in the numerical resolution of the nonlinear minimization problem (Q) is to have a function that evaluates the Jacobian associated to the state constraints (71) of our problem. For this, there are two possibilities, to use the linearized equations or to use the adjoint state equations. The choice of one method or another depends on the relationship between the size of the control and the constraints. That is, if we want to use the linearized equations to calculate the Jacobian matrix, we must solve N CT × N times the linearized equations, whereas if we employ the adjoint state to calculate the Jacobian matrix, we must solve 5N times the adjoint state system. If we take into account that the computational cost associated with the resolution of the linearized equations and the adjoint state is similar, we have that if N CT < 5 it is more convenient to use the linearized equations, whereas if N CT > 5, it is more advantageous to employ the adjoint state. In the example presented in the next section, we have considered four pairs collector-injector, therefore, we will use the linearized equations for computing the Jacobian of the discretized state constraints.
So, using the linearized state equations, the Jacobian associated to the state constraints can be computed by: where δv 0 = 0, δu 0 = 0 and, given a direction δg, the linearized state is obtained by: (i) The linearized hydrodynamic model: For each n = 0, 1, . . . , N − 1, the pair is the solution of: (ii) The linearized eutrophication model: For each n = 0, . . . , N , the element δu n+1 ∈ Z h , with is the solution of: 6. Numerical results. In this final section we will present a few numerical results that we have obtained by solving our problem in a realistic scenario that, for the sake of simplicity, corresponds to a two-dimensional case. For this purpose, we have considered a rectangular domain Ω = [0, 14] × [0, 16] (measured in meters), corresponding to a reservoir, in which we have distributed N CT = 4 pairs collector-injector with the geometric configuration shown in Figure 1. The control domain we have chosen is the lower rectangle (shaded in green in Figure 1) of dimensions Ω C = [0, 14] × [0, 2] ⊂ Ω. The parameters used for the numerical resolution of the state equations are given in Table 1.
For time discretization we have considered a time step of ∆t = 3600 seconds (1 hour), and for space discretization we have used a regular mesh formed by triangles of characteristic size h = 0.25 meters (corresponding to 3075 vertices). The finite element spaces we have employed for space discretizations have been the minielement (P 1 -bubble/P 1 ) for the hydrodynamic model, and the Lagrange P 1 element for the eutrophication model. In Figure 2 we compare the evolution of the state 308 A. MARTÍNEZ, F. J. FERNÁNDEZ AND L. J. ALVAREZ-VÁZQUEZ 3.00e +01 m 1 ' 2 0.00e +00 m 1 por mgC/l I 0 1.00e +02 Cal/m 2 por s I S 1.00e +02 Cal/m 2 por s C fz 1.10e +00 -C T 1.06e +00µ i 1.00e 05 m 2 /s Para las discretizaciones temporales hemos considerado un paso de tiempo de 3600 segundos (1 hora) y, para la discretización espacial, hemos considerado un mallado regular formado por triángulos de longitud característica 1/4 metros (3075 vértices). Los espacios de elementos finitos que hemos empleado para las discretizaciones espaciales han sido el minielemento (P 1b , P 1 ) para el modelo hidrodinámico y P 1 para 37 Figure 1. Physical domain Ω for the numerical example, showing the control domain Ω C (shaded area) and the four pumping groups connecting collectors C k with injectors T k , for k = 1, . . . , 4.  This decay in concentration is mainly due to the decomposition of organic detritus which, by the effects of sedimentation, are deposited in the bottom. In order to mitigate this oxygen decrease, it is observed that the pumping of water from upper to lower layers relieves the effects discussed above. However, a side effect of this turbination is the increase in the rate of deposition of the organic detritus in the lower layers, which may aggravate the problem when we stop acting on the system. In Figures 3 and 4 we show the concentrations of dissolved oxygen and of organic detritus in the last step of time (n = 49), where previous behaviour can be observed (that is, with mechanical mixing of water, stratification is completely broken for both species, and aeration of the bottom layer is accomplished).

Parameters Values Units
In this first approach to solving the problem, we are only interested in controlling the concentration of dissolved oxygen in the lower layer, for which we consider the state constraints rectangle given by K C = (−∞, ∞) × (−∞, ∞) × (−∞, ∞) × (−∞, ∞)×, [λ m 5 , ∞), where λ m 5 ∈ R N is such that λ n,m 5 = min{4.6, G 5,n (g ref )}, for n = 1, . . . , N , (this is, we consider as a lower threshold for the dissolved oxygen that one obtained by turbinating water at medium power until reaching a concentration    4.6 mg/l). It is then expected that the optimal control g opt is lower, from the point at which the threshold of 4.6 mg/l is exceeded, than that used to generate previous threshold, with the consequent, desired energy saving. Now, let us present the results we have obtained by optimizing the process during the first 12 steps of time. In Figure 5 we can see the optimal control g opt obtained for the four groups against the reference control g ref = 5.0 10 4 m/s, and the evolution of the constraint associated with the dissolved oxygen for the optimal case, for the reference case, and for the uncontrolled case. We observe that the control is optimal in the sense that it saturates the constraint associated with the dissolved oxygen. On the other hand, it can be also noted that, with this optimal configuration, groups      Figure 6. Concentration of dissolved oxygen at final time for the optimal control g opt .