Actual and Forecasted Vulnerability Assessment to Seawater Intrusion via GALDIT-SUSI in the Volturno River Mouth (Italy)

: Coastal areas have become increasingly vulnerable to groundwater salinization, especially in the last century, due to the combined effects of climate change and growing anthropization. In this study, a novel methodology named GALDIT-SUSI was applied in the ﬂoodplain of the Volturno River mouth for the current (2018) and future (2050) evaluation of seawater intrusion accounting for the expected subsidence and groundwater salinization rates. Several input variables such as digital surface model, land use classiﬁcation, subsidence rate and drainage system have been mapped via remote sensing resources. The current assessment highlights how areas affected by salinization coincide with the semiperennial lagoons and inland depressed areas where paleosaline groundwaters are present. The future assessment (2050) shows a marked increase of salinization vulnerability in the coastal strip and in the most depressed areas. The results highlight that the main vulnerability driver is the Revelle index, while predicted subsidence and recharge rates will only slightly affect groundwater salinization. This case study indicates that GALDIT-SUSI is a reliable and easy-to-use tool for the assessment of groundwater salinization in many coastal regions of the world.


Introduction
Coastal zones (CZ) represent exceptionally sensitive environments worldwide which have been deeply influenced by natural processes and human activities capable of modifying the landscape pattern structure and the ecological processes in the related ecosystems [1,2]. Climate change, especially sea-level rise (SLR) [3], severe drought and anomalous waves [4] are the natural processes that most affect coastal systems [5], and anthropic pressures, namely high population density and widespread agricultural and touristic activities [6,7], have been identified as important drivers responsible for coastal systems degradation, especially in deltaic coastal settings.
In particular, seawater intrusion (SI) and more generally groundwater salinization, are serious phenomena that could negatively affect these delicate settings, being the main discriminating factor for the groundwater quality in coastal aquifers worldwide. In the Mediterranean coastal areas, anthropogenic activities requiring large quantities of water, which is usually met by abstracting groundwater resources, have led to continuous groundwater depletion, triggering SI. Moreover, among all factors, SLR in Mediterranean coastal aquifers has been identified as one of the main reasons responsible for actual and future extensional sedimentary basin formed during the Quaternary, between the western broadside of the Southern Apennines and the Eastern Tyrrhenian margin [28]. During the Pleistocene, the CP was characterized by the onset of a rapid tectonic subsidence and consequently it was filled with continental and marine deposits (alluvial, fluvial and lagoon deposits), interfingered with volcanic deposits related to the explosive eruptions of Somma-Vesuvius, Roccamonfina and Campi Flegrei districts [29]. Among the volcanic events that characterized the area, the 39 ka pyroclastic flow deposits vented by the Campi Flegrei (Campania Grey Tuff-CGT) blanketed the entire plain and represented the substrate for the Holocene sedimentation that followed [30]. These deposits were deeply eroded by the paleo-Volturno River in response to the Last Glacial Maximum sea-level drop [31,32]. During the Pleistocene-Early Holocene (ca. 6-15 ka) the post-glacial rising of sea level caused a broadening of the inner shelf with a rapid flooding of the lower part of the Volturno plain. Since 6.5 ky, a coastal pro-gradation phase, established at the termination of the transgressive phase, allowed the formation of a wave-dominated delta system with flanking strand plains forming beach-dune systems that partially enclose shallow lagoons and marshes ( Figure 2). The present-day morphology is largely inherited by the Holocene sedimentary evolution. The Volturno River mouth is a low-lying area, with mean ground elevation above sea level (a.s.l.) of approximately −2.0 to 4.0 m and characterized by semipermanent water bodies (lagoons and wetlands). Here, the aquifer is hosted in sandy dune deposits interlayered with silt lens. During the Pleistocene, the CP was characterized by the onset of a rapid tectonic subsidence and consequently it was filled with continental and marine deposits (alluvial, fluvial and lagoon deposits), interfingered with volcanic deposits related to the explosive eruptions of Somma-Vesuvius, Roccamonfina and Campi Flegrei districts [29]. Among the volcanic events that characterized the area, the 39 ka pyroclastic flow deposits vented by the Campi Flegrei (Campania Grey Tuff-CGT) blanketed the entire plain and represented the substrate for the Holocene sedimentation that followed [30]. These deposits were deeply eroded by the paleo-Volturno River in response to the Last Glacial Maximum sea-level drop [31,32]. During the Pleistocene-Early Holocene (ca. 6-15 ka) the post-glacial rising of sea level caused a broadening of the inner shelf with a rapid flooding of the lower part of the Volturno plain. Since 6.5 ky, a coastal pro-gradation phase, established at the termination of the transgressive phase, allowed the formation of a wave-dominated delta system with flanking strand plains forming beach-dune systems that partially enclose shallow lagoons and marshes ( Figure 2). The present-day morphology is largely inherited by the Holocene sedimentary evolution. The Volturno River mouth is a low-lying area, with mean ground elevation above sea level (a.s.l.) of approximately −2.0 to 4.0 m and characterized by semipermanent water bodies (lagoons and wetlands). Here, the aquifer is hosted in sandy dune deposits interlayered with silt lens.
has been reduced, resulting in erosional processes that have led to the loss of hectares of beach [35].
The site selected for the present study corresponds to the area of maximum ingression during the Holocene transgression (Figures 1 and 2), characterized by marine and transitional environments and currently resting largely below the sea level ( Figure 3). The area has a typical Mediterranean climate with wet winters and dry summers, an average annual temperature of 21.0 °C and precipitation of 800-1000 mm/year [35].  A pro-delta formation made of silty clay deposits is present starting from 16-20 m below ground level (b.g.l.); it represents the impermeable bed of the coastal aquifer [30] ( Figure 3). High subsidence rates characterize the entire coastal environments (from −5 to −25 mm/year) with the most affected area corresponding to the ancient buried lagoon, which was filled by highly compactable peat and clay-rich sediments [9].
The dynamics of the whole coastal system dramatically changed in the last century following reclamation interventions that completed the canalization of most secondary streams in the lower alluvial plain [33]. The availability of reclaimed lands along the coastal alluvial plain promoted the development of agriculture and farming, together with strong coastal urbanization. Therefore, landscape fragmentation increased significantly between the 1960s and the 1990s, as well as built-up land area that overgrew to the sea. An overall reduction of high-quality ecosystems (humid coastal setting, lacustrine/marshy back-dune area and beach-dune system) occurred, resulting in biodiversity loss and a dramatic reduction of environmental quality [34]. Moreover, sediment supply to the coast has been reduced, resulting in erosional processes that have led to the loss of hectares of beach [35].
The site selected for the present study corresponds to the area of maximum ingression during the Holocene transgression (Figures 1 and 2), characterized by marine and transitional environments and currently resting largely below the sea level ( Figure 3). The area has a typical Mediterranean climate with wet winters and dry summers, an average annual temperature of 21.0 • C and precipitation of 800-1000 mm/year [35].

