Evidence and Implications of Hydrological and Climatic Change in the Reno and Lamone River Basins and Related Coastal Areas (Emilia-Romagna, Northern Italy) over the Last Century

Abstract

Climate change and human activities have consequences on coastal areas as they affect hydrological processes in the related river basins. The riverine sediment supply to the beaches of the Emilia-Romagna coast, a highly urbanized area with high economic and naturalistic value, has been heavily impacted by human activities throughout the catchment, reducing solid transport to the coast and increasing the threat of coastal erosion and flooding. Despite the introduction of safeguard policies in the early 1980s and the consequent stoppage of such activities, the expected return in solid transport has not yet been reflected at the coast. To better understand the various processes acting at the river basin scale, we utilized empirical mode decomposition to analyze the variability in different parameters (river discharge, rainfall, air temperature, and sea level) from the headwaters to the coast of the Reno and Lamone rivers over the last century. The anthropogenic footprint, linked to the large-scale dimming/brightening phenomenon, is visible in the long-term trends. Moreover, natural signals with variable periodicity are evident and partially correlated with two major climate modes (North Atlantic Oscillation and Atlantic Multidecadal Oscillation). The coupled interactions among these processes, combined with the changes in land use and evapotranspiration during the last century, have resulted in the prolonged scarcity of river sediment supply and a long-term trend of erosion of the coastal area.

1. Introduction

Riverine sediment discharge, the foremost feeding source for low-elevation coastal zones (LECZ) [1,2], has declined drastically worldwide over the last century (e.g., [3,4,5,6,7,8]) and will likely continue to decline [9,10]. The amount and role of sediment supplies at the coast varies with the local climate, geological setting, and degree of human impact within each catchment area [11,12,13]. LECZs worldwide have reached a critical tipping point, shifting from a natural pseudo-stable state since the Holocene to a state controlled by human dynamics [14] owing to land use modifications throughout their catchment areas [15,16]. Furthermore, rising sea levels and other climate-driven effects are expected to exacerbate such impacts on LECZs [17]. Despite climate-driven modifications remaining influential, they may locally be of secondary importance [18] since human effects may more heavily interact with hydrological processes. Anthropogenic modifications over catchment areas, together with the hydrological and climatic changes, and their related connection with atmospheric dynamics, should be thus taken into account in order to better understand which processes are most likely to drive the evolution of coastal areas.

The Emilia-Romagna (ER) coast is a 130 km long LECZ south of the Po River Delta, Northern Italy, facing the northern Adriatic Sea (Figure 1a), which has undergone severe anthropization since the 1950s [19]. Increased subsidence due to underground fluid exploitation [20,21,22,23], land reclamation, and the widespread use of coastal defense structures [24,25,26] have affected the ER coastal plain, making it sensitive to floods [27]. Furthermore, the heavy, anthropogenically driven reduction in sediment delivery by rivers (mainly due to sediment digging and river regulation), mostly between the 1940s and 1980s, has led to generalized erosion and shoreline retreat [28,29,30]. Although safeguard policies for riverbed excavation, introduced by the ER regional government during the early 1980s, were expected to gradually reverse this, only the amount of suspended material has slowly increased, whereas fluvial bedload has decreased [31,32]. Granulometric analysis on samples collected since the 1970s in the submerged beach of the ER coast, at or below 3 m depth, has shown a generalized overall decrease in grain size until the present day [33,34,35,36]. Indeed, a recovery in bedload sediment supply has not yet occurred since the introduction of the safeguard policies; thus, most of the ER shoreline has become threatened by erosion and flooding [34]. Previous analyses have been conducted over the regional scale in order to better elucidate the hydrological, environmental, and climatic evolution of the area, finding a significant decrease in river discharge [37] and indicating that sediments are currently not being transferred downstream as they were previously. While a clear increase in air temperature has also been observed for the last 70 years (e.g., [38]), trends in rainfall have been found to be statistically significant only at certain locations [39,40,41], highlighting a great spatial variability for this parameter. This effect has primarily been linked to the frequencies of diverse weather regimes which, locally, have led to differentiation in both total annual rainfall amount and its variability through time [42]. Generally speaking, significant/non-significant decreases in rainfall are a matter of debate on both the Italian [43,44,45,46] and Mediterranean scales [47,48], for which the contrast has been mainly attributed to the length of the time interval considered in the analysis [45].

In this study, the Reno and Lamone river basins (Figure 1a), which feed an approximately 20 km stretch of the ER coast between the Porto Garibaldi and Porto Corsini harbors (Figure 1b), have been considered. The Reno River, with a length of 210 km and a basin area of approximately 4630 km², is the main river of the ER region after the Po River, and the sixth nationwide in terms of basin size. It drains both mountains and plains and is joined by many tributaries from the Apennines and artificial channels in the plains (Figure 1c). The basin of the Lamone River (approximately 520 km², Figure 1a) is much smaller than the Reno River basin; nevertheless, its sediment supply also provides an input for the coastal budget of this area. Coastal prograding at the mouth of the Reno River during the early 20th century was previously associated with a sustained sediment supply [33], whereas from the late 1930s onward, progressive erosion has occurred [34]. The coastal areas of these river basins and adjacent wetlands (the Comacchio Valleys and the Bellocchio marsh; Figure 1b), represent a significant economic and naturalistic resource (Ramsar Convention site); however, these areas are among the most threatened by shoreline retreat and flooding along the Mediterranean and European coasts [49]. Furthermore, these coasts are significantly impacted by economic activities that cause land subsidence, such as gas exploitation by ENI (Ente Nazionale Idrocarburi) in the “Dosso degli Angeli” reservoir (Figure 1b) [50].

