2.1.1. The Mazzocchio Reclamation Area and the Drainage Hydraulic Network

The Mazzocchio basin, with a surface about equal to 103 km2, occupies the largest depression of the Pontina Plain, which lies between the Appia and the Ufente River. The drainage of the basin is provided by the pumping plant of Mazzocchio (see Figure 1a) located at the downstream end of the Selcella river. The Selcella river is 18.9 km long with an average slope (h/L) of 0.00035. The terrain of the Mazzocchio basin is shown in Figure 1b; the ground height above mean sea level ranges from −2 m in the central stretch of the Selcella river to +8 m in the upstream region. From Figure 1b, we can observe the presence of a large area with about the same or lower height than the mean sea level. Furthermore, a dense channel network hierarchically structured, can be observed in Figure 1a collecting water to the Selcella River. All the water that flows to the Mazzocchio plant is lifted and through a short channel, conveyed into the Ufente river. The Mazzocchio pumping system (Figure 2) at the downstream end of the river Selcella consists of six pump groups of 6 m3/s each with a total capacity of about 36 m3/s. The pumping system was built on 20 July 1934 and inaugurated on 19 December of the same year; it is equipped with six engines and six water-immersion pumps with a modern conception for the era, each with a capacity of 6000 L per second and powered by an electric motor of about 600 hp. During the Second World War, German troops sabotaged the plant, taking away the engines. Flooding of the surrounding land meant hindering the entry by the allies. In 1948, near the Brenner railway station, located in Northern Italy, the engines were found and immediately reinstalled.

**Figure 2.** Hydraulic pumps in Mazzocchio's pumping station.

The starting and stopping of each pump is ruled by floating switches, which depend on the prefixed free surface levels in the collection pond immediately upstream from the pumping station. The depth of the pond close to the floating switches is −4 m.a.s.l. The pump starting and stopping levels are reported in Table 1.

**Table 1.** Starting and stopping levels for each pump (data provided by Agro Pontino Reclamation Consortium (APRC)).


Downstream of the pumping system of Mazzocchio, the waters flow into the Ufente River, which flows into the hydraulic node of Ponte Maggiore. Ponte Maggiore also receives water from the Amaseno River and the Pio Line Channel flow. From the Ponte Maggiore node, through the Portatore Channel, all the water flows into the sea (see Figure 1a). The main characteristics of the hydraulic network downstream from the pumping station are reported in Table 2.

**Table 2.** Hydraulic characteristics of channels, rivers, and basins (data provided by Agro Pontino Reclamation Consortium (APRC)).


The geometric and morphological characteristic of the entire hydraulic network (river slope, size and shape of cross-sections, etc.) as well as maintenance and operation of the pumping systems described above were provided by the Agro Pontino Reclamation Consortium (APRC).

#### 2.1.2. Precipitation and Pumping Rate Data

The closest station to the study area is the Borgo San Michele station (Lat. 41◦25 12 N Lon. 12◦58 12 E), from which 56 years (1950–2005) of hourly rainfall data has been collected. From this record, intensity and duration curves, with return periods equal to 10 and 100 years are shown in Figure 3. Hourly data of pump operation for the period 1950–2016, recorded in the archive of the Mazzocchio pumping station and provided by APRC, are also collected. From the time series of hourly rainfall amount, a number of heavy rainfall events that occurred in the past are identified, and the corresponding temporal trends of pumping discharge during such events are reconstructed. The latter data are used to calibrate the hydraulic model. The data related to the ground surface level, the bed profiles of rivers, and the geometry of river cross sections were also provided by the APRC.

**Figure 3.** Intensity–duration curve from hourly rainfall records at the rain gauge located in Borgo San Michele station (Lat. 41◦25 12 N Lon. 12◦58 12 E).

#### *2.2. Methods*

