Exploring the Relationship between River Discharge and Coastal Erosion: An Integrated Approach Applied to the Pisa Coastal Plain (Italy)

Abstract

Coastal erosion coupled with human-induced pressure has severely affected the coastal areas of the Mediterranean region in the past and continues to do so with increasing intensity today. In this context, the Pisa coastal plain shows a long history of erosion, which started at the beginning of the nineteenth century. In this work, shoreline positions derived from historical maps as well as airborne and DGPS (Differential Global Positioning System) surveys were analyzed in a GIS (Geographic Information System) environment to identify the main changes that have occurred in the last 142 years. These analyses were compared with 100 years of discharge data measured at the S. Giovanni alla Vena gauge to identify a possible correlation between the two sets of information. Finally, Sentinel-2 and Landsat images were studied to identify the dispersion of sediments transported by the Arno River. In particular, we found a minimum of fluvial discharge in the years 1954, 1978, and 2012 corresponding to a peak of erosion, while the reduced erosion rate and the fluvial discharge increased in the years 1928–1944, 1954–1975, and after 2012. The qualitative anticorrelation between discharge and erosion is particularly true if we take into account flood events with a value of discharge greater than 700 m³/s, which are those able to transport suspended sand. The remote sensing analyses of Sentinel-2 images acquired during the floods of 6 February 2019 and 3 December 2019, under the most typical wind and sea state conditions for this area (wind coming from SW and storms coming from W/SW and SW) show that during these events a consistent amount of sediment was transported by the river. However, the majority of these sediments are not deposited along the coastline but are dispersed offshore. Grain-size analyses on the transported sediment show that plumes are formed by coarse-to-medium sand, suitable for coastal nourishment, but the reconstructed sediment dispersion lines show that some sectors of the coastline are constantly in the shade. These areas are the most affected by erosion.

1. Introduction

In the last decades coastal erosion has become one of the main environmental threats worldwide [1,2,3,4]. Approximately 28,000 km² of the global coastline was eroded between 1984 and 2015 and about twice as many as those formed by accumulation processes [5]. According to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) [6] the foreseen future scenario for coastal zones will worsen as a result of the gradual rise in sea level and of the possible increase of extreme events due to global warming linked to the use of fossil fuels and human pressure [6]. Specifically, the Mediterranean region is severely affected by the impact of extreme climatic events (e.g., storm surges) coupled with human-induced pressures (e.g., poorly planned buildings on the coast, dam construction, land use changes inland), resulting in growing vulnerability of the coastal areas [7]. In these areas, the number of coexisting socio-economic activities (urbanization, tourism, natural protected areas) makes it imperative to understand and further monitor coastal dynamics [8,9,10,11]. The change in land cover dynamics is considered one of the most important variables of global change affecting the coastal systems, especially for the effects on river sediment supply [12,13,14]. A decrease in solid load documented for several river-systems of the Mediterranean basin (e.g., Nile, Ebro, Rhône), in particular after the 1970s [8,14,15,16], has been considered one of the dominant causes of coastal erosion, together with subsidence and sea level rise [9,17,18,19].

Despite the importance of a correct estimation of solid load, for many Mediterranean rivers solid load measurements are inadequate, probably because they are cost- and time-consuming. For these reasons the contribution of solid load in countering coast erosion is largely speculative. Even fewer data are available regarding the dynamics of long-shore and off-shore dispersion of sediments carried by rivers, which is another key point for the understanding of coastal erosion dynamics.

Thanks to the availability of a long series of data that make it possible to better understand the current evolution and to predict future trends, in this paper we propose an integrated approach to the study of coastal erosion along the Pisa coastal plain (central Italy). This area has been affected by marked erosive coastal processes since the late nineteenth century (1850–1870) [8,20,21,22,23,24,25,26]. Here, the temporal evolution of erosion and the possible effects of the construction of anthropogenic defenses have been extensively studied [22,27,28], while the primary causes of this phenomenon have been less analyzed.

