Tidal inlet seafloor changes induced by recently built hard structures

Tidal inlets are extremely dynamic environments that are often strongly modified by anthropogenic intervention. In this study, we describe the rapid evolution of a highly human-impacted tidal inlet, studied through repeated high-resolution multibeam surveys and geomorphometric analysis. We document the rapid change induced by new hard coastal structures built to protect the historical city of Venice (Italy). A new breakwater erected between 2011 and 2013 induced the formation of large scour holes with the consequent erosion of about 170 · 103 ± 15.6% m3 of sediment until 2016. The construction of a new island in the middle of the inlet and the restriction of the inlet channel caused a general change of the inlet sedimentary regime from depositional to erosive with a net sediment loss of about 612 · 103 ± 42.7% m3, a reduction of the dune field area by more than 50% in about five years, and a coarsening in the sediment distribution. Our results give new insight on the tidal inlet resilience to changes, distinguishing two different phases in its recent evolution: (i) a very rapid response (from 2011 to 2013) of the seafloor morphology with scour-hole erosion at the new breakwater tips at a rate of about 45⋅103 m3/year and the disappearing of dune fields at a rate of 104⋅103 m2/year; and (ii) a general slowdown of the erosive processes from 2013 to 2016. Nevertheless, the erosion continues at the breakwater, though at a reduced rate, possibly representing a threat to the hard structure. In view of global mean sea level rise and consequent proliferation of hard structures along the coast all over the world, the combined use of very high resolution multibeam surveys and repeatable geomorphometric analysis proposed in this study will be crucial for the monitoring and future management of coastal environments.