The simulation–optimization model combines a hydraulic model and a multiobjective optimization model. The hydraulic model calculates the free surface level and the flowrate along rivers and channels upstream and downstream of the Mazzocchio pumping station, the flooding areas over the basins, and the pumped flowrate at the pumping station. The multiobjective optimization genetic algorithm, using the outputs from the hydraulic model, calculates the objective functions, verifies the constraints violations, and by an iterative procedure based on tournament selection of nondominant solutions and generation of new populations by crossover and mutation, identifies the set of Pareto optimal solutions. A sketch of the flowchart of the combined simulation–optimization model is shown in Figure 4.

**Figure 4.** Optimization model's flowchart.

#### 2.2.1. Hydraulic Model

For given input variables—rainfall amount over the basins and sea levels—as well as for given pumping switching levels in the collection pond upstream of the Mazzocchio pumping station, the hydraulic model calculates the outputs related to free surface levels and flowrates along the entire hydraulic network, flooding surface, water volumes, and depths over the basins. To simulate the temporal and spatial distribution of these quantities, 2D and 3D hydraulic models are generally used, as proposed by Cioffi and Gallerano [26] and Orton et al. [27]. Such models require a computation time which conflicts with the very high number of simulations necessary for the multiobjective optimization algorithm to calculate the Pareto front. Therefore, in this paper, a simplified hydraulic model has been constructed. This model is able to calculate the above listed hydraulic quantities after calibration. The scheme of such models, for the hydraulic networks upstream and downstream from the Mazzocchio pumping station, is shown in Figure 5.

Three different hydraulic elements may be recognized in the figure: the main river hydraulic networks upstream and downstream from the pumping station, the storage areas representing the basins, and the ideal channels connecting the storage areas to points of the rivers belonging to the hydraulic networks. As seen in the figure, the basin of Mazzocchio, upstream from the pumping station, is represented by seven storage areas connected to the river Selcella by ideal channels. Figure 6 shows how the storage areas refer to parts of the basin with homogeneous morphology and altimetry, in order to provide a sufficiently accurate representation of the spatial distribution of ground surface elevation.

**Figure 5.** Sketch of hydraulic model.

**Figure 6.** Storage areas partition of the Mazzocchio basin.

The ideal channels mimic the secondary drainage network collecting the rainfall water to the Selcella river. For the hydraulic network downstream from the pumping station, storage areas and related ideal channels are located upstream of the Amaseno, Ufente, and Linea Pio rivers. The combination of storage areas and ideal channels allows for a representation of both the rainfall–runoff processes on the basins and the flooding over the riverbanks due to levee overflowing.

The 1D de Saint-Venant equations are used to simulate the flow along the river which belongs to the hydraulic networks:

$$\frac{\partial \mathcal{Q}}{\partial \mathbf{x}} + \frac{\partial \mathcal{S}}{\partial t} = q\_l \tag{1}$$

$$\frac{\partial \mathcal{Q}}{\partial t} + \frac{\partial}{\partial x} \left( \beta \frac{\mathcal{Q}^2}{S} \right) + \mathcal{g} \mathcal{S} \left( \frac{\partial H}{\partial x} + I \right) = \gamma\_I \tag{2}$$

where *Q* is the flow rate, *S* the cross-section area, *H* the free surface elevation from a reference plane, and *ql* represents the lateral inflow, *γ<sup>l</sup>* represents the sum of the forces applied, and *J* is a dimensionless number, representing the average rate of energy dissipation:

$$J = \frac{Q^2}{K\_m^2 S^2 R^{4/3}}\tag{3}$$

where *Km* is the Strickler's roughness coefficient and *R* the hydraulic radius.

A number of river junctions (nodes) are present in the hydraulic network downstream from the pumping station (Figure 5). The equations of the junctions are derived from the equality of the elevations and the conservation of the discharges at the junction following the approach of 1D–2D coupling proposed by Goutal et al. [28]. The temporal trend of water level over the storage areas is calculated by the continuity equation, which is a function of the flows entering or going out from the ideal channels and of the rainfall amount directly falling over the area. It is assumed that the water free surface of the storage areas remains horizontal while moving vertically. The link between the storage area and the river is seen in Figure 7.

