Flood Modeling and Groundwater Flooding in Urbanized Reclamation Areas: The Case of Rome (Italy)

: Coastal regions below the sea level, subject to reclamation, are becoming more and more exposed to ﬂooding following increasing urbanization and hydrological changes. In these areas, groundwater and water table dynamics during intense rainfall events can be an important component of ﬂooding and inundation, leading to groundwater ﬂooding. Thus, the commonly employed hydrological models based on only the surface component of ﬂow may result in a poor estimation of the extension and persistence of inundation events. We introduce here a simple and parsimonious approach for handling the groundwater contribution to ﬂooding in such areas, which can be easily implemented and introduced into surface hydraulic models for ﬂood management and the delineation of inundation maps. The approach involves few relevant parameters, requiring a minimum of information regarding the hydrogeological setup. The method is exempliﬁed through the ﬂood analysis of a wide reclamation area located in the southern part of Rome, Italy. The introduction of the groundwater component could explain the large water volumes pumped by the stations, which are much larger than excess rainfall. The application conﬁrmed the validity of the proposed approach, emphasizing the important role played by groundwater to ﬂooding in areas similar to the one considered here.


Introduction
Low-elevation coastal areas are home to more than 600 million people, that is, roughly one in 10 people in the world live in such territories [1]. It turns out that two thirds of world's largest cities, i.e., those with more than five million people, are located (at least partially) in these low areas. According to the United Nations [2], nine of the ten largest urban agglomerations in the world are located in deltas or floodplain areas.
Floodplains usually offer favorable conditions for urban and economic development, and they have therefore attracted human settlements over the years. Small urban settlements evolve into dense cities, and are then vulnerable to high-consequence and low-probability events [3,4]. In such cases, structural and/or non-structural measures are required to control floods and reduce flood losses. Flooding is estimated to affect more people globally than any other natural hazard, with an average of 250 million people affected annually [5]. Overland flooding, e.g., caused by insufficient river capacity, has long been recognized as a risk, and is commonly included in land-use policies and insurance risk considerations [6,7]. Hazard maps and policies, generally based on flood recurrence intervals for design river discharges, are typically based on overland flooding, which is mainly studied using the Water 2020, 12 tools available in surface hydrology, leading to either pluvial (i.e., direct runoff from rainfall) and/or river flooding (due to overflow from rivers and channels).
Of particular interest are the coastal regions below sea level that were originally wetlands and have been reclaimed by lowering the groundwater level by pumping. Such coastal areas are widespread and have been the subject of intensive urbanization over the years, and are becoming more and more exposed to flooding. In those areas, groundwater and the water table dynamics during intense rainfall events can be an important component of flooding and inundation. Such behavior is called groundwater flooding, which is defined as the "emergence of groundwater at the ground surface away from perennial river channels or the rising of groundwater into man-made ground, under conditions where the 'normal' ranges of groundwater level and groundwater flow are exceeded" [8][9][10]. Flood damage can occur before the rise of the water table to the ground surface, e.g., in the case of inundation of subsurface structures (i.e., basements) by groundwater [7], or the rise of the water table due to riverbank filtration [11]. When compared to overland flooding, groundwater flooding is inadequately recognized and poorly understood worldwide [12], with the possible exception of the EU Flood Directive in 2007 [13]. In comparison to fluvial and pluvial flooding, flooding from groundwater has only recently been given adequate consideration following groundwater floods in Northern and Central Europe in the past two decades [14][15][16][17]. The available studies on groundwater flooding are rather scattered, with many dealing with different groundwater processes and case studies (e.g., References [18][19][20][21][22]).
In this work, we dealt with groundwater flooding determined by the sudden rise of the water table due to intense direct precipitation, as is characteristic in low areas, and particularly in coastal zones subject to reclamation. The latter are areas below sea level and are typically characterized by a network of natural and artificial channels, usually connected to the groundwater table, that drain water to a pumping station or a network of stations. Such systems keep the groundwater levels below the ground, and the drained water is further delivered to the sea by a system of channels. The pumps usually operate automatically as a function of the water levels in the pump chambers. The water level in the channels displays a sudden increase after rainfall due to runoff, which in turn recharges the aquifer by infiltration. Thus, when rainfall is intense and the water table is high, the dewatering stations receive both surface water and groundwater feeding the channels. As a consequence, the volume pumped is potentially larger than the total rain that has fallen on the surface catchment, being also impacted by the size of the hydrogeological (groundwater) basin.
The correct assessment of the flooding risk conditions in such areas, and the subsequent design of any mitigation measures, therefore requires the implementation of suitable hydrological and hydraulic models capable of representing in a robust, possibly parsimonious way (i) the joint surface-subsurface hydrological processes and (ii) the specific role played by the drainage network and the pumping stations, of which the (often limited) capacity determines the effective drainage of rainwater and the extension of the inundated areas. This is in contrast with the common practice of employing modeling approaches for flooding based on surface hydrological/hydraulic dynamics, typically neglecting the groundwater component [23].
In the present work, we introduce a simple approach for handling the groundwater contribution to flood in reclaimed areas, which can be easily implemented and introduced into surface hydraulic models for flood management and the delineation of inundation maps. The approach aims at integrating the groundwater contribution to flooding in a parsimonious manner, with few relevant parameters, requiring a minimum of information regarding the hydrogeological setup. The approach differs from other approaches involving direct groundwater simulations (e.g., References [7,15]) or from combined, fully numerical surface/groundwater modeling (e.g., Reference [18]); such approaches are more complex than the one employed here, requiring significant information, data, and computational resources. The proposed method was exemplified through the flood analysis and delineation of the inundation maps of a wide reclamation area located in the southern part of Rome, Italy, emphasizing the important role played by groundwater in flooding.

