Figure 2: Morphological evolution of the Venice Lagoon. Bathymetry of the Venice Lagoon in 1887 (A), 1901 (B), 1932 (C), 1970 (D), 2003 (E), and 2014 (F), as reconstructed from available historical topographic and bathymetric data. G,H,I,J) Morphology of the Venice Lagoon according to the hypothetical scenarios of marsh erosion analyzed in the present study. These scenarios assume different rates of marsh erosion equal to 25% (2014-E25), 50% (2014-E50), 75% (2014-E75), and 100% (2014-E100), respectively. Temporal variation of salt-marsh extent (K) and spatially-averaged bottom elevation of tidal flats (L) between 1887 and 2014 are also shown.
Further manmade modifications of the inlet morphologies were carried out between 2006 and 2014 to accommodate the mobile floodgates of the Mo.S.E. (“Modulo Sperimentale Elettromeccanico”) system (Figure 1 B,C,D), designed to protect the city of Venice and other lagoon settlements from extensive floodings (Mel, Viero, et al., 2021). These interventions slightly increased hydraulic resistances and led to both a reduction of tidal amplitude and an increase in tidal-phase delays within the lagoon (Ghezzo et al., 2010; Matticchio et al., 2017), thus partially mitigating the dominance of ebb-tidal flows (e.g., Finotello et al. 2019). In spite of this, however, salt-marsh erosion within the lagoon is still ongoing, though at reduced rates compared to the last century (Figure 2 K) (Finotello et al., 2020; Tommasini et al., 2019). Moreover, operations of the Mo.S.E. floodgates will further reduce the resilience of salt marshes to rising relative sea levels by preventing inorganic deposition during storm-surge events which, though episodic, critically contribute to marsh accretion in sediment-starving systems such as the Venice Lagoon (Tognin et al., 2021, 2022).
3 Methods

Numerical model