**Figure 7.** Sketch of the link between a river and a storage cell by an ideal channel. Hc is the bottom elevation of the ideal channel, h1 and h2 are free surface elevation in the river and in the storage cell respectively.

The flow rate through the ideal channels is calculated by the Manning–Stricker uniform free surface flow formula and depends on the difference of the free surface levels in the storage areas and in the river according to the following equations:

$$Q = K l h^{5/3} \sqrt{\frac{\Delta H}{L}} \tag{4}$$

where, as shown in Figure 7, Δ*H* is the difference between the free surface elevation *h*<sup>2</sup> in the storage cell and *h*<sup>1</sup> in the river, *L* is the length of the channel, *l* is the channel width, and *h* is the water depth in the channel, which is equal to the average between *h*<sup>2</sup> and *h*<sup>1</sup> (see Figure 7). Both *h*<sup>2</sup> and *h*<sup>1</sup> depend on the bottom elevation *Hc* of the ideal channel.

#### 2.2.2. Multiobjective Optimization Problem Formalization

In formalizing the multiobjective problem, we assume the levels switching on and off the pumps in the Mazzocchio station (*Hs* ) as decisional variables. Temporal and spatial distribution of the rainfall amount (*I*) and the sea levels (*Sl*) are the input variables. The state variables are the free surface levels - *H<sup>d</sup>* , the flow rates (*Qd*) along the hydraulic networks, and the surface of the flooding areas (*A<sup>m</sup> i* ). State and decisional variables range within imposed constraints. Specifically, we impose that the free surface level in the hydraulic network downstream from the pumping station cannot be more than the height of the river levees - *H<sup>d</sup> max* . The multiobjective optimization problem is formalized defining two optimality criteria: (a) to minimize the maximum surface of flooding areas (*Am*) over the Mazzocchio basin during a flood event caused by heavy precipitation; (b) to minimize the pumping power necessary for flooding control over the Mazzocchio basin. The pumping power here is defined as the energy consumed to lift water from the beginning of the event to the moment in which the maximum flooding surface over the Mazzocchio basin is reached. The second criterion is aimed to identify pumping schedules that limit, as much as possible, the power used to reduce the flooding surface. Such a criterion also responds to the need to limit the number of pumps working simultaneously in order to allow a more efficient and robust maintenance and operation of the plant.

The multiobjective algorithm is aimed to identify the levels of switching on and off the pumps in the Mazzocchio station which satisfy the optimality criteria described above. The mathematical formulation of the multiobjective optimization problem and the Pareto set calculation algorithms are as follows. The objective functions may be expressed as:

$$Of1 = \left( max\left(A\_1^m, \ A\_2^m, \dots, A\_n^m\right) \right) \tag{5}$$

$$
\Delta O \mathbf{f} \mathbf{2} = \left( \sum\_{i=1}^{T\_{\text{Amax}}} \gamma Q\_i^p \left( H\_i^m - H\_i^p \right) \Delta \mathbf{t} \right) \quad (i = 1, n) \tag{6}
$$

where *A<sup>m</sup> <sup>i</sup>* , *<sup>Q</sup><sup>p</sup> <sup>i</sup>* , *<sup>H</sup><sup>m</sup> <sup>i</sup>* , and *<sup>H</sup><sup>p</sup> <sup>i</sup>* at each time step *i* (*i* = 1, *n*) are the surface of the flooding area over the Mazzocchio basin, the pumping flowrate at the Mazzocchio pumping station, the free surface levels (from the reference level) in the pond of the Mazzocchio pumping station, and the free surface levels in the basin downstream from the pumping station, respectively; *TAmax* is the temporal step at which the flooding surface in Mazzocchio reaches the maximum value; *n* is the number of time steps in which the period of hydraulic simulation *Ts* is divided; and Δ*t* the temporal step.