Site Description
The study area is located in southern Rome, Italy (Figure 1), along the Tyrrenian coast. The region includes a wide catchment that drains water to zones with elevations below the sea level; in the past, it was a wetland, with groundwater levels partially above the ground. The land was reclaimed in the previous century via the realization of a network of channels connected to a series of pumping stations. The area has seen a very intense urbanization in recent years, with an increasing population during the past century, and it is now incorporated into the municipality of Rome.
The study area is located in southern Rome, Italy (Figure 1), along the Tyrrenian coast. The region includes a wide catchment that drains water to zones with elevations below the sea level; in the past, it was a wetland, with groundwater levels partially above the ground. The land was reclaimed in the previous century via the realization of a network of channels connected to a series of pumping stations. The area has seen a very intense urbanization in recent years, with an increasing population during the past century, and it is now incorporated into the municipality of Rome.
The site has experienced several floods due to the insufficient transport capacity of the drainage channels and the pumping stations, and it is frequently inundated, with a relatively small return period of the events. The overall area was divided for simplicity into two major zones, as depicted in Figure 1c (the elevation map is provided in Figure 1d, where the areas below sea level are evident): (i) the areas at higher altitudes (hereinafter high areas, or HA), draining water by gravity directly into the Tyrrhenian Sea through a system of natural and artificial channels; and (ii) the part of the catchment draining water into areas below the sea level (low areas, or LA), for which relies on a system of pumps that deliver the excess water collected by the channel network to the sea. The study area (c), with its location within the Italian territory (a) and the Lazio region (b); the high and low areas and their drainage networks are represented by solid lines; the dots display the dewatering pumping stations. The shaded high areas were studied and modeled in a previous work [24]. The topography of the study area (d) clearly shows the zones with elevation below sea level. The red star in panel (d) is the center position of the sub-area represented in Figure 6.
The area under consideration covers an area of approximately 36 km 2 , of which about 13 km 2 flows towards the sea or is moved by gravity (HA) and the remaining 23 km 2 requires pumping, which is performed by four main pumping stations: Ostiense (catchment of about 9 km 2 ), Stagni (10 km 2 ), Bagnolo (2 km 2 ), and Tor San Michele (about 2 km 2 ). The water dynamics in the HA of higher elevation (shaded area in Figure 1), with an area 57 km 2 , were studied and modeled in a previous work [24].
The hydrological and hydraulic models employed in the present study required a preliminary campaign to acquire the data needed for an accurate description of the transport capacity of the existing drainage network and the hydraulic characterization of the area. Such information included the characteristics of the urbanization, the structure of the drainage networks, the various interconnections, the hydraulic cross-sections, and the state of conservation of the canalizations, as The study area (c), with its location within the Italian territory (a) and the Lazio region (b); the high and low areas and their drainage networks are represented by solid lines; the dots display the dewatering pumping stations. The shaded high areas were studied and modeled in a previous work [24]. The topography of the study area (d) clearly shows the zones with elevation below sea level. The red star in panel (d) is the center position of the sub-area represented in Figure 6.
The site has experienced several floods due to the insufficient transport capacity of the drainage channels and the pumping stations, and it is frequently inundated, with a relatively small return period of the events. The overall area was divided for simplicity into two major zones, as depicted in Figure 1c (the elevation map is provided in Figure 1d, where the areas below sea level are evident): (i) the areas at higher altitudes (hereinafter high areas, or HA), draining water by gravity directly into the Tyrrhenian Sea through a system of natural and artificial channels; and (ii) the part of the catchment draining water into areas below the sea level (low areas, or LA), for which relies on a system of pumps that deliver the excess water collected by the channel network to the sea.
The area under consideration covers an area of approximately 36 km 2 , of which about 13 km 2 flows towards the sea or is moved by gravity (HA) and the remaining 23 km 2 requires pumping, which is performed by four main pumping stations: Ostiense (catchment of about 9 km 2 ), Stagni (10 km 2 ), Bagnolo (2 km 2 ), and Tor San Michele (about 2 km 2 ). The water dynamics in the HA of higher elevation (shaded area in Figure 1), with an area 57 km 2 , were studied and modeled in a previous work [24].
The hydrological and hydraulic models employed in the present study required a preliminary campaign to acquire the data needed for an accurate description of the transport capacity of the existing drainage network and the hydraulic characterization of the area. Such information included the characteristics of the urbanization, the structure of the drainage networks, the various interconnections, the hydraulic cross-sections, and the state of conservation of the canalizations, as well as the characteristics of the intersections with the other civil infrastructures together with the hydraulic characteristics of the draining systems. This intense campaign covered the entire network of the reclamation channels, a total of 100 channels, and overall involved the survey of • more than 95 km of drainage channels (riverbed and bank profiles), • more than 1250 sections of the remediation channels, • more than 500 structures (bridges, manholes, weirs, etc.), including more than 300 bridges.
The network of the main and secondary channels previously illustrated has a total extension of 97.56 km and is in several cases interconnected.
As previously stated, the study area has been subject to several floods, mainly due to the absence of rain sewers and the insufficient capacity of the reclamation channels and the four major pumping stations. In particular, frequent floods have been reported in the past fifteen years, including the most relevant ones of 1 November 2002, 20 October 2011, 31 January 2014, and the more recent event of 10 September 2017. For all of these intense events, the recording of rainfall intensity (from four stations: Ostia, Acilia, Ponte Galeria, and Isola Sacra, time resolution of 10 min) and the operation data of some or all of the pumping stations (in terms of discharge and water levels in the system) were available; such information was necessary for the calibration and verification of the hydraulic simulations of the events and, for some of them, the detailed description of the hydraulic effects on the territory. Unfortunately, no quantitative information on the piezometry during the flood events was available as no monitoring is routinely performed in the area.

