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

,

1 The Wind-Wave-Tidal Model

Hydrodynamic model
The hydrodynamic module solves the 2D depth-integrated shallow water equations, phase averaged over a representative elementary area in order to deal with wetting and drying processes in very shallow and irregular domains (Defina, 2000): where t is time, D/Dt is the material time derivative, ∇ and ∇• denote the 2D gradient and divergence operators, respectively.
In the continuity equation (Eq.1), η is the free surface elevation, ϑ is the wet fraction of the computational domain, which is a function of water depth and local topographic unevenness (Defina, 2000), q = (q x , q y ) is the depth-integrated velocity (i.e., discharge per unit width).In the momentum equations (Eq.2), Y is the effective water depth (i.e., the water volume per unit area), τ t and τ s are the shear stresses at the bottom and at the free surface, respectively, ρ is water and g is gravity acceleration.
The Reynolds stresses Re are computed using a depth-averaged version of Smagorinsky's model (Smagorinsky, 1963) and they read: with i, j denoting either the x or y coordinate, u = q/Y and the eddy viscosity ν e is computed as S A e 2(u x,x ) 2 + (u x,y + u y,x ) 2 + 2(u y,y ) 2 (4) where the Smagorinsky coefficient C S is set equal to 0.2 and A e is the area of the computational element.
1 Supplement In the numerical scheme, the material time derivative in Eq. 2 is solved with the method of characteristics by expressing it as the finite difference in time.This allows solving the continuity equation (Eq. 1) with a semi-implicit scheme, resulting in a self-adjoint spatial operator, which is solved on a staggered triangular grid with the finite element method of Galerkin and flow rates are obtained by back substitution.

Wind wave model
The wind-wave module (Carniello et al., 2011) solves the wave action conservation equation using the same computational grid of the hydrodynamic module, which provides water depths and depth-averaged flow velocities, used to propagate the wind-wave field.In the frequency domain, the wave action density, N 0 , evolves according to (Carniello et al., 2005) where c ′ gx and c ′ gy are the wave group celerity components used to approximate the speed propagation of N 0 (Holthuijsen et al., 1989), S 0 denotes the wind-wave source terms, accounting both for positive (wind energy input) and negative (bottom friction, whitecapping, and depth-induced breaking) contributions.
The wave-action conservation equation (Eq.5) is solved with an upwind finite volume scheme based on the same computational grid of the hydrodynamic model.In each element at each time step, the wind-wave model computes the wave action, from which the significant wave height is obtained by applying the linear theory.

Model calibration
The calibration of the model was performed for different periods, but here we focus as an example on the period 8-13 November 2004, for which not only water levels, but also flow rates at the inlets are available.This period was initially characterized by quite low wind speeds (U wind =7.0 m/s between 8 and 10 November) followed by an intense Bora event with wind speeds up to 18 m/s (Figure S2).Model capability to capture the process is evaluated by means of the Nash-Sutcliffe Model Efficiency (NSE) (Allen et al., 2007).Figure S2 shows the comparison between measured and computed water levels at three tide-gauge stations: Pagliaga, Punta della Salute and Brondolo, located in the northern, central and southern portions of the lagoon, respectively.The comparison highlights a very good agreement between computed and measured water levels (NSE mean = 0.970, NSE median = 0.984, NSE std = 0.040) Figure S3 compares computed and measured flow rates at the three inlets in the same period.The comparison is quite favourable (NSE mean = 0.853, NSE median = 0.931, NSE std = 0.184) and only when the wind speed reaches the highest value of about 18 m/s the model slightly underestimates the maximum water discharge.The windwave measurements within the Venice lagoon were performed only in 2002 and 2003, therefore, we focus as an example on the period 11-14 February 2003. Figures S4 and S5 show the comparison between measured and computed wave peak periods and wave heights, respectively, at the 1BF station, located in the northern lagoon, and the 2BF station located in the central lagoon (see Carniello et al. (2011) for the location of the stations).The comparison between measured and computed wave heights and periods confirms the capability of the coupled WWTM to reproduce the modulation of the wave height induced by water-level oscillations and wind-speed variations (NSE mean = 0.627, NSE median = 0.756, NSE std = 0.357).In particular, Figure S4 shows and water depth (Carniello et al., 2011).

Boundary conditions selection
To simulate the typical wave-driven BSS conditions in the Venice Lagoon, we carefully analyse the wind climate measured in the period 2000-2020.In particular, we compared the wind velocity measured in each year with that of the whole period 2000-2020 by means of the two-sample Kolmogorov-Smirnov (KS) test and the Wilcoxon (W) test (Table S6).We selected the 55 wind velocity measured in 2005 because it does not significantly differ from that of the whole considered period according to both the KS and W tests.A visual comparison of the cumulative distribution frequency F (x) of wind velocity also confirms this choice (Figure S6).

Figure S1 .
Figure S1.Meshes used in the numerical model.Meshes of the six different configurations of the Venice Lagoon: (a) 1611, (b) 1810, (c) 1901, (d) 1932, (e) 1970, and (f) 2012.For the sake of clarity, the mesh portion representing the sea is shown only in panel (e) and (f), but it is present in all meshes.

Figure S2 .
Figure S2.Model calibration: water level.Comparison between observations and model results for the period 8-13 November 2004: (a) Measured wind speed and direction; comparison of measured (circles) and computed (solid lines) water levels at the Pagliaga (b); at the Punta della Salute (c); and at the Brondolo station (d) (adapted from Carniello et al., 2011).

Figure S3 .
Figure S3.Model calibration: flow rate.Comparison between observations and model results for the period 8-13 November 2004: (a) Measured wind speed and direction; comparison of measured (circles) and computed (solid lines) water discharges at the Lido Inlet (b); at the Malamocco Inlet (c); and at the Chioggia Inlet (d).Positive values refer to the food phase, negative values to the ebb phase (adapted from Carniello et al., 2011).

Figure S7 .Figure
Figure S7.Number of upcrossings of the erosion threshold.Spatial distribution of the number of upcrossings of the threshold for erosion τc = 0.4 Pa for the six different configurations of the Venice Lagoon: (a) 1611, (b) 1810, (c) 1901, (d) 1932, (e) 1970, and (f) 2012.

Figure
Figure S9.Over-threshold BSS events at the Chioggia inlet.Statistical analysis at SC station in the Chioggia inlet: time series of the computed BSS (a-f); probability distributions of the interarrival times (circles) and exponential distributions (dashed lines) (g-l).

Figure S10 .
Figure S10.Spatial probability density function of interarrival time, intensity and duration of BSS over-threshold events.Probability density function (left), mean (mean ± standard deviation) and median value (right) of interarrival times t (a), intensity e (b) and duration d (c) of BSS over-threshold events.

Figure S14 .
Figure S14.Spatial probability density function of cross-correlation between interarrival time, intensity and duration of BSS overthreshold events.Probability density function (left) and mean value (mean ± standard deviation, right) of cross-correlation between intesity and duration ρ(e − d) (a), interarrival time and duration ρ(t − d) (b) and interarrival time and intensity ρ(t − e) (c).

Table S1 .
Results of the Kolmogorv-Smirnov (KS) and Wilcoxon (W) tests on wind velocity.