The objective functions *Of* 1 and *Of* 2 are constrained by:

$$Q\_i^p = n\_k^p \cdot \Delta Q \text{ if } \ H\_i^m \ge H\_k^s \ k = 1, n\_p \text{ } i = 1, n \tag{7}$$

$$H\_{i,j}^d \le H\_{\max,j}^d \; j = 1, n\_c \; i = 1, n \tag{8}$$

where *n<sup>p</sup> <sup>k</sup>* is the number of pumps working simultaneously, Δ*Q* is the flowrate of each single pump, *Hs <sup>k</sup>* is the level in the collecting pond upstream of the pumping station switching on the *k*th pump, and *H<sup>d</sup> <sup>i</sup>*,*<sup>j</sup>* and *<sup>H</sup><sup>d</sup> max*,*<sup>j</sup>* are the free surface level and the levee top level, respectively, in the *j*th point along the river stretches of the hydraulic network downstream from the Mazzocchio pumping station.

Equation (7) imposes a constraint on the maximum number of pumps that may work simultaneously and the pumping switching on/off levels. Equation (8) constraints mean that the free surface levels along the hydraulic network downstream from the pumping station must be lower or equal to the top level of the levees of the Ufente, Amaseno, and Linea Pio Rivers and in the Portatore Channel.

The variables *A<sup>m</sup> <sup>i</sup>* , *<sup>Q</sup><sup>p</sup> <sup>i</sup>* , *<sup>H</sup><sup>m</sup> <sup>i</sup>* , *<sup>H</sup><sup>p</sup> <sup>i</sup>* and *<sup>H</sup><sup>d</sup> <sup>i</sup>*,*<sup>j</sup>* are calculated by the hydraulic model as a function of the decision variable *H<sup>s</sup> <sup>k</sup>* and of the input variables *Ii* and *Sli*, that is, the temporal trend of the rainfall over the storage areas and the sea level; the latter varies in time with the tidal cycle.

*Water* **2018**, *10*, 820

All the hydrodynamic variables of the river flows, that is, free surface levels and flow rates along the two hydraulic networks, are subject to the constraints imposed by the continuity and momentum equations previously described (see Equation (1)).

For given inputs *Ii* and *Sli*, the resolution of the multiobjective problem requires the minimization of the two objectives *Of* 1 and *Of* 2 in Equations (5) and (6), within the constraints imposed by Equations (7) and (8), as a function of the decision vector H<sup>s</sup> <sup>k</sup> with *k* = 1, *np*.

Given two decision vectors *H<sup>s</sup> <sup>k</sup>*,1 and *<sup>H</sup><sup>s</sup> <sup>k</sup>*,2 with *k* = 1, *np*, solutions of the multiobjective problem, the identification of the Pareto front is obtained applying Pareto dominance criteria.

In accordance with such criteria, an objective vector <sup>→</sup> *y*<sup>1</sup> ≡ / *O f* 12,*O f* 22 0 is said to dominate another objective vector <sup>→</sup> *y*<sup>2</sup> ≡ / *O f* 12,*O f* 22 0 (i.e., <sup>→</sup> *<sup>y</sup>*<sup>1</sup> < <sup>→</sup> *y*2), if no component of <sup>→</sup> *y*<sup>1</sup> is greater than the corresponding components of <sup>→</sup> *y*<sup>2</sup> and at least one component of <sup>→</sup> *y*<sup>2</sup> is greater; consequently, the solution *Hs <sup>k</sup>*,1 dominates *<sup>H</sup><sup>s</sup> <sup>k</sup>*,2. The nondominant solutions are optimal solutions of the problem whose ensemble identifies the set of Pareto optimal solutions.