In this paper, the linear and nonlinear evolution of environmental-climatic parameters for both river basins have been analyzed by considering river discharge, rainfall, and air temperature datasets, representative of the river basins over the last century. The main goal of the study is to better understand the implications of the hydrological and climatic variations occurring over the river basins and at the coast. According to Preciso et al. [37], the Reno River catchment no longer represents the primary sediment source for the coast, primarily because of reforestation within the catchment since the 1950s. Moreover, in several Mediterranean locations, changes in land use, together with increased evapotranspiration and soil erosion control measures, have reduced runoff, and thus, sediment supply [51,52,53,54]. Due to this, land use changes and evapotranspiration datasets were also taken into account in our analyses in order to obtain a more holistic view of the changes occurring throughout the catchment areas. Moreover, the local sea level, here affected by both natural and anthropogenic components (see [55,56] for details), was considered, as it is a critical factor in coastal evolution.

In detail, this study addresses the following questions: (1) How have anthropogenic and natural signals acted on hydrological processes at the river-basin scale, from the headwater to the coast, over the last century? (2) Why, despite safeguards for river basin management, has the solid river supply to coastal sites not yet recovered? (3) How have climate change and anthropogenic footprints affected sea-level changes in the studied LECZ?

2. Materials and Methods

2.1. Data

The average monthly river discharge, rainfall, and air temperature data were obtained from the Dext3r-SIMC platform [57] and the annals of the National Hydrological Service [58]. Only long-term data within both the river basins were obtained, in order to consider a wider, historical view of the surface processes. Furthermore, liquid transport time series were reconstructed from three selected river gauges (Figure 1c); however, these datasets contained gaps (Figure 2). To rebuild a representative time series, thereby increasing the signal-to-noise ratio, data from the twelve longest and more complete individual rain gauges (out of a total of 92 stations available throughout the two basins) were considered and stacked. Although these datasets were affected by a common gap in the year 1924, minor gaps in each time series were filled using records from the closest station (excluded from the analysis) within 10 km after careful comparison of the overlapping data. To further examine rainfall variability, the number of wet days (≥1 mm/day) were also considered. The minimum and maximum temperatures were used to produce a monthly climatological assessment of the river basins. Additionally, a stacked curve was computed to achieve a basin-representative climate time series that considered data from nine stations (Figure 1c).

The relative sea level data (RSL, which measures the sea level with respect to a benchmark on land) were obtained from the Marina di Ravenna tide gauge (RTG; Figure 1b), which is currently managed by ISPRA (Istituto Superiore per la Protezione e la Ricerca Ambientale). The gauge has recorded data since 1873, with a gap from 1922 to 1933 and discontinuities since 2016 (Figure 2); hence, data after 2015 were discarded. The RTG time series was assembled using data from various sources (1873–1979 and 2001–2016 from the Permanent Service for Mean Sea Level (PSMSL) data bank [59,60]; 1971–2001 from Romagnoli et al. [61], and homogenized considering data overlaps and following Bruni et al. [56]).

River discharge, rainfall, temperature, and RTG time series were compared with three climate indices, namely the North Atlantic Oscillation (NAO), Atlantic Multidecadal Oscillation (AMO), and Western Mediterranean Oscillation (WeMO). The NAO is among the most significant modes of climate variability in the Northern Hemisphere, with seasonally fluctuating manifestations, especially during the boreal winter [62,63]. The winter (December–March) station-based NAO dataset was used; this represents the normalized sea-level pressure between Lisbon (Portugal) and Stykkishólmur (Iceland) since 1864, by the National Centers for Environmental Prediction/National Center for Atmospheric Research [64], based on Hurrell et al. [63]. The AMO [65] is another prominent climate variability mode in the Northern Hemisphere, defined by variability in sea surface temperature (SST) in the North Atlantic Ocean and characterized by a dominant periodicity of approximately 60–70 years. The AMO is the area-weighted average SST (between 0° and 70° N) with the linear trend removed. The annual AMO index, available since 1856, was provided by the National Oceanic and Atmospheric Administration/Physical Sciences Division [66]. The WeMO index [67] is defined as a dipole structure formed (in its positive phase) by an anticyclone over the Azores and a depression over north-western Italy; its index is the result of the difference between the standardized values of surface atmospheric pressure in San Fernando (Spain) and Padua (Italy). In this study, we also accounted for the WeMO since it is considered by some authors more significant than the NAO to explain rainfall anomalies over the north-western Mediterranean (e.g., [68] and references within).