River Flooding and Groundwater Flooding
The joint analysis of rainfall and the water discharge at the pumping stations, which is described and discussed in the next sections, revealed that in most of the events, the runoff component of rainfall alone could not justify the high water volumes that are pumped during and after the rain events. As shown below, the volume of water pumped was often larger than the total average rainfall in the catchment. The reason for this is simple and understandable in view of the nature of the reclamation area, and it relates to the interaction between the drainage channels and the aquifer, such that the hydrological behaviors of the high and low areas are different. The high areas are poorly connected to the aquifer, and the flow discharge in the channels is determined by the surface hydrological component alone (river flooding), which in turn derives from the excess rainfall. In contrast, the low areas are mostly connected to the aquifer, and the channels are fed by both surface (river flooding) and subsurface (groundwater flooding) waters.
Thus, groundwater flooding derives from the increase of the water table in the area, which in turn is due to infiltration. Such increase of the water table is relatively fast because the aquifer is shallow and highly permeable (mostly alluvial sands). The mechanism is described in Figure 2. The hydrological conditions in absence of rain (Figure 2a), which are determined by the water table as lowered by the pumps, change after rainfall as a function of its intensity, which is usually measured in terms of the return period T (Figure 2b). Rainfall determines an increase of the water table which in turn increases the contributing areas near the channels (Figure 2c), leading to a groundwater flooding component that adds to the surface one. This mechanism is in fact very similar to that of the contributing areas in the runoff generation mechanism originally envisioned by Dunne and extensively studied in the literature [25]; because of the low elevation and shallow water table, it represents a fundamental contribution to floods in the region under examination and in similar coastal regions.
Water 2020, 12, x FOR PEER REVIEW 5 of 16 quantitatively similar to or greater than that determined by surface runoff via the standard Hortonian mechanism. The combined effect of flooding due to excess runoff and groundwater can in principle be studied using a coupled hydrological surface-subsurface modeling approach, e.g., as was done in Reference [18]. However, such an approach may easily require significant modeling effort, especially if a two-dimensional approach is pursued; the latter is necessary in view of the complexity of the hydraulic infrastructures, which includes the network of channels, the various structures present (weirs, manholes, small and large bridges, etc.) that have a paramount impact on the dynamics of the floods, and the resulting inundated areas. Additionally, a combined surface-subsurface model would require adequate characterization, especially concerning the groundwater component and the aquifer, which are typically difficult to monitor and characterize.
The objective of the proposed method is to integrate the groundwater contribution to flooding in a parsimonious manner, with few relevant parameters, requiring a minimum of information regarding the hydrogeological setup. The approach is able to be easily implemented and introduced into fully-fledged two-dimensional hydraulic models, where the flood dynamics are modeled along the channel network, together with the pumping sequences of the pumping stations. The approach envisioned here is to model groundwater flooding by a simple groundwater component, along the instantaneous unit hydrograph (IUH) approach, which is distributed along each reach of the channel network. The method requires the determination of the groundwater catchment area drained by the channels, which can be much larger than the surface drainage area and makes use of the same runoff coefficients employed in the excess runoff hydrological component. The method involves only one parameter which is determined after calibration.
In the next section, we detail the methodological aspects of the proposed approach, with application to the case study under examination.

Materials and Methods
We describe in the following the hydrological components providing the inflow to each channel of the network and the 2D hydraulic model of the area; the model includes the pumping stations and their operation as a function of the water levels in the pump chambers.
The hydrological component represents the channel contributions coming both from surface runoff and from the aquifer; it is usually the most critical and uncertain component in the modeling chain [30]. The two contributions, surface and subsurface, reach the drainage channels at different times, quickly in the case of the surface component and more slowly for the groundwater component; The rapidity of groundwater flooding is usually explained by the celerity of water waves in groundwater, which depend on transmissivity and storativity [26], which justifies the rapid streamflow generation from groundwater [27][28][29]. In the case of reclaimed areas, the rapidity of streamflow generation also depends on the density of the channel network that is connected to the aquifer. The combination of these effects may be such that it generates streamflow that is quantitatively similar to or greater than that determined by surface runoff via the standard Hortonian mechanism.
The combined effect of flooding due to excess runoff and groundwater can in principle be studied using a coupled hydrological surface-subsurface modeling approach, e.g., as was done in Reference [18]. However, such an approach may easily require significant modeling effort, especially if a two-dimensional approach is pursued; the latter is necessary in view of the complexity of the hydraulic infrastructures, which includes the network of channels, the various structures present (weirs, manholes, small and large bridges, etc.) that have a paramount impact on the dynamics of the floods, and the resulting inundated areas. Additionally, a combined surface-subsurface model would require adequate characterization, especially concerning the groundwater component and the aquifer, which are typically difficult to monitor and characterize.
The objective of the proposed method is to integrate the groundwater contribution to flooding in a parsimonious manner, with few relevant parameters, requiring a minimum of information regarding the hydrogeological setup. The approach is able to be easily implemented and introduced into fully-fledged two-dimensional hydraulic models, where the flood dynamics are modeled along the channel network, together with the pumping sequences of the pumping stations. The approach envisioned here is to model groundwater flooding by a simple groundwater component, along the instantaneous unit hydrograph (IUH) approach, which is distributed along each reach of the channel network. The method requires the determination of the groundwater catchment area drained by the channels, which can be much larger than the surface drainage area and makes use of the same runoff coefficients employed in the excess runoff hydrological component. The method involves only one parameter which is determined after calibration.
In the next section, we detail the methodological aspects of the proposed approach, with application to the case study under examination.

