Optimizing the Design of an Estuarine Water Quality Monitoring Network by Optimal Control Techniques

In this work, we propose a novel methodology in order to automatically optimize the location of the sampling points for a water quality monitoring network in an estuary, in such a way that any unknown pollution source can be identified (both in intensity and location) from the data supplied by those sampling points. In the central part of the article, after a rigorous mathematical formulation of the environmental problem, the full details of its numerical implementation are given. Finally, we present and analyze the results when applying the above proposed technique to study a real case in Ría of Vigo (northwestern Spain).


Introduction
Environmental monitoring of estuarine waterbodies is a fundamental tool to assure the fulfillment of water quality standards in these ecosystems. Data obtained from the monitoring network -consisting of a set of measures of different chosen parameters, as pollutants or nutrients concentrations, salinity, temperature, pH, chlorophyll, dissolved oxygen, and so on -can be used for general management and/or restoration of water quality in these areas.
One of most important issues in the design of a monitoring network is the number and the location of the sampling stations. Their number is usually limited by the available budget, but the determining of the monitoring locations -in past times mainly fixed by the intuitive experiences of stakeholders and decision-makers -needs to be systematically and scientifically chosen in order to optimize the effective performance of the network. Scientific studies show that the accuracy of identifying pollution sources is highly dependent on the location of these monitoring stations [1]. Therefore, finding a set of optimal locations for the set of sampling points is essential to correctly characterize pollution sources (wastewater discharges, accidental spills, runoffs, etc.).
In fact, as it has been remarked by several authors [1,2], the placement of the sampling stations can be considered the most critical factor in the design of any water quality monitoring network. The selection of these optimal sampling points has been addressed by several authors, but mainly from a statistical viewpoint (a geostatistical approach combined with simulated annealing [3,4], fuzzy logic based on a geographic information system [5], multivariate statistical techniques [6], cellular automata-Markov chain models [7], graphical optimization by interpolation via correlation coefficients and standard deviations [8,9], Kriging variance combined with simulated annealing [10], a profile likelihood approach [11], etc.).
The aim of our research is to present a novel and effective approach to the problem of the optimal sampling points allocation within a simulation-based optimization framework, in the spirit of previous works of the authors for the case of a river water quality monitoring system [12][13][14], although this previous one-dimensional issue was a much simpler problem, both from the simulation viewpoint and from the optimization one. Here, we formulate the problem as a two-level optimization problem, where the upper level problem is the optimal fixing of the sampling points locations -given by their ability to capture the correct information on intensity and location of possible pollution releases -and the lower level problems are related to the optimal determination of these pollution sources.
In a specific way, the upper level optimization problem concerns the finding of the optimal sampling locations which best determine a large number of random point source pollution episodes. This problem can be formulated as an optimization problem where the objective function -measuring the global accuracy of the set of sampling stations -is given by the sum of the optimal approximation errors at the set of sampling points for all the different source pollution cases considered in the above formulation. In our study, the minimization process is executed via a controlled random search procedure for global optimization [15] in order to try to avoid the possibility of being trapped in local minima. The (subsidiary) lower level problems are related to the optimal identification (location and intensity) of the numerous random pollution sources, that is a critical step in managing the quality of estuarine waters. This problem has been much more extensively studied (see, for instance, a general survey in the recent review [16] for surface and groundwaters, and the numerous references therein). However, in our bidimensional case, this inverse problem can be mathematically ill-posed (in this sense, several authors have proved the well-posedness of the problem if we have data for the whole boundary of the domain [17,18], or if the known data from three sampling points in the case of the unbounded whole bidimensional space [19,20]; but the case of a finite number of sampling points in a bounded domain remains, as far as we know, as an open problem). This identification problem can be also reformulated as an optimization problem where the cost function depends on the differences between the observed and the predicted values for the different pollutant concentrations at the sampling points. There exists a large variety of proposed methods for solving the pollution source identification problem, but they can be categorized into three main groups, according to their approach: the probability-based approach (including Bayesian inference [21], backward probability method [22], the minimum relative entropy [23], and many others), the classification approach [24,25], and, finally, the optimization-based approach where the differences between simulated and observed pollutant concentrations at several points -obtained by solving a numerical model for pollutant concentrations -are minimized by means of a large range of optimization algorithms of derivative, derivative-free, or hybrid type [26][27][28]. In our case, we have chosen this linked simulation-optimization option, employing the gradient-free Nelder-Mead algorithm [29] for the minimization process and a convection-diffusionreaction equation for the simulation step.
This article is divided into five sections. After this introduction follows a second section devoted to the rigorous formulation of the problem, whose computational implementation is detailed in Section 3. Final sections are devoted to present the numerical results for a case study posed in Ría of Vigo (NW Spain), showing several discussions and conclusions.