The Standardized Precipitation Evapotranspiration Index (SPEI) [69] was used for long-term drought analysis (1921–2019). The SPEI allows for the identification of drought severity using its length and intensity at various timescales, considering the effects of both temperature and precipitation. This index combines the strength of the multi-temporal Standardized Precipitation Index [70] and the sensitivity of the Palmer Drought Severity Index [71]. The Global SPEI database [72] used in this study covers the globe at a spatial resolution of 0.5 degrees and offers timescales from 1 to 48 months. The potential evapotranspiration in the SPEI database is based on the FAO-56 Penman-Monteith method [73]. Despite its coarser spatial resolution with respect to local datasets, the SPEI has been adopted as it has been thoroughly validated by several authors (e.g., [74,75,76,77,78,79,80,81]) to reveal moisture anomalies for environmental and agricultural applications. Since our study primarily considers the long-term (decadal to centennial) alternation of drought and wet conditions rather than short-term (monthly) events, the 48 month-average time series have been used. The Bologna grid cell was considered as representative of the Lamone and Reno river basins because the variability between adjacent cells at this spatial scale was negligible.

Snow cover percentage over the period 1950–2022 was obtained from the ERA5-Land reanalysis dataset [82,83,84] and limited to the area under analysis. This parameter represents the fraction (0–100%) of the grid cell occupied by snow and is provided with a spatial resolution of 0.1 × 0.1 degrees.

Land use was evaluated using land use maps for 1976, 1994, 2003, 2008, 2014, and 2017, produced by the GIS and Statistic Service of the ER Region [85]. These maps were developed utilizing aerial photography and are characterized by their minimal detectable area: 0.375 (1976), 1.56 (1994 and 2003), and 0.16 (2008, 2014, and 2017) hectares. The land use for these products was mapped using the Corine Land Cover legend until level three.

2.2. Methods

All datasets in this study were homogenized into monthly means. The time series were plotted using a 1-year moving average low-pass filter to highlight the interannual variability. Modified non-parametric Mann-Kendall (MK) tests [86,87,88] were performed to assess temporal trends (statistically significant at the 95% confidence interval). MK analysis was implemented using the Theil-Sen estimator [89,90] to evaluate linear trends. This method is insensitive to outliers and more accurate than simple linear regression [91]. To account for serial autocorrelation and properly estimate trend-related errors, the procedure described in Zervas [92] has been followed.

Considering the stacked rainfall time series and the monthly sum of wet days, the Simple Daily Intensity Index (SDII), which is the ratio of total rainfall to the number of wet days, was computed based on Frich et al. [93]. Only the SDII values above the 95th and 99th percentiles were considered.

For RSL analysis, several authors have addressed the removal or assessment of vertical crustal movements (VLM) from tide gauge data [55,94,95,96,97]. Optimally, VLM contributions can be distinguished by Global Navigation Satellite System monitoring of tide gauge data [98]; however, this began in the late 1990s [99,100]. To remove the VLM components, the RTG time series was subtracted from the Trieste Tide Gauge (TTG; location in Figure 1a) dataset (data 1875–2021 from the PSMSL archive); from this, a locally weighted polynomial regression was created and subsequently removed from the original RTG datasets. This procedure relies on the assumption that TTG location is considered a relatively stable area in terms of vertical land movements, with a weak uplift in the range of a few tenths of mm·year⁻¹ or negligible (e.g., [101,102,103,104,105,106]).

To assess nonlinear and non-stationary signals, empirical mode decomposition (EMD) [107] was applied to decompose monthly mean river discharge, rainfall, air temperature, and sea level into finite empirical orthogonal intrinsic mode functions (IMFs). Each IMF describes cyclic variations, though not necessarily constant phases and amplitudes, representative of the oscillating mode from the highest to the lowest frequency. The lasting, non-oscillatory mode is the residual (RES), usually associated with the long-term trend of the signal, which is typically monotonic and provides information about the reliability of the linear model for nonperiodic components. Previous studies have employed EMD to correlate nonlinear variations among different phenomena [108,109,110,111,112,113,114]; moreover, it is well suited for analyzing non-stationary time series.

The non-parametric Kendall rank correlation was employed to examine statistical relationships for NAO, AMO, and WeMO with our time series. This method has significantly smaller gross error sensitivity and asymptotic variance than Spearman rho or Pearson r and provides significantly more accurate p-values [115]. To evaluate autocorrelation, the time series were linearly detrended, and their effective sample size, based on the lag-1 autocorrelation coefficient, was considered [116,117]. The statistical significance of the correlation results was set at 95% CI.

Periodic signals were characterized by computing Lomb-Scargle periodograms [118,119] for each time series and both the NAO and AMO indices, providing a description of the overall frequencies. The periodograms were confined to the standard normalized power interval 0–1, according to Baluev [120], and limited to >4 years to focus on low-frequency signals.

In the land use analysis, only polygons within the Reno and Lamone river basins were considered. A further subdivision was made between the polygons within the alluvial plain and Apennine portions of the maps to highlight the various changes between the two sectors.

3. Results

3.1. Long-Term Linear Analysis

