An optimization problem related to water artificial recirculation for controlling eutrophication

TL;DR

This paper introduces an optimal control framework for water recirculation aimed at preventing eutrophication in large waterbodies, combining nonlinear hydrodynamics and ecological models with numerical solutions.

Contribution

It presents a novel coupled PDE model for water recirculation and eutrophication control, along with existence analysis and an efficient numerical algorithm.

Findings

Numerical results demonstrate the effectiveness of the proposed control strategy.

The model successfully captures complex interactions in water quality dynamics.

The approach provides a practical tool for environmental management.

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.

Figures

Tables

Table 1. Table 1: Functional cost evaluated in the Reference Control ( 𝐠 ~ ~ 𝐠 \widetilde{\mathbf{g}} ) vs. Optimal Control ( 𝐠 ^ ^ 𝐠 \widehat{\mathbf{g}} ).

	$\Delta t=3600$ s	$\Delta t=1800$ s	$\Delta t=900$ s	$\Delta t=450$ s
$J\left(\widetilde{𝐠}\right)$	$1.2000e-07$	$2.4000e-07$	$4.8000e-07$	$9.6000e-07$
$J\left(\widehat{𝐠}\right)$	$1.0973e-07$	$2.1865e-07$	$4.3195e-07$	$8.7104e-07$

Equations

$$\mathbf{v}={\mathbf{\varphi }}_{\mathbf{g}}\text{mbox}on\partial \Omega \times (0,T),$$

$${\mathbf{\varphi }}_{\mathbf{g}}(\mathbf{x},t)=\sum _{k=1}^{N_{CT}}{g}^k(t)\left[\frac{{\phi }^k(\mathbf{x})}{\mu (T^k)}-\frac{{\tilde{\phi }}^k(\mathbf{x})}{\mu (C^k)}\right]\mathbf{n}(\mathbf{x})$$

$${\int }_{\partial \Omega }{\mathbf{\varphi }}_{\mathbf{g}}\cdot \mathbf{n}d\gamma =0.$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial\mathbf{v}}{\partial t}+\nabla\mathbf{v}\,\mathbf{v}-\text{div}(\Xi(\mathbf{v}))+\nabla p=\alpha^{0}(\theta-\theta^{0})\,\mathbf{a}_{g}\quad\mbox{in}\;\Omega\times(0,T),\\ \displaystyle\nabla\cdot\mathbf{v}=0\quad\mbox{in}\;\Omega\times(0,T),\\ \displaystyle\mathbf{v}=\boldsymbol{\phi}_{\mathbf{g}}\quad\mbox{on}\;\partial\Omega\times(0,T),\\ \displaystyle\mathbf{v}(0)=\mathbf{v}^{0}\quad\mbox{in}\;\Omega,\end{array}\right.$$

$$\Xi (\mathbf{v})={\frac{\partial D(e)}{\partial e}|}_{e=e(\mathbf{v})},\text{mbox}withe(\mathbf{v})=\frac{1}{2}\left(\nabla \mathbf{v}+\nabla {\mathbf{v}}^t\right),$$

$$D(e)=\nu \left[e:e\right]+\frac{2}{3}{\nu }_{tur}{\left[e:e\right]}^{3/2},$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial\theta}{\partial t}+\mathbf{v}\cdot\nabla\theta-\nabla\cdot(K\nabla\theta)=0\quad\mbox{in}\;\Omega\times(0,T),\\ \displaystyle\theta=\phi_{\theta}\quad\mbox{on}\;\Gamma_{T}\times(0,T),\\ \displaystyle K\frac{\partial\theta}{\partial\mathbf{n}}=0\quad\mbox{on}\;\Gamma_{C}\times(0,T),\vspace{0.1cm}\\ \displaystyle K\frac{\partial\theta}{\partial\mathbf{n}}=b_{1}^{N}(\theta_{N}-\theta)\quad\mbox{on}\;\Gamma_{N}\times(0,T),\vspace{0.1cm}\\ \displaystyle K\frac{\partial\theta}{\partial\mathbf{n}}=b_{1}^{S}(\theta_{S}-\theta)+b_{2}^{S}(T_{r}^{4}-|\theta|^{3}\theta)\quad\mbox{on}\;\Gamma_{S}\times(0,T),\\ \theta(0)=\theta^{0}\quad\mbox{in}\;\Omega,\end{array}\right.$$

$${\varphi }_{\theta }(\mathbf{x},t)=\sum _{k=1}^{N_{CT}}{\phi }^k(\mathbf{x}){\int }_{-T}^T{\rho }_{ϵ}(t-ϵ-s){\gamma }_{\theta }^k(s)ds$$

