Numerical Modeling of the Hydro-Morphodynamics of a Distributary Channel of the Po River Delta (Italy) during the Spring 2009 Flood Event

: One-dimensional (1D) numerical models generally provide reliable results when applied to simulate river hydraulics and morphodynamics upstream of the tidal inﬂuence, given the predominantly unidirectional ﬂow conditions. Such models, however, can also be used to reproduce river hydraulics across the ﬂuvial to marine transition zone when speciﬁc conditions occur, as during high discharge events, and the results obtained via these simple modeling tools can provide indicative trends that may guide more structured and detailed modeling of a particularly critical area. In this study, the application of a 1D model setup with hydrologic engineering centers river analysis system (HEC-RAS) for simulating the hydro-morphodynamic conditions of a distributary channel of the Po River Delta (Italy) during a ﬂooding event that occurred in Spring 2009 is presented. The channel bathymetry and the grainsize composition was taken from ﬁeld measurements, while the dimension of the plume o ﬀ shore the delta was derived from a MODIS image acquired at the peak of the ﬂood. The comparison between the numerical outcomes and the ﬁeld evidence shows the reliability of the proposed 1D modeling approach in representing the delta dynamics at a large scale, as well as in showing locations where more spatially detailed studies are needed. The code was also able to adequately reproduce the channel hydro-morphodynamics and the sediment data as derived from a core sample taken a few km o ﬀ shore during the ﬂooding event of April–May 2009. Through a sensitivity analysis, it is also proven that the dimension of the river plume can inﬂuence the evolution of the prodelta, while having a rather negligible e ﬀ ect inland, because of the major stresses induced by the high river discharge during the ﬂood event.


Introduction
Nearly 10% of the world's population lives in coastal areas that are less than ten meters above sea level, along river deltas and adjacent coastlines [1]. Such transitional environments are under threat due to relative sea-level rise and increasing anthropogenic pressure [2] and it is thus fundamental to understand how they will respond to external perturbations [3,4]. To reach this goal, and to adequately simulate future delta dynamics resulting from natural processes and human activities, it is fundamental to quantify their present state and the fundamental mechanisms that control their growth. Likewise, reconstructions of the past evolution of river deltas rely on the knowledge of deltas' morphodynamics that is gained by studying modern systems, as well as by the modeler´s ability to numerically reproduce The annual hydrograph at the gauging station of Pontelagoscuro shows two peaks in discharge, normally in autumn and spring, generated by rainfall and snowmelt, respectively [Error! Reference source not found.]. The mean annual discharge measured at this gauging station during the period 1984-2018 was around 1400-1500 m 3 /s, with a peak discharge higher than 12,000 m 3 /s. The total annual sediment and freshwater yield delivered to the Northern Adriatic Sea were about 13 × 10 9 kg and 40-50 km 3 , respectively [Error! Reference source not found.-Error! Reference source not found.]. Before reaching the Adriatic Sea, the main stem of the Po River forms five major distributary channels within the delta; from north to south, these are the Maistra, Pila, Tolle, Gnocca (or Donzella), and Goro. Among these five, the Pila carries the majority of water and sediment (about 60%), with secondary contributions from the other distributary channels [Error! Reference source not found.,Error! Reference source not found.]. Spring and neap tides in the northern Adriatic Sea have amplitudes of ~86 cm and 20 cm, respectively. At the Pila mouth, the river flow is reversed during flood tide, when the river discharge is less than 1000 m 3 /s [Error! Reference source not found.].

