**Table 4.** Expected accumulated damage and adaptation strategies linked to sea-level rise percentiles (in million EUR) in 2100 for the 10 Iberian coastal cities with the largest populations under the RCP 8.5 scenario.

#### 4.2.2. Adapting by Fixing a Certain Height

If defence heights were to be decided exogenously, for instance, at 2 m and 3 m by 2100, the average adaptation costs for the top 10 cities of the ranking would range between 282 and 1225 million euros for a protection infrastructure of 2 m (Table 5). In the case of a 3 m defence, average costs could reach 400–2400 million euros, but even in this case, the investment can be well justified to avoid the accumulated damage. Note that if standard height defences of 2 m or 3 m are to be implemented, the probability of damage in the year 2100 is rather small.

**Table 5.** Expected accumulated damage and adaptation strategies with standard height defences in 2100 (costs shown in million euros) for the 10 largest Iberian coastal cities, under RCP 8.5.

#### **5. Discussion**

#### *5.1. The Economic Damage of Sea-Level Rise*

In this paper, we propose a model for the calculation of accumulated expected damage under uncertainty. The model uses an expanded version of the projected sea-level rise percentiles [9] and a deterministic function that depends on flood height (in our case, sea-level height) [19]. The combination of both provides a group of stochastic damage distributions that allow us to calculate annual expected damage and, consequently, also the expected accumulated damage. Accumulated average damage in all 62 cities of our sample exceed 1700 million euros in 2030, and the difference between the results under RCP 4.5 and RCP 8.5 is 45 million euros. However, by the end of the century, accumulated damage reach 41,200 million euros under RCP 4.5 and more than 50,400 million euros in the highest emission scenarios (RCP 8.5). Accumulated damage increase by more than 20 times by 2100 and the difference between RCP 4.5 and RCP 8.5 could reach 10 billion euros (Table S8 and Table S9).

Furthermore, our results show that limiting the assessments to median values (50th percentile) significantly underestimates the potential damage. While differences between median and average values might be small by 2030, they increase considerably by one order of magnitude by the end of the century. For example, in Barcelona the difference between median and average damage is 0.5 million euros in 2030, but increases to almost 4 million euros in 2100 (Table 2). This difference occurs because the shape of the sea-level rise distribution is positively skewed, i.e., the distribution is not symmetrical but has a heavier upper tail The shape of the distribution of the damage is also another cause of important underestimations in the medium term when calculating expected damage. These findings are consistent with previous studies that warned about the larger potential risk of climate-induced impacts.

Damage costs do not account for coastal erosion or ecosystem loss, impacts that are both expected to be relevant in the Iberian Peninsula. Additionally to their intrinsic value, ecosystems provide a number of services to people that can be monetised [50]. The value of services such as storm protection or erosion control provided by salt marshes and mangroves worldwide is estimated to be 194,000 \$/ha/year [51], which has not been accounted for in this study. Of course, one should acknowledge that other important factors such as shoreline evolution have to be integrated into a cost-benefit analysis to account for the complexities of coastal defence solutions, as argued earlier by other authors [52]. In this case, we have simplified this part of the analysis to focus on the impact of uncertainty on the economic estimates but by no means do we neglect the need to integrate those factors.

#### *5.2. Comparing the Costs of Inaction to Protection Investment Costs*

Two different adaptation options have been explored in this study: the first is based on tailoring the level of risk cities might decide to protect themselves from. This is illustrated through the assumption that coastal or city managers decide to protect each city for 95% or 99.5% of the cases under RCP 8.5 by 2100. In the first case, with protection strategies to face 95% of the cases, investment needs in the largest coastal cities by population would range between 87 (Porto) and 1015 (Bilbao) million euros. Investment costs almost double if a higher protection threshold (99.5%) is considered, and these would vary between 164 (Porto) and 1015 (Bilbao) million euros. The second adaptation option explored the investment costs needed for certain protection heights (2 and 3 m). In the case of 3 m protection infrastructures, construction costs range between 404 and 2368 million euros.

An important finding of this analysis is that adaptation costs in the long term are much smaller than the increasingly expensive costs of inaction, and thus, investing in adaptation is a good decision when comparing costs and benefits of the actions, a finding in line with previous research carried out in this area of study [53].

Of course, one should acknowledge that many other adaptation options that have not been considered in this paper may exist, and that many other criteria other than economic costs and benefits may be worth being taken into account when making such decisions. Another clear finding is that adaptation investment needs will be smaller under more favourable scenarios (RCP 2.6 and RCP 4.5). However, even under the most optimistic scenario (RCP 2.6), expected damage will be very important in the second half of the century. In other words, despite strong emission reduction efforts being really effective, adaptation will be inevitable. This is to highlight the importance of the need to implement both mitigation and adaptation policies.

Our findings show that not adapting to climate change is not, by any means, a good strategy in the medium and long term. Note that, even if many types of adaptation options should be considered when protecting coasts from climate risks, the building of defence infrastructures that require strong investments would be justified based on the accumulated damage by 2100, which are much greater than the investment costs of protection. Moreover, the adaptation strategies considered in this paper are shown to be very cost-effective in every city studied and reduce considerably the probability of experiencing high economic damage by the end of the century.

#### *5.3. Policy Implications*

In this paper we propose a method to avoid the underestimation of climate risk and estimate what the size of it may be for the case of 62 cities in the Iberian peninsula. We have used two infrastructure-based adaptation options for illustrative purposes. This has been done because the data needed for the analysis is available in these cases, and not because we argue that these solutions should be the ones implemented. Results show that the underestimation of damage can be very important if we consider median values instead of the average (expected) damage. Therefore, we argue that strong efforts should be made to understand the full probability distributions of sea-level rise and the consequent economic damage, for each city, emission scenario and year. Not accounting for all these may lead to investment decisions (and even maladaptation) that underestimate future risks and cannot respond to them.

The method can also be applied to adequately assess the cost effectiveness of other adaptation options when data becomes available. The purpose of the paper has been to propose a methodology and illustrate the case for sea-level rise and hard adaptation measures for which good data was available. Other pieces of research should help us to better consider alternative adaptation options or even other impacts related to the adaptation options considered. Nonetheless, the main message of this paper remains, not accounting for the full distribution of sea level rise as well as damage costs is a clear and very significant underestimation of climate risks that may lead to inadequate policy decisions.

**Supplementary Materials:** A dataset with additional results for all 62 cities, scenarios and years is available online at http://www.mdpi.com/2073-4441/12/4/1220/s1.

**Author Contributions:** Conceptualization, L.M.A.; formal analysis, L.M.A. and E.S.d.M.; methodology, L.M.A.; writing—original draft, L.M.A., E.S.d.M. and I.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is supported by the Basque Government through the BERC 2018-2021 program and by the Spanish Ministry of Economy and Competitiveness MINECO through BC3 María de Maeztu excellence accreditation MDM-2017-0714. Additionally, Luis M. Abadie is grateful for financial support from The Spanish Ministry of Science and Innovation (ECO2015-68023). E. Sainz de Murieta acknowledges funding from the Basque Government (Postdoctoral Fellowship grant no. POS\_2018\_2\_0027). Ibon Galarraga received financial support from the European Union's Horizon 2020 research and innovation programme under grant agreement No 776479 for the project CO-designing the Assessment of Climate CHange costs.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Mitigating the E**ff**ects of Sea-Level Rise on Estuaries of the Mississippi Delta Plain Using River Diversions**

**Eric D. White 1,2,\* , Ehab Meselhe 2 , Denise Reed 3 , Alisha Renfro 4 , Natalie Peyronnin Snider 5 and Yushi Wang 6**

Received: 7 August 2019; Accepted: 23 September 2019; Published: 28 September 2019

**Abstract:** Using the Mississippi River as a tool for restoration has been a key element of restoration planning in Louisiana for decades. The results of allowing river water and sediment back into the coastal system are manifested in a number of places in present day Louisiana, with additional plans for large scale sediment and water diversions from the Mississippi River. Many previous numerical modeling studies have focused on sediment delivery to Louisiana estuaries. This study examines the effects of river diversions on salinity gradients in receiving estuarine basins. The Integrated Compartment Model, a planning-level model that simulates multi-decadal change in estuarine hydrodynamics and wetland systems under assumed sea-level rise scenarios, was used to assess the estuarine salinity gradient under potential management regimes. The simulations for current conditions are compared to a future 50-year simulation with additional diversions, as well as cases with a variety of diversion options. This modeling analysis shows that without additional action, 50-years of sea-level rise could result in substantial increases in salinity throughout the Mississippi Delta Plain estuaries. This can be largely offset with additional large river diversions which can maintain variable salinity gradients throughout the estuary basins.

**Keywords:** Mississippi River Delta; Louisiana; salinity; sediment; wetland loss; estuary; diversion; sea-level rise; environmental planning; coastal restoration

#### **1. Introduction**

The impacts from projected future sea-level rise are expected to be extensive for the World's River Deltas, with the loss of land and habitat diversity threatening the inhabitants, resources and ecology of these low-lying coastal areas [1,2]. Many deltas have already experienced broad change, resulting in a loss of habitat extent and diversity due to human activities [3]. The Mississippi River Delta is a river-dominated delta that was built in a series of successive lobes over the last 7500 years [4]. However, over the past 80 years, Louisiana has lost 3800 km2 of land on the Mississippi River's deltaic coastal plain [5] due to a variety of human-induced and natural factors, such as rising sea levels, levee construction along the Mississippi River, sediment compaction due to pressure

from overlying sediments and organic decomposition, fluid withdrawal, decreased fluvial sediment discharge, hurricanes, and widespread hydrologic alteration [6,7].

Over time, the interaction between the Mississippi River and the delta plain has been progressively managed [8,9] such that there are presently a few uncontrolled connections upstream of the Birdsfoot Delta [10]. Reconnecting the river with the delta plain is now widely recognized as an essential tool in the future sustainability of the wetland-dominated coastal ecosystem [11–14], and this approach has been central to restoration planning in Louisiana for decades [15]. Two freshwater diversions, Caernarvon and Davis Pond (Table 1) were built in the late 20th century specifically for the management of the estuarine gradient to support conditions for fisheries in the receiving basins. Future managed reconnections are usually in the form of sediment diversions. These are structures with gates built in the existing river levee system that would be operated to allow river water, sediment and nutrients to flow into the wetlands adjacent to the river to build new land and sustain existing wetlands. The results of allowing river water and sediment back into the coastal system are manifest in a number of places including the delta of the Wax Lake Outlet, mining of accumulated sediment in the Bonnet Carre Spillway, crevasse splays in the Birdsfoot Delta and revitalized marshes adjacent to freshwater diversions and siphons [14]. At the same time, adjustments in the estuarine gradient as freshwater is introduced are expected to have some effect on the habitat for important commercial fisheries, including oysters [16], brown shrimp [17,18], and other commercial species [19].

**Table 1.** Locations of existing or planned connections between the Mississippi River and estuarine basins. Barataria (BA), Breton (BR) and Pontchartrain (PO) identify the receiving basin. ID tracks the locations on Figure 1.

**Table 1.** *Cont*.

Louisiana has a long history of coastal restoration and protection planning. In response to the past and the growing threat of future coastal land loss due to rising sea levels, subsidence, and altered hydrology, Louisiana has developed a coastwide plan for coastal restoration and storm surge-based risk reduction over the next 50 years. The implementation of this plan would reduce land loss and expected annual damages, compared to a future with no action [20–22]. Through the 2017 Coastal Master Plan process, the restoration projects were analyzed and selected based on their ability to build or sustain land against a potential future landscape that incorporated up to 0.8 m of eustatic sea-level rise over the next 50 years. With the effort and investment envisaged in the 2017 plan, nearly 3000 km2 of land is projected to be built or sustained that may have otherwise been lost. The expected annual damages from storm surge-based flooding are projected to be reduced by \$12.2 billion at year 50 [22]. The 2017 plan selected a variety of restoration project types including river sediment diversions.

A number of sediment diversions from the Mississippi River have been identified to support land building and wetland sustainability [21,22]. The focus on sediment recognizes the value of that resource. Freshwater diversions of various scales have been used for decades in Louisiana to manage estuarine salinity gradients. The 2017 Coastal Master Plan also introduced the concept of using diversions at strategic and suitable locations that operate only under low-moderate river discharges to provide supplemental freshwater to estuarine wetlands, especially freshwater wetlands, during the times of the year when salinities often increase. The ecosystem diversity can be important outcomes of strategically placed riverine diversions, as well as land-building and land maintenance.

Moving water and sediment from the river into the estuarine basins results in ecosystem change relative to current or past conditions. However, the forward-looking analysis that considers future change, especially that associated with sea-level rise is needed.

The analysis described here explores how the use of river freshwater and sediment resources at multiple locations, individually or in combinations, year-round or seasonally, and in different magnitudes can be used to sustain the coastal ecosystem. It uses numerical experiments to consider several strategies for diversion operation and shows how the effects of diversion vary over time as the sea-level rises, between the years as the river fluctuates, and among estuarine basins with different landscape settings.

#### **2. Methods**

