The Holocene influence on the future evolution of the Venice Lagoon tidal marshes

The resilience of marsh ecosystems to expected sea-level rise is determined by a complex interplay of organic and inorganic sedimentation dynamics. Marshes have formed over past centuries to millennia and consist of extremely reactive bodies with sediments that can experience high compaction. Here we provide a quantification of the degree to which the past history of a salt marsh can affect its long-term evolution. A dataset of elevation dynamics was established in the Venice Lagoon (Italy) and interpreted using a physics-based model of deposition and large consolidation of newly deposited material. We found that the fate of low-lying tidal landscapes over the next century of accelerating sea-level rise will be highly dependent on compaction of soft, recently deposited soils. Our results imply that a sedimentation rate twice the present rate will be needed to counterbalance the expected sea-level rise. Holocene compaction in the Venice Lagoon, Italy, highlights the importance of soil properties and deposition rate in predicting the evolution of tidal marshes, according to numerical model simulations calibrated with 20 years of observations.

T idal marshes are one of the most dynamic ecosystems. They have formed over the last millennia as the superposition of interrelated biotic and abiotic processes 1 . The interactions and feedbacks among these mechanisms exert a significant influence on the evolution of marsh elevations in relation to mean sea level (MSL), which is crucial to maintain the delicate equilibrium that these complex systems provide to wildlife and plants [2][3][4] . The increasing marsh loss experienced in many coastal areas worldwide 5,6 indicates that an investigation of these aspects with a more complete perspective is needed for a better understanding of the ongoing processes. This new knowledge is especially needed for conservation and restoration policy planning to preserve coastal wetlands, which are recognized as crucial ecosystems that purify water; provide protection from floods, droughts and other disasters; provide food and livelihoods to millions of people; support rich biodiversity; and store more carbon than any other ecosystem 7 . The notion that salt marshes are carbon-sink ecosystems is widely recognized worldwide (e.g. 8,9 ) and great efforts are currently put into finding blue carbon offset methodologies 10 . However, wetlands are globally threatened by various factors such sea-level rise (SLR), local subsidence (shallow compaction and deep subsidence), decrease in sediments availability, edge erosion by ship waves, and other anthropogenic activities (e.g., uncontrolled tourism), interacting at various spatio-temporal scales [11][12][13][14][15][16][17] . In particular, a reduction in land elevation may decrease carbon-storage capacity for future CO 2 sequestration, becoming a potential source of carbon previously buried in the sediments.
In this contribution we study the coupling dynamics between surface and subsurface processes occurring over a millennial time scale in the San Felice salt marsh located in the Venice Lagoon (Italy) (Fig. 1) that is representative of a micro tidal marsh environment (tidal range~1.0 m). The main focus of this research is to investigate how the Holocene strata formed in the last few millennia consolidates and influences the current surface elevation changes. This information is fundamental to providing reliable predictions of salt marsh resilience under mutable forcing factors, such as those related to climate change.
The city of Venice and its lagoon is a critical coastal area due to the peculiar morphological setting. Indeed, the map in Fig. 1a highlights the large variability in terms of subsidence suggesting that even a few centimeters of land elevation loss relative to the Adriatic MSL may have critical effects amplified by the impacts associated with SLR. During autumn-winter and under specific meteorological conditions, the lagoon is exposed to the largest inflows of water from the sea through three inlets, causing the socalled "Aqua Alta" (high tide > 110 cm above MSL). Depending on the duration and intensity of these events, the city and the lagoon salt marshes undergo partial or total submersion in relation to their elevation above MSL. Clearly, this scenario will become worse in the case of increasing subsidence rates and SLR. Recent studies in this area show that salt marshes are characterized by a wide range of subsidence rates with displacements ranging from small uplifts to subsidence rates of more than 20 mm year −1 , with average values of 3-4, 1-2 and 2-3 mm year −1 in the northern, central and southern portions of the lagoon, respectively, as detected by SAR (Synthetic Aperture Radar) interferometry on a network of artificial radar reflectors 18,19 . The elevation gain/loss in natural landforms is mainly controlled by the availability of organic/inorganic sediments deposition, the consolidation of deposited material due to sediments' weight, the erosion caused by currents and waves, and external forcings such as increasing MSL and deep subsidence related to the movements of the Pleistocene basement [20][21][22][23][24] . Usually, the formation and evolution of salt marshes are studied with models that do not explicitly account for the consolidation of the marsh sediments. Conversely, this process is roughly included in relative SLR 25,26 or calculated through simple  18 ). Positive values indicate uplift, negative land subsidence. The enlargement in panel b shows the SEC behavior of the San Felice salt marsh whereas panel c shows the frequency distribution of the SEC rates (i.e., the percentage of point targets (PTs) in each SEC class-rate over the total amount of PTs detected on the marsh) (data from 33 ). Panels d-f show typical marsh environments with creeks, vegetated and muddy areas. empirical relationships [27][28][29] . A physics-based numerical model of natural consolidation may be critical to accurately quantifying the dynamics of the surface elevation at centennial to millennial time scales 30 , allowing for specific estimations of deposition rates required to counteract future scenarios of relative SLR given by the superposition of regional subsidence, shallow compaction and SLR. Adequate quantification of compaction has been shown as significant for relative sea level reconstructions using salt marsh sediment sequences 31 , providing valuable insights into the mechanisms affecting the resilience of marshes to SLR.
We propose a first-of-its-kind dataset interpretation with a physics-based model of soil consolidation originally developed by Zoccarato and Teatini 32 . The elevation dynamics of the San Felice (Fig. 1b, c) salt marsh is accurately described by taking advantage of high-resolution remote sensing and on site observations performed using a comprehensive monitoring network of subsidence and accretion rates at different spatio-temporal scales. The model is calibrated using data collected in the period 1993-2013 through a network of sensors and instruments deployed in San Felice salt marsh (Figs. 1d-f and 2), even though not all instruments worked simultaneously. The resilience of the San Felice salt marsh was investigated considering current and accelerating SLR rates under RCP-derived scenarios. Preliminary predictions of salt marsh losses by the end of 2100 in the Venice Lagoon are given based on sediment availability up to twice the present accretion rate. However, under SLR rates much higher than present rates, we claim the need for a comprehensive morphodynamicgeomechanical model to properly estimate the accretion rates accordingly to the MSL over the marsh surface.