MK tests on all three river gauge stations indicated a statistically significant decrease in river liquid discharge (Figure 3a and Figure A1) over the last century. The Casalecchio di Reno time series, which is the most complete, showed a decrease of −0.18 ± 0.08 m³/s·year⁻¹ during 1921–2021 (Figure 3a). This station was more reliable than that at Ponte Bastia (Figure 1c) since, although the latter could depict the behavior of the lower Reno River basin, it was affected by massive anthropogenic contamination, mainly due to the Cavo Napoleonico (Figure 1c), one of the most significant hydraulic structures in the Po Plain, whose function since the late 1960s has been that of diverting the Reno River water into the Po River during floods and vice versa during drought periods. The river gauge at Pracchia also provides a nearly complete time series of the Reno River runoff (Figure 2). Despite its location in the high Apennine and the relatively limited drainage area upstream, a significant, slight decrease in river discharge occurred. In terms of seasonality (by considering the sole Casalecchio di Reno time series) the monthly mean data show clear and significant decreasing trends, with the lowest discharge values during summer (namely June, July, and August), marked decrease trends for both winter (December, January, and February) and autumn (September, October, and November), and a stronger decrease during spring (March, April, and May) (Figure A2a).

The stacked monthly rainfall time series (Figure 3b) showed a negative but not significant trend (−0.07 ± 0.08 mm·year⁻¹); the same was true for most of the individual rain gauges (Figure A3). Seasonality within the rainfall stacked curve shows the same behavior as that described by Pavan et al. [40], with the minimum values observed during summer and the maximum in autumn (Figure A2b); however, only the other two seasons (winter and spring), comparable in terms of both the amount of precipitation and trends, denote a significant decreasing trend.

The monthly mean temperature trends over the Reno and Lamone river basins during the last century almost all showed statistically significant increases (Figure A4). The stacked temperature curve (Figure 3c) shows a statistically significant increase of 0.016 ± 0.007 °C·year⁻¹, which corresponds to a warming of over 1.5 °C during the last century. The variability of temperature highlights an increasing and significant trend for all seasons (Figure A2c). During autumn, however, the trend is very slight, while stronger warming for the last century is observed in winter. Monthly mean temperatures ranged from 9 to over 14 °C at the highest station in the Apennines (Maresca) and urban areas in the plains (Bologna, Ferrara, and Imola), respectively, showing thermal variability due to both elevation and urbanization. Moreover, the latter has caused an urban heat island effect, which affects local and regional air temperatures near densely settled areas [121]; its main impact is a large increase in nighttime temperatures [122]. In addition, the local topography and land use have led to large variability in temperature, with the minimum (5.5 °C at Maresca) and maximum (18.6 °C at Bologna) averages differing by over 13 °C. However, grouping datasets from similar geographic contexts allows for a gross assessment of the variability induced by land use and local topography. Long-term minimum and maximum averages from the Apennine datasets were 0.2 and 2.4 °C lower than those from the rural plains, respectively. In contrast, maximum averages in the rural plains were not substantially different from those in urban areas (0.16 °C), whereas minimum averages were approximately 1.6 °C higher in urban areas.

The four-year low-pass filtered time series of wet days in Figure 3d shows a nonlinear, cyclic pattern, with minima during the 1920s, 1940s, and the early 1980s–2010. However, the >95th and >99th percentile values of the SDII do not give a clear pattern throughout the century; rather, they suggest an intensification of their occurrence, beginning in the 1970s.

Figure 3e shows the evolution of the SPEI during 1921–2019 in 48-month intervals. Paulo et al. [123] categorized the SPEI drought classes as non-drought (greater than −0.5), mild (−0.5 to −1), moderate (−1.5 to −1), severe (−1.5 to 2), and extreme (less than −2) conditions. In our data, moist and dry periods fluctuated throughout the century, in agreement with the wet days (Figure 3d); however, the wet periods have declined in recent decades. Since the 1980s, drought conditions have been increasingly persistent. Furthermore, extreme drought conditions only occurred during the 1940s and recently.

Trends in ground snow cover percentage (Figure A2d) over the area under study denote a marked and significant decrease in winter (−0.58 ± 0.18 %·year⁻¹) over the period 1950–2022, corresponding to a decrease in the winter snow surface of about 63% in the last 72 years; no significant trends emerged from both autumn and spring.

RSL data from the RTG time series (Figure 4a) clearly indicates a higher rate of sea-level rise during 1940–1980, compared to the prior and subsequent periods. This behavior witnesses the well-known subsidence phenomena due to local human activities (exploitation of gas and water from the underground) which overlapped with the natural land subsidence of the area [124,125]. However, subsidence has slowed significantly since the 1980s, mainly due to safeguarding policies introduced by the regional administration. Accounting for the TTG time series (Figure 4b) as a stable reference at the century scale, Zerbini et al. [126] evaluated the natural subsidence at the RTG site at around 1.88 mm·year⁻¹, by considering the period prior to the start of the mining activity (1873–1922); this value is in agreement with previous geological investigations [127,128,129] and with the first portion of the TTG-RTG time series of Figure 4c. The estimated trend of the RTG time series (1875–2016), detrended for VLMs (Figure 4d) is 1.3 ± 0.2 mm·year⁻¹, which is comparable with the TTG estimate of 1.33 ± 0.17 mm·year⁻¹ (1875–2021).