Introduction
Coastal systems are among the most productive and yet most threatened ecosystems in the world [1]. In addition, it is estimated that 40% of the world population lives within the coastal a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 The processes induced by the construction of coastal hard structures are well studied in the scientific literature (e.g. [9], and references therein for a review). These studies are mainly based on laboratory experiments (e.g. [20]) and/or numerical models (e.g. [21]).
Shallow-water Multi-Beam Echo Souder System (MBES) has been widely used in mapping seabed morphology and composition (e.g. [ In the last decade, geomorphometric methods to characterize quantitatively terrestrial landscapes (e.g. terrain attributes, feature extraction, automated classification) became well established and have been extensively applied to the marine environment ([42], [43]). The most recent studies based on high resolution repeated multibeam surveys focused on assessment of change in patterns of habitat distribution or on natural spatio-temporal morphological evolution (see for example [ [48]; [49]; [50]; etc.). However, as far as we know, the combination of high resolution data from MBES surveys and geomorphometric and sedimentological analysis has not been used to assess quantitatively and monitor the seafloor changes over time induced by newly built hard defence structures and more generally by coastal infrastructures.
This study investigated the mid-term effects of the construction of hard structures related to the MoSE project, employing a set of three repeated MBES bathymetric surveys, carried out over a period of five years. The aim was to assess the inlet evolution through the monitoring of the seabed sedimentary regime of the inlet system. To achieve this aim, we classified the inlet seafloor morphologies with a repeatable semi-automatic geomorphometric analysis of the digital elevation model of each survey ([43]). Moreover, considering the global increase of hard coastal defence structures to protect the shoreline from global mean sea level rise ( [51]; [52]), our methodology can be applied to other coastal areas around the world for the future management of new coastal defence structures.

Study area
The lagoon of Venice (Fig 1) is located at the northern tip of the Adriatic Sea, along the eastern coast of Italy (45˚N, 12˚E). This is the largest coastal lagoon in the Mediterranean area ( [53]), covering a surface of 550 km 2 (it stretches for 50 km along the coastline, with a mean width of 15 km and an average depth of 1.5 m) ( [54]). The tide is semidiurnal with an average range of 55 cm increasing to 110 cm during spring tides ( [55]) and the residence time varies between 24 hours close to the inlets ( [56]; [57]) and 30 days in the internal lagoon ( [58]; [54]). The lagoon today exchanges water and sediment with the Adriatic Sea through three wide inlets ( [59]): Lido, Malamocco and Chioggia inlet from north to south (Fig 1 left). The total exchange of water with the sea is about 350 � 10 6 m 3 during spring tides and 175 � 10 6 m 3 during neap tides ( [60]).
The Lido inlet (45˚25 0 18@N,12˚26 0 0@E- Fig 1) is the northernmost inlet and is crossed by cruise ship and ferry traffic. The tidal prism at Lido inlet is about 145 � 10 6 m 3 ([61]). Like the rest of the Lagoon, the Lido inlet is exposed to two major wind events; the scirocco, an autumnal/spring wind that blows from south-east, and the bora, that prevails in winter and blows from north-east ( [55]). Until 1882, there were three channels in this area ( [62]), but between 1882 and 1910 the building of the jetties still present today merged these three channels into one and fixed the inlet position ( [63]) (Fig 1 right). In the last 15 years, the inlet configuration has changed again, due to the construction of the MoSE structures (Fig 1 right). The most relevant intervention were the construction of the island (from 2004 to 2008) and of the breakwater (from 2010 to the end of 2012) ( Table 1).  Table 2. Despite some technical differences in the instruments adopted during the successive surveys and in the accuracy of the positioning systems, the three surveys can be reliably compared to support quantitative assessments on the modifications of the sea floor in the study area. The processing of data was the same for all the three datasets: the software CARIS Hydrographic and Side Scan Information Processing System ( [64]) was used to take into account sound velocity variations, tides, and basic quality control of bathymetric data. All corrections were referred to the ZMPS. Backscatter data were processed with the Fledermaus Geocoder Toolbox (FMGT) software. Mosaics and terrain digital models, with a 0.5 m resolution, were created using the software ArcGIS v 10.2 ( [65]).

Sediment sampling
A total of 20 sediment samples were collected with a Van Veen grab in May 2016 and their position was strategically selected based on the main morphological features (Fig 5), according to the 2016 bathymetric map. In this operation a Seapath 300 system supplied by a Fugro HP DGPS was used. The samples were all acquired at slack water, to reduce any positioning error. At the same location of the sediment samples, a Qumox SJ4000 camera installed on an aluminium frame was dropped on the seafloor to shoot underwater images of each sediment sample site. The fine fraction (<1 mm) was analysed by laser diffraction (Mastersizer 3000, Malvern). The measurement range of the analyser is 0.03-1000 μm. The coarse fraction (>1 mm), mainly composed of shell detritus, was analysed with a mechanical sieve. The samples were classified according to the Folk and Ward method ( [66]) using Gradistat statistical package ( [67]) and EntropyMax ( [68]) software. The photoshoots were all visualised with the VLC software to extract a representative image of the bottom.

Geomorphometric analysis
The morphological features of the Lido inlet were identified visually within the ArcGIS platform following the classification of shallow coast landforms of [69]. For all the three datasets we generated contour lines with a spacing of 0.5 m. Bathymetric (with a 50% transparency), hill shade and contour layers were overlapped to help the interpretation. Every feature was classified and inserted in a Geodatabase with its dimensions for further analysis. The classification separated erosive and depositional bedforms and, in particular, identified scour holes, dunes and anthropogenic features on the sea-floor. Scour holes can be defined as localised erosional features generated on a sediment surface by a turbulent current ( [69]; [70]). Dunes are defined as flow-transverse repetitive bedforms that develop when a sediment bed is subjected to a current ( [71]). Anthropogenic features include MoSE main structures, rip-rap areas and dredging marks. A visual classification is subjective and non-repeatable. Therefore, we implemented semiautomatic and repeatable protocols to identify the morphological features for the three datasets. In this way, we could compare the identified morphologies and assess their changes over time.
Scour holes and big dune. Bathymetric Position Index (BPI) from the ArcGis toolbox Benthic Terrain Modeler (BTM) ( [72]) is a second order derivative of the bathymetry ( [73]). The BPI algorithm compares each cell's elevation to the mean elevation of the surrounding cells within a user defined annulus. In the case of crests, a cell will be higher than the surrounding cells in the annulus, giving positive BPI values. Similarly, for depressions the BPI values will be negative ( [74]). The BPI tool was used to automatically identify scour holes and a large dune located in the study area. To run the BPI function, the three bathymetric rasters were resampled at 2 m resolution using the Resample tool. This operation reduced details of rasters but did not affect the broad-scale morphological identification. The BPI inner and outer radius of the annulus were selected for every scour hole for the dataset 3 by comparing the BPI results with the results of the visual classification (Table 3). In this way the BPI algorithm identified concave and convex areas, respectively.
After comparing the BPI values with the results of the visual classification, the broad BPI values were divided into 5 classes (Fig 2 and Table 3): Class A (broad BPI � -3) identified the scour holes S1, S3 and S4; class B (-3 < broad BPI � -1) identified the scour holes S2, S4 and S6; class C (-1 < broad BPI � 0) identified the scour hole S5; class D and E (0 < broad BPI � 2 and broad BPI>2) identified the large dune close to the inlet mouth.
Once the broad BPI class was defined, the following steps were repeated for each dataset: first we applied the Reclassify tool to the BPI raster, keeping only the class values that correctly identified the selected feature. This raster class was then converted into a shapefile by Raster to polygon tool (Fig 2c). Calculate areas was used to calculate the areas of the selected features (see [70]). Scour hole S6 was identified with a semi-automatic procedure: the BPI tool was not able to distinguish between the scour hole and adjacent dredged channel, so this separation was made following the -12 m contour. Table 3. BPI identification protocol with the inner and outer radius for each scour and the resulting classes that identify them: Class A = (broad BPI � -3); class B = (-3 < broad BPI � -1); class C = (-1 < broad BPI � 0); class D = (0 < broad BPI � 2) and class E = (2 < broad BPI). Dune fields. After several tests with different choices of parameters and classes, by comparing with the visual interpretation, we found that the tool that was better suited for the purpose of automatic identification of dune fields was the Vector Ruggedness Measure (VRM or Terrain Ruggedness), that can be found inside the Benthic Terrain Modeler (BTM) toolbox of ArcGIS ( [72]). The VRM calculates the dispersion of vectors perpendicular to each grid cell of the surface and therefore calculates the ruggedness of the sea-floor ( [75]). VRM shows low values both in steep and flat areas, but when the floor is both steep and rugged, its values increase. For the seafloor, VRM measures the complexity of the terrain and, in this way, it highlights the presence of an irregular profile ( [75]).

Morphology
With a bathymetric raster resolution of 0.5 m, VRM was computed by means of a 11 x 11-pixel window (Fig 3). The continuous values from the resulting raster were divided into 2 classes representing minimum (< 0.005) and maximum (>0.005) ruggedness. In this way, we automatically extracted the dune field that was also visually identified (red polygon in Fig 3b).
VRM raster was then reclassified maintaining only the class that represented the dune field. Then the one-class raster was converted into a shapefile with the tool Raster to polygon. From this shapefile, it was possible to measure the dune field area with the tool Calculate areas.
Another two layers were created to help in the bedform interpretation: Slope and Aspect, both deriving from the bathymetric layer (with a 0.5 m resolution).

Volume calculations of bathymetric variations
We calculated the volume differences by integrating the difference of the DTMs over a certain area S. We computed the net sediment volume displaced between any pair of surveys: where V is the volume, S is the area that we analysed, i indicates a grid cell, ΔZ(i) is the difference of depth value for the grid cell i (i.e. the difference between the depth from survey Z 1 (i) and the survey Z 2 (i), in m) and A is the surface area of one grid cell (in m 2 ), i.e. 0.25 m 2 .
The 2016 bathymetric raster lacked data for the scour hole S4 and for other small areas in the whole bathymetry. To overcome this problem, data were interpolated with the tool Focal statistics and so data gaps were substituted using interpolated data.

Error computation
As pointed out by [76], it is crucial to estimate the uncertainty connected with every areal and volumetric measurement. For this estimate, we selected a common reference area for all datasets. The software CARIS calculated the Total Propagated Uncertainty (TPU) for the reference area, both vertical and horizontal for every dataset, with a confidence level of 95%. The error σ z12 related to the difference of two DEMs (e.g. 1 and 2) was computed as the quadratic propagated error of each DEM s z12 ¼ p ðs 2 z1 þ s 2 z2 Þ, where σ z1 and σ z2 are the vertical TPU values of each DEM used to calculate the difference. To calculate the error in the volume measurements, we applied the formula: where N is the total number of pixels contained in the morphological feature under consideration.
For the measurements of areas, the formula applied for the error was: where P is the perimeter of the morphological feature under investigation. The vertical and horizontal TPUs for each dataset are collected in Table 4.

Backscatter analysis
Various methods have been proposed to classify backscatter intensity maps (e.g. [77]; [78]; [47]) in order to identify sub-regions with similar surficial seafloor composition. In this study, given the effectiveness proved by this method in other areas of the Venice lagoon ([40]; [79]), we decided to follow an unsupervised methodology ( [77]) and the implementation of the Jenks' Optimization clustering technique, a tool that can be found in the ArcGis software. Given a user-defined number of classes, the Jenks' algorithm provides the classification of the backscatter intensity map reducing the variance within classes and maximizing the variance between classes ( [80]). We first applied this classification for the 2016 dataset, the only year in which sediment samples were collected.

Results
We visually identified the different morphological features for each dataset ( The presence of a dredged channel is also evident in all three pictures on the south-western side of the inlet, crossing it entirely from the lagoon to the sea. In the north-eastern side of the inlet, instead, the water is shallower (Fig 4 left).
In the classification, we mapped some dredging marks (purple), mainly set close to the MoSE structures (Fig 4 right). The number of scour holes (orange) increased over time: from the 4 depressions visually identified in 2011 (Fig 4d), the number increased to 6 in 2013 and 2016 (Fig 4e and 4f, respectively).
The number of dune fields (red) observable inside the inlet decreased throughout the years (Fig 4 right and S1 Fig); in 2011, 20 dune fields were identified (Fig 4d), in 2013 they fell to 6 and in 2016 only 5 were left (Fig 4e and 4f, respectively). The average properties of each dune field were measured manually and are collected in S1 Table in the Supporting Information. The large dune at the seaward end of the inlet was recorded in all the three surveys (light blue in Fig 4 right). Given that scour holes and dune fields appear as the most dynamic features, we applied the geomorphometric analysis described before to map them semi-automatically: the scour holes and the large dune using the broad BPI; the dune field depicted in Fig 4d (1) using the VRM classification. In fact, we found that the VRM terrain attribute can identify univocally only the dune fields with wavelength (λ) between 5 and 6 m. For this reason, to compare the dune fields recognised with this method, we needed this wavelength to be  Tidal inlet seafloor changes induced by recently built hard structures maintained throughout the whole period (2011-2016). The only dune field that satisfied this was the one represented in Fig 4(1). After the comparison of the inlet morphological features, in order to understand the evolution of the whole tidal inlet over time, we first compared the satellite images of the inlet of 2004 and 2017 (Fig 1 right) to estimate the surface occupied by the new structures (in green in Fig 5). The comparison shows that the surface of the inlet inside the channel was reduced by about 437200 m 2 , whereas the breakwater occupies 34013 m 2 . Then, we qualitatively compared the latest complete low-resolution bathymetry of the Venice Lagoon collected in 2002 ([81]) before the MoSE works even started (before MoSE) and our MBES bathymetric map of 2011 when the construction of the MoSE structures was all finished (after MoSE). The residual bathymetry is shown in Fig 5. The North-Eastern channel, which is shallower than the South-Western challange (average depth of 6 m versus 12 m), was eroded more, as testified by the prevalent blue colour in Fig 5. The trenches corresponding to the lodgement of the mobile barriers were deepened more than 10 m. The area that became deeper (more than 3 m) is close to the rip-rap around the lodgement of the mobile barriers.
Overall, the inlet seems to be in deposition at the northern and southern flanks of the navigation channel, that became deeper.
To quantitatively estimate the more recent changes of the inlet, we compared the DEMs of each MBES survey from 2011 to 2016, assuming that areas with a bathymetric difference outside the error interval can be considered in deposition or erosion, depending on the sign of the difference. The areas with difference values inside the error interval can be considered stable. During the first two years (2011-2013- Fig 6a), we observe a prevalent erosion with an overall net sediment loss of 726.6 � 10 3 ± 46.3% m 3 ( Table 5). The deposition (138.1 � 10 3 ± 39.0% m 3 ) predominantly occurred at the lagoon-side of the newly-built offshore breakwater. From 2013 to 2016 (Fig 6b), instead, the deposition process was dominant (549.3 � 10 3 ± 33.1% m 3 ) with a net sediment gain of 207.3 � 10 3 ± 100.7% m 3 (Table 5). Over about five years, from June 2011 to September 2016, there was a net sediment erosion of 612.2 ± 42.7% m 3 (Fig 6c  and Table 5). Also in this case, like in the period 2002-2011, the North-East channel seems to experience a more severe erosion than the South-West channel.

Scour holes
Scour holes S1, S2 and S3 are located inside the inlet (Fig 4f); scour hole S4 is located at the seaward end of the southern jetty; scour holes S5 and S6 are sited at the south and north ends of the breakwater located outside the inlet, respectively. Fig 7 shows the green, blue and red polygons that identify the scour holes in 2011, 2013 and 2016, respectively. Starting from the DEMs we computed the areas of the scour holes and the volume differences over the years, with relative errors ( Table 6).
The only scour hole whose shape remained constant over the studied interval was S3 ( Fig  7). All the other depressions showed a change in size and shape. In detail, scour hole S5 and S6 did not exist when the first dataset was collected (2011), they appeared in 2013, changing consistently their shape and extent over time. Their areas doubled from 2013 to 2016 (S5: from 32.3 � 10 3 ± 1.5% m 2 to 66.3 � 10 3 ± 0.5% m 2 and S6: from 12.3 � 10 3 ± 2.1% m 2 to   27.6 � 10 3 ± 0.6% m 2 ) ( Table 6). We computed that the newly formed scour holes (S5 and S6) increase their area at a speed of 14215 m 2 /year and 5923 m 2 /year respectively. The differences in volume during the five years show a pronounced erosion of the scour holes S5 and S6 (Table 6), recording a sediment loss of 114.3 � 10 3 ± 12.1% m 3 from 2011 to 2016 in the scour hole S5 and of 58.2 � 10 3 ± 9.9% m 3 in the scour hole S6. The other scour holes (S1, S2, S3 and S4) show irrelevant changes in their volumes, below the estimated error (Table 6).

Dune fields and large dune
The area of dune fields shows a drastic decrease. In 2016, it is less than half of the surface of 2011: from 437.9 � 10 3 ± 0.1% m 2 in 2011 to 256.0 � 10 3 ± 0.2% m 2 in 2013 and, finally, to 212.5 � 10 3 ± 0.1% m 2 in 2016 ( Table 7). The average reduction rate is of 48298 m 2 /year. In particular, the dune field identified with the VRM method (Fig 8a) reduces its area from 122.7 � 10 3 ± 0.2% m 2 in 2011 to 92.4 � 10 3 ± 0.5% m 2 in 2013 and to 32.5 � 10 3 ± 0.3% m 2 in 2016 (Table 7). Its original area shrunk by about 75% at a speed of 19337 m 2 /year. Fig 8b shows the variation in extent of the large dune, identified with the BPI procedure. The large dune stepped backwards and rotated, changing its main axis from West-East to South West-North East from 2011 to 2016. In Table 7 the area and volume differences between surveys are reported. In the period 2011-2013, the large dune area decreases from 27.5 � 10 3 ± 0.7% m 2 to 20.7 � 10 3 ± 1.5% m 2 . In the following three years, it increases from 20.7 � 10 3 ± 1.5% m 2 to 22.7 � 10 3 ± 0.7% m 2 .
The volume differences show no relevant change and remains below the estimated error during the five years (Table 7).  Tidal inlet seafloor changes induced by recently built hard structures one at the lagoon-side of the breakwater and the remaining 11 samples on the dune fields ( Fig  5). The Gradistat analysis classified the samples according to the Folk and Ward diagram ( [66]) as: Sand, Slightly Gravelly Sand, Gravelly Muddy Sand, Gravelly Sand, Muddy Sandy, Sandy Gravel, Gravel and only one as Sandy Mud (S1 Table). To better enhance the grain size variation between samples, a cluster analysis (EntropyMax) was applied. The 20 samples were grouped in 8 groups. Backscatter classification. It was already demonstrated that backscatter intensity is related to sediment particles dimension: coarse sediment is characterised by higher backscatter intensity values, while fine sediment correspond to lower backscatter intensity values ( [82]; [83]; [84]; [26]). By comparing the backscatter distribution and the grain size analysis, we found a good agreement between samples and backscatter classification by rearranging the 8 groups of samples in three classes, corresponding to the three backscatter classess obtained thanks to the Jenks' optimization clustering technique. The three backscatter classes (Fig 9)   • Class II (mainly sandy sediment): this class is characterised predominantly by well sorted sand (72.5%-100% of sand, main mode * 250 μm).

Seafloor sediment distribution
• Class III (gravelly sediment): this class corresponds to an high content of shells (> 33%). In this case the grain size curve has two modes at * 250 μm and * 8000 μm. The analysis of the images shows the seafloor completely covered by shells and shell detritus.
We found that this unsupervised Jenks' classification shows an overall accuracy of 75% ( [85]) obtained deriving the confusion matrix (S2 Table) and counting the percentage of correctly allocated cases ( [86]). The same Jenks' classification was applied by analogy to the 2011 and 2013 datasets. This procedure allowed a qualitative comparison of the classified backscatter images obtained for each survey. Overall we observed the following main changes in the seafloor sediment distribution: at the lagoon-side of the breakwater, from 2011 to 2016, a muddy area has become bigger, while inside the inlet we observed a change in sediment composition, from mainly sandy to mainly gravelly, as shown in Fig 9c).

Discussion
Fontolan et al. [62] defined the Lido inlet as stable, whereas Tambroni and Seminara [87] stated that the inlet was prone to deposition. Helsby [55] demonstrated that from 1930 to 2000 the northern lagoon channels were experiencing higher rates of deposition than of erosion. Lido inlet, in particular, was considered depositional. Defendi et al. [88] measured the total transport in the Lido inlet based on calibrated ADCP measurements in 2005-2006 and by coupling the transport model SEDTRANS96 ( [89], [90]) to the SHYFEM hydrodynamic model ( [91]; [92]), showing a sediment loss of 381 � 10 3 m 3 /year for the lagoon and of 256 � 10 3 m 3 / year for the Lido inlet, with a bed-load transport of 29 � 10 3 m 3 /year. These studies referred to the inlet before the most recent modifications. In Fig 5, we show that in the period between 2002 (before MoSE) and 2011 (after MoSE) in the Lido inlet the deposition was still prevalent, except for the area close to the new island and the navigation channel, that is periodically dredged. The sediment transport observed in [88] could be related mainly to the erosion of these areas.
Using high-resolution MBES, we then quantitatively estimated changes in the time span of almost five years (2011-2016) immediately after the construction of the island in the inner part of the inlet and contemporary to the building of the offshore breakwater (2010-2012).
By comparing the bathymetric survey of 2011 and 2016 (Fig 6c), we found a dominant erosive trend with a net volume loss of 612. However, if we just consider the last interval 2013-2016, the volume of deposited sediment is almost four times higher with respect to the previous interval 2011-2013 and the volume of eroded sediment decreases by 2.5 times (Table 5). This shift from erosion to deposition could be related to the construction of the offshore breakwater, since the seasonal changes in erosion and deposition are unlikely to play an important role: in fact, all the surveys were carried out in the same season (summer). This hypothesis seems to be confirmed by the fact that the deposition concentrates close to the hard-coastal structure (Fig 6). Fig 9 shows that in this area mainly mud deposition occurred.
The construction of the offshore breakwater led also to the formation of the two scour holes (S5 and S6) located at its south and north edges respectively (Fig 4f). These scour holes eroded the ebb-tidal delta, as the breakwater was built on top of it. Ebb-tidal deltas are lobes of sediment that accumulates seaward of the inlet and form because of the interaction of tidal and wave-generated currents ( [94]; [95]). They are usually found in tidal inlets of most barrier island complexes around the globe ( [94]). This feature was completely mapped in the 2011 survey, but not entirely included in the 2013 and 2016 surveys.
Scour holes S5 and S6 increased their maximum depth of more than 2 m over five years. The scouring process is highlighted in Table 6, where the volume differences show a sediment loss both in S5 and S6. The scouring process took place with different velocities, higher for the years 2011-2013, lower for the years 2013-2016: in the two-year-period 2011-2013 the erosion proceeded at 33205 m 3 /year and 12405 m 3 /year for S5 and S6, respectively. In the first-time interval, the breakwater was still under construction, substantially altering the system equilibrium, with consequently very high erosion rates. These rates reduced in the following three years reaching almost half of the previous values (8158 m 3 /year for S5 and 6471 m 3 /year for S6). In 2013 the building of the breakwater was completed, and the scour holes seemed to have adjusted to the new configuration with a consequent slowing down of the erosion process. This is in agreement with observations of the scouring process in correspondence of offshore structures, where there is an exponential deepening of the scours after the structure installation ( [96]; [97]).
Scour holes occurring around breakwaters have been observed globally, like for example in Japan ( [98,99]), in The Netherlands ( [100]) and in the U.S. ( [101]). Processes leading to the formation of scour holes around hard coastal structures have been extensively studied mainly on the basis of tank experiments ([20]; [102]; [103]; [104]). Fredsøe and Sumer ([102]) investigated the scouring at the round head of a rubble-mound breakwater by using regular waves. They found that the major mechanism responsible for the scouring is the formation of leewake vortices in each half period of the waves. The scouring process is governed by the Keulegan-Carpenter number, KC, which depends on the base width of the breakwater head. Larger values of KC imply the formation of larger scour holes. In our case, we determined KC from Eq 6 of [20]: KC = 1 + ((L/1.75B)) 2 , obtaining KC = 1.16 with the width of the breakwater B = 70 m and the width of protection layer (rip-rap) on the seafloor L = 50 m. This value of KC corresponds to a separated flow regime with no horse-shoe-vortex formation in front of the breakwater.
The shape of the scour holes we found is very similar to the scour holes due to non-breaking waves described by [20] at the head of a vertical breakwater. From 2012 to 2016, the monthly averaged wave height and period in the study area ranged from 0.3 to 1.55 m and from 3 to 6.2 s ( [105]). The presence of co-directional currents likely contributed to the wave action enhancing the depth of the scour holes, given that large-scale vortices generated at the breakwater tip can increase the transport capacity of the flow.
We estimated that the ratio d max /h (relative maximum scour depth/water depth) is equal to 0.5 for the scour hole S5 and 0.2 for the scour hole S6. These values were found for wave-dominated scour holes ( [106]). For these reasons, assuming that this ratio will not increase substantially in the near future and being aware that this is a qualitative comparison, it is reasonable to consider that these scour holes are wave-dominated, even though other factors like currents and local bathymetry may play an important role. However, a 3D hydrodynamic modelling analysis would be required to fully understand the role of currents and waves in the scouring process. This process may lead to the gradual dislocation of the rubble mound foundation at the breakwater toe and could progressively endanger the stability of the structures ( [107]; [103]).
Several studies applied geomorphometric semi-automated techniques to describe dune fields parameters as wavelength, height, crest orientation ( [108]; [109]) or to evaluate their migration ( [110]; [111]; [112]). Repeated multibeam bathymetric mapping was used to measure the bed elevation variations caused by dune migration and to estimate the sand transport rates in different parts of the world (e.g. Duffy et al., [110] and [111], in the Bay of Fundy, Canada; Ernsten et al. [113] in the Danish Wadden Sea; Nittrouer et al. [114] in the Mississippi river). More recently, Fraccascia et al. ([45]) studied the morphology and hydrodynamics of the natural tidal inlet of Knudedyb, Danish Wadden Sea, relating the bedforms migration patterns to residual currents. In some cases, geomorphometric analysis of dune fields was crucial to prevent hazards related to human activities (e.g. navigation, construction of pipelines) ( [115]; [112]).
In the Lido Inlet, we classified the dune fields by measuring their average properties from the bathymetry (see S3 Table in Supporting Information). From the backscatter classification (Fig 9) we observed that their composition was mainly gravelly in the throughs and sandy on the crests, in agreement with [79]. The direct comparison of the dune properties, however, was not possible since many dune fields disappeared from one survey to the next. In 2011 there were 20 dune fields, whereas in 2016 only 5 were left, with an overall shrinking of the dune field extent of more than half (51.5%) over about five-years (Table 6).
This shrinking of the dune fields occurred at different speeds: between September 2011 and June 2013, the dune fields diminished at a rate of 103916 m 2 /year; between June 2013 and May 2016 the rate was of 14927 m 2 /year, considerably lower than before. All the surveys were carried out in summer, so it is unlikely that these changes are due to seasonality. They could be related instead to the construction of the island inside the inlet that, as observed before, likely increased the current velocity ( [93]). Higher currents and higher bottom shear stresses possibly induced a rapid reduction of the dune fields. This reduction was accompanied by a coarsening of the sediment in the inlet channel. The same happened in the southernmost inlet of the Venice Lagoon ( [79]. When the breakwater was built (2013), the overall dune fields surface reduction slowed down. For the dune field shown in Fig 8a, however, the reduction speed increased from 17309 m 2 /year in the interval 2011-2013, to 20554 m 2 /year in the following three years. At this rate of reduction, this dune field will eventually disappear in less than two years.
The disappearance of the dunes can be explained by the change of shape of the inlet related to the MoSE defence: a) the construction of the island at the land-ward tip of the inlet channel determined a bifurcation of the flux entering the lagoon and a confluence of two channel exiting the lagoon with probably high current velocities that contributed to the erosion of the dune fields close to the mobile barrier structures and in the centre of the inlet; b) the construction of the breakwater determined another flow separation: part of the flow deviates south of the inlet mouth while the other part continues straight.
Consequently, the dunes reduction decelerated, and the orientation of the dune fields changed over time (Figs 4 and 8). In particular, the dune field immediately outside the inlet mouth (indicated by the number 2 in Fig 4d), in 2011 extended parallel to the jetties as was the flux without the breakwater and was about to disappear in 2013. After the construction of the breakwater it was recreated in a south-west direction, with the result that in 2016 a dune field area is still present with an orientation rotated by more than 125˚clockwise with respect to the one of 2011.
The general shrinking of the dune fields suggests that the sediment transport regime should have changed radically with the completion of the construction of the artificial island for the MOSE project with an enhanced seaward flux of sand. This is evident by comparing the sediment distribution over time (Fig 9). This is in agreement with the findings of [116], that analyzed the tidal components in the Northern Adriatic Sea and in the Venice Lagoon during the last 70 years. They observed an increase amplitude of the major tidal components and a shift of the Venice Lagoon tidal asymmetry towards ebb dominance. Particularly, in the last few years, they observed that the recent reduction of the inlets cross-sectional area further enhanced the ebb dominance over the whole lagoon and probably increased the seaward flux of sand, as well. The sand flux through the inlets is dominated by bed-load transport as it was shown by the direct measurements in the Lido Inlet in 2006 (before the construction of the island), by [117]. The disappearing of dune fields may imply the overall reduction of sand transport inside the inlet.

Conclusion
Tidal inlets are extremely dynamic environments governed by natural processes which control their morphological variations. Often these changes represent a problem to human activities and many actions are taken to limit this natural phenomenon. This is the case of the Lido inlet (Venice, Italy) that went through several human-induced modifications in the last century.
Our In detail: 1. Two new scour holes formed at the tips of the recurved breakwater positioned on the seaward side of the inlet. From the comparison with literature on scour holes at breakwaters, we found that these scour holes are more likely wave-dominated than current dominated. During this 5-years research, both these depressions deepened and widened: S5 area went from 32.3 � 10 3 ± 1.5% m 2 to 66.3 � 10 3 ± 0.5% m 2 ; S6 area from 12.3 � 10 3 ± 2.1% m 2 to 27.6 � 10 3 ± 0.6% m 2 . Scouring took place with two different speeds: faster in the first two year, slower in the following three years. These different velocities are explained by the fact that in the first two-year-period the breakwater was under construction, probably leading to drastic and rapid modifications of the seafloor. In the following three years the breakwater was fully built, and the rate of system alteration slowed down. To fully understand whether this slower variation implies a return of the system to a different equilibrium, it will be necessary to plan new surveys over the same area and to compare the results with 3D hydrodynamic models including the effect of currents and waves. However, if the scouring will not stop, the new-built structures are probably going to face problems connected with their structural stability.
2. Dune fields drastically shrunk, going from an area of 437.9 � 10 3 ± 0.1% m 2 in 2011 to one of 212.5 � 10 3 ± 0.1% m 2 in 2016, with a reduction rate of 103916 m 2 /year in the first twoyear-period. In the following three years the mean rate of disappearance decreased to 14927 m 2 /year. The reduction of dune fields corresponds to the reduction of the sand content in the seafloor sediment composition in favour of shell and shell detritus.
In view of global mean sea level rise and probable future construction of new "hard" defences against it, the combined use of very high resolution multibeam surveys and repeatable geomorphometric analysis proves to be an effective tool for knowledge-based coastal monitoring and management.