This analysis utilizes tools developed for the 2017 Coastal Master Plan to further explore the potential of the Mississippi River as a tool to support the sustainability of the Mississippi Delta Plain and the economic activities that depend on it. The Integrated Compartment Model (ICM) is a planning-level model that was developed by integrating into a single modeling platform several models that had previously been used for coastal zone planning and research in Louisiana [23]. The models that are included as subroutines within the ICM framework include a hydrologic and hydraulic model [24], a vegetation dynamics model [25], a wetland morphology/elevation change model [26], as well a barrier island morphology model and several receptor models that summarize the hydrologic and landscape conditions into numerous habitat indices and decision-making metrics. The hydraulic and hydrologic model, ICM-Hydro, simulates on a sub-daily timestep: water level (stage), flow rate, salinity, water temperature, suspended sediment concentration, sediment deposition and resuspension within open water areas, sediment deposition on the marsh surface, and a variety of water quality/nutrient constituents. The seasonal hydrologic conditions and salinity are then utilized by the vegetation dynamics model, ICM-LAVegMod, to predict the relative likelihood of wetland vegetation species being present under such conditions. The patterns of inundation, salinity penetration, sediment deposition and organic matter accretion are all taken into account within the wetland morphology/elevation change model, ICM-Morph, to predict whether wetland areas remain as viable emergent vegetation or collapse into open water due to either persistent inundation stress or acute salinity stress. Finally, as the coastal wetlands change in vegetative cover, and the resultant changes in elevation and land/water composition, the hydraulic network of the estuary is updated prior to the next model year in ICM-Hydro, to simulate the changing hydraulic and hydrologic conditions that occur as a result of collapsing wetland areas.

Prior to any modeling analysis, the ICM was calibrated using observational data from 2010 to 2013, and then validated on observations from 2006 to 2009. These time periods correspond to the development and widespread availability of data from the Coastwide Reference Monitoring System (CRMS), which collects hourly water level and salinity data across coastal Louisiana [27]. The water level predictions were calibrated with a goal that 80% of compartments should have a prediction bias of less than 0.15 m in the daily mean water level prediction. The magnitude of this bias corresponds approximately to the error in the underlying topographic digital elevation model (DEM) used in this analysis, which varies from a RMSE of 0.07 to 0.3 m [28]. The hydraulic link capacity and bed roughness were adjusted if the observed stage signal and amplitude were not in satisfactory agreement with the model. Due to the apparent datum inconsistencies in certain observed stations (as well as clear patterns of hydraulic controls influencing observed water levels), certain observed data were excluded when assessing the overall model fit. These inconsistent datasets were still used to visually compare modeled and observed hydrographs, but they were excluded from any aggregate model performance statistics.

Salinity calibration was also conducted for the same 2010 through 2013, time frame. The model agreement was defined as 80% of the compartments with a bias in salinity prediction of less than 1 ppt. Coastwide, 76% of the model-observation pairs met this goal. In the Mississippi River Delta (Pontchartrain, Breton and Barataria basins, Figure 1), 87% of the model-observation pairs met the bias target of 1 ppt or less (refer to Tables S2 and S3 and Figures S2–S5 in Supplementary Material). To attain this, a model parameter representing a combined dispersion diffusion coefficient, Exy, was adjusted if the observed salinity signal and amplitude were not matched by the model.

90

*Water* **2019**, *11*, x FOR PEER REVIEW 4 of 26

65 **Figure 1.** Location map showing the estuarine basins and the location of existing or planned 66 connections between the Mississippi River and estuarine basins. See Table 1 for more information. **Figure 1.** Location map showing the estuarine basins and the location of existing or planned connections between the Mississippi River and estuarine basins. See Table 1 for more information.

67 Louisiana has a long history of coastal restoration and protection planning. In response to the 68 past and the growing threat of future coastal land loss due to rising sea levels, subsidence, and altered 69 hydrology, Louisiana has developed a coastwide plan for coastal restoration and storm surge-based 70 risk reduction over the next 50 years. The implementation of this plan would reduce land loss and 71 expected annual damages, compared to a future with no action [20–22]. Through the 2017 Coastal 72 Master Plan process, the restoration projects were analyzed and selected based on their ability to 73 build or sustain land against a potential future landscape that incorporated up to 0.8 meters of 74 eustatic sea-level rise over the next 50 years. With the effort and investment envisaged in the 2017 plan, nearly 3000 km2 75 of land is projected to be built or sustained that may have otherwise been lost. 76 The expected annual damages from storm surge-based flooding are projected to be reduced by \$12.2 The higher Exy values allowed for more exchange between compartments and lower values allowed for less exchange between compartments. In addition to adjusting the dispersion-diffusion coefficient, the salinity predictions were improved by the addition of a term in the hydrodynamic code that replaced the original central-difference method used for salinity convection with a first-order upwinding scheme [29]. This scheme was used for compartments with typical river/canal flows, where upwinding was an appropriate approximation, as compared to slower estuarine/marsh flow regimes where the upwinding scheme was not applied. This addition of the upwinding technique greatly increased the stability of salinity predictions. Further discussion and analysis of the ICM methodology, calibration, validation, model performance and analysis are provided in the Supplementary Materials as well as several peer reviewed articles [24–26,30–34] and technical reports [35–41].

77 billion at year 50 [22]. The 2017 plan selected a variety of restoration project types including river 78 sediment diversions. 79 A number of sediment diversions from the Mississippi River have been identified to support 80 land building and wetland sustainability [21,22]. The focus on sediment recognizes the value of that 81 resource. Freshwater diversions of various scales have been used for decades in Louisiana to manage 82 estuarine salinity gradients. The 2017 Coastal Master Plan also introduced the concept of using 83 diversions at strategic and suitable locations that operate only under low-moderate river discharges 84 to provide supplemental freshwater to estuarine wetlands, especially freshwater wetlands, during 85 the times of the year when salinities often increase. The ecosystem diversity can be important 86 outcomes of strategically placed riverine diversions, as well as land-building and land maintenance. This paper focuses on six sets of analyses. Five of the six simulations were analyzed using the 50-year flow hydrograph used in the 2017 Coastal Master Plan analysis, whereas the current conditions simulation used the observed flow hydrograph for 2006 through 2013. Three distinct Mississippi River flow conditions were identified from the historical river discharge data from 1964 to 2013 and were used in single year model runs for some of the analysis. Figure 2 shows the daily mean flow in the Mississippi River at Tarbert Landing, from 1964 to 2013. Figure 3 shows the annual mean flow anomaly (i.e., deviation from the long-term mean) for the Mississippi River for this 50-year flow hydrograph. A low river year (1988, elapsed year 25), a high river year (1993, elapsed year 30), and an average river year (1992, elapsed year 29) were selected between year 25–35. These river years were used to examine the effects of project configurations under a range of river flow conditions for some of the simulations described here.

87 Moving water and sediment from the river into the estuarine basins results in ecosystem change 88 relative to current or past conditions. However, the forward-looking analysis that considers future Table 1 shows the operation regimes for diversions and other river outlets to estuarine basins. The operational flow rates for most passes, cuts, and diversions on the Mississippi River were based

89 change, especially that associated with sea-level rise is needed.

152 [35–41].

164

170

upon the flow in the river at each respective pass/diversion location. Therefore, the operation of any diversion was impacted by all operated diversions located upstream. 161 30), and an average river year (1992, elapsed year 29) were selected between year 25–35. These river 162 years were used to examine the effects of project configurations under a range of river flow conditions 163 for some of the simulations described here.

*Water* **2019**, *11*, x FOR PEER REVIEW 6 of 26

142 The higher Exy values allowed for more exchange between compartments and lower values 143 allowed for less exchange between compartments. In addition to adjusting the dispersion-diffusion 144 coefficient, the salinity predictions were improved by the addition of a term in the hydrodynamic 145 code that replaced the original central-difference method used for salinity convection with a first-146 order upwinding scheme [29]. This scheme was used for compartments with typical river/canal 147 flows, where upwinding was an appropriate approximation, as compared to slower estuarine/marsh 148 flow regimes where the upwinding scheme was not applied. This addition of the upwinding 149 technique greatly increased the stability of salinity predictions. Further discussion and analysis of the 150 ICM methodology, calibration, validation, model performance and analysis are provided in the 151 Supplementary Materials as well as several peer reviewed articles [24–26,30–34] and technical reports

153 This paper focuses on six sets of analyses. Five of the six simulations were analyzed using the 154 50-year flow hydrograph used in the 2017 Coastal Master Plan analysis, whereas the current 155 conditions simulation used the observed flow hydrograph for 2006 through 2013. Three distinct 156 Mississippi River flow conditions were identified from the historical river discharge data from 1964 157 to 2013 and were used in single year model runs for some of the analysis. Figure 2 shows the daily 158 mean flow in the Mississippi River at Tarbert Landing, from 1964 to 2013. Figure 3 shows the annual