Our proposed replicable approach—based on a strong integration of data coming from the analysis of shoreline position, river discharge and sediment dispersion, applied to the study of the Pisa coastal plain—could also be applied to other coastal areas affected by coastal erosion. This would provide guidance on some basic elements that need to be taken into account and on appropriate planning choices necessary for the protection of the territory.

2. Study Area

The Pisa coastal plain has been progressively shaped by the Arno river since the Late Holocene [29,30,31,32] (Figure 1a,c–e). In this area SW winds have the highest frequency in each velocity range, while NW winds are frequent but subordinate to the former [33]. Major storm events come from W/SW and SW, even if high energy events can also be related to NW winds [33]. Detailed analyses [34] show that more than 90% of storms (on average 48 events per year) originate in the 220–260° N sector, while less than three events per year generally occur in the 210–180° N sector. Finally, only less than one event per year comes from the other sectors (Figure 1b) [34]. In the coastal sector located in the hydrographic left of the Arno River the littoral drift is mainly oriented towards the south, with the exception of the area between Calambrone and the Scolmatore Channel, where the littoral drift is oriented towards the north. Conversely, in the coastal sector located in the hydrographic right, the littoral drift is oriented towards the north [11,30] (Figure 1c,d). The coast experiences a microtidal regime, where the spring tide is only ca. 30 cm [35].

The coastal plain studied has been affected by a prolonged period of coastal erosion [8,20,21,22,23,24,25,26]. This process is thought to be mainly linked to anthropogenic activities (riverbed dredging, mountain reforestation, diffused riverbed reinforcement, river barriers), which have caused a drastic decrease in the sedimentary load of the Arno river [11,33]. Some key areas for the local touristic economy have been particularly affected by erosion (area A in Figure 1d and area B in Figure 1e). In these two areas, separated from the mouth of the Arno river with its engineering structures built in 1926, the local authorities have applied different approaches in the management of the territory. Specifically, the Area A located in the Migliarino San Rossore Massaciuccoli Natural Park, characterized by a highly natural environment and few anthropogenic structures, was allowed to evolve without the construction of coastal protection structures until 2001. The only exception was five detached breakwaters, built between 1962 and 1968 to defend il Gombo beach. Here, in 1984, after limited results, four out of five breakwaters were elevated and extended, resulting in the formation of tombolos [28]. More recently (2001 and 2003), 10 emerged groins with submerged extensions were built between il Gombo beach and the Arno River mouth [27,28,34]. In 2009 an artificial tombolo was constructed to connect the fourth breakwater and two submerged groins. These last structures were quickly destroyed by coastal erosion [28].

The area B, characterized by the presence of the Marina di Pisa village (founded around 1872), and of tourist facilities and seaside settlements, has been constantly subjected to anthropogenic interventions aimed at mitigating the erosion [35,36]. In this sector, the battle against erosion started in the early twentieth century [28], as documented by the construction of perishable structures designed to protect the beach in front of the recently built settlement of Marina di Pisa [28], which led to the construction of a complex defense system. Before 1928 several short groins were constructed along the first kilometer southward of the river mouth [27,28]. Later, a detached breakwater was constructed and connected to an older groin [28]. Before 1954 another breakwater was constructed, making the first kilometer of coastline totally artificial. A further improvement in defense construction is observed in the following years. Specifically, 10 detached breakwaters were constructed, some of which were connected to the coast with groins [27,34,37] (Figure 1d).

In this way the shoreline position was fixed but erosion continued on the submerged profile [28]. Up until now, the structures have been modified in different ways and recently four of the ten breakwaters have been converted into submerged structures. In the last years, without a general development project, the construction of shoreline protection has undergone a southward migration towards the village of Tirrenia, in order to maintain the touristic potential of the beach. However, the touristic value of this coast has been strongly reduced, as documented by [28]. In the last 20 years this area has received several artificial nourishments [38,39]. This different evolution of sectors A and B has determined the asymmetrical development of the Arno river delta [27]. Coastal erosion is still active in both sectors, despite the major efforts to counter this process [22,24,34].