$$\gamma_{\theta}^{k}(s)=\left\{\begin{array}[]{lcl}\displaystyle\frac{1}{\mu(C^{k})}\int_{C^{k}}\theta^{0}\,d\gamma&\mbox{ if }&s\leq 0,\\ \displaystyle\frac{1}{\mu(C^{k})}\int_{C^{k}}\theta(s)\,d\gamma&\mbox{ if }&s>0,\displaystyle\end{array}\right.$$

$$\rho_{\epsilon}(t)=\left\{\begin{array}[]{lcl}\displaystyle\frac{c}{\epsilon}\,\exp\left(\frac{t^{2}}{t^{2}-\epsilon^{2}}\right)&\mbox{ if }&|t|<\epsilon,\\ \displaystyle 0&\mbox{ if }&|t|\geq\epsilon,\end{array}\right.$$

$${\int }_{\mathbb{R}}{\rho }_1(t)dt=1.$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial u^{i}}{\partial t}+\mathbf{v}\cdot\nabla u^{i}-\nabla\cdot(\mu^{i}\nabla u^{i})=A^{i}(\mathbf{x},t,\theta,\mathbf{u})\quad\mbox{in}\;\Omega\times(0,T),\\ \displaystyle u^{i}=\phi_{u^{i}}\quad\mbox{on}\;\Gamma_{T}\times(0,T),\\ \displaystyle\mu^{i}\frac{\partial u^{i}}{\partial\mathbf{n}}=0\quad\mbox{on}\;(\Gamma_{S}\cup\Gamma_{N}\cup\Gamma_{C})\times(0,T),\\ \displaystyle u^{i}(0)=u^{0,i}\quad\mbox{in}\;\Omega,\qquad i=1,\ldots,5,\end{array}\right.$$

$${\varphi }_{u^i}(\mathbf{x},t)=\sum _{k=1}^{N_{CT}}{\phi }^k(\mathbf{x}){\int }_{-T}^T{\rho }_{ϵ}(t-ϵ-s){\gamma }_{u^i}^k(s)ds,$$

$$\gamma_{u^{i}}^{k}(s)=\left\{\begin{array}[]{lcl}\displaystyle\frac{1}{\mu(C^{k})}\int_{C^{k}}u^{0,i}\,d\gamma,&\mbox{if}&s\leq 0,\vspace{0.1cm}\\ \displaystyle\frac{1}{\mu(C^{k})}\int_{C^{k}}u^{i}(s)\,d\gamma,&\mbox{if}&s>0.\displaystyle\end{array}\right.$$

$$\mathbf{A}(\mathbf{x},t,\theta,\mathbf{u})=\left[\begin{array}[]{l}\displaystyle-\frac{C_{nc}L(\mathbf{x},t,\theta)}{K_{N}+|u^{1}|}u^{1}u^{2}+C_{nc}K_{r}u^{2}+C_{nc}K_{rd}D(\theta)u^{4}\vspace{0.1cm}\\ \displaystyle\frac{L(\mathbf{x},t,\theta)}{K_{N}+|u^{1}|}u^{1}u^{2}-K_{r}u^{2}-K_{mf}u^{2}-\frac{K_{z}}{K_{F}+|u^{2}|}u^{2}u^{3}\vspace{0.1cm}\\ \displaystyle\frac{C_{fz}K_{z}}{K_{F}+|u^{2}|}u^{2}u^{3}-K_{mz}u^{3}\vspace{0.1cm}\\ \displaystyle K_{mf}u^{2}+K_{mz}u^{3}-K_{rd}D(\theta)u^{4}\vspace{0.1cm}\\ \displaystyle\frac{C_{oc}L(\mathbf{x},t,\theta)}{K_{N}+|u^{1}|}u^{1}u^{2}-C_{oc}K_{r}u^{2}-C_{oc}K_{rd}D(\theta)u^{4}\end{array}\right]$$

$$D(\theta )={\Theta }^{\theta -{\theta }^0},$$

$$D(\theta )=1+\log (\Theta )(\theta -{\theta }^0)$$

$$L(\mathbf{x},t,\theta )=\mu C_t^{\theta -{\theta }^0}\frac{I^0(t)}{I_s}{e}^{-{\phi }_1{x}_3},$$

$$L(\mathbf{x},t,\theta )=\mu \left(1+\log (C_t)(\theta -{\theta }^0)\right)\frac{I^0(t)}{I_s}{e}^{-{\phi }_1{x}_3}$$

$$(\mathcal{P})\min \{J(\mathbf{g}):\mathbf{g}\in {\mathcal{U}}_{ad},\frac{1}{\mu ({\Omega }_C)}{\int }_{{\Omega }_C}{u}^5(t)d\mathbf{x}\in [{\lambda }^m,{\lambda }^M]\},$$

$$\begin{array}[]{l}\mathcal{U}_{ad}=\{\mathbf{g}\in[H^{1}(0,T)]^{N_{CT}}\,:\,\mathbf{g}(0)=\boldsymbol{0},\,\text{and}\\ \displaystyle\qquad\|g^{k}\|_{H^{1}(0,T)}\leq c,\,\forall k=1,\ldots,N_{CT}\}\end{array}$$

$$J(\mathbf{g})=\frac{1}{2}\sum _{k=1}^{N_{CT}}{\int }_0^T{g}^k(t{)}^{2}dt+\frac{1}{2}\sum _{k=1}^{N_{CT}}{\int }_0^T\frac{d{g}^k}{dt}(t{)}^{2}dt,$$

$$W^{1,p,q}(0,T;V_1,V_2)=\left\{u\in L^p(0,T;V_1):\frac{du}{dt}\in L^q(0,T;V_2)\right\},$$

$$\begin{array}[]{rcl}\displaystyle\mathbf{X}_{1}&=&\left\{\mathbf{v}\in[W^{1,3}(\Omega)]^{3}\,:\,\nabla\cdot\mathbf{v}=0,\,\mathbf{v}_{|_{\partial\Omega\setminus(\Gamma_{C}\cup\Gamma_{T})}}=\boldsymbol{0}\right\},\\ \displaystyle\widetilde{\mathbf{X}}_{1}&=&\left\{\mathbf{v}\in[W^{1,3}(\Omega)]^{3}\,:\,\nabla\cdot\mathbf{v}=0,\,\mathbf{v}_{|_{\partial\Omega}}=\boldsymbol{0}\right\}.\end{array}$$

$$\begin{array}[]{rcl}\displaystyle\mathbf{W}_{1}&=&\displaystyle W^{1,\infty,2}(0,T;\mathbf{X}_{1},[L^{2}(\Omega)]^{3})\cap\mathcal{C}([0,T];\mathbf{X}_{1}),\\ \displaystyle\widetilde{\mathbf{W}}_{1}&=&\displaystyle W^{1,\infty,2}(0,T;\widetilde{\mathbf{X}}_{1},[L^{2}(\Omega)]^{3})\cap\mathcal{C}([0,T];\widetilde{\mathbf{X}}_{1}).\end{array}$$

$$\begin{array}[]{rcl}X_{2}&=&\{\theta\in H^{1}(\Omega)\,:\,\theta_{|_{\Gamma_{S}}}\in L^{5}(\Gamma_{S})\},\\ \widetilde{X}_{2}&=&\{\theta\in X_{2}\,:\,\theta_{|_{\Gamma_{T}}}=0\}.\end{array}$$

$$\Vert \theta {\Vert }_{X_2}=\Vert \theta {\Vert }_{H^1(\Omega )}+\Vert \theta {\Vert }_{L^5({\Gamma }_S)},$$

$$\begin{array}[]{rcl}\displaystyle W_{2}&=&\{\theta\in W^{1,2,5/4}(0,T;X_{2},X_{2}^{\prime}):\theta_{|_{\Gamma_{S}}}\in L^{5}(0,T;L^{5}(\Gamma_{S}))\}\cap L^{\infty}(0,T;L^{2}(\Omega)),\\ \widetilde{W}_{2}&=&\{\theta\in W^{1,2,5/4}(0,T;\widetilde{X}_{2},\widetilde{X}_{2}^{\prime}):\theta_{|_{\Gamma_{S}}}\in L^{5}(0,T;L^{5}(\Gamma_{S}))\}\cap L^{\infty}(0,T;L^{2}(\Omega)).\end{array}$$

$$\begin{array}[]{rcl}\displaystyle\mathbf{X}_{3}&=&[H^{1}(\Omega)]^{5},\\ \displaystyle\widetilde{\mathbf{X}}_{3}&=&\{\mathbf{u}\in\mathbf{X}_{3}\,:\,\mathbf{u}_{|_{\Gamma_{T}}}=\mathbf{0}\},\end{array}$$

$$\begin{array}[]{rcl}\displaystyle\mathbf{W}_{3}&=&W^{1,2,2}(0,T;\mathbf{X}_{3},\mathbf{X}_{3}^{\prime}),\\ \displaystyle\widetilde{\mathbf{W}}_{3}&=&W^{1,2,2}(0,T;\widetilde{\mathbf{X}}_{3},\widetilde{\mathbf{X}}_{3}^{\prime}).\end{array}$$

Full text

††journal: JOTA

1 Introduction: The environmental problem

Eutrophication is one of the most important problems of large masses of water (estuaries, lakes, reservoirs, etc.) and it is caused by undue high levels of nutrients (usually nitrogen and phosphorus) reaching the water. These nutrients mainly come from human activities (resulting in the discharge of sewage, detergents, fertilizers and so on, very rich in phosphate or nitrate), and can cause an excessive phytoplankton growth that lead to undesirable effects like algal blooms. This abnormal growth of algae directly affects the concentration of dissolved oxygen, mainly in the deeper layers, since the processes of remineralization of organic detritus (accumulated in the bottom due to the effects of sedimentation) consumes oxygen, which can lead to oxygen depletion of the body of water intro1 . In Figure 1 (left) we can find a schematic representation of the problem and its consequences.

Artificial circulation is a management technique for oxygenating eutrophic water bodies subject to quality problems, such as loss of oxygen, sediment accumulation and algal blooms. It disrupts stratification and minimizes the development of stagnant zones that may be subject to above commented water quality problems. In our particular case we are only interested in increasing the dissolved oxygen concentration in the bottom layers (but our methodology could be extended in a straightforward way to any phenomenon and any region). In the process of artificial recirculation, a set of flow pumps takes water from the well aerated upper layers by means of a collector and injects it into the poorly oxygenated bottom layers, through a pipeline, setting up a circulation pattern that prevents stratification. Then, oxygen-poor water from the bottom is circulated to the surface, where oxygenation from the atmosphere and photosynthesis can naturally occur intro2 . In Figure 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 intro3 ; intro4 ; intro5 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 a recent paper of the authors intro2 , where a simplified preliminary formulation of the problem is posed and briefly analyzed). Thus, 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, in the central part of the paper, we analyze the wellposedness of the corresponding state system, in order to 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.

2.1 The physical domain

We consider a domain $\Omega\subset\mathbb{R}^{3}$ corresponding, for instance, to a reservoir. In order to promote the artificial circulation of water inside the domain $\Omega$, we suppose the existence of a set of $N_{CT}$ pairs collector-injector $\{(C^{k},T^{k})\}_{k=1}^{N_{CT}}\subset\partial\Omega$ 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 $\partial\Omega$, such that it can be split into four disjoint subsets $\partial\Omega=\Gamma_{S}\cup\Gamma_{C}\cup\Gamma_{T}\cup\Gamma_{N}$, where $\displaystyle\Gamma_{C}$ corresponds to the part of the boundary where the water collectors are located ($\Gamma_{C}=\cup_{k=1}^{N_{CT}}C^{k}$), $\displaystyle\Gamma_{T}$ corresponds to the part of the boundary where the water injectors are located ($\Gamma_{T}=\cup_{k=1}^{N_{CT}}T^{k}$), $\Gamma_{S}$ is the top part of the boundary in contact with air, and $\displaystyle\Gamma_{N}=\partial\Omega\setminus\left(\Gamma_{S}\cup\Gamma_{C}\cup\Gamma_{T}\right)$ corresponds to the rest of the boundary. In particular, we suppose the boundary $\partial\Omega$ regular enough to assure the existence of elements $\varphi^{k},\,\widetilde{\varphi}^{k}\in H^{3/2}(\partial\Omega)$, for $k=1,\ldots,N_{CT}$, satisfying the following assumptions (mainly corresponding to suitable regularizations of the indicator functions of $T^{k}$ and $C^{k}$, respectively):

• •

$\varphi^{k}(\mathbf{x}),\,\widetilde{\varphi}^{k}(\mathbf{x})\geq 0$, a.e. $\mathbf{x}\in\partial\Omega$,

• •

$\varphi^{k}(\mathbf{x})=0$, a.e. $\mathbf{x}\in\partial\Omega\setminus T^{k}$, and $\int_{T^{k}}\varphi^{k}\,d\gamma=\mu(T^{k})$,

• •

$\widetilde{\varphi}^{k}(\mathbf{x})=0$, a.e. $\mathbf{x}\in\partial\Omega\setminus C^{k}$, and $\int_{C^{k}}\widetilde{\varphi}^{k}\,d\gamma=\mu(C^{k})$,

where $\mu(S)$ represents the $n-1$ dimensional measure of a generic set $S$ and $\beta_{0}:u\in H^{3/2}(\partial\Omega)\rightarrow\beta_{0}(u)\in H^{2}(\Omega)$ denotes the right inverse of the classical trace operator $\gamma_{0}$, i.e., $(\gamma_{0}\circ\beta_{0})(u)=u$ (cf. Theorem 8.3. in Chapter 1 of magenes1 ). Finally, we also consider a subdomain $\Omega_{C}\subset\Omega$, corresponding to the part of the domain where we want to increase the dissolved oxygen concentration (denoted as control domain in Figure 2).

2.2 The control variable

As above commented, our control will be the volumetric flow rate $({\rm m}^{3}\,{\rm s}^{-1})$ by pump $k$ at each time $t$, $g^{k}(t)\in H^{1}(0,T)$, for $k=1,\ldots,N_{CT}$, where $T$ ($\rm{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:

[TABLE]

where $\mathbf{v}(\mathbf{x},t)$ will denote the water velocity, and where:

[TABLE]

represents the given Dirichlet condition for the hydrodynamic system. It is immediate that, thanks to the regularity of the control $\mathbf{g}$ and of the functions $\{(\varphi^{k},\widetilde{\varphi}^{k})\}_{k=1}^{N_{CT}}$, we have that $\boldsymbol{\phi}_{\mathbf{g}}\in W^{1,2,2}(0,T;H^{3/2}(\partial\Omega),H^{3/2}(\partial\Omega))$ (cf. expression (20) below for a detailed definition of this Sobolev-Bochner space), and also that

[TABLE]

2.3 The thermo-hydrodynamic model

We denote by $\mathbf{v}(\mathbf{x},t)\ ({\rm m}\,{\rm s}^{-1})$ the solution of the following modified Navier-Stokes system with a Smagorinsky model of turbulence:

[TABLE]

where $\mathbf{a}_{g}$ $({\rm m}\,{\rm s}^{-2})$ is the gravity acceleration, $\alpha^{0}=-\frac{1}{\rho}\frac{\partial\rho}{\partial\theta}$ $({\rm K}^{-1})$ is the thermic expansion coefficient, $\rho$ is the density ,$\mathbf{v}^{0}$ is the initial velocity, and the boundary field $\boldsymbol{\phi}_{\mathbf{g}}$ is the element given by (2). The diffusion term $\Xi(\mathbf{v})$ is given by:

[TABLE]

where $D$ is a potential function (for instance, in the standard case of the classical Navier-Stokes equations, $D(e)=\nu\left[e:e\right]$, with $\nu$ $({\rm m}^{2}\,{\rm s}^{-1})$ the kinematic viscosity of the water, and, consequently, $\Xi(\mathbf{v})=2\nu\,e(\mathbf{v})$). However, in our case, the Smagorinsky model, the potential function is defined as in lady1 :

[TABLE]

where $\nu_{tur}$ $({\rm m}^{2})$ is the turbulent viscosity.

Regarding thermic effects, water temperature $\theta(\mathbf{x},t)$ ($\rm{K}$) is the solution of the following convection-diffusion partial differential equation with nonhomogeneous, nonlinear, mixed boundary conditions:

[TABLE]

where Dirichlet boundary condition $\phi_{\theta}$ is given by expression:

[TABLE]

with, for each $k=1,\ldots,N_{CT}$,

[TABLE]

representing the mean temperature of water in the collector $C_{k}$, and with the weight function $\rho_{\epsilon}$ defined by:

[TABLE]

for $c\in\mathbb{R}$ the positive constant satisfying the unitary condition:

[TABLE]

In other words, we are assuming that the mean temperature of water at each injector $T_{k}$ is a weighted average in time of the mean temperatures of water at its corresponding collector $C_{k}$. 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\epsilon,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\epsilon,t)$. Parameter $0<\epsilon<T$ represents, 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).

Moreover, for the other terms appearing in the formulation of problem (6) we have that:

• •

$\mathbf{n}$ is the unit outward normal vector to the boundary $\partial\Omega$.

• •

$K>0\ ({\rm m}^{2}\,{\rm s}^{-1})$ is the thermal diffusivity of the fluid, that is, $K=\frac{\alpha}{\rho\,c_{p}}$, where $\alpha\ ({\rm W}\,{\rm m}^{-1}\,{\rm K}^{-1})$ is the thermal conductivity, $\rho\ ({\rm g}\,{\rm m}^{-3})$ is the density, and $c_{p}\ ({\rm W}\,{\rm s}\,{\rm g}^{-1}\,{\rm K}^{-1})$ is the specific heat capacity of water.

• •

$b_{1}^{K}\geq 0\ ({\rm m}\,{\rm s}^{-1})$, for $K\in\{N,S\}$, are the coefficients related to convective heat transfer through the boundaries $\Gamma_{N}$ and $\Gamma_{S}$, obtained from the relation $\rho\,c_{p}\,b_{1}^{K}=h^{K}$, where $h^{K}\geq 0\ ({\rm W}\,{\rm m}^{-2}\,{\rm K}^{-1})$ are the convective heat transfer coefficients on each surface.

• •

$b_{2}^{S}>0\ ({\rm m}\,{\rm s}\,{\rm K}^{-3})$ is the coefficient related to radiative heat transfer through the boundary $\Gamma_{S}$, given by $b_{2}^{S}=\frac{\sigma_{B}\,\varepsilon}{\rho\,c_{p}}$, where $\sigma_{B}\ ({\rm W}\,{\rm m}^{-2}\,{\rm K}^{-4})$ is the Stefan-Boltzmann constant and $\varepsilon$ is the emissivity.

• •

$\theta^{0}\geq 0\ ({\rm K})$ represents the initial temperature.

• •

$\theta_{N},\,\theta_{S}\geq 0\ ({\rm K})$ are the temperatures related to convection heat transfer on the surfaces $\Gamma_{N}$ and $\Gamma_{S}$.

• •

$T_{r}\geq 0\ ({\rm K})$ is the radiation temperature on the surface $\Gamma_{S}$, derived from the expression $\sigma_{B}\,\varepsilon\,T_{r}^{4}=(1-a)R_{sw,net}+R_{lw,down}$, where $a$ is the albedo, $R_{sw,net}\ ({\rm W}\,{\rm m}^{-2})$ denotes the net incident shortwave radiation on the surface $\Gamma_{R}$, and $R_{lw,down}\ ({\rm W}\,{\rm m}^{-2})$ denotes the downwelling longwave radiation.

2.4 The eutrophication model

We consider the following system for modelling the eutrophication processes, based in Michaelis-Menten kinetics (further details can be found, for instance, in fran1 ; DRAGO200117 and the references therein), where we consider the concentrations of five different species: $u^{1}(\mathbf{x},t)\ ({\rm mg}\,{\rm l}^{-1})$ stands for the nutrient (nitrogen in this case), $u^{2}(\mathbf{x},t)\ ({\rm mgC}\,{\rm l}^{-1})$ for the phytoplankton, $u^{3}(\mathbf{x},t)\ ({\rm mgC}\,{\rm l}^{-1})$ for the zooplankton, $u^{4}(\mathbf{x},t)\ ({\rm mgC}\,{\rm l}^{-1})$ for the organic detritus, and $u^{5}(\mathbf{x},t)\ ({\rm mg}\,{\rm l}^{-1})$ for the dissolved oxygen:

[TABLE]

where, for $i=1,\ldots,5$,

[TABLE]

and, for $k=1,\ldots,N_{CT}$, and $i=1,\ldots,5$,

[TABLE]

Finally, the reaction term $\mathbf{A}=(A^{i}):\Omega\times(0,T)\times\mathbb{R}^{6}\rightarrow\mathbb{R}^{5}$ is defined by the following expression:

[TABLE]

where:

• •

$C_{oc}\geq 0\ ({\rm mg}\,{\rm mgC}^{-1})$ is the oxygen-carbon stoichiometric relation,

• •

$C_{nc}\geq 0\ ({\rm mg}\,{\rm mgC}^{-1})$ is the nitrogen-carbon stoichiometric relation,

• •

$C_{fz}\geq 0$ is the zooplankton grazing efficiency factor,

• •

$K_{rd}\geq 0\ ({\rm s}^{-1})$ is the detritus regeneration rate,

• •

$K_{r}\geq 0\ ({\rm s}^{-1})$ is the phytoplankton endogenous respiration rate,

• •

$K_{mf}\geq 0\ ({\rm s}^{-1})$ is the phytoplankton death rate,

• •

$K_{mz}\geq 0\ ({\rm s}^{-1})$ is the zooplankton death rate (including predation),

• •

$K_{z}\geq 0\ ({\rm s}^{-1})$ is the zooplankton predation (grazing),

• •

$K_{F}>0\ ({\rm mgC}\,{\rm l}^{-1})$ is the phytoplankton half-saturation constant,

• •

$K_{N}>0\ ({\rm mg}\,{\rm l}^{-1})$ is the nitrogen half-saturation constant,

• •

$\mu^{i}\geq 0\ ({\rm m}^{2}\,{\rm s}^{-1})$, $i=1,\ldots,5$, are the diffusion coefficients of each species,

• •

$D$ is the thermic regeneration function for the organic detritus, defined as:

[TABLE]

with $\log(\Theta)\ ({\rm K}^{-1})$ the thermic regeneration constant for the reference temperature $\theta^{0}$. In order to simplify the mathematical analysis of the state equations we will consider the following linear approximation:

[TABLE]

if $\Theta>0$, and $D(\theta)=1$ if $\Theta=0$.

• •

$L$ is the luminosity function, given by:

[TABLE]

with $I^{0}\ ({\rm W}\,{\rm m}^{-2})$ the incident light intensity, $I_{s}\ ({\rm W}\,{\rm m}^{-2})$ the light saturation, $\log(C_{t})\ ({\rm K}^{-1})$ the phytoplankton growth thermic constant for the reference temperature $\theta^{0}$, $\varphi_{1}\ ({\rm m}^{-1})$ the light attenuation due to depth, and $\mu\ ({\rm s}^{-1})$ the maximum phytoplankton growth rate. Again, for the sake of simplicity, we will consider the following linear approximation:

[TABLE]

if $C_{t}>0$, and $L(\mathbf{x},t,\theta)=\mu\,\frac{I^{0}(t)}{I_{s}}\,e^{-\varphi_{1}x_{3}}$ if $C_{t}=0$.

2.5 The optimal control problem

Our main objective is to ensure that the concentration of dissolved oxygen in the bottom layer is in an admissible range 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

[TABLE]

where

[TABLE]

is the admissible set, $c>0$ is a constant related to technological limitations of the pumps, $J(\mathbf{g})$ is the cost function:

[TABLE]

and $\lambda^{m},\lambda^{M}>0$ represent, respectively, minimum and maximum permissible concentrations in the control domain $\Omega_{C}$. Finally, $(\mathbf{v},\theta,\mathbf{u})$ are the solutions of the coupled state systems (3), (6) and (10).

3 Mathematical analysis of the state equations

In order to establish the appropriate framework for mathematically analyzing the coupled state systems (3), (6) and (10), we consider, for a Banach space $V_{1}$ and a locally convex space $V_{2}$ such that $V_{1}\subset V_{2}$, and for $1\leq p,q\leq\infty$, the following Sobolev-Bochner space (cf. Chapter 7 of Roubicek1 for further details):

[TABLE]

where $\frac{du}{dt}$ denotes the derivative of $u$ in the sense of distributions. It is well known that, if both $V_{1}$ and $V_{2}$ are Banach spaces, then $W^{1,p,q}(0,T;V_{1},V_{2})$ is also a Banach space endowed with the norm $\|u\|_{W^{1,p,q}(0,T;V_{1},V_{2})}=\|u\|_{L^{p}(0,T;V_{1})}+\left\|\frac{du}{dt}\right\|_{L^{q}(0,T;V_{2})}$.

So, for the modified Navier-Stokes system (3) we consider the following spaces:

[TABLE]

Then, associated to the previous spaces, we define:

[TABLE]

Now, for the water temperature system (6), we consider the following spaces:

[TABLE]

If we define the following norm associated to above space $X_{2}$:

[TABLE]

we have that $X_{2}$ is a reflexive separable Banach space (cf. Lemma 3.1 of Zolesio1 ), and that $\widetilde{X}_{2}\subset L^{2}(\Omega)\subset\widetilde{X}_{2}^{\prime}$ is an evolution triple. So, we consider:

[TABLE]

Finally, for the eutrophication system (10), we define:

[TABLE]

and we consider the following spaces associated to them:

[TABLE]

From this section we will assume the following hypotheses for coefficients and data in the analytical study of the problem:

• •

$g^{k}\in H^{1}(0,T)$, with $g^{k}(0)=0$, $\forall k=1,\ldots,N_{CT}$,

• •

$\mathbf{v}^{0}\in\left[H_{\sigma}^{2}(\Omega)\right]^{3}=\{\mathbf{v}\in[H^{2}(\Omega)]^{3}\,:\,\nabla\cdot\mathbf{v}=0,\,\mathbf{v}_{|_{\partial\Omega}}=\boldsymbol{0}\}\subset\widetilde{\mathbf{X}}_{1}$,

• •

$\theta_{S}\in L^{2}(0,T;L^{2}(\Gamma_{S}))$,

• •

$\theta_{N}\in L^{2}(0,T;L^{2}(\Gamma_{N}))$,

• •

$T_{r}\in L^{5}(0,T;L^{5}(\Gamma_{S}))$,

• •

$\theta^{0}\in X_{2}$,

• •

$I_{0}\in L^{\infty}(0,T)$,

• •

$\mathbf{u}^{0}\in\mathbf{X}_{3}$.

Under these hypotheses we will state now two lemmas (whose demonstrations can be found in fran7 and fran8 , respectively), which will allow us to reformulate the state systems (3), (6) and (10) as homogeneous Dirichlet problems.

Lemma 1

There exists a linear continuous extension:

[TABLE]

such that ${\boldsymbol{\zeta}_{\mathbf{g}}}_{|_{\partial\Omega}}=\boldsymbol{\phi}_{\mathbf{g}}$, where $\boldsymbol{\phi}_{\mathbf{g}}$ is defined by (2), and $H_{\sigma}^{2}(\Omega)=\{\mathbf{u}\in[H^{2}(\Omega)]^{3}\,:\,\nabla\cdot\mathbf{u}=0\}$. ∎

Remark 1

It is worthwhile emphasizing here that, thanks to the construction done in the proof of Lemma 1, we have

[TABLE]

and, consequently, this term will not appear in the corresponding variational formulation. ∎

Lemma 2

We have that the following operator is compact

[TABLE]

where:

[TABLE]

with $\gamma^{k}_{\mathbf{h}}(s)\in L^{2}(-T,T)$, for $k=1,\ldots,N_{CT}$, defined by:

[TABLE]

and $\beta_{0}:u\in H^{3/2}(\partial\Omega)\rightarrow\beta_{0}(u)\in H^{2}(\Omega)$ the right inverse of the classical trace operator $\gamma_{0}$, i.e., such that $(\gamma_{0}\circ\beta_{0})(u)=u$ (cf. Theorem 8.3. of magenes1 ).

We also have the existence of a constant $C$, that depends continuously on the space-time computational domain and the initial temperature $\theta^{0}$, such that:

[TABLE]

∎

Now, we will establish the following notations, in order to stablish the homogeneous Dirichlet systems. Given elements $(\mathbf{z},\xi,\mathbf{w})\in\widetilde{\mathbf{W}}_{1}\times\widetilde{W}_{2}\times\widetilde{\mathbf{W}}_{3}$, we define $(\mathbf{v},\theta,\mathbf{u})\in\mathbf{W}_{1}\times W_{2}\times\mathbf{W}_{3}$ in the the following way:

• •

$\mathbf{v}=\mathbf{z}+\boldsymbol{\zeta}_{\mathbf{g}}\in\mathbf{W}_{1}$, with $\boldsymbol{\zeta}_{\mathbf{g}}\in W^{1,2,2}(0,T;[H_{\sigma}^{2}(\Omega)]^{3},[H_{\sigma}^{2}(\Omega)]^{3})$ the extension of control $\mathbf{g}$ given by Lemma 1.

• •

$\theta=\xi+\zeta_{\mathbf{h}_{\theta}}\in W_{2}$, with $\displaystyle\zeta_{\mathbf{h_{\theta}}}\in W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$ the extension of $\mathbf{h}_{\theta}$ obtained from Lemma 2, where:

[TABLE]

• •

$u^{i}=w^{i}+\zeta_{\mathbf{h}^{i}_{\mathbf{u}}}$, with $\zeta_{\mathbf{h}^{i}_{\mathbf{u}}}\in W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$ the extension of $\mathbf{h}^{i}_{\mathbf{u}}$ obtained from Lemma 2 with obvious modifications, where:

[TABLE]

As it is immediate, $\mathbf{w}\in\mathbf{W}_{3}$.

Thus, using above notations, we can reformulate the state systems (3), (6) and (10) in the following way:

[TABLE]

[TABLE]

[TABLE]

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 (3), (6) and (10) in terms of the modified state systems (34), (35) and (36). It should be also noted that, in the case of systems (6) and (10), the coupling terms in the Dirichlet boundary conditions are now transferred to the partial differential equations in systems (35) and (36).

Definition 1 (The concept of solution for the state systems)

An element $(\mathbf{v},\theta,\mathbf{u})\in\mathbf{W}_{1}\times W_{2}\times\mathbf{W}_{3}$ is a solution for the state systems (3), (6) and (10), if there exists an element $(\mathbf{z},\xi,\mathbf{w})\in\widetilde{\mathbf{W}}_{1}\times\widetilde{W}_{2}\times\widetilde{\mathbf{W}}_{3}$ such that:

• •

$\mathbf{v}=\mathbf{z}+\boldsymbol{\zeta}_{\mathbf{g}}$, with $\boldsymbol{\zeta}_{\mathbf{g}}\in W^{1,2,2}(0,T;[H_{\sigma}^{2}(\Omega)]^{3},[H_{\sigma}^{2}(\Omega)]^{3})$ as given by Lemma 1, $\mathbf{z}(0)=\mathbf{v}^{0}$, a.e. $\mathbf{x}\in\Omega$, and $\mathbf{z}\in\widetilde{\mathbf{W}}_{1}$ the solution of the following variational formulation:

[TABLE]

where:

[TABLE]

• •

$\theta=\xi+\zeta_{\mathbf{h}_{\theta}}$, with $\displaystyle\zeta_{\mathbf{h_{\theta}}}\in W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$ obtained from Lemma 2 for $\mathbf{h}_{\theta}\in[L^{2}(0,T)]^{N_{CT}}$ defined by (32), $\xi(0)=\theta^{0}-\zeta_{\mathbf{h}_{\theta}}(0)$, a.e. $\mathbf{x}\in\Omega$, and $\xi\in\widetilde{W}_{2}$ the solution of the following variational formulation:

[TABLE]

where:

[TABLE]

• •

$u^{i}=w^{i}+\zeta_{\mathbf{h}^{i}_{\mathbf{u}}}$, with $\zeta_{\mathbf{h}^{i}_{\mathbf{u}}}\in W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$ obtained from Lemma 2 for $\mathbf{h}^{i}_{\mathbf{u}}\in[L^{2}(0,T)]^{N_{CT}}$ defined by (33), $\mathbf{w}(0)=\mathbf{u}_{0}-\boldsymbol{\zeta}_{\mathbf{h}_{\mathbf{u}}}(0)$, a.e. $\mathbf{x}\in\Omega$, and $\mathbf{w}\in\widetilde{\mathbf{W}}_{3}$ the solution of the following variational formulation:

[TABLE]

where $\Lambda_{\mu}={\rm diag}(\mu^{1},\ldots,\mu^{5})\in M_{5\times 5}(\mathbb{R})$ is a diagonal matrix with diffusion coefficients, and:

[TABLE]

∎

Remark 2

We have the following dependence scheme between the elements of state system:

[TABLE]

Therefore, we can separate the mathematical analysis of systems (3)-(6) from system (10). The coupled system (3)-(6) has been fully analyzed by the authors in fran7 and fran8 . Thus, following the results there, we can assure that, for each control $\mathbf{g}\in[H^{1}(0,T)]^{N_{CT}}$, there exists a solution $(\mathbf{v},\theta)\in\widetilde{\mathbf{W}}_{1}\times\widetilde{W}_{2}$ of the thermo-hydrodynamic system (3)-(6). We must remark here that, due to the complexity of this nonlinear system, we cannot obtain a uniqueness result for the thermo-hydrodynamic solution $(\mathbf{v},\theta)$ 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 in analyzing the solution $\mathbf{u}$ of the eutrophication system (10) or, equivalently, in studying the solution $\mathbf{w}$ of the modified system (36). ∎

Thus, in order to analyze the existence of a solution $\mathbf{u}$ by means of a fixed point technique, we consider the operator:

[TABLE]

where $\mathbf{u}^{*}=({u^{1}}^{*},\ldots,{u^{5}}^{*})$, with ${u^{i}}^{*}\in L^{2}(0,T;L^{2}(\Omega))$, for $i=1,\ldots,5$, $\mathbf{h}^{*}_{\mathbf{u}}=({\mathbf{h}^{1}_{\mathbf{u}}}^{*},\ldots,{\mathbf{h}^{5}_{\mathbf{u}}}^{*})$, with ${\mathbf{h}^{i}_{\mathbf{u}}}^{*}\in[L^{2}(0,T)]^{N_{CT}}$, for $i=1,\ldots,5$, $\mathbf{u}=(u^{1},\ldots,u^{5})\in\mathbf{W}_{3}$, $\mathbf{h}_{\mathbf{u}}=(\mathbf{h}^{1}_{\mathbf{u}},\ldots,\mathbf{h}^{5}_{\mathbf{u}})\in[L^{2}(0,T)]^{5\times N_{CT}}$, such that:

• •

$\zeta_{{\mathbf{h}^{i*}_{\mathbf{u}}}}\in W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$, for $i=1,\ldots,5$, is defined by Lemma 2.

• •

$\mathbf{u}\in{\mathbf{W}}_{3}$ is the solution, in the sense of Definition 1 with the obvious modifications, of the following decoupled problem with resolution order $i=3\rightarrow 2\rightarrow 4\rightarrow 1\rightarrow 5$:

[TABLE]

where the Caratheodory function $\widehat{\mathbf{A}}=(\widehat{A}^{i}):\Omega\times(0,T)\times\mathbb{R}^{6}\times\mathbb{R}^{6}\rightarrow\mathbb{R}^{5}$ is defined by:

[TABLE]

• •

$\mathbf{h}^{i}_{\mathbf{u}}\in[L^{2}(0,T)]^{N_{CT}}$, for $i=1,\ldots,5$, is such that:

[TABLE]

Remark 3

All the five equations of the decoupled system (44) can be expressed in the way of the following generic equation:

[TABLE]

where $\Gamma_{1}=\Gamma_{S}\cup\Gamma_{C}\cup\Gamma_{N}$, $\Gamma_{2}=\Gamma_{T}$, $w^{0}\in H^{1}(\Omega)$, $k_{1}\in L^{4}(0,T;L^{3}(\Omega))$, $k_{2}\in L^{2}(0,T;L^{3/2}(\Omega))$ and $k_{3}\in L^{2}(0,T;L^{2}(\Gamma_{1}))$.

In one hand, $k_{3}\in L^{2}(0,T;L^{2}(\Gamma_{1}))$ thanks to Lemma 2 and Theorem 3.1 of fursikov1 for the trace operator $\gamma_{1}=\partial^{1}/\partial\mathbf{n}^{1}$, as well $\|\gamma_{1}(\zeta_{{\mathbf{h}^{i*}_{\mathbf{u}}}})\|_{L^{2}(0,T;L^{2}(\Gamma))}\leq C_{3}\|{\mathbf{h}^{i*}_{\mathbf{u}}}\|_{[L^{2}(0,T)]^{N_{CT}}}$, $\forall i=1,\ldots,5$.

In the other hand, for the coefficients $k_{1}$ and $k_{2}$, we need to study the particular case for each one of the five species. We have the following lemma.

Lemma 3

We have the following estimates for the coefficients $k_{1}$ and $k_{2}$ associated to each species:

• •

Species $u^{3}$:

[TABLE]

• •

Species $u^{2}$:

[TABLE]

• •

Species $u^{4}$:

[TABLE]

• •

Species $u^{1}$, $k_{1}=0$ and:

[TABLE]

• •

Species $u^{5}$, $k_{1}=0$ and:

[TABLE]

Where $C_{1}$ and $C_{2}$ are positive constants that depend on coefficients and data associated to problem (10).

Proof

We will follow the same order of resolution of the decoupled problem:

• •

Equation for $u^{3}$:

[TABLE]

So, we have that $k_{1}\in L^{\infty}(0,T;L^{\infty}(\Omega))$ and $k_{2}\in L^{2}(0,T;L^{2}(\Omega))$. We also have the following estimates:

[TABLE]

where we have used $\|\zeta_{{\mathbf{h}_{\mathbf{u}}^{3*}}}\|_{W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))}\leq C(u^{0,3})\|{\mathbf{h}_{\mathbf{u}}^{3*}}\|_{[L^{2}(0,T)]^{N_{CT}}}$, in particular, $\|\nabla\zeta_{{\mathbf{h}_{\mathbf{u}}^{3*}}}\|_{L^{\infty}(0,T;L^{6}(\Omega))}\leq C(u^{0,3})\|{\mathbf{h}_{\mathbf{u}}^{3*}}\|_{[L^{2}(0,T)]^{N_{CT}}}$