GALDIT-SUSI
GALDIT-SUSI, proposed by Kazakis et al. [23], considers the six standard GALDIT parameters, such as the characteristics of the aquifer, the morphology and the vadose zone, together with the presence/absence of superficial water bodies (lagoon, torrents and river) to determine the SI vulnerability of coastal aquifers. The SUSI modification allows accounting for the saline water that can indirectly intrude into the aquifer flowing from surface water bodies towards groundwater. Moreover, the introduction of the elevation parameters also allows accounting for the effect of naturally depressed areas along with water bodies bathymetry.
The six additional parameters used to implement the GALDIT-SUSI method are (1) vadose zone hydraulic conductivity (K), (2) elevation, (3) torrents, (4) rivers, (5) wetlands and (6) lagoons (Table S1). The initial weights of the overall twelve parameters were defined using the analytic hierarchy process [36]. Table S1 shows the weight and class range for each parameter. The final SI vulnerability map of GALDIT-SUSI was produced using the following formula: where the parameters are R and W their corresponding weights from Table S1. All twelve parameters were constructed and classified in a GIS environment retaining a cell resolution of 2 × 2 m. The existing status of SI was calculated in 44 groundwater and surface water samples using the Revelle index (RI) [37].

Data Collection and Analysis
To obtain the main geomorphological features characterizing the study area, extensive cartographic documentation (topographic maps, ortho images, aerial photographs) was acquired covering a period of about 200 years; the maps were acquired as raster files in tiff format, 300 and 1200 dpi and georeferenced through ground control points (GCPs)

GALDIT-SUSI
GALDIT-SUSI, proposed by Kazakis et al. [23], considers the six standard GALDIT parameters, such as the characteristics of the aquifer, the morphology and the vadose zone, together with the presence/absence of superficial water bodies (lagoon, torrents and river) to determine the SI vulnerability of coastal aquifers. The SUSI modification allows accounting for the saline water that can indirectly intrude into the aquifer flowing from surface water bodies towards groundwater. Moreover, the introduction of the elevation parameters also allows accounting for the effect of naturally depressed areas along with water bodies bathymetry.
The six additional parameters used to implement the GALDIT-SUSI method are (1) vadose zone hydraulic conductivity (K), (2) elevation, (3) torrents, (4) rivers, (5) wetlands and (6) lagoons (Table S1). The initial weights of the overall twelve parameters were defined using the analytic hierarchy process [36]. Table S1 shows the weight and class range for each parameter. The final SI vulnerability map of GALDIT-SUSI was produced using the following formula: where the parameters are R and W their corresponding weights from Table S1. All twelve parameters were constructed and classified in a GIS environment retaining a cell resolution of 2 × 2 m. The existing status of SI was calculated in 44 groundwater and surface water samples using the Revelle index (RI) [37].

