Controlling eutrophication by means of water recirculation: An optimal control perspective

Journal of Computational and Applied Mathematics 421 (2023) 114886 Contents lists available at ScienceDirect Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam Controlling eutrophication by means of water recirculation: An optimal control perspective Francisco J. Fernández a,∗, Aurea Martínez b, Lino J. Alvarez-Vázquez b aUniversidade de Santiago de Compostela, Instituto de Matemáticas, 15782 Santiago, Spain bUniversidade de Vigo, E.I. Telecomunicación, 36310 Vigo, Spain article info Article history: Received 6 July 2021 Received in revised form 14 July 2022 MSC: primary 49J20 65K10 secondary 49M25 Keywords: Optimal control Numerical optimization Artificial recirculation Eutrophication abstract In this work, the artificial recirculation of water is presented and analyzed, from the perspective of the optimal control of partial differential equations, as a tool to prevent eutrophication effects in large waterbodies. A novel formulation of the environmental problem, based on the coupling of nonlinear models for hydrodynamics, water temperature and concentrations of the different species involved in the eutrophication processes, is introduced. After a complete and rigorous analysis of the existence of optimal solutions, a full numerical algorithm for their computation is proposed. Finally, some numerical results for a realistic scenario are shown, in order to prove the efficiency of our approach. ©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction: The environmental problem One of the most important environmental problems related to human activities near large masses of water is Eutrophication. Eutrophication is caused by an excessive supply of nutrients (normally nitrogen and phosphorus) to the water. These nutrients cause undesirable effects like algal blooms that directly affects the concentration of dissolved oxygen in the deeper layers since the processes of remineralization of organic detritus consume oxygen, which can lead to oxygen depletion of the body of water [1]. In Fig. 1 (left) we can find a schematic representation of the problem and its consequences. Artificial circulation is a management technique that increases the dissolved oxygen concentration in the bottom layers. Water from the well aerated upper layers is taken by means of a set of collectors and injects it into the bottom layers that are poorly oxygenated. Then, oxygen-poor water from the bottom is circulated to the surface, where oxygenation from the atmosphere and photosynthesis can naturally occur [2,3]. In Fig. 1 (right) we can find a representation of the main idea of water artificial circulation. Although eutrophication has received some attention from the mathematical viewpoint in last decade (see, for instance, the recent publications [4–6] and the references therein), the study of artificial circulation as a eutrophication control tool has remained unaddressed in the mathematical literature up to now, as far as we know. We can only mention two recent papers of the authors [2,3], where a simplified preliminary formulation of the problem is posed and briefly analyzed. The main contribution of this work with respect to [3] is that in this case, we are considering the effect that solar radiation has on the increase in water temperature and, therefore, on the evolution of the different species. The incorporation of ∗Corresponding author. E-mail addresses: [email protected] (F.J. Fernández), [email protected] (A. Martínez), [email protected] (L.J. Alvarez-Vázquez). https://doi.org/10.1016/j.cam.2022.114886 0377-0427/©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons. org/licenses/by-nc-nd/4.0/). F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 Fig. 1. On the left side, a diagram representing algal blooms caused by eutrophication and its consequences. On the right side, a basic scheme depicting the water artificial circulation process [2]. solar radiation to the model represents a notable increase in complexity in the mathematical analysis of the existence of a solution to the state equations and also in mathematical analysis of the control problem and its numerical resolution. This complexity is revealed in [2], where the authors show that, under certain assumptions about the regularity of the data, it is possible to find a solution to the state equations. We must mention that in this work we will also present a complete algorithm for the numerical resolution of the control problem and we will present results that we have obtained using realistic data. In next section we present a detailed mathematical formulation of the physical problem as a control/state constrained optimal control problem of nonlinear partial differential equations. Then, we briefly recall the existence results for the state system proved in [2]. Next, in the central part of the paper, we demonstrate in a rigorous way the existence of optimal solutions. Finally, we present the numerical resolution of the problem, introducing a full computational algorithm and a realistic numerical example, showing the efficiency of our approach. 2. Mathematical formulation of the control problem In this section we will formulate the environmental problem in the framework of optimal control of partial differential equations. For a better understanding of this novel mathematical formulation, we will divide this section into five subsections: in the first subsection we will introduce and describe the physical domain; in the second one, the control variables (in our case, the volumetric flow rate for each pump); in the third subsection we will establish the mathematical formulation for the thermo-hydrodynamic model; in the fourth one we will present the eutrophication model that will be used (the core of our model) and, finally, in the fifth subsection we will formulate the optimal control problem. In order to establish the appropriate framework for mathematically analyzing the control problem, we consider, for a Banach space V1and a locally convex space V2such that V1⊂V2, and for 1 ≤p,q≤ ∞, the following Sobolev–Bochner space (cf. Chapter 7 of [7] for further details): W1,p,q(0,T;V1,V2)={u∈Lp(0,T;V1):du dt ∈Lq(0,T;V2)},(1) where du dt denotes the derivative of uin the sense of distributions. It is well known that, if both V1and V2are Banach spaces, then W1,p,q(0,T;V1,V2) is also a Banach space endowed with the norm ∥u∥W1,p,q(0,T;V1,V2)= ∥u∥Lp(0,T;V1)+du dt Lq(0,T;V2). 2.1. The physical domain We consider a domain Ω⊂R3corresponding, for instance, to a reservoir. In order to promote the artificial circulation of water inside the domain Ω, we suppose the existence of a set of NCT pairs collector–injector {(Ck,Tk)}NCT k=1⊂∂Ωin such a way that each water collector is connected to its corresponding injector by a pipe with a pumping group. We also assume a smooth enough boundary ∂Ω, such that it can be split into four disjoint subsets ∂Ω=ΓS∪ΓC∪ΓT∪ΓN, where ΓCcorresponds to the part of the boundary where the water collectors are located (ΓC= ∪NCT k=1Ck), ΓTcorresponds to the part of the boundary where the water injectors are located (ΓT= ∪NCT k=1Tk), ΓSis the top part of the boundary in contact 2 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 Fig. 2. Geometrical configuration of an example domain Ωwith NCT =4 collector/injector pairs, showing the different boundary sections: ΓS, ΓC= ∪4 k=1Ck,ΓT= ∪4 k=1Tkand ΓN, and also the control domain ΩC. with air, and ΓN=∂Ω\(ΓS∪ΓC∪ΓT)corresponds to the rest of the boundary. In particular, we suppose the boundary ∂Ωregular enough to assure the existence of elements ϕk,˜ϕk∈H3/2(∂Ω), for k=1,...,NCT , satisfying the following assumptions (mainly corresponding to suitable regularizations of the indicator functions of Tkand Ck, respectively): •ϕk(x),˜ϕk(x)≥0, a.e. x∈∂Ω, •ϕk(x)=0, a.e. x∈∂Ω\Tk, and ∫Tkϕkdγ=µ(Tk), •˜ϕk(x)=0, a.e. x∈∂Ω\Ck, and ∫Ck˜ϕkdγ=µ(Ck), where µ(S) represents the (n−1) dimensional measure of a generic set S, and β0:u∈H3/2(∂Ω)→β0(u)∈H2(Ω) denotes a right inverse of the classical trace operator γ0, i.e., (γ0◦β0)(u)=u(cf. Theorem 8.3. in Chapter 1 of [8]). Finally, we also consider a subdomain ΩC⊂Ω, corresponding to the part of the domain where we want to increase the dissolved oxygen concentration (denoted as control domain in Fig. 2). 2.2. The control variable As above commented, our control will be the volumetric flow rate (m3s−1) by pump kat each time t,gk(t)∈H1(0,T), for k=1,...,NCT , where T(s) denotes the length of the time interval. We will suppose that the control acts over the system through a Dirichlet boundary condition on the hydrodynamic model: v=φgon ∂Ω×(0,T),(2) where v(x,t) will denote the water velocity, and where: φg(x,t)= NCT ∑ k=1 gk(t)[ϕk(x) µ(Tk)−˜ϕk(x) µ(Ck)]n(x) (3) represents the given Dirichlet condition for the hydrodynamic system. It is immediate that, thanks to the regularity of the control gand of the functions {(ϕk,˜ϕk)}NCT k=1, we have that φg∈W1,2,2(0,T;H3/2(∂Ω),H3/2(∂Ω)) (cf. expression (1) for a detailed definition of this Sobolev–Bochner space), and also that ∫∂Ω φg·ndγ=0. 2.3. The thermo-hydrodynamic model We denote by v(x,t) (m s−1) the solution of the following modified Navier–Stokes system with a Smagorinsky model of turbulence: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ∂v ∂t+ ∇v v −div(Ξ(v)) + ∇p=α0(θ−θ0)agin Ω×(0,T), ∇ · v=0 in Ω×(0,T), v=φgon ∂Ω×(0,T), v(0) =v0in Ω, (4) 3 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 where ag(m s−2) is the gravity acceleration, α0= − 1 ρ ∂ρ ∂θ (K−1) is the thermic expansion coefficient (see, for instance, section 10.7 of [9]), ρ(θ) (g m−3) is the density, v0is the initial velocity, and the boundary field φgis the element given by (3). The diffusion term Ξ(v) is given by: Ξ(v)=∂D(e) ∂e⏐⏐⏐⏐e=e(v) ,with e(v)=1 2(∇v+ ∇vt),(5) where Dis a potential function (for instance, in the standard case of the classical Navier–Stokes equations, D(e)=ν[e:e], with ν(m2s−1) the kinematic viscosity of the water, and, consequently, Ξ(v)=2νe(v)). However, in our case, the Smagorinsky model, the potential function is defined as in [10]: D(e)=ν[e:e]+2 3νtur [e:e]3/2,(6) so, Ξ(v)=∂D(ϵ) ∂ϵ ⏐⏐⏐⏐ϵ=ϵ(v) =2ν ϵ(v)+2νtur [ϵ(v):ϵ(v)]1/2ϵ(v) =(2ν+2νtur [ϵ(v):ϵ(v)]1/2)ϵ(v)=β(ϵ(v)) ϵ(v), (7) where β(ϵ(v)) =2ν+2νtur [ϵ(v):ϵ(v)]1/2and νtur (m2) is the turbulent viscosity. Regarding thermic effects, water temperature θ(x,t) (K) is the solution of the following convection–diffusion partial differential equation with nonhomogeneous, nonlinear, mixed boundary conditions: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂θ ∂t+v· ∇θ− ∇ · (K∇θ)=0 in Ω×(0,T), θ=φθon ΓT×(0,T), K∂θ ∂n=0 on ΓC×(0,T), K∂θ ∂n=bN 1(θN−θ) on ΓN×(0,T), K∂θ ∂n=bS 1(θS−θ)+bS 2(T4 r− |θ|3θ) on ΓS×(0,T), θ(0) =θ0in Ω, (8) where Dirichlet boundary condition φθis given by expression: φθ(x,t)= NCT ∑ k=1 ϕk(x)∫T −T ρϵ(t−ϵ−s)γk θ(s)ds (9) with, for each k=1,...,NCT , γk θ(s)=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1 µ(Ck)∫Ck θ0dγif s≤0, 1 µ(Ck)∫Ck θ(s)dγif s>0, (10) representing the mean temperature of water in the collector Ck, and with the weight function ρϵdefined by: ρϵ(t)=⎧ ⎨ ⎩ c ϵexp (t2 t2−ϵ2)if |t|< ϵ, 0 if |t| ≥ ϵ, (11) for c∈Rthe positive constant satisfying the unitary condition: ∫R ρ1(t)dt =1. In other words, we are assuming that the mean temperature of water at each injector Tkis a weighted average in time of the mean temperatures of water at its corresponding collector Ck. In order to obtain the mean temperature at each injector, we convolute the mean temperature at the collector with a smooth function with support in (t−2ϵ, t). In this way, we have that the temperature in the injector only depends on the mean temperatures in the collector along the time interval (t−2ϵ, t). Parameter 0 < ϵ < Trepresents, in a certain sense, the technical characteristics of the pipeline that define the stay time of water in the pipe. We also suppose that there is not heat transfer through the walls of the pipelines (that is, they are isolated). 4 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 Moreover, for the other terms appearing in the formulation of problem (8) we have that: •nis the unit outward normal vector to the boundary ∂Ω. •K>0 (m2s−1) is the thermal diffusivity of the fluid, that is, K=α ρcp, where α(W m−1K−1) is the thermal conductivity, ρ(θ) (g m−3) is the density, and cp(W s g−1K−1) is the specific heat capacity of water. •bK 1≥0 (m s−1), for K∈ {N,S}, are the coefficients related to convective heat transfer through the boundaries ΓNand ΓS, obtained from the relation ρcpbK 1=hK, where hK≥0 (W m−2K−1) are the convective heat transfer coefficients on each surface. •bS 2>0 (m s K−3) is the coefficient related to radiative heat transfer through the boundary ΓS, given by bS 2=σBε ρcp, where σB(W m−2K−4) is the Stefan–Boltzmann constant and εis the emissivity. •θ0≥0 (K) represents the initial temperature. •θN, θS≥0 (K) are the temperatures related to convection heat transfer on the surfaces ΓNand ΓS. •Tr≥0 (K) is the radiation temperature on the surface ΓS, derived from the expression σBεT4 r=(1 −a)Rsw,net + Rlw,down, where ais the albedo, Rsw,net (W m−2) denotes the net incident shortwave radiation on the surface ΓR, and Rlw,down(W m−2) denotes the downwelling longwave radiation. 2.4. The eutrophication model We consider the following system for modeling the eutrophication processes, based in Michaelis–Menten kinetics (further details can be found, for instance, in [11,12] and the references therein), where we consider the concentrations of five different species: u1(x,t) (mg l−1) stands for the nutrient (nitrogen in this case), u2(x,t) (mgC l−1) for the phytoplankton, u3(x,t) (mgC l−1) for the zooplankton, u4(x,t) (mgC l−1) for the organic detritus, and u5(x,t) (mg l−1) for the dissolved oxygen: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ui ∂t+v· ∇ui− ∇ · (µi∇ui)=Ai(x,t, θ, u) in Ω×(0,T), ui=φuion ΓT×(0,T), µi∂ui ∂n=0 on (ΓS∪ΓN∪ΓC)×(0,T), ui(0) =u0,iin Ω,i=1,...,5, (12) where, for i=1,...,5, φui(x,t)= NCT ∑ k=1 ϕk(x)∫T −T ρϵ(t−ϵ−s)γk ui(s)ds,(13) and, for k=1,...,NCT , and i=1,...,5, γk ui(s)=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1 µ(Ck)∫Ck u0,idγ , if s≤0, 1 µ(Ck)∫Ck ui(s)dγ , if s>0. (14) Finally, the reaction term A=(Ai):Ω×(0,T)×R6→R5is defined by the following expression: A(x,t, θ, u)= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ −Cnc L(x,t, θ) KN+ |u1|u1u2+Cnc Kru2+Cnc KrdD(θ)u4 L(x,t, θ) KN+ |u1|u1u2−Kru2−Kmf u2−Kz KF+ |u2|u2u3 CfzKz KF+ |u2|u2u3−Kmzu3 Kmf u2+Kmzu3−KrdD(θ)u4 Coc L(x,t, θ) KN+ |u1|u1u2−Coc Kru2−Coc KrdD(θ)u4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (15) where: •Coc ≥0 (mg mgC−1) is the oxygen–carbon stoichiometric relation, •Cnc ≥0 (mg mgC−1) is the nitrogen–carbon stoichiometric relation, •Cfz ≥0 is the zooplankton grazing efficiency factor, •Krd ≥0 (s−1) is the detritus regeneration rate, 5 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 •Kr≥0 (s−1) is the phytoplankton endogenous respiration rate, •Kmf ≥0 (s−1) is the phytoplankton death rate, •Kmz ≥0 (s−1) is the zooplankton death rate (including predation), •Kz≥0 (s−1) is the zooplankton predation (grazing), •KF>0 (mgC l−1) is the phytoplankton half-saturation constant, •KN>0 (mg l−1) is the nitrogen half-saturation constant, •µi≥0 (m2s−1), i=1,...,5, are the diffusion coefficients of each species, •Dis the thermic regeneration function for the organic detritus, defined as: D(θ)=Θθ−θ0,(16) with log(Θ) (K−1) the thermic regeneration constant for the reference temperature θ0. In order to simplify the mathematical analysis of the state equations we will consider the following linear approximation: D(θ)=1+log(Θ)(θ−θ0) (17) if Θ>0, and D(θ)=1 if Θ=0. •Lis the luminosity function, given by: L(x,t, θ)=µCθ−θ0 t I0(t) Is e−ϕ1x3,(18) with I0(W m−2) the incident light intensity, Is(W m−2) the light saturation, log(Ct) (K−1) the phytoplankton growth thermic constant for the reference temperature θ0,ϕ1(m−1) the light attenuation due to depth, and µ(s−1) the maximum phytoplankton growth rate. Again, for the sake of simplicity, we will consider the following linear approximation: L(x,t, θ)=µ(1+log(Ct)(θ−θ0))I0(t) Is e−ϕ1x3(19) if Ct>0, and L(x,t, θ)=µI0(t) Ise−ϕ1x3if Ct=0. 2.5. The optimal control problem Our main objective is to ensure that the concentration of dissolved oxygen in the bottom layer lies within an admissible range of values by means of an optimal artificial circulation of water from the well-aerated upper layer. So, we want to solve the following optimal control problem: (P) min{J(g):g∈Uad,1 µ(ΩC)∫ΩC u5(t)dx∈ [λm, λM]}, where Uad = {g∈ [H1(0,T)]NCT :g(0) =0,∥gk∥H1(0,T)≤c,∀k=1,...,NCT }(20) is the admissible set, c>0 is a constant related to technological limitations of the pumps, J(g) is the cost function: J(g)=1 2 NCT ∑ k=1∫T 0 gk(t)2dt +1 2 NCT ∑ k=1∫T 0 dgk dt (t)2dt,(21) and λm, λM>0 represent, respectively, minimum and maximum permissible concentrations in the control domain ΩC. Finally, (v, θ, u) are the solutions of the coupled state systems (4),(8) and (12). 3. Mathematical analysis of the state equations The mathematical analysis of the state equations (4),(8) and (12) have been previously studied by authors in [2]. We briefly recall the main results for convenience of the reader. For the modified Navier–Stokes system (4) we consider the following spaces: X1={v∈ [W1,3(Ω)]3: ∇ · v=0,v|∂Ω\(ΓC∪ΓT)=0}, ˜ X1={v∈ [W1,3(Ω)]3: ∇ · v=0,v|∂Ω=0}. (22) Then, associated to the previous spaces, we define: W1=W1,∞,2(0,T;X1,[L2(Ω)]3)∩C([0,T]; X1), ˜ W1=W1,∞,2(0,T;˜ X1,[L2(Ω)]3)∩C([0,T];˜ X1).(23) 6 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 Now, for the water temperature system (8), we consider the following spaces: X2= {θ∈H1(Ω):θ|ΓS∈L5(ΓS)}, ˜ X2= {θ∈X2:θ|ΓT=0}.(24) If we define the following norm associated to above space X2: ∥θ∥X2= ∥θ∥H1(Ω)+ ∥θ∥L5(ΓS), we have that X2is a reflexive separable Banach space (cf. Lemma 3.1 of [13]), and that ˜ X2⊂L2(Ω)⊂˜ X′ 2is an evolution triple. So, we consider: W2= {θ∈W1,2,5/4(0,T;X2,X′ 2):θ|ΓS∈L5(0,T;L5(ΓS))} ∩ L∞(0,T;L2(Ω)), ˜ W2= {θ∈W1,2,5/4(0,T;˜ X2,˜ X′ 2):θ|ΓS∈L5(0,T;L5(ΓS))} ∩ L∞(0,T;L2(Ω)).(25) Finally, for the eutrophication system (12), we define: X3= [H1(Ω)]5, ˜ X3= {u∈X3:u|ΓT=0},(26) and we consider the following spaces associated to them: W3=W1,2,2(0,T;X3,X′ 3), ˜ W3=W1,2,2(0,T;˜ X3,˜ X′ 3).(27) From this section we will assume the following hypotheses for coefficients and data in the analytical study of the problem: •gk∈H1(0,T), with gk(0) =0, ∀k=1,...,NCT , •v0∈[H2 σ(Ω)]3= {v∈ [H2(Ω)]3: ∇ · v=0,v|∂Ω=0} ⊂ ˜ X1, •θS∈L2(0,T;L2(ΓS)), •θN∈L2(0,T;L2(ΓN)), •Tr∈L5(0,T;L5(ΓS)), •θ0∈X2, •I0∈L∞(0,T), •u0∈X3. Under these hypotheses we will state now two lemmas (whose demonstrations can be found in [14,15], respectively), which will allow us to reformulate the state systems (4),(8) and (12) as homogeneous Dirichlet problems. Lemma 1. There exists a linear continuous extension: Rv: [H1(0,T)]NCT →W1,2,2(0,T; [H2 σ(Ω)]3,[H2 σ(Ω)]3) g→Rv(g)=ζg (28) such that ζg|∂Ω=φg, where φgis defined by (3), and H2 σ(Ω)= {u∈ [H2(Ω)]3: ∇ · u=0}.□ Remark 1. It is worthwhile emphasizing here that, thanks to the construction done in the proof of Lemma 1, we have ν∫Ω e(ζg):e(η)dx=0,∀η∈˜ X1, and, consequently, this term will not appear in the corresponding variational formulation. □ Lemma 2. We have that the following operator is compact Rh: [L2(0,T)]NCT →W1,2,2(0,T;H2(Ω),H2(Ω)) h→Rh(h)=ζh,(29) where: ζh(x,t)= NCT ∑ k=1 β0(ϕk(x)) ∫T −T ρϵ(t−ϵ−s)γk h(s)ds,(30) with γk h(s)∈L2(−T,T), for k =1,...,NCT , defined by: γk h(s)=⎧ ⎨ ⎩ 1 µ(Ck)∫Ck θ0dγif s ≤0, hk(s)if s >0, (31) 7 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 and β0:u∈H3/2(∂Ω)→β0(u)∈H2(Ω)the right inverse of the classical trace operator γ0, i.e., such that (γ0◦β0)(u)=u (cf. Theorem 8.3. of [8]). We also have the existence of a constant C, that depends continuously on the space–time computational domain and the initial temperature θ0, such that: ∥ζh∥W1,2,2(0,T;H2(Ω),H2(Ω)) ≤C(θ0) (1 + ∥h∥[L2(0,T)]NCT ).□(32) Now, we will establish the following notations, in order to consider the homogeneous Dirichlet systems. So, given elements (z, ξ, w)∈˜ W1ט W2ט W3, we define (v, θ, u)∈W1×W2×W3in the following way: •v=z+ζg, with ζg∈W1,2,2(0,T; [H2 σ(Ω)]3,[H2 σ(Ω)]3) the extension of control ggiven by Lemma 1. •θ=ξ+ζhθ, with ζhθ∈W1,2,2(0,T;H2(Ω),H2(Ω))the extension of hθobtained from Lemma 2, where: hk θ(s)=1 µ(Ck)∫Ck θ(s)dγ , k=1,...,NCT .(33) •ui=wi+ζhi u, with ζhi u∈W1,2,2(0,T;H2(Ω),H2(Ω)) the extension of hi uobtained from Lemma 2 with obvious modifications, where: hi,k u(s)=1 µ(Ck)∫Ck ui(s)dγ , k=1,...,NCT ,i=1,...,5.(34) Thus, using above notations, we can reformulate the state systems (4),(8) and (12) in the following way: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂z ∂t+ ∇(ζg+z)z+ ∇zζg −div (2νe(z)+2νtur ∫Ω[e(ζg+z):e(ζg+z)]1/2e(ζg+z)) +∇p=α0(θ−θ0)ag−∂ζg ∂t− ∇ζgζg+2ν∇ · e(ζg) in Ω×(0,T), z=0on ∂Ω×(0,T), z(0) =v0in Ω. (35) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ξ ∂t+v· ∇ξ− ∇ · (K∇ξ) = −∂ζhθ ∂t−v· ∇ζhθ+ ∇ · (K∇ζhθ) in Ω×(0,T), ξ=0 on Tk×(0,T),for k=1,...,NCT , K∂ξ ∂n= −K∂ζhθ ∂non Ck×(0,T),for k=1,...,NCT , K∂ξ ∂n=bN 1(θN−ζhθ−K bN 1 ∂ζhθ ∂n−ξ)on ΓN×(0,T), K∂ξ ∂n=bS 1(θS−ζhθ−K bS 1 ∂ζhθ ∂n−ξ) +bS 2(T4 r− |ξ+ζhθ|3(ξ+ζhθ))on ΓS×(0,T), ξ(0) =θ0−ζhθ(0) in Ω. (36) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂wi ∂t+v· ∇wi− ∇ · (µi∇wi)=Ai(x,t, θ, ζhu+w) −∂ζhi u ∂t−v· ∇ζhi u+ ∇ · (µi∇ζhi u) in Ω×(0,T), ∂wi ∂n= −µi∂ζhi u ∂non (ΓS∪ΓN∪ΓC)×(0,T), wi=0 on Tk×(0,T),for k=1,...,NCT , wi(0) =u0,i−ζhi u(0) in Ω,i=1,...,5. (37) It is worthwhile noting here that all three previous systems show homogeneous Dirichlet boundary conditions and, consequently, we will be able to define the concept of solution of the original state systems (4),(8) and (12) in terms of the modified state systems (35),(36) and (37). It should be also noted that, in the case of systems (8) and (12), the coupling terms in the Dirichlet boundary conditions are now transferred to the partial differential equations in systems (36) and (37). Definition 1 (The Concept of Solution for the State Systems).An element (v, θ, u)∈W1×W2×W3is a solution for the state systems (4),(8) and (12), if there exists an element (z, ξ, w)∈˜ W1ט W2ט W3such that: 8 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 •v=z+ζg, with ζg∈W1,2,2(0,T; [H2 σ(Ω)]3,[H2 σ(Ω)]3) as given by Lemma 1,z(0) =v0, a.e. x∈Ω, and z∈˜ W1the solution of the following variational formulation: ∫Ω ∂z ∂t·ηdx+∫Ω ∇(ζg+z)z·ηdx+∫Ω ∇zζg·ηdx+2ν∫Ω e(z):e(η)dx +2νtur ∫Ω[e(ζg+z):e(ζg+z)]1/2e(ζg+z):e(η)dx =∫Ω Hg·ηdx,a.e. t∈(0,T),∀η∈˜ X1, (38) where: Hg=α0(θ−θ0)ag−∂ζg ∂t− ∇ζgζg∈L2(0,T; [L2(Ω)]3).(39) •θ=ξ+ζhθ, with ζhθ∈W1,2,2(0,T;H2(Ω),H2(Ω)) obtained from Lemma 2 for hθ∈ [L2(0,T)]NCT defined by (33), ξ(0) =θ0−ζhθ(0), a.e. x∈Ω, and ξ∈˜ W2the solution of the following variational formulation: ∫Ω ∂ξ ∂tηdx+∫Ω v· ∇ξη dx+K∫Ω ∇ξ· ∇ηdx+bN 1∫ΓN ξη dγ +bS 1∫ΓS ξη dγ+bS 2∫ΓS |ξ+ζhθ|3(ξ+ζhθ)ηdγ=∫Ω Hhθηdx +∫ΓC gC hθηdγ+bN 1∫ΓN gN hθηdγ+bS 1∫ΓS gS hθηdγ +bS 2∫ΓS T4 rηdγ , a.e. t∈(0,T),∀η∈˜ X2, (40) where: Hhθ=∂ζhθ ∂t−v· ∇ζhθ+ ∇ · (K∇ζhθ)∈L2(0,T;L2(Ω)), gC hθ= −K∂ζhθ ∂n∈L2(0,T;L2(ΓC)), gN hθ=θN−ζhθ−K bN 1 ∂ζhθ ∂n∈L2(0,T;L2(ΓN)), gS hθ=θS−ζhθ−K bS 1 ∂ζhθ ∂n∈L2(0,T;L2(ΓS)). (41) •ui=wi+ζhi u, with ζhi u∈W1,2,2(0,T;H2(Ω),H2(Ω)) obtained from Lemma 2 for hi u∈ [L2(0,T)]NCT defined by (34), w(0) =u0−ζhu(0), a.e. x∈Ω, and w∈˜ W3the solution of the following variational formulation: ∫Ω ∂w ∂t·ηdx+∫Ω ∇wv ·ηdx+Λµ∫Ω ∇w: ∇ηdx =∫Ω A(θ, ζhu+w)·ηdx+∫Ω Hu·ηdx +∫ΓS∪ΓN∪ΓC gu·ηdγ , a.e. t∈(0,T),∀η∈˜ X3, (42) where Λµ=diag(µ1, . . . , µ5)∈M5×5(R) is a diagonal matrix with diffusion coefficients, and: Hi u= −∂ζhi u ∂t−v· ∇ζhi u+ ∇ · (µi∇ζhi u)∈L2(0,T;L2(Ω)), gi u= −µi∂ζhi u ∂n∈L2(0,T;L2(ΓS∪ΓN∪ΓC)),i=1,...,5. □(43) Remark 2. We have the following dependence scheme between the elements of state system: g→(v←→ θ)→u. Therefore, we can separate the mathematical analysis of systems (4)–(8) from system (12). The coupled system (4)– (8) has been fully analyzed by the authors in [14,15]. Thus, following the results there, we can assure that, for each control g∈ [H1(0,T)]NCT , there exists a solution (v, θ)∈˜ W1ט W2of the thermo-hydrodynamic system (4)–(8). We must remark here that, due to the complexity of this nonlinear system, we cannot obtain a uniqueness result for the thermo-hydrodynamic solution (v, θ) under our general hypotheses. However, this property will not be necessary in our approach, and previous existence result will be sufficient for our argumentation. So, we can focus now all our attention 9 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 Fig. 3. Dependence scheme for the discretized variables. with un+1∈Zhthe solution of the discretized eutrophication model. Thus, we can express: U= {g∈Uad :Gn(g)∈ [λm, λM],∀n=1,...,N}.(78) It is worthwhile remarking here that, due the type of time discretization considered for the material derivatives (64), the control gNacts over the species and the temperature at time tN+1. So, it will be necessary to compute one additional time step in the case of temperature and species in order to take into account this control. This fact can be more clearly noticed in the dependence scheme shown in Fig. 3. 5. Water velocity and pressure: Given v0∈Wh, the pair velocity/pressure (vn+1,pn+1)∈Vh×Mh, for each n=0,1,...,N−1, with: vn+1 |Tk= −gn+1,k µ(Tk)n,vn+1 |Ck=gn+1,k µ(Ck)n,∀k=1,...,NCT ,(79) is the solution of the fully discretized system: α∫Ω vn+1·ηdx+∫Ω β(vn)e(vn+1):e(η)dx−∫Ω pn+1∇ · ηdx −∫Ω ∇ · vn+1q dx−λ∫Ω pn+1q dx=α∫Ω (vn◦Xn −)·ηdx +∫Ω α0(θn−θ0)ag·ηdx,∀η∈Wh,∀q∈Mh, (80) where λ > 0 is the penalty parameter and β(vn)=2ν+2νtur [e(vn):e(vn)]1/2. 6. Water temperature: Given θ0∈Kh, the temperature θn+1∈Kh, for each n=0,...,N, with θn+1 |Tk=1 µ(Ck)∫Ck θndγ , ∀k=1,...,NCT ,(81) is the solution of the fully discretized system: α∫Ω θn+1ηdx+K∫Ω ∇θn+1· ∇ηdx+bN 1∫ΓN θn+1ηdγ +bS 1∫ΓS θn+1ηdγ=α∫Ω (θn◦Xn −)ηdx+bN 1∫ΓN θn+1 Nηdγ +bS 1∫ΓS θn+1 Sηdγ+bS 2∫ΓS (Tn r 4− |θn|3θn)ηdγ , ∀η∈Hh. (82) 7. Eutrophication species concentration: Given u0∈Xh, the species concentration un+1∈Xh, for each n=0,...,N, with: un+1 |Tk=1 µ(Ck)∫Ck undγ , ∀k=1,...,NCT ,(83) 16 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 is the solution of the fully discretized system: α∫Ω un+1·ηdx+∫Ω Λµ∇un+1: ∇ηdx+∫Ω An(θn,un)un+1·ηdx =α∫Ω(un◦Xn −)·ηdx,∀η∈Zh, (84) where An(θn,un)∈R5×5is the following matrix: ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0Cnc Ln(θn)un,1 KN+ |un,1|−Cnc Kr0 0 −Cnc KrdD(θn) 0Kr−Ln(θn)un,1 KN+ |un,1|+Kmf Kzun,2 KF+ |un,2|0 0 0 0 −CfzKzun,2 KF+ |un,2|+Kmz 0 0 0−Kmf −Kmz 0KrdD(θn) 0Coc Kr−Coc L(θn)un,1 KN+ |un,1|0 0 Coc KrdD(θn) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . 5.2. Numerical resolution of the optimization problem Once developed above space–time discretization, as introduced in previous section, we obtain the following discrete optimization problem: (P) min{J(g):g∈U} In order to solve this nonlinear optimization problem, we will use the interior point algorithm IPOPT [21] interfaced with the FreeFem++ code that we have developed. One of the requirements for using previous algorithm is the knowledge of functions that compute the gradient of the cost functional and the Jacobian matrix of the constraints. In the case of the cost functional, we have that its differential δgJ(g)∈L(RN×NCT ,R) is such that, for any δg= (δg1, . . . , δgN)∈RN×NCT : δgJ(g)(δg)=σ1 N ∑ n=1 NCT ∑ k=1 gn,kδgn,k +σ2 N−1 ∑ n=1 NCT ∑ k=1 (gn+1,n−gn,k)(δgn+1,n−δgn,k). (85) Therefore, [∇gJ(g)]i=δgJ(g)(ei), where ei,i=1,...,N×NCT , is the ith vector of the canonical basis in RN×NCT . In the case of the Jacobian matrix of the constraints, we know that the differential associated to the application G:Uad ⊂RN×NCT →RNis such that δgG(g)∈L(RN×NCT ,RN). So, given any element δg∈RN×NCT , we have that δgG(g)(δg)∈RN, and the Jacobian matrix JgG(g)∈MN×(NCT ×N)is such that [JgG(g)]j,i= ⟨δgG(g)(ei),˜ ej⟩, where ˜ ej, j=1,...,N, is the jth vector of the canonical basis in RN. As above commented, for computing previous matrix we can use either the linearized state systems or the adjoint state systems. The choice of one method or another depends on the relation between the dimension of the space of controls (N×NCT ) and the dimension of the space where the application Gtakes values (N). •When using the linearized systems, we would have to solve NCT ×Ntimes these systems (in this case, we would compute the Jacobian matrix column by column): JgG(g)=(δgG(g)(e1)δgG(g)(e2)· · · δgG(g)(eNCT −1)δgG(g)(eNCT )), where δgG(g)(ek)∈MN×1(R), for k=1,...,NCT . •When employing the adjoint state systems, we would have to solve Ntimes these systems (now, we would compute the Jacobian row by row): JgG(g)(δg)=⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∇gG1(g) ∇gG2(g) . . . ∇gGN−1(g) ∇gGN(g) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , where ∇gGi(g)∈RNCT is such that [∇gGj(g)]i=δgGj(g)(ei), j=1,...,N,i=1,...,N×NCT . 17 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 In our case, NCT >1. So, it is more advantageous to employ the adjoint state systems and computing the Jacobian matrix row by row (j=1,...,N). However, in order to obtain a computational expression for the Jacobian matrix using the adjoint systems it will be necessary deriving first a theoretical expression using the linearized systems and then applying a transposition procedure. Lemma 7 (Computing the Jacobian Matrix Using Linearized Systems).Within the framework introduced in this Section, we have the following expression for the Jacobian matrix of the constraints using the linearized equations: Given an element δg∈RN×NCT , then δgG(g)(δg)=⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 µ(ΩC)∫ΩC δu2,5dx . . . 1 µ(ΩC)∫ΩC δuN+1,5dx ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , where {(δvn, δpn)}N n=0⊂Vh×Mh,{δθn}N+1 n=0⊂Khand {δun}N+1 n=0⊂Xhare, respectively, the solutions of the linearized hydrodynamic model, the linearized thermic model and the linearized eutrophication model, defined as: •Linearized system for water velocity and pressure: Given (δv0, δp0)=(0,0), for each n =0,...,N−1,(δvn, δpn)∈ Vh×Mh, with δvn+1 |Tk= −δgn+1,k µ(Tk)n, δvn+1 |Ck=δgn+1,k µ(Ck)n,∀k=1,...,NCT ,(86) is the solution of: α∫Ω δvn+1·ηdx+∫Ω β(vn)e(δvn+1):e(η)dx−∫Ω δpn+1∇ · ηdx −∫Ω ∇ · δvn+1q dx−λ∫Ω δpn+1q dx=α∫Ω (δvn◦Xn −)·ηdx −∫Ω (∇vn◦Xn −)δvn·ηdx+∫Ω α0δθnag·ηdx −∫Ω γ(vn)e(vn):e(δvn)e(vn+1):e(η)dx,∀η∈Wh,∀q∈Mh, (87) where γ(vn)=2νtur [e(vn):e(vn)]−1/2. •Linearized system for water temperature: Given δθ0=0, for each n =0,...,N, δθn+1∈Kh, with: δθn+1 |Tk=1 µ(Ck)∫Ck δθndγ , ∀k=1,...,NCT ,(88) is the solution of: α∫Ω δθn+1ηdx+K∫Ω ∇δθn+1· ∇ηdx+bN 1∫ΓN δθn+1ηdγ +bS 1∫ΓS δθn+1ηdγ=α∫Ω (δθn◦Xn −)ηdx −∫Ω (∇θn◦Xn −)·δvnηdx−4bS 2∫ΓS |θn|3δθnηdγ , ∀η∈Hh. (89) •Linearized system for eutrophication model: Given δu0=0, for each n =0,...,N, δun+1∈Xh, with: δun+1 |Tk=1 µ(Ck)∫Ck δundγ , ∀k=1,...,NCT ,(90) is the solution of: α∫Ω δun+1·ηdx+∫Ω Λµ∇δun+1: ∇ηdx +∫Ω An(θn,un)δun+1·ηdx=α∫Ω (δun◦Xn −)·ηdx −∫Ω (∇un◦Xn −)δvn·ηdx−∫Ω δθAn(θn,un)(δθn)un+1·ηdx −∫Ω δuAn(θn,un)(δun)un+1·ηdx,∀η∈Zh. (91) 18 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 Proof. The proof is straightforward, where the only drawback is related to the computation of terms of the type δg(ϕ(g,x−∆tv(g,x)))(δg), where ϕ(g,x) and v(g,x) are vector functions smooth enough (the scalar case would be analogous). Nevertheless, using the chain rule, we can easily obtain that: δg(ϕ(g,x−∆tv(g,x)))(δg) =δgϕ(g,x−∆tv(g,x))(δg)−∆tδxϕ(g,x−∆tv(g,x))(δgv(g,x)(δg)) ≡(δϕ◦X−)−∆t(∇ϕ◦X−)δv. (we must note here that, in our specific formulation, we deal with the function b:x∈R→b(x)=x|x|3, that is differentiable in R, with b′(x)=4|x|3). Lemma 8 (Computing the Jacobian Matrix Using the Adjoint Equations).Within the framework introduced in this Section, we have the following expression for the Jacobian matrix of the constraints using the adjoint systems: For each row k =1,...,N, the matrices {∇gnGk(g)}k n=1⊂M1×NCT (R)can be computed using the following expressions: •If n ∈ {1,...,k}\{N}, δgnGk(g)(δgn)= NCT ∑ i=1 δgn,i µ(Ti)∫Ti β(vn−1)e(wn−1)n·n−qn−1dγ + NCT ∑ i=1 δgn,i µ(Ci)∫Ci qn−1−β(vn−1)e(wn−1)n·ndγ + NCT ∑ i=1 δgn,i µ(Ti)∫Ti γ(vn)e(vn+1):e(wn)e(vn)n·ndγ − NCT ∑ i=1 δgn,i µ(Ci)∫Ci γ(vn)e(vn+1):e(wn)e(vn)n·ndγ . •If n =N, δgnGk(g)(δgn)= NCT ∑ i=1 δgn,i µ(Ti)∫Ti β(vn−1)e(wn−1)n·n−qn−1dγ + NCT ∑ i=1 δgn,i µ(Ci)∫Ci qn−1−β(vn−1)e(wn−1)n·ndγ , where if we introduce, for each row k =1,...,N, the following vector (defined from the usual Kronecker delta δij and the indicator function of subset ΩC): Hn+1 k=(0,0,0,0,1 µ(ΩC)χΩCδkn)∈M1×5(R),n=0,...,N, then the adjoint states associated to the eutrophication system {zn}N+1 n=0⊂Zh, to the hydrodynamic system {(wn,qn)}N n=0⊂ Wh×Mh, and to the temperature system {ξn}N+1 n=0⊂Hhare, respectively, the solution of the following systems: •Adjoint system for eutrophication model: –For n =N+1,zn=0. –For n =N, zn∈Zhis such that: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ αzn− ∇ · (Λµ∇zn)+An(θn,un)Tzn=α(zn+1◦Xn+1 +)+Hn+1 kin Ω, zn=0on ΓT, Λµ∇znn=0on ∂Ω\(ΓT∪ΓC) Λµ∇znn= − 1 µ(Ck)∫Tk Λµ∇zn+1ndγ′on Ck,k=1,...,NCT , (92) where Xn+1 +(x)=x+∆tvn+1. 19 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 –For n =N−1,...,0,zn∈Zhis such that: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ αzn− ∇ · (Λµ∇zn)+An(θn,un)Tzn=α(zn+1◦Xn+1 +)+Hn+1 k − 5 ∑ l=1 [∇uAn+1 l(θn+1,un+1)]Tun+2zn+1,lin Ω, zn=0on ΓT, Λµ∇znn=0on ∂Ω\(ΓT∪ΓC), Λµ∇znn= − 1 µ(Ck)∫Tk Λµ∇zn+1ndγ′on Ck,k=1,...,NCT , (93) •Adjoint system for water temperature: –For n =N+1,ξn=0. –For n =N, ξn∈Hhis such that: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ αξn− ∇ · (K∇ξn)=α(ξn+1◦Xn+1 +)in Ω, ξn=0on ΓT, K∇ξn·n= −bS 1ξn−4bS 2|θn+1|3ξn+1on ΓS, K∇ξn·n= −bN 1ξnon ΓN, K∇ξn·n= − K µ(Ck)∫Tk ∇ξn+1·ndγ′on Ck,k=1,...,NCT . (94) –For n =N−1,...,0,ξn∈Hhis such that: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ αξn− ∇ · (K∇ξn)=α(ξn+1◦Xn+1 +) −d dθAn+1(θn+1,un+1)un+2·zn+1+α0ag·wn+1in Ω, ξn=0on ΓT, K∇ξn·n= −bS 1ξn−4bS 2|θn+1|3ξn+1on ΓS, K∇ξn·n= −bN 1ξnon ΓN, K∇ξn·n= − K µ(Ck)∫Tk ∇ξn+1·ndγ′on Ck,k=1,...,NCT . (95) •Adjoint system for water velocity and pressure: –For n =N, (wn,q0)=(0,0). –For n =N−1,(wn,qn)∈Wh×Mhis such that: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ αwn−div(β(vn)e(wn)) + ∇qn=α(wn+1◦Xn+1 +) −(∇vn+1◦Xn+1 −)Twn+1−(∇un+1◦Xn+1 −)Tzn+1 −(∇θn+1◦Xn+1 −)Tξn+1in Ω, ∇ · wn+λqn=0in Ω, wn=0on ∂Ω. (96) –For n =N−2,...,0,(wn,qn)∈Wh×Mhis such that: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ αwn−div(β(vn)e(wn)) + ∇qn =α(wn+1◦Xn+1 +)−(∇vn+1◦Xn+1 −)Twn+1 −(∇un+1◦Xn+1 −)Tzn+1−(∇θn+1◦Xn+1 −)Tξn+1 +∇ · (γ(vn+1)e(vn+2):e(wn+1)e(vn+1)) in Ω, ∇ · wn+λqn=0in Ω, wn=0on ∂Ω. (97) Remark 5. In order to simplify the proof of above Lemma, we have established the adjoint systems (96)–(97),(94)–(92) and (95)–(93) in a strong formulation (contrary to the case of the linearized systems (87),(89) and (91), where we have proposed a variational formulation). It is also clear that these adjoint systems easily admit a variational formulation, but we have chosen to formulate them in a strong form for a better understanding of the demonstration. □ Proof. Let us consider as a test functions in the linearized systems (87),(89) and (91), respectively, the nth component of the sequences {(wn,qn)}N n=0⊂Wh×Mh,{ξn}N+1 n=0⊂Hhand {zn}N n=0⊂Zh, such that wN=0,qN=0, ξN+1=0 and zN+1=0, and let us sum in nfrom 0 to N. Then, after some straightforward computations, taking into account the final conditions for the adjoint systems and the initial conditions for the linearized ones, we have: 20 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 •For eutrophication model: N ∑ n=0[α∫Ω zn·δun+1dx−∫Ω ∇ · (Λµ∇zn)·δun+1dx +∫Ω An(θn,un)Tzn·δun+1dx]= N ∑ n=0[α∫Ω (zn+1◦Xn+1 +)·δun+1dx −∫Ω (∇un+1◦Xn+1 −)Tzn+1·δvn+1dx − NCT ∑ k=1∫Ck(Λµ∇znn+1 µ(Ck)∫Tk Λµ∇zn+1ndγ′)·δun+1dγ] − N−1 ∑ n=0[∫Ω(d dθAn+1(θn+1,un+1)un+2·zn+1)δθn+1dx +∫Ω(5 ∑ l=1 [∇uAn+1 l(θn+1,un+1)]Tun+2zn+1,l)·δun+1dx], (98) with Xn+1 +(x)=x+∆tvn+1, and where we are assuming δvN+1=0in order to simplify the notation. •For water temperature: N ∑ n=0[α∫Ω ξnδθn+1dx−∫Ω ∇ · (K∇ξn)δθn+1dx +∫ΓN(bN 1ξn+K∇ξn·n)δθn+1dγ +∫ΓS(bS 1ξn+4bS 2|θn+1|3ξn+1+K∇ξn·n)δθn+1dγ] = N ∑ n=0[α∫Ω (ξn+1◦Xn+1 +)δθn+1dx −∫Ω (∇θn+1◦Xn+1 −)Tξn+1·δvn+1dx − NCT ∑ k=1∫Ck(K∇ξn·n+1 µ(Ck)∫Tk K∇ξn+1·ndγ′)δθn+1dγ], (99) where, for the sake of simplicity, we have also assumed δvN+1=0. •For water velocity: N−1 ∑ n=0[α∫Ω wn·δvn+1dx−∫Ω div (β(vn)e(wn))·δvn+1dx −∫Ω ∇ · wnδpn+1dx+∫Ω ∇qn·δvn+1dx−λ∫Ω qnδpn+1dx] = N ∑ n=0[∫Ω α0agwn+1δθn+1dx]+ N−1 ∑ n=0[α∫Ω (wn+1◦Xn+1 +)·δvn+1dx −∫Ω (∇vn+1◦Xn+1 −)Twn+1·δvn+1dx] + N−2 ∑ n=0[∫Ω div(γ(vn+1)e(vn+2):e(wn+1)e(vn+1)) ·δvn+1dx] + N−1 ∑ n=0 NCT ∑ k=1 δgn+1,k[1 µ(Tk)∫Tk(β(vn)e(wn)n·n−qn)dγ −1 µ(Ck)∫Ck(β(vn)e(wn)n·n−qn)dγ]21 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 Fig. 4. Triangular mesh of the domain Ωfor the numerical tests. + N−2 ∑ n=0 NCT ∑ k=1 δgn+1,k[1 µ(Tk)∫Tk γ(vn+1)e(vn+2):e(wn+1)e(vn+1)n·ndγ −1 µ(Ck)∫Ck γ(vn+1)e(vn+2):e(wn+1)e(vn+1)n·ndγ], (100) where we have assumed wN+1=0. Thus, if we define {(wn,qn)}N n=0⊂Wh×Mh,{ξn}N+1 n=0⊂Hhand {zn}N+1 n=0⊂Zh, such that wN=0,qN=0, ξN+1=0 and zN+1=0, as the solutions of the adjoint system (96)–(97),(94)–(92) and (95)–(93), respectively, we obtain, after summing above expressions (98),(99) and (100), that: N ∑ n=0∫Ω Hn+1 k·δun+1dx= N−1 ∑ n=0 NCT ∑ k=1 δgn+1,k[1 µ(Tk)∫Tk(β(vn)e(wn)n·n−qn)dγ −1 µ(Ck)∫Ck(β(vn)e(wn)n·n−qn)dγ] + N−2 ∑ n=0 NCT ∑ k=1 δgn+1,k[1 µ(Tk)∫Tk γ(vn+1)e(vn+2):e(wn+1)e(vn+1)n·ndγ −1 µ(Ck)∫Ck γ(vn+1)e(vn+2):e(wn+1)e(vn+1)n·ndγ]. And, finally, from the definition: N ∑ n=0∫Ω Hn+1 k·δun+1dx=1 µ(ΩC)∫ΩC δuk+1,5dx.(101) 5.3. Numerical results In order to simplify the graphical representation of the computational results for the numerical tests developed in this study, we will present here only the case of a two dimensional domain Ω. So, we consider a space configuration similar to that presented in Fig. 1, with NCT =4 collector/injector pairs, in a rectangular domain of 20 m×16 m. We suppose that the diameter of each collector is 1 m and the diameter of each injector is 2 m. For the coefficients of the eutrophication model (12), we have used the same values as those appearing in [12], and for the thermo-hydrodynamic system (4),(8) we have employed the same values as in [15]. For the space discretization we have generated a regular mesh of 2989 vertices, as shown in Fig. 4. 22 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 Fig. 5. Standard profile for radiation temperature Tr. Fig. 6. Comparison of the mean concentrations of dissolved oxygen in ΩC, taking a time step length ∆t=450 s, for a constant flow rate of 1.0×10−4m3s−1in all the pumps, and for the case without pumping. The control domain ΩCcorresponds to a 3 m strip at the bottom of the domain, and all the numerical tests have been performed in a temporal horizon of 12 h (T=43200 s). Finally, in order to simulate the effects of solar radiation for the heat Eq. (8), we consider the standard function Trdepicted in Fig. 5. We must remark that our main goal in this first approximation to the numerical resolution of the problem is trying to understand if we can improve the management of the pumps with respect to a constant operating regime. So, given a constant reference control ˜ g, with ˜ gn,k=C(constant), for n=1,...,N,k=1,...,NCT , we will solve the following modification of the original optimization problem (P): (ˆ P) min{J(g):g∈Uad,G(g)≥G(˜ g)}. In other words, we want to find an optimal control ˆ g∈Uad that supplies us with a higher concentration of dissolved oxygen than that obtained with the constant control ˜ g, and that minimizes the energy cost functional J. As an illustration to this behavior, in Fig. 6 we can see the evolution of the mean concentration of dissolved oxygen in the control domain ΩC considering a constant reference control˜ gn,k=1.0×10−4m3s−1,∀n=1,...,N,∀k=1,...,NCT , compared to the mean concentration assuming that all the pumps are out of service (that is, ˜ gn,k=0.0 m3s−1,∀n=1,...,N,∀k=1,...,NCT ). We observe how, if the pumps are out of service, the mean concentration of dissolved oxygen in the control domain decays gradually but, nevertheless, if we consider a constant flow rate (not necessarily large), this mean concentration of dissolved oxygen increases in a significant way. In this final part of the Section we present several numerical results that we have obtained using different choices of the time step length ∆t. We must mention that in the numerous numerical tests developed, we have always obtained that G(ˆ g)=G(˜ g), and also a reduction in the value of the cost functional J(ˆ g)<J(˜ g). So, in Fig. 7 we can see the optimal control that we have obtained taking σ1=0.5 and σ2=1−σ1=0.5, for time steps of ∆t=3600 s and ∆t=1800 s (corresponding to N=12 and N=24, respectively). In Fig. 8 we can find the optimal control corresponding to time steps of ∆t=900 s and ∆t=450 s (N=48 and N=96, respectively), showing the robustness of our methodology. 23 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 Fig. 7. Evolution of the optimal flow rates for the four pumps taking weights σ1=σ2=0.5, and ∆t=3600 s (left) or ∆t=1800 s (right). Fig. 8. Evolution of the optimal flow rates for the four pumps taking σ1=σ2=0.5, and ∆t=900 s (left) or ∆t=450 s (right). Table 1 Functional cost evaluated in the reference control (˜ g) vs. optimal control (ˆ g). ∆t=3600 s ∆t=1800 s ∆t=900 s ∆t=450 s J(˜ g) 1.2000 ×10−72.4000 ×10−74.8000 ×10−79.6000 ×10−7 J(ˆ g) 1.0973 ×10−72.1865 ×10−74.3195 ×10−78.7104 ×10−7 We observe that the flow rates associated to the two upper collectors (g1and g3) are significantly higher than the corresponding to lower collectors (g2and g4). This is caused by the fact that the photosynthesis is more intense in the superficial layers and, consequently, the presence of dissolved oxygen is higher there. In Table 1 we can see the comparison between the functional cost evaluated in the reference control and in the optimal control. We can observe that, as we decrease the time step, the difference between the reference cost and the optimal cost increases. As it seems obvious, this is a straightforward consequence of the fact that, as we decrease the time step, we can act in a more precise way over the system and achieve better results. In Fig. 9 we can see the evolution of the constraints for the choice of the time step length ∆t=450 s. We can verify there that the optimal constraint G(˜ g) and the reference constraint G(ˆ g) are virtually indistinguishable, that is, with optimal strategy ˜ gwe obtain the same water quality in the control region as with the constant reference flow rate ˜ g, but with a significative decrease in energy cost. Finally, in Fig. 10 we show the concentration of dissolved oxygen in the whole domain Ωassociated to the optimal control solution for ∆t=450 s (left), and the concentration of dissolved oxygen when all the pumps are off (right), both in the last time step (corresponding to N=96). We can easily notice here the pumping effects associated to the optimal control in the bottom layer, with an evident improvement of water quality in the region. 24 F.J. Fernández, A. Martínez and L.J. Alvarez-Vázquez Journal of Computational and Applied Mathematics 421 (2023) 114886 Fig. 9. Evolution of the constraints, for ∆t=450 s, in the controlled and uncontrolled cases. Fig. 10. Concentration of dissolved oxygen in the last time step corresponding to the optimal solution (left), and without control (right). Data availability No data was used for the research described in the article. Acknowledgments The authors would like to thank the anonymous referees for their interesting comments, suggestions and corrections, as they have greatly contributed to improve the quality of the manuscript. The authors were partially supported by Xunta de Galicia, Spain, project ED431C 2019/02, and by the Agencia Estatal de Investigación (AEI) of Spain under grant MTM2016-75140-P, co-financed by the European Community fund FEDER. References [1] L.J. Alvarez-Vázquez, F.J. Fernández, A. Martínez, Optimal control of eutrophication processes in a moving domain, J. Franklin Inst. 351 (2014) 4142–4182. [2] F.J. Fernández, A. Martínez, L.J. Alvarez-Vázquez, A nonlinear problem related to artificial circulation in a lake, in: Nonlinear Analysis and Boundary Value Problems. NABVP 2018, Vol. 292, Springer Proceedings in Mathematics & Statistics, 2019. [3] A. Martínez, L.J. Alvarez-Vázquez, F.J. Fernández, Water artificial circulation for eutrophication control, Math. Control Relat. Fields 8 (2018) 277–313. [4] T.-S. Yeh, Bifurcation curves of positive steady-state solutions for a reaction–diffusion problem of lake eutrophication, J. Math. Anal. Appl. 449 (2017) 1708–1724. [5] B. Wang, Q. Qi, Modeling the lake eutrophication stochastic ecosystem and the research of its stability, Math. Biosci. 300 (2018) 102–114. [6] J. Jiao, P. Li, D. Feng, Dynamics of water eutrophication model with control, Adv. Diff. Equ. Pap. (383) (2018). [7] T. Roubíček, Nonlinear Partial Differential Equations with Application, Birkhäuser, Basel, 2013. 25