The Spring 2009 Flood
The April-May 2009 river flood was triggered by the combination of spring snowmelt and a significant atmospheric perturbation coming from the Atlantic Ocean [Error! Reference source not found.]. In fact, the winter 2008-2009 was particularly wet, allowing for significant accumulation of snow on the Alps (Figure 1b), which was then released starting from mid-April, because of an The annual hydrograph at the gauging station of Pontelagoscuro shows two peaks in discharge, normally in autumn and spring, generated by rainfall and snowmelt, respectively [21]. The mean annual discharge measured at this gauging station during the period 1984-2018 was around 1400-1500 m 3 /s, with a peak discharge higher than 12,000 m 3 /s. The total annual sediment and freshwater yield delivered to the Northern Adriatic Sea were about 13 × 10 9 kg and 40-50 km 3 , respectively [23][24][25]. Before reaching the Adriatic Sea, the main stem of the Po River forms five major distributary channels within the delta; from north to south, these are the Maistra, Pila, Tolle, Gnocca (or Donzella), and Goro. Among these five, the Pila carries the majority of water and sediment (about 60%), with secondary contributions from the other distributary channels [26,27]. Spring and neap tides in the northern Adriatic Sea have amplitudes of~86 cm and 20 cm, respectively. At the Pila mouth, the river flow is reversed during flood tide, when the river discharge is less than 1000 m 3 /s [12].

The Spring 2009 Flood
The April-May 2009 river flood was triggered by the combination of spring snowmelt and a significant atmospheric perturbation coming from the Atlantic Ocean [28]. In fact, the winter 2008-2009 was particularly wet, allowing for significant accumulation of snow on the Alps (Figure 1b), which was then released starting from mid-April, because of an increment of the temperature in the basin. At the end of April, heavy precipitation occurred in the entire watershed, but mainly concentrated in the western part of the Po River catchment. The co-occurrence of rainfall and snowmelt resulted in increasing the river flow, and the peak discharge recorded at the beginning of May 2009 was significantly over 8000 m 3 /s at Pontelagoscuro (Figure 2). addition, whereas the 2000 Po River flood was a single event that reached a maximum discharge of 9650 m 3 /s at Pontelagoscuro, the 2009 flood occurred at the end of a long rainy period which led to a series of smaller events exceeding 4000 m 3 /s prior to the Spring flood. Furthermore, the 2009 river flood occurred during marine fair weather conditions (wave height < 1 m) and relatively small tidal oscillations (∼50 cm) (Figure 2d,e), suggesting possible deposition in the shallow part of the prodelta, as observed for the October 2000 flood [27].

Numerical Model Setup
To simulate the changes of the water profile and the bed elevation in the period of April 10-May 10, 2009, the well-known code HEC-RAS version 5.0.7 [28] was used, imposing the input data The most notable Po River flood event before the one that occurred in Spring 2009 was the flood of October 2000, the second-largest flood of the last 150 years [27]. As pointed out by Tesi et al. [29], these two major flooding events exhibited interesting contrasts and similarities. In detail, the two floods occurred during different seasons and were triggered by different events: exclusively exceptional

Numerical Model Setup
To simulate the changes of the water profile and the bed elevation in the period of 10 April-10 May 2009, the well-known code HEC-RAS version 5.0.7 [28] was used, imposing the input data reported in Section 2.3.1. The calibration and validation procedures show that the code reproduces in a reliable manner the water levels measured in an intermediate section.
HEC-RAS (www.hec.usace.army.mil/software/hec-ras) allows the user to perform 1D steady flow, 1D and 2D unsteady flow calculations, sediment transport/mobile bed computations, and water temperature/water quality modeling, along a network of natural and man-made channels. In the present application, 1D unsteady flow calculations were run for computing the sediment transport and the morphological changes due to a flooding event, using the data reported below, and following the calibration and validation procedure described in Section 2.3.2. It is worth to note that, being HEC-RAS a 1D model, the bed shear stress values are cross-section averages, computed as a function of the average depth (total flow area divided by the top width) and the slope of the computed local energy grade line.