Mathematical Setting of the Problem
We consider a domain Ω ⊂ ℝ 2 occupied by shallow waters, for instance an estuary or a ría (river end flowing into the sea), and we are interested in determining the optimal locations of a (usually small) number N of sampling points in a water quality monitoring network, that is, we want to find the best locations p i ∈ Ω, i = 1, … , N, for the sampling points, with the only constraint that each We understand as the best locations a vector So, if we denote by c (m,b) (x, t) the concentration at point x ∈ Ω and at time t ∈ (0, T) of an undesired pollutant (say, for instance, coliform bacteria Escherichia coli) coming from a discharge of intensity m at point b, then its evolution along Ω × (0, T) can be obtained as the solution of the following initial/boundary value problem (see, for instance, [30]): where Γ is the boundary of Ω (assumed to be smooth enough), m(t) is the mass flow rate of E. coli discharged in b, and (x − b) denotes the Dirac measure at point b. Experimentally known parameters and correspond to horizontal viscosity and decay rate, respectively. Finally, h(x, t) denotes the water depth, and (x, t) is the depth-averaged horizontal velocity of water, which can be measured in situ or estimated as the solution of the classical shallow water equations. Then, for each set of sampling points p ∈ ∏ N i=1 U i and for each set of sampling data d(t) ∈ [C(0, T)] N , we define the cost function: which measures the difference between the pollution levels caused by a discharge of intensity m at point b, and the levels collected in the samples taken at points p i , for i = 1, … , N . The objective is to determine the vector p so that each pollution discharge (m,b) can be recovered by numerically solving the inverse problem: being d (t) the pollution level samples, observed at monitoring points p, caused by discharge (m,b).
It is worthwhile remarking here that, for any set of locations p, if for a discharge (m,b) we solve the problem Eq. (1) and consider the synthetic samples given by is a solution of the inverse problem Eq. (3). However, we need to assure a correct numerical resolution of this inverse problem, which strongly depends on the chosen vector p, since an unsuitable choice of the set of sampling locations might lead to an inaccurate solution, posibly very different to the theoretical solution (m,b) . Our methodology aims to determine the best set p of sampling points location, which can recognize all possible discharges.
So, in order to determine the good quality of a particular set of sampling points locations p = (p 1 , … , p N ) , we consider a (large) number M of random pollution sources (located at points b j ∈ Ω and with constant intensities m j ∈ [m min , m max ] , for j = 1, … , M ), which must be identified by the monitoring network from data given by this particular set of sampling points.
After solving problem Eq. (1) for the different pollution scenarios (m j , b j ), j = 1, … , M, we can obtain the synthetic data d (2) Thus, we can define the function: which determines the goodness of the particular set of sampling points locations p, that is, function J j (p) measures the accurateness given by the set of sampling points locations p = (p 1 , … , p N ) when identifying the j-th pollution source (m j , b j ).
Finally, in order to find the best locations for the sampling points, we need to solve the following global optimization problem ( P ): find p = (p 1 , … ,p N ) , with p i ∈ U i , ∀i ∈ {1, … , N}, minimizing the objective function J given by: that is, we look for the optimal set of sampling points locations that identifies in the most accurate way all the random pollution scenarios chosen at the beginning of the monitoring network design process. Last but not least, we must note that, due to the strong nonlinearities of the global problem, uniqueness for solution p is not expected. Nevertheless, this does not represent any difficulty for our approach, since any of these possible global minima is good enough and suitable for our purposes. So, it suffices to compute one of them.
For readers' convenience, a detailed flowchart corresponding to the global process for the optimal network design can be found in Fig. 1.