• •

Equation for $u^{2}$:

[TABLE]

In this case, the regularity of the term $k_{1}$ is imposed by the regularity of the term $L(\mathbf{x},t,\theta)$. In particular, we have $W_{2}\subset L^{2}(0,T;L^{6}(\Omega))\cap L^{\infty}(0,T;L^{2}(\Omega))$ $\subset L^{4}(0,T;L^{3}(\Omega))$ and then, $k_{1}\in L^{4}(0,T;L^{3}(\Omega))$. In the other hand, it is clear that $k_{2}\in L^{2}(0,T;L^{2}(\Omega))$. Finally, we have the following estimates:

[TABLE]

where we have used $\|\zeta_{{\mathbf{h}_{\mathbf{u}}^{2*}}}\|_{W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))}\leq C(u^{0,2})\|{\mathbf{h}_{\mathbf{u}}^{2*}}\|_{[L^{2}(0,T)]^{N_{CT}}}$, in particular, $\|\zeta_{{\mathbf{h}_{\mathbf{u}}^{2*}}}\|_{L^{\infty}(0,T;L^{\infty}(\Omega))},\,\|\nabla\zeta_{{\mathbf{h}_{\mathbf{u}}^{2*}}}\|_{L^{\infty}(0,T;L^{6}(\Omega))}\leq C(u^{0,2})\|{\mathbf{h}_{\mathbf{u}}^{2*}}\|_{[L^{2}(0,T)]^{N_{CT}}}$.