Input Data: Channel Geometry, Flow Parameters, and Grainsize Composition
The geometry of the Po River reach downstream of Pontelagoscuro was described by combining 57 cross-sections measured in 2005 during a LiDAR campaign (cross-sections available at geoportale.agenziapo.it) with cross-sections of the delta front and prodelta derived from bathymetric data [30] and from satellite images to reproduce the whole plume up to the location of the sediment core E16 (Figure 1c,d). For the offshore part, three sections spacing ca. 1000 m and having a rather regular shape were assumed to adequately represent the plume. Aiming to allow for a better spatial resolution, all the cross-sections were linearly interpolated imposing a maximum separation of 50 m. Different interpolation distances were tested to find a compromise between the spatial accuracy and the computational effort.
The hydrological boundary conditions were defined through the 30 min averaged discharge at the upstream end of the simulated reach (Pontelagoscuro), as derived from the discharge-stage relationship provided by the AIPO (Agenzia Interregionale per il fiume Po)-Po Water Authority. Following a previous report [26], the river discharge was reduced to 60% to account for the reduced transport along the Pila distributary channel. The downstream boundary conditions were set using the water levels measured at the river mouth and supposing that these levels do not change on the prodelta plume. All the water levels were derived from the AIPO website (www.agenziapo.it/content/monitoraggioidrografico-0). To adequately reproduce the hydrological variations, unsteady flow simulations were performed, imposing an integration time-step of 30 s, which was chosen after a sensitivity analysis for code stability.
The bed grain size composition was derived from field samples, collected and analyzed by AIPO in the early 2000s, along the main channel of the Po River and the sediment sample in the prodelta (location E16 in Figure 1c,d) [31]. For computing the sediment transport, the Meyer Peter Müller formula was adopted, imposing equilibrium conditions between sediments and water flow at the upstream end. The reliability of the proposed approach was tested in previous applications made on the same river, also using other solvers [20,32,33].

Model Calibration and Validation
The model was calibrated and validated, in terms of bed roughness, for the period 10 April-10 May 2009, using the first two weeks as a calibration period and the rest for the validation. The 30-minute-averaged computed values were compared to the ones measured every 30 min at the gauging station of Polesella, located around 20 km downstream of Pontelagoscuro ( Figure 1c).
As visible from Figure 3a, the calibration period covered relatively low water levels, while the flooding event characterized the validation timespan. Despite this difference, a roughness Manning coefficient n of 0.020 s/m 1/3 in the main channel provides reliable values for all the water levels observed in April-May 2009, with a slight underestimation of the high water values during the rising limb of both floods (Figure 2b). Such values of roughness agree rather well with literature evidence [1,34].
Being offshore and representing the river plume over the prodelta, for representing the last to km, a very high Manning value of 0.10 s/m 1/3 was assumed. Between the river mouth and the prodelta (around 500 m), the roughness values were linearly interpolated, as made also for the bed geometry. Being offshore and representing the river plume over the prodelta, for representing the last to km, a very high Manning value of 0.10 s/m 1/3 was assumed. Between the river mouth and the prodelta (around 500 m), the roughness values were linearly interpolated, as made also for the bed geometry.

Longitudinal Profile
The numerical model shows that a drawdown M2 profile characterizes the water surface for all discharges during the flood and that water profiles with gentler slopes occur during high tides ( Figure 4). Regardless of the tide, the drawdown profile starts from the Pila mouth and extends for around 15-20 km in the case of high discharges, while during low flow the influence is reduced to only 8-10 km, and the tide is not changing the water surface in a significant manner.

Longitudinal Profile
The numerical model shows that a drawdown M2 profile characterizes the water surface for all discharges during the flood and that water profiles with gentler slopes occur during high tides ( Figure 4). Regardless of the tide, the drawdown profile starts from the Pila mouth and extends for around 15-20 km in the case of high discharges, while during low flow the influence is reduced to only 8-10 km, and the tide is not changing the water surface in a significant manner. Being offshore and representing the river plume over the prodelta, for representing the last to km, a very high Manning value of 0.10 s/m 1/3 was assumed. Between the river mouth and the prodelta (around 500 m), the roughness values were linearly interpolated, as made also for the bed geometry.