3.2. Nonlinear Signal

The EMD-decomposed time series of river discharge, rainfall, and temperature are shown in Figure 5a–c respectively, along with the sea level time series from the corrected RTG (Figure 5d). EMD decomposition produced seven IMFs for each of these four time series; however, to focus on long-term variability, only the lowest frequency IMFs (IMF4, IMF5, and IMF6) and RES were considered.

The RES component in river discharge was not perfectly linear but decreased over time; this decrease was 90% and 190% greater during 1993–2021 compared to that during 1950–1980 and 1921–1950, respectively. The RES in both rainfall and temperature show non-monotonic behavior, suggesting the presence of cycles longer than the time series themselves. In detail, rainfall decreased during 1950–2000 but otherwise increased. In contrast, temperature has increased since the 1920s; although it dropped slightly during the 1960s and 1970s, it has strongly increased since the 1980s. Conversely, the RES in sea level consistently increased over the last 140 years. RES trends (

Table 1. Datasets considered in this study and related residual trends.

Station	Variable	Residual Rate (*·Year−1)	Period
Casalecchio di Reno (Reno River)	River discharge (m3/s)	−0.202 ± 0.002	1921–2021
		−0.101 ± 0.004	1921–1950
		−0.190 ± 0.001	1950–1980
		−0.268 ± 0.002	1993–2021
Stacked rainfall	Rainfall(mm)	−0.101 ± 0.005	1921–2021
		0.157 ± 0.004	1921–1950
		−0.263 ± 0.005	1950–2000
		0.200 ± 0.006	2000–2021
Stacked temperature	Temperature (°C)	0.032 ± 0.002	1925–2021
		0.040 ± 0.001	1925–1950
		−-0.007 ± 0.001	1950–1980
		0.071 ± 0.001	1980–2021
Marina di Ravenna (VLM corrected)	Sea level (mm)	0.786 ± 0.008	1875–2016

Note: * symbol: unit of measurement of the specific variable.

) were statistically significant in all cases, even for rainfall, in contrast to the original time series (see Section 4.1).

The time series compared with the NAO and AMO climate indices (Figure 5e,f) show significant (95% CI) negative correlations of NAO phases with river discharge, rainfall, and sea level of −0.31, −0.16, and −0.34, respectively. In contrast, air temperature was positively correlated with both NAO (0.26) and AMO (0.21) fluctuations, and the sea level correlated with the AMO (0.14) index. However, only the air temperature-NAO correlation was statistically significant at the 95% CI, while the correlations with the AMO index become significant only if the 90% CI is considered. No correlations were found between the AMO with either river discharge or rainfall, nor was the WeMO index with any variable considered. Since no correlation was found on any annual time series representative of this area, the latter is not further discussed in this paper.

Several periodicities “contaminate” the time series, as detected by the LSP analysis (

Table 2. Most evident periodicities that affect the variables considered in this study, and their related Standard Normalized Power (Snp, according to Baluev [120]). Indications in parentheses for some periodicities highlight their dominance within a specific IMF.

Variable	Period	Snp
NAO	4–5	0.041
	~6	0.042
	7–8	0.087
	8–9	0.023
	~12	0.012
	~14	0.022
	16–17	0.011
	22–23	0.04
	37–40	0.04
AMO	4–5	0.015
	~6	0.018
	7–8	0.025
	8–9	0.05
	10–11	0.047
	16–17	0.018
	20–22	0.012
	~67	0.5
River discharge	4–5	0.056
	~6 (IMF 4)	0.073
	8–9	0.037
	10–11	0.072
	~12	0.068
	~14 (IMF 5)	0.058
	16–17	0.028
	22–23 (IMF 6)	0.056
	33–35	0.052
	53–56	0.097
Rainfall	4–5 (IMF 4)	0.094
	~6	0.038
	8–9	0.035
	10–11 (IMF 5)	0.052
	~12	0.027
	16–17	0.019
	23–26	0.021
	38–41 (IMF 6)	0.061
Temperature	4–5 (IMF 4)	0.04
	8–9	0.03
	~12	0.042
	~14 (IMF 5)	0.06
	16–17	0.02
	22–23 (IMF 6)	0.07
	30–33	0.022
	53–56	0.19
Sea level	4–5	0.039
	8–9 (IMF 4)	0.043
	10–11	0.02
	~12	0.038
	~14	0.076
	16–17 (IMF 5)	0.079
	22–23	0.06
	28–30	0.13
	40–42 (IMF 6)	0.29
	~67	0.16

). Most are common to all the datasets, albeit with different values. Periodicities of 4–5 and 8–9 years were detected in all the time series and are dominant within the IMF4 of rainfall and temperature and the IMF4 of sea level, respectively (Figure 5). Additionally, an approximately 6-year periodicity, which is dominant in the IMF4 of river discharge, is not present in the temperature and sea level. At lower frequencies, only the approximately 12-year period is common among all the time series, while periods of 10–11, approximately 14, and 16–17 years are dominant within IMF5 of rainfall, river discharge and temperature, and sea level, respectively. A periodicity of 22–23 years is common to all time series (dominant in the IMF6 of river discharge and temperature) when considering >20 years, whereas in rainfall and sea level, a powerful periodicity of approximately 40 years characterizes IMF6. Furthermore, two periodicities prevail at the longest wavelength: 53–56 years for river discharge and temperature, and approximately 67 years for sea level. All these periodicities occur in the NAO and AMO time series (Table 2), often with values common to both indices, and fluctuate in phase or anti-phase (as described above) with the IMFs of all the decomposed time series.