• •

Equation for $u^{4}$:

[TABLE]

The regularity of the term $k_{1}$ is given by the regularity of the water temperature. So, $k_{1}\in L^{4}(0,T;L^{3}(\Omega))$. In the case of $k_{2}$ we are in the same situation as above and then $k_{2}\in L^{2}(0,T;L^{2}(\Omega))$. We also have the following estimates:

[TABLE]

• •

Equation for $u_{1}$:

[TABLE]

In the term $k_{2}$, the most restrictive regularity is determined by the product of two functions ($u^{4}$ and $\theta$) lying in the space $L^{4}(0,T;L^{3}(\Omega))$. So, $k_{2}\in L^{2}(0,T;L^{3/2}(\Omega))$. Besides, we have the following estimate:

[TABLE]

• •

Equation for $u_{5}$:

[TABLE]

We are in the same situation as above, in particular, $k_{2}\in L^{2}(0,T;L^{3/2}(\Omega))$, and we also have the following estimate:

[TABLE]

∎

Now we can state the following existence result for the generic equation (47). The proof of this result can be done using techniques analogous to those presented in fran1 :

Theorem 3.1

For any given elements $\mathbf{v}\in\mathbf{W}_{3}$, $w_{0}\in H^{1}(\Omega)$, $k_{1}\in L^{4}(0,T;L^{3}(\Omega))$, $k_{2}\in L^{2}(0,T;L^{3/2}(\Omega))$ and $k_{3}\in L^{2}(0,T;L^{2}(\Gamma_{1}))$, there exists a unique element $w\in W^{1,2,2}(0,T;H^{1}_{0,\Gamma_{2}}(\Omega),{H^{1}_{0,\Gamma_{2}}(\Omega)}^{\prime})\cap L^{\infty}(0,T;L^{2}(\Omega))$, with $w(0)=w^{0}$ a.e. $x\in\Omega$, that satisfies the following variational formulation:

[TABLE]

where $H^{1}_{0,\Gamma_{2}}(\Omega)=\{\eta\in H^{1}(\Omega):\;\eta_{|_{\Gamma_{2}}}=0\}$. The solution also verifies the following estimate:

[TABLE]

where $C$ is a positive constant depending on $k_{1}$, $w^{0}$ and $\mathbf{v}$. ∎

Thus, thanks to Remark 3 and Theorem 3.1, we have the following estimates:

Lemma 4

A solution $\mathbf{u}$ of the uncoupled system (44) verifies the following:

• •

Estimates for $u^{3}$ and $\mathbf{h}^{3}_{\mathbf{u}}$:

[TABLE]

where the constants $C_{3}$ and $C_{8}$ also depend on the initial condition $\|u^{0,3}\|_{H^{1}(\Omega)}$.

• •

Estimates for $u^{2}$ and $\mathbf{h}^{2}_{\mathbf{u}}$:

[TABLE]

where constants $C_{2}$ and $C_{7}$ also depend on the initial conditions $\|u^{0,2}\|_{H^{1}(\Omega)}$ and $\|u^{0,3}\|_{H^{1}(\Omega)}$.

• •

Estimates for $u^{4}$ and $\mathbf{h}^{4}_{\mathbf{u}}$:

[TABLE]

where the constants $C_{4}$ and $C_{9}$ also depend on the initial conditions $\|u^{0,2}\|_{H^{1}(\Omega)}$, $\|u^{0,3}\|_{H^{1}(\Omega)}$ and $\|u^{0,4}\|_{H^{1}(\Omega)}$.

• •

Estimates for $u^{1}$ and $\mathbf{h}^{1}_{\mathbf{u}}$:

[TABLE]

where constants $C_{1}$ and $C_{6}$ also depend on the initial conditions $\|u^{0,1}\|_{H^{1}(\Omega)}$, $\|u^{0,2}\|_{H^{1}(\Omega)}$, $\|u^{0,3}\|_{H^{1}(\Omega)}$ and $\|u^{0,4}\|_{H^{1}(\Omega)}$.

• •

Estimates for $u^{5}$ and $\mathbf{h}^{5}_{\mathbf{u}}$:

[TABLE]

where the constants $C_{5}$ and $C_{10}$ also depend on the initial conditions $\|u^{0,2}\|_{H^{1}(\Omega)}$, $\|u^{0,3}\|_{H^{1}(\Omega)}$, $\|u^{0,4}\|_{H^{1}(\Omega)}$ and $\|u^{0,5}\|_{H^{1}(\Omega)}$. ∎

Remark 4

We must to note here that all above estimates do not depend on the variable $\mathbf{u}^{*}$, since the dependence on $\mathbf{u}^{*}$ appears within terms of the form:

[TABLE]

with $K>0$, and those terms are bounded a.e. $(\mathbf{x},t)\in\Omega\times(0,T)$ by a constant independent on $\mathbf{u}^{*}$. ∎

Now, we will prove the main result of this Section:

Theorem 3.2 (Existence of solution for the eutrophication system)

If there exist coefficients and data such that:

[TABLE]

for all $\mathbf{g}\in\mathcal{U}_{ad}$, then there will exist positive constants $\widetilde{C}_{i}$, $i=1,\ldots,10$, such that the operator $M_{\mathbf{u}}:\mathbf{B}_{\mathbf{u}}\rightarrow\mathbf{B}_{\mathbf{u}}$ defined by (43) has a fixed point, which is solution of the state system (3), (6) and (10) in the sense of Definition 1, where:

[TABLE]

Proof

In order to apply the Schauder fixed point Theorem (see, for instance, Theorem 9.5 of conway1 ), we will prove that the operator $M_{\mathbf{u}}$ is compact and that there exist positive constants $\widetilde{C}_{i}$, $i=1,\ldots,10$, such that the operator $M_{\mathbf{u}}$ maps elements from the set $\mathbf{B}_{\mathbf{u}}$ (which is closed and convex) into itself.

• •