Longitudinal Profile
The numerical model shows that a drawdown M2 profile characterizes the water surface for all discharges during the flood and that water profiles with gentler slopes occur during high tides ( Figure 4). Regardless of the tide, the drawdown profile starts from the Pila mouth and extends for around 15-20 km in the case of high discharges, while during low flow the influence is reduced to only 8-10 km, and the tide is not changing the water surface in a significant manner.

Gauging Station at Pila
Generally, at the gauging station of Pila (Figure 1a), the tidal oscillations have a very low influence on the water levels, as reported in Figure 5a. Indeed, here, the system is controlled by the river flow, therefore the discharge is the main driver, as also recognizable looking at the average flow velocity (Figure 5b). Despite the distance between the gauging station and the sea, however, tides can have an influence in the case of relatively low river discharges, as visible from the high water levels measured around April 27 (Figure 5a,b). On the other part, such an oscillation does not force significant changes in terms of bed shear stress and median diameter, given that these quantities are driven by the flow velocity (i.e., by the river discharge). A similar behavior is observable also at the E16 core location, as shown in Figure 6.
April can mobilize sediments having a diameter around 140-150 μm, while, during the flooding event at the beginning of May, the model simulates the transport of particles with smaller diameters. This could be correlated to the presence of finer sediments along the banks, which will be moved in the case of floods, occupying the floodplains.
To summarize, the model points out that, at the Pila gauging station, the hydromorphodynamics are controlled by the river flow in the case of discharges higher than 1000 m 3 /s, in accordance with Maicu et al. [Error! Reference source not found.].

Location of the Core Sample, Cross-Section E16
As reported in Figure 6a, the effects of the tidal oscillations are recognizable in the river discharge computed at the cross-section E16: generally, the lower the discharge, the higher the influence of the tide on the prodelta, even if during flooding conditions. However, given that the water level is mostly controlled by the tidal variations, oscillations of the river discharge have a lower  impact on it. The proposed 1D model computes the depth-averaged velocity along the longitudinal direction (Figure 6b), therefore it is not possible to derive a velocity field from this simple model. Nonetheless, it is possible to recognize the tidal oscillations also from the computed velocities, even if with lower resolution. Lower discharges are associated with a lower stream power and a smaller bed shear stress, resulting in a lower capacity for moving sediments (Figure 6c). Furthermore, the shear stress, and consequently the sediment transport and the erosion/deposition processes, are clearly controlled by the tidal oscillations, even during high flow conditions. For high discharges (flooding conditions), the influence of the downstream boundary conditions (tide) is less significant, with respect to low discharges, therefore these oscillations are smoothed.
The modest shear stress results in low transport capacity, and, therefore, only the finest fractions of sediments can be moved offshore (Figure 6d). In particular, the flooding event contributed to transport of sediments having a median diameter >30 μm, while during low flow conditions, only very fine material can be moved [Error! Reference source not found.].

Variation of the Plume Dimensions
The extension of the plume was derived from a satellite image (250 m resolution) acquired on 2 May 2009 (Figure 1) by MODIS and available through the NASA Earth Observing System Data and Information System (EOSDIS). To test the influence of such dimensions on the final results, two additional simulations were performed assuming smaller or larger plume cross-section, but maintaining the depths as derived from the bathymetry [Error! Reference source not found.]. In designing the new geometries, the section widths were changed of around 30% with respect to the Being correlated to the flow velocity, the bed shear stress has a similar trend (Figure 5c,d), while the active mean diameter has a rather dissimilar behavior. In this case, the discharges observed in April can mobilize sediments having a diameter around 140-150 µm, while, during the flooding event at the beginning of May, the model simulates the transport of particles with smaller diameters. This could be correlated to the presence of finer sediments along the banks, which will be moved in the case of floods, occupying the floodplains.
To summarize, the model points out that, at the Pila gauging station, the hydro-morphodynamics are controlled by the river flow in the case of discharges higher than 1000 m 3 /s, in accordance with Maicu et al. [12].