Results
Measurement dataset at the San Felice salt marsh. At the San Felice salt marsh land subsidence and uplift were measured with the monitoring network as shown in Fig. 2. It consists of (i) artificial radar trihedral corner reflectors (TCRs) 19 , (ii) a group of deep benchmarks monitored by very high-precision levelling (HPL), (iii) two different surface elevation table (SET) apparatuses (RSET and SET in the following), (iv) feldspar marker horizons (FMH), and (v) point targets (PTs) consisting of natural reflectors scattered over the marsh surface detected by SAR interferometry. Figure 2 shows the reference depths and the data collection time window for each instrument as well as the images of the typical set up of the monitoring stations and their geographic locations. The main characteristics of the measurement techniques are summarized in Table 1. Each instrument/sensor monitors different variables as described in the following section (see also the sketch in Fig. 3a):  • High-precision levelling was used to measure the relative elevation between three benchmarks, whose reference depths are 8.0, 2.4 and 0.5 m. Records processing allows to quantify the deformation (DE) of the two layers between the benchmarks, that is, the relative displacements between 0.5 and 2.4 m and 8.0 and 2.4 m. DE can be either positive (expansion) or negative (compaction).
• FMH consists of feldspar minerals upon which sediments naturally settle 34 . They were placed on the marsh surface to measure its vertical accretion (AC). Accretion corresponds to the thickness of the sediment layer deposited over the feldspar layer.
• Two SET apparatuses 35 allowed to measure the relative surface elevation change (rSEC) with respect to the reference depth of the instrument. SET and RSET mainly differed in terms of the reference depths that are 0.5 m and 2.0 m from the marsh surface, respectively. SET data combined with FMH enabled the quantification of the deformation caused by shallow compaction of the surface layer ranging from the SET-rod depth to the marsh surface. DE can be easily estimated by the simple relationship DE = rSEC − AC (modified after 36 where negative DE values indicate compaction).
• SEC values were obtained by satellite measurements on point targets (PTs) consisting of natural radar reflectors scattered over the marsh and detected by a specific PSI approach, adopting a "short-term candidate pixel" criterion suitable for salt marsh surfaces 33 . A total of 151 PTs (Fig. 1b) were retrieved at the San Felice salt marsh. From a physical point of view these radar targets are portions of the marsh surface keeping a coherent response or, for example, large wood remnants deposited on the ground (Fig. 2f). Positive and negative values indicate land surface uplift and subsidence, respectively. We underline that SEC values by InSAR-derived PTs data provide a comprehensive information on all processes contributing to the displacements of the marsh surface, i.e., below-ground deformation plus accretion.
The measurement values are reported in the plots in Fig. 3b. They represent the average values over the total number of available measurements for each instrument/sensor. Absolute displacement on the order of −1 ± 0.1 mm year −1 was detected by the 4 m-depth TCR, providing comparable values with the regional subsidence of this area 37,38 . The HPL values confirmed minimal deformations below a 0.5 m-depth with +0.3 ± 0.  (Fig. 1b, c). Combining the FHM and SET measurements, the mean rate of shallow compaction was −1.2 mm year −1 (1.8 ± 3.7 mm year −1 minus 3.0 ± 1.8 mm year −1 ) for the 0.5 m topsoil and −1.0 mm year −1 (2.0 ± 2.5 mm year −1 minus 3.0 ± 1.8 mm year −1 ) for the 2 m topsoil, with the major soil consolidation occurring in the shallowest 0.5 m depth interval. Notice that HPL and SET measurements show discrepancies probably due to the high heterogeneity of these environments and the consequent monitoring complexity.
Salt marsh vertical dynamics: contribution of below-ground deformations and accretion rates. Integrated analysis of available observations at the San Felice salt marsh (Fig. 3b) enables quantification of the belowground deformations at different depths and net accretion rates of the marsh surface (Fig. 4a, b). The thickness of the Holocene sediments plays a significant role in the variability of the measured subsidence rates in the Venice Lagoon and in areas of increased sediment deposition 6,39,40 . In the study area, the Holocene deposits are~10 m thick according to the map of the Holocene deposits published by Tosi et al. 41 and obtained by combining data from stratigraphic, sedimentological, and radiochronological cores with high-resolution seismic surveys 42 .  (Table 1).
However, for the purpose of this analysis, we considered only depths in the range of 0-4 m rather than the entire Holocene thickness because we consider a negligible compaction of the layer between 4 and 8 m. Since the 4 m-depth TCR quantifies deep (regional) subsidence values, the deformation measured by the HPL between 8 and 2.4 m is entirely attributed to the layer between 4 and 2.4 m. Two profiles of the below-ground vertical dynamics are provided in Fig. 4a, b using data from the 4 m-depth TCR, the HPL, the RSET (only Fig. 4a) and the SET (only Fig.4b).
To calculate the land-surface vertical displacements (VD surface ) of the profiles in Fig. 4a, b, the following relationship was applied: Fig. 4b. VD at a certain depth below the marsh surface was simply obtained by the above relationships by adding to the 4 m-depth TCR term the deformation measured below the depth of interest. Notice that, assuming deformation uniformly distributed within each monitored depth interval, the displacements change linearly. By using the profiles of Fig. 4a, b, VD surface = −2.4 mm year −1 and VD surface = −2.3 mm year −1 , respectively. Employing data from the FMH, with accretion rates on the order of 3 ± 1.7 mm year −1 , surface-elevation gain less than 1 mm year −1 was found for this area, confirming the analysis obtained with Persistent Scatterer Interferometry on natural PTs. It is clear that shallow compaction plays a major role in the net accretion of such ecosystems that consist of relatively recent deposits undergoing substantial consolidation. To further investigate the relative importance of the driving mechanisms on marsh surface vertical movements we employed a physics-based model 32 , where the reconstruction of the transitional landform during the Holocene followed the accretion/enlargement of the sedimentary body due to deposition and consolidation of the newly deposited material. The accretion term considers material sedimentation over the body surface, whereas compaction reflects the progressive consolidation of the porous medium under the overburden load of overlying younger deposits. The modeling approach is based on a 2D groundwater flow simulator coupled to a 1D vertical geomechanical module. The vertical dynamics were studied within a temporal scale of 7000 years covering marsh formation and evolution of the upper Holocene strata 43 above the Pleistocene basement. The model also considered the geometric non-linearity arising from the consideration of large solid grain movements by using a Lagrangian approach with an adaptive mesh. The dynamic mesh allowed the grid nodes to follow the soil grains in their accretion and consolidation movements. The sediment properties vary with the intergranular effective stress, according to the depth of the buried material. We used data from the literature to characterize the hydrogeomechanical properties of the San Felice salt marsh. In particular, laboratory tests conducted by Cola et al. 44 using conventional oedometers allowed us to characterize relatively shallow samples collected from the marsh surface up to a 1 m depth.
We performed a number of 100 model runs (Monte Carlo simulations) by varying the input parameter ω, which defines the model accretion rate function sampled from a lognormal distribution (Fig. 4c). A uniform spatiotemporal distribution of ω was assumed as a preliminary hypothesis without any loss in generality for further discussion. The vertical dynamics were highly influenced by accretion rates, providing a quasi-linear response of the system to ω variations. The model calibration against vertical displacement reconstructions in Fig. 4a, b provided the best VD surface match with ω =~2.4 mm year −1 . Model surface VDs were computed within the interval VD = −2.5 ± 0.8 mm year −1 due to the uncertainty associated ω. On average, 60% of the total deposition was offset by shallow compaction of highly compressible sediments.
Projections of the marsh resilience to sea-level rise. We used the model calibration obtained in the previous analysis to predict the possible evolution of the tidal marsh landscape under increasing SLR. The prediction was performed over a period of 80 years, covering the interval of 2020-2100. Various scenarios (SC 1-3) were analyzed based on different accretion rates possibly depositing over the marsh surface. SC-1, SC-2 and SC-3 were characterized by accretion rates of 0 mm year −1 , 2.5 mm year −1 and 5 mm year −1 , respectively (Table 2), representing the cases where sediment transport is completely prevented in the area and equal to or twice the current deposition. Model compaction of the salt marsh in 80 years is estimated to be equal to 0 mm (SC-1), 138 mm (SC-2) and~228 mm (SC-3). By adding the contribution of 1 mm year −1 of deep subsidence to consider the regional pattern related to long-term geological processes 37,45 , the salt marsh elevation between 2020 and 2100 will lose −80 mm in SC-1 and −18 mm in SC-2, whereas it will gain 92 mm in SC-3 ( Table 2). To compare the salt marsh elevation to MSL, we first considered a reference scenario of SLR characterized by a constant value of 1.2 mm year −1 in the period 2020-2100. This value is obtained from data collected by the tide gauge station of Trieste, Northern Adriatic Sea, Italy, in the period 1890-2007 46-48 and confirmed by Zanchettin et al. 49 with a value of 1.2 ± 0.1 mm year −1 for the period from 1872 to 2019, corresponding to~100 mm up to 2100. The result will be a complete submersion of the salt marsh in SC-1 and SC-2 with accretion not able to counterbalance the elevation loss due to compaction and deep subsidence. We estimated land-elevation annual loss equal to~2.2 mm year −1 in SC-1 and~1.4 mm year −1 in SC-2, possibly leading to different stable equilibria such as intertidal flats or subtidal open water 26 . The salt marsh will keep peace with SLR in SC-3 where almost null relative elevation change is expected between 2020 and 2100. To account for the global mean sea level projections under the Representative Concentration Pathway RCP2.6 and RCP8.5 scenarios 50 , a median SLR of 430 mm and 840 mm is expected by 2100 relative to 1986-2005. In scenario SC-2, the salt marsh will loose elevation with respect to the MSL of about 420 mm and 830 mm in RCP2.6 and RCP8.5, respectively (Table 2). However, these projections need two major discussions. The first is related to the correction of the global mean SLR with regional variations due to specific physical processes dominating the causes of sea-level changes. In the case of the Northern Adriatic Sea, arm of the Mediterranean Sea, projections need to consider the particular configuration as semiclosed basin, exchanging water with the Atlantic Ocean through only the Gibraltar Strait. Indeed, these water masses exchanges are a great source of uncertainty for estimating future regional SLR with respect to global projections, with a likely range of MSL rise between 210 mm (low RCP2.6) and 1000 mm (high RCP8.5) by 2100 in Venice relative to 1986-2005 49,51,52 . The second important factor is the relationship between mean SLR and accretion rates. Indeed, the deposition dynamics is strictly related to the relative elevation of the marsh surface with respect to the water level. For high rates of SLR, a more complex model is needed including the dynamics of the surface processes 26,53,54 .