3. Materials and Methods

The changes in the coastline position of the Pisa coastal plain were studied for a chronological interval of 142 years (data from 1878 to 2020) by using a heterogeneous database consisting of topographic maps, airborne images and DGPS (Differential Global Position System) measurements (

Table 1. Census of shoreline data used in this study.

Year	Ownership Organization	Data Type	Properties	Source
1878	I.G.M	Cartography	Scale 1:25,000	This work
1907	I.G.M	Cartography	Scale 1:25,000	This work
1928	I.G.M	Cartography	Scale 1:25,000	This work
1944	R.A.F.	Aerial photographs	Black and white film	[22]
1954	Tuscany Region	Aerial photographs	Black and white film	[22]
1965	I.G.M	Aerial photographs	Black and white film	[22]
1975	Tuscany Region	Aerial photographs	Black and white film	[22]
1978	Tuscany Region	Aerial photographs	Color film	[22]
1982	Tuscany Region	Aerial photographs	Color film	[22]
1982	I.G.M	Aerial photographs	Black and white film	[22]
1986	I.G.M	Aerial photographs	Black and white film	[22]
1990	I.G.M	Aerial photographs	Infrared	[22]
1996	Tuscany Region	Aerial photographs	Black and white film	[22]
2003	Orthophoto	Aerial photographs	Black and white film	[22]
2004	Pisa Province	Aerial photographs	Color film	[22]
2008		DGPS measurements		[34]
2009		DGPS measurements		[34]
2010		DGPS measurements		[34]
2011		DGPS measurements		[34]
2012		DGPS measurements		[34]
2013		DGPS measurements		[34]
2014		DGPS measurements		[34]
2015		DGPS measurements		This work
2020		DGPS measurements		This work

). The data were managed using the open source QGIS 3.10 software. A linear shapefile with the position of the shoreline (derived by manual operator delimitation and digitalization) was created for each analyzed year until 2008. The data of the last 12 years were derived from dedicated DGPS surveys performed following the methods proposed by [37]. The differences in the area between the shorelines acquired in different years were calculated in terms of annual gains and losses.

All the data available were processed using GIS tools in order to identify the coastal areas characterized by the advancing shoreline (gain of territory) and those characterized by the retreating shoreline (loss of territory) for the entire period considered. Specifically, a comparison of the shoreline in the successive years made it possible to identify coastal sectors that had never experienced erosion and others that had never experienced advancement in the last 142 years. The area constantly in advancement and the area constantly in erosion were quantified and plotted on a histogram. The sectors that experienced a mixed behavior were subdivided into areas mainly in erosion and areas mainly in advancement, according to the prevailing behavior.

