## 1 Introduction

Usually, wastewater originated from urban areas or industry undergoes a physical/chemical/biological treatment in a purifying plant. From there, wastewater is discharged through a submarine outfall into an aquatic media like a lake, a river or a coastal area, at a suitable distance from protected areas like beaches, fisheries or marine recreation zones. When several purifying plants are going to discharge wastewater into the same domain, the problem of design and management of the whole treatment system arises. To solve this problem, optimal control theory and optimization methods can be very useful: they can help decision makers in formulating rational policies in order to minimize costs while keeping the prescribed levels of water quality. Sections 4 and 5 are devoted to present two typical examples: the optimal management of a wastewater treatment system and the optimal location of the submarine outfalls.

Obviously, the knowledge of mathematical models for the evolution of pollutant concentration is a unavoidable first step if one wants to use optimal control techniques. So, the first part of this work is devoted to study in detail one of them, related with Biochemical Oxygen Demand (BOD) and Dissolved Oxygen (DO), which is frequently used in the case of domestic discharges. In Section 2 we deal with the complete two-dimensional model consisting of two uncoupled system of partial differential equations which give us the height of water, the depth-averaged horizontal velocity of water and the depth-averaged concentration of BOD and DO. In Section 3 we obtain a new zero-dimensional model, by integrating the previous one, which can be very useful to study the global performance of the system. We apply this model in the ría of Arousa (Spain) to obtain the global concentration of BOD and DO along the time.

## 2 A complete two-dimensional model for Biochemical Oxygen Demand and Dissolved Oxygen

We consider a domain Ω occupied by shallow water (as can be a ría or an estuary), where polluting wastewater are discharged through ${N}_{\mathrm{E}}$ submarine outfalls located at points ${b}_{j}\in \Omega $ and connected to their respective purifying plants located at points ${a}_{j}\in \Gamma $ (see Fig. 1). Moreover, we assume the existence of several areas ${A}_{i}\subset \Omega ,\phantom{\rule{0.25em}{0ex}}i=1,\dots ,{N}_{\mathrm{Z}}$, representing fisheries, beaches or marine recreation zones where it is necessary to guarantee the water quality with pollution levels lower than some allowed threshold levels.

The flow of an effluent from a submarine outfall in a shallow water domain is mainly governed by horizontal transport due to currents (produced by tides, wind …) and turbulent diffusion. It allows uncouple hydrodynamical equations and transport equations. The former give the height and the velocity field which are used in latter to obtain the pollution levels.

### 2.1 Hydrodynamic model: the shallow water equations

In this section we recall the shallow water equations which constitute an useful mathematical model for hydrodynamic flows in shallow regions. If we consider Γ, the boundary of Ω, divided into three parts: ${\Gamma}^{-}$ (corresponding to the effluent), ${\Gamma}^{+}$ (corresponding to open sea) and ${\Gamma}^{0}$ (corresponding to the cost), then the shallow water equations can be written in the following way:

$$\begin{array}{cc}\frac{\partial h}{\partial t}+\overrightarrow{\mathrm{\nabla}}\cdot (h\overrightarrow{u})=0\text{\hspace{1em}}\hfill & \text{in}\phantom{\rule{0.25em}{0ex}}\Omega \times (0,T)\hfill \\ \multicolumn{2}{c}{\frac{\partial \overrightarrow{u}}{\partial t}+(\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla}})\overrightarrow{u}-\nu \mathrm{\Delta}\overrightarrow{u}+g\overrightarrow{\mathrm{\nabla}}h=\overrightarrow{F}}\\ \text{\hspace{1em}}\hfill & \text{in}\phantom{\rule{0.25em}{0ex}}\Omega \times (0,T)\hfill \\ h=\eta \text{\hspace{1em}}\hfill & \text{on}\phantom{\rule{0.25em}{0ex}}{\Gamma}^{-}\times (0,T)\hfill \\ h=\varphi \text{\hspace{1em}}\hfill & \text{on}\phantom{\rule{0.25em}{0ex}}{\Gamma}^{+}\times (0,T)\hfill \\ h(0)={h}_{0}\text{\hspace{1em}}\hfill & \text{in}\phantom{\rule{0.25em}{0ex}}\Omega \hfill \\ \overrightarrow{u}\cdot \overrightarrow{n}=q\text{\hspace{1em}}\hfill & \text{on}\phantom{\rule{0.25em}{0ex}}{\Gamma}^{-}\times (0,T)\hfill \\ \overrightarrow{u}\cdot \overrightarrow{n}=0\text{\hspace{1em}}\hfill & \text{on}\phantom{\rule{0.25em}{0ex}}{\Gamma}^{0}\times (0,T)\hfill \\ \overrightarrow{\mathrm{\nabla}}\cdot \overrightarrow{u}=0\text{\hspace{1em}}\hfill & \text{on}\phantom{\rule{0.25em}{0ex}}{\Gamma}^{+}\times (0,T)\hfill \\ \overrightarrow{u}(0)={\overrightarrow{u}}_{0}\text{\hspace{1em}}\hfill & \text{in}\phantom{\rule{0.25em}{0ex}}\Omega \hfill \end{array}\}$$ | (1) |