We employed the bidimensional, depth-averaged, finite element numerical model developed by Carniello et al. (2005, 2011), which is suitable to reproduce the hydrodynamics of shallow tidal basins driven by tidal flows and wind fields. In the following, we report a brief description of the model and refer the reader to Carniello et al. (2011) for further details. The model consists of two coupled modules, namely a hydrodynamic and a wind-wave module, and will be referred to as WWTM (Wind-Wave Tidal Model) hereinafter.
The hydrodynamic model solves the bidimensional, depth-averaged shallow water equations, suitably modified to reproduce wetting and drying processes in very shallow and irregular domains, which read (after Defina 2000):
\begin{equation} \frac{\partial q_{x}}{\partial t}+\frac{\partial}{\partial x}\left(\frac{q_{x}^{2}}{Y}\right)+\frac{\partial}{\partial y}\left(\frac{q_{x}q_{y}}{Y}\right)-\left(\frac{\partial R_{\text{xx}}}{\partial x}+\frac{\partial R_{\text{xy}}}{\partial y}\right)+\frac{\tau_{\text{bx}}}{\rho}-\frac{\tau_{\text{wx}}}{\rho}+gY\frac{\partial H}{\partial x}=0\nonumber \\ \end{equation}\begin{equation} \frac{\partial q_{y}}{\partial t}+\frac{\partial}{\partial x}\left(\frac{q_{x}q_{y}}{Y}\right)+\frac{\partial}{\partial y}\left(\frac{q_{y}^{2}}{Y}\right)-\left(\frac{\partial R_{\text{xy}}}{\partial x}+\frac{\partial R_{\text{yy}}}{\partial y}\right)+\frac{\tau_{\text{by}}}{\rho}-\frac{\tau_{\text{wy}}}{\rho}+gY\frac{\partial H}{\partial y}=0\nonumber \\ \end{equation}\begin{equation} \eta\frac{\partial H}{\partial t}+\frac{\partial q_{x}}{\partial x}+\frac{\partial q_{y}}{\partial y}=0\nonumber \\ \end{equation}
where \(t\) is time, the \(x\) and \(y\) subscripts denote the directions of a given variable in a Cartesian reference system, \(q\) is the flow rate per unit width, \(R\) represents the depth-averaged Reynolds stresses, \(\tau_{b}\) is the bottom shear stress produced by tidal currents, \(\tau_{w}\) is the wind-induced shear stress at the free surface whose elevation is \(H\), \(\rho\) stands for the fluid density, \(g\) is the gravitational acceleration, \(Y\) is the water volume per unit area ponding the bottom (i.e., the equivalent water depth) and \(\eta\) is the wet fraction of the computational domain which accounts for surface irregularities during the wetting and drying processes (see Defina, 2000, for a detailed description of the hydrodynamic equations). A semi-implicit staggered finite element method based on discontinuous Galerkin’s approach is used to solve the governing equations (Defina, 2003). In order to solve the closure of the Reynolds stresses which appear in the momentum equations in the two horizontal directions, a suitable eddy viscosity is as customary introduced and evaluated by means of the Smagorinsky (1963)’s model. The horizontal components of the Reynolds stresses then read
\begin{equation} \begin{matrix}R_{\text{xx}}=2\nu_{e}\left(\frac{\partial q_{x}}{\partial x}\right)\\ R_{\text{xy}}=\nu_{e}\left(\frac{\partial q_{x}}{\partial y}+\frac{\partial q_{y}}{\partial x}\right)\\ \end{matrix}\nonumber \\ \end{equation}
The eddy viscosity \(\nu_{e}\) differs from the standard one as it also encloses the contribution from the stresses produced by the subgrid momentum exchange, which is yet difficult to evaluate owing to its dependence on the full three-dimensional morphology of the bottom surface (see Defina, 2000 and D’Alpaos and Defina, 2007 for further details). The hydrodynamic module provides the wind-wave module with water levels and depth-averaged velocities that are employed for calculating wave group celerity and bottom shear stresses (induced by both wind waves and tidal currents), as well as for evaluating the influence of flow depth on wind-wave propagation.
The wind-wave module employs the same computational grid of the hydrodynamic model to solve the wave action conservation equation (Holthuijsen et al., 1989). The latter is simplified by assuming that the direction of wave propagation instantaneously readjusts to match the wind direction (i.e., neglecting refraction). The module describes the evolution of the wave action density (\(N_{0}\)) in the frequency domain and it reads (Carniello et al., 2011):
\begin{equation} \frac{\partial N_{0}}{\partial t}+\frac{\partial}{\partial x}c_{\text{gx}}^{{}^{\prime}}N_{0}+\frac{\partial}{\partial y}c_{\text{gy}}^{{}^{\prime}}N_{0}=S_{0}\nonumber \\ \end{equation}
where \(c_{\text{gx}}^{{}^{\prime}}\) and \(c_{\text{gy}}^{{}^{\prime}}\) represent the wave group celerity in the \(x\) and \(y\) direction, respectively, and are used to approximate the propagation speed of \(N_{0}\) (Carniello et al., 2005; Holthuijsen et al., 1989), while \(S_{0}\) represents all the source terms describing the external phenomena contributing to wave energy variations, which can be either positive (wind energy input) or negative (bottom friction, white capping, and depth-induced breaking). Based on the relationship between peak-wave period and local wind speed and water depth (Young & Verhagen, 1996), the model can compute both the spatial and temporal distribution of the wave period. The linear wave theory also allows one to relate the local significant wave height to the horizontal orbital velocity at the bottom and, therefore, to the wind-wave induced bottom shear stress (\(\tau_{\text{ww}}\)). The nonlinear interactions between the latter and the current-induced bottom shear stress \(\tau_{b}\) are accounted for through the empirical Soulsby’s (1995) formulation, which enhances the value of the bottom shear stress \(\tau_{\text{wc}}\) beyond the mere sum of \(\tau_{b}\)and \(\tau_{\text{ww}}\) (see detailed descriptions in Carniello et al. (2005), their equations (26) and (27)). The WWTM has been extensively tested and calibrated against both hydrodynamic and wind-wave field data from the Venice Lagoon (Carniello et al., 2005, 2011; L. D’Alpaos & Defina, 2007; Tognin et al., 2022; Tommasini et al., 2019).

Numerical simulations