One of the longest and most complete sets of daily discharge data of the Arno River was also analyzed. The daily discharge data were recorded at the S. Giovanni alla Vena gauge (data from 1924 to 2020; https://www.sir.toscana.it/consistenza-rete), 35 km inland from the river mouth (Figure 1b). The discharge value measured at this gauge is significant because it could be considered close to the discharge value at the river mouth, considering that there are no important tributaries downstream from the gauge.

The trend of the daily discharge data since 1924 has been analyzed to identify the main variations of that period.

The return times for flood events were calculated as opposites of probability:

$$p=\frac{1}{T},$$

(1)

where $p$ is the probability and $T$ the return time. This equation can be applied for frequent events, but it can be difficult to determine the return time of an extreme event if the occurrence of these events is too small in the series history [40]. For this statistical reason, we investigated the role of flow events with a maximum return time of about eight years.

In order to identify possible correlations, the set of discharge data were analyzed and compared with the erosion trends, as documented by the study of the shoreline positions. Since a rigorous statistical treatment of the data was not possible because of the limited number of shoreline measurements (despite the long period considered), a qualitative comparison was made between the trend of erosion and the trend of discharge.

The trend of the selected flow range was compared with erosion over time by using a moving average calculated over ten years. The use of the moving average to understand the trend of river discharge is a validated method, as demonstrated by several authors (e.g., [41,42,43]). To support the analyses, we investigated Sentinel-2 images acquired during two recent flood events occurred on 3 February 2019 and on 3 December 2019, in the most frequent wind and sea state (winds and waves from NW and SW, respectively; Figure 1d) in order to identify the presence of sediments and their longshore and offshore dispersion. The selected images were analyzed by means of the tone mapping technique using Photomatix, a program that allows the user to choose among different HDR profiles. After numerous tests, the use of the profile called “Painterly” was chosen. This profile gives the image a very photographic look and a well-exposed dynamic range. The image acquires an aspect that is very similar to the one perceived by the human eye and is often referred to “as the eye sees” in the field of photography. After the application of this profile, the parameters of gamma, luminance and contrast varied, so as to make more readable all the information contained in the images. This method makes it possible to qualitatively identify the area subject to solid transport (according to the different shades and gradations of color). Moreover, the method allows the identification of flow lines of solid transport direction that show the dispersion of sediments out of the mouth towards the coast.

In addition to the study of these two recent events, in order to extend the analyses of the plumes generated by the Arno River, we considered 151 frames from the Landsat (5, 7 and 8) and Sentinel-2 satellites. All the images were filtered at a threshold of less than 20% cloudiness during the download procedures and were then checked manually. Specifically, 51 Sentinel-2 frames were chosen from 2015 to 2020 and 101 Landsat frames from 1984 to 2020. In particular, in this second case 68 frames were from Landsat 5, 14 from Landsat 7, and 19 from Landsat 8. The red band was analyzed for both satellites, corresponding to bands number 4 for Sentinel-2 and Landsat 8 and to number 3 for Landsats 5 and 7. A good correlation was found between the red band (630–690 nm) reflectance values and low-to-moderate turbidity [44,45,46,47]. The images were resampled on the same matrix. Two grids were extracted. The first provides the maximum value recorded in the analyzed images for each pixel, while the second provides the mean value calculated using the analyzed images for each pixel.

In an attempt to understand whether the types of sediment transported by the river were suitable to counter coastal erosion, we performed grain size analyses on sediments sampled on the riverbed near the city of Pisa (where the sediments left by the river were easy to sample) and on the beach (area C, Figure 1) immediately after a flooding event that occurred on 3 December 2019, which was one of most representative in this area, in terms of both wave direction and fluvial discharge. The samples were dried in an oven at 40 °C for 48 h. Each sample was then weighed and analyzed by means of standard sieving with mesh size 1 mm, 0.5 mm, 0.315 mm, 0.250 mm, 0.125 mm and 0.063 mm. After weighing the fraction left in each sieve a grain size curve was constructed using GRADISTAT [48]. Finally, to check the effective ability of solid load to counter coastal erosion, we undertook a dedicated DGPS survey in a small stretch of coast in the southern sector (part of area C) about one month after the flood event of 3 February 2019.

4. Results

4.1. Shoreline GIS Analysis

The analysis of shoreline positions (Figure 2a–d) shows that more than 2.5 km² of the studied area experienced constant erosion, while only 0.5 km² experienced constant advances (Figure 3). The areas most affected by erosion are located in the hydrographic right (Area A, Figure 2c and the area northward of the Morto Nuovo River, Figure 2a). In the hydrographic left (where erosion has been countered by several engineering interventions such as groins and breakwaters [22,27,28,49]) the total area eroded is less in extent compared with Area A, but the erosion is persistent in the area marked B in Figure 2c. Coastal progradation has been documented mainly in the C area of the southern sector (Figure 2d), while minor local advances have also been identified in the northern sector immediately south of the mouth of the Morto Nuovo River (Figure 2b). From a chronological point of view, the highest values of erosion are documented around the years 1954, 1978 and 2012 (Figure 4), while the lowest values of erosion are reported for the periods 1881–1944, 1965–1975, and 2013–2020. It is worth noting that the erosion rates for the period before 1944 are similar to those recorded in the last 10 years (Figure 4).

4.2. Discharge Data Analysis

Discharge analysis (Figure 4) shows the highest values until 1944, followed by a reduction in the period at the end of the 1950s and around the year 1980. These two periods are separated by a peak in the values of discharge. An increase in discharge values, which are still well below the values recorded before 1944, is documented from the end of the 1980s to the year 2005. Despite some fluctuations, a further increase in the trend is documented after this year. We are aware that the lack of data in the 1950s and 1970s may have affected the results, even if a clear decreasing trend, probably not so marked, can be identified from the available data. Return time analysis of the historical series of S. Giovanni alla Vena (Figure 5) with a 100 m³/s discretization shows an exponential trend with a one-year return time corresponding to a flow rate of 700/800 m³/s.

4.3. Remote Sensing Analysis

The sentinel-2 image recorded at 10:24 a.m. on 3 December 2019 shows the Arno River sediment dispersion along the coast during a flood event characterized by a peak discharge of 902.81 m³/s at the S. Giovanni alla Vena gauge (return time of 1.5 years, Figure 5). This event occurred with a much more frequent wave direction varying from 240°N to 180°N (Tuscany Region, wave buoy of Gorgona: https://www.cfr.toscana.it/index.php?IDS=42&IDSS=282). The image processed by tone mapping (Figure 6) shows a wide plume formed by sediments transported by the river during this event and dispersed at the river mouth. Plume n. 1 in Figure 6 covers a total area of about 100 km². A minor plume (identified as 1A in Figure 6) covering a total area of about 10 km² is also well visible in the image. The dispersion of sediments shows that a large amount of transported sediment was dispersed offshore. In agreement with the wind and wave directions, the plume shows a deviation towards NNW, together with the other minor plumes of the Scolmatore channel located further south. The reconstruction of the flow line allows the identification of two shadow areas in sediment dispersion: one located in the northern sector between the Lame della Gelosia and il Gombo (area A in Figure 6), and the other one situated in the southern sector in front of the village of Marina di Pisa (area B in Figure 6).

The Sentinel-2 image recorded at 10:22 a.m. on 6 February 2019 shows the final phase of the flood event of 3 February, characterized by a wave direction from 300° N (which is the second in occurrence), and a fluvial discharge of 1131 m³/s at the S. Giovanni alla Vena gauge (return time of 3 years, Figure 7). In this case a significant volume of sediments was dispersed offshore, while the entire solid longshore transport faced south (Figure 6). Apart from the first kilometer to the north of the river mouth, the northern sector does not receive any sediments. In the southern sector, the shadow area in front of Marina di Pisa is still visible; instead, according to wave direction analyses, a part of the solid load reaches the southernmost area in correspondence to the village of Tirrenia and the borough of Calambrone (area C, Figure 1d).

A more generalized analysis of sediment dispersion of the Arno River along the coast was performed considering several (Landsat and Sentinel-2) images for the period 1984–2020. Specifically, Figure 8 shows the result of processing the 50 Sentinel-2 and the 101 Landsat 5, 7 and eight images. The analysis was performed using the red band. In the first processing data the final images (Figure 8a,c) were obtained considering in each pixel the mean value recorded in all the analyzed images. Specifically, the mean value on the red band of the 50 analyzed Sentinel-2 images is recorded in each pixel of Figure 8a, while the mean value of the 101 analyzed Landsat images is recorded in each pixel of Figure 8c. In the second processing data the final images are obtained considering in each pixel the maximum value recorded on the red band in all the analyzed images (Figure 8b for the Sentinel-2 Images and 8d for the Landsat images.). All the images obtained show that the area to the north of the river mouth is characterized by a greater reflection on the red band than the southern area. This suggests a different distribution of the sediment between the northern part and the southern part. In particular, the southern area near the mouth of the river has the lowest reflection values for the entire stretch of coast studied. Moreover, the shadow areas described in the events previously described (Figure 6 and Figure 7) in the northern and southern sectors are confirmed also in Figure 8b,d, which consider the maximum pixel value. The images based on the mean values (Figure 8a,c) allow identification of the main distribution of the sediment in the period studied (considering only the mean value of the differences among the various sectors are minimized). The images based on the maximum values (Figure 8b,d) accentuate the differences between the areas that received a greater supply of sediment and those with the opposite behavior.

4.4. Post-Flood Field Investigations

In C area between the village of Tirrenia and the borough of Calambrone, new DGPS shoreline measurements acquired before and about one month after the flood of 3 February 2019 document an advance of the coastline (Figure 9). More precisely, on a linear stretch of coast of about 6.5 km, there was an increase of the territory of about 16,400 m² following this event, which occurred under more favorable conditions of the wind and sea states (wind and wave directions from NW) for nourishment of the southern sector. Grain size distribution of sediment acquired in two points (the first on the sediments left by the river on the city center embankment and the second on the beach of the area C: Figure 1b) shows an evident analogy between the samples and a predominance of coarse sand (50% ca) and medium sand (45% ca) in both cases (Figure 10), excluding the plumes identified by remote sensing that were formed by fine material unsuitable for beach nourishment.

5. Discussions

The analysis of the coastline position during the last 142 years represents the most comprehensive study ever performed on the history of coastline evolution of the Pisa coastal plain. Our results actually expand the work of [8], both in terms of time ([8] considers only the period between 1944 and 2015) and resolution, since they mainly consider satellite images. Analysis of the data shows a predominance of areas permanently affected by erosion, corresponding to 60% of the total area investigated. Following a period of progradation due to the phases of Arno delta construction beginning at about 3000 ka BP [29,30,31,32,50], different parts of the Pisa plain experienced marked coastal erosion, which started at the end of the 19th century and increased after the construction of the river mouth jetty, especially on the hydrographic right [20,22,28,34].

However, with the exception of area A (Figure 1), which was particularly affected by erosion, the overall value of erosion remained low until the 1950s when there was a rapid documented increase of the process. Erosion was particularly severe at the end of the 1980s, probably caused, among other things, by the effects of dredging/damming [51,52]. The following period was characterized by an increase in erosion around 2010, while a reduction in the erosion rate has been documented in the last eight years (since 2012). The period most affected by erosion was around the 1980s, while the areas that most experienced this phenomenon are the area northwards of the Morto Nuovo River, area A where the Arno mouth jetty limits sediment flow that evolved naturally until 2001, and area B area, where several human interventions designed to mitigate the erosion effects have been constructed since the beginning of the twentieth century [28,36].

In the absence of precise modern measures of solid load, [53,54] provide an estimate of 1.5 t/yrs. In addition, only eight sporadic measurements were made far from the river mouth by [55]. It is possible to indirectly define the role of solid load in countering erosion, in an attempt to create a correlation with fluvial discharge [8,14]. The sporadic data obtained by [55] made it possible to highlight the presence of sand in the suspended sediments transported by the river, starting from the discharge value of about 500 m³/s at the Nave di Rosano gauge (Figure 1). Figure 11a shows the sections of S. Giovanni alla Vena and of Nave di Rosano, while Figure 11b shows the relationship between the discharges of these two stations. The linear regression of the two datasets is represented by the following equation:

$$y=0.60x-1.86,$$

(2)

where $y$ is the discharge of Nave di Rosano and $x$ the discharge of the S. Giovanni alla Vena. The equation highlights that the discharge at Nave di Rosano is about 60% the discharge at S. Giovanni alla Vena. R² of the linear regression is 0.79 and RMSE (Root Mean Square Error) is 40.25 m³/s. By using the outflow scale provided by the Regional Hydrologic Service of the two investigated stations, we can derive velocity as a function of discharge (Figure 11c). A discharge of 500 m³/s at Nave di Rosano corresponds to a flow rate of ca 1 m/s. The same flow rate at S. Giovanni alla Vena corresponds to a little less than 700 m³/s discharge.

For this reason, it may be hypothesized that values of discharge starting from 700 m³/s at S. Giovanni alla Vena are the most suitable to counter coastal erosion transported by suspended sand. These data are consistent with those observed in the Tiber River [56] in a climatic and environmental setting not very different from the Arno River catchment.

The role of solid load in countering the coastal erosion of this territory is documented by the qualitative anticorrelation between fluvial discharge and erosion rate (Figure 12). In particular, Figure 12 shows a minimum of fluvial discharge during the years 1954, 1978, 2012 corresponding to a peak of erosion, while in the years 1928–1944, 1954–1975, and after 2012 the erosion rate (despite some fluctuations) diminished, and the fluvial discharge increased. In the period 1960–2012 the river discharge was significantly low. This was particularly true for events >700 m³/s (purple bars in Figure 12). It is logical to presume that the same occurred with the transported sediment load. It is perhaps for this reason that several authors (e.g., [11,22,57]), on the basis of the data available until 2010, considered the Pisa coastal plain to have evolved as a relict beach. However, this general trend seems to have changed over the last decade, with a slight decrease in erosion possibly caused by a slight increase of the solid load related to an increase in events greater than 700 m³/s. In more recent years the data concerning the solid transport of sand in suspension starting from 700 m³/s have been more reliable, and large environmental changes (in sections of the riverbed, availability of sediments resulting from different soil use in the catchment basin, etc.) can be excluded. Moreover, the quality of the shoreline data in the last 12 years (through high frequency DGPS acquisition) has made it possible to better understand the relationship between river discharge and the trend of erosion.

The study of Sentinel-2 images related to the flood event of 3 December 2019 with fluvial discharge of ca 900 m³/s (a little higher than the threshold of 700 m³/s identified for the transport of sand), allowed us to identify the presence of a large plume of sediments in front of the mouth of the Arno River. Moreover, grain size analysis of the transported sediments highlighted the presence of well-sorted medium sand (>0.315 mm), perfectly consistent with the sediments characterizing different parts of the Pisa coastal plain [37] and suitable for coast nourishment. However, the sediments forming the plume were largely dispersed offshore.

Analysis of the images presented in Figure 6, Figure 7 and Figure 8 also made it possible to infer that strong erosion still characterizing some sectors of the Pisa coastal plain is related to the mode of longshore dispersion of the sediments that reached the coastline. Indeed, the analyzed flood event can be considered representative of most of the floods characterizing the Pisa coastal plain, which took place under the most common wind and sea state (more than 90% of events occurred with wind and wave directions from W/SW). The reconstruction of the flow lines (Figure 6) shows that areas A and B were not enhanced by the dispersion of sediments transported during the floods that occurred under these wind and sea state conditions.

The study of the events of 3 February (higher in discharge: ca 1100 m³/s), which arose in less frequency but had more suitable sea state conditions in the southern sector (wind and wave directions from NW) underlines the presence of the same shadow areas. The most favorable condition for nourishing the northern beaches is the occurrence of floods with winds blowing and waves coming from SW, while the most favorable condition for nourishing the southern sector is the occurrence of floods with waves from the W and/or NW (Figure 6). Thanks to these studies it is possible to hypothesize that coastal sectors that are not fed by sediments during these events generally lack sedimentary supply. These interpretations have been also confirmed by analyses on several Landsat and Sentinel-2 images showing the process of sediment dispersion over the chronological interval 1984–2020. In all the cases analyzed, many Arno River sediments do not reach the coastline but are dispersed offshore. This aspect is likely to have more relevant effects on coastal erosion than the decrease in the solid load of sediments transported by the Arno River, as hypothesized by several authors [11,22,30]. However, our work also suggests that sediment load may have had an opposite trend in the last decade, as shown by the analysis of fluvial discharge analysis, with a slight increase in the number of events higher than 700 m³/s, suitable for sand transport.

As highlighted by Sentinel-2 analyses, the shape of the jetty at the river mouth also played a key role in the dispersion of sediments offshore and in the creation of shadow areas longshore. As documented in Figure 6 and Figure 7, some sediments reach the northern sector passing through an opening in the jetty (see plume 1a, Figure 6), but this intervention is still insufficient to limit the formation of the shadow zone in the northern sector. In this respect, it may be useful to seriously consider the maintenance/change of shape of this engineering structure, built at the beginning of the 1900s, which was supposed to avoid the silting up of the river mouth and to reduce the flood hazard that is creating significant coastal erosion problems. In more recent years, a similar crucial role has been played by the jetty at the Morto Nuovo River. It is important to reflect on the extensive use of expensive and impacting engineering defenses in this area, which do not always give the expected protection results. In some cases the defenses have been destroyed by coastal erosion as documented by [28], and in other cases the defenses have shifted the coastal erosion from one coastal sector to another (e.g., areas of the village of Marina di Pisa, where erosion and coastal defenses are moving southwards and seaward acting on the submerged profile) [28,34].

The possibility for the solid load to reach the coastline to counter erosion was confirmed by using a DGPS to measure a small stretch of coast in the southern sector in front of the village of Tirrenia (area C, which exhibits a general trend of stability, reached by sediments during the described events). We documented an advance by comparing the shoreline position obtained by DGPS before the event of 3 February 2019 with the shoreline measured about one month after this event. Since grain size analysis of the sediments of the riverbed and of the southern coastal area are consistent and correspond to medium sand, we can deduce that part of the southern coast receives sediments from the Arno River. Furthermore, the grain size range of the transported sediments is useful to counter erosion. On the contrary, areas A and B are affected by erosion because they are not influenced by the dispersion of the sediments of the plume under different marine weather conditions.

6. Conclusions

Coastal erosion in the Pisa coastal plain at the end of the nineteenth century increased according to our review, in particular after the 1950s, at the end of the 1970s and around 2012. After 2012 there is a slight decrease in the erosion rate, which reaches the same values observed at the beginning of the nineteenth century. On the basis of the erosion rate and the remote sensing analyses on sediment dispersion, we can state that the particular shape of the jetty and of other engineering defenses at the mouth of the Arno River favors the dispersion offshore of the sediments carried by the watercourse and determines the formation of shaded areas in longshore sediment distribution. Despite the construction of several engineering defenses, particularly active erosional trends still appear in two coastal areas (in front of the small town of Marina di Pisa in the southern sector and between the locations of Lame della Gelosia and il Gombo in the northern sector). Erosion is particularly evident in the first location, while it has been countered by numerous, repeated, impacting and expensive engineering interventions for coastal protection in the second. The fluvial discharge data of the Arno River highlight a discharge decrease that roughly matches the periods that correspond to an increase in erosion. On the other hand, the increase in the flow rate recorded in the last decade can be considered a key element for a reduction in the overall erosion rate recorded in this last period. This correlation is particularly true if we take into account flood events with a value of discharge greater than 700 m³/s, which are those able to transport suspended sand.

The amount of sediment transported by the watercourse has increased in recent years and could certainly counter erosion more effectively if it were possible to reduce the amount of sediment dispersed offshore. Furthermore, it would be important that the sediments were distributed more uniformly longshore, so as to avoid the formation of shadow areas in the distribution of the sedimentary load transported by the watercourse. Although modest, the changes of the erosive trend that have been observed in the last eight years represent an important signal for the development of this territory. It will be necessary to monitor this coastal stretch with DGPS at high temporal resolution in order to understand whether this documented decrease will continue over time.

Future sporadic measurements of the solid load conducted near the river mouth during flood events could allow for calibration of the satellite observations (plume amplitude versus measured solid load data) to derive the solid load data directly from the study of satellite images, with considerable saving of time and money. A similar integrated approach could be easily used in other contexts affected by coastal erosion, where a holistic approach of this type could help identify unclear causes of erosion and support future development of these sensitive areas.