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).