Data Collection and Analysis
To obtain the main geomorphological features characterizing the study area, extensive cartographic documentation (topographic maps, ortho images, aerial photographs) was acquired covering a period of about 200 years; the maps were acquired as raster files in Remote Sens. 2021, 13, 3632 6 of 18 tiff format, 300 and 1200 dpi and georeferenced through ground control points (GCPs) technique, using as a base map the georeferenced Topographic Regional Map (CTR) of the Campania Region (1:25.000) and GeoMedia ® Pro 6.1 with "Image Registration" tool, which allowed for resampling the maps by bilinear interpolation. All new data were reprojected to UTM Zone 33 (WGS84). The accuracy of the maps (Table 1) could be affected by errors induced either by digitization and georeferencing, the latter being more common if the original map was realized in different coordinate systems. In this study, the root mean square error (RMSe) calculated for all the maps used to draw the land use/land cover (LULC) polygons is lower than the pixel resolution. Multitemporal data were topologically overlapped to facilitate spatial analysis (change, persistence and trend). Table 1. Input data used to obtain land use and land cover (LULC), hydrographic network (HN), coastline (C), subsidence trends (ST), and digital terrain model (DTM). The year of data acquisition, its resolution and root mean square error (RMSe) are also reported.
All the data for the GALDIT-SUSI application were collected from previous literature and online regional databases. The information to define the aquifer characteristics, such as (i) aquifer type, (ii) vadose zone materials, (iii) aquifer thickness and (iv) vadose zone/aquifer K values, were extracted from the interpretation of several geological boreholes collected from a local institutional database and literature [9,35]. Major ions in surface waters and groundwater samples were collected in a sampling campaign conducted in September 2016 on 44 monitoring stations ( Figure 1) and analyzed using an isocratic dual pump ion chromatograph ICS-1000 Dionex.
The main data processing and study workflow are shown in Figure 4.  . Workflow for spatial data processing and SI vulnerability.

Land Use Map
The land use map was reconstructed for 2012. The land cover categories were identified by visual on-screen photointerpretation and then digitized and classified. All the activities were performed in GeoMedia ® ; the polygons were digitized in a vectorial format and classified with the aim of field surveys. The land use map CORINE Land Cover (CLC, 2012) was used to check and verify the results and obtain the land use maps. The related land use classes were combined or renamed, according to the CLC nomenclature, to produce the final outputs. They were classified into nine categories that could be easily recognized in each time period. To outline the drainage system evolution over time, the topographic maps from 1957 to 1987 were scanned and imported into the GIS environment; the maps were analyzed by manual on-screen interpretation and the main linear elements (i.e., rivers and canals) were digitized. No changes occurred after the beginning of the 1990s.

Coastline Changes
To document the processes and geomorphological changes that occurred in the study area, 150 years of cartographic sources were acquired in raster format with 300 and 1200 dpi resolution and georeferenced through the GCPs technique using as a base map the georeferenced CTR of the Campania Region (1:5000 scale, 2006). Approximately 100 GCPs were identified on each document in order to achieve the most reliable accuracy. These douments were managed in GIS using the software Intergraph Geomedia Pro 6.1 and georeferenced in the UTM Zone 33 (WGS84) coordinate system.
The uncertainties in mapping the shoreline position owing to short-term dynamic cyclic fluctuations [38] were considered, although the available records from which longterm shoreline changes were determined span over 100 years and thus considered reliable, whereas an accuracy of few meters is limited only to the most recent datasets [39]. The coastline position from 1957 to 2018 varied from a maximum of 90 m landward in the central part of the study area at the Volturno River mouth to a maximum of 65 m seaward in the northern part of the study area between the Savone and Agnena rivers. Assuming that the coastline changes in this area will occur at an almost steady trend and rate in the near future, as testified by several long-term monitoring studies in the area [33,39], the maximum expected shift for 2050 will be approximately 48 m landward and 35 m seaward.

Land Use Map
The land use map was reconstructed for 2012. The land cover categories were identified by visual on-screen photointerpretation and then digitized and classified. All the activities were performed in GeoMedia ® ; the polygons were digitized in a vectorial format and classified with the aim of field surveys. The land use map CORINE Land Cover (CLC, 2012) was used to check and verify the results and obtain the land use maps. The related land use classes were combined or renamed, according to the CLC nomenclature, to produce the final outputs. They were classified into nine categories that could be easily recognized in each time period. To outline the drainage system evolution over time, the topographic maps from 1957 to 1987 were scanned and imported into the GIS environment; the maps were analyzed by manual on-screen interpretation and the main linear elements (i.e., rivers and canals) were digitized. No changes occurred after the beginning of the 1990s.