Location of the Core Sample, Cross-Section E16
As reported in Figure 6a, the effects of the tidal oscillations are recognizable in the river discharge computed at the cross-section E16: generally, the lower the discharge, the higher the influence of the tide on the prodelta, even if during flooding conditions. However, given that the water level is mostly controlled by the tidal variations, oscillations of the river discharge have a lower impact on it. The proposed 1D model computes the depth-averaged velocity along the longitudinal direction (Figure 6b), therefore it is not possible to derive a velocity field from this simple model. Nonetheless, it is possible to recognize the tidal oscillations also from the computed velocities, even if with lower resolution.
Lower discharges are associated with a lower stream power and a smaller bed shear stress, resulting in a lower capacity for moving sediments (Figure 6c). Furthermore, the shear stress, and consequently the sediment transport and the erosion/deposition processes, are clearly controlled by the tidal oscillations, even during high flow conditions. For high discharges (flooding conditions), the influence of the downstream boundary conditions (tide) is less significant, with respect to low discharges, therefore these oscillations are smoothed.
The modest shear stress results in low transport capacity, and, therefore, only the finest fractions of sediments can be moved offshore (Figure 6d). In particular, the flooding event contributed to transport of sediments having a median diameter >30 µm, while during low flow conditions, only very fine material can be moved [31].

Variation of the Plume Dimensions
The extension of the plume was derived from a satellite image (250 m resolution) acquired on 2 May 2009 ( Figure 1) by MODIS and available through the NASA Earth Observing System Data and Information System (EOSDIS). To test the influence of such dimensions on the final results, two additional simulations were performed assuming smaller or larger plume cross-section, but maintaining the depths as derived from the bathymetry [30]. In designing the new geometries, the section widths were changed of around 30% with respect to the values adopted in [31].
The computed values at the Pila mouth ( Figure 7) show that the dimensions of the plume have practically no influence in the inland. There are only small changes in terms of median diameter ( Figure 7b) and bed shear stress (Figure 7c) during the flooding event, while no variations are visible for lower discharges during April 2009.
On the other part, a higher influence of the plume dimensions is recognizable looking at the prodelta. Figure 8 reports four parameters computed at the core location E16, ca. 2 km downstream of the Pila channel mouth (see Figure 1). Comparing the calibrated simulation (black line) with the one performed assuming narrower sections (blue line), it is possible to see that smaller dimensions drive to higher velocities (Figure 8d) and consequently higher bed shear stresses (Figure 8c), which are more prone to move coarser sediments (Figure 8b), especially during flooding conditions. On the contrary, by imposing larger sections (red line), an opposite behavior can be observed.  On the other part, a higher influence of the plume dimensions is recognizable looking at the prodelta. Figure 8 reports four parameters computed at the core location E16, ca. 2 km downstream of the Pila channel mouth (see Figure 1). Comparing the calibrated simulation (black line) with the one performed assuming narrower sections (blue line), it is possible to see that smaller dimensions drive to higher velocities ( Figure 8d) and consequently higher bed shear stresses (Figure 8c), which are more prone to move coarser sediments (Figure 8b), especially during flooding conditions. On the contrary, by imposing larger sections (red line), an opposite behavior can be observed.

Discussion
Being a major source of freshwater, the Po River plays a fundamental role in driving the physical and biogeochemical processes of the Adriatic Sea. For this reason, in the past, several researchers focused on representing the river plume variability, accounting for both the inland drivers (sediments, water discharge) and the sea-related variables (wind, currents, and wave resuspension).
Even if derived from a relatively simple 1D modeling approach, the results reported here are in