Discussion and conclusion
The city of Venice and its lagoon are seriously threatened by accelerated rates of SLR. In this contribution, we studied salt marsh vertical dynamics based on SLR projections, sediment availability and shallow compaction over the coming century. A physics-based model combined with an integrated analysis of data from a comprehensive monitoring system allowed an investigation of the response of the San Felice salt marsh under three different scenarios, including the potential impact of the MOSE system. MOSE consists of flap-gate barriers installed at the bottom of the three tidal inlets allowing for temporary closure of the lagoon from the sea during high tide events (>110 cm above MSL). The barriers are currently under construction and it is estimated they will be operational by the end of 2021. Repercussions of the MOSE closures might influence lagoon-sea water exchanges, possibly diminishing the sediment supply from the Adriatic Sea into the Venice Lagoon 55,56 . Based on the hypothetical assumption that the sediment flux is totally precluded, we estimated a progressive marsh elevation loss rate of~2.2 mm year −1 (SC-1, Table 2). This is a pessimistic scenario because no permanent closure of the barriers is planned in the future, and sediment supply is also available from resuspension from channels and tidal flats. Moreover, organic sediments accumulate over the marsh surface contributing to accretion rates. At the San Felice salt marsh, organic matter percentages evaluated using Loss On Ignition amounts to 6% along the marsh edge and increases to 17.3% at about 5-7.5 m from the channel 57 . To estimate the amount of organic production, a biomorphological model is needed to account for the biomass productivity of different vegetation species and the twoway feedback on marsh elevation 1 . A recently implemented coupled model between surface and subsurface processes 53 shows that in typical Spartina-dominated environment of the Venice Lagoon, the organic contribution over the total amount of sediments is minimal at the margin of the salt marsh and increases up to 40% moving towards the marsh center. We emphasize that, even though we do not distinguish between inorganic and organic accretion rates, the hydrogeomechanical properties used to characterize the sediments are those typical of the studied salt marsh, representing a spatially average behavior. Further analysis could improve our projection estimates by including the coupling dynamics between biogeomorphological and geomechanical models, and considering the spatial distribution of different properties characterizing the soil deposited on the marsh surface, Scenarios differ for the accretion rate over the marsh surface and the SLR, considering a constant deep subsidence of 1 mm year −1 37 and SLR of 1.2 mm year −1 46-48 . Parameter δ is the marsh elevation loss/gain between 2020 and 2100; parameter Δ is the net elevation loss/gain of the marsh surface with respect to the 2020-MSL. Negative Δ indicates that the vegetated marsh is losing elevation with respect to MSL possibly leading to a different stable equilibrium, e.g., intertidal flat or subtidal open water 26 . a Note that these scenarios need correction for the possible increase of accretion rate as consequence of SLR.
e.g., inorganic sediments or organic matter provided by biomass degradation.
Our results emphasize that the vulnerability of these landforms to SLR is highly dependent on their past formation history. We show that a large influence on loss of elevation relative to MSL is due to compaction of the subsoil, thus neglecting this component in the modeling framework may lead to incorrect predictions of salt marshes resilience.
Holocene compaction may highly impact San Felice salt marsh resilience if sedimentation is maintained at the current rate with an expected 100 mm of SLR by 2100. The scenario becomes worse if sediments are drastically prevented from entering the Venice Lagoon. To counterbalance shallow compaction, deep subsidence, and SLR, a sediment accretion rate of~5 mm year −1 , i.e., twice the one recorded over the last decade, will be required. Our results confirm those of previous studies that demonstrate the powerful influence of peat compaction on the evolution of the late Holocene coastal wetlands in southeastern England, UK 58 . Similarly, compaction of Holocene strata has been found to contribute significantly to high rates of relative SLR in the Mississippi Delta, Louisiana (US) 6 , in the eastern margin of the North Sea 59 , and in the Mekong Delta (Vietnam) 60 . We also highlight that these results, under the reference scenario of SLR rate of 1.2 mm year −1 , should be taken as a starting point to demonstrate the importance of considering shallow subsidence to study longterm marsh evolution. To be able to predict the future state of these environments under the RCP-derived sea-level scenarios 50 , we claim the need of involving a more sophisticated model which also account for the possible variation of the accretion rate, which depends on sediment availability but also on the water level above the marsh surface.
Although our model uses accretion rates as time-constant input value, this is a first approximation that can be released as new information on the sediment spatio-temporal variability in the Venice Lagoon will become available. However, we provide a quantification of uncertainty of the model response with a Monte Carlo-based analysis accounting for the variability in the measured deformation and accretion rates. We found that the most likely average accretion rate over the Holocene was to 2.4 mm year −1 with a ±0.8 standard deviation. This uncertainty, intrinsically related to the consolidation state of the recently deposited strata, will obviously affect the capability of the salt marsh to be resilient against SLR.
Our research emphasizes that the features of the Holocene deposits that constitute the marsh body play a major role in the future resilience of these transitional landforms. It follows the importance of characterizing the hydrogeomechanical behavior of highly compressible sediments, in particular within the shallowest meters of soil where major consolidation takes place. Laboratory experiments carried out with oedometer compression tests on cores sampled at the Cowpen Marsh, Tees Estuary, UK shows that measured material properties are largely influenced by local factors, including climatic, sedimentological, geomorphic and hydrographic conditions 61 . We highlight the importance of considering the past history of the marsh to study its long-term evolution in combination with other aspects controlling coastal marsh dynamics, such as the regional hydrological setting 62 . Moreover, high levels of disturbance when core sampling are usually associated with these soil types due to the presence of very soft material and the complex root system in the first 30 cm of topsoil. This suggests the need for additional appropriate investigations to determine the compression behaviors of salt marshes using, for example, in situ loading tests 63 .
Based on our findings, planning and designing salt marsh restoration should also consider the high probability of compaction experiencing the Holocene sequence, as result of adding "extra" sediment load on top of existing marshes. Depending on the characteristic properties of the nourishment material and the marsh body, the accretion rates, which account for the combination of accretion and compaction, can be disparate and the elevation of the marsh surface may vanish if tidal silty-sand sediments are placed over low-density organic soils that are prone to compaction.