A genetic algorithm is applied to solve the multiobjective optimization problem formalized above. Genetic algorithms are based on Darwin's theory of natural selection, which involves the language of microbiology, and in developing new potential solutions, mimics genetic operations. A population represents a group of solution points, and a generation represents the algorithm iteration, while a chromosome is equivalent to a component of the design vector. In accordance with these definitions, a genetic algorithm deals with a population of points, and hence multiple Pareto optimal solutions can be obtained from a population in a single run. Random numbers and information from previous iterations are combined to evaluate and improve a population of points, and then to select nondominant solutions. In this paper, the nondominant-sorting genetic algorithm II described by Deb et al. [29], NSGA2, is used, which has been applied successfully to many optimization problems. This algorithm uses tournament selection [30], simulated binary crossover (SBX) [31], a mutation operator, and crowding distance for diversity preservation.

#### **3. Results**

#### *3.1. Hydraulic Model Calibration*

As mentioned in Section 2.2.1, in order to limit the computation time needed to run the combined optimization–hydraulic models, the simplified hydraulic model previously described has been constructed. Such a model has a number of parameters that have to be identified: in particular, the length *L*, average heights *h*<sup>1</sup> and *h*2, and width *l* of the ideal channels linking the rivers to the storage areas, in which the different basins have been divided, as well as the runoff coefficient φ of the storage areas and the Strickler's roughness coefficient *Km* in the de Saint-Venant equations.

To identify such parameters, a calibration procedure was applied. Such a procedure consists of carrying out a single-objective optimization, through which the parameter values that minimize the sum of squares of the difference between the simulated and observed hydraulic data are identified.

Since the hydraulic network upstream of the Mazzocchio pumping station is hydraulically disconnected by the pumping from the downstream network, calibration of the upstream hydraulic network can be performed separately. For this network, two datasets were selected: the first is the precipitation trend recorded at the Borgo San Michele station, and the second is the pumping discharge trend at the Mazzocchio pumping station. Collected data from five different extreme events that occurred in the past were used to calibrate the hydraulic model. Then, the set of hydraulic parameters in the Mazzocchio area were identified, which minimizes the sum of squares of the difference between the simulated and the recorded pumping discharges during that period. The values of the identified parameters are shown in Table 3. For an extreme event, during the time period 16–20 March 2011, recorded pumping discharge and simulated pumping discharge was compared to validate the model. Figure 8a,b shows the precipitation trend used as input as well as the comparison between the

simulated and recorded pumping discharges. Model accuracy was tested using the root mean square error (RMSE) between observed and simulated data, giving ~70% of accuracy.

**Figure 8.** (**a**) Comparison between simulated and observed pumping discharges. (**b**) Rainfall intensity trend during the extreme event in March 2011.

**Table 3.** Hydraulic parameters after the calibration of the model upstream from the Mazzocchio pumping system (data provided by Agro Pontino Reclamation Consortium (APRC)).


For the hydraulic network downstream from the pumping station, we used a different heavy rainfall event (from 12–17 December 2008), since this was the only event where there were simultaneous measurements of both free surface level at Ponte Maggiore's hydrometric station and rainfall amount trend in the Pontinia rain gauge. In this case, the optimization was performed by running the simulation model of both the hydraulic networks forced by the same rainfall amount trend. In this last case, the objective function was formalized using the measured free surface levels at the hydrometric gauge of Ponte Maggiore. The values of the parameters that minimized the objective functions and referred to the ideal channels connecting the storage areas of Linea Pio, Ufente, and Amaseno Rivers are reported in Table 4.

**Table 4.** Hydraulic parameters after the calibration of the model downstream from the Mazzocchio pumping system.


#### *3.2. Optimal Pareto Set for the Event in March 2011*

An optimal Pareto curve was calculated by using the rainfall temporal distribution of the event on 16–20 March 2011 as an input (see Figure 8). The main aim is to verify whether and how much the pumping schedule applied during such an event was far from the Pareto front. Furthermore, we also test the process of convergence to the Pareto front of the multiobjective optimization algorithm, finally selecting 100 populations and 20 generations. In the analysis, the discrete levels of free surface *Hs*, regulating the switching on or off of the pumps, was varied within the prefixed range reported in Table 5.