The operator $M_{\mathbf{u}}$ is compact in the sense that it is continuous and $\overline{M_{\mathbf{u}}(A)}$ is compact whenever $A$ is a bounded subset of $[L^{2}(0,T;L^{2}(\Omega)]^{5}\times[L^{2}(0,T)]^{5\times N_{CT}}$:

In fact, given a convergent sequence $\{(\mathbf{u}^{*}_{n},\mathbf{h}_{n}^{*})\}_{n\in\mathbb{N}}\subset[L^{2}(0,T;L^{2}(\Omega)]^{5}\times[L^{2}(0,T)]^{5\times N_{CT}}$ such that $\mathbf{u}^{*}_{n}\rightarrow\mathbf{u}^{*}$ in $[L^{2}(0,T;L^{2}(\Omega))]^{5}$ and $\mathbf{h}_{n}^{*}\rightarrow\mathbf{h}_{\mathbf{u}}^{*}$ in $[L^{2}(0,T)]^{5\times N_{CT}}$, we have that $(\mathbf{u}_{n},\mathbf{h}_{n})=M_{\mathbf{u}}(\mathbf{u}^{*}_{n},\mathbf{h}_{n}^{*})\in\mathbf{W}_{3}\times[H^{1}(0,T)]^{5\times N_{CT}}$ is such that $\mathbf{u}_{n}=\mathbf{w}_{n}+\boldsymbol{\zeta}_{\mathbf{h}_{n}^{*}}$, with $\mathbf{w}_{n}\in\widetilde{\mathbf{W}}_{3}$ the solution of the following variational formulation:

[TABLE]

where:

[TABLE]

Using Lemma 2, we know that $\boldsymbol{\zeta}_{\mathbf{h}_{n}^{*}}\rightarrow\boldsymbol{\zeta}_{\mathbf{h}^{*}}$ in $W^{1,2,2}(0,T;[H^{2}(\Omega)]^{5},[H^{2}(\Omega)]^{5})$, $\mathbf{H}_{n}\rightarrow\mathbf{H}_{\mathbf{u}}$ in $[L^{2}(0,T;[L^{2}(\Omega))]^{5}$, and $\mathbf{g}_{n}\rightarrow\mathbf{g}_{\mathbf{u}}$ in $[L^{2}(0,T)]^{5\times N_{CT}}$. Thus, taking subsequences if necessary, we have that $\mathbf{u}_{n}^{*}\rightarrow\mathbf{u}^{*}$ a.e. $(x,t)\in\Omega\times(0,T)$, $\frac{\partial\mathbf{w}_{n}}{\partial t}\rightharpoonup\frac{\partial\mathbf{w}}{\partial t}$ weakly in $L^{2}(0,T;\widetilde{\mathbf{X}}_{3}^{\prime})$, $\mathbf{w}_{n}\rightharpoonup\mathbf{w}$ weakly in $L^{2}(0,T;\widetilde{\mathbf{X}}_{3})$, and $\mathbf{w}_{n}\rightarrow\mathbf{w}$ strongly in $[L^{10/3-\epsilon}(0,T,L^{10/3-\epsilon}(\Omega))]^{5}$, for all $\epsilon>0$. It is straightforward to prove, using previous convergences, that we can pass to the limit in variational formulation (72) obtaining that $M_{\mathbf{u}}(\mathbf{u}_{n}^{*},\mathbf{h}^{*}_{n})\rightarrow M_{\mathbf{u}}(\mathbf{u}^{*},\mathbf{h}^{*})$ in $[L^{2}(0,T;L^{2}(\Omega))]^{5}\times[L^{2}(0,T)]^{5\times N_{CT}}$, with $(\mathbf{u},\mathbf{h}_{\mathbf{u}})=M_{\mathbf{u}}(\mathbf{u}^{*},\mathbf{h}^{*})$. We also derive that $\mathbf{u}=\mathbf{w}+\boldsymbol{\zeta}_{\mathbf{h}_{\mathbf{u}}^{*}}\in\mathbf{W}_{3}$, with $\mathbf{w}\in\widetilde{\mathbf{W}}_{3}$ the solution of the following variational formulation:

[TABLE]

where:

[TABLE]

Finally, the compactness of operator $M_{\mathbf{u}}$ is a direct consequence of the compact embedding of space $\mathbf{W}_{3}\times[H^{1}(0,T)]^{5\times N_{CT}}$ into space $[L^{2}(0,T;L^{2}(\Omega))]^{5}\times[L^{2}(0,T)]^{5\times N_{CT}}$.

• •

There exist positive constants $\widetilde{C}_{i}$, $i=1,\ldots,10$, such that the operator $M_{\mathbf{u}}$ applies elements from the set $\mathbf{B}_{\mathbf{u}}$ into itself:

We only need to prove that, if we define the following constants:

[TABLE]

then $(\mathbf{u},\mathbf{h}_{\mathbf{u}})=M_{\mathbf{u}}(\mathbf{u}^{*},\mathbf{h}_{\mathbf{u}}^{*})\in\mathbf{B}_{\mathbf{u}}$, for all $(\mathbf{u}^{*},\mathbf{h}_{\mathbf{u}}^{*})\in\mathbf{B}_{\mathbf{u}}$. Indeed, given an element $(\mathbf{u}^{*},\mathbf{h}_{\mathbf{u}}^{*})\in\mathbf{B}_{\mathbf{u}}$, we have that, thanks to the estimates (60)-(69), we can easily obtain by simple (but tedious) algebraic computations that $(\mathbf{u},\mathbf{h}_{\mathbf{u}})=M_{\mathbf{u}}(\mathbf{u}^{*},\mathbf{h}_{\mathbf{u}}^{*})$ satisfies the following estimates:

[TABLE]

[TABLE]

Thus, $M_{\mathbf{u}}(\mathbf{u}^{*},\mathbf{h}_{\mathbf{u}}^{*})\in\mathbf{B}_{\mathbf{u}}$

Then, as a direct consequence of the Schauder Theorem, we obtain the existence of a fixed point $(\mathbf{u},\mathbf{h}_{\mathbf{u}})\in\mathbf{B}_{\mathbf{u}}$, which is, from the construction of operator $M_{\mathbf{u}}$, a solution of problem (10). ∎

4 Mathematical analysis of the optimal control problem

In this section we will prove the existence of solution for the optimal control problem $(\mathcal{P})$. It is important to remark here that, since we have not demonstrated the uniqueness of solution for the state systems (3), (6) and (10), we will treat the problem as a multistate control problem (cf. fran9 ). Thus, we define the set:

[TABLE]

where the set of admissible controls $\mathcal{U}_{ad}$ is bounded, convex and closed (in particular, $\mathcal{U}_{ad}$ is weakly closed). We observe that the constraints in the sets $\mathcal{U}$ and $\mathcal{U}_{ad}$ are well defined since $g^{k}\in H^{1}(0,T)\subset\subset\mathcal{C}([0,T])$, $k=1,\ldots,N_{CT}$, and ${u}^{i}\in\mathcal{C}([0,T];L^{2}(\Omega))$, $i=1,\ldots,5$. Then, we prove the following property for the set $\mathcal{U}$.

Lemma 5

The set $\mathcal{U}$ is weakly closed.

Proof

Let us consider a sequence of elements $\{(\mathbf{v}_{n},\theta_{n},\mathbf{u}_{n},\mathbf{g}_{n})\}_{n\in\mathbb{N}}\subset\mathcal{U}$ such that $(\mathbf{v}_{n},\theta_{n},\mathbf{u}_{n},\mathbf{g}_{n})\rightharpoonup(\mathbf{v},\theta,\mathbf{u},\mathbf{g})$ in $L^{3}(0,T;\mathbf{X}_{1})\times L^{2}(0,T;X_{2})\times L^{2}(0,T;\mathbf{X}_{3})\times\mathcal{U}_{ad}.$ In particular, the sequence $\{\mathbf{g}_{n}\}_{n\in\mathbb{N}}$ is bounded in $[H^{1}(0,T)]^{N_{CT}}$ and then, thanks to the estimates obtained in Lemma 7 of fran7 , in Theorem 9 of fran8 , and in above Lemmas 1, 2 and 4, we have that the sequence $\{(\mathbf{z}_{n},\boldsymbol{\zeta}_{\mathbf{g}_{n}},\xi_{n},\zeta_{\mathbf{h}_{\theta_{n}}}\mathbf{w}_{n},\boldsymbol{\zeta}_{\mathbf{h}_{\mathbf{u}_{n}}})\}_{n\in\mathbb{N}}\subset\widetilde{\mathbf{W}}_{1}\times W^{1,2,2}(0,T;[H_{\sigma}^{2}(\Omega)]^{3},[H_{\sigma}^{2}(\Omega)]^{3})\times\widetilde{W}_{2}\times W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))\times\widetilde{\mathbf{W}}_{3}\times W^{1,2,2}(0,T;[H^{2}(\Omega)]^{5},[H^{2}(\Omega)]^{5})$ induced by Definition 1 is bounded, where, for all $n\in\mathbb{N}$, $\mathbf{v}_{n}=\mathbf{w}_{n}+\boldsymbol{\zeta}_{\mathbf{g}_{n}}$, $\theta_{n}=\xi_{n}+\zeta_{\mathbf{h}_{\theta_{n}}}$ and $\mathbf{u}_{n}=\mathbf{w}_{n}+\boldsymbol{\zeta}_{\mathbf{h}_{\mathbf{u}_{n}}}.$

Now, if we denote by $\mathbf{w}=\mathbf{v}-\boldsymbol{\zeta}_{\mathbf{g}}$, $\xi=\theta-\zeta_{\mathbf{h}_{\theta}}$ and $\mathbf{w}=\mathbf{u}-\boldsymbol{\zeta}_{\mathbf{h}_{\mathbf{u}}}$, we have (taking subsequences if necessary) the following convergences for the elements associated to the sequence of controls:

• •

$g_{n}(t)\rightarrow g(t)$ strongly for all $t\in[0,T]$ (so, in particular, $g_{n}(0)\rightarrow g(0)$ and, consequently, $g(0)=0$),

• •

$\displaystyle\boldsymbol{\zeta}_{\mathbf{g}_{n}}\rightharpoonup\boldsymbol{\zeta}_{\mathbf{g}}$ weakly in $W^{1,2,2}(0,T;[H_{\sigma}^{2}(\Omega)]^{3},[H_{\sigma}^{2}(\Omega)]^{3})$,

• •

$\displaystyle\zeta_{\mathbf{h}_{\theta_{n}}}\rightarrow\zeta_{\mathbf{h}_{\theta}}$ strongly in $W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$,

• •

$\displaystyle\boldsymbol{\zeta}_{\mathbf{h}_{\mathbf{u}_{n}}}\rightarrow\boldsymbol{\zeta}_{\mathbf{h}_{\mathbf{u}}}$ strongly in $W^{1,2,2}(0,T;[H^{2}(\Omega)]^{5},[H^{2}(\Omega)]^{5})$,

where the first convergence is a direct consequence of compactness of $H^{1}(0,T)$ in $\mathcal{C}([0,T])$ and the two last convergences are consequence of Lemma 2. In a similar way we also have, for the sequence $\{\mathbf{z}_{n}\}_{n\in\mathbb{N}}$, the following convergences:

• •

$\mathbf{z}_{n}\rightarrow\mathbf{z}$ strongly in $L^{p}(0,T;[L^{q}(\Omega)]^{3})$ for all $1<p<\infty$ and $2\leq q<\infty$,

• •

$\displaystyle\mathbf{z}_{n}\rightharpoonup\mathbf{z}$ weakly in $L^{3}(0,T;\widetilde{\mathbf{X}}_{1})$,

• •

$\displaystyle\frac{d\mathbf{z}_{n}}{dt}\rightharpoonup\frac{d\mathbf{z}}{dt}$ weakly in $L^{2}(0,T;[L^{2}(\Omega)]^{3})$,

• •

$\displaystyle\nabla\mathbf{z}_{n}\rightharpoonup^{*}\nabla\mathbf{z}$ weakly-$*$ in $L^{\infty}(0,T;[L^{3}(\Omega)]^{3})$,

• •

$\displaystyle\beta(e(\boldsymbol{\zeta}_{\mathbf{g}_{n}}+\mathbf{z}_{n}))\,e(\boldsymbol{\zeta}_{\mathbf{g}_{n}}+\mathbf{z}_{n})\rightharpoonup\widehat{\beta}$, weakly in $L^{3/2}(0,T;\widetilde{\mathbf{X}}_{1}^{\prime})$.

Moreover, for the sequence $\{\xi_{n}\}_{n\in\mathbb{N}}$ we have:

• •

$\displaystyle\xi_{n}\rightharpoonup\xi$ in $L^{2}(0,T;\widetilde{X}_{2})$,

• •

$\displaystyle\xi_{n}\rightharpoonup^{*}\xi$ in $L^{\infty}(0,T;L^{2}(\Omega))$,

• •

$\displaystyle\xi_{n}\rightarrow\xi$ in $L^{10/3-\epsilon}(0,T;L^{10/3-\epsilon}(\Omega))$, for all $\epsilon>0$ small enough,

• •

$\displaystyle\xi_{n}\rightarrow\xi$ in $L^{2}(0,T;L^{2}(\Gamma_{C}))$,

• •

$\displaystyle\xi_{n}\rightarrow\xi$ in $L^{4}(0,T;L^{4}(\Gamma_{S}))$.

Finally, for the sequence $\{\mathbf{w}_{n}\}_{n\in\mathbb{N}}$, we have the following convergences:

• •

$\displaystyle\mathbf{w}_{n}\rightharpoonup\mathbf{w}$ weakly in $L^{2}(0,T;\widetilde{\mathbf{X}}_{3})$,

• •

$\displaystyle\frac{d\mathbf{w}_{n}}{dt}\rightharpoonup\frac{d\mathbf{w}}{dt}$ weakly in $L^{2}(0,T;\widetilde{\mathbf{X}}_{3}^{\prime})$,

• •

$\displaystyle\mathbf{w}_{n}\rightarrow\mathbf{w}$ strongly in $[L^{10/3-\epsilon}(0,T;L^{10/3-\epsilon}(\Omega))]^{5}$, for all $\epsilon>0$ small enough,

• •

$\displaystyle\mathbf{w}_{n}\rightarrow\mathbf{w}$ strongly in $[L^{2}(0,T;L^{2}(\Gamma_{C}))]^{5}$.

So, we are able to pass to the limit in the corresponding variational formulations using the same arguments that we have employed for proving the compactness of operator $M_{\mathbf{u}}$, and in the Galerkin approximations for systems (3) and (6) (cf. fran7 and fran8 ). The only difficulty here is to prove that $\widehat{\beta}=\beta(e(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z}))\,e(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z})$. However, by Lemma 4.2 of fran7 we have that

[TABLE]

for all $\boldsymbol{\eta}\in L^{3}(0,T;\widetilde{\mathbf{X}}_{1})$, and then, using similar techniques that we can find in the proof of Theorem 4.3 of fran7 , we can prove that

[TABLE]

for all $\boldsymbol{\eta}\in L^{3}(0,T;\widetilde{\mathbf{X}}_{1})$. Finally, choosing $\boldsymbol{\eta}=\mathbf{z}\pm\lambda\boldsymbol{\zeta}$, with $\boldsymbol{\zeta}\in L^{3}(0,T;\widetilde{\mathbf{X}}_{1})$ and $\lambda$ an arbitrary positive number, and multiplying both sides of the inequality by $\lambda^{-1}$, we obtain

[TABLE]

Now, letting $\lambda$ tend to zero, we deduce that, for all $\boldsymbol{\zeta}\in L^{3}(0,T;\widetilde{\mathbf{X}}_{1})$:

[TABLE]