3.3. Land Use Changes

Recent land use changes within the Reno and Lamone River basins were categorized using the first level of the Corine Land Cover code as artificial areas, agricultural areas, forests and semi-natural areas, wetlands, and water bodies (see also Heymann et al. [130]). The land use changes in the Apennines for both basins (Figure 6) highlight that agricultural areas have been replaced since the 1970s, mainly by forests and semi-natural areas (approximately 70%) and some artificial surfaces (30%). Recent reforestation in the upper ER river basins has been driven by the migration of rural populations to urban areas and consequent cropland abandonment [32,34].

Trends in the plain sector have been nearly constant since the second postwar period, when the landscape began changing, following a rapid shift from a largely peasant economy to an industrial and touristic economy [131,132]. This has caused a gradual reduction in agricultural areas and intensive urban sprawl. This “uncontrolled spread of towns and villages into undeveloped areas” [133] has increased soil sealing due to impermeable artificial materials, which are considered the largest issue for pervasive soil degradation [134]. Moreover, the area of sealed soils (excluding areas such as gardens, urban wastelands, and parks) has increased by approximately 123% since 1976. Thus, sealed surfaces accounted for approximately 83% of the total artificial areas in all six analyzed maps. Figure 6 shows that forests and wetlands represent very little of the land use in the plains, which is a consequence of human activities beginning with progressive deforestation by the Romans during the first century B.C. and land reclamation of swamps and salt marshes, begun in the 1800s and completed in 1964. Finally, the slight increase in water bodies in the plains was primarily due to progressive increases in artificial basins (over 120%) and, to a lesser extent, in waterways and irrigation canals (approximately 45% increase).

4. Discussion

4.1. Climatic and Hydrological Variability Induced by the Anthropogenic Dimming/Brightening Phenomenon

Surface water runoff has decreased in several areas of the ER region in recent decades [135]. Local water scarcity has been driven by the onset of persistent drought conditions, which began in the 1980s and have worsened, as also suggested by the SPEI (Figure 3e); conditions are projected to further deteriorate over the next decades [136,137].

In our analysis, air temperature (Figure 3c) showed an increasing trend over the last century, as previously observed throughout the ER region by Antolini et al. [38]. However, the RES curve of temperature (Figure 5c) showed an increase from the 1920s to 1950s, a generally stable (or slight cooling) period until the late 1970s, and rapid warming thereafter. This general behavior in temperature, as previously observed throughout Italy [138,139,140], Europe [141,142], and worldwide [143], has been linked to a dimming/brightening phenomenon [144,145] from surface solar radiation (SSR) changes over the last century [146,147] driven by anthropogenic air pollution, which altered the transparency of the atmosphere [148,149,150]. The marked change around 1980 led to a “global warming-type drought” [151,152]; its effects on water evaporation and local temperature have been further accentuated by the presence of artificial surfaces.

Despite the lack of statistical significance (Figure A3), the decrease in rainfall over the last century over large portions of the ER was mainly attributed to a decrease in spring and winter rainfall [38,39,40,153,154], as confirmed by our analysis (Figure A2b). However, when the oscillations are removed, a significant decrease and non-monotonic behavior are visible (Figure 5b). The latter is likely due to the anthropogenic dimming/brightening process described above, as any change in SSR leads to large modifications in the water cycle [155,156]. Wild [157] have found evidence of concurrent decreases in rainfall and SSR from the 1950s to 1980s across the Northern Hemisphere, followed by rainfall intensification from the 1980s onward because of subsequent brightening and increased evapotranspiration [158]. This is consistent with our results, although rainfall resumed long after the 1980s (Figure 5b). This may be due to climatological moisture divergence over the Mediterranean, which transports water vapor outside the region; therefore, moisture does not fully translate to regional increases in rainfall [159,160,161]. Moreover, this phenomenon is likely connected to a possible increase in rainfall extremes [161,162], as also suggested by the increase in the >95% SDII values from the 1970s (Figure 3d). This could be strictly linked to the concordant reduction in rainy days, and in this case, not necessarily associated to a substantial decrease in rainfall totals. It can be argued that, following this decline in the 1970s, the amount of total monthly rainfall has not changed considerably (Figure 3b and Figure 5b), and therefore falls in a shorter time.

4.2. Anthropogenic Influences at the River-Basin Scale

