Statistical characterization of erosion and sediment transport mechanics in shallow tidal environments – Part 1: Erosion dynamics

Reliable descriptions of erosion events are foundational to effective frameworks relevant to the fate of tidal landscape evolution. Besides the rhythmic, predictable action of tidal currents, erosion in shallow tidal environments is strongly influenced by the stochastic wave-induced bottom shear stress (BSS), mainly responsible for sediment resuspension on tidal flats. However, the absence of sufficiently long, measured time series of BSS prevents a direct analysis of the combined tide- and wave-driven erosion dynamics and its proper representation in long-term morphodynamic models. Here we test the hypothesis of describing erosion dynamics in shallow tidal environments as a Poisson process by analysing, with the peak-over-threshold theory, the BSS time series computed using a fully coupled, bi-dimensional numerical model. We perform this analysis on the Venice Lagoon, Italy, taking advantage of several historical surveys done in the last 4 centuries, which allow us to investigate the effects of morphological modifications on spatial and temporal erosion patterns. Our analysis suggests that erosion events on intertidal flats can effectively be modelled as a marked Poisson process in different morphological configurations because the interarrival times, durations, and intensities of the over-threshold exceedances are always well described by exponentially distributed random variables. The resulting statistical characterization allows a straightforward computation of morphological indicators, such as the erosion work, and paves the way for a novel synthetic, yet reliable, approach for the long-term morphodynamic modelling of tidal environments.

1 Introduction

Understanding the morphodynamic evolution of shallow tidal systems is crucial to suggest long-term sustainable management strategies under the increasing pressure of relative sea level rise (Nicholls et al., 2021) and human interventions (Tognin et al., 2021; Orton et al., 2023). Several process-based models have been developed to describe the morphodynamics of shallow tidal basins. Although these models were originally developed to deal with short-term timescales, various techniques were proposed to accelerate bed evolution and upscale the results at much longer timescales. The most commonly adopted upscaling techniques are the tide-averaging procedure and the online update with a morphological factor (Latteux, 1995; Roelvink, 2006). The tide-averaging approach assumes that the timescale at which morphological changes take place is far longer than that typical of hydrodynamic changes, and therefore, the bottom can reasonably be considered fixed over a tidal cycle so that bed level changes can be upscaled offline from gradients in the tidally averaged sediment transport (e.g. Lanzoni and Seminara, 2002; Sgarabotto et al., 2021). On the contrary, to continuously update the morphology describing short-term interactions between hydrodynamics and sediment transport (e.g. Defina, 2003; Lesser et al., 2004; Carniello et al., 2012) but, at the same time, to overcome the prohibitively time-consuming computation associated with the direct upscaling of this fully coupled approach, the so-called “morphological factor” can be adopted (Lesser et al., 2004; Roelvink, 2006). The morphological factor N is essentially an online multiplication factor for the bed change so that the results of a simulation over one tidal cycle in fact depict the morphological changes over N cycles. Fully coupled models upscaled using a morphological factor are increasingly adopted to investigate not only decadal (e.g. Gourgue et al., 2022) but also centennial to millennial morphological evolution of tidal and estuarine systems (e.g. Braat et al., 2017; Pinton et al., 2023; Baar et al., 2023).

These upscaling techniques are based on the underlying assumption that the actual morphological evolution of a system is equivalent to that resulting from a repetitive pattern representative of the dominant forcing conditions. This hypothesis is reasonable when the main hydrodynamic forcing is represented by tidal oscillation, although attention should also be paid when selecting representative boundary conditions for astronomic tidal patterns (Schrijvershof et al., 2023). Instead, also taking into account merely stochastic processes, such as wind waves and storm surges, is far less straightforward. When wind climate may be reduced to a limited set of representative conditions, multiple simulations can be run and then a weighted average of the different results can be determined to estimate the upscaled morphological evolution (Roelvink, 2006). However, when representative wind and storm climate cannot be reduced to a limited batch of boundary conditions or, more importantly, not only the magnitude but also the temporal succession of these events is likely to strongly affect the morphological evolution of the system, then these upscaling techniques cannot be properly applied.

A different perspective would be to directly consider the stochasticity of morphodynamic processes. From this point of view, the first step is to test the possibility of setting up a statistically based framework in order to generate synthetic, yet reliable, time series to model the morphodynamic evolution on long-term timescales and compare possible scenarios in a computationally effective way through the use of independent Monte Carlo realizations. Although the statistical characterization of the long-term behaviour of several geophysical processes is becoming increasingly popular in hydrology and geomorphology (e.g. Rodriguez-Iturbe et al., 1987; D'Odorico and Fagherazzi, 2003; Botter et al., 2007; Park et al., 2014), applications to tidal landscapes are still quite rare (D'Alpaos et al., 2013; Carniello et al., 2016).

The morphological evolution of tidal systems can be described by Exner's equation:

$$\begin{array}{}\text{(1)}& (\mathrm{1}-n){\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial t}}+\mathrm{\nabla }{q}_{\mathrm{b}}=D-E,\end{array}$$

where n is the bed porosity, z_b is the bed elevation, q_b is the bedload, and D and E are the deposition and entrainment rates of sediment, respectively. In mud-dominated tidal systems, sediment is primarily transported in suspension, and the bedload is negligible; hence, the bed level changes can be determined by accurately describing erosion and deposition. Erosion, E, is directly influenced by the local bottom shear stress (BSS), which results from the interaction between tidal currents and wind waves in shallow tidal systems (Green and Coco, 2014). Instead, deposition, D, is linked to the suspended sediment concentration (SSC). However, SSC is largely affected by advection and dispersion processes at a larger scale and, therefore, cannot be solely determined by local resuspension. Consequently, to effectively model bed level variations, it is essential to accurately describe both BSS and SSC. This contribution focuses on characterizing BSS, while the analysis of SSC is presented in the companion paper (Tognin et al., 2024).