Coastline Changes
To document the processes and geomorphological changes that occurred in the study area, 150 years of cartographic sources were acquired in raster format with 300 and 1200 dpi resolution and georeferenced through the GCPs technique using as a base map the georeferenced CTR of the Campania Region (1:5000 scale, 2006). Approximately 100 GCPs were identified on each document in order to achieve the most reliable accuracy. These douments were managed in GIS using the software Intergraph Geomedia Pro 6.1 and georeferenced in the UTM Zone 33 (WGS84) coordinate system.
The uncertainties in mapping the shoreline position owing to short-term dynamic cyclic fluctuations [38] were considered, although the available records from which longterm shoreline changes were determined span over 100 years and thus considered reliable, whereas an accuracy of few meters is limited only to the most recent datasets [39]. The coastline position from 1957 to 2018 varied from a maximum of 90 m landward in the central part of the study area at the Volturno River mouth to a maximum of 65 m seaward in the northern part of the study area between the Savone and Agnena rivers. Assuming that the coastline changes in this area will occur at an almost steady trend and rate in the near future, as testified by several long-term monitoring studies in the area [33,39], the maximum expected shift for 2050 will be approximately 48 m landward and 35 m seaward.

DTM
The digital terrain model (DTM) was realized by merging the DTM FIRST with 2 × 2 m resolution for the coastal stretch (750 m wide) and 1 x 1 m resolution for the inland area; they were finally resampled at 2 × 2 m resolution for the whole area ( Figure 3b). The data were acquired from LIDAR scanning available from the Ministry of Environment and Protection of the Sea and Territory (http://www.pcn.minambiente.it/mattm/accessed, accessed on 10 January 2020).

Subsidence and Geology
The assessment of the subsidence trends was based on a temporal analysis and mapping of persistent scatterers (PS) data, obtained from interferometric processing of ERS-1/2 (1992-2001), RADARSAT (2003RADARSAT ( -2007 and ENVISAT (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) radar satellite scenes of the study area using permanent scatterers (PS-InSAR) and persistent scatterers pairs (PSP) techniques [9]. The information concerning post-processing, analysis and mapping of the available PS interferometry datasets, derived from the combination of both ascending and descending orbits of three different radar satellite systems during an observation period of almost two decades (1992-2010), are available in Matano et al. [9].
According to the identified subsidence rate, the predicted elevation for 2050 was calculated by subtracting from the DTM the product between the subsidence for the time frame considered, 30 years (2020-2050), using the raster calculator.
To assess the relationship between lithology and subsidence, the stratigraphic architecture of the subsoil was reconstructed on the basis of 400 shallow borehole stratigraphic logs (mostly up to 20-30 m in depth) located along the lower Volturno plain, and on previous studies results [30][31][32][33][34].
Borehole stratigraphic data were interpreted and grouped into three main units: 1.
Pre-CGT-comprises all the lithofacies deposited before the Phlegraean volcanic eruption (39 Ky) and settled below the tuff unit.

2.
CGT-consists of the Campania Grey Tuff volcanic deposits.

3.
Post-CGT-comprises all the sediments deposited after the CGT eruption and that likely constituting Holocene sedimentation.
Among this last unit, the thickness of each lithology characterizing the alluvial/ transitional deposits was also recorded to understand their impact in the subsidence processes. Top and bottom depths of each lithological unit surface (peat, clay, silt, sand and gravel) were stored in a relational database. The ratio between the thickness of peat, clay and silt (the most compressible materials) and the whole post-CGT thickness was then considered. Reference geological sections were realized to compare thicknesses and lithology of the post-CGT deposits and the subsidence rates. Cumulative subsidence profiles were drawn on each profile.

Future Groundwater Quality and Head Variations
The future status of groundwater quality was assessed on the basis of a linear increasing trend of the RI in the whole area, +0.011 per year, derived by analyzing the observed RI in area wells from 2002 to 2018. The trend obtained was also in accordance with the groundwater salinization trend highlighted by previous research in the same study area [27].
Similarly, groundwater head values for the next 30 years were evaluated according to the piezometric records from 2002 to 2018, available for five monitoring wells located inside the study area and showing no increasing or decreasing trends, and thus kept fixed at actual conditions for the future SI assessment.

Geochemical Background
Descriptive statistics of all physicochemical parameters measured within the 44 water samples is presented in Table 2. Both surface water and groundwater samples inside the Remote Sens. 2021, 13, 3632 9 of 18 study area mainly belong to the carbonate-alkali and sodium chloride facies, showing a very high salinity content. The high correlation found between Cl − , Na + , Br − and SO 4 2− (r > 0.9) further confirms the presence of salinization phenomena. The ternary diagrams representing the main anions ( Figure 5a) and cations (Figure 5b) indicate that water mineralization is mainly shifted towards the predominance of Cl − in solution, followed by Na + , K + and Mg 2+ , highlighting the dissolution of saline sediment and recent SI only near the coast [26]. Cl − concentrations in the area range from a minimum of 9 mg/L to as much as 19,400 mg/L for the Tyrrhenian Sea. Another characteristic of groundwater samples in the area is the elevated NO 3 − pollution. The average detected value of this pollutant is 33.5 mg/L, and in some samples it exceeds the WHO limit of 50 mg/L, reaching a maximum concentration six times higher than the threshold limit. The presence of NO 3 − pollution can be directly connected to the intensive agricultural activities that widely characterize the CZ of the entire CP [40].

Geochemical Background
Descriptive statistics of all physicochemical parameters measured within the 44 water samples is presented in Table 2. Both surface water and groundwater samples inside the study area mainly belong to the carbonate-alkali and sodium chloride facies, showing a very high salinity content. The high correlation found between Cl − , Na + , Br − and SO4 2− (r > 0.9) further confirms the presence of salinization phenomena. The ternary diagrams representing the main anions ( Figure 5a) and cations (Figure 5b) indicate that water mineralization is mainly shifted towards the predominance of Clin solution, followed by Na + , K + and Mg 2+ , highlighting the dissolution of saline sediment and recent SI only near the coast [26]. Clconcentrations in the area range from a minimum of 9 mg/L to as much as 19,400 mg/L for the Tyrrhenian Sea. Another characteristic of groundwater samples in the area is the elevated NO3 -pollution. The average detected value of this pollutant is 33.5 mg/L, and in some samples it exceeds the WHO limit of 50 mg/L, reaching a maximum concentration six times higher than the threshold limit. The presence of NO3 -pollution can be directly connected to the intensive agricultural activities that widely characterize the CZ of the entire CP [40].

Land Use Changes
The comparison between pre-urbanization and actual land use maps shows that the landscape composition went through major changes from the beginning of 1900s, when land reclamation works resulted in the canalization of most secondary streams in the lower alluvial plain, and croplands and/or grassland covered most of the territory [35]. Since the second half of the XX century, massive urbanization of the coastal area replaced most of the former land cover types. Salt marshes, lagoons and wetland, beach-dune systems and a conspicuous fraction of agricultural land were nearly obliterated and converted into settlement areas in response to population growth and tourism expansion. Unfortunately, the land use change induced by the socioeconomic development exposed the residential areas to a high vulnerability in a few decades since the disappeared pristine environments acted as natural buffers to reduce wave action, storm surges, shoreline erosion and, ultimately, SI [41,42].

Geological Framework and Subsidence Rate
The analysis of the borehole data allowed reconstruction of the upper surface of the GCT in the whole CP, and to recognize paleovalley morphology in correspondence to the course of the modern Volturno River [9]. The vertical ground deformation assessed by Matano et al. [9] documented a continuous subsidence, mostly confined in the paleovalley perimeter ( Figure S1). Inside this perimeter, the vertical ground movements could be considered stable with values within ±10 mm. Looking in detail at the study site, high subsidence rates (>15 mm/y) are distributed in the external sectors of the alluvial plain, and in the back-dune depressions, extending about 8 km 2 (Figure 6a). The coastal sector near the Volturno River mouth, extending for about 10 km 2 , is characterized by moderate subsidence values (−5 to −7 mm/y), while the dune system sector shows very slight subsidence or stability. Considering the near future (2050) (Figure 6b), inside this perimeter, the vertical ground movements can be regarded as stable in most of the territory, with values within ±10 mm. In correspondence with the more subsided areas identified previously, a cumulated subsidence of -418 mm can be reached at the upstream stretch of the Volturno River course.

Land Use Changes
The comparison between pre-urbanization and actual land use maps shows that the landscape composition went through major changes from the beginning of 1900s, when land reclamation works resulted in the canalization of most secondary streams in the lower alluvial plain, and croplands and/or grassland covered most of the territory [35]. Since the second half of the XX century, massive urbanization of the coastal area replaced most of the former land cover types. Salt marshes, lagoons and wetland, beach-dune systems and a conspicuous fraction of agricultural land were nearly obliterated and converted into settlement areas in response to population growth and tourism expansion. Unfortunately, the land use change induced by the socioeconomic development exposed the residential areas to a high vulnerability in a few decades since the disappeared pristine environments acted as natural buffers to reduce wave action, storm surges, shoreline erosion and, ultimately, SI [41][42].

Geological Framework and Subsidence Rate
The analysis of the borehole data allowed reconstruction of the upper surface of the GCT in the whole CP, and to recognize paleovalley morphology in correspondence to the course of the modern Volturno River [9]. The vertical ground deformation assessed by Matano et al. [9] documented a continuous subsidence, mostly confined in the paleovalley perimeter ( Figure S1). Inside this perimeter, the vertical ground movements could be considered stable with values within ±10 mm. Looking in detail at the study site, high subsidence rates (>15 mm/y) are distributed in the external sectors of the alluvial plain, and in the back-dune depressions, extending about 8 km 2 (Figure 6a). The coastal sector near the Volturno River mouth, extending for about 10 km 2 , is characterized by moderate subsidence values (−5 to −7 mm/y), while the dune system sector shows very slight subsidence or stability. Considering the near future (2050) (Figure 6b), inside this perimeter, the vertical ground movements can be regarded as stable in most of the territory, with values within ±10 mm. In correspondence with the more subsided areas identified previously, a cumulated subsidence of -418 mm can be reached at the upstream stretch of the Volturno River course.  Two lithostratigraphic sea-land transects were constructed, focusing on the stratigraphic position of the compressible layers (clays and peat) in relation to the different average subsidence rates (Figure 7). They show anomalous increases in annual average subsidence in correspondence to large thicknesses of peat and clay. The inclusion of significant amounts of peat and organic matter clearly reflects high values of the secondary compression coefficient (creep), that is, the continuous soil deformation under constant effective stresses [43]. To identify the possible origin of the increase of the vertical displacement observed in the abovementioned sectors of the study area and to validate the described trends, the lithology of the post-CGT deposits was considered. In particular, the thickness of materials highly subject to secondary consolidation (i.e., clay, silt and peat) was compared to the subsidence values (Figure 5a). The map shows that zones with high values of subsidence (not corresponding to the absolute maximum values) correspond to a large thickness of peat and clay lithologies in the subsoil. From a geotechnical point of view, these soils can be classified as fine-grained soil characterized by high values of the secondary compression coefficient, which can contribute to soil deformation.
Two lithostratigraphic sea-land transects were constructed, focusing on the stratigraphic position of the compressible layers (clays and peat) in relation to the different average subsidence rates (Figure 7). They show anomalous increases in annual average subsidence in correspondence to large thicknesses of peat and clay. The inclusion of significant amounts of peat and organic matter clearly reflects high values of the secondary compression coefficient (creep), that is, the continuous soil deformation under constant effective stresses [43].

Actual Seawater Intrusion Assessment
The GALDIT-SUSI was applied in the coastal plain of the Volturno River to assess the actual SI vulnerability. In accordance with the available hydrogeological information and the GALDIT-SUSI method ( Figure S2), the local aquifer is classified as (i) unconfined within sand deposits with a bottom impermeable layer of pro-delta silty clay sediments along the coast, and so classified as highly vulnerable (rating 7), (ii) confined inside clay and silty clay deposits, which constitute the alluvial plain of the Volturno River and characterized by an average vulnerability (rating 5) and (iii) bounded in those areas adjoining the Volturno and Savone rivers, with a low vulnerability (rating 2). The K value of the sandy formations, which constitute the main aquifer body for the whole study area, in some places interbedded with silt lenses, varies between 15-30 m/day. Following Table S1, the main aquifer body is classified as highly vulnerable (rating 10). The water table ranges from −2.0 to +8.0 m a.s.l. and all the vulnerability classes are represented in the study area with bands subparallel to the coastline. The most vulnerable zone (rating 10) extends from the north coast to the central eastern part of the study area, while the less vulnerable area (rating 2.5) lies in the northern part of the study area; the remaining territory, more than 50% of the plain, is characterized by an average to high vulnerability to SI (rating 5 and 7.5) due to the low piezometric values, ranging from +1.0 to +3.0 m. a.s.l. The calculated RI ranges from 0.02 to over 200. Highly positive values (>2), corresponding to high and very high vulnerability to SI, characterize the whole coastal zone and some internal areas where paleoseawater trapped in Holocene sediments contributes to increased overall salinity. The remaining alluvial plain is instead characterized by a medium to low vulnerability (negative RI). The aquifer thickness is commonly utilized in many rating methodologies to express pollutant dilution [44]. Generally high values of aquifer thickness express a high dilution capacity, thus the vulnerability of the aquifer to SI is higher where the aquifer is thinner. In the study area, the aquifer thickness, calculated through the interpretation of available lithological boreholes and according to previous studies, is not greater than 30 m [27], which is sufficient to be considered highly vulnerable to SI (rating 10). In the GALDIT-SUSI methodology, additional parameters were included with the standard ones, such as vadose zone K and ground elevation. The elevation is classified using the previously created DTM. In detail, the morphology ranges from −0.05 m to +12.0 m a.s.l. All the areas close to the river mouth up to 6 km inland range from −0.05 to +2.0 a.s.l., while the elevation reaches its maximum in the remaining relict dunes and in the northeastern part of the study area. The depressed morphology, coincident with the ancient reclamation works of Lake Patria and Clanio River, is classified as high (rating 7.5) and very high (rating 10) vulnerability, while the remaining valley is characterized by medium (rating 5) to low (rating 2.5) vulnerability. Another important parameter introduced by the GALDIT-SUSI modification is vadose zone K, which plays a key role in determining the interactions between groundwater and surface water [45]. Where low values of vadose K occur, the interaction is limited and vice versa. In the coastal area, the vadose zone is mainly made of coarse-grained material with a K ranging between 5 and 10 m/d, thus corresponding to medium vulnerability (rating 5), while all the unsaturated materials far from the coastline are silt and silty clay deposits characterized by a low vulnerability (rating 2.5). Regarding the distance from lagoons, rivers, torrents, wetlands and shoreline, the SI vulnerability is considered higher in the areas that lie nearer water bodies or courses. All the distances were calculated using a buffer zone of 75, 150, 300 and greater than 300 m for lagoons, rivers and torrents, while a buffer of 500, 750, 1000 and greater than 1000 m was used for the shoreline. All the parameters previously discussed are represented in the maps of Figure S2.
The composite SI vulnerability map was constructed by using Equation (1) and the raster calculator tool in ArcGIS 10.2. The actual SI assessment map (Figure 8a) subdivides the area into five classes of vulnerability, classified using a geometrical interval, and ranging from very low to very high vulnerability to SI. operations, where the presence of drainage canals contributes to worsen the situation. Moreover, note that this highly vulnerable zone coincides with the greatest thickness of peat and clay sediment (Figure 7) suggesting a controlling role of lithology towards the SI phenomena. In this portion of the study area, a wide part of the territory is classified as highly vulnerable.
The low and very low vulnerability areas are present only in the northern part of the study area, at least 5 km from the coastline, and are divided from the most vulnerable areas by a thin stripe where the vulnerability to SI is medium. The area located between the Regi Lagni River and Lago Patria, characterized by lacustrine sediments and where the subsidence rates are small, is also classified as medium vulnerable.  The highest vulnerability to SI spans the entire coastline, widening inland towards the north of the study area. Very high vulnerability also characterizes the terminal stretches of the Agnena, Volturno and Regi Lagni rivers, up to 4.0 km from the shoreline. In these zones, the rivers act as a preferential way for SI due to their riverbed elevation largely below ground level, especially close to the coastal zone [46]. The low elevation, together with the high salinity index (RI > 2), and the unconfined aquifer condition are responsible for the natural high vulnerability to salinization. In the southeastern portion of the study area, a wide area of very high vulnerability to SI also corresponds with depressed areas affected by a high rate of subsidence and recently involved in reclamation operations, where the presence of drainage canals contributes to worsen the situation. Moreover, note that this highly vulnerable zone coincides with the greatest thickness of peat and clay sediment ( Figure 7) suggesting a controlling role of lithology towards the SI phenomena. In this portion of the study area, a wide part of the territory is classified as highly vulnerable.
The low and very low vulnerability areas are present only in the northern part of the study area, at least 5 km from the coastline, and are divided from the most vulnerable areas by a thin stripe where the vulnerability to SI is medium. The area located between the Regi Lagni River and Lago Patria, characterized by lacustrine sediments and where the subsidence rates are small, is also classified as medium vulnerable.

Future Seawater Intrusion Assessment
The future vulnerability to SI was assessed considering the future impact of subsidence, groundwater heads and salinization rate upon the whole study area. According to the previous elaboration, the subsidence rate in the area ranges from a minimum of 0.0 to a maximum of 20 mm/year. These rates could significantly impact the study area in 2050, especially within those areas slightly above the sea level (0.0-2.0 m), also increasing the erosional processes through the years if the trend remains constant. This assumption is confirmed from the study conducted by Aucelli et al. [47], who predicted the topographic model for the Volturno River plain using the actual subsidence rate. Considering this scenario and applying the predicted trends for 30 years into the future, a maximum of 0.45 m of subsidence (Figure 5b) will be registered in 2050, which will become almost 1.0 m in 2100. Together with subsidence, the increasing trend of salinization, which will affect the whole Campanian coastal zone [27], will surely influence the future vulnerability to SI. The predicted RI and elevation are represented in Figure S3. Concerning the future groundwater heads, despite slight piezometric changes, which will affect only the shoreline strip [8], the SUSI classification (Table S1) for the whole area will not show significant changes in the parameter classification. Thus, the same map was used for the future evaluation of SI. A similar result was obtained for the coastline changes.
The final step was the assessment of SI for 2050 using GALDIT-SUSI retaining the same weight and rate used for the actual SI evaluation. The composite vulnerability map was so calculated using Equation (1) and represented in Figure 8b. The future SI vulnerability map divided the entire valley in five classes of vulnerability: very high, high, medium, low and very low.
From a first glance at the map, it is already clear that the very high and high vulnerability classes will characterize a wider area at the expense of the medium vulnerability. Specifically, the very high vulnerability areas that line the Agnena, Volturno and Regi Lagni rivers further expand laterally, while the areas of high vulnerability expand to the detriment of the areas of medium vulnerability, the latter having almost completely disappeared in the map of future vulnerability to SI. The low and very low vulnerability will remain confined in those areas characterized by groundwater levels markedly above sea level and higher topographical elevations.
The changes recorded in the distribution of the vulnerability to SI, and in particular its exacerbation, are mainly attributable to the increase in RI, which in 2050 will be up to three time stronger than 2018 [27] and, even if to a lesser extent, to the impact of subsidence. On the other hand, the decrease in recharge [48][49][50] with the relative lowering of the water table and coastal erosion will only marginally affect the enlargement of highly vulnerable areas. The predominant role of RI ultimately derives from the expected sea-level rise, which will push highly saline sea waters farther into the coastal aquifers, but also along the terminal stretches of the rivers that in turn will let the saline water percolate downward to the unconfined aquifers. This effect, specifically addressed by the GALDIT-SUSI methodology, will be proportional to the increased difference in hydraulic head between the sea and the continental waters (both superficial and underground). On the other hand, coastal erosion, the increased frequency of anomalous waves and storm surges, despite being highly destructive phenomena on a local scale, will not lead to an evident worsening of the vulnerability to SI as they are limited to a thin coastal stripe already classified as highly vulnerable currently, while their negative effects cannot spread much inland. Figure 9 shows a comparison between the distribution of the high and very high SI vulnerability with the study area's land use. The actual distribution of the vulnerability to SI (Figure 9a) affects several typologies of land use, such as urban area, forest, greenhouses and a large portion of arable land, for a total extension of 6000 ha. The situation worsens considering the 2050 scenario (Figure 9b), where the extension of SI could increase, influencing a larger portion of arable land. In detail, approximately 4000 ha more could face salinization issues compared to 2020; in fact, in 2050 the total extension of arable land with high and very high vulnerability will be 10,260 ha. Whereas more than 70% of these territories are currently utilized for agricultural production, the predicted scenario for 2050 could lead to a reduction in suitable land for agriculture [50] due to soil salinization. Here, the shallow depth to the water table together with the depressed morphology, which will enhance water salinization, could directly influence the quality of water taken up by roots, thus limiting plant growth [51].

Actual and Future Seawater Intrusion Comparison
are not always available, thus limiting the applicability of this method to highly characterized coastal aquifers. A future hybridization of this method with numerical models could allow for consideration of a wider number of scenarios and performance of sensitivity analyses on important forcing conditions [54], such as (i) climate change, (ii) increasing groundwater overexploitation, (iii) sea-level rise and (iv) severe drought conditions. Finally, a holistic assessment could also be generated accounting for soil salinization risk, which represents an important issue for agricultural areas [55]. Thus, the methodology applied represents an easy and reliable tool for screening vulnerability to SI, especially for studies at the watershed rather than plot scale, for which necessary data for a numerical model application are often missing.  This study highlights how the use of GIS and remote sensing resources can allow for (i) optimizing the input data resolution, (ii) producing new inputs for future assessment and (iii) increasing the overall result's quality. Nevertheless, it is worth clarifying that this rating methodology cannot overcome the reliability of a process-based numerical model [52,53]. In fact, the rating methods can only generate static snapshots of vulnerability, according to data availability, without highlighting the evolution and trend of the phenomenon over the time. Despite this drawback, the proposed assessment demonstrates how the exploitation of a multidisciplinary approach integrating geological, hydrogeological, geochemical and morphological analysis can return highly reliable results and provide a forecast of the future evolution of SI. The use of high-resolution data contributes to increase the overall assessment reliability, with an accuracy of 7 × 7 m and an error of 5.2 m, calculated by summing the lowest resolution and highest error of the input data (Table 1) with the accuracy of georeferenced points used to create the raster maps.
Furthermore, it must be stressed that GALDIT-SUSI has already been proved to be a reliable assessment tool for different of hydrogeological settings [23]; the future seawater intrusion assessment scenario needs subsidence and groundwater quality trend data that are not always available, thus limiting the applicability of this method to highly characterized coastal aquifers. A future hybridization of this method with numerical models could allow for consideration of a wider number of scenarios and performance of sensitivity analyses on important forcing conditions [54], such as (i) climate change, (ii) increasing groundwater overexploitation, (iii) sea-level rise and (iv) severe drought conditions. Finally, a holistic assessment could also be generated accounting for soil salinization risk, which represents an important issue for agricultural areas [55]. Thus, the methodology applied represents an easy and reliable tool for screening vulnerability to SI, especially for studies at the watershed rather than plot scale, for which necessary data for a numerical model application are often missing.

Conclusions
In this study, the GALDIT-SUSI methodology is combined with remote sensing techniques and applied to a low-lying coastal area of the Volturno River mouth in Campania (Italy) to investigate the vulnerability to water resource salinization. In the study area, different subsidence rates are predicted for 2050, along with a significant increase of the Revelle index. The GALDIT-SUSI mapping showed a higher vulnerability of coastal and inland depressed areas, which will widen appreciably in 2050. Furthermore, the study highlighted the controlling role of lithology both on subsidence and groundwater salinization. The large-scale vertical displacements are linked to the land use. At a closer scale, however, they appear largely controlled by the sediment nature, which determines the compaction rate. The 2050 map shows a worrying scenario depicting the possible degradation of approximately 4000 ha of agricultural land due to salinization. The study demonstrates the suitability of the proposed integration of the GALDIT-SUSI method with remote sensing resources as a screening tool for the assessment of the vulnerability to salinization in coastal areas, emphasizing how classical hydrological studies may be drastically improved by the employment of radar satellite scenes and land use classification from remote sensing for land and water resources management.