The river discharge time series in this study shows a significant decreasing trend since the 1920s (Figure 3a and Figure A1), in agreement with Preciso et al. [37]. Unfortunately, data on bedload materials are lacking, and the relationship between them and river flow discharge cannot be verified [32]. As stated in Section 1, the Reno catchment no longer represents the primary sediment source for the coast, primarily because of reforestation in the catchment since the 1950s [37]. This was confirmed by our analysis (Figure 6); increased forest cover represents the prevailing land use change in the Apennines portion of both basins (see [163]), which likely influenced the reduction in discharge and sediment production. The decrease in discharge showed a negative acceleration in the RES curve (Figure 5a) during 1950–1980, compared to the previous 30 years, which is due to the influence of both river regulation and reforestation, as pastures, shrubs, and crops, consume much less water than forests [164,165,166]. However, the huge decrease during the following years (1980–2021; Table 1 and Section 4.2) cannot be attributed to river regulation. First, the increase in reforested areas has remained constant, although it has slowed since the mid-1990s. Second, the notable increase in temperature since the 1980s (Figure 5c) has significantly increased the atmospheric evaporative demand and led to persistent local drought conditions (Figure 3e), which, coupled with evapotranspiration, may have driven further declines in discharge. The non-monotonic variation in rainfall (Figure 5b) also contributed to the decreased river discharge since the 1950s; however, its significant but slight changes played a secondary role to the above processes. Furthermore, another consequence of the increasing temperatures may be the decreasing extent of snow cover and number of days with snow present on the Italian peninsula [41] and throughout the Northern Hemisphere [167] over the last century, especially since the late 1970s. Declining snow cover extent throughout the upper catchment areas during winter (Figure A2d) could be a further driver of the reduced river discharge by weakening or even depriving the fluvial systems of water supplied by snowmelt [168,169]. The decline in snow cover extent throughout the upper portion of the studied catchments area can be ascribed as the main driver for the deterioration in spring trends of river discharge (Figure A2a), by weakening, or even depriving, the fluvial systems of the water supply due to the melting of snow.

Once eroded, soil particles can undergo transportation-deposition cycles throughout their travel downstream and either remain within the basin or outflow from it. Thus, processes throughout the plains, usually of anthropogenic origin, are crucial for sediment supply. Our analysis noted that agricultural areas are substantially decreasing in the plains of both basins (Figure 6) and are being replaced by artificial surfaces and, to a lesser extent, water bodies. The increase in the latter, despite the small area, may have significantly influenced discharge by withholding water from the system. Artificial basins (+120% during 1976–2017) and irrigation channels (+45%) have been continuously implemented throughout the plains to satisfy increasing water demand. Intensive groundwater withdrawals in these areas, mainly for industrial and domestic needs [170], also contributed significantly to the decreased discharge, as over-exploitation of aquifers weakens (and often prevents) natural river recharge from the aquifers. The increased temperatures and UHI effect also directly influence the evaporative processes and water demand for both agriculture and domestic needs in the plains. These processes have been magnified during recent decades owing to urban sprawl, the increase in artificial surfaces (Figure 6), and soil sealing (+123% during 1976–2017). The latter reduces and sometimes prevents water infiltration into soils and, therefore, aquifer recharge, increasing both water scarcity and dependence on irrigation and artificial reservoirs [171].

4.3. Natural Contributions from Coupled Ocean/Atmosphere Processes

The NAO influences the weather and hydrology by controlling the wind intensity and direction and the interactions among air masses over much of Europe, and the AMO influences several hydroclimatic variables, including Mediterranean air temperature (e.g., [172,173,174,175]). Additionally, the sea level in the Mediterranean is influenced by the NAO, which drives atmospheric sea-level pressure changes [63,176] and the net water flux at Gibraltar by directly altering wind and oceanic circulation near the strait [177,178] and indirectly by impacting regional river runoff, evaporation, and rainfall [179]. Moreover, the AMO influences the SST and sea evaporation [174,179].

Previous studies have found significant correlations between the NAO and Mediterranean surface air temperatures [62,174,180,181,182,183]. Our findings of a positive correlation between temperature and AMO cyclicity (Table 2) agree; however, they cannot be confirmed because of the statistical non-significance of the correlation at the 95% CI. Mariotti and Dell’Aquila [174] have also found no correlation between the AMO variability and rainfall over the Mediterranean; our study confirms this apparent lack of influence.

Several authors have found significant negative correlations between rainfall and NAO over Southern Europe and the Mediterranean, with the positive phases of the NAO linked to drier conditions and vice versa [62,63,139,159,181,184,185]. Generally, positive and negative NAO phases decrease and increase rainfall, respectively, over Southern Europe through anomalous behavior of the Atlantic westerlies [184,186], thereby severely impacting regional hydrological cycles [159]. This is reflected by river discharge variability, as indicated by our analysis, which shows a significant negative correlation between NAO phases and river runoff (Figure 5a,e). Negative correlations between river discharge and NAO have been found throughout Europe [187,188,189,190]. However, López-Moreno et al. [191] note that, unlike rainfall, river discharge is characterized by a nonlinear response of short, but intense, discharge anomalies during negative NAO phases and persistent, but weaker, anomalies during positive phases. The same authors highlight that during positive NAO years, surface runoff and aquifer recharge are severely reduced, whereas during negative NAO years, river basins receive more rainfall, leading to soil saturation and rapid rainfall transfer to runoff. Positive anomalies, such as the strong peak observed in river discharge in 2010 (Figure 5a), may be a consequence of the negative NAO phase that occurred around 2010, which might have left a record in the river basins and coastal systems. Granulometric analysis performed by ARPAE (Agenzia Regionale per la Prevenzione, l’Ambiente e l’Energia—Emilia Romagna) on samples collected in 2012 in the submerged beach of the studied coastal tract [25] showed larger grain sizes than those from later [36] or prior [33,34,35] regional sampling. Excluding the possible effects of local interventions (i.e., beach replenishment) [36], these variations in grain size may be related to the large, positive anomaly in river discharge that occurred in 2010. This natural anomaly could have momentarily strengthened the solids’ contribution to the coast by enhancing the ability of rivers to carry coarser sediment across the submerged beach through the river mouth bypassing and subsequent redistribution by longshore currents. Thus, the generalized increase in grain size observed in the studied coastal stretch and elsewhere in the nearshore ER coast in 2012 could have been associated with this short-term event, which temporarily disrupted the tendency towards gradually finer and weaker sediment supply observed in recent decades.