Thus, $\widehat{\beta}=\beta(e(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z}))e(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z})$ a.e. $(\mathbf{x},t)\in\Omega\times]0,T[$, and then, $(\mathbf{v},\theta,\mathbf{u})$ is a solution of the systems (3), (6) and (10) associated to the control $\mathbf{g}$.

Finally, by the strong convergence of $\{\mathbf{u}_{n}\}_{n\in\mathbb{N}}$ in $[L^{2}(0,T;L^{2}(\Gamma_{C}))]^{5}$, we have

[TABLE]

and, consequently, the element $(\mathbf{v},\theta,\mathbf{u},\mathbf{g})\in\mathcal{U}$. ∎

Theorem 4.1 (Existence of optimal solution)

The optimal control problem $(\mathcal{P})$ has, at least, a solution.

Proof

Let us consider a minimizing sequence $\{(\mathbf{v}_{n},\theta_{n},\mathbf{u}_{n},\mathbf{g}_{n})\}_{n\in\mathbb{N}}\subset\mathcal{U}$. Then, $\{\mathbf{g}_{n}\}_{n\in\mathbb{N}}$ is bounded in $[H^{1}(0,T)]^{N_{CT}}$, which implies, thanks to the estimates (66), (62), (60)and (64), and to the Hypotheses of Theorem 3.2, that the sequence $\{\mathbf{u}_{n}\}_{n\in\mathbb{N}}$ is bounded in $\mathbf{W}_{3}$. We also have, thanks to estimates obtained in fran7 and fran8 that the sequence $\{(\mathbf{v}_{n},\theta_{n})\}_{n\in\mathbb{N}}$ is also bounded in $\mathbf{W}_{1}\times W_{2}$. Thus, we can take a subsequence of $\{(\mathbf{v}_{n},\theta_{n},\mathbf{u}_{n},\mathbf{g}_{n})\}_{n\in\mathbb{N}}\subset\mathcal{U}$, still denoted in the same way, such that $(\mathbf{v}_{n},\theta_{n},\mathbf{u}_{n},\mathbf{g}_{n})\rightharpoonup(\widetilde{\mathbf{v}},\widetilde{\theta},\widetilde{\mathbf{u}},\widetilde{\mathbf{g}})$ in $L^{3}(0,T;\mathbf{X}_{1})\times L^{2}(0,T;X_{2})\times L^{2}(0,T;\mathbf{X}_{3})\times\mathcal{U}_{ad}.$ Moreover, from previous Lemma, we have that $(\widetilde{\mathbf{v}},\widetilde{\theta},\widetilde{\mathbf{u}},\widetilde{\mathbf{g}})\in\mathcal{U}$.

Finally, due to the continuity and the convexity of the cost functional $J$ (in particular, $J$ is weakly lower semicontinuous), we deduce that:

[TABLE]

Thus, $\widetilde{\mathbf{g}}\in\mathcal{U}_{ad}$ is a solution of the optimal control problem $(\mathcal{P})$. ∎

Remark 5

It is worthwhile remarking here that, using standard techniques in the spirit of those presented in below section, it is possible to obtain a formal optimality system for the characterization of the optimal solutions of the control problem $(\mathcal{P})$. However, since this is not the main aim of this paper, and for the sake of brevity, we will not present here this optimality system, focusing our attention on the numerical computation of these optimal solutions. ∎

5 Numerical resolution of the control problem

In this section we will present a numerical approximation for the optimal control problem $(\mathcal{P})$. So, we will discretize the state systems (3), (6) and (10) using a standard finite element method, and we will compute the numerical approximation of the resulting nonlinear optimization problem (that appears after the full space-time discretization of the control problem) using an interior point algorithm. In this particular case, due to the specific relations between the dimensions of the control and the constraint variables, the numerical approximation of the Jacobian matrix of the constraints will be performed using the discretized adjoint system (row by row) instead of the linearized systems (column by column). In addition, the computation of each row of the Jacobian matrix will be parallelized.

5.1 Space-time discretization

Let us consider a regular partition $0=t_{0}<t_{1}<\ldots<t_{N}=T$ of the time interval $[0,T]$ such that $t_{n+1}-t_{n}=\Delta t=\frac{1}{\alpha}$, $\forall n=0,\ldots,N-1$, and recall the material derivative of a generic scalar field $\phi$ defined as:

[TABLE]

where $\mathbf{X}$ represents the characteristic line, that is, verifies the equation:

[TABLE]

So, we can approximate the material derivative in the following way:

[TABLE]

where $\phi^{n}$ represents an approximation to $\phi(t_{n})$, and $\mathbf{X}^{n}_{-}(\mathbf{x})=\mathbf{X}(\mathbf{x},t_{n+1};t_{n})$ (i.e., the position at time $t_{n}$ of a particle that at time $t_{n+1}$ was located at point $\mathbf{x}$) is the solution of the following trajectory equation:

[TABLE]

approached by the Euler scheme, that is, we consider the following approximation (see further details, for instance, in cita1 ; cita2 ):

[TABLE]

For the space discretization, we take a family of meshes $\tau_{h}$ for the domain $\Omega$ with characteristic size $h$ and, associated to this family of meshes, we define the following finite element spaces (cf. Section 4.1 of girault1 ):

• •

$\mathbf{V}_{h}$ ($\mathbb{P}_{1b}$ FEM space) for the water velocity $\mathbf{v}$:

[TABLE]

and, for the test functions and the adjoint state, the subspace:

[TABLE]

• •

$M_{h}$ ($\mathbb{P}_{1}$ FEM space) for the water pressure $p$:

[TABLE]

• •

$K_{h}$ ($\mathbb{P}_{1}$ FEM space) for the water temperature $\theta$:

[TABLE]

and, for the test functions and the adjoint state, the subspace:

[TABLE]

• •

$\mathbf{X}_{h}$ ($\mathbb{P}_{1}$ FEM space) for the concentration $\mathbf{u}$ of the species involved in eutrophication process:

[TABLE]

and, for the test functions and the adjoint state, the subspace:

[TABLE]

With respect to the computational treatment of the problem, we have used the open code FreeFem++ HECHT1 for the space-time discretizations of the problem. We have also employed a penalty method (cf. Section 4.3 of girault1 ) for computing the solution of the Stokes problems that appear after discretization. Finally, in order to reduce the CPU time necessary for computing the solution of the state systems, we have applied an explicit scheme (evaluations in previous time step) for the nonlinearities and the coupled terms of the discretized problem.

So, we consider the following space-time discretization for the optimal control problem $(\mathcal{P})$ where, for the sake of simplicity, we will use the same notations for the discrete problem as in the case of the continuous one:

Coupling of temperature/species in collectors and injectors:

We denote by $\theta^{n}\in K_{h}$ and $\mathbf{u}^{n}\in\mathbf{X}_{h}$, respectively, the water temperature and the species concentration at time step $n=0,\ldots,N$. Then, we consider the following approximation for functions $\gamma_{\theta}^{k}$, $k=1,\ldots,N_{CT}$, defined in (8), (analogously for functions $\gamma_{u^{i}}^{k}$, $k=1,\ldots,N_{CT}$, $i=1,\ldots,5$, defined in (12)):

[TABLE]

Moreover, if we assume the value $\epsilon=\frac{\Delta t}{2}$ in the definition (9) of function $\rho_{\epsilon}$ we have that the support of $\rho_{\Delta t/2}(t^{n}-\frac{\Delta t}{2}-s)$ is contained in $(t^{n}-\Delta t,t^{n})=(t^{n-1},t^{n})$, for all $n=1,\ldots,N$, and then:

[TABLE]

Finally, we approximate each element $\varphi^{k}$ by the indicator function of the injector $T^{k}$, $k=1,\ldots,N_{CT}$, and each element $\tilde{\varphi}^{k}$ by the indicator function of the collector $C^{k}$, $k=1,\ldots,N_{CT}$. Thus, the temperature in each injector at time step $t_{n}$ is the mean temperature in the corresponding collector at time step $t_{n-1}$ (analogously for the species of the eutrophication model). 2. 2.

Discretized control:

We consider the following discretization of the admisible set (18) (we will also denote by ${\mathcal{U}}_{ad}$ the set of admissible discrete controls):

[TABLE]

where $c_{1}\,,c_{2}>0$ technological bounds related to mechanical characteristics of pumps and they are chosen so that $\|g^{k}\|_{H^{1}(0,T)}\leq c$, $\forall k=1,\ldots,N_{CT}$. So, if we consider the standard basis for the previous finite element space, we can consider the following discrete control:

[TABLE] 3. 3.

Discretized cost functional:

In order to simplify the numerical resolution of the control problem, we will consider the following modification of the cost functional restriction to the previous admissible set:

[TABLE]

where $\sigma_{1}$ and $\sigma_{2}$ are positive weights that we will take into account in the numerical tests. 4. 4.

Discretized state constraints:

We consider the function:

[TABLE]

where, for each $n=1,\ldots,N$,

[TABLE]

with $\mathbf{u}^{n+1}\in\mathbf{Z}_{h}$ the solution of the discretized eutrophication model. Thus, we can express:

[TABLE]

It is worthwhile remarking here that, due the type of time discretization considered for the material derivatives (81), the control $\mathbf{g}^{N}$ acts over the species and the temperature at time $t_{N+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 Figure 3.

Water velocity and pressure:

Given $\mathbf{v}^{0}\in\mathbf{W}_{h}$, the pair velocity/pressure $(\mathbf{v}^{n+1},p^{n+1})\in\mathbf{V}_{h}\times M_{h}$, for each $n=0,1,\ldots,N-1$, with:

[TABLE]

is the solution of the fully discretized system:

[TABLE]

where $\lambda>0$ is the penalty parameter and $\beta(\mathbf{v}^{n})=2\nu+2\nu_{tur}[e(\mathbf{v}^{n}):e(\mathbf{v}^{n})]^{1/2}$. 6. 6.

Water temperature: Given $\theta^{0}\in K_{h}$, the temperature $\theta^{n+1}\in K_{h}$, for each $n=0,\ldots,N$, with

[TABLE]

is the solution of the fully discretized system:

[TABLE] 7. 7.

Eutrophication species concentration:

Given $\mathbf{u}^{0}\in\mathbf{X}_{h}$, the species concentration $\mathbf{u}^{n+1}\in\mathbf{X}_{h}$, for each $n=0,\ldots,N$, with:

[TABLE]

is the solution of the fully discretized system:

[TABLE]

where $\mathbf{A}^{n}(\theta^{n},\mathbf{u}^{n})\in\mathbb{R}^{5\times 5}$ is the following matrix:

[TABLE]

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:

[TABLE]

In order to solve this nonlinear optimization problem, we will use the interior point algorithm IPOPT Biegler1 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 $\delta_{\mathbf{g}}J(\mathbf{g})\in\mathcal{L}(\mathbb{R}^{N\times N_{CT}},\mathbb{R})$ is such that, for any $\delta\mathbf{g}=(\delta\mathbf{g}^{1},\ldots,\delta\mathbf{g}^{N})\in\mathbb{R}^{N\times N_{CT}}$:

[TABLE]

Therefore, $[\nabla_{\mathbf{g}}J(\mathbf{g})]_{i}=\delta_{\mathbf{g}}J(\mathbf{g})(\mathbf{e}_{i})$, where $\mathbf{e}_{i}$, $i=1,\ldots,N\times N_{CT}$, is the $i$-th vector of the canonical basis in $\mathbb{R}^{N\times N_{CT}}$.

In the case of the Jacobian matrix of the constraints, we know that the differential associated to the application $\mathbf{G}:\mathcal{U}_{ad}\subset\mathbb{R}^{N\times N_{CT}}\rightarrow\mathbb{R}^{N}$ is such that $\delta_{\mathbf{g}}\mathbf{G}(\mathbf{g})\in\mathcal{L}(\mathbb{R}^{N\times N_{CT}},\mathbb{R}^{N})$. So, given any element $\delta\mathbf{g}\in\mathbb{R}^{N\times N_{CT}}$, we have that $\delta_{\mathbf{g}}\mathbf{G}(\mathbf{g})(\delta\mathbf{g})\in\mathbb{R}^{N}$, and the Jacobian matrix $J_{\mathbf{g}}\mathbf{G}(\mathbf{g})\in\mathcal{M}_{N\times(N_{CT}\times N)}$ is such that $[J_{\mathbf{g}}\mathbf{G}(\mathbf{g})]_{j,i}=\langle\delta_{\mathbf{g}}\mathbf{G}(\mathbf{g})(\mathbf{e}_{i}),\widetilde{\mathbf{e}}_{j}\rangle$, where $\widetilde{\mathbf{e}}_{j}$, $j=1,\ldots,N$, is the $j$-th vector of the canonical basis in $\mathbb{R}^{N}$. 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\times N_{CT}$) and the dimension of the space where the application $\mathbf{G}$ takes values ($N$).

• •

When using the linearized systems, we would have to solve $N_{CT}\times N$ times these systems (in this case, we would compute the Jacobian matrix column by column):

[TABLE]

where $\delta_{\mathbf{g}}\mathbf{G}(\mathbf{g})(\mathbf{e}_{k})\in\mathcal{M}_{N\times 1}(\mathbb{R})$, for $k=1,\ldots,N_{CT}$.

• •

When employing the adjoint state systems, we would have to solve $N$ times these systems (now, we would compute the Jacobian row by row):

[TABLE]

where $\nabla_{\mathbf{g}}G^{i}(\mathbf{g})\in\mathbb{R}^{N_{CT}}$ is such that $[\nabla_{\mathbf{g}}G^{j}(\mathbf{g})]_{i}=\delta_{\mathbf{g}}G^{j}(\mathbf{g})(\mathbf{e}_{i})$, $j=1,\ldots,N$, $i=1,\ldots,N\times N_{CT}$.

In our case, $N_{CT}>1$. So, it is more advantageous to employ the adjoint state systems and computing the Jacobian matrix row by row ($j=1,\ldots,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 6 (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 $\delta\mathbf{g}\in\mathbb{R}^{N\times N_{CT}}$, then

[TABLE]

where $\{(\delta\mathbf{v}^{n},\delta p^{n})\}_{n=0}^{N}\subset\mathbf{V}_{h}\times M_{h}$, $\{\delta\theta^{n}\}_{n=0}^{N+1}\subset K_{h}$ and $\{\delta\mathbf{u}^{n}\}_{n=0}^{N+1}\subset\mathbf{X}_{h}$ are, 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 $(\delta\mathbf{v}^{0},\delta p^{0})=(\mathbf{0},0)$, for each $n=0,\ldots,N-1$, $(\delta\mathbf{v}^{n},\delta p^{n})\in\mathbf{V}_{h}\times M_{h}$, with

[TABLE]

is the solution of:

[TABLE]

where $\gamma(\mathbf{v}^{n})=2\nu_{tur}[e(\mathbf{v}^{n}):e(\mathbf{v}^{n})]^{-1/2}$.

• •

Linearized system for water temperature: Given $\delta\theta^{0}=0$, for each $n=0,\ldots,N$, $\delta\theta^{n+1}\in K_{h}$, with:

[TABLE]

is the solution of:

[TABLE]

• •

Linearized system for eutrophication model: Given $\delta\mathbf{u}^{0}=\mathbf{0}$, for each $n=0,\ldots,N$, $\delta\mathbf{u}^{n+1}\in\mathbf{X}_{h}$, with:

[TABLE]

is the solution of:

[TABLE]

Proof

The proof is straightforward, where the only drawback is related to the computation of terms of the type $\delta_{\mathbf{g}}\left(\boldsymbol{\varphi}(\mathbf{g},\mathbf{x}-\Delta t\,\mathbf{v}(\mathbf{g},\mathbf{x}))\right)(\delta g)$, where $\boldsymbol{\varphi}(\mathbf{g},\mathbf{x})$ and $\mathbf{v}(\mathbf{g},\mathbf{x})$ are vector functions smooth enough (the scalar case would be analogous). Nevertheless, using the chain rule, we can easily obtain that:

[TABLE]

(We must note here that, in our specific formulation, we deal with the function $b:x\in\mathbb{R}\rightarrow b(x)=x\,|x|^{3}$, that is differentiable in $\mathbb{R}$, with $b^{\prime}(x)=4\,|x|^{3}$). ∎

Lemma 7 (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,\ldots,N$, the matrices $\{\nabla_{\mathbf{g}^{n}}G^{k}(\mathbf{g})\}_{n=1}^{k}\subset\mathcal{M}_{1\times N_{CT}}(\mathbb{R})$ can be computed using the following expressions:

• •

If $n\in\{1,\ldots,k\}\setminus\{N\}$,

[TABLE]

• •

If $n=N$,

[TABLE]

where if we introduce, for each row $k=1,\ldots,N$, the following vector (defined from the usual Kronecker delta $\delta_{ij}$ and the indicator function of subset $\Omega_{C}$):

[TABLE]

then the adjoint states associated to the eutrophication system $\{\mathbf{z}^{n}\}_{n=0}^{N+1}\subset\mathbf{Z}_{h}$, to the hydrodynamic system $\{(\mathbf{w}^{n},q^{n})\}_{n=0}^{N}\subset\mathbf{W}_{h}\times M_{h}$, and to the temperature system $\{\xi^{n}\}_{n=0}^{N+1}\subset H_{h}$ are, respectively, the solution of the following systems:

• •

Adjoint system for eutrophication model:

• –

For $n=N+1$, $\mathbf{z}^{n}=\mathbf{0}$.

• –

For $n=N$, $\mathbf{z}^{n}\in\mathbf{Z}_{h}$ is such that:

[TABLE]

where $\mathbf{X}_{+}^{n+1}(\mathbf{x})=\mathbf{x}+\Delta t\,\mathbf{v}^{n+1}$.

• –

For $n=N-1,\ldots,0$, $\mathbf{z}^{n}\in\mathbf{Z}_{h}$ is such that:

[TABLE]

• •

Adjoint system for water temperature:

• –

For $n=N+1$, $\xi^{n}=0$.

• –

For $n=N$, $\xi^{n}\in H_{h}$ is such that:

[TABLE]

• –

For $n=N-1,\ldots,0$, $\xi^{n}\in H_{h}$ is such that:

[TABLE]

• •

Adjoint system for water velocity and pressure:

• –

For $n=N$, $(\mathbf{w}^{n},q^{0})=(\mathbf{0},0)$.

• –

For $n=N-1$, $(\mathbf{w}^{n},q^{n})\in\mathbf{W}_{h}\times M_{h}$ is such that:

[TABLE]

• –

For $n=N-2,\ldots,0$, $(\mathbf{w}^{n},q^{n})\in\mathbf{W}_{h}\times M_{h}$ is such that:

[TABLE]

Remark 6

In order to simplify the proof of above Lemma, we have established the adjoint systems (113)-(114), (111)-(112) and (109)-(110) in a strong formulation (contrary to the case of the linearized systems (104), (106) and (108), where we have proposed a variational formulation). It is also clear that these adjoint systems easily admits 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 (104), (106) and (108), respectively, the $n$-th component of the sequences $\{(\mathbf{w}^{n},q^{n})\}_{n=0}^{N}\subset\mathbf{W}_{h}\times M_{h}$, $\{\xi^{n}\}_{n=0}^{N+1}\subset H_{h}$ and $\{\mathbf{z}^{n}\}_{n=0}^{N}\subset\mathbf{Z}_{h}$, such that $\mathbf{w}^{N}=\mathbf{0}$, $q^{N}=0$, $\xi^{N+1}=0$ and $\mathbf{z}^{N+1}=\mathbf{0}$, and let us sum in $n$ from [math] 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:

• •

For eutrophication model:

[TABLE]

with $\mathbf{X}_{+}^{n+1}(\mathbf{x})=\mathbf{x}+\Delta t\,\mathbf{v}^{n+1}$, and where we are assuming $\delta\mathbf{v}^{N+1}=\mathbf{0}$ in order to simplify the notation.

• •

For water temperature:

[TABLE]

where, for the sake of simplicity, we have also assumed $\delta\mathbf{v}^{N+1}=\mathbf{0}$.

• •

For water velocity:

[TABLE]

where we have assumed $\mathbf{w}^{N+1}=\mathbf{0}$.

Thus, if we define $\{(\mathbf{w}^{n},q^{n})\}_{n=0}^{N}\subset\mathbf{W}_{h}\times M_{h}$, $\{\xi^{n}\}_{n=0}^{N+1}\subset H_{h}$ and $\{\mathbf{z}^{n}\}_{n=0}^{N+1}\subset\mathbf{Z}_{h}$, such that $\mathbf{w}^{N}=\mathbf{0}$, $q^{N}=0$, $\xi^{N+1}=0$ and $\mathbf{z}^{N+1}=\mathbf{0}$, as the solutions of the adjoint system (113)-(114), (111)-(112) and (109)-(110), respectively, we obtain, after summing above expressions (115), (116) and (117), that:

[TABLE]

And, finally, from the definition:

[TABLE]

∎

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 $\Omega$. So, we consider a space configuration similar to that presented in Figure 1, with $N_{CT}=4$ collector/injector pairs, in a rectangular domain of $20\,{\rm m}\times 16\,{\rm m}$. We suppose that the diameter of each collector is $1\,{\rm m}$ and the diameter of each injector is $2\,{\rm m}$. For the coefficients of the eutrophication model (10), we have used the same values as those appearing in DRAGO200117 , and for the thermo-hydrodynamic system (3), (6) we have employed the same values as in fran8 . For the space discretization we have generated a regular mesh of $2989$ vertices, as shown in Figure 4.

The control domain $\Omega_{C}$ corresponds to a $3\,{\rm m}$ strip at the bottom of the domain, and all the numerical tests have been performed in a temporal horizon of 12 hours ($T=43200\,{\rm s}$). Finally, in order to simulate the effects of solar radiation for the heat equation (6), we consider the standard function $T_{r}$ depicted in Figure 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 $\widetilde{\mathbf{g}}$, with $\widetilde{g}^{n,k}=C$ (constant), for $n=1,\ldots,N$, $k=1,\ldots,N_{CT}$, we will solve the following modification of the original optimization problem $(\mathcal{P})$:

[TABLE]

In other words, we want to find an optimal control $\widehat{\mathbf{g}}\in\mathcal{U}_{ad}$ that supplies us with a higher concentration of dissolved oxygen than that obtained with the constant control $\widetilde{\mathbf{g}}$, and that minimizes the energy cost functional $J$. As an illustration to this behaviour, in Figure 6 we can see the evolution of the mean concentration of dissolved oxygen in the control domain $\Omega_{C}$ considering a constant reference control $\widetilde{g}^{n,k}=1.0\times 10^{-4}\,{\rm m}^{3}\,{\rm s}^{-1}$, $\forall n=1,\ldots,N$, $\forall k=1,\ldots,N_{CT}$, compared to the mean concentration assuming that all the pumps are out of service (that is, $\widetilde{g}^{n,k}=0.0\,{\rm m}^{3}\,{\rm s}^{-1}$, $\forall n=1,\ldots,N$, $\forall k=1,\ldots,N_{CT}$). 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 $\Delta t$. We must mention that in the numerous numerical tests developed, we have always obtained that $\mathbf{G}(\widehat{\mathbf{g}})=\mathbf{G}(\widetilde{\mathbf{g}})$, and also a reduction in the value of the cost functional $J(\widehat{\mathbf{g}})<J(\widetilde{\mathbf{g}})$. So, in Figure 7 we can see the optimal control that we have obtained taking $\sigma_{1}=0.5$ and $\sigma_{2}=1-\sigma_{1}=0.5$, for time steps of $\Delta t=3600\,{\rm s}$ and $\Delta t=1800\,{\rm s}$ (corresponding to $N=12$ and $N=24$, respectively). In Figure 8 we can find the optimal control corresponding to time steps of $\Delta t=900\,{\rm s}$ and $\Delta t=450\,{\rm s}$ ($N=48$ and $N=96$, respectively), showing the robustness of our methodology.

We observe that the flow rates associated to the two upper collectors ($g^{1}$ and $g^{3}$) are significantly higher than the corresponding to lower collectors ($g^{2}$ and $g^{4}$). 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. This is because as we decrease the time step, we can act more precisely over the system and achieve better results.

In Figure 9 we can see the evolution of the constraints for the choice of the time step length $\Delta t=450\,{\rm s}$. We can verify there that the optimal constraint $\mathbf{G}(\widetilde{\mathbf{g}})$ and the reference constraint $\mathbf{G}(\widehat{\mathbf{g}})$ are virtually indistinguishable, that is, with optimal strategy $\widetilde{\mathbf{g}}$ we obtain the same water quality in the control region as with the constant reference flow rate $\widetilde{\mathbf{g}}$, but with a significative decrease in energy cost.

Finally, in Figure 10 we show the concentration of dissolved oxygen in the whole domain $\Omega$ associated to the optimal control solution for $\Delta t=450\,{\rm 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.