**Table 5.** Variation range for activation levels.

In Figure 9, in the space of the objective functions, the points belonging to the optimal Pareto set as well as the points obtained from the different generations show that 20 generations were sufficient to allow for the convergence of the algorithm at the Pareto front. From Figure 9, it is evident how the two objectives are conflicting, since one objective cannot be improved without sacrificing the other one; i.e., if less pumping power is consumed during the event, a greater surface of the Mazzocchio basin is flooded.

**Figure 9.** Pareto front obtained for the rainfall trend of the event that occurred in March 2011.

In the figure, the green-filled square refers to the values of the two objective functions obtained assuming the real operative pumping schedule applied during the rainfall event by the pumping station operators. From Figure 9, it is evident that such a point is far from the set of optimal Pareto solutions. This suggests that a more efficient pump schedule could have been identified to manage the extreme event of March 2011 between the solutions belonging to the Pareto front. As suggested by Kurek and Ostfeld [32], we used a procedure based on an utopian solution [33] to identify the optimal solution, circled in blue on Figure 9. Therefore, other comparison methods could be employed, such as game theory techniques [34]. The utopian solution has been identified as minimizing the objective functions:

$$m = \left\{ \min\{Of\_1\}, \min\{Of\_2\} \right\} \tag{9}$$

Then, the span of the front Δ*O fi* for each of the considered objectives is computed:

$$
\Delta Of\_i = \max\{Of\_i\} - \min\{Of\_i\} \quad i = 1, 2 \tag{10}
$$

Subsequently, the fronts are normalized using their span in the objective space; finally, a solution having the minimum distance to the utopian solution is selected for comparison:

$$\min\_{p \in \Omega} \left\{ \left[ \left( \frac{Of\_1 - m\_1}{\Delta O f\_1} \right)^2 + \left( \frac{Of\_2 - m\_2}{\Delta O f\_2} \right)^2 \right]^{\frac{1}{2}} \right\} \tag{11}$$

where Ω is the Pareto set resulting from solving the problem (5) or (6) and *p* is the selected "balanced" solution.

#### *3.3. Sets of Pareto Optimal Solutions for Heavy Rainfall Events and Different Sea Level Rises*

Running the combined simulation–optimization model, sets of Pareto optimal solutions were obtained for different design hyetographs associated to prefixed return periods of daily rainfall amount and different mean sea level. In this study, we used daily rainfall amount because previous studies [11,35] have shown that 24-h-long heavy rainfall (with a daily rainfall amount greater than 100 mm) has in the past induced serious flash flooding in the examined site. It should also be underlined that downscaling models aimed to project future changes in the precipitation regime [10] generally refer to the daily rainfall amount, and therefore it is reasonable to use the return period associated to such a quantity. In order to construct the design hyetograph, that is, the artificial rainfall temporal distribution, having a given 24-h rainfall amount return period, the dimensionless approach suggested by Kimura et al. [35], also called the modified ranking method, was used.

Such an approach can be summarized in the following steps:


Using the rainfall data recorded in the Borgo San Michele rain gauge, the dimensionless hyetograph shown in Figure 10 was calculated.

Then, the hyetographs referring to the daily rainfall amount with return periods (TR ) equal to 10 and 100 years were assumed as inputs of the simulation–optimization model, and different sets of Pareto optimal solutions were calculated. While calculating the sets, we impose a further constraint on the maximum number of simultaneously working pumps. In Figure 11, two sets of Pareto optimal solutions are compared for a return period (TR) equal to 100 years.