There are several works related to the existence, uniqueness and regularity of solution of (1) in particular cases (see, for instance, [1–5]) but in the general case it is still an open problem. A numerical approximation of h and $\overrightarrow{u}$ in $\Omega \times (0,T)$ can be obtained by finite difference, finite element or finite volume methods (see, for instance, [6–9]) in order to be employed in the pollutant dispersion model.

### 2.2 Pollutant dispersion: the BOD-DO model

Firstly, in order to simulate the water quality in Ω, we have to choose some indicators of pollution levels. Two of the most important (especially in the case of domestic discharges) are the Dissolved Oxygen (DO) and the organic matter, which can be measured in terms of the need of oxygen to decompose it, the so-called Biochemical Oxygen Demand (BOD). If the pollution level is not too high the BOD can be satisfied by the DO. However, if the organic matter increases beyond a maximum value the DO is not enough for its decomposition, leading to important modifications (anaerobic processes) in the ecosystem. To avoid them a threshold value of BOD may not be exceeded and a minimum level of DO must be guaranteed.

The evolution of the BOD and the DO in the domain $\Omega \subset {\mathbb{R}}^{2}$ is governed by a system of partial differential equations (cf. [10]). Let us denote by ${\rho}_{1}(x,t)$ and ${\rho}_{2}(x,t)$ the concentrations of BOD and DO at point $x\in \Omega $ and at time $t\in [0,T]$, respectively. Then these concentrations are obtained as the solution of the following two initial-boundary value problems:

$$\begin{array}{c}\frac{\partial {\rho}_{1}}{\partial t}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla}}{\rho}_{1}-{\beta}_{1}\mathrm{\Delta}{\rho}_{1}\text{\hspace{1em}}\hfill \\ \phantom{\rule{1em}{0ex}}=-{\kappa}_{1}{\rho}_{1}+\frac{1}{h}\sum _{j=1}^{{N}_{\mathrm{E}}}{m}_{j}\delta (x-{b}_{j})\text{\hspace{1em}}\hfill & \text{in}\phantom{\rule{0.25em}{0ex}}\Omega \times (0,T)\phantom{\rule{0.25em}{0ex}}\hfill \\ \frac{\partial {\rho}_{1}}{\partial n}=0\text{\hspace{1em}}\hfill & \text{on}\phantom{\rule{0.25em}{0ex}}\Gamma \times (0,T)\hfill \\ {\rho}_{1}(x,0)={\rho}_{10}(x)\text{\hspace{1em}}\hfill & \text{in}\phantom{\rule{0.25em}{0ex}}\Omega \hfill \\ \frac{\partial {\rho}_{2}}{\partial t}+\overrightarrow{u}\cdot \overrightarrow{\mathrm{\nabla}}{\rho}_{2}-{\beta}_{2}\mathrm{\Delta}{\rho}_{2}\text{\hspace{1em}}\hfill \\ \phantom{\rule{1em}{0ex}}=-{\kappa}_{1}{\rho}_{1}+\frac{1}{h}{\kappa}_{2}({d}_{\mathrm{s}}-{\rho}_{2})\text{\hspace{1em}}\hfill & \text{in}\phantom{\rule{0.25em}{0ex}}\Omega \times (0,T)\hfill \\ \frac{\partial {\rho}_{2}}{\partial n}=0\text{\hspace{1em}}\hfill & \text{on}\phantom{\rule{0.25em}{0ex}}\Gamma \times (0,T)\hfill \\ {\rho}_{2}(x,0)={\rho}_{20}(x)\text{\hspace{1em}}\hfill & \text{in}\phantom{\rule{0.25em}{0ex}}\Omega \hfill \end{array}\}$$ | (2) |