Materials and Methods
We describe in the following the hydrological components providing the inflow to each channel of the network and the 2D hydraulic model of the area; the model includes the pumping stations and their operation as a function of the water levels in the pump chambers.
The hydrological component represents the channel contributions coming both from surface runoff and from the aquifer; it is usually the most critical and uncertain component in the modeling chain [30]. The two contributions, surface and subsurface, reach the drainage channels at different times, quickly in the case of the surface component and more slowly for the groundwater component; therefore, the two contributions are identified below as the "fast" and "slow" components of the system response, respectively. The fast component was represented by a standard kinematic IUH model; instead, the slow component was modeled using a linear reservoir model, properly adapted to the channel network. While the first component (surface runoff) is rather standard, the second one (groundwater) was novel, at least in the current implementation of flooding in reclaimed areas, in combination with a surface hydraulic model, and represents the major contribution of this paper. The total flow feeding each channel was obtained by adding the slow and fast components. For channels that are not connected to the aquifer (high areas) only the runoff surface component was considered.
Both models, fast and slow, fall into the category of linear rainfall-runoff models, based on IUH theory, which represents the system response (i.e., the flood hydrograph) consequent to an instantaneous unit excess rainfall. Following the IUH approach, the relationship between effective rainfall and flow rate at the closure of the catchment is given by the well-known convolution where q is the discharge at the outlet, A is the area of the (surface) catchment pertaining to each channel, p is the excess rainfall (total rainfall minus the hydrological losses), and u is the IUH. The IUH corresponds to the probability density function of the travel time of water particles from a random location within the catchment to the outlet.

Surface Hydrological Component
The surface (fast) component of flooding, which is determined by the excess water runoff, was calculated using the standard linear kinematic model, such that the IUH u in Equation (1) was univocally determined by the catchment travel time τ b , as follows.
The rainfall input p in Equation (1) is given by either • the rainfall intensity measured during the past flood events for the calibration and verification of the model, or • a synthetic design rainfall, with given return period T, for the assessment of the inundation maps and the flooding scenarios; the design rainfall adopted was the Chicago one [31], which was built on the intensity-duration-frequency (IDF) curves derived from the VAPI project Valutazione delle Piene in Italia, [32] and which provides the IDF for the Italian territory for any choice of T.
The return periods adopted in this study were 10, 30, 50, 200, and 500 years.
The excess rainfall was calculated using the simple runoff coefficient φ, such that the excess rainfall p was equal to the product of the total rainfall by φ < 1. The values of the runoff coefficients were those generally employed in the area of interest [24] and are represented in Table 1. The averaged values were inferred and delimited on the basis of recent aerial photos of the study area, and the runoff coefficient for each contributing basin was therefore calculated as the weighted average based on the fractions of the area occupied by the different types of soil cover. The catchment travel time τ b is the sum of an average travel time to reach the channel under consideration by surface runoff and the travel time through the channel itself. The former is determined based on the size of the drained catchment and a velocity parameter that varies with land cover and use and with the average slope of the catchment, in between 0.25 and 0.5 m/s for HA and 0.01 and 0.035 for LA; the latter is equal to the ratio between the channel length and a travel velocity, which is inferred from the hydraulic model (illustrated later).
The hydrological input, in terms of discharge q(t), is introduced in each of the channels of the network represented in the hydraulic model. The total number of sub-basins considered in the simulations is 265, with an average size of 13 hectars.

Subsurface Hydrological Component
The groundwater (slow) component of channel flow was analyzed using a simple and parsimonious model, as emphasized in the previous sections. Adopting the same linear approach as Equation (1), we modeled the groundwater discharge to the channels using a linear reservoir IUH, leading to the exponential model where K is the storage coefficient. Since the above model reflects the dynamics of the aquifer, the precipitation p that feeds the basin in Equation (1) is that falling over the whole groundwater basin, which generally differs from the surface catchment [33]. For the case under consideration, the extension of the groundwater basin was derived from Reference [34], and it was approximately four times larger the catchment basin of the LA (Figure 3). Thus, the rainfall p was that recharging the aquifer over the groundwater catchment. Consistent with the analysis of runoff, the recharge to the aquifer was obtained by multiplying the total rainfall by the complement to unity of the runoff coefficient, with the same values adopted for the surface component (Table 1). We implicitly neglected the storage in the shallow unsaturated area and evapotranspiration, which are both reasonable assumptions when analyzing intense flood events with high return periods. The storage coefficient K represents the average response time of the slow component in the system; its value was calculated for each sub-basin as the ratio between the half-length of the path (already used to calculate the travel time of the surface, fast response) and the underground flow celerity, not to be confused with velocity, as previously discussed. In the conceptual scheme adopted here, celerity in the subsoil is a parameter that represents the complex physical processes occurring in the subsurface, which was, for many reasons (model complexity, lack of information, etc.), not represented in detail. The groundwater celerity depends on the aquifer characteristics, mainly transmissivity, which in turn is affected by the hydraulic conductivity of the formation and the aquifer thickness and the derived piezometric gradients. The latter two quantities are highly variable during the year and during intense rainfall events. As a consequence, water celerity is bound to vary among events. In order to ensure the robustness of the model, the value of the subsurface celerity was reasonably set as equal for all the basins identified in the study area. Such a constant value was meant to be an average, "effective" value along each rainfall event. In fact, the derived storage coefficient K is bound to change between sub-basins as the average length is different.
pumped and, where available, of recorded levels. In particular, the celerity calibrated for the 2014 event was used for the simulation of the synthetic events, which among all the simulated real events was the one characterized by a higher return period (around 30 years) and occurred at a time of the year (late January) preceded by significant rainfall that contributed to an increase of the aquifer levels; the resulting calibrated celerity was 6.5 × 10 −4 m/s. Adopting such a value for the high scenarios allowed a more precautionary and safe estimation of the inundated areas.