3.2.1 Computational Grids
Numerical simulations were performed considering ten different morphological configurations of the Venice Lagoon (Figure 2 ). Six of these configurations represent actual past-lagoon morphologies, reconstructed from available topographic and bathymetric data (Figure 2 A-F), while the additional four configurations consist of hypothetical scenarios characterized by additional marsh loss relative to the present-day lagoon morphology (Figure 2 G-J). Specifically, to understand how marsh loss affected the hydrodynamics of the Venice Lagoon in the past, we utilized six already existing WWTM computational grids representing the morphological configurations of the lagoon from 1887 to 2014 (Figure 2). The 1887 and 1901 grids were constructed based on “Topographic/hydrographic map of the Venice Lagoon” produced by the Genio Civile of Venezia in 1901, and are identical to each other except for the different morphology of the Lido inlet, where only the northern jetty was present in 1887 while both the jetties were completed in 1901. In contrast, different topographic surveys carried out in 1932, 1970, and 2003 by the Venice Water Authority (Magistrato alle Acque di Venezia) were employed to create the computational grids relative to the years 1932, 1970, 2003, and 2014, respectively (Carniello et al., 2009) (Figure 2A-E). The 2014 computational grid is based on the most recent, 2003, bathymetric survey and accounts also for the anthropogenic modifications at the three inlets related to the Mo.S.E. system, which were completed in 2014 (Figure 2F). All these grids have been accurately calibrated and utilized in several previous studies (Carniello et al., 2009, 2016; L. D’Alpaos, 2010; Finotello et al., 2019, 2020; Finotello, Capperucci, et al., 2022; Silvestri et al., 2018; Tommasini et al., 2019). Because cell-bed elevations in each computational grid refer to the mean sea level recorded when each survey was performed, and numerical simulations were carried out by forcing the model with tidal waves imposed at the open sea boundary and oscillating around the mean sea level (see next paragraph), historical rises in relative sea level are implicitly accounted for.
In contrast, hydrodynamic effects of additional marsh losses were investigated based on four possible scenarios characterized by different degrees of marsh-area loss equal to 25% (E25, Figure 2G), 50% (E50, Figure 2H), 75% (E75, Figure 2I), and 100% (E100, Figure 2J) relative to the present-day marsh extent (see Donatelli et al. 2018). These scenarios were modeled utilizing the 2014 computational grid as a baseline, and gradually removing marsh areas following the approach adopted by Donatelli et al. (2018). Specifically, computational elements located along eroding marsh margins are selected and their characteristics in terms of elevation and roughness are modified to match those of the surrounding tidal flats. The major shortcoming of this approach is that the hydrodynamic effects of marsh erosion may be largely dependent on the location of the eroded marshes. To overcome this limitation and provide plausible results, marsh elements are removed preferentially where larger wind-induced erosion rates are observed based on combined remote sensing analyses and numerical modeling of wave power (see Tommasini et al. 2019; Finotello et al. 2020). The mean sea level was held constant in all the simulations involving additional loss of salt marsh area. It is worthwhile noting that this approach implicitly assumes that sediments eroded from salt marshes are almost instantaneously removed from the basin, so that they can no longer contribute to active sedimentation on either tidal flats or salt marshes and potentially mitigate the negative effects of marsh loss (Elsey-Quirk et al., 2019; Mariotti & Carr, 2014). Hence, even though our analyses can hardly be used to draw inferences on the long-term effects of marsh loss on the lagoon net sediment budget, they can still be considered representative of possible scenarios of marsh loss in the Venice Lagoon, especially in view of the reduced sediment input to the marshes expected after the activation of the Mo.S.E. barrier system (Tognin et al., 2021).
3.2.2 Boundary Conditions
In the numerical model, water levels are imposed at the seaward boundary of the computational domain, the latter including also a portion of the northern Adriatic Sea in front of the Venice Lagoon (see Figure 3 B). Water level data are measured at the CNR Oceanographic Platform, which is located in the Adriatic Sea, about 15 km away from the coastline. In contrast, data concerning wind speeds and directions are measured at the “Chioggia Diga Sud” anemometric station (Figure 3 B) and are applied to the whole lagoonal basin, as they nicely correlate with wind measurements in other locations within the lagoon (see Carniello et al., 2005 for details).