For this system we can prove the existence and uniqueness of solution (cf. [11]). A numerical solution for (2) can be obtained by a method which combines characteristics for time discretization with finite elements for space discretization (cf. [12]).

## 3 The global performance: a zero-dimensional model

The main aim in this paper is to study the global performance of the BOD-DO model in a ría, an estuary or a lake. In order to do it we integrate the system (2) in the region occupied by water. Thus, we obtain a new simple zero-dimensional model.

Let ${M}_{1}(t)$ and ${M}_{2}(t)$ be, respectively, the total mass of BOD and DO in the region at time t, that is to say,

$${M}_{1}(t)=\underset{\Omega}{\int}h(x,t){\rho}_{1}(x,t)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x$$ |

$${M}_{2}(t)=\underset{\Omega}{\int}h(x,t){\rho}_{2}(x,t)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x$$ |

$$\{\begin{array}{c}\frac{\mathrm{d}{M}_{1}}{\mathrm{d}t}=-{\kappa}_{1}{M}_{1}+Q+{F}_{1}\hfill \\ \frac{\mathrm{d}{M}_{2}}{\mathrm{d}t}=-{\kappa}_{1}{M}_{1}+A{\kappa}_{2}({d}_{\mathrm{s}}-\frac{{M}_{2}}{V})+{F}_{2}\hfill \end{array}$$ | (3) |

If we study a stationary region (as can be a closed lake) A and V are constant along the time. In this case, if we also suppose that Q, ${F}_{1}$ and ${F}_{2}$ are constant, the exact solution of (3) is given by:

(4) |

In order to develop an asymptotic study, for $t\gg \frac{1}{{\kappa}_{1}}$ and $t\gg \frac{1}{a}$, the terms with ${\mathrm{e}}^{-{\kappa}_{1}t}$ and ${\mathrm{e}}^{-at}$ can be negligible in (4) and the asymptotic values of BOD and DO are ${M}_{1}(\infty )=\raisebox{1ex}{$(Q+{F}_{1})$}\!\left/ \!\raisebox{-1ex}{${\kappa}_{1}$}\right.$ and ${M}_{2}(\infty )=\raisebox{1ex}{$(b-Q-{F}_{1})$}\!\left/ \!\raisebox{-1ex}{$a$}\right.$. However, it is as well to point out that, in a shallow water domain, the typical values of ${\kappa}_{1}$ and a are, respectively, of the order of 10^{−5} and ${10}^{\mathrm{-13}}\text{}{\text{s}}^{\mathrm{-1}}$, and then, the previous asymptotic values would be reached 300,000 years later!

For $\frac{1}{{\kappa}_{1}}\ll t\ll \frac{1}{a}$ (only about some days later), we can neglect ${\mathrm{e}}^{-{\kappa}_{1}t}$ and approximate ${\mathrm{e}}^{-at}\simeq 1-at$. Then we obtain an asymptotic value for the BOD and a linear approximation for the DO:

$${M}_{1}(t)\simeq \frac{Q+{F}_{1}}{{\kappa}_{1}}$$ |

$${M}_{2}(t)\simeq ({M}_{2}(0)+\frac{Q+{F}_{1}-{\kappa}_{1}{M}_{1}(0)}{{\kappa}_{1}-a})+(b-Q-{F}_{1}-a({M}_{2}(0)+\frac{Q+{F}_{1}-{\kappa}_{1}{M}_{1}(0)}{{\kappa}_{1}-a}))t$$ |

When we are working in an open region (as can be a ría or an estuary), we must take into account the water renewal. For example, in a ría, we can assume that inflow water is more pure than outflow water (in the sense that DO concentration is greater and BOD concentration is lesser) and, moreover, we know that DO concentration in inflow water is always lesser than oxygen saturation density (${d}_{\mathrm{s}}$). In this case, if $\varphi (t)$ denotes the height over a fixed level (half tide) at the mouth of the ría and $\tilde{V}$ denotes the water volume corresponding to half tide, then, by geometrical reasons, we have

$$V(t)=A\varphi (t)+\tilde{V}$$ |

Moreover, taking into account the purification of the water in the ría because of the seawater renewal, we obtain the following expressions for the mass flow rates across the boundary:

$${F}_{1}(t)={\varphi}^{\prime}(t)A\frac{\alpha {M}_{1}(t)}{V(t)}$$ |

$${F}_{2}(t)={\varphi}^{\prime}(t)AC(t)$$ |

$$\{\begin{array}{c}\frac{\mathrm{d}{M}_{1}(t)}{\mathrm{d}t}=-{\kappa}_{1}{M}_{1}(t)+Q+{\varphi}^{\prime}(t)A\frac{\alpha {M}_{1}(t)}{V(t)}\hfill \\ \frac{\mathrm{d}{M}_{2}(t)}{\mathrm{d}t}=-{\kappa}_{1}{M}_{1}(t)+A{\kappa}_{2}({d}_{\mathrm{s}}-\frac{{M}_{2}(t)}{V(t)})\hfill \\ +{\varphi}^{\prime}(t)AC(t)\hfill \end{array}$$ | (5) |

This system can be easily solved by finite difference methods. We have used the backward Euler scheme to solve it for data corresponding to the ría of Arousa (Spain). Moreover, in order to study the weight of the water renewal in the ría, we have compared the result with the exact solution of (3) with $V(t)=V$ and ${F}_{1}={F}_{2}=0$ (in this case the exact solution is given by (4)). As one can see in Figs. 2 and 3, for data given in Table 1, the BOD concentration reaches an asymptotic value in both cases, but the DO concentration observes a linear decrease without water renewal and reaches an asymptotic value (very close to saturation) when the water renewal is considered.

– Total time | T=4⋅10^{6} s |

– Initial mass of BOD | ${M}_{1}(0)=0\text{}\text{kg}$ |

– Initial mass of DO | ${M}_{2}(0)=1.16311\cdot {10}^{7}\text{}\text{kg}$ |

– Mass flow rate of BOD discharged in the region | $Q=0.256\text{}\text{kg}\phantom{\rule{0.2em}{0ex}}{\text{s}}^{\mathrm{-1}}$ |

– Oxygen saturation density | ${d}_{\mathrm{s}}=8.98\cdot {10}^{\mathrm{-3}}\text{}\text{kg}\phantom{\rule{0.2em}{0ex}}{\text{m}}^{\mathrm{-3}}$ |

– Kinetic coefficients ${\kappa}_{1}$ and ${\kappa}_{2}$ | ${\kappa}_{1}=1.15\cdot {10}^{\mathrm{-5}}\text{}{\text{s}}^{\mathrm{-1}}$ |

${\kappa}_{2}=9\cdot {10}^{\mathrm{-12}}\text{}\text{m}\phantom{\rule{0.2em}{0ex}}{\text{s}}^{\mathrm{-1}}$ | |

– Area of water surface in contact with air | $A=6.77399\cdot {10}^{7}\text{}{\text{m}}^{2}$ |

– Water volume corresponding to half tide | $V=1.43897\cdot {10}^{9}\text{}{\text{m}}^{3}$ |

Exclusive data for model (5) | |

– Height over a fixed level at the mouth of the ría | $\varphi (t)=1.4(\mathrm{sin}\left(\frac{8\pi t}{178560-\pi}\right)+1)$ |

– Parameters α and β if ${\varphi}^{\prime}(t)<0$ | α=1, β=1 |

– Parameters α and β if ${\varphi}^{\prime}(t)>0$ | α=0.95, β=1.02 |

## 4 Problem 1: Optimal management of a wastewater treatment system

In this section we pose an optimal control problem related to the management of a wastewater treatment system. We recall the situation is Section 2: we consider a domain occupied by shallow water where polluting wastewater is discharged through ${N}_{\mathrm{E}}$ submarine outfalls (in this section we suppose that these outfalls are located at fixed points ${b}_{j}\in \Omega $). Moreover, there exist several sensitive areas, ${A}_{i}\in \Omega $, in the domain, where it is necessary to guarantee the water quality with pollution levels lower than an allowed threshold. We suppose that wastewater arrives to the purifying plants with a certain BOD concentration. Before discharging it into the sea, its BOD concentration can be reduced in the plants by different biological or biochemical treatments. From the ecological point of view the depuration in each plant must be as high as possible but, from the economical point of view, there is a cost proportional to the developed depuration. Then, the optimal management problem is determining the depuration at each plant along the time, in such a way that the global depuration cost is minimized and the above constraints on the water quality are satisfied.

### 4.1 Mathematical formulation

In order to formulate this problem we need to take into account some issues. Firstly, if ${\overline{m}}_{j}$ denotes the BOD of wastewater arriving to the jth plant and ${\underline{m}}_{j}$ is the BOD corresponding to the maximum depuration at that plant, then determining the depuration at the jth plant is equivalent to finding the mass flow rate of BOD, ${m}_{j}(t)$, discharged through the corresponding outfall. We assume that they satisfy the constraints

$$\underline{m}\u2a7d{m}_{j}(t)\u2a7d\overline{m},\phantom{\rule{1em}{0ex}}j=1,2,\dots ,{N}_{\mathrm{E}}$$ | (6) |

$$\{\begin{array}{cc}{\rho}_{1}{|}_{{A}_{i}}\u2a7d\sigma \hfill & i=1,\dots ,{N}_{\mathrm{Z}}\hfill \\ {\rho}_{2}{|}_{{A}_{i}}\u2a7e\zeta \hfill & i=1,\dots ,{N}_{\mathrm{Z}}\hfill \end{array}$$ | (7) |

Finally, we suppose that the cost of the depuration process at the jth plant is known and it is a strictly convex ${C}^{2}$-function of the BOD discharged through the corresponding outfall. Hence, if ${f}_{j}$ denotes the cost function at jth plant, the cost of the whole depuration is given by,

$${J}_{1}(m)=\sum _{j=1}^{{N}_{\mathrm{E}}}\underset{0}{\overset{T}{\int}}{f}_{j}({m}_{j}(t))\phantom{\rule{0.2em}{0ex}}\mathrm{d}t$$ | (8) |

This is an optimal control problem with pointwise state constraints and with pointwise control. We can obtain the existence and uniqueness of solution and an optimality system in order to characterize it (see [11] for further details).

### 4.2 Numerical resolution

The numerical solution of the optimal control problem $({\mathcal{P}}_{1})$ requires a discretization of the state system (2). Moreover, in practice, because of the particular shape of the functions ${f}_{j}$, the constraints ${m}_{j}(t)\u2a7d\overline{m}$ can be suppressed. Thus, we are led to consider the discretized constraints:

$${g}_{1}\phantom{\rule{0.2em}{0ex}}:m\to {g}_{1}(m)=({\rho}_{1}-\sigma ,\phantom{\rule{0.25em}{0ex}}\zeta -{\rho}_{2})$$ |

$${g}_{2}\phantom{\rule{0.2em}{0ex}}:m\to {g}_{2}(m)=\underline{m}-m$$ |

Then, the optimal control problem $({\mathcal{P}}_{1})$ can be written as follows,

$$({\mathcal{P}}_{1\mathrm{D}})\phantom{\rule{0.25em}{0ex}}\{\begin{array}{c}\mathrm{min}{J}_{1}(m)\hfill \\ \text{such that}{g}_{i}(m)\u2a7d0,\hfill & i=1,2\hfill \end{array}$$ |

### 4.3 Numerical results

In this work we resolve this problem in a realistic situation posed in the ría of A Coruña (Spain). We consider one purifying plant and one protected area (see Fig. 4). We assume pollutant concentration of the wastewater arriving to the purifying plant is $150\text{}\text{kg}\phantom{\rule{0.2em}{0ex}}{\text{m}}^{\mathrm{-3}}$ (then the depuration cost above this value is constant and equal to 100) and we suppose that the complete depuration is not possible. Thus, we take the following cost function

$${f}_{1}(x)=\{\begin{array}{cc}\frac{100\times {150}^{3}}{{x}^{3}-3\times 150{x}^{2}+3\times {150}^{2}x}\hfill & \text{if}\phantom{\rule{0.25em}{0ex}}x\u2a7d150\hfill \\ 100\hfill & \text{if}\phantom{\rule{0.25em}{0ex}}x\u2a7e150\hfill \end{array}$$ |

In order to make sure that the previous interpretation is correct we have found the optimal discharge during several tidal cycles. In Fig. 6, we show the optimal discharge during three cycles. In effect, we observe some periodicity, however, because of BOD accumulation, the discharge is greater in the next cycle than in the previous one. The DO concentration at the last point of time, when the saturation of the constraints takes place, can be seen in Fig. 7.

## 5 Problem 2: Optimal location of wastewater outfalls

The second problem is connected with the optimal design of a wastewater treatment system. Particularly, it consists of finding the optimal location of the submarine outfalls.

We consider a similar situation to the previous one: a domain occupied by shallow water where we are going to discharge polluting wastewater through submarine outfalls and where there exist several sensitive areas in which we have to guarantee the water quality in terms of BOD and DO. Moreover, we also suppose that there are ${N}_{\mathrm{E}}$ purifying plants (located at points ${a}_{j}\in \Gamma $) but, unlike the previous problem, we now assume that the depuration in every plant is fixed (the functions ${m}_{j}(t)$ are known beforehand) and our goal is to determine the points ${b}_{j}$ where wastewater will be discharged. These points must be determined in order to minimize the construction cost of the submarine outfalls while guaranteeing the water quality at the protected areas.

### 5.1 Mathematical formulation

The constraints on the water quality are given in (7). Moreover, taking into account technological limitations, the jth outfall must be placed in a suitable region ${U}_{j}$, where ${U}_{j}\subset \Omega \backslash {\bigcup}_{i=1}^{{N}_{\mathrm{Z}}}{\overline{A}}_{i}$ is a compact convex polyhedral set representing all the admissible points where outfalls can be located. Thus, the optimal locations must verify ${b}_{j}\in {U}_{j},\phantom{\rule{0.25em}{0ex}}\forall j=1,\dots ,{N}_{\mathrm{E}}$. If we define ${U}_{\mathrm{ad}}={\prod}_{j=1}^{{N}_{\mathrm{E}}}{U}_{j}$, this constraint can be written in the simpler way

$$b\in {U}_{\mathrm{ad}}$$ | (9) |

$${J}_{2}(b)=\sum _{j=1}^{{N}_{\mathrm{E}}}\frac{1}{2}\Vert {b}_{j}-{a}_{j}{\Vert}^{2}$$ | (10) |

This is a control problem with quadratic cost but with nonconvex pointwise state constraints which makes difficult its analysis and resolution. A complete theoretical analysis of this problem can be seen in [14].

### 5.2 Numerical resolution

Now, in order to solve the problem $({\mathcal{P}}_{2})$, we introduce a discretization of the control problem in the same way that the previous problem. Firstly, the function collecting the discretized state constraints is denoted by ${\overline{g}}_{1}$,

$${\overline{g}}_{1}\phantom{\rule{0.2em}{0ex}}:b\to {\overline{g}}_{1}(b)=({\rho}_{1}-\sigma ,\phantom{\rule{0.25em}{0ex}}\zeta -{\rho}_{2})$$ |

Then the optimal control problem $({\mathcal{P}}_{2})$ is approximated by the following discrete optimization problem,

$$({\mathcal{P}}_{2\mathrm{D}})\phantom{\rule{0.25em}{0ex}}\{\begin{array}{c}\mathrm{min}{J}_{2}(b)\hfill \\ \text{such that}{\overline{g}}_{i}(b)\u2a7d0,\hfill & i=1,2\hfill \end{array}$$ |

In [15], three different algorithms are used to obtain the numerical solution of $({\mathcal{P}}_{2\mathrm{D}})$ namely an admissible points algorithm, the Nelder–Mead simplex method and a duality method. As we can see in that paper, due to the geometric nature of the problem, the three algorithms present a good performance, specially the Nelder–Mead method.

### 5.3 Numerical results

In this section we present the numerical results obtained when solving the problem $({\mathcal{P}}_{2})$ for a realistic situation posed in the ría of Vigo (Spain) during a complete tidal cycle. We have taken two purifying plants, located near the coast at points ${a}_{1}=(0,11000)$ and ${a}_{2}=(6630,7200)$, and two protected areas (see Fig. 8).

The state constraints for both protected areas corresponds to ${\sigma}_{1}=0.0002$, ${\sigma}_{2}=0.000135$, ${\zeta}_{1}=0.008067$, ${\zeta}_{2}=0.0000805$. The admissible set ${U}_{\mathrm{ad}}$ and the optimal locations ${b}_{1}=(69,10224)$ and ${b}_{2}=(5535,8278)$, given by the Nelder–Mead method can be seen in Fig. 8. Moreover, this figure shows the BOD concentration at high tide, at the end of the tidal cycle that we have simulated.

## 6 Conclusions

In this work, mathematical modelling and optimal control theory have been successfully applied to an interesting ecological problem: the management and design of a wastewater treatment system for coastal areas. In the last sections of the paper we present numerical results for two different cases related to this ecological problem.

However, the main contribution of this paper is the obtaining of a zero-dimensional model, which can be used in the study of the global performance of the system; we present an application of this simple model in a realistic situation posed in the ría of Arousa (Spain).

Obtained results indicate the good performance of our models and show how mathematical tools can be very useful in the study of a wide range of ecological problems.

## Acknowledgments

This work was supported by Project BFM2003-00373 of Ministerio de Ciencia y Tecnología (Spain).