Hydraulic Modeling
The hydrological loads calculated via the above procedure were introduced in all channels of the hydraulic model. The contribution of local drainage networks (sewage), when present, was also considered. The hydraulic model employed was HEC-RAS [35], which is an integrated, 1D-2D model used to solve shallow water equations in an interconnected channel network-floodplain environment.
The 1D calculation domain included the channel network and it covered the channel sections from the bottom to the top of the embankment. As is well known, HEC-RAS is able to model head losses due to the presence of obstructions (bridges, manholes, etc.) or sudden changes in the riverbed sections (weirs, obstructions, etc.), which was a necessary requirement in view of the complexity of the hydraulic system under consideration. The model also included the pump operation, such that the discharge was automatically set as a function of the calculated water levels in the pump chambers, which in turn depended on the water discharge from the channel networks.
The 2D solver employed a parabolic approximation of the two-dimensional shallow water equations, thus neglecting the inertial terms. The hydraulic connections between the 1D and 2D domains were carried out via bidirectional side spillways. About 3600 lateral spillways were introduced and adapted to the bank profiles inferred from the topographic surveys specifically created for this study and verified against the digital elevation map of the Ministry of the Environment (Lidar 1k). Altogether, the following elements were introduced and modeled: 100 The value of the groundwater celerity, and consequently the values of the storage coefficients of the sub-basins, were determined during the calibration of the hydrological and hydraulic model, after reproducing the behavior of the pumping stations during the observed events in terms of volumes pumped and, where available, of recorded levels. In particular, the celerity calibrated for the 2014 event was used for the simulation of the synthetic events, which among all the simulated real events was the one characterized by a higher return period (around 30 years) and occurred at a time of the year (late January) preceded by significant rainfall that contributed to an increase of the aquifer levels; the resulting calibrated celerity was 6.5 × 10 −4 m/s. Adopting such a value for the high T scenarios allowed a more precautionary and safe estimation of the inundated areas.

Hydraulic Modeling
The hydrological loads calculated via the above procedure were introduced in all channels of the hydraulic model. The contribution of local drainage networks (sewage), when present, was also considered. The hydraulic model employed was HEC-RAS [35], which is an integrated, 1D-2D model used to solve shallow water equations in an interconnected channel network-floodplain environment.
The 1D calculation domain included the channel network and it covered the channel sections from the bottom to the top of the embankment. As is well known, HEC-RAS is able to model head losses due to the presence of obstructions (bridges, manholes, etc.) or sudden changes in the riverbed sections (weirs, obstructions, etc.), which was a necessary requirement in view of the complexity of the hydraulic system under consideration. The model also included the pump operation, such that the discharge was automatically set as a function of the calculated water levels in the pump chambers, which in turn depended on the water discharge from the channel networks.
The 2D solver employed a parabolic approximation of the two-dimensional shallow water equations, thus neglecting the inertial terms. The hydraulic connections between the 1D and 2D domains were carried out via bidirectional side spillways. About 3600 lateral spillways were introduced and adapted to the bank profiles inferred from the topographic surveys specifically created for this study and verified against the digital elevation map of the Ministry of the Environment (Lidar 1k). Altogether, the following elements were introduced and modeled: 100 remediation channels, 5800 hydraulic sections, and 500 local structures (bridges, manholes, weirs, etc.), for a total of about 98 km.
The 2D integration domain consisted of mixed-type meshes, both structured (square) of 5 × 5 m mesh, and unstructured (triangles or polygons up to eight faces), depending on the topographical condition. Overall, more than 35 km 2 of land was modeled, made up of approximately 55 2D interconnected areas, or connected to the 1D domain, for a total of approximately 1,870,000 cells (see Figure 4). The pumping stations were also modeled, using the dedicated module available in HEC-RAS.
Water 2020, 12, x FOR PEER REVIEW 9 of 16 remediation channels, 5800 hydraulic sections, and 500 local structures (bridges, manholes, weirs, etc.), for a total of about 98 km. The 2D integration domain consisted of mixed-type meshes, both structured (square) of 5 × 5 m mesh, and unstructured (triangles or polygons up to eight faces), depending on the topographical condition. Overall, more than 35 km 2 of land was modeled, made up of approximately 55 2D interconnected areas, or connected to the 1D domain, for a total of approximately 1,870,000 cells (see Figure 4). The pumping stations were also modeled, using the dedicated module available in HEC-RAS. Manning roughness coefficients in the channels were assessed based on the morphological characteristics of the canals, inferred from the photographs taken along the canals during the survey campaign. Manning coefficients ranging from 0.020 to 0.035 m −1/3 s were adopted for the concretecovered canals, depending on the amount and state of the vegetation present, while for the uncoated riverbeds they were assumed, always based on the vegetation present, to vary between 0.035 and 0.050 m −1/3 s. The boundary conditions downstream were selected as supercritical for the LA, such as to avoid any backwater effect, and the maximum storm surge of the sea level in the case of HA was equal to 0.75 m a.s.l (impact of sea level rise, as considered for instance in Reference [22], was not considered here). A constant Manning coefficient = 0.060 m −1/3 s, was employed in the 2D subdomains of the model. The Manning coefficients and the runoff coefficients were calibrated and verified in the aforementioned previous study focused on the higher areas, after comparison of the model results with the flood areas that actually occurred during the major flood events in the period from 2002 to 2014 [24]. Such values have been adopted by the local reclamation authority.
The input of the model were the discharges calculated with the previously illustrated hydrological analysis. The initial condition was a minimum flow rate in all channels, similar to what is usually found in normal conditions (a few tens of liters per second).