Materials and methods
Estimation of long-term shallow subsidence. Modeling approach and notation is here described following previous works 32, 64 . We assume a vertical-cross section of a marshland and denote it as the domain Ω t & R 2 at some instant t 2 I ¼0; T½, with T being the final simulation time, in a x − z reference system. This section has fixed length L in the x-coordinate but varies in time t along the z-coordinate. The boundary of the domain Ω t , Γ t , is subdivided into four portions: Γ N,b is the bottom boundary representing the impervious marshland basement, Γ t N;r is the right boundary and is assumed to be a symmetry axis for the cross-section (no cross-flow is allowed), Γ t D;l is the left boundary representing the marshland border that is in equilibrium with the tidal creek. This means that no overpressure with respect to the hydrostatic distribution is allowed. Γ t D;t is the top boundary representing conditions of hydrostatic pressure. We denote as Γ t N ¼ Γ N;b ∪ Γ t N;r the boundary portion where Neumann conditions are imposed and Γ t D ¼ Γ t D;l ∪ Γ t D;t the boundary portions with Dirichlet conditions. These conditions are prescribed at every instant t ∈ I to follow the mesh updates. Moreover, let conditions Γ t D ∪ Γ t N ¼ Γ t and Γ t D \ Γ t N ¼ ; hold true. The governing equations of the pore pressure variation with respect to hydrostatic conditions, p, in Ω t ¼ Ω t I is described by the fluid mass balance in a deforming and saturated porous medium. Note that the bar above a set identifies the union of the set with its boundary, i.e., Ω t ¼ Ω t ∪ Γ t . Following the original rigorous theory of 1D flow in an elastic saturated compacting porous medium 65,66 and further development to 2D domains 32 , the mass balance equation is written under the assumptions of (i) large solid grain motion by dropping the hypothesis of infinitesimal displacement of the grains, (ii) Darcy's law applied to the relative velocity between fluid and solid grains and (iii) incompressible solid grains. Marshland accretion depends on the sediment accretion rate, ω(x, t), which causes a variation in time of the total vertical stress,σ z ðx; tÞ, relative to the initial total vertical stress. The vertical displacement, u(x, t), is obtained as a function of the variation in the vertical effective stress, σ z (x, t), which is related tô σ z ðx; tÞ by Terzaghi's principle 67 . Hence, the strong form of the initial/boundary value problem may be stated as follows. Given ω : Ω t I ! R, find p : Ω t ! Rσ z : Ω ! R, and u : Ω t ! R such that: where ψ is the storage coefficient, κ is the rank-2 hydraulic conductivity tensor, γ the specific weight of water, ζ the soil oedometric compressibility; ϕ 0 is the porosity at the initial reference vertical stress σ z,0 , γ s the specific weight of the solid grains; α is the classical compressibility as a function of the soil Young's modulus and Poisson's ratio, E and ν, respectively;∇⋅ and ∇ are the divergence and gradient operators, respectively; and the superposed dot, _ ðÞ, denotes a derivative with respect to time. The use of the Lagrangian approach allows us to replace the total Eulerian derivative on a moving particle with partial derivatives in time. The governing PDE is nonlinear when equipped with appropriate constitutive laws for the material under consideration. Having defined the classical compressibility α as (1 − 2ν)(1 + ν)/[(1 − ν)E], the following relationships hold for ψ and ζ 32,67 : where β represents the groundwater compressibility and ϕ represents the soil porosity.
Hydrogeomechanical properties and model set-up of San Felice salt marsh. Hydrogeomechanical properties of a salt marsh are not easy to determine. The high-porosity and low-density sediments that usually characterize these marshlands are difficult to sample for laboratory compression tests. At the San Felice salt marsh some attempts to sample and test the compressibility of shallow soil has been conducted in the past 44 . The mechanical behavior of the material is assumed to follow virgin loading conditions. The void ratio, e, is a function of σ z , following with C c the compression index, e 0 the initial void ratio, and σ z,0 the initial vertical effective stress. The parameters e 0 and C c ( Table 3) are used to characterize the Holocene unit according to previous laboratory results 44 . The values of σ z are computed by the model depending on depth z. As consequence, e and ϕ vary with σ z according to the empirical relationships that depend on the sediment type and the loading conditions of the sediments (Supplementary Fig. S1). κ is constant with σ z ( Table 3).
The initial domain is a vertical cross-section discretized by 200 nodes and 198 triangular elements. We use a horizontal scaling factor of 3. The horizontal and vertical discretizazion is 0.1-m in x and z directions. The new elements are added on top of the existing domain are 0.1 m-high. A constant time-step size Δt = 10 years is prescribed (Table 3).

Data availability
The authors declare that the data supporting the findings of this study are available on Figshare repository with the identifier https://doi.org/10.6084/m9.figshare.14192480.

Code availability
The code NATSUB2D that support the findings of this study is available from the corresponding author upon request.