The first sets refer to the configuration of the system with unlimited flow carrying capacity of the hydraulic network downstream from the pumping station; that is, the Equation (6) constraint was removed. The second sets refer to the case in which the constraint of Equation (6) remains. As expected, within the admissible region defined by the constraints of Equation (6), the Pareto optimal solutions of the two sets are coinciding. The comparison in Figure 11 shows that the constraints of Equation (6) (the free surface levels must be lower than the height of levees in the downstream channel network) significantly reduces the discharge from the Mazzocchio basin that can be pumped into the river Ufente. For the reconstructed hyetograph with a return period (TR) equal to 100 years, no more than four pumps can simultaneously work without violating the constraints. Therefore, due to the limited flow carrying capacity of the hydraulic network downstream from the pumping system, not all the potential of the pumping system to mitigate the flooding in the Mazzocchio basin can be achieved. It would be important to investigate in detail why this occurs.

**Figure 10.** Design hyetograph.

**Figure 11.** Pareto sets for daily rainfall amount with return period of 100 years. (**a**) Pareto set with all the pumping configurations. (**b**) Optimal Pareto set with a maximum of six pumps simultaneously working; and with (**c**) maximum of five pumps; (**d**) maximum of four pumps; (**e**) maximum of three pumps; (**f**) maximum of two pumps.

Figure 12 shows the temporal trend of the flow rates in significant cross sections of the Ufente, Amaseno, and Portatore rivers, close to the Ponte Maggiore node and the free surface level at the junction. These trends refer to the configuration of pump switching on/off corresponding to the Pareto optimal solution in the middle of the range, with a maximum number of working pumps equal to five. From Figure 12, it can be observed that a large increase in discharge, up to 250 m3/s in the Amaseno river, causes a significant increase in the free surface level at the Ponte Maggiore node. Such an increase

in the free surface level produces a backwater with a consequent inversion of the flow direction in the Ufente river (as seen by the negative values of the flow rate in the Figure 12), with a consequent blocking effect for the pumped flow rate from Mazzocchio. The consequences of this phenomenon can cause overflows and levee collapse along the Ufente River, which has already occurred in the recent past. From the results of hydraulic simulation, it is evident that the area included in a 3-km radius from Ponte Maggiore's junction, as seen in Figure 13, is at high risk of embankment overtop.

**Figure 12.** River discharges and free surface level trends at Ponte Maggiore's junction: rivers (**a**) Ufente, (**b**) Amaseno, and (**c**) Portatore, and (**d**) water level at node.

**Figure 13.** High risk of embankment overtop near Ponte Maggiore's junction using a 100-year return period hyetograph as the boundary condition for the hydraulic simulation.

In Figure 14, the sets of Pareto optimal solutions obtained for a hyetograph with the daily rainfall amount corresponding to a return period of T<sup>R</sup> = 10 years and referring to three different maximum sea levels (0, +0.5, +1.0 m) at the outlet of the Portatore river are shown. These levels are typical of the shore examined and are due to the combined effects of the astronomical tide, storm, and barometric surge.

**Figure 14.** Pareto sets for a likely return period of 10 years. (**a**) Pareto set with all the pumping configurations. (**b**) Optimal Pareto set with a maximum of six pumps simultaneously working; and with (**c**) maximum of five pumps; (**d**) maximum of four pumps; (**e**) maximum of three pumps; and (**f**) maximum of two pumps.

As expected, the mean sea level rise modifies the boundary condition at the outlet of the network, causing an increase in free surface level in a larger part of the hydraulic network, making the violation of level constraints more likely. As a consequence of such a violation of level constraints, there is a limitation of the number of the pumps that can simultaneously work. As the figure shows, for a sea level rise equal to 0.5 m, the maximum number of pumps that may simultaneously work is equal to four, while in the case of a sea level rise equal to 1.0 m, no more than three pumps may be used simultaneously.

#### **4. Discussion and Conclusions**