160 year flow hydrograph. A low river year (1988, elapsed year 25), a high river year (1993, elapsed year

165 **Figure 2.** Daily mean flow in the Mississippi River at Tarbert Landing, 1964 to 2013. *Water* **2019**, *11* **Figure 2.** , x FOR PEER REVIEW Daily mean flow in the Mississippi River at Tarbert Landing, 1964 to 2013.

7 of 26

171 **Figure 3.** Annual mean flow anomaly from 1964 to 2013. **Figure 3.** Annual mean flow anomaly from 1964 to 2013.

172 All of the analyses presented here (with the exception of the current conditions run, which used 173 observed boundary conditions) were conducted using the medium scenario for the 2017 Coastal 174 Master Plan [35]. The medium scenario uses 0.63 meters of eustatic sea-level rise (ESLR) from 2015 175 through the end of 2064. The acceleration term increased over time, and the rate corresponded to 1.5 176 meters of ESLR by 2100 compared to the 1992 sea level. The relative sea-level rise rates were assigned 177 based upon these ESLR rates and subsidence values, which varied spatially across the model domain 178 [42]. For the area included in this study, the subsidence ranges from 2.6 mm/yr. in the Upper 179 Pontchartrain basin to 17 mm/yr. in the Birdsfoot Delta (refer to Figure S1 in Supplementary 180 Materials: ICM Overview, for a map showing subsidence zones for coastal Louisiana used in this All of the analyses presented here (with the exception of the current conditions run, which used observed boundary conditions) were conducted using the medium scenario for the 2017 Coastal Master Plan [35]. The medium scenario uses 0.63 m of eustatic sea-level rise (ESLR) from 2015 through the end of 2064. The acceleration term increased over time, and the rate corresponded to 1.5 m of ESLR by 2100 compared to the 1992 sea level. The relative sea-level rise rates were assigned based upon these ESLR rates and subsidence values, which varied spatially across the model domain [42]. For the area included in this study, the subsidence ranges from 2.6 mm/year in the Upper Pontchartrain basin to 17 mm/year in the Birdsfoot Delta (refer to Figure S1 in Supplementary Materials: ICM Overview, for a map showing subsidence zones for coastal Louisiana used in this analysis).

181 analysis). 182 A set of simulations was designed to explore a range of uses of Mississippi River water and 183 sediment. The locations of existing and planned connections from the Mississippi River to the 184 estuarine basins are shown in Figure 1 and their configurations are shown in Table 1. A set of simulations was designed to explore a range of uses of Mississippi River water and sediment. The locations of existing and planned connections from the Mississippi River to the estuarine basins are shown in Figure 1 and their configurations are shown in Table 1.

185 The six simulations conducted for this analysis were: The six simulations conducted for this analysis were:

193 rise can be compared.

186 1. Current Conditions. Before the ICM was used to assess future hydrologic and 187 landscape conditions throughout coastal Louisiana, the ICM-Hydro subroutine was 188 run from 2006 to 2013 using observed salinity and water level data collected across 189 coastal Louisiana. The descriptions of the data used for this run and discussions on 190 general model performance and limitations can be found in previously published 191 literature [28,30–32]. This run is used to illustrate the current pattern of salinity 1. Current Conditions. Before the ICM was used to assess future hydrologic and landscape conditions throughout coastal Louisiana, the ICM-Hydro subroutine was run from 2006 to 2013 using observed salinity and water level data collected across coastal Louisiana. The descriptions of the data used for this run and discussions on general model performance and limitations can be found in previously published literature [28,30–32]. This run is used to illustrate the current pattern of salinity fluctuations within the estuaries, against which simulations that include sea-level rise can be compared.

192 fluctuations within the estuaries, against which simulations that include sea-level

195 a. Future Without Action Run (FWOA). The ICM was run for a future 50-year 196 period without any potential future projects included. This simulation is 197 termed the future without action (FWOA). In the FWOA, all sediment and 198 freshwater diversions that were already built were implemented in the 199 model. A full discussion of all boundary conditions assumed for this 50-year 200 future condition is available in Brown et al., [28]. This run is used to 201 demonstrate the effects of relative sea-level rise on the estuarine salinity

203 b. Base Run with Both Mid-Basin Sediment Diversions (Base Run). This 50-year 204 model run includes all diversions that are currently in some phase of the 205 engineering and design process (e.g. Mid-Barataria sediment diversion,

202 gradient in the absence of additional freshwater inflow.

	- a. Future Without Action Run (FWOA). The ICM was run for a future 50-year period without any potential future projects included. This simulation is termed the future without action (FWOA). In the FWOA, all sediment and freshwater diversions that were already built were implemented in the model. A full discussion of all boundary conditions assumed for this 50-year future condition is available in Brown et al., [28]. This run is used to demonstrate the effects of relative sea-level rise on the estuarine salinity gradient in the absence of additional freshwater inflow.
	- b. Base Run with Both Mid-Basin Sediment Diversions (Base Run). This 50-year model run includes all diversions that are currently in some phase of the engineering and design process (e.g., Mid-Barataria sediment diversion, Mid-Breton sediment diversion, and Maurepas freshwater diversion) and represents a likely future condition against which the need for additional management of river-estuarine can be assessed. In the base run, these diversions were implemented in addition to the existing diversion included in FWOA (Table 1). This run is used to demonstrate the effects of the diversion projects which are expected to be built within the next 5 years in the context of future sea-level rise.
	- c. All 2017 Master Plan Diversions (MP17). This simulation shows the potential effect of all 2017 Coastal Master Plan diversions working together over a 50-year period. This includes 5 diversions in addition to those included in the base run (Table 1).

To examine the effect of diversions over time as the relative sea-level rise progresses, the salinity gradients within the receiving basins plots were constructed showing the area of the basin (wetland and open water combined) subject to different salinity levels. Six bins were identified representing different salinity ranges. These salinity ranges were selected on the basis of salinity tolerance of seven wetland species which have previously been used to assess the effects of Mississippi River diversions [44]. The range in salinity values in which these species are optimally located (as defined by the LAVegMod model rules [25]) are shown in Figure 4. The daily salinity values for the growing season for the Louisiana coastal wetland vegetation were used to calculate the area of the basin, experiencing the salinity range for each bin in each year of the simulation. The model output for each ICM-Hydro compartment, which included both open water and marsh areas, was summarized based on three ecoregions shown in Figure 1.

*Water* **2019**, *11*, x FOR PEER REVIEW 9 of 26

250 **Figure 4.** Salinity ranges considered optimal for seven indicator wetland vegetation species. **Figure 4.** Salinity ranges considered optimal for seven indicator wetland vegetation species.

#### 251 **3. Results 3. Results**

249

#### 252 *3.1. Current Conditions 3.1. Current Conditions*

253 A comparison across the basins (Figure 5) shows the differences in the penetration of salinity 254 into the basins. Both Pontchartrain and Barataria show extensive areas with salinity 16–35 ppt in all 255 years, while for four of the years this salinity range does not occur in Breton. This is due to the 256 proximity of Breton to the Mississippi River and the many small connections (see further analysis of 257 the effects of these connections below). A comparison across the basins (Figure 5) shows the differences in the penetration of salinity into the basins. Both Pontchartrain and Barataria show extensive areas with salinity 16–35 ppt in all years, while for four of the years this salinity range does not occur in Breton. This is due to the proximity of Breton to the Mississippi River and the many small connections (see further analysis of the effects of these connections below).

258 Breton is also different at the lower end of the salinity range. Both Pontchartrain and Barataria 259 extend inland some distance, are larger in size (increasing freshwater inputs through precipitation) 260 and receive some municipal stormwater runoff. Pontchartrain also received freshwater inflow from 261 tributary streams. Both Barataria and Pontchartrain include areas <0.6 ppt in all years in this run. 262 However, the conditions in 2006, 2007 and 2012 eliminate this zone in Breton. Breton is also different at the lower end of the salinity range. Both Pontchartrain and Barataria extend inland some distance, are larger in size (increasing freshwater inputs through precipitation) and receive some municipal stormwater runoff. Pontchartrain also received freshwater inflow from tributary streams. Both Barataria and Pontchartrain include areas <0.6 ppt in all years in this run. However, the conditions in 2006, 2007 and 2012 eliminate this zone in Breton.

263 The period shown in Figure 5 included both a high flood year (2011) and a low river flow year 264 (2012) as shown in Figure 3. However, there is a substantial change in the salinity distribution in only 265 one of the three basins. The high flow in 2011 does not increase the area <0.6 ppt in any basin 266 indicating the limited connections between the Mississippi River and estuarine basins in the current 267 condition. However, Breton which has several small connections, does show a substantial decrease 268 in the area with salinities <2 ppt in 2012 during the low flow year. The period shown in Figure 5 included both a high flood year (2011) and a low river flow year (2012) as shown in Figure 3. However, there is a substantial change in the salinity distribution in only one of the three basins. The high flow in 2011 does not increase the area <0.6 ppt in any basin indicating the limited connections between the Mississippi River and estuarine basins in the current condition. However, Breton which has several small connections, does show a substantial decrease in the area with salinities <2 ppt in 2012 during the low flow year.

#### 269 *3.2. Simulations of Future Conditions with Varying Levels of Freshwater Inflow 3.2. Simulations of Future Conditions with Varying Levels of Freshwater Inflow*

270 Three of the simulations progressively increase the exchange between the Mississippi River and 271 estuarine basins (Table 1). Table 2 summarizes the differences in the distribution of salinity for 272 decadal intervals during the 50-year model runs, showing the relative area of the basins for three 273 different salinity ranges (based on the six bins described above) for ease of comparison. Three of the simulations progressively increase the exchange between the Mississippi River and estuarine basins (Table 1). Table 2 summarizes the differences in the distribution of salinity for decadal intervals during the 50-year model runs, showing the relative area of the basins for three different salinity ranges (based on the six bins described above) for ease of comparison.

275 (FWOA), the base run and MP17 simulations. **FWOA Table 2.** Relative area of each basin subject to three salinity ranges by decade for future without action (FWOA), the base run and MP17 simulations.

274 **Table 2.** Relative area of each basin subject to three salinity ranges by decade for future without action

**Table 2.** *Cont*. *Water* **2019**, *11*, x FOR PEER REVIEW 10 of 26

277 **Figure 5.** April to September salinity region variation in Pontchartrain (top), Breton (middle) and 278 Barataria (lower) for the current conditions simulation. **Figure 5.** April to September salinity region variation in Pontchartrain (**top**), Breton (**middle**) and Barataria (**lower**) for the current conditions simulation.

The FWOA run maintains the same connections between the basins and Mississippi River as the current conditions simulation. However, the 50-year simulation includes the effects of the relative sea-level rise and the changing wetland landscape within the basins [33]. Under FWOA, all basins show an increase over time in the proportion of the basin in the 2–16 ppt range, and both Pontchartrain and Breton show an increase in the area >16 ppt by year 50. The rate of eustatic sea-level rise accelerates over time in the simulation and is assumed to be the lead cause of that significant increase in salinities in the last decade.

The base run shows the same increase in salinity levels in Pontchartrain, Breton, and Barataria for April–September as found in FWOA with some individual years, fresher or more saline. The simulation including all the 2017 Master Plan diversions does not result in a substantially different salinity regime compared to the base run, mostly because both simulations include several diversion projects that divert large freshwater volumes to the basins. The operation plans are different between the simulations but not sufficiently to create major differences in the salinity regime. Table 2 shows the least amount of area >16 ppt at year 50 in this run compared to the other runs in Barataria, with the other basins being similar.

The annual snapshots shown in Table 2 mask the extensive interannual variation, and the differences in the response among the basins to events included in the 50-year boundary conditions. Of note, here is a severe drought (low rainfall) in the simulation in years 42–43. Figure 6 shows the effect of this drought in Pontchartrain in years 42–43 for the three simulations. For each simulation, these years show the greatest extent >16 ppt and the least extent <0.6 ppt.

Figure 6 shows a few major trends over time in Pontchartrain with substantial interannual variability, as might be expected where there are connections to riverine inflows. The last 15 years of the simulations show a generally smaller area of <0.6 ppt. In Pontchartrain, there appears to be an increase in the frequency of years with large areas >16 ppt in the later years which may be associated with the progressive influence of sea-level rise. This occurs in all three simulations and may be more related to the changes in freshwater from tributary inflows to the basin from inland streams, which remains the same for all simulations, rather than the river diversion inflows.

The changes in the Breton basin are shown in Figure 7. This basin has no inputs of freshwater apart from the precipitation and multiple relatively small controlled or uncontrolled connections to the Mississippi River (Table 1). The last 15 years of the simulations show a generally smaller area of <0.6 ppt in all simulations shown in Figure 7, and in FWOA there is also a decrease in the area <2 ppt. The FWOA also shows some incursion of water >16 ppt in isolated years in the last decade in contrast to earlier years of the simulation.

With the mid-basin diversion, the area in Breton <0.6 ppt increases compared to FWOA, showing that the large diversion has a greater impact on salinity than the existing connections (see next section). However, there are more years in the base and MP17 runs compared to FWOA, especially early in the simulations, when Breton includes areas >16 ppt. This seems counterintuitive when more freshwater is entering the basin through the diversion. However, as diversions into Barataria are also operating in these simulations, the flow toward the mouth of the Mississippi is diminished, reducing the discharge through distributaries in the Birdsfoot Delta, and increasing salinity in the lower parts of the Breton basin (Figure 3).

*Water* **2019**, *11*, x FOR PEER REVIEW 12 of 26

321 **Figure 6.** April to September salinity region variation in Pontchartrain for runs with varying levels of 322 freshwater inflow **Figure 6.** April to September salinity region variation in Pontchartrain for runs with varying levels of freshwater inflow

Barataria shows a very different pattern from the other basins, especially under FWOA (Figure 8). The FWOA run maintains the same connections between the basins and Mississippi River as current conditions. The 50-year simulation includes the effects of the relative sea-level rise and the changing wetland landscape and these effects are most dramatically shown in FWOA. There is a major increase in the extent of the basin with 16–35ppt salinity after year 30 (Figure 8). This period also shows a major decrease in the area <0.6 ppt. In the last few years of the simulation, only a small area of the basin experiences salinities <2 ppt.

*Water* **2019**, *11*, x FOR PEER REVIEW 13 of 26

324 **Figure 7.** April to September salinity region variations in Breton for runs with varying levels of 325 freshwater inflow. **Figure 7.** April to September salinity region variations in Breton for runs with varying levels of freshwater inflow.

326 Barataria shows a very different pattern from the other basins, especially under FWOA (Figure 8). 327 The FWOA run maintains the same connections between the basins and Mississippi River as current 328 conditions. The 50-year simulation includes the effects of the relative sea-level rise and the changing 329 wetland landscape and these effects are most dramatically shown in FWOA. There is a major increase 330 in the extent of the basin with 16–35ppt salinity after year 30 (Figure 8). This period also shows a 331 major decrease in the area <0.6 ppt. In the last few years of the simulation, only a small area of the 332 basin experiences salinities <2 ppt. With the diversions operating, Barataria does not experience the dramatic increase in the areas >16 ppt found in FWOA. Figure 8 shows remarkably little change over time. There is an increase for most of the 50-year simulation in the extent <0.6 ppt in Barataria in the base run but during the drought period, a large area <2 ppt remains. The MP17 simulation includes an additional diversion into the Barataria basin at Ama (Table 1, Figure 1) which appears to have little effect on the salinity distribution for most of the simulation, but reduces the effect of the drought. Even with this diversion, there are still conditions in specific years later in the simulation, that dramatically reduce the area of the basin <0.6 ppt.

*Water* **2019**, *11*, x FOR PEER REVIEW 14 of 26

334 **Figure 8.** April to September salinity region variations in Barataria for runs with varying levels of 335 freshwater inflow. **Figure 8.** April to September salinity region variations in Barataria for runs with varying levels of freshwater inflow.

#### 336 With the diversions operating, Barataria does not experience the dramatic increase in the areas *3.3. Modified Exchange between the Mississippi River and Breton Sound*

350 Master Plan FWOA and the base run (Table 1).

337 >16 ppt found in FWOA. Figure 8 shows remarkably little change over time. There is an increase for 338 most of the 50-year simulation in the extent <0.6 ppt in Barataria in the base run but during the 339 drought period, a large area <2 ppt remains. The MP17 simulation includes an additional diversion 340 into the Barataria basin at Ama (Table 1, Figure 1) which appears to have little effect on the salinity 341 distribution for most of the simulation, but reduces the effect of the drought. Even with this diversion, 342 there are still conditions in specific years later in the simulation, that dramatically reduce the area of The currently active exchange points (distributaries) between the Mississippi River and Breton Sound (e.g., Caernarvon Freshwater Diversion, Mardi Gras Pass, Bohemia Spillway, Fort St. Philip, etc.) were modeled as inactive with the flow being maintained in the main channel of the Mississippi River and eventually reaching the Birdsfoot Delta. All other simulation conditions (including diversion operations) and environmental assumptions were as represented in the 2017 Coastal Master Plan FWOA and the base run (Table 1).

343 the basin <0.6 ppt. 344 *3.3. Modified Exchange Between the Mississippi River and Breton Sound*  345 The currently active exchange points (distributaries) between the Mississippi River and Breton 346 Sound (e.g., Caernarvon Freshwater Diversion, Mardi Gras Pass, Bohemia Spillway, Fort St. Philip, 347 etc.) were modeled as inactive with the flow being maintained in the main channel of the Mississippi Figure 9 shows the salinity difference maps for the run with modified exchange compared to the base run at year 25 and 50, respectively. The simulation results show a salinity increase that extends east into the Lower Pontchartrain basin due to the reduced freshwater into Breton. The April-September salinity distribution plots in Figure 10 show the salinity gradient in Breton due to the disconnection from the Mississippi River, which shows a greater extent of saline bins compared to the base run (Figure 7) while little change is observed in Barataria or Pontchartrain.

348 River and eventually reaching the Birdsfoot Delta. All other simulation conditions (including 349 diversion operations) and environmental assumptions were as represented in the 2017 Coastal

357

362

362

365

365

351 Figure 9 shows the salinity difference maps for the run with modified exchange compared to the 352 base run at year 25 and 50, respectively. The simulation results show a salinity increase that extends 353 east into the Lower Pontchartrain basin due to the reduced freshwater into Breton. The April-354 September salinity distribution plots in Figure 10 show the salinity gradient in Breton due to the

351 Figure 9 shows the salinity difference maps for the run with modified exchange compared to the 352 base run at year 25 and 50, respectively. The simulation results show a salinity increase that extends 353 east into the Lower Pontchartrain basin due to the reduced freshwater into Breton. The April-354 September salinity distribution plots in Figure 10 show the salinity gradient in Breton due to the

*Water* **2019**, *11*, x FOR PEER REVIEW 15 of 26

356 the base run (Figure 7) while little change is observed in Barataria or Pontchartrain.

356 the base run (Figure 7) while little change is observed in Barataria or Pontchartrain.

358 **Figure 9.** The salinity difference when exchanges with Breton Sound through distributary channels 359 downstream of the Mid-Basin diversions are reduced (modified exchange run compared to the base 360 run with both Mid-Basin diversions) at year 25 (**left**) and year 50 (**right**). Eustatic sea level rise by year 361 25 was 0.33 m and 0.63 m by year 50. **Figure 9.** The salinity difference when exchanges with Breton Sound through distributary channels downstream of the Mid-Basin diversions are reduced (modified exchange run compared to the base run with both Mid-Basin diversions) at year 25 (**left**) and year 50 (**right**). Eustatic sea level rise by year 25 was 0.33 m and 0.63 m by year 50. 358 **Figure 9.** The salinity difference when exchanges with Breton Sound through distributary channels 359 downstream of the Mid-Basin diversions are reduced (modified exchange run compared to the base 360 run with both Mid-Basin diversions) at year 25 (**left**) and year 50 (**right**). Eustatic sea level rise by year 361 25 was 0.33 m and 0.63 m by year 50.

363 **Figure 10.** April to September salinity region variation in Breton Sound eco-region with minimal 363 **Figure 10.** April to September salinity region variation in Breton Sound eco-region with minimal 364 exchange between the Mississippi River and Breton Sound. **Figure 10.** April to September salinity region variation in Breton Sound eco-region with minimal exchange between the Mississippi River and Breton Sound.

364 exchange between the Mississippi River and Breton Sound.

The vegetation coverage for the base run and minimal exchange for year 25 and 50 are shown in Figure 11. The vegetation in lower Breton Sound and along the east bank side of the Mississippi River converts towards a more saline type due to the freshwater reduction with the modified exchange as the direct freshwater input through Mardi Gras Pass and Fort St. Philip (Figure 3) are reduced. While less vegetated land is shown in year 50 as the effects of the sea-level rise and subsidence have resulted in land loss, a similar effect is apparent. At year 50, with 0.63 m of additional sea-level rise, the diversion alone (Modified Exchange in Figure 11) maintains fresh marsh in the upper part of the Breton basin but is insufficient to prevent saline marsh along the east bank of the river.

378

366 The vegetation coverage for the base run and minimal exchange for year 25 and 50 are shown in 367 Figure 11. The vegetation in lower Breton Sound and along the east bank side of the Mississippi River 368 converts towards a more saline type due to the freshwater reduction with the modified exchange as 369 the direct freshwater input through Mardi Gras Pass and Fort St. Philip (Figure 3) are reduced. While 370 less vegetated land is shown in year 50 as the effects of the sea-level rise and subsidence have resulted 371 in land loss, a similar effect is apparent. At year 50, with 0.63 m of additional sea-level rise, the

373 Breton basin but is insufficient to prevent saline marsh along the east bank of the river.

375 **Figure 11.** Vegetation coverage for base run with both of the Mid-Basin diversions (left) and the 376 modified exchange run with both diversions active but other connections to Breton Sound reduced in 377 flow capacity (right). Eustatic sea level rise by year 25 (top) was 0.33 m and 0.63 m by year 50 (bottom). **Figure 11.** Vegetation coverage for base run with both of the Mid-Basin diversions (**left**) and the modified exchange run with both diversions active but other connections to Breton Sound reduced in flow capacity (**right**). Eustatic sea level rise by year 25 (**top**) was 0.33 m and 0.63 m by year 50 (**bottom**).

There is also a marked difference in the extent of the brackish marsh zone between year 25 and year 50. At year 25, this is a very narrow band in both simulations. The freshwater from the diversion and saline water from the Gulf squeeze out the brackish zone. However, as the sea-level rise progresses, saltwater penetrates further, mixing with the fresh conditions maintained by the diversion to produce a more extensive zone of brackish marsh. It should be noted that in addition to the vegetation shifts over time and under different freshwater exchange scenarios, there is a marked decrease in the amount of wetland area throughout lower Barataria and lower Breton basins under this assumed relative sea level rise scenario. This conversion of wetland area to open water is a consequence of persistent inundation and salinity penetration due to the combined effects of eustatic sea level rise and subsidence [33].

#### *3.4. Operational Regime and Sediment Delivery*

The purpose of most planned river diversions in Louisiana is to provide sediment for rebuilding of coastal wetlands and their maintenance under a future sea-level rise [45]. Their location and design are often optimized based on the potential for delivery of sand from the river through the diversion structure [46,47]. Table 3 presents the cumulative suspended sand load (not including bedload or fines in suspension) diverted for the simulations described above. These observations also point to the need for careful management of the operation of proposed restoration strategies (with dynamic components).

These simulations reflect the combined influence of different combinations of diversions with and without various existing connections between the river and estuary, over 50 years. Table 3 shows greater potential for sand delivery into Breton, which may be the result of the Mid-Breton diversion being upstream of the Mid-Barataria diversion, reducing the flow available for distribution when they are both operating simultaneously.

Varying the operational regimes shown in Table 1 would modify these sand loads. However, as fine sediments are important for coastal marsh maintenance in the face of sea-level rise [45], an additional simulation was conducted to assess the effects of an operational regime and the river flow on fine sediment delivery. This simulation used the Central Wetlands diversion as an example (see Figure 1 for location). The single-year runs were conducted to examine the different configurations and operational rules. All other simulation conditions (including diversion operations—see Table 1) and environmental assumptions were as represented in the 2017 Coastal Master Plan FWOA and the base run.

Three Central Wetlands diversion configurations were evaluated:

The diversion discharge for these three different operational configurations in the selected low, medium, and high river years (see Methods) are shown in Figure S1 in Supplementary Materials: ICM Simulation Settings. The three diversion operations were evaluated in these pre-selected river flow conditions and the sand concentration in the diverted water was set to zero in order to evaluate sediment accumulation from fine particles.

The fine sediment accumulation in the receiving basin for the Central Wetlands diversion for different diversion configurations in different river years are shown in Figure 12.

428 sediment accumulation from fine particles.

417 3. The diversion opened when river flow was above 12,742 cms with a discharge of 418 0.83% of the Mississippi River flow. This diversion to river ratio was determined as 419 the sum of the total diversion discharge if it operated with a constant 141 cms 420 flowrate in 50-years and divided by the total river flow when the diversion is opened 421 (river flow above 12,742 cms) in 50 years. The idea is to operate the diversion only 422 during the high river condition while retaining the same total annual discharge

424 The diversion discharge for these three different operational configurations in the selected low, 425 medium, and high river years (see Methods) are shown in Figure S1 in Supplementary Materials*:*  426 ICM Simulation Settings. The three diversion operations were evaluated in these pre-selected river 427 flow conditions and the sand concentration in the diverted water was set to zero in order to evaluate

423 volume as with a constant flowrate (configuration 1).

432 **Figure 12.** Sediment accumulation in the receiving compartment for the three Central Wetlands 433 diversion configurations: A constant year-round flowrate of 141 cms as modeled for the 2017 Master 434 Plan (**1**), a flowrate activated when the river flow was above a given threshold with a maximum 435 diverted rate of 141 cms (**2**), and a diversion such that the annual volume diverted was equal to the 436 annual volume diverted in configuration 1, but only during periods of peak flow in the river (**3**). **Figure 12.** Sediment accumulation in the receiving compartment for the three Central Wetlands diversion configurations: A constant year-round flowrate of 141 cms as modeled for the 2017 Master Plan (**1**), a flowrate activated when the river flow was above a given threshold with a maximum diverted rate of 141 cms (**2**), and a diversion such that the annual volume diverted was equal to the annual volume diverted in configuration 1, but only during periods of peak flow in the river (**3**).

The receiving compartment is a 17 km2 tract of marsh (16% of which is open water) on the east bank of the Mississippi River. The Central Wetlands tract is bounded on all sides by hydrologic barriers (flood walls, levees, roadways) and currently receives freshwater inflows from several pump station discharge lines which drain the eastern portions of New Orleans and neighboring municipalities. The sediment accumulation for diversion configuration 3 is higher compared to diversion configuration 1 and 2 in all simulated river conditions. Diversion configuration 3 performs better in low and medium river years when compared to the high river year.

The fifty-year simulations with diversion operational configuration 1 and 3 were conducted. Salinity in the immediate receiving area is highly influenced by the different operational configurations. As illustrated in Figure 13, operational configuration 1 resulted in a sustained low salinity in the diversion receiving area, while in operational configuration 3, salinity is only slightly different from the no diversion condition.

448 different from the no diversion condition.

The receiving compartment is a 17 km2 437 tract of marsh (16% of which is open water) on the east 438 bank of the Mississippi River. The Central Wetlands tract is bounded on all sides by hydrologic 439 barriers (flood walls, levees, roadways) and currently receives freshwater inflows from several pump 440 station discharge lines which drain the eastern portions of New Orleans and neighboring 441 municipalities. The sediment accumulation for diversion configuration 3 is higher compared to 442 diversion configuration 1 and 2 in all simulated river conditions. Diversion configuration 3 performs

444 The fifty-year simulations with diversion operational configuration 1 and 3 were conducted. 445 Salinity in the immediate receiving area is highly influenced by the different operational 446 configurations. As illustrated in Figure 13, operational configuration 1 resulted in a sustained low

443 better in low and medium river years when compared to the high river year.

450 **Figure 13.** Salinity of the diversion receiving compartment for different operational configurations. **Figure 13.** Salinity of the diversion receiving compartment for different operational configurations.

#### 451 **4. Discussion 4. Discussion**

449

#### 452 *4.1. Future Effects of Relative Sea-Level Rise 4.1. Future E*ff*ects of Relative Sea-Level Rise*

453 The effect of 50 years of sea-level rise on estuarine environments is often focused on the fate of 454 the coastal wetlands and the impact of their loss of habitat for fish and wildlife [48–50]. The 455 simulations described here show that the effects of the sea-level rise on these estuarine basins can be 456 largely offset by more extensive use of river resources. Barataria basin (Figure 8) shows the most 457 variation in salinity during the FWOA with a marked increase in later years of the simulation. It has 458 a few existing connections to the Mississippi River, other than the Davis Pond freshwater diversion 459 (Table 1), and no other freshwater inputs other than precipitation, stormwater runoff from local 460 communities. Under FWOA, the fresher parts of the estuarine gradient, i.e., <2 ppt, are greatly 461 reduced and almost eliminated during drought years with low rainfall (Figure 8, Table 2). The effect of 50 years of sea-level rise on estuarine environments is often focused on the fate of the coastal wetlands and the impact of their loss of habitat for fish and wildlife [48–50]. The simulations described here show that the effects of the sea-level rise on these estuarine basins can be largely offset by more extensive use of river resources. Barataria basin (Figure 8) shows the most variation in salinity during the FWOA with a marked increase in later years of the simulation. It has a few existing connections to the Mississippi River, other than the Davis Pond freshwater diversion (Table 1), and no other freshwater inputs other than precipitation, stormwater runoff from local communities. Under FWOA, the fresher parts of the estuarine gradient, i.e., <2 ppt, are greatly reduced and almost eliminated during drought years with low rainfall (Figure 8, Table 2).

462 Barataria also has a relatively intact chain of barrier islands, separated by tidal inlets, at the 463 Gulfward margin (Figure 3). Barrier islands have been shown to modulate salinity in Louisiana 464 estuaries [51]. The simulations described here do not discuss barrier island dynamics over time. 465 Further, how the barrier shoreline is maintained in the future, as it also faces impacts from the sea-466 level rise, will be important for future salinities and the relative influence of diversion in the Barataria Barataria also has a relatively intact chain of barrier islands, separated by tidal inlets, at the Gulfward margin (Figure 3). Barrier islands have been shown to modulate salinity in Louisiana estuaries [51]. The simulations described here do not discuss barrier island dynamics over time. Further, how the barrier shoreline is maintained in the future, as it also faces impacts from the sea-level rise, will be important for future salinities and the relative influence of diversion in the Barataria basin.

467 basin. 468 As this modeling has shown, the interaction of the future sea-level rise with other climatic factors 469 is more complex. The inclusion of the drought years in year 42–43 of the simulations showed that 470 even in Pontchartrain where sea-level rise did not cause a major progressive increase in salinity > 16 471 ppt, the effect of a higher sea-level and reduced precipitation caused a temporary increase (Figure 6). 472 The potential for adaptive operations for diversions to mitigate such events is discussed below. 473 However, one other key influence on the issues discussed here is the effect of future climate change 474 on Mississippi River discharge. The simulations included here are all based on the same 50-year 475 record of the Mississippi River flow in order that the effects of diversions on the change in the As this modeling has shown, the interaction of the future sea-level rise with other climatic factors is more complex. The inclusion of the drought years in year 42–43 of the simulations showed that even in Pontchartrain where sea-level rise did not cause a major progressive increase in salinity > 16 ppt, the effect of a higher sea-level and reduced precipitation caused a temporary increase (Figure 6). The potential for adaptive operations for diversions to mitigate such events is discussed below. However, one other key influence on the issues discussed here is the effect of future climate change on Mississippi River discharge. The simulations included here are all based on the same 50-year record of the Mississippi River flow in order that the effects of diversions on the change in the estuaries due to the sea-level rise can be specifically explored. Studies have shown that trends in flooding on the Mississippi River are influenced by climate variability including the El Niño–Southern Oscillation and the Atlantic Multidecadal Oscillation [52], although predicting future river flow is challenging in large river systems [53].

#### *4.2. E*ff*ects of River Diversions on Estuarine Conditions*

This study has shown the influence of different magnitudes of river inflow on salinity gradients within the estuaries. The comparison of the future with only minor connections between the river and estuarine basin (FWOA) with two major diversions (base run) and then with additional diversions (MP17) for salinity distribution shows differences among the basins (Figures 5–7). For Pontchartrain (Figure 6) which does not directly receive additional inflows under the base run, the difference is the greatest in the last decade in the extent of the basin <0.6 ppt. The incremental addition of diversions increases the area <0.6 ppt. The area <0.6 ppt, however, is smaller in extent than experienced in earlier decades, reflecting the progressive influence of the sea-level rise on the system. Even though the diversions are not directly influencing Pontchartrain basin, the extent of the freshwater influence from large diversions into Breton basin extends into lower Pontchartrain (Figure 9) perhaps ameliorating the penetration of higher salinities from the Gulf. There is much less difference in Pontchartrain among the higher salinity levels.

In contrast, in FWOA Breton (Figure 7) has fluctuating areas of <0.6 ppt, likely due to fluctuations in the river flow (Figure 2) and the many small existing connections between the river and the basin (Table 1). Breton Sound basin both experienced fewer periods of the freshest salinity range, with three of the years having no portion of the basin with a mean salinity <0.6 ppt. While in four other years, the basin had no portion with the mean salinities greater than 16 ppt. These freshening events occur during the wetter-than-average years of the flow in the Mississippi River. This was the case in five of the last six years in this simulation, 2008 through 2013 with the exception of 2012 (Figure 2). All five of these years were years with higher-than-average flowrates in the Mississippi River. The sixth year, 2012, was well below the average for the river flow. The sensitivity to the Mississippi River flow is due to the lack of any other pluvial source into Breton Sound. Additionally, the most Gulfward portions of the Breton Sound Basin are somewhat buffered from the open waters of the Gulf by both the distal barrier Chandeleur Islands enhancing the influence of Mississippi River fluctuations on salinities in the lower Breton basin [54]. This results in a basin that has less overall range to salinity in most years. However, the current conditions simulations ended in 2013, and did not include the flowrates for Mardi Gras Pass which initially formed during the 2011 flood. While one flow rating curve was used in this modeling (Table 1), the observed flows leaving this diversion and other connections to Breton had a noticeable impact on salinities in Breton, which was also seen in future simulations with these connections turned off (Figure 9).

Thus, the addition of the Mid-Breton diversion in the base run increases the freshwater area <0.6 ppt. However, Figure 7 also shows little change in the salinity zones between 2 and 16 ppt. This is shown in the vegetation patterns for year 25 (Figure 11) and the greater brackish area shown in Figure 11 at year 50 comes at the expense of the freshwater areas, the latter of which is limited in all three simulations in the last decade (Figure 7). Even the addition of the Lower Breton diversion (Table 1, Figure 1) in the MP17 simulation does not substantially change the basin scale conditions in the last decade, although there is likely localized freshening.

#### *4.3. E*ff*ect on Estuarine Fisheries*

Estuaries are an important habitat for the Gulf of Mexico's commercial and recreational fisheries with 97% (by weight) depending on estuaries for some part of their life cycle [55]. The importance of the freshwater inflow to estuarine ecosystems is widely recognized [56–58]. Often the focus of concern is the upstream impoundment of streams for water supply or flood management [59,60] but in the estuarine basins discussed here, which have been essentially isolated from riverine inflows for more than a century due to river management in the delta plain, the key issue is where, when, and how much inflow is required to meet certain goals and avoid undesirable outcomes for other parts of the ecosystem [61].

The existing diversions have changed the aquatic ecology of their receiving area. Studies of Davis Pond and Caervarvon have examined the effects of diversion influxes of freshwater on the ecology of the current system [62,63], including the effects on important commercial and recreational species, such as spotted seatrout *Cynoscion nebulosus*, red drum *Sciaenops ocellatus*, Gulf menhaden *Brevoortia patronus*, blue crab *Callinectes sapidus*, and penaid shrimp. These effects vary with discharge and are often localized [64]. Modeling has also been used as a tool in such evaluations. For example, the simulations of fish response to the operation of the Caernarvon freshwater diversion showed movement down the estuary of Bay Anchovy *Anchoa mitchilli*, which agrees with field observations [63]. The long-term simulations of fish communities in Barataria, Breton, and Pontchartrain basins in response to the

sea-level rise and the inflow from four large river diversions (Mid-Barataria, Mid-Breton, Lower Breton in Table 1 and an additional diversion into Lower Barataria) showed local shifts in the distribution of six species of ecological and economic importance including the Eastern oyster *Crassostrea virginica* [65]. Some changes were the direct effects of decreased salinity while other effects were related to reduced Chl *a* due to increased turbidity, a trend also identified in the field and literature-based studies [66]. The simulations also showed that local biomass reductions were mostly the result of redistribution, as there was minimal impact on the total biomass compared to a future without action.

The ecological effect of the larger scale riverine inflow is likely complex. An examination of the salinity sensitivity of 44 estuarine fish species, including many of economic importance, in the Gulf of Mexico estuaries found many species showing low or moderate sensitivity to the salinity change [67]. However, freshwater fishes currently in the extensive areas of the basin >0.6 ppt, may be impacted by increased salinity through the impacts on physiological costs (e.g., osmoregulation) which may be offset through increased productivity and greater food availability [68]. Further, experimental studies using coastal marsh fish which are found across the salinity gradient showed that descendants of individuals from brackish populations had better survival when exposed to salinity than fish with no historical exposure to salinity, suggesting genetic adaptation [69].

#### *4.4. Lessons for Future Management*

It is not a new concept that large-scale river diversions into the Mississippi River estuary could utilize an operational regime which mimics the natural fluctuations of riverine flood events [70]. Crevasses and small uncontrolled diversions show such a flow regime [71,72] in contrast to the near continuous flow of a delta distributary channel. This analysis examined three distinct operational regimes to explore how a diversion could be operated to mimic a seasonal high flow and sediment pulse that is typically seen in the springtime Mississippi River flows [10,45,70]. Capturing an equivalent amount of diverted water, the flood-wave specific operations are able to capture and deliver a greater load of sediment from the river and distribute it to the estuary (Figure 12). This regime has the largest impact during river years with below-average flowrates when the total annual load must be scheduled during shorter periods (Figure S1 in Supplementary Material). This operation regime, while simple enough to implement in a model with assigned boundary conditions, the real-world implementation would be challenging, given the inability for accurate seasonal forecasts of the flow in the Mississippi River. In practice, such a regime would likely operate as soon as the river is high enough (typically in the spring) and continue until the annual sediment/freshwater load mass/volume were met.

An additional benefit from flood-mimicking operational regimes is that the concentrations of suspended sediments within the Mississippi River waters (both fines and sands) vary in relation to flood events and annual cycles [73]. The current flow-to-sediment rating curves for the lowermost Mississippi River show a hysteresis effect where the rising limb of the annual river flood has a flushing effect and carries a greater suspended sediment load than an equivalent flowrate on the falling limb of the flood hydrograph [10]. The scheduling diversion operations for the arrival of the floodwave allow for the most efficient delivery of suspended sediments from the Mississippi River into the estuary (Figure 12).

The simulations presented in this report provided insights on how the system responds to a variety of operation plans. Some of the responses are not intuitive, emphasizing the value of using such holistic landscape models that capture (to the extent possible) feedbacks and interactions among various processes. The FWOA simulation shows a dramatic change in salinity in parts of the system due to the sea-level rise. The rate of the eustatic sea-level rise accelerates over time in the simulation and is assumed to be the lead cause of that significant increase in salinities in the last decade. These changes have been shown to have dramatic effects on vegetation [32] and, in combination with the changes in temperature, can dramatically influence the distribution of aquatic species [34].

One important aspect of using the resources of the Mississippi River to manage the coastal ecosystem in the face of the sea-level rise is that it allows the management of conditions to achieve a variety of goals, e.g., land building or positioning suitable salinity conditions for fishery benefits, as well as the ability to adjust to changing conditions. The simulations here all apply a fixed set of rules over time. In reality, the diversions are likely to be managed in response to estuarine (and river) conditions. The effects of the drought event in the last decade of the simulation were in some way ameliorated by the various diversion configurations (Figures 5–7) even with operational rules which were not designed to specifically address such situations. As the implementation of diversions proceeds, near-term forecasting models of estuarine conditions and adaptable operational rules for diversion flows may be useful to further protect against massive interannual fluctuations on estuarine conditions. However, this study has shown that planning how many diversions are needed and at which locations within the estuary must focus on the long-term in order that the effects of the sea-level rise can be integrated into decision making.

#### **5. Conclusions**

The simulations described here show, for the estuarine basins of the Mississippi Delta Plain, that the effects of the sea-level rise decades into the future are potentially dramatic, but the effects on the estuarine salinity gradient varies according to the landscape setting. The analysis demonstrates that the ecological character of the future system is dependent on both the sea-level rise and the management of river resources. The pronounced incursion of higher salinity water into Barataria in later years of the simulation associated with the relative sea-level rise was largely offset by the freshwater inputs from the planned Mid-Barataria diversion. While sediment delivery is frequently the focus of future river diversions in the Mississippi River Delta, this analysis has shown how large-scale river diversions also affect spatiotemporal patterns in the salinity gradient throughout the estuarine basins, combatting the effects of a future sea-level rise by maintaining present-day salinity variations throughout the basin. The spatial extent of these salinity gradients varies annually due to the variability in riverine freshwater flows reaching the delta and are also susceptible to the number (and operational regime) of engineered diversions and natural distributary channels into the estuary.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/11/10/2028/s1: *ICM Overview* and *ICM Simulation Settings*. Please contact the corresponding author for access to model source code repositories on GitHub.

**Author Contributions:** The co-authors contributed to the research and preparation of this manuscript in the following ways: conceptualization, E.M., D.R., N.P.S., A.R., and E.D.W.; methodology, E.D.W., Y.W., E.M., and D.R.; software, E.D.W., Y.W., E.M., and others listed in acknowledgements; investigation, E.D.W., Y.W., D.R., and E.M.; writing—original draft preparation, E.D.W., Y.W., D.R., E.M. and A.R.; writing—review and editing, D.R., A.R., E.D.W., N.P.S. and E.M.; visualization, E.D.W., Y.W. and A.R.; supervision, D.R., E.M., A.R. and N.P.S.; project administration, N.P.S. and A.R.; funding acquisition, N.P.S.

**Funding:** This analysis was funded by Environmental Defense Fund (EDF) (grant number 1048-000005-10650; #7557790/#7930492) through support from the Gilder Foundation.

**Acknowledgments:** The Integrated Compartment Model development and calibration discussed in this work was initially funded by the Louisiana Coastal Protection and Restoration Authority; the authors are indebted to the many researchers and engineers who participated in the 2017 Coastal Master Plan analyses and assisted in those efforts. We would like to acknowledge fellow model developers for the various subroutines of the ICM including (among many others): Alex McCorquodale, Jenneke Visser, Scott Duke-Sylvester, Brady Couvillion and Zhifei Dong.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**

1. Nicholls, R.J.; Wong, P.P.; Burkett, V.; Codignotto, J.; Hay, J.; McLean, R.; Ragoonaden, S.; Woodroffe, C.D.; Abuodha, P.A.O.; Arblaster, J.; et al. Coastal systems and low-lying areas. In *Climate Change 2007: Impacts, Adaptation and Vulnerability*; Parry, M.L., Canziani, O.F., Palutikof, J.P., van der Linden, P.J., Hanson, C.E., Eds.; Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK, 2007; pp. 315–356.

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Natural Variability and Vertical Land Motion Contributions in the Mediterranean Sea-Level Records over the Last Two Centuries and Projections for 2100**

#### **Antonio Vecchio 1,2 , Marco Anzidei 3,\*, Enrico Serpelloni 4 and Fabio Florindo 3**

Received: 14 June 2019; Accepted: 5 July 2019; Published: 16 July 2019

**Abstract:** We analyzed a set of geodetic data to investigate the contribution of local factors, namely the sea level natural variability (SLNV) and the vertical land motion (VLM), to the sea-level trend. The SLNV is analyzed through the Empirical Mode Decomposition (EMD) on tidal data (>60 years of recordings) and results are used to evaluate its effects on sea levels. The VLM is measured at a set of continuous GPS (cGPS) stations (>5 years of recordings), located nearby the tide gauges. By combining VLM and SLNV with IPCC-AR5 regional projections of climatic data (Representative Concentration Pathways (RCP) 2.6 and 8.5), we provide relative sea-level rise projections by 2100. Results show that the combined effects of SLNV and VLM are not negligible, contributing between 15% and 65% to the sea-level variability. Expected sea levels for 2100 in the RCP8.5 scenario are between 475 ± 203 (Bakar) and 818 ± 250 mm (Venice). In the Venice Lagoon, the mean land subsidence at 3.3 ± 0.85 mm a−1 (locally up to 8.45 ± 1.69 mm a−1 ) is driving the local sea-level rise acceleration.

**Keywords:** sea-level rise; Mediterranean Sea; tide gauges; natural variability; vertical land motion

#### **1. Introduction**

In situ and remote sensing data show that the rate of global sea-level rise over the past two centuries has increased at faster rates than in the last two or three millennia [1–4]. Satellite radar altimetry and tide gauge data indicate a rise of the global mean sea-level rate at 1.7 mm a−1 and 3.2 mm a−1 for the 20th century and the last two decades, respectively [5–7]. The Mediterranean Sea has been rising at the mean rate of ~1.8 mm a−1 in the last two-three centuries [8–12]. At global scales, the sea-level change is described as the sum of eustatic, glacio-hydro-isostatic, and land-hydrology components [13], which are independently modeled to provide future sea-level predictions. Local sea-level change can significantly differ from mean global sea-level rise because local factors become relevant and their contribution adds up to the global components mentioned above [14]. Oceanographic and climatic effects, glacial isostatic adjustment (GIA), and continuous vertical land motion, especially coastal subsidence due to tectonics (volcanic included) or anthropic origins, are some of the factors that mainly affect sea-level changes at local scales [15]. Therefore, depending on the geographical location they cause spatial variability in the sea-level change pattern. While local effects have little influence at global scale, since they are smoothed out by averaging operations, they become critical in adaptation planning and risk management, where localized assessments are needed.

At time scales of years or decades, oceanographic and climatic factors drive non-stationary sea-level fluctuations, the so-called sea-level natural variability, whose assessment is important for future sea-level projections since their contribution is not negligible [16–26], (http://www.ipcc.ch). For example, in the Indian Ocean the effects of the natural variability, mainly driven by the decadal El Niño–Southern Oscillation (ENSO), Indian Ocean Dipole (IOD) and monsoon, dominate the external forcing contributing up 85% of the observed trend [27]. In the Mediterranean basin, the glacio-isostatic contribution, acting from the Last Glacial Maximum, has a predictable pattern that has been estimated and compared with direct observational data in deforming zones [4,9,15]. The vertical land motion due to natural (tectonics, isostasy and ground compaction) [4,9,13], and/or anthropogenic factors (fluid extraction) [28,29], produces changes in the shape of the basins and changes in the height of the land with respect to the sea (see e.g., [9,13] and references therein). Although this effect does not directly change the quantity of water and its physical properties it results in a net increase (in case of subsidence) or decrease (in case of uplift) of the sea level, when evaluated in restricted geographical areas (as for example the punctual measurements by tide gauges), since the land is moving. If the vertical land motion operates over time it locally results in a net sea-level rate that adds up to the rate induced by GIA, thus affecting the local sea-level change. Tide gauge records are strongly sensitive to the vertical land motion since these instruments are tied to land. Although the role and rates of the main contributors to the sea-level rise are still debated [30], it is a matter of fact that most of the coasts of the world are undergoing the effects of the continuing sea-level rise. The vertical motion due to tectonic effects is strongly variable in the Mediterranean region (as shown in [9] and references therein) and in some areas it can strongly enhance the effect of the sea-level rise. In particular, along subsiding coasts, the rising sea is expected to flood in advance, with respect to stable zones, over the coming decades [1,15,31].

In this framework, the Northern Mediterranean region, being characterized by coastal plains with low elevation (even placed below sea level as some areas of the Adriatic Sea) is therefore prone to marine flooding [9,32–34], (www.savemedcoasts.eu), and it represents an interesting case study. Several historical maritime densely-inhabited cities (i.e., Venice and Trieste) and valuable environmental and industrial areas (i.e., Genoa and Marseille) are located in this area. Therefore, detailed studies of local sea-level variations and future projections are relevant to support marine flooding scenarios for a cognizant coastal management.

This study aims to characterize the effect of local factors, such as the sea level natural variability (SLNV) and vertical land motion (VLM), on the sea-level changes at nine locations in the Northern Mediterranean Sea (Figure 1). The SLNV was inferred from tide gauge data through the Empirical Mode Decomposition (EMD) while VLM was measured by means of continuous GPS (cGPS) data. Our approach to estimate SLNV is comprised of two steps: (i) the EMD is applied to tide gauge data to characterize and isolate the natural variability; (ii) a simple theoretical model, based on EMD results, is developed to reproduce the observed SLNV. By combining the SLNV and VLM contributions with the local time-dependent trend based on regional IPCC- fifth Assessment Report (AR5) estimations, the expected sea levels for 2050 and 2100 at each tide gauge location were estimated.

**Figure 1.** Location of the nine analyzed tidal stations: Marseille (5.35◦–43.32◦ , 1885–2012), Trieste (13.76◦–45.57◦ , 1875–2012), Genova (9.92◦–44.41◦ , 1884–1997), Venice Punta Salute (12.32◦–45.39◦ , 1909–2000), Bakar (14.5◦–45.30◦ , 1930–2011), Dubrovnik (19.13◦–42.65◦ , 1956−2009), Rovinj (13.64◦–45.06◦ , 1955–2011), Split RT Marjana (16.39◦–43.51◦ , 1952–2011), Split Gradska Luka (16.44◦–43.51◦ , 1954–2011). Coordinates of the station and time span of recordings are indicated in parentheses.

#### **2. Materials and Methods: Data Set and Natural Variability**

In this study we analyzed the longest available monthly tide gauge records (>60 years of Revised Local Reference (RLR) data collected in the 1888–2008 time span) at nine tide gauge stations located in the Mediterranean Sea, belonging to the Permanent Service for Mean Sea Level (www.psmsl.org). The location of the stations, the time span of the records and the analyzed time series are shown in Figures 1 and 2. The choice of data sets with length exceeding 60 years aims to avoid the effect of inter-decadal variations affecting tidal recordings shorter than ~50 years [35,36]. In addition, to avoid the introduction of artificial signals in the records, interpolations on data gaps were not applied.

The identification of sea-level trends and large timescales variability from sea-level records is a debated argument and the best methods to discriminate between acceleration and intrinsic climate variability remains an open issue. Indeed, since in tide gauge data the long-term trend is embedded in fluctuations at many different time scales, from daily tides to the seasonal cycle as well as interannual and decadal variations, the filtering method used to remove these variations may impact the accuracy of the calculated trend. Different approaches to this problem have been proposed (for a detailed review see [37]). Among these approaches, the EMD has been used to characterize the intrinsic fluctuations and the long-term trend in sea-level data [38–45]. This technique has been developed to analyze non-stationary and non-linear data, such as tide gauge records, by decomposing the signal in modes that provide a description of the intrinsic timescales present in the time series. Being "adaptive"" the modes are empirically obtained from the original time series according to its local characteristics. The EMD reveals the proper components of the phenomenon under study better than other methods, such as Fourier analysis, whose modes with predefined functional form are often far from being eigenfunctions of the phenomenon at hand. This makes the EMD very powerful for identifying and isolating, for example, deformed wave profiles, associated with nonlinearities in the data, preserving their local properties. For the applications of this paper, two main reasons guided us in using the EMD. Firstly, the sea level natural variability is described by a few numbers of modes with a quite simple functional form (oscillating functions). This will be exploited in the following for building up a simple low dimensional mode for sea-level estimations. Secondly, the EMD does not require a predetermined

timescale to carry out the filtering operation (as it happens e.g., in moving averages commonly used as low-pass filter to recover the large timescale variability in sea-level records): this has an important rational basis, since in a nonstationary process the local timescale is unknown a priori.

**Figure 2.** Time series of the monthly mean tidal records from stations in Figure 1 (data from www.psmsl.org). Red lines show linear regressions on the data (see the text for details).

In the EMD framework the measured sea level *L*(*t*) is decomposed into a finite number *n* of band-limited oscillating components (called intrinsic mode functions, IMFs) *Lj* (*t*), whose functional shape is not fixed a priori, but obtained from the data, and a residue *rn*(*t*) providing the nonlinear trend [40]:

$$L\_{\cdot}(t) = \sum\_{j=0}^{n} l\_{j}(t) + r\_{n}(t). \tag{1}$$

The procedure of IMF extraction consists of three steps: (i) local maxima and minima of *L*(*t*) are identified; (ii) these are connected through a cubic spline identifying the lower and upper envelope; (iii) the signal *h*(*t*), the difference between *L*(*t*) and the mean of lower and upper envelopes, are calculated. The signal *h*(*t*) represents the first IMF, *Lj* (*t*), if it has exactly one zero crossing between two consecutive local extrema and its mean is zero. The first condition ensures a global narrowband requirement, namely periods which are too different are not mixed together into an IMF; the second requirement ensures getting oscillations about the zero level. Both these conditions ensure that *Lj* (*t*) is an oscillating function with time-dependent amplitude and phase, namely it can be expressed as *Lj* (*t*) = *Bj* (*t*)cos[φ*j* (*t*)]. If none of these conditions are fulfilled, then step (iii) is repeated by using *h*(*t*) as input data. This procedure is iterated until the resulting signal fulfills the IMF properties and it is repeated to find the next IMF by using *L*(*t*) − *Lj* (*t*) as new input data. When no more modes can be extracted, the difference between the original *L*(*t*) and the sum of all IMFs results in the residue *rn*(*t*). Technical details on the procedure for IMF extraction are reported in [41]. The EMD allows calculating a meaningful instantaneous phase for each IMF by using the Hilbert transform

$$L\_j^\*\left(t\right) = \pi^{-1} \times P \times \int\_{-\infty}^{\infty} \frac{L\_j\left(t'\right)}{t - t'} \times dt' \tag{2}$$

where P denotes the Cauchy principal value and *Lj*\* represents the complex conjugate of *L*(*t*) as φj (*t*) = arctan (*Lj*\*/*Lj* ). Finally, the instantaneous frequency, ω*j* (*t*), is given by *d*φ*j* (*t*)/*dt*. For each *Lj* (*t*) we can define a characteristic timescale τ*j* , as the time average <*2*π/ω*j* (*t*)>, that represents an average period of the IMF and provides an estimate of the timescale characterizing the EMD mode. We remark that τ*j* is not to be intended as the Fourier one. It just gives an indication of the timescale characterizing the EMD mode for which it is computed, although many modes with different average periods may contribute to the variability of the actual signal at a particular timescale. The IMF statistical significance with respect to white noise can be evaluated through a test based on the comparison between IMF amplitude from the real signal and from a white noise series [46] at different confidence levels.  Moreover, since a low number of modes are produced, this property can be exploited to build up low order theoretical models to describe the dynamics of the investigated systems.

#### **3. Results: The EMD Analysis and a Generalized Model for Long-Term Sea-Level Variations**

When the EMD is used to analyze the nine tide gauge time series, the number of IMFs, *n,* varies between 10 and 15, depending on the station. This is mainly related to the different lengths of the time series which vary between 137 and 53 years: for the longest time series more components of the natural variability are recovered.

Let *S* be the set of IMF with timescale >15 years (IMFs timescales τ*j* are shown in Table 1) and significant with respect to a white noise at the 90th significance level. The long-term variations of the measured sea level, *LLT*, can be defined as:

$$L\_{LT}\left(t\right) = \sum\_{j \in \mathcal{S}} L\_j(t) \, | \, + \, r\_n\left(t\right), \tag{3}$$

where IMFs in *S* and the residue are added together.

**Table 1.** Characteristic periods (years) of the significant intrinsic mode functions (IMFs) with a characteristic timescale >15 years, obtained by the Empirical Mode Decomposition (EMD) for the analyzed data set. Errors are estimated as the standard deviation.

The choice of a 15-year threshold follows the approach used in other studies to remove the short-term variability from sea-level data [1,18]. We remark that, as pointed out by [45], the results of single IMFs from tide gauge data need of a careful interpretation. However, a discussion about the physical meaning of the single IMF is out of the scope of this paper since we use the EMD as a low-pass filter by adding up several IMFs. In this case, the EMD has shown to be a useful tool to reduce the impact of high-frequency variability and noise [45].

The contribution *LLT*, for the analyzed tide gauges, is shown in Figure 3. All records show a long-term variability in a well-defined range of timescales of about 20–30 and, for the longer records, >45 years, in the typical ranges where internal modes of natural variability are detected [47]. Observational datasets and simulations of the typical scales of Atlantic internal modes of natural variability show variations at decadal timescales that are linked to the ocean-atmosphere coupling [47–50]. In particular, the 20–30 year variability has been related to the Atlantic Meridional Overturning Circulation (MOC), while the >45 year fluctuations are plausibly due to salinity and matter exchange processes between the Atlantic and Arctic Ocean causing the Atlantic Multidecadal Oscillation (AMO) pattern. Both patterns are quasi-persistent and present up to 8000 years ago [51,52]. The similarity of timescales we found in the tide gauge data with those of the Atlantic modes of variability suggests a connection between climate variability and sea-level variations in the Mediterranean Sea as

also indicated by [53,54]. The combined analysis of global oceanic circulation, climate models and observations, indeed, indicates that the main physical mechanism driving the connection between sea level in the Mediterranean Sea and internal modes of natural variability is the mass exchange through the Strait of Gibraltar [53]. This mechanism generates a coherent and uniform signal in the whole Mediterranean at decadal time scales [54,55]. The detected sea-level variations at 20–30 years can be thus related to the combined effect of the North Atlantic Oscillations (NAO) and AMO [26]. An additional contribution could come from the Inter-decadal Pacific Oscillation (IPO), causing variations of the global surface temperature at these timescales, which affects also the climate dynamics in the Atlantic Ocean [56,57]. Likewise, major volcanic eruptions, affecting the global temperature by acting on the radiative forcing [57,58], could also contribute to the SLNV. Fluctuations at longer timescales can be due to both atmospheric forcing [59] and exchange processes [47].

**Figure 3.** Function *LLT*(*t*) (black line), obtained by summing up IMFs at timescales >15 years and the residue, and the best-fit function *LM*(*t*) (red line) described by Equation (4).

A first check of the EMD results can be performed by comparing the SLNV obtained for the four closest stations in the North Adriatic sea, namely Venice, Trieste, Bakar, and Rovinij that are separated by a minimum and maximum distance of ~65 and ~150 km, respectively. Figure 4 shows a comparison of the EMD reconstructions, through IMF with timescale >15 years (the residue *rn*(*t*) has not been added up), for these records. For the four stations the variability is very similar with maxima and minima almost in phase. Small differences can be due to border effects for shorter series (e.g., Rovinij) and/or presence of data gaps affecting the IMF calculation or natural reasons such as small local differences in the signal variability. We note that the Venice station shows a different behavior with a flatter signal between 1945 and 1980 due to temporary anthropogenic subsidence induced by groundwater [28].

Results show that, for all the stations, the SLNV is described by a small number of EMD modes (from 3 to 6). Since IMFs are oscillating functions, we modelled the long-term behavior of the tide gauge signal as the superposition of a linear trend and the sum of oscillating functions of given amplitude, frequency and phase:

$$L\_M(t) = (rt + c) + \sum\_{i=j}^k A\_j \times \cos[\Omega\_j t + \Phi\_j]. \tag{4}$$

**Figure 4.** EMD reconstruction by using IMFs with timescale >15 years (the residue *rn*(*t*) has not been added up) for the four closest stations in the North Adriatic Sea.

In Equation (4) the trend is linear, and the long-term fluctuations are accounted by oscillating functions whose frequency Ω*j* is the time average of the instantaneous frequency of IMFs in *S*. This represents the most elementary approximation to reproduce the long-term variability inferred by the EMD through modes with time-dependent amplitude and phase. The free parameters in Equation (4) (rate *r*, intercept *c*, *Aj* and Φ*j* ) are calculated by fitting *LM*(*t*) to *LLT*(*t*). To avoid possible end effects due to the procedure of IMF extraction and the computation of derivative at borders, we cut three years of data at both boundaries of *LLT*(*t*) before performing the fit. In addition, since data from Trieste, Genoa, Venice, and Bakar show gaps in the initial part of their records, for these stations the fit starts from year 1901, 1931, 1914, and 1949, respectively, when recordings begin to be continuous for many decades. Two sources of uncertainty were taken into account. Firstly, the uncertainty due to the use of constant frequencies, instead of time varying ones, is evaluated by building 1000 *LM*(*t*) realizations whose Ω*j* are randomly extracted from a Gaussian distribution with mean < ω*j* (*t*) > and standard deviation σ [ω*j* (*t*)]. Secondly, effects of possible undetected modes of climate variability, with periods of the order and/or longer than the record lengths, are taken into account as an increase of the estimated uncertainty. Indeed, since these very long period modes are not included in the model, the accuracy for estimating the sea level is reduced. This additional uncertainty is accounted for by adding to each realization *LM*(*t*) a normally distributed zero average white noise. Since we expect that modes with periods of the order and/or longer than the record lengths are settled on the EMD nonlinear residue, the standard deviation of the white noise signal was taken as the standard deviation of the EMD residue from which we subtracted a linear trend.

Figure 3 shows the best solutions obtained by fitting *LM*(*t*) to *LLT*(*t*). Fit parameters and adjusted R-square statistics are shown in Table 2. R-square values close to 1 indicate that a good fit is achieved for the analyzed records. Rates *r* obtained by a fit through Equation (4) are compatible or slightly higher than rates estimated by a linear fit to the raw data for the same time intervals shown in Figure 2 and varying between 0.35 ± 0.38 mm a−1 (Split M.) and 2.43 ± 0.23 mm a−1 (Venice P.S.). Only for Split Gradska Luka, the linear rate from Equation (4) is higher than the linear fit. This is possibly due to a steeper sea-level rise observed in *LLT* at this station after 2006. In the raw data, this feature is hidden by short-term fluctuations, thus resulting in a lower rate when a linear fit is computed on the original time series. The origin of this steepening has been proposed to be the coincidence of AMO-NAO phase opposition and warm AMO phase [60]. The same feature is also present in *LLT* for the Trieste and Rovinj data. However, for Trieste the linear rate from Equation (4) does not show significant differences

with respect to the rate of the linear fit likely due to the longer duration of recordings for which the fitting procedure is very robust and not affected by variations in the last part of the data.

**Table 2.** Fit parameters and adjusted R-square values for the analyzed data set. Uncertainties define the 90% confidence limits.

#### **4. Discussion: The Vertical Land Motion and Relative Sea Levels by 2050 and 2100**

The VLM trends were derived from the processing of raw data collected by regional continuous GPS networks operating in the Euro-Mediterranean region for more than 7 years [61] and located near the tide gauges. The geodetic rates of VLM were obtained from the analysis of position-times series, following the procedures described in [62], according to the IGS recommendations (http://acc.igs.org/reprocess2), to produce solutions with the most updated standards and models. GPS velocities were estimated and represented in the absolute geocentric reference frame (the IGb08 realization of the global ITRF08 [63]. Measured VLM and general information about GPS stations are shown in Figure 5 and Table 3. Due to the poor quality of the data, the VLM at ROGS (Bakar) and RIGS (Rovinj) were not used in this study.

**Figure 5.** Map of the absolute (i.e., IGb08 frame) vertical velocities for the continuous GPS (cGPS) stations used in this work. Black dots refer to the stations ROGS (Bakar) and RIGS (Rovinj) for which, due to the limited available data, calculated vertical land motion (VLM) are not reliable.

**Table 3.** VLM measured at the cGPS station closest to the tide gauges used in this work. The cGPS station name, position, rate of VLM, time interval of the time series, mean VLM, and the closest tide gauges to the cGPS stations are indicated. Negative sign in the VLM corresponds to land subsidence.

The pattern of vertical GPS rates is variable: the VLM (the average value is considered when more than one cGPS station is available) is almost null in the Gulf of Trieste; it is ≤1 mm a−1 at Marseille, Genova, and Split and it is about 1.5 mm a−1 at Dubrovnik. Sudden temporary tectonic episodes, which can cause permanent signatures during the tidal recordings at these stations, have not occurred (www.psmsl.org).

On the other hand, vertical GPS velocities in the North Adriatic coast around the Venice lagoon, display significant subsidence rates, with values in agreement with previous studies [32,64]. The average VLM in the region, calculated as the average over 14 GPS stations, is −3.30 ± 0.85 mm a−1 . Part of the detected subsidence is due to the continuing global glacio-hydroisostatic signal, estimated in this region to be between −0.12 and −0.21 mm a−1 [15], and the remaining part can be related to natural and anthropogenic effects [64,65]. We can reasonably assume that the VLM rate will remain unchanged until 2100 AD, in the absence of additional episodes of land movement that may result from eventual earthquakes with significant magnitude in the area [66] (http://www.isc.ac.uk/). The remarkable subsidence in the Venice lagoon results in a steeper sea-level rate with respect to those provided by IPCC estimations in the area [29,67].

Equation (4) can be used to evaluate sea levels at future times by including both the sea-level trend and VLM, included in *r*, and the effects of the natural long-term variability, reproduced by the oscillating contributions. By keeping the linear rate derived from a linear fit on the raw data, we implicitly assume a constant trend scenario for which the sea-level trend remains unchanged for the next decades. This represents, of course, a very unlikely scenario since there is high confidence that the rate of sea-level rise is increasing (IPCC – AR5). Therefore, to evaluate more realistic sea levels, a time-dependent rate, obtained from the local IPCC projections and including the VLM contribution is substituted to the constant rate *r* in Equation (4). Concerning the modeling of the SLNV, since the Atlantic climatic variability is quasi-persistent on secular times [51,52] and since it is connected to the sea-level variations in the Mediterranean sea (see Section 3), we can reasonably assume that the modes of oscillation associated with the natural variability observed in the past, can be considered as reliable for the next ~90 years. The contribution of the natural variability is also significant in a system dominated by a forced CO2 emission, as the Earth nowadays and possibly in future years. Studies on the effects of the long-term variability on regional sea levels from different climatic models indicate that the internal variability contribution is already relevant when 100 year integrations are considered [68–70]. Moreover, the estimated time of emergence, namely the time at which the CO2-forced sea level signal starts to dominate the internal variability, is close to or longer than a century in most of the Oceans. Other studies [71,72] indicate that the interplay of natural and forced variability in sea-level changes is expected to persist in time throughout the 21st century. In addition, possible variations of the modes of natural variability, between 2016 and 2100, were taken into account in our analysis by using modes with several frequencies for the uncertainty calculation.

To check the robustness of our approach in estimating future sea levels and the accuracy of the SLNV evaluation, we carried out a simple test using data from Marseille, the station with the longest and almost complete time series. Firstly, we applied the EMD analysis on sea-level data collected in the interval 1888–1957 (70 years). Two IMFs with τ*j* > 15 years (τ9 = 33.4 ± 5.4 years and τ10 = 58.1 ± 7.6 years) were found. Then, sea levels for the next 50 years (1957–2007) were evaluated through the model, by using the linear trend from the fit to Equation (4) and were compared to the real complete data set.

Figure 6 shows real data, both reduced at the 1 year of resolution (red dashed line) and smoothed by a 15-years sliding window (red full line), and the model curve (black line) with the 90% confidence interval (enclosed by black dashed lines). The model reproduces well variations in both curves derived from the full time series. For a quantitative comparison:

**Figure 6.** Sea levels in the period 1888–2007 at Marseille as estimated by deriving the model from the first 70 years (1888–1957) of tidal data. Fit parameters are: *r* = 1.57 ± 0.04; *c* = 3882 ± 75; *A* = 13.8 ± 1.0, 9.8 ± 1.01; Φ = 1.22 ± 0.06, 1.89 ± 0.09, *R* 2 = 0.91. The model (black line) with 90% confidence interval (enclosed by black dashed lines), sea-level data at 1 year of resolution (red dashed line) and smoothed by a 15-year sliding window (red line), are shown. The vertical bar marks the boundary between sea level data used for the model (left) and estimated sea levels (right).

The test shows that the modeling of the last 50 years of sea levels, based on the previous 70 years of data, fits the real sea levels within a 90% confidence interval well, thus indicating that the approach used is reliable.

To get sea-level projections at future times, time dependent rates from regional IPCC AR5 sea-level projections (spatial resolution: 1◦ × 1 ◦ ), discussed in the Fifth Assessment Report of the IPCC-AR5 [73,74], and available from the Integrated Climate data Center-ICDC of the University of Hamburg (http://icdc.cen.uni-hamburg.de/1/daten/ocean/ar5-slr.html), were considered. Rates obtained under the Representative Concentration Pathways RCP2.6 and RCP8.5, providing the least and most amounts of future sea level rise respectively, were used. These data consist of mean values and upper/lower 90% confidence bounds obtained by summing up the contributions of several geophysical sources [73,74]: the thermosteric/dynamic contribution, from 21 CMIP5 coupled atmosphere-ocean general circulation models (AOGCMs); the surface mass balance and dynamic ice sheet contributions from Greenland and Antarctica; the glacier and land water storage contributions; the GIA and the inverse barometer effect. For each tide gauge record the IPCC regional sea-level rate was evaluated at the grid point closest to the location of the corresponding station. To take into account the VLM, we substituted the modeled GIA contribution to the IPCC rate with the measured cGPS rate that includes both real GIA and tectonics/anthropic vertical motion. When more than one cGPS station was present, the average VLM value was used (see Table 3). The IPCC rate, corrected for the measured VLM, represents the function *r* in Equation (4). Uncertainties for the sea-level estimations were then calculated by combining lower and upper sea level bounds from IPCC projection, errors from GPS measurements, and uncertainties from the model. Note that for the Bakar and Rovinj sea levels we used the modeled GIA component calculated from the mean of the ICE-5G model [75] and the Australian National University's (ANU) ice sheet model [76] and subsequent improvements, since the measured VLM is not reliable.

Computed sea levels for the two considered IPCC scenarios are shown in Figure 7 and values at 2050 and 2100 (with respect to 2005) are indicated in Table 4. Let us first discuss the results without considering Venice P.S. A mean sea-level rise of 198 ± 102 mm by 2050 and 570 ± 255 mm by 2100 in the most severe RCP8.5 scenario, and 174 ± 97 mm and 342 ± 211 mm by 2050 and 2100, respectively, for the RCP2.6 scenario results from our modelling. In the most favorable scenario (RCP2.6) a sea-level rise between 142 ± 82 (Trieste) and 191 ± 106 (Split M.) mm by 2050, and between 259 ± 165 (Bakar) and 445 ± 200 (Dubrovnik) mm by 2100 is expected.

**Table 4.** Sea levels for 2050 and 2100, with respect to 2005, at the tide gauge for the two scenarios RCP2.6 and RCP8.5 including both the VLM and sea level natural variability (SLNV) contribution. Uncertainties define the 90% confidence limits. Time span used for the modeling and duration of raw records (in parenthesis) are also indicated. Since cGPS data are not reliable, the contribution of VLM is not included for Bakar and Rovinj.

The most severe scenario projections indicate a sea-level rise between 150 ± 86 (Trieste) and 246 ± 95 (Dubrovnik) mm by 2050 and between 475 ± 203 (Bakar) and 681 ± 246 (Dubrovnik) mm by 2100. The combined effect of the VLM and SLNV produces an average increase of the sea level, with respect to pure IPCC projections, of 23/58 mm in RCP2.6 scenario and 37/65 mm in RCP8.5 scenario for 2050/2100 corresponding to a variation of about the 20% (2050) and 13% (2100). The maximum sea-level variation, with respect to pure IPCC projections, of 146/205 mm in RCP2.6 and 150/222 mm in RCP2.6 for 2050/2100 is found at Dubrovnik where the highest VLM rate is measured. The natural variability produces an average variation on the estimated sea level of about 20 (2050) and 10 (2100) mm with maximum values reaching about 65 mm with respect to the IPCC projections. Its average contribution, for the analyzed stations, is about the 9% in 2050 and 2.5% in 2100 (note that the SLNV contribution, in percentage, decreases when the effects of global changes and VLM increase). We remark that, although the contribution of the natural variability in the analyzed stations is of the order of few centimeters, its contribution is not negligible especially for coastlines that are more prone to marine flooding. In these areas, indeed, changes of some centimeters in projected sea levels can change the flooding scenarios [33].

**Figure 7.** Sea-level height up to 2100 for the analyzed stations, as obtained by the model of Equation (4) that includes VLM and AR5 RCP2.6 (blue line) and RCP8.5 (red line) rate. Color bands represent the 90% confidence interval, obtained by including uncertainties from the model, VLM, and AR5 trend. The black line represents sea levels when the constant trend obtained from fits is used. Since cGPS data are not reliable, the contribution of VLM is not included for Bakar and Rovinj.

The sea level at Venice deserves a separate discussion. At Venice P.S., future sea levels are strongly driven by the land subsidence that accounts alone up to 68% (2050) and 60% (2100) of the sea-level rise. VLM, together with the natural variability, accounts for up to 78% (2050) and 64% (2100) of the level rise. In the RCP2.6/RCP8.5 scenarios, the expected sea levels are 283 ± 103 mm/311 ± 114 mm (2050) and 603 ± 217 mm/818 ± 258 mm (2100) well above the values found for the other analyzed records. These values exceed, by about 220 (2050) and 385 mm (2100), the values obtained for Venice when pure IPCC scenarios are considered. These results have considerable implications for the area of the Venice lagoon, characterized by high-density of population, valuable historical setting, residential and industrial infrastructures, besides areas for agricultural use and national parks, for which sea-level rise can become a serious hazard.

The sea-level projections obtained in this study were finally compared with previous investigations in the same geographical areas. References [15,32] estimated the sea-level rise in the Venice area by making use of IPCC predictions, GIA modeling, and including the VLM contribution from geological and archeological data. By using the global IPCC 2007, in [15] the authors found a sea-level rise in 2100 (relative to the sea level average 1986–2005) of about 315 mm; an average sea level between 580 (minimum) and 996 mm (maximum) was estimated by [32], using the global IPCC AR5 RCP8.5 scenario. By modeling and combining the contributions of terrestrial ice melt, GIA and steric sea-level components, in [26] the authors elaborated a minimum and maximum sea-level rise scenario for sea level in 2050: sea levels at Marseille, Venice P.S., Trieste, Bakar were 60, 113, 113, 113 mm and 221, 227, 227, 229 mm respectively. Reference [14] presented a worldwide set of local sea-level projections up to 2100 and 2200 by modeling the contribution to the sea level through three ice sheet components, glacier and ice cap, oceanographic processes, and land water storage. A long-term, local, non-climatic sea-level change calculated through a Gaussian process model, based on historical tide gauge data, was also included to estimate the non-climatic contribution. A quite different approach consists of using a dynamical system model with spatial analysis capability [77–79] to assess the spatial variability of both temperature and sea-level changes over the 21st century. In this framework, temperature and

sea level are considered as coupled variables whose behavior is described by two ordinary differential equations. After calibration through historical data, the model is able to provide reliable predictions of the two variables for the next century. For the Mediterranean Sea region, in particular, an average sea-level rise below 600 mm, relative to 1990, is expected.

The sea-level value from [15] (315 mm) is significantly lower than our RCP2.6 result for 2100 (603 ± 217 mm). This can be ascribed to the use of a global not upgraded IPCC model and an estimation of the VLM from geological and archeological data instead of GPS measurements. On the other hand, our RPC 8.5 results for 2100 are in agreement, within the error limits, with the results by [32]. Additionally, in this case small differences can arise from the use of a global IPCC model (instead of a local one) and the VLM estimation from geological data. Concerning the analysis proposed by [26], sea levels from their maximum scenario are in agreement with our predictions for both RCP2.6 and 8.5. Differences (mainly with their minimum scenario) possibly arise from differences in modeling ice melt, GIA, and steric components with respect to IPCC projections and from the lack of inclusion of the VLM contribution in their study.

Our results are in agreement with the mean sea-level values of [14], within the uncertainty limit, for both 2050 and 2100 at the nine analyzed locations. Differences between the two results are smaller (at most 10%) for the records spanning a longer period (Marseille, Trieste, and Genova) which indicates that both analyses seem to converge when longer time series are considered. We recall that both the approaches use observationally based rates to account for local effects in the sea-level predictions and both the EMD as well as the statistical approach from [14] are more accurate when longer data points are available. Small differences can also be due to differences in modeling ice sheet components, glacier and ice cap, oceanographic processes, and land water storage with respect to IPCC projections. Regarding the results from the Venice station, differences reaching 27% were found. Due to the strong subsidence affecting only the Venice area, and not the entire North Adriatic Sea, we believe that an approach that takes into account the punctual vertical land motion measured closer to the tide gauge station is more suitable for modeling this effect on the sea-level. A comparison of the sea level in 2100 in the Mediterranean Sea area from [77] shows very good agreement with our results when the RCP85 scenario is considered. Indeed, when averaged over all stations (excluding Venice) and reported to the same time basis (year 1990) as [77] we found an average sea-level of 626 ± 255 mm that fits well, within the uncertainty limits, the value found by [77].

Finally, the modeling of the effect of the natural variability on future sea-level projections, present in this study, can also introduce, to a lesser extent, differences with sea-level estimations from the works mentioned above.

#### **5. Conclusions**

Sea levels for 2050 and 2100 at nine locations in the Central-Northern Mediterranean coasts were estimated from tide gauge data by including the contribution of both the vertical land motion and the natural variability. The VLM was inferred at a set of locations nearby tide gauges through continuous GPS measurements. The SLNV was estimated by the EMD analysis on the tide gauge data. The EMD approach revealed modes at characteristic time scales in the ranges 20–30 and >45 years, associated with the long-term natural variability, resulting from the combined effect of teleconnection patterns. Expected sea levels in 2050 and 2100 were then estimated by using IPCC projection, for the time dependent sea-level rate, the measured VLM rate, and a superposition of periodic functions, whose properties were obtained from the EMD results, to model the contribution of the long-term variability. Results (excluding the Venice P.S. tide gauge station) provided a mean sea-level rise of 198 ± 102 mm by 2050 and 570 ± 255 mm by 2100 in the most severe RCP8.5 scenario, while a lower bound is estimated to be 174 ± 97 and 342 ± 211 mm by 2050 and 2100, respectively, for the RCP2.6 scenario. The contribution of natural variability to local sea-level variations is significant and can reach, on average, up to the 9% with respect to the sea-level changes induced only by global climate change. VLM and SLNV together account, on average, for about the 15% of the sea-level variation.

In the Venice lagoon, the land subsidence is increasing the local sea-level rise. By taking into account VLM rates measured by the cGPS stations in this area, relative sea levels are expected to reach 603 ± 217 and 818 ± 258 mm in 2100 for the RCP2.6 and RCP8.5 scenarios, respectively. In particular, at the Venice P.S. tidal station, VLM and SLNV contribute more than the 60% to the expected sea-level rise.

Finally, our analysis shows that the contribution of local effects, such as the natural variability and the vertical land motion, play a key role in local sea-level rise projections. Effects are even more relevant for low elevated coastal areas where marine ingression may represent a potential hazard for the environment and human activities. Even in the more optimistic scenario, the estimated sea-level rise will have an important impact along the coasts, causing diffuse erosion of the shorelines, hence determining serious impacts on many coastal areas, increasing the hazard related to flooding events, storm surges, and possibly tsunamis. In subsiding areas, like the Venice Lagoon, severe environmental impact with subsequent loss of economic value of shores and coastal infrastructures can be reasonably expected before 2100. In this regard, land planners and decision makers should take into account the rising sea levels reported in this study for cognizant coastal management.

**Author Contributions:** Conceptualization, A.V. and M.A.; Formal analysis, A.V. and E.S.; Funding acquisition, M.A.; Methodology, A.V.; Project administration, M.A.; Supervision, M.A. and F.F.; Validation, A.V.; Writing – original draft, A.V.; Writing—review & editing, M.A., E.S. and F.F.

**Funding:** This work was carried out under the umbrella of the SAVEMEDCOASTS project "Sea level rise scenarios along the Mediterranean coasts", funded by the EU (www.savemedcoasts.eu; agreement number ECHO/SUB/2016/742473/PREV16).

**Acknowledgments:** We acknowledge two anonymous reviewers for the useful comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Beach-Foredune Sediment Budget Response to Sea Level Fluctuation. Curonian Spit, Lithuania**

#### **Darius Jarmalaviˇcius \*, Donatas Pupienis, Gintautas Žilinskas, Rasa Janušaite and ˙ Viktoras Karaliunas ¯**

Laboratory of Geoenvironmental Research, Nature Research Centre, Akademijos 2, LT-08412 Vilnius, Lithuania<EMAIL_ADDRESS>(D.P<EMAIL_ADDRESS>(G.Ž<EMAIL_ADDRESS>(R.J<EMAIL_ADDRESS>(V.K.)

**\*** Correspondence<EMAIL_ADDRESS>
Received: 14 January 2020; Accepted: 17 February 2020; Published: 20 February 2020