Model Verification against Observed Events
The hydrological-hydraulic model was verified against the available data (mainly rainfall and discharge data from the pumping stations, including levels in the pump chambers) for some of the Manning roughness coefficients in the channels were assessed based on the morphological characteristics of the canals, inferred from the photographs taken along the canals during the survey campaign. Manning coefficients ranging from 0.020 to 0.035 m −1/3 s were adopted for the concrete-covered canals, depending on the amount and state of the vegetation present, while for the uncoated riverbeds they were assumed, always based on the vegetation present, to vary between 0.035 and 0.050 m −1/3 s. The boundary conditions downstream were selected as supercritical for the LA, such as to avoid any backwater effect, and the maximum storm surge of the sea level in the case of HA was equal to 0.75 m a.s.l (impact of sea level rise, as considered for instance in Reference [22], was not considered here). A constant Manning coefficient n = 0.060 m −1/3 s, was employed in the 2D subdomains of the model. The Manning coefficients and the runoff coefficients were calibrated and verified in the aforementioned previous study focused on the higher areas, after comparison of the model results with the flood areas that actually occurred during the major flood events in the period from 2002 to 2014 [24]. Such values have been adopted by the local reclamation authority.
The input of the model were the discharges calculated with the previously illustrated hydrological analysis. The initial condition was a minimum flow rate in all channels, similar to what is usually found in normal conditions (a few tens of liters per second).