Numerical Implementation
In this section, we present the full details for the computational resolution of the problem by means of a suitable discretization process, addressing the numerical resolution of the boundary value problem Eq. (1), the minimization of the intermediate optimization problems ( P j p ), and the resolution of the global optimization problem ( P).
In particular, for solving problem Eq. (1), we consider the standard variational formulation of the problem, and apply finite element method techniques for its resolution on a triangular mesh Ω h of the domain. To do this, we use the open-source finite element software Freefem++ [31], through a full programming of the associated formulation. Moreover, in order to assure the robustness of our approach, we have compared our achieved results to those obtained by the 2D finite volume hydrodynamic model MIKE 21 [32], developed by the Danish Institute of Technology (DHI), showing a good agreement, both from a qualitative and a quantitative viewpoint.
When solving the intermediate minimization problems ( P j p ), for any particular p and for each j ∈ {1, … , M} , given the essentially geometric nature of the problem, the authors propose to use a direct search algorithm: the Nelder-Mead simplex method [29]. This gradient-free algorithm has been successfully used by the authors in other related environmental problems (see, for instance, their previous work [33], where a short description of the method can be also found), and presents good convergence properties in low dimensions (in our particular case, we are dealing with a three-dimensional design variable (m j , b j ) ∈ [m min , m max ] × Ω ⊂ ℝ 3 ). In addition, the classical Nelder-Mead algorithm can be effectively modified with an oriented restarting when stagnation at a nonoptimal point is detected. However, since the Nelder-Mead algorithm was originally designed for unconstrained minimization problems, in order to apply it to the constrained optimization problem ( P p j ), we need first to modify the corresponding cost function by adding a penalty term related to the fulfilling of the constraints (m j , b j ) ∈ [m min , m max ] × Ω , which can be made in a simple and straightforward way: in case that m j ∉ [m min , m max ] or b j ∉ Ω we add to the original cost funtion a high penalty value, what makes vector (m j , b j ) inadmissible to be a minimum value.
Finally, for solving the global minimization problem ( P ), we use a controlled random search procedure for global optimization [15] combined with a multi-start strategy in order to assure a better performance of the algorithm. Again, as in the intermediate optimization problems, constraints related to p i ∈ U i , i = 1, … , N need to be penalized in cost function Eq. (6) by adding a penalty term.

A Case Study
This section presents some numerical tests for a real-world scenario posed in the estuary Ría of Vigo (Galicia, NW Spain). This shallow water region, whose finite element mesh Ω h is depicted in Fig. 2 [34], with water depth h and horizontal velocity computed by means of our own Fortran code.
In our case, due to budget constraints, only N = 5 monitoring stations will be allocated. Then, in order to avoid an incorrect accumulation of monitors in too limited areas, we have decided to divide the estuary into five vertical stripes and place a monitoring station in each of them. In particular, the case shown here corresponds to the following five admissible areas: y) ∈ Ω ∶ 520.9 ≤ x ≤ 526.6} , and U 5 = {(x, y) ∈ Ω ∶ 526.6 ≤ x}.
To determine the goodness of each set of sampling locations p, we employ following M = These achieved optimal and sub-optimal locations for sampling points are shown in Figs. 3 and 4, respectively. Finally, for both solutions, the identified locations and intensities for the synthetic discharges can be seen in Table 1.
By a straightforward analysis of above results, we can see how both the optimal and the sub-optimal solution are able to identify in a very accurate way all the random synthetic discharges -employed to calibrate the goodness of the set of sampling points locations -giving exact intensities and locations in the optimal case, and almost exact results in the sub-optimal one.
We must also note that, contrary to the suboptimal case where the five monitors are completely separated, in the optimal distribution case the monitors p 2 and p 3 are practically stuck together, which could indicate that the number of sampling stations could be maybe reduced from five to four without a loss of quality.  B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 U U U U U U U U U U U U U U U U U U U U U U U U U U U U U U 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 U U U U U U U U U U U U U U U U U U U U U U U U U U U U 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3  Galicia Gal Gal G Ga Gal Gal G G l l l G G Ga Ga Ga G G Gal Gal Gal Ga G G G G Gal Gal G G Ga G G Ga Ga Ga Ga Gal Gal G Ga Ga Ga a Ga a Ga Gal G Ga a Ga a a Gal a Gal a Ga a a G G G G G Ga a a a Gal G G G G Ga a a a a a a a  G Ga a a a a a ic c a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a

Conclusions
This paper proposes a new technique to automate the design of the sampling points for an estuarine water quality monitoring network by means of a linked simulation-optimization algorithm. After presenting a detailed and rigorous formulation of the problem, including its whole computational details, we have studied a real-world case posed in Rìa of Vigo (NW Spain), where the achieved optimal solutions show a very good ability to capture both the locations and the intensities of a large amount of possible discharges in the estuary. Moreover, although we have formulated our problem for the particular case of the concentration of E. coli in an estuary, our methodology can be immediately extended with the minimal changes to the analysis of any other water quality indicator -or indicators -in any type of 1D, 2D, or 3D domains.
The novel methodology introduced here represents not only an advance towards the scientifical rationalization of the design of water quality monitoring systems, but it also shows its wide possibilities in other different fields of application (atmospheric contamination, groundwater pollution...), and for other different interests from the stakeholders and decision-makers (pollution detection in minimal time, minimization of the number of monitoring stations...).