Discussion
Being a major source of freshwater, the Po River plays a fundamental role in driving the physical and biogeochemical processes of the Adriatic Sea. For this reason, in the past, several researchers focused on representing the river plume variability, accounting for both the inland drivers (sediments, water discharge) and the sea-related variables (wind, currents, and wave resuspension).
Even if derived from a relatively simple 1D modeling approach, the results reported here are in line with what was computed by Bever et al. [16] using a more complex 3D model. Indeed, by combining the river hydrodynamics with sea currents, waves, and wind, they have shown that the sediment deposition occurred directly offshore of the distributary mouths, with the highest accumulation near the largest source, the Pila distributary, as a combination of the buoyant river plume, wave resuspension, and suspended transport due to the currents of the Adriatic Sea. With respect to previous works [17,35,36], in the present study, the focus was on understanding the general dynamics of a Po River distributary channel. Therefore, gravity flows were not considered here for representing the channel hydro-morphodynamics during the Spring 2009 flood event. Indeed, as pointed out in a previous report [16], except on the steepest portions of the delta front, suspended transport and plume settling are likely to dominate sediment flux during floods, and therefore representing high-density gravity flows is not strictly needed.
Falcieri et al. [37] applied a 3D hydrodynamic model for reproducing the Po River plume pattern variability at different temporal scales, showing that, on long time scales, river discharges represent the dominant forcing in defining the plume size and surface pattern, while on time scales of few days the plume dynamics are modulated mostly by the wind structure. In particular, they demonstrated that, during the period May 2008-July 2009, the river discharges were the dominant driver in shaping the Po River plume, mostly because the weak wind forcing, according to the results reported here.
During high flow conditions, the Po River is a river-dominated system, and the resulting high-water profile can be adequately simulated via the classic backwater profile [4,38], similarly to other fluvial systems worldwide. As an example, Leonardi et al. [39] monitored the flow velocity and the water levels of a distributary channel of the Apalachicola Delta (Florida, USA) by using an acoustic Doppler current profiler (ADCP), pointing out an M2 profile for relatively high tides. Based on such measures, the authors showed that, in this micro-tidal system, a drop in sea level during low tides forces the drawdown of the backwater profile of the river, thus increasing the river flow velocity. The HEC-RAS model applied in the present study shows a similar result, with steeper backwater profiles during low tides (Figure 4b, red line), corresponding to peaking flow velocity.
On the contrary, the code implemented in this application cannot be applied to deltas having very complex geometries or tidal-dominated estuaries like the Fly River [21,40,41], where neap tides generate a bed shear stress higher than the critical one, therefore moving large quantities of sediments that redeposit, preventing the erosion of islands and bars [22].
Focusing on the location of the core sample, as the computed sediment transport is correlated with the velocity, in particular the suspended load, it is possible to conclude that the tidal oscillations have a major influence on the sediment transport during low flow conditions, while, in the case of high flow conditions, it is the river discharge that drives the overall process. Despite small, the oscillations of the river profile due to the tides during the peak of the flood can still affect the sediment transport capacity of the river and the sediment in suspension, leading to the formation of a tidal-modulated flood deposit in the prodelta [31].

Conclusions
This study shows that 1D numerical models can be adopted for describing the hydro-morpho dynamics of a distributary channel during a river flood event, in particular in the case of a fluvial-dominated deltaic system, as the Italian Po River. Using a combination of data acquired at different times and with multiple techniques, this study highlights the reliability of a 1D HEC-RAS model in reproducing the channel hydro-morphodynamics and the sediment data derived from a core sample taken a few km offshore during the flooding event of April-May 2009.
The use of simple codes with a reduced mathematical complexity, combined with data having a rather coarse scale, allows for quick reproduction of observed phenomena, especially in the case of relatively simple geometry, and for describing the effects of small tidal changes on a fluvial-dominated environment, like the Po River mouth. The results obtained with this simplistic approach can be used for implementing more complete models, in particular in the case of environments characterized by complex 2D or even 3D dynamics. In essence, to adequately describe riverine morphodynamics in tidal environments, 1D models can be used to provide the general picture, while providing more complex codes to allow for a detailed description of the involved mechanisms.