Model Verification against Observed Events
The hydrological-hydraulic model was verified against the available data (mainly rainfall and discharge data from the pumping stations, including levels in the pump chambers) for some of the most intense rainfall events that occurred in the study area in the last twenty years. The events included those of 1 November 2002, 20 October 2011, and 31 January 2014, and the more recent event of 10 September 2017; unfortunately, a complete dataset was not available for all the events. The calibration involved only the groundwater celerity, while the other relevant parameters were adopted from a previous study focused on the higher areas, as previously discussed.
The event of 31 January, 2014 (return period T = 30 years) was the most intense and critical event in this area among those examined. Figure 5a Figure 5b shows the temporal behavior of the water levels in the pumping chamber together with the flow rates pumped by the Stagni plant. In all cases, good behavior of the model was observed, with an overall agreement between observations and simulations. The persistence of high levels in the pump chambers, higher than the bank levels of the channels, was due to the limited lifting capacity of the system; this caused significant flooding in the upstream areas, consistent with the photographic documentation of the event, as discussed in the following section.
On the basis of the information available for the 2014 event, it was possible to further check and verify the model capabilities and accuracy for reproducing the dynamics and extension of the inundated areas at the local scale. This was done for several check-points within the drainage network. As an example, Figure 6 display a series of photos taken near the Stagni area during the event, compared with the flood dynamics simulated at the same times as the shots were taken; the comparisons show that the flood advancement and levels compared favorably with what was documented in the field while the event was taking place. The good accuracy of the model in reproducing the propagation dynamics was also confirmed by comparison between simulations and photos along a few relevant channels and in proximity of the Stagni pumping plant (not shown in figure).  Figure 5b shows the temporal behavior of the water levels in the pumping chamber together with the flow rates pumped by the Stagni plant. In all cases, good behavior of the model was observed, with an overall agreement between observations and simulations. The persistence of high levels in the pump chambers, higher than the bank levels of the channels, was due to the limited lifting capacity of the system; this caused significant flooding in the upstream areas, consistent with the photographic documentation of the event, as discussed in the following section.
On the basis of the information available for the 2014 event, it was possible to further check and verify the model capabilities and accuracy for reproducing the dynamics and extension of the inundated areas at the local scale. This was done for several check-points within the drainage network. As an example, Figure 6 display a series of photos taken near the Stagni area during the event, compared with the flood dynamics simulated at the same times as the shots were taken; the comparisons show that the flood advancement and levels compared favorably with what was documented in the field while the event was taking place. The good accuracy of the model in reproducing the propagation dynamics was also confirmed by comparison between simulations and photos along a few relevant channels and in proximity of the Stagni pumping plant (not shown in figure). The numerical model was further verified for the event that occurred in 1 November 2002. The event was less intense than the others, being the one characterized by the smallest return period. Nevertheless, the event caused several floods, mainly in the most depressed areas. Unfortunately, data for this event were available only for the Stagni station. It was seen (Figure 7) that the model behaved quite well in terms of pumped volumes and water level at the pumps. The large volumes pumped compared to the total rainfall confirmed the importance of the groundwater component to this flood. As the total rainfall was smaller than during 2014 event (about 100 mm in 12 h compared to about 150 mm in 24 h), the 2002 event was characterized by a lower level and discharge at the Stagni pumping station; this notwithstanding, the pumped volume at Stagni was greater than the total rainfall, highlighting again the critical role of groundwater flooding. The observed and simulated behaviors also depend on the antecedent conditions (number of rainy antecedent days and the alternation of wet and dry periods).
The importance of the groundwater component in pluvial flooding was further tested by analyzing the event that occurred on 10 September 2017, after a prolonged summer period characterized by an absence of rainfall, such that the groundwater levels in the area were presumably very low. In Figure 8, we show that, about 24 h after the end of the rainfall, the volume pumped at  As the total rainfall was smaller than during 2014 event (about 100 mm in 12 h compared to about 150 mm in 24 h), the 2002 event was characterized by a lower level and discharge at the Stagni pumping station; this notwithstanding, the pumped volume at Stagni was greater than the total rainfall, highlighting again the critical role of groundwater flooding. The observed and simulated behaviors also depend on the antecedent conditions (number of rainy antecedent days and the alternation of wet and dry periods).
The importance of the groundwater component in pluvial flooding was further tested by analyzing the event that occurred on 10 September 2017, after a prolonged summer period characterized by an absence of rainfall, such that the groundwater levels in the area were presumably very low. In Figure 8, we show that, about 24 h after the end of the rainfall, the volume pumped at the Stagni, Ostiense, and Bagnolo pumping stations tended to remain constant and much lower than the total, in contrast to what happened during the previously examined events. This can be explained As the total rainfall was smaller than during 2014 event (about 100 mm in 12 h compared to about 150 mm in 24 h), the 2002 event was characterized by a lower level and discharge at the Stagni pumping station; this notwithstanding, the pumped volume at Stagni was greater than the total rainfall, highlighting again the critical role of groundwater flooding. The observed and simulated behaviors also depend on the antecedent conditions (number of rainy antecedent days and the alternation of wet and dry periods).
The importance of the groundwater component in pluvial flooding was further tested by analyzing the event that occurred on 10 September 2017, after a prolonged summer period characterized by an absence of rainfall, such that the groundwater levels in the area were presumably very low. In Figure 8, we show that, about 24 h after the end of the rainfall, the volume pumped at the Stagni, Ostiense, and Bagnolo pumping stations tended to remain constant and much lower than the total, in contrast to what happened during the previously examined events. This can be explained by the low groundwater levels at the beginning and during the events, confirming again the importance of this component in the flooding dynamics and the water volume pumped. local aquifer features like piezometric surge or groundwater discharge at all locations within the catchment. If such quantities are also of interest, a more integrated approach is needed, including a detailed and distributed analysis of groundwater flow (e.g., References [7,15,18]).

Application of the Model to Synthetic Rainfall Events for Inundated Area Delineation
Thus, the groundwater component is a very important component in the determination of inundated areas in reclaimed zones, like the one under examination here, after intense rainfall events. Neglecting this component may lead to a severe underestimation of water volumes and the extension of inundated areas, and an underestimation of the pumping plants that are necessary to deliver the excess water to the final recipient. We briefly discuss in the following the main results of the pumped water volumes and inundation areas after the considered synthetic rainfall events, as a function of a few return periods. The calibrated and verified 2D numerical model was employed for the simulations, and we applied the groundwater celerity pertaining to the 2014 event, as previously discussed.  Figure 9 shows the rainfall volume (solid lines) and the pumped volumes at the Bagnolo station (dashed lines) as a function of time and the five return periods considered. The great impact of the groundwater component, with volumes always larger than those related to rainwater within the catchment, was clearly visible and consistent with previous results. The differences between solid and dashed lines increased with the return period T, i.e., with the intensity and rarity of the event; in fact, more intense rainfall events lead to a more significant saturation of the surface and groundwater components, leading to a larger runoff and water volumes to be pumped. While the modeling of the groundwater component via a lumped approach approximated the water volumes delivered in the channels and the pumps well, it is worth noting that it cannot capture local aquifer features like piezometric surge or groundwater discharge at all locations within the catchment. If such quantities are also of interest, a more integrated approach is needed, including a detailed and distributed analysis of groundwater flow (e.g., References [7,15,18]).