The statistical characterization of erosion events requires a sufficiently long BSS time series. For this purpose, direct, continuous BSS measurements are rarely available, and thus, process-based numerical models can represent a useful tool to provide sufficiently long and spatially distributed time series to characterize erosion dynamics at the basin scale. Moreover, to test the hypothesis that such a simplified approach holds independently on the specific morphological configuration, this characterization should ideally be performed in different tidal systems or in a single tidal system, where information on the morphological evolution is available for a sufficiently long time. From this point of view, the Venice Lagoon, Italy (Fig. 1), represents an almost unique opportunity to test long-term models because several bathymetric surveys are available for the last 4 centuries (Carniello et al., 2009; D'Alpaos, 2010a; Tommasini et al., 2019; Finotello et al., 2023).

Figure 1 Morphological features and wind conditions characterizing the Venice Lagoon. (a) Spatial distribution of the morphological features characterizing the Venice Lagoon. The locations of the anemometric (Chioggia) and oceanographic (CNR Oceanographic Platform) stations are also shown, together with the locations of the three stations at the inlets (SL, SM, and SC) and two stations (S1 and S2) for which we provide a detailed statistical characterization of over-threshold events. (b) Wind rose for the data recorded at the Chioggia station in 2005. The dashed line shows the wind rose for the period 2000–2020.

Toward the goal of developing a synthetic theoretical framework to represent wind-wave-induced erosion events and accounting for their influence on the long-term morphodynamic evolution of tidal systems, we applied a two-dimensional finite-element model for reproducing and analysing the combined effect of wind waves and tidal currents in generating BSSs in several historical configurations of the Venice Lagoon. More in detail, in the present study, we used the fully coupled Wind Wave–Tidal Model (WWTM) (Carniello et al., 2005, 2011) to investigate the hydrodynamic behaviour of six historical configurations of the Venice Lagoon, namely 1611, 1810, 1901, 1932, 1970, and 2012. For each configuration, we run a 1-year-long simulation, considering representative tidal and meteorological boundary conditions. The resulting spatial and temporal dynamics of BSSs for the six selected configurations were analysed on the basis of the peak-over-threshold (POT) theory once a critical shear stress for bed sediment erosion was chosen.

The main goal of this analysis is to find whether, in line with previous results on the present-day configuration of the Venice Lagoon (D'Alpaos et al., 2013), wave-induced erosion events can also be modelled as marked Poisson processes in different morphological settings. The relevance of this result lies in the possibility to describe erosion processes as a Poisson process over time, which represents a promising framework for long-term studies. Our analyses provide a temporally and spatially explicit characterization of wind-induced erosion events for the Venice Lagoon, starting from the beginning of the 17th century, thus allowing us to investigate and understand the main features of the erosive trends the lagoon has been experiencing and to provide predictions on future scenarios.

2 Materials and methods

2.1 Geomorphological setting

The Venice Lagoon, located in the northern Adriatic Sea and characterized by an area of 550 km², formed over the last 7500 years covering alluvial Late Pleistocene, silty–clayey deposits, locally known as caranto (Zecchin et al., 2008). In the present-day morphology, the lagoon is connected to the sea with three inlets, namely Lido, Malamocco, and Chioggia (Fig. 1), through which the tide propagates within the back-barrier system. The tidal regime is semidiurnal, with a maximum tidal amplitude of about 0.75 m, which is typical of the northern Adriatic Sea (D'Alpaos et al., 2013; Valle-Levinson et al., 2021).

Figure 2 Historical bathymetries of the Venice Lagoon. Colour-coded bathymetries of the six different configurations of the Venice Lagoon: (a) 1611, (b) 1810, (c) 1901, (d) 1932, (e) 1970, and (f) 2012.

Meteorological conditions also importantly affect the hydrodynamics of the Venice Lagoon. In particular, storm surges generated by the southeasterly sirocco wind (Fig. 1b) often overlap astronomical tides, thus leading to a large increase in water levels (Mel et al., 2014). Instead, the northeasterly bora wind (Fig. 1b) is mainly responsible for the generation of wind waves and water level set-up, especially in the central and southern portions of the lagoon (Carniello et al., 2009).

The morphology of the lagoon has changed deeply over the last 4 centuries (Fig. 2), especially owing to anthropogenic modifications (Carniello et al., 2009; D'Alpaos, 2010a; Finotello et al., 2023). By the end of the 16th century, all the major rivers flowing into the lagoon were diverted to debouch directly into the open sea, thus dramatically decreasing fluvial sediment input. Furthermore, between the 1900s and 1950s, the exploitation of groundwater for industrial purposes accelerated the local subsidence rates, with anthropogenically induced subsidence reaching values of about 10 to 14 cm in the area of the city of Venice (Carbognin et al., 2004; Zanchettin et al., 2021). During the same period, the total area open to the propagation of tides was largely reduced, due to extensive land reclamation carried out to accommodate industrial, agricultural, and aquacultural activities, especially along the landward margin of the lagoon (Fig. 2c–e). On the seaward side, massive jetties were built between 1839 and 1934 to stabilize the sections of the three inlets and maintain the water depths required for increasingly bigger commercial ships (Fig. 2c and d). For the same reason, navigation channels were excavated in the central part of the lagoon to connect the inner harbour with the sea (Fig. 2e). As a result, these human interventions, together with eustatic sea level rise (average value 1.23±0.13 mm yr⁻¹ between 1872 and 2019; 2.76±1.75 mm yr⁻¹ between 1993 and 2019; see Zanchettin et al., 2021), played a primary role in affecting the morphological evolution of the lagoon. Only in the northern lagoon was the morphological degradation less pronounced because this area is sheltered by the mainland against the northeasterly bora wind and because the human pressure in that area was less intense than in other areas of the lagoon. Therefore, the northern basin displays, also in the present-day configuration, relatively shallow intertidal flats and larger salt marsh areas compared to the central and southern lagoon (Fig. 2f).

We considered six different configurations of the Venice Lagoon (Figs. 2 and S1 in the Supplement), covering a time span of 4 centuries, in order to assess the evolution through time of the feedback mechanisms between morphology and wave-induced erosion. The oldest three configurations (1611, 1810, and 1901) were reconstructed by using historical maps, while the more recent ones make use of the topographic surveys carried out by the Venice Water Authority (Magistrato alle Acque di Venezia) in 1932, 1970, and 2003. The updated description of the more recent morphological modifications, which mainly occurred at the three inlets in the context of the Mo.S.E. project for the safeguarding of the city of Venice by high tides (almost completed in 2012), was included in the 2003 configuration, so we refer to the latter configuration as to the 2012 configuration.

In each bathymetry and, hence, in each computational grid, bed elevation refers to the local mean sea level at the time when each survey was performed. We refer the reader to Tommasini et al. (2019) for a detailed description of the methodology adopted to reconstruct the historical configurations and for information on the bathymetric data of the Venice Lagoon. The computational grids reproducing all six of the considered configurations of the lagoon are shown in Fig. S1 and were calibrated in previous studies, namely in 1611 by Tommasini et al. (2019); in 1810 by D'Alpaos and Martini (2005) and D'Alpaos (2010b); and in 1901, 1932, 1970 and 2012 by Carniello et al. (2009) and Finotello et al. (2023).

2.2 Numerical model and simulations

To compute the hydrodynamic and the wind wave fields in the six selected configurations of the Venice Lagoon, we used the two-dimensional (2-D) fully coupled Wind Wave–Tidal Model (WWTM) (Carniello et al., 2005, 2011). The numerical model, coupling a hydrodynamic module and a wind wave module, describes the hydrodynamic flow field, together with the generation and propagation of wind waves using the same computational grid.

The hydrodynamic module uses a semi-implicit staggered finite-element method, based on Galerkin's approach to solving the 2-D shallow water equations, suitably rewritten in order to deal with partially wet and morphologically irregular domains (Defina, 2000; Martini et al., 2004, Supplementary Information). The bottom shear stress induced by currents, τ_tc, is evaluated using the Strickler equation, considering the case of a turbulent flow over a rough wall, which can be written as (Defina, 2000)

$$\begin{array}{}\text{(2)}& {\mathit{\tau }}_{\mathrm{tc}}=\mathit{\rho }gY\left({\displaystyle \frac{\left|q\right|}{K_{\mathrm{s}}^{\mathrm{2}}{H}^{\mathrm{10}/\mathrm{3}}}}\right)\mathit{q},\end{array}$$

where ρ is water density; g is the gravity acceleration; Y is the effective water depth, defined as the volume of water per unit area actually ponding the bottom; q is the flow rate per unit width; K_s is the Strickler roughness coefficient; and H is an equivalent water depth that accounts for the typical height of ground irregularities. The hydrodynamic module also provides the water levels that are used by the wind wave module to assess the wave group celerity and the bottom influence on wind wave propagation.

For the wind wave module (Carniello et al., 2005, 2011), the wave action conservation equation is parameterized using the zero-order moment of the wave action spectrum in the frequency domain (Holthuijsen et al., 1989, Supplement). An empirical correlation function relating the peak wave period to the local wind speed and water depth determines the spatial and temporal distribution of the wave period (Young and Verhagen, 1996; Breugem and Holthuijsen, 2007; Carniello et al., 2011). The wind wave module computes the bottom shear stress induced by wind waves as (Carniello et al., 2005)

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{\mathrm{ww}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{f}_{\mathrm{w}}{u}_{\mathrm{m}}^{\mathrm{2}},\end{array}$$

where u_m is the maximum horizontal orbital velocity associated with wave propagation, and f_w is the wave friction factor. According to the linear theory, the bottom velocity u_m can be evaluated as

$$\begin{array}{}\text{(4)}& u_{\mathrm{m}}={\displaystyle \frac{\mathit{\pi }{H}_{\mathrm{w}}}{T\sinh \left(kh\right)}},\end{array}$$

where H_w is the wave height, T denotes the wave period, k is the wave number, and h is the water depth. The wave friction factor can be approximated as (Soulsby, 1997)

$$\begin{array}{}\text{(5)}& f_{\mathrm{w}}=\mathrm{1.39}{\left[{\displaystyle \frac{u_{\mathrm{m}}T}{\mathrm{2}\mathit{\pi }\left(D_{\mathrm{50}}/\mathrm{10}\right)}}\right]}^{-\mathrm{0.52}},\end{array}$$

where D₅₀ is the median grain size.

The total bottom shear stress, τ_wc, resulting from the combined effect of tidal currents and wind waves, is enhanced beyond the sum of the two contributions because of the non-linear interaction between the wave and the current boundary layer. In the WWTM, this is accounted for by using the empirical formulation, suggested by Soulsby (1995, 1997), as follows:

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{\mathrm{wc}}={\mathit{\tau }}_{\mathrm{tc}}+{\mathit{\tau }}_{\mathrm{ww}}\left[\mathrm{1}+\mathrm{1.2}\left({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{ww}}}{{\mathit{\tau }}_{\mathrm{ww}}+{\mathit{\tau }}_{\mathrm{tc}}}}\right)\right].\end{array}$$

Even if BSSs induced by the tidal currents are typically smaller than those produced by wind waves, they are of fundamental importance in modulating the temporal evolution of the total BSSs and can increase the peak BSS values by up to 30 % (Mariotti et al., 2010; D'Alpaos et al., 2013).

The WWTM has been widely tested against field observations not only in the Venice Lagoon (e.g. Carniello et al., 2005, 2011; Tognin et al., 2022a) but also in other shallow microtidal environments worldwide, for example, in the back-barrier lagoons of the Virginia Coast Reserve (Mariotti et al., 2010) and the Cádiz Bay (Zarzuelo et al., 2018, 2019). Concerning the Venice Lagoon, model calibration and testing have been performed only in the most recent configuration, i.e. when field data are available. For the older configurations of the lagoon, no hydrodynamic measurements are available, and consequently, the roughness coefficient values have been derived analogously to those selected for the calibrated grid in the most recent configuration by comparing local sediment grain size, bed elevation, and morphological classes (e.g. channel, tidal flat, and salt marsh), which also take into account the possible presence of vegetation (Finotello et al., 2023).

We summarize here the model performance in reproducing tidal levels, significant wave heights, and flow rates at the inlets by reporting the standard Nash–Sutcliffe model efficiency (NSE) parameter computed when field data are available and refer the interested reader to the Supplement (Figs. S2–S5) and the literature (Carniello et al., 2005, 2011; Tognin et al., 2022a) for a more detailed comparison. Adopting the classification proposed by Allen et al. (2007), the model performance can be rated from excellent to poor (i.e. NSE > 0.65 excellent; 0.5 < NSE < 0.65 very good; 0.2 < NSE < 0.5 good; NSE < 0.2 poor). The WWTM model is excellent in reproducing tidal levels (NSEmean=0.970; NSE_median=0.984; NSE_SD=0.040), very good to excellent in reproducing significant wave heights (NSE_mean=0.627; NSE_median=0.756; NSE_SD=0.357), and excellent in replicating flow rates at the inlets (NSE_mean=0.853; NSE_median=0.931; NSE_SD=0.184) (statistics are derived from the calibration reported in Carniello et al., 2011, their Tables 1–3, and Tognin et al., 2022a, their Table S2).

We applied the numerical model to the six computational domains representing the Venice Lagoon and a portion of the Adriatic Sea in front of it in order to perform 1-year-long simulations (Fig. S1). The boundary conditions of the model are the hourly tidal levels measured at the Consiglio Nazionale delle Ricerche (CNR) Oceanographic Platform, located in the Adriatic Sea approximately 15 km offshore of the coastline, and wind velocities and directions are recorded at the Chioggia anemometric station, for which a quite long data set was available (Fig. 1a). In particular, we selected as boundary conditions the measurements recorded in 2005 because they can be considered representative of the wind climate of the Venice Lagoon. Indeed, the probability distribution of wind speeds in 2005 shows the minimum difference with the mean probability distribution computed for the period 2000–2020, compared to any other year within the same time interval (Fig. S6 and Table S1 in the Supplement). Indeed, a visual comparison between the wind roses for 2005 and the entire period 2000–2020 supports this choice (Fig. 1b). Because bed elevation in each computational grid refers to the coeval mean sea level, by using the same water levels as boundary conditions, we implicitly take into account the effect of historical relative sea level variations. Considering the same wind and tidal forcing in each historical configuration of the Venice Lagoon allowed us to isolate the effects of the morphology on the hydrodynamics and wind wave fields.

2.3 Peak-over-threshold analysis of BSS

The morphodynamic evolution of tidal environments is controlled by the complex interaction among hydrodynamic, biologic, and geomorphologic processes, which include both deterministic and stochastic components. As an example, it was shown that sediment transport dynamics in the Venice Lagoon is mostly linked to some limited and severe events induced by wind waves (Carniello et al., 2011), whose dynamics are markedly stochastic in the present-day configuration of the lagoon (D'Alpaos et al., 2013; Carniello et al., 2016). In this work, at any location within each considered configuration of the Venice Lagoon, we used the peak-over-threshold (POT) theory (Balkema and de Haan, 1974) to analyse the temporal and spatial evolution of the total BSS, τ_wc.

In general, the selection of the threshold for the POT method must satisfy two contrasting requirements. On the one hand, the threshold must be large enough to discern stochastic events from the deterministic background. On the other hand, the threshold should not be too high to avoid the loss of important information and the need for a much longer time series to compute meaningful statistics because of the lower number of threshold exceedances. Moreover, the extreme value theory postulates the general emergence of Poisson processes whenever the censoring threshold is high enough (Cramér and Leadbetter, 1967). To comply with these requirements, in the present analysis, the threshold is maintained well below the maximum observed values in order to remove only the background modulation induced by tidal currents without losing significant information on the stochastic wave-driven erosion process.

In applying the POT method to BSS time series, setting the threshold equal to a critical BSS for erosion, τ_c, presents the non-trivial advantage of also preserving the physical meaning of the erosion mechanism. Values of critical BSS for erosion for fine, cohesive mixtures typical of shallow tidal settings largely vary in the literature and are affected by multiple physical and biotic factors (Mehta et al., 1989). Erosion shear stress from in situ measurements on the tidal flats of the Venice Lagoon ranges between 0.2 and 2.3 Pa (0.7±0.5 Pa − median ± standard deviation), with values higher than 0.9 Pa usually recorded within densely vegetated patches (Amos et al., 2004). In the present analysis, we could not take into account the role of the biotic component because of the impossibility to reconstruct the spatial distribution of vegetated tidal flats in the oldest configurations of the Venice Lagoon. For all the above reasons, and following the approach suggested by D'Alpaos et al. (2013), we set the critical shear stress, τ_c, equal to 0.4 Pa for all the selected historical configurations of the Venice Lagoon.

Before performing the POT analysis, the time series of BSSs were processed by applying a moving average filter in order to remove spurious up-crossings and down-crossings of the prescribed threshold. This low-pass filter with a time window of 6 h removes short-term fluctuations, preserving the modulation given by the semidiurnal tidal oscillation. Thanks to this preprocessing procedure, over-threshold events satisfy the independence assumption required by the statistical analysis applied.

The POT method allowed us to identify:

1. the interarrival time of over-threshold events, defined as the time between two consecutive up-crossings of the threshold;

2. the duration of over-threshold events (that is, the time elapsed between any up-crossing and the subsequent down-crossing of the threshold);

3. its intensity, calculated as the largest exceedance of the threshold in the time elapsed between an up-crossing and the following down-crossing.

These three random variables synthetically characterize the over-threshold erosion events and can be combined to obtain further metrics to describe the erosion process (e.g. see the erosion work defined later on).

Once the probability density functions and the corresponding moments of these variables were defined, a statistical analysis was performed for each location in all the considered configurations of the Venice Lagoon in order to provide an accurate description of the BSS evolution through the last 4 centuries. This enabled us to highlight the feedback between morphology and resuspension events over long-term timescales.

We performed the non-parametric Kolmogorov–Smirnov (KS) goodness-of-fit test to verify the hypothesis that the interarrival time of over-threshold events is an exponentially distributed random variable. The interarrival probability distribution plays an important role because if interarrival times between subsequent exceedances of the threshold τ_c are independent and exponentially distributed random variables, then the mechanics of erosion events can be mathematically described as a 1-D marked Poisson process characterized by a vector of random marks (intensity and duration of each over-threshold event) associated with a sequence of random events along the time axis (Cramér and Leadbetter, 1967; Gallager, 2013). Memorylessness is one of the most interesting mathematical features of Poisson processes, since it allows one to set the probability of observing a certain number of events in a pre-established time interval, dependent only on its duration, regardless of its position along the time series (Gallager, 2013). Therefore, the description of over-threshold BSS events as a Poisson process will allow one to immediately identify the probabilities of observing a certain number of resuspension events in a year or during a season because all the sources of stochasticity in the physical drivers are described by a single parameter (i.e. the mean frequency of the process). This suggests the possibility of setting up a synthetic theoretical framework to model the wave-induced events through the use of Monte Carlo realizations, bearing important consequences for the long-term evolution of tidal landscapes.

The result of modelling erosion events as a Poisson process stands, regardless of the specific value of the censoring threshold selected for the POT analysis, provided that it is high enough to exclude deterministic exceedances, and this is confirmed also by the sensitivity analysis performed by D'Alpaos et al. (2013) on the present-day configuration of the Venice Lagoon. Indeed, when considering values of the threshold that are too low (e.g. τ_c=0.2 Pa), deterministic exceedances driven by tidal currents occur and make the interarrival time not exponentially distributed. On the contrary, as the threshold value increases (e.g. τ_c≥0.6 Pa), the KS test is still verified, thus confirming that the process remains Poisson for increasing censoring thresholds (see Fig. 6 in D'Alpaos et al., 2013, for further details).

3 Results and discussion

We analysed the time series of computed total BSSs, τ_wc, at any element of the computational grids reproducing the six selected configurations of the Venice Lagoon on the basis of the POT method in order to characterize the over-threshold erosion events in terms of interarrival times, peak excess, and duration. The KS test is then performed in each element of the six domains in order to test where interarrival times can be described by an exponential distribution, and thus, the over-threshold erosion events can be modelled as a Poisson process. We also performed the KS test on peak excess and duration to test if these marks of the process can also be described by exponential distributions. Whether peak excess and duration can be described by exponential distributions does not affect the chance to model erosion as a Poisson process, which indeed relies only on the exponentiality of interarrival times, but it can simplify the set-up of the final stochastic framework. Therefore, in the spatial distribution of the KS test results (Fig. 3), we distinguished:

1. the dark blue area, where the KS test is not verified for the interarrival time (i.e. wave-induced erosion events cannot be described as a Poisson process);

2. the red area, where the KS test is verified for all the three considered stochastic variables, namely interarrival times, intensity, and duration (i.e. wave-induced erosion events are indeed a marked Poisson process where its marks (intensity and duration) are also exponentially distributed random variables);

3. the yellow area, where the KS test is verified for the interarrival time, but it is not verified for the intensity and/or duration (i.e. wave-induced erosion events are a marked Poisson process but at least one between intensity and duration is not an exponentially distributed random variable).

Figure 3 Kolmogorov–Smirnov test for over-threshold erosion events. Spatial distribution of Kolmogorov–Smirnov (KS) test at significance level α=0.05 for the six different configurations of the Venice Lagoon: (a) 1611, (b) 1810, (c) 1901, (d) 1932, (e) 1970, and (f) 2012. We can distinguish areas where the KS test is not verified (dark blue); is verified for all the considered stochastic variables (interarrival time, intensity over the threshold, and duration) (red); and is verified for the interarrival time and not for intensity and/or duration (yellow).

The mean interarrival times (Fig. 4), mean peak excesses (Fig. 5), and mean durations of over-threshold erosion events (Fig. 6) in the six selected configurations of the Venice Lagoon are shown in every location where the KS test is satisfied for interarrival times (Fig. 3), and thus, erosion events can be described as a Poisson process. Mean peak excess and mean duration are shown also where at least the interarrival times are exponentially distributed (i.e. yellow areas in Fig. 3) because mean values are considered to be informative anyway, and erosion events can still be modelled as a Poisson process, although the marks are described by a distribution that is different from the exponential one.

Figure 4 Mean interarrival time of over-threshold erosion events. Spatial distribution of mean interarrival times of over-threshold exceedances at sites where bed shear stress can be modelled as a marked Poisson process, as confirmed by the KS test (α=0.05) for the six different configurations of the Venice Lagoon: (a) 1611, (b) 1810, (c) 1901, (d) 1932, (e) 1970, and (f) 2012.

Figure 5 Mean intensity of over-threshold erosion events. Spatial distribution of mean intensity of peak excesses of over-threshold exceedances at sites where bed shear stress can be modelled as a marked Poisson process, as confirmed by the KS test (α=0.05) for the six different configurations of the Venice Lagoon: (a) 1611, (b) 1810, (c) 1901, (d) 1932, (e) 1970, and (f) 2012.

Figure 6 Mean durations of over-threshold erosion events. Spatial distribution of mean durations of over-threshold exceedances at sites where bed shear stress can be modelled as a marked Poisson process, as confirmed by the KS test (α=0.05) for the six different configurations of the Venice Lagoon: (a) 1611, (b) 1810, (c) 1901, (d) 1932, (e) 1970, and (f) 2012.

Wind wave generation is determined by energy transfer from the wind to the water surface, and thus, it clearly depends on wind characteristics, namely wind intensity and duration, as well as on fetch length and water depth (Fagherazzi et al., 2006; Fagherazzi and Wiberg, 2009). As a consequence, the spatial distribution and morphological characteristics of channels, tidal flats, and, more importantly, salt marshes and islands strongly influence the response of a shallow tidal basin to wind forcing and the resulting distribution of BSSs (Fagherazzi et al., 2006; Defina et al., 2007). Large portions in the oldest configurations of the lagoon were occupied by salt marsh areas that were continuously interrupting the fetch and thus reducing the exceedances of the critical threshold. As a result, in the four oldest configurations, the characteristics of erosion events globally display a more complex spatial pattern, which conversely tends to be more uniform in the most recent configurations, due to the reduction in salt marsh areas, the increase in fetch length, and the deepening of tidal flats.

In all the selected configurations, salt marshes and tidal channel networks mostly represent the portion of the lagoon where wave-induced erosion events cannot be modelled as a Poisson process (dark blue area in Fig. 3). Over salt marsh platforms, almost no exceedances of the prescribed threshold, τ_c, tend to occur (Fig. S7) because of the low water depth that prevents the formation of significant waves (e.g. Möller et al., 1999). Moreover, the colonization of the salt marsh surface by halophytic vegetation almost completely prevents any vertical erosion (Christiansen et al., 2000; Temmerman et al., 2005). On the contrary, exceedances of the threshold can be detected along the channel network and at the three inlets (Fig. S7), but these are mostly associated with shear stresses produced by tidal currents, especially after the construction of the jetties at the inlets. Consequently, at these points, the KS test is not satisfied, and erosion events cannot be modelled as a Poisson process because of the strictly deterministic nature of tide-induced shear stress.

Particularly interesting is the temporal evolution of BSS at the inlets. For instance, the time series and the interarrival time probability distribution of over-threshold erosion events at the SM station in the Malamocco inlet (see Fig. 1a for the location) clearly show how the morphological modifications affected the BSS (Fig. 7). In the 1611 and 1810 simulations, in the absence of jetties at the inlets, the BSS was very small, so the number of exceedances of the threshold was too low to be representative (Fig. 7a and b). After the construction of the jetties at the Malamocco inlet in 1872, erosion mechanics abruptly changed; BSS considerably increased, but it was driven by tidal forcing, and thus, interarrival times were not exponentially distributed, since the erosion threshold was exceeded on average once per day because of tidal fluxes (Fig. 7c–f). The BSS analysis at the SL station in the Lido inlet, where the construction of the jetties ended in 1892, provides analogous results (Fig. S8). Instead, at the SC station in the Chioggia inlet, BSS still does not systematically exceed the threshold also in the 1901 configuration, since the construction of the jetties at the Chioggia inlet took place between 1930 and 1934 (D'Alpaos, 2010b) (Fig. S9).

Figure 7 Over-threshold BSS events at the Malamocco inlet. Statistical analysis at SM station in the Malamocco inlet. (a–f) Time series of the computed BSS. (g–l) Probability distributions of the interarrival times (circles) and exponential distributions (dashed lines).

The KS test is verified over subtidal platforms and tidal flats, where current-induced BSSs are typically below the critical value, but wave-induced BSSs mainly contribute to the total BSS. Locations where interarrival time, duration, and intensity follow an exponential distribution (see red areas in Fig.3) remain the vast majority of the tidal basin in all of the configurations. As a result, a synthetic framework that models erosion as a Poisson process is deemed to be suitable for wide tidal-flat areas. More generally, the chance to model erosion as a Poisson process lies in the intrinsic nature of BSS drivers. Wherever the stochastic action of wind waves and storm surges plays a prominent role in generating BSS compared to the deterministic tidal component, erosion is likely to be properly described by a Poisson process. This is the case of shallow tidal environments, where the open-water surface allows for the generation of wind waves, such as in back-barrier lagoons. On the contrary, the chance to use the Poisson-process-based approach diminishes where tidal currents substantially modulate BSS dynamics and mask the signature of stochastic processes, such as in tidal inlets and narrow meso- or macrotidal estuaries.

In almost all configurations, large interarrival times (Fig. 4) are essentially found in sheltered areas, where only particularly intense events are able to generate BSSs large enough to exceed τ_c. A clear example is provided by the area protected by marsh platforms and by the mainland in the northeastern and in the western portion of the lagoon, sheltered from the northeasterly bora wind, which is the main morphologically significant wind in the Venice Lagoon (Fig. 1b). This pattern becomes even more evident in the configurations of 1611, 1810, and 1901, where portions of the lagoon occupied by salt marshes are wider than in the more recent configurations and display interarrival times longer than 30 d. Large interarrival times can also be observed close to the three inlets, where the water depth is such that, only during intense events, the bottom can be affected by wave oscillations and the total BSSs can exceed the threshold. Globally speaking, in the four oldest configurations, we found relatively short (about 5 d) interarrival times spread all around the lagoon, while the present configuration, characterized by a more uniform and larger water depth (in some areas deeper than 1.5 m), displays longer interarrival times, e.g. between 10 and 15 d for the tidal flats located in the central-southern portion of the lagoon (Figs. 4 and S10a). This is mainly due to the relationship existing between τ_ww and water depth that, for a prescribed wind velocity, decreases as the water depth increases (Defina et al., 2007). Indeed, in the historical configurations, large areas occupied by tidal flats are characterized by lower water depth (≤0.5 m), and as a result, τ_ww is higher also for weak wind speeds, thus increasing the number of exceedances of the threshold.

A point-by-point comparison among different configurations provides further insight into the effects of morphological changes on interarrival times (Fig. 8). On a tidal flat in the northern lagoon named Palude Maggiore (see station S1 in Fig. 1a), as in most areas of the lagoon, the mean interarrival time ${\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_{\mathrm{t}}$ between two subsequent over-threshold events increases through time (Fig. 8a). This is because this area preserved the same morphological features, i.e. relatively shallow tidal flats protected by the mainland and salt marshes, over the last four centuries. On the contrary, the tidal flat in the watershed divide area between the Chioggia and the Malamocco inlets, named Fondo dei Sette Morti (see point S2 in Fig. 1a), shows a reverse trend; the interarrival times decrease in time from 1611 to the present (i.e. wave-induced erosion events are more frequent; Fig. 8d). Although the almost constant, relatively deep bottom elevation that characterized this area through centuries (Carniello et al., 2009; D'Alpaos, 2010b) prevents the exceedance of the threshold τ_c during less intense erosion events, the generalized deepening experienced by the surrounding portion of the lagoon in the most recent configurations promotes more frequent and less intense events within this area and, therefore, a decrease in the interarrival times.

Figure 8 Over-threshold erosion events at stations S1 and S2. Statistical characterizations of over-threshold events at two stations S1 Palude Maggiore and S2 Fondo dei Sette Morti (see Fig. 1a for locations) in the six configurations of the Venice Lagoon. (a, b) Probability distributions of interarrival times, t. (c, d) Probability distributions of intensities of over-threshold exceedances, e. (e, f) Probability distributions of durations of over-threshold events, d. ${\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_{\mathrm{t}}$ denotes the mean interarrival time, ${\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_{\mathrm{e}}$ the mean peak excess intensity, and ${\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_{\mathrm{d}}$ the mean duration.

The over-threshold peak intensities generally strongly increased during the last 4 centuries (Figs. 5 and S10b). For all the selected configurations, intensities are lower in the more pristine northern part of the lagoon, which is sheltered from the dominant bora wind by the mainland and by preserved salt marsh areas, interrupting the fetch. The central and southern portions of the lagoon are characterized by much larger-intensity values, which more rapidly increased over the last few decades. In particular, in the central part of the lagoon, the mean intensities increased from around 0.13 to 0.25 Pa above the threshold, due to the flattening and deepening of this area. A rather similar situation characterizes also the southern part of the Venice Lagoon between the Malamocco and Chioggia inlets.

For all the configurations investigated, the durations of over-threshold events (Figs. 6 and S10c), as intensities, present much lower values in the areas sheltered by salt marshes (i.e. the northern lagoon and the western portion of the southern lagoon) than in the fetch-unlimited central-southern portion of the lagoon. In the latter area, indeed, over-threshold events last more than 15 h, compared to a duration of about 5 h in the more sheltered areas. The increase in the peak intensities and durations of erosion events over time is also clearly shown by the probability distributions computed at points S1 and S2 (Fig. 8).

The larger over-threshold peak intensities, as well as the longer durations characterizing the central-southern portion of the lagoon and increasing from the past to the present, are in agreement with recent observations highlighting a critical erosive trend for the tidal flats and subtidal platforms in this area (Carniello et al., 2009; Molinaroli et al., 2009; D'Alpaos, 2010b; Defendi et al., 2010; Sarretta et al., 2010).

To investigate the relationship among the three random variables, the temporal cross-correlation is computed for each location and for all six configurations (Figs. S11–S13). In particular, the temporal cross-correlation between the intensity of peak excesses and the duration of over-threshold exceedances displays values very close to 1 for all the lagoon morphologies, thus suggesting a pseudo-deterministic link between peak intensities and the corresponding durations (Figs. S11 and S14a). On the contrary, almost no correlation exists between durations and interarrival times (Figs. 12 and S14b); the same applies to the correlation between intensities and interarrival times (Figs. S13 and S14c). These results, in line with the temporal cross-correlation obtained for the statistical analysis of suspended sediment concentration for the present lagoon by Carniello et al. (2016), suggest that resuspension events can be modelled as a 3-D Poisson process in which the marks (duration and intensity) are mutually dependent but independent from the interarrival time between two subsequent over-threshold events.

In order to provide a more quantitative estimation of the spatial heterogeneity of interarrival times, duration, and intensities of the critical BSS exceedances, we computed the “erosion work” (geomorphic work sensu Wolman and Miller, 1960; see also Mariotti and Fagherazzi, 2013), which represents the total amount of sediment resuspended during a selected time interval and, thus, the potential erosion because it does not consider any possible subsequent deposition. The erosion work [$E_{\mathrm{w}}^{*}$] experienced by a single point during the time interval (t₂−t₁) can be computed as

$$\begin{array}{}\text{(7)}& \left[E_{\mathrm{w}}^{*}\right]=\underset{t_{\mathrm{1}}}{\overset{t_{\mathrm{2}}}{\int }}{\displaystyle \frac{e}{{\mathit{\rho }}_{\mathrm{b}}}}\left({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{wc}}-{\mathit{\tau }}_{\mathrm{c}}}{{\mathit{\tau }}_{\mathrm{c}}}}\right)\mathrm{d}t,\end{array}$$

where e is the value of the erosion coefficient, which depends on the sediment properties, and ${\mathit{\rho }}_{\mathrm{b}}={\mathit{\rho }}_{\mathrm{s}}(\mathrm{1}-n)$ is the sediment bulk density, with the porosity n. We set ρ_s=2650 kg m⁻³ and n=0.4, in agreement with Carniello et al. (2012). In principle, when computing the actual erosion, thus, taking into account both erosive and depositional processes, the parameter e could be calibrated by comparing the modelled erosion with that retrieved from the comparison among subsequent surveys, provided that other non-natural processes (e.g. boat wave and dredging) do not strongly affect the local morphological evolution. However, because here we are considering the total potential rather than the actual erosion, such a calibration would not be correct. For this reason, we set e equal to $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ kg m⁻² s⁻¹, as usually suggested for sand–mud mixtures (van Ledden et al., 2004; Le Hir et al., 2007).

Figure 9 Erosion work. Spatial distribution of erosion work for the six different configurations of the Venice Lagoon: (a) 1611, (b) 1810, (c) 1901, (d) 1932, (e) 1970, and (f) 2012. Black identifies sites where the bottom shear stress cannot be modelled as a marked Poisson process (i.e. the KS test is not verified for the interarrival time).

Using the mean values of the stochastic variables considered herein (i.e. interarrival time, intensity, and duration), they can be modelled as a Poisson process once verified, and we can simplify Eq. (7) as follows:

$$\begin{array}{}\text{(8)}& \left[E_{\mathrm{w}}\right]={\displaystyle \frac{e}{{\mathit{\rho }}_{\mathrm{b}}}}\left({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{wc}}-{\mathit{\tau }}_{\mathrm{c}}}{{\mathit{\tau }}_{\mathrm{c}}}}\right)\left(t_{\mathrm{2}}-t_{\mathrm{1}}\right),\end{array}$$

where we assume (t₂−t₁) to be the mean duration of over-threshold events and (τ_wc−τ_c) their mean intensity. In order to estimate the erosion work for 1 year, E_w, we multiplied the result obtained with the Eq. (8) for the number of events, computed as 365 (days per year) divided by the mean interarrival time at each point within the lagoon. Instead, using the complete formulation in Eq. (7), the erosion work over 1 year, $E_{\mathrm{w}}^{*}$, can be simply computed simply by extending the integration period to the entire year.

We computed the spatial distribution of the annual erosion work E_w using the synthetic approach for the six configurations of the Venice Lagoon (Fig. 9). We also computed the erosion work according to Eq. 7, in order to compare differences between the complete formulation based on the computed BSS time records and the synthetic approach exploiting the possibility of describing resuspension events as marked Poisson processes (Fig. S15). The erosion work computed following the two approaches is quite similar, as shown by the map of the relative error (Fig. S16) and by the computed values of the spatially averaged relative error, which varies between 10 % and 14 %, considering all the analysed configurations of the lagoon (Table 1). Such an agreement between the two estimates of the erosion work supports the validity of the provided statistical characterization of resuspension events.

Table 1 Spatially averaged relative error between erosion work computed with Eqs. (7) and (8).

As it is defined, the erosion work represents the total potential erosion, rather than the effective erosion, because it neglects the subsequent possible deposition. However, a comparison between the erosion work and the actual erosion, which can be retrieved from the comparison between surveys, can still provide some interesting insights and highlight erosive trends.

Focusing in particular on the central-southern lagoon, the erosion work is almost constant between 1611 and 1932; it reaches its maximum in 1970, and then it decreases in the present configuration (Fig. 9). The four oldest configurations (i.e. 1611, 1810, 1901, and 1932) display a more complex spatial pattern of the computed erosion work because of the wider presence of salt marshes and islands distributed throughout the basin and because of the shallower and more irregular bathymetry characterizing the tidal flats. This morphology is such that the fetch is continuously interrupted, and wind waves are prevented from fully developing, while generating and propagating over areas whose bathymetry is continuously varying. Interestingly, even if the present configuration of the lagoon displays larger mean intensities and longer mean durations than in 1970 (see Figs. 5 and 6), the combination with generally longer mean interarrival times (Fig. 4) affects the erosion work. Indeed, the erosion work is at its maximum in the 1970 configuration, when it reaches a peak of more than 4.0 cm yr⁻¹. This promoted an intense and uniform erosion of the lagoon, thus leading to the present morphology and bathymetry characterized by less complex erosion patterns and a roughly constant erosion work on the tidal flats in the central-southern lagoon of about 2.5 cm yr⁻¹. Our results quantitatively support previous studies (Carniello et al., 2009; Molinaroli et al., 2009) that identified two different evolutionary trends in the northern lagoon and in the central-southern part, with the northern lagoon displaying, on average, much lower erosion rates.

4 Conclusions

Our results provide a statistical characterization of sediment erosion dynamics, aimed at testing the possibility for describing erosion events as a Poisson process in a synthetic modelling framework able to reproduce the long-term evolution of shallow tidal systems. The proposed approach aims to better describe erosion events in shallow tidal environments, where BSS dynamics are strongly affected by wind conditions.

In the present study, the approach is applied to the specific case of the Venice Lagoon for which six morphological configurations along the last 4 centuries are available. We applied the extensively calibrated and tested WWTM to the six historical configurations of the Venice Lagoon, in order to perform a spatially explicit analysis of the BSS time series under the same wind- and water level forcing. We analysed the computed BSS temporal evolution following the peak-over-threshold theory. We verified whether wind wave erosion events could be modelled as a marked Poisson process by performing the non-parametric Kolmogorov–Smirnov goodness-of-fit test to confirm the hypothesis that the interarrival time of over-threshold BSS events together with their durations and intensities are exponentially distributed random variables.

Statistical analyses of the wave-driven erosion processes suggest that interarrival times between two consecutive over-threshold events and their durations and intensities can be described as exponentially distributed random variables over wide areas in all the selected configurations of the Venice Lagoon. As a consequence, the wave-induced erosion can be represented by a marked Poisson process through centuries.

Furthermore, we observed that the durations and intensities of over-threshold BSS exceedances are highly correlated, while almost no correlation exists between the duration and interarrival time, as well as between intensity and interarrival time. These observations indicate that a 3-D Poisson process, in which the marks (duration and intensity of the over-threshold events) are mutually dependent but independent from the interarrival time, provides a suitable description of the wave-induced erosion processes.

Moreover, we showed that in the last 4 centuries, the interarrival times of erosion events generally increased everywhere within the lagoon, as well as their intensities and durations, thus leading to less frequent but more intense wave-induced erosion events. These modifications in the bottom shear stress field are generated by, but at the same time they are also responsible for, the morphological modifications of the Venice Lagoon, in particular the generalized deepening of tidal flats and reduction in the salt marsh area. Only in the Fondo dei Sette Morti, located close to the watershed divide between the Malamocco and the Chioggia inlets, did interarrival times decrease in the last 4 centuries. Such an opposite trend is associated with the relatively deep and constant bottom elevation characterizing this area, combined with the generalized deepening experienced by the surrounding areas that allows more frequent events reaching the Fondo dei Sette Morti.

The erosion work, computed as a combination of interarrival times, durations, and intensities, remained almost constant and characterized by an irregular spatial pattern until the beginning of the 20th century, when it rapidly increased and reached a peak in 1970. In the last few decades, the erosion work decreased, presenting a more uniform pattern and suggesting that the quite intense erosive trend the Venice Lagoon has been experiencing since the beginning of the last century is, at present, slowing down as a consequence of the generalized deepening and flattening of the lagoonal bed. Owing to the choice of forcing the domain with the same conditions, these changes in the erosive trend are, in fact, only due to morphological modifications experienced by the tidal basin.

The present findings, together with the statistical characterization of suspended sediment dynamics (Tognin et al., 2024), represent a step towards the set up of a synthetic statistically based framework which can be used to model the long-term morphodynamic evolution of shallow tidal systems through the use of independent Monte Carlo realizations, thus possibly exploring a large set of equally likely lagoonal configurations.