This study assessed in the Mazzocchio area whether pumping systems and related channel hydraulic networks are reliable in managing the risks of flooding associated to extreme rainfall events and sea level rise. To conduct this research, we developed and proposed an approach that takes into account the different pumping schedules which are possible at parity of hydrological inputs. To compare the different states of the hydraulic system as a function of the different hydrological inputs, we referred to sets of Pareto optimal solutions calculated by a simulation–optimization model, which combined a multiobjective evolutionary algorithm (the nondominating sorting genetic algorithm, NSGA2) and a hydraulic model in order to satisfy two optimality criteria: minimize the consumption of pumping power and minimize the flooding surface over the Mazzocchio basin. The application to the study case shows how the increase of extreme rainfall amount, as well as sea level rise, affects the reliability and efficiency of the pumping system and hydraulic channel networks to mitigate the flooding in the study site. The use of the model allows identification of pumping system management solutions that are more efficient than those currently used, which are based on empirical assumptions. As seen in Figure 9, the current operating configuration of Mazzocchio's pumping plant should be

more efficient in terms of energy consumption, with significant economic benefit. The Pareto set includes the optimal solutions, then through decision-making techniques, the preferred solutions can be identified. Pareto sets provide decision makers with a fundamental tool for defining and choosing the actions to be taken to reduce the hydraulic risk. For example, improving safety also means increasing costs. At present, the model can be part of a flood management system, included in a model system in which a submodel makes rainfall and sea level forecasts, starting from large-scale meteorological configurations. In this context, the developed model can be part of an optimal control system to manage both the pumping plant and the hydraulic network. The Pareto sets obtained from the optimization model (Figures 11 and 14) show that the hydraulic network will not be able to manage the increase in river flows due to increasingly intense rainfall and rising sea level, phenomena generated by climate changes. Therefore, the study suggests that in order to exploit the entire pumping capacity of the plant to mitigate the flooding in the Mazzocchio area, the flow carrying capacity of the downstream hydraulic channel network should be increased by designing new interventions; for instance, reshaping of the riverbeds and levees and designing structures for flow diversion or compensating–balancing reservoirs; all of those measures should be designed using a multiobjective optimization approach in order to minimize construction costs, but guaranteeing adequate protection from future extreme events. The identification of such technical solutions is beyond the goal of this paper, but the simulation–optimization model proposed here could still be used in order to explore such flood risk intervention strategies. This study could be further developed by taking into account the uncertainty in future projections of heavy rainfall and sea level rise, as suggested by Woodward et al. [25]. It is possible to take into account future uncertainties regarding sea level and rainfall amount and duration with the help of flexible sets of defense strategies linked to predefined situations. To assess such uncertainties, we are developing rainfall downscaling models which use General Circulation Models (GCM) ensemble simulations as inputs within a more general project aimed to assess and to identify possible strategies to preserve the coastal areas from the negative effect of climate change. Sea level rise will also be calculated using hydraulic models that will simulate storm surges caused by waves and wind, thus providing more accurate values.

**Author Contributions:** Conceptualization, F.C., A.D.B.T. and F.R.C.; Methodology, F.C., A.D.B.T. and F.R.C.; Software, F.C., A.D.B.T. and F.R.C.; Validation, F.C., A.D.B.T. and F.R.C.; Formal Analysis, F.C., A.D.B.T. and F.R.C.; Investigation, F.C., A.D.B.T. and F.R.C.; Resources, F.C., A.D.B.T. and F.R.C.; Data Curation, F.C., A.D.B.T. and F.R.C.; Writing-Original Draft Preparation, F.C., A.D.B.T. and F.R.C.; Writing-Review & Editing, F.C., A.D.B.T. and F.R.C.; Visualization, F.C., A.D.B.T. and F.R.C.; Supervision, F.C.; Project Administration, F.C.; Funding Acquisition, F.C.

**Funding:** This research was funded by Consorzio di Bonifica dell'Agro Pontino (APRC).

**Conflicts of Interest:** The authors declare no conflict of interest

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Water* Editorial Office E-mail: water@mdpi.com www.mdpi.com/journal/water

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-03943-387-2