Application of the Model to Synthetic Rainfall Events for Inundated Area Delineation
Thus, the groundwater component is a very important component in the determination of inundated areas in reclaimed zones, like the one under examination here, after intense rainfall events. Neglecting this component may lead to a severe underestimation of water volumes and the extension of inundated areas, and an underestimation of the pumping plants that are necessary to deliver the excess water to the final recipient. We briefly discuss in the following the main results of the pumped water volumes and inundation areas after the considered synthetic rainfall events, as a function of a few return periods. The calibrated and verified 2D numerical model was employed for the simulations, and we applied the groundwater celerity pertaining to the 2014 event, as previously discussed. Figure 9 shows the rainfall volume (solid lines) and the pumped volumes at the Bagnolo station (dashed lines) as a function of time and the five return periods considered. The great impact of the groundwater component, with volumes always larger than those related to rainwater within the catchment, was clearly visible and consistent with previous results. The differences between solid and dashed lines increased with the return period T, i.e., with the intensity and rarity of the event; in fact, more intense rainfall events lead to a more significant saturation of the surface and groundwater components, leading to a larger runoff and water volumes to be pumped. Unfortunately, the capacity of the existing pumping stations is not adequate to drain the significant water volumes coming from the surface and groundwater hydrological components, resulting in large inundated areas, as a function of the return period. The overlap of the flooding maps for the different return period is displayed in Figure 10. It was observed that the area is prone to flooding even for low to moderate return periods, e.g., = 10 years. This is indeed what usually happens in this area, as explained in Section 2. The differences in the inundated areas reflect the wide range of synthetic rainfall as a function of the return periods considered, ranging from 10 to 500 years. Additionally, the flatter zones generally led to larger inundated areas for the same water volumes. Regarding the latter observation, we remark again that neglecting groundwater contributions severely underestimates these Interestingly, for the larger values of , the inundated area reproduced in the Ostiense/Stagni region almost exactly the Levante lake that was present in the area before the reclamation works started, until the end of the 19th century. Figure 10. Inundated areas as function of the return period T. The map includes the inundated areas evaluated in a previous study [24] (shaded area of Figure 1). Unfortunately, the capacity of the existing pumping stations is not adequate to drain the significant water volumes coming from the surface and groundwater hydrological components, resulting in large inundated areas, as a function of the return period. The overlap of the flooding maps for the different return period is displayed in Figure 10. It was observed that the area is prone to flooding even for low to moderate return periods, e.g., T = 10 years. This is indeed what usually happens in this area, as explained in Section 2. The differences in the inundated areas reflect the wide range of synthetic rainfall as a function of the return periods considered, ranging from 10 to 500 years. Additionally, the flatter zones generally led to larger inundated areas for the same water volumes. Regarding the latter observation, we remark again that neglecting groundwater contributions severely underestimates these Interestingly, for the larger values of T, the inundated area reproduced in the Ostiense/Stagni region almost exactly the Levante lake that was present in the area before the reclamation works started, until the end of the 19th century.
Water 2020, 12, x FOR PEER REVIEW 13 of 16 Figure 9. Synthetic design events: comparison between the pumped water volumes (per unit catchment area) and the total design rainfall at the Bagnolo pumping station for different return periods T.
Unfortunately, the capacity of the existing pumping stations is not adequate to drain the significant water volumes coming from the surface and groundwater hydrological components, resulting in large inundated areas, as a function of the return period. The overlap of the flooding maps for the different return period is displayed in Figure 10. It was observed that the area is prone to flooding even for low to moderate return periods, e.g., = 10 years. This is indeed what usually happens in this area, as explained in Section 2. The differences in the inundated areas reflect the wide range of synthetic rainfall as a function of the return periods considered, ranging from 10 to 500 years. Additionally, the flatter zones generally led to larger inundated areas for the same water volumes. Regarding the latter observation, we remark again that neglecting groundwater contributions severely underestimates these Interestingly, for the larger values of , the inundated area reproduced in the Ostiense/Stagni region almost exactly the Levante lake that was present in the area before the reclamation works started, until the end of the 19th century. Figure 10. Inundated areas as function of the return period T. The map includes the inundated areas evaluated in a previous study [24] (shaded area of Figure 1). Figure 10. Inundated areas as function of the return period T. The map includes the inundated areas evaluated in a previous study [24] (shaded area of Figure 1).

Summary and Conclusions
We analyzed flood and inundation areas in a reclaimed coastal zone, i.e., a low-elevation area where the groundwater level is kept low by a system of reclamation channels and pumping stations. An important component of flooding in such areas is groundwater flooding, complementary to surface runoff, which is determined by the sudden rise of the water table due to intense direct precipitation. When rainfall is intense and the water table is high, the channel network drains both the surface (runoff) water and groundwater components, resulting in a greater pumping effort by the stations. As a consequence, the volume pumped is potentially larger than the total volume of rain that has fallen on the surface catchment, being impacted by the size of the groundwater basin.
Given that most flood models are based on surface hydrological/hydraulic methods and tools, e.g., HEC-RAS, we introduced in this work a simple lumped approach for handling the groundwater contribution to flood in reclaimed areas. The approach is parsimonious (only one calibration parameter), requiring little information on the groundwater setup, and it is easy to implement with existing surface hydraulic models for flood management and the delineation of inundation maps. The method was applied to the analysis of inundation maps of a wide reclamation area located in the southern part of Rome, Italy.
The results obtained confirmed that the simple modeling scheme adopted, although simplified and parsimonious, was able to represent the complex dynamics of the surface water drained by the system of channels with enough accuracy for the scopes at hand. Hence, groundwater can be an important component in the delineation of flood-prone areas, especially for high return periods. Neglecting this component may lead to a severe underestimation of water volumes and the extension of inundated areas, and an underestimation of the pumping plants that are necessary to deliver the excess water to the final recipient.
The proposed approach represents an approximation of the complex groundwater dynamics that occur during and after an intense rainfall event, and it clearly cannot capture local aquifer features like piezometric surge or groundwater discharge at all locations within a catchment. Nevertheless, it may provide a simple and useful procedure to be implemented with existing surface flood models applied to reclaimed areas similar to the one examined in this study.