At the ER coast, persistent sea-level rise (Figure 5d), periodically altered by various cyclicities, inevitably impacted shoreline evolution. The large positive sea-level anomaly shown in the IMF4 in Figure 5d, has also been observed in previous studies throughout the Mediterranean [178,192,193] and has been attributed to a strongly negative NAO [194]. The dominant periodicity observed in sea-level fluctuations was approximately 8–9 years, which is comparable to the approximately 10-year periodicity found by Bonaduce et al. [112]. Additionally, Galassi and Spada [111] associated the 2010 sea-level anomalies with the coincidence of both negative NAO and positive AMO phases, which has a periodicity of approximately 20 years.

Since the 1990s, sea level has risen rapidly throughout the Mediterranean, which has been linked to increased SST [195], mostly controlled by increasing air temperatures. This also happened during the early 20th century, the other major period of constant rise in global sea and air temperatures [196]. The sea level oscillations at the RTG during these periods are described by IMF6 (Figure 5d), which has a dominant periodicity of approximately 40 years (Table 2). Indeed, major sea-level changes at the RTG over the last 140 years (excluding the effects of local VLMs) were described by multiple phases of acceleration, during the early 20th century, the 1970s–1980s, and 1990s to mid-2000s. Moreover, this behavior is highly phase-correlated with the rainfall IMF6 (Figure 5b), which is also characterized by a similarly dominant periodicity. The NAO may be the driver of this process, because it has a similar periodicity (approximately 37–40 years) and a good negative correlation with both sea level and rainfall. For the same study area, Meli et al. [114] noticed a slight negative acceleration in sea level change since the mid-2000s and attributed it to the possible presence of periodic signals at periods longer than the 26-year period they considered. Therefore, this recent change in the local sea level trend could be associated with the above periodicities, as the final portion of the sea-level IMF6 denotes a temporally concordant, progressive flattening of the rising trend that began in the early 1990s.

5. Conclusions

The environmental and climatic data considered in this study provide a holistic view of the hydrological and climatic processes that have affected the Reno and Lamone river basins of the ER coastal areas in Northern Italy over the last century. The analysis provided responses to the questions posed in Section 1.

• The anthropogenic footprint, attributed to effects such as land use changes and the large-scale anthropogenic dimming/brightening phenomenon, acting at both the river basin and regional-to-global scales, profoundly impacts the catchment dynamics by driving long-term, nonlinear trends, upon which natural oscillations are superimposed. Interactions with major climate modes, in fact, affect the signals over various periodicities. Both positive and negative correlations among some of the studied parameters and the main climatic indexes (NAO, AMO, WeMo) are evidenced;
• The marked negative acceleration in river discharge provides an indication of the lack of recovery in terms of sediment supply to the coast, despite safeguard policies introduced by the ER regional administration in the early 1980s. This decline resulted from river regulation and land use changes during 1950–1980 and related implications. Moreover, since the 1980s, local air temperatures have increased significantly, leading to persistent drought conditions contributing to the drastic reduction in river discharge. However, periodic natural signals can significantly restore the river discharge, as observed during the strongly negative NAO event in 2010;
• Persistent sea-level rise has affected the coastal site under study over the last 140 years. Apart from the extremely high rates during the 1950s–1980s due to anthropogenic-induced land subsidence, sea level is periodically amplified or reduced by natural fluctuations, such as that observed in 2010, and lower frequency fluctuations in the local sea level, as during the early 20th century, 1970s–1980s, and 1990s to mid-2000s. This nonlinear behavior, locally enhanced by subsidence that overwhelms the river sediment inputs, notably impacted the ER coast and may continue to do so.

Overall, the interconnection found between various climatic, hydrological, and anthropogenic processes with individual nonlinear behaviors illustrates the need for a holistic approach when considering complex environments such as LECZs and heavily anthropized areas. Expanded observation and monitoring at the scale of entire river basins is necessary to better understand the natural and anthropogenic influences on coastal behavior. This approach should be thoroughly considered in coastal management and in adaptation strategies, especially when facing a changing climate and increased vulnerability in the most threatened environments, such as LECZs.