Insights into exceptional freshening events in the northern Adriatic Sea throughout 2023–2024

The northern Adriatic Sea experienced multiple extreme freshening events in 2023 and 2024, underscoring the intensifying influence of heavy precipitation, river discharge, and coastal flooding in the region. This study combines high-frequency in-situ observations from two autonomous meteo-marine monitoring platforms offshore of Fano and Senigallia with river discharge records, satellite imagery, and Copernicus Marine Service products to examine the evolution, drivers, and ecological implications of these events.

The first freshening event in early summer 2023 was linked to consecutive heavy rainfall episodes, including Storm Minerva, which triggered catastrophic flooding across the Emilia-Romagna region in northern Italy. Near-surface salinity declined from ∼ 36 g kg⁻¹ to below 24 g kg⁻¹ within two weeks, corresponding to an estimated freshwater input of ∼ 42.5 m³ h⁻¹ in the vicinity of the buoy. Events in late 2023 and 2024 were even more abrupt and intense, culminating in late October 2024 with salinity values dropping to ∼ 15 g kg⁻¹ following Storm Boris. A key distinction between the two years was the marked recovery of Po River discharge in 2024 after a prolonged drought, which significantly amplified freshwater inputs. These events triggered cascading environmental responses, including stratification, turbidity peaks, and phytoplankton blooms including summertime mucilage-forming aggregations, posing risks to aquaculture, water quality, and coastal tourism.

By combining autonomous in-situ measurements, satellite observations, and Copernicus reanalysis products, this study demonstrates a robust and cost-effective approach to tracking extreme hydrological events in coastal zones. The integrated dataset not only captures rapid thermohaline and biogeochemical changes near-real time but also places them in a broader environmental context. These findings underscore the value of sustained and interoperable observing systems in supporting early warning capabilities, scientific research, and climate adaptation in vulnerable marine regions such as the northern Adriatic Sea.

1 Introduction

Coastal regions like the shallow Northern Adriatic are complex and densely-populated socio-ecological systems that provide high-value services (Costanza et al., 2014) which are susceptible to the impacts of climate change and other anthropogenic and natural pressures.

Table 1 Products from the Copernicus Service and other complementary datasets used in this study, including the Product User Manual (PUM) and Quality Information Document (QUID). For complementary datasets, the links to the product description, data access and/or references are provided.

Among these, flooding events, triggered by factors such as high tides, storm surges and precipitation whiplashes can lead to large adverse impacts on human and natural systems (Zscheischler et al., 2018; Wright and Nichols, 2019; Tan et al., 2023). According to the latter, precipitation whiplashes, including abrupt shifts between dry and wet extremes, favor the generation of significant surface runoff and flooding, concurrent landslides and erosion that can also amplify negative trends in water quality. Overall, river discharge plays a crucial role in coastal ecosystems by providing freshwater and nutrients that impact hydrology, biogeochemistry, and productivity. Along the Northwestern Adriatic coast, several small rivers and the major Po river discharge nutrient-rich freshwater which is mainly transferred southward along the Italian coast by the Western Adriatic Current (Goudeau et al., 2013), while the Po plume is known to remain close to the coast under stratified conditions (Cozzi and Giani, 2011). This, in turn, affects nutrient cycling, species growth and survival rates (Ludwig et al., 2009; Cozzi and Giani, 2011) as well as water mass formation processes (Falcieri et al., 2014).

The coastal impacts of flooding events extend beyond environmental concerns, affecting local communities and economies that depend on infrastructure and marine resources (e.g., IPCC, 2022, 2023). The need for coastal observations and monitoring efforts is therefore increasingly recognized as essential for guiding policy and public safety needs. The Copernicus Marine Service plays a key role in this context by offering a wide range of physical and biogeochemical near-real time data. However, the location, rapidity and complexity of floods and the related coastal events present specific challenges. Accurate detection and monitoring require an integrated methodology that combines local high-resolution in-situ measurements, satellite data and modeling approaches. This study leverages the high-resolution data from a coastal meteo-marine buoy and several Copernicus sources to provide a comprehensive analysis of a series of exceptional freshening events in the northern Adriatic Sea throughout the last 2 years.

2 Methodology

In order to assess the spatio-temporal characteristics, environmental drivers, and ecological impact of extreme freshening events in the northwestern Adriatic Sea during 2023 and 2024, we integrated a set of in-situ, satellite-derived, reanalysis, and model-based datasets.

2.1 Local meteo-marine conditions

High-frequency in-situ measurements of thermohaline, atmospheric and biogeochemical properties were obtained from two autonomous meteo-marine monitoring systems, positioned ∼ 3.7 km offshore of Fano (43°53^′30.0^′′ N, 13°00^′30.0^′′ E) and Senigallia (43°45^′20.9^′′ N, 13°12^′32.4^′′ E) on the western coast of the northern Adriatic Sea (product ref. 1 and 2; Table 1).

Raw data from these platforms (mentioned as Buoy Fortunae offshore of Fano and Meda Station offshore of Senigallia henceforth) were first converted with the calibration coefficients in the configuration file provided by Sea-Bird Electronics and processed according to the typical SBE data processing sequences. The necessary calculation and analysis of thermodynamic properties was based on the Gibbs Seawater (GSW) Oceanographic Toolbox (IOC et al., 2010), which contains the TEOS-10 subroutines (http://www.TEOS-10.org, last access: 12 February 2025) for the evaluation of seawater properties (McDougall and Barker, 2011). Detailed sensor specifications and processing steps are available in Appendix A1.1–A1.3.

2.2 Precipitation and river discharge

Hourly precipitation data were retrieved from the ERA5 reanalysis dataset, produced by the Copernicus Climate Change Service, which provides gridded atmospheric fields at 0.25° spatial and 1 h temporal resolution (product ref. 3; Table 1). In parallel, 6-hourly river discharge estimates were extracted from the European Flood Awareness System (EFAS) model, based on the LISFLOOD hydrological model forced with gridded meteorological observations (product ref. 4; Table 1). Considering the dominant southeastward flow of the Western Adriatic Current, freshwater input from major (Po, Adige, Reno) and minor regional rivers (e.g., Lamone, Savio, Metauro, and 14 others) was assessed at their upstream mouths along the northwestern Adriatic coast (see Fig. 3 for river mouth locations and characteristic circulation patterns).

2.3 Freshwater rate estimation

To estimate the required freshwater inflow necessary to account for the observed decrease in surface seawater salinity over a given time period, a simplified salt mass conservation approach was applied. This method was based on in-situ salinity data from the Fortunae buoy (product ref. 1; Table 1) and the estimated area of the affected region, derived from satellite SST and ocean-colour observations (product refs. 5 and 6; Table 1) capturing the cold and turbid coastal plume. Assuming that salinity and density are conserved during mixing and proportional to the volumes involved, the following salt mass balance equation was applied:

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_{\mathrm{1}}{S}_{\mathrm{1}}{V}_{\mathrm{1}}={\mathit{\rho }}_{\mathrm{2}}{S}_{\mathrm{2}}{V}_{\mathrm{2}}\end{array}$$

where ρ₁ and ρ₂ are the initial and final density in kg m⁻³; S₁ and S₂ are the initial and final salinity in g kg⁻¹; V₁ is the initial considered volume of seawater in m³; V₂ is the final volume after mixing with the added freshwater in m³.

And the freshwater flow (Q) in m³ s⁻¹, required to achieve the observed salinity decrease over the time period Δt is calculated as:

$$\begin{array}{}\text{(2)}& Q={\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }t}}\end{array}$$

A practical example of this calculation, including measured values from the Fortunae buoy and assumptions on the affected volume based on the buoy's position and local bathymetry, is provided in Appendix A2, along with a schematic diagram illustrating the considered water volume (Fig. A1).

2.4 Complementary satellite and model-based datasets

To complement the in-situ observations and support the freshwater rate estimation (Sect. 2.3), we incorporated satellite-derived and model-based datasets (products 5–8; Table 1) providing full spatial coverage of the northwestern Adriatic Sea. These included daily fields of sea surface temperature (SST; 0.01° × 0.01°), chlorophyll a concentration (Chl-a, 1 km × 1 km), light attenuation coefficient at 490 nm (KD490; 1 km × 1 km), surface salinity (SSS; 0.042° × 0.042°), and horizontal current components (u and v).

These datasets served three main purposes:

1. to delineate the surface extent of freshwater-affected waters and define the horizontal area used in freshwater input calculations (see Appendix A2);

2. to fill temporal gaps in the buoy and coastal station records, ensuring continuous environmental context during key periods of data loss; and

3. to extend the analysis beyond the nearshore region, offering insight into the broader evolution of surface thermohaline and biogeochemical properties across the Adriatic Sea.

The satellite products provided gap-free, multi-sensor interpolated surface fields, while model outputs offered physically consistent estimates of key oceanographic variables across space and time.

2.5 Statistical analysis

Basic statistics were computed following Emery and Thomson (1998). Root Mean Square Difference (RMSD), bias, and cross-correlation were applied to evaluate temporal coherence between thermohaline, biogeochemical and environmental conditions. Detailed equations and assumptions are reported in Appendix A3.

3 Results

Over the past two years, the northern Adriatic Sea experienced a series of hydro-meteorological extremes, resulting in rapid and recurrent freshening episodes with cascading environmental and societal implications.

3.1 Surface freshening

Between May 2023 and November 2024, five distinct surface freshening events were identified in the northern Adriatic Sea from high-frequency salinity observations collected by the Fortunae buoy and the Senigallia Meda Station. These events varied in duration (12 to 53 d), intensity (6.7 to 20.1 g kg⁻¹ salinity drop), and temporal structure, shaped both by prolonged and abrupt freshwater input dynamics (Tables A4 and A5).

Event I (2 May–24 June 2023) was the longest, spanning 53 d. It was characterized by a multi-phase freshening, including an extended period of anomalously low salinity and a gradual recovery. The most pronounced salinity decrease occurred over 13 d (11–24 May), with salinity dropping by 12.77 g kg⁻¹, from 36.19 to 23.42 g kg⁻¹. This prolonged freshening suggests sustained riverine discharge into the northwestern Adriatic. Event II (29 September–16 October 2023, 17 d) showed a more abrupt signature. Over just 5 d, salinity fell by 6.72 g kg⁻¹ (from 37.17 to 30.45 g kg⁻¹), indicating a rapid freshwater input followed by a relatively short recovery. Event III (12 November–2 December 2023, 20 d) included the highest initial salinity observed (37.63 g kg⁻¹), followed by a significant 14 d drop to 26.53 g kg⁻¹, which corresponds to a total decrease of 11.10 g kg⁻¹. Event IV (9–21 April 2024, 12 d) was short but intense, featuring a 6 d salinity decrease of 13.89 g kg⁻¹, from 36.16 to 22.27 g kg⁻¹, representing one of the steepest gradients recorded. Event V (20 October–22 November 2024, 33 d) was the most extreme in terms of salinity reduction. A 10 d drop of 20.12 g kg⁻¹ brought surface salinity from 35.09 to 14.97 g kg⁻¹, making this the most intense freshening event of the record.

These events illustrate the basin's dual sensitivity to sustained meteorological forcing and intense short-term hydrological extremes. Cross-event comparisons emphasize differing recovery behaviors and highlight the coupling between catchment saturation and downstream marine responses (Figs. A2–A6).

3.2 Environmental connections

Hydrological analysis revealed that major freshening events coincided with intense precipitation and river discharge anomalies in the Emilia-Romagna region, with events driven by both catchment saturation and localized convective storms.

A sequence of five rainfall episodes culminated in catastrophic flooding from 1–24 May 2023 (Fig. 2b). The first rainstorm (1–3 May) brought over 200 mm of precipitation to central-eastern Emilia-Romagna. Initial flooding affected the Reno and Lamone basins. Two subsequent moderate episodes (9–10 and 12–14 May) further saturated soils. The most severe storm, associated with Storm Minerva (16–18 May), delivered over 300 mm in areas such as Forlì, with cumulative precipitation more than 5 times the monthly climatological average.

During this period, the Po River remained well below the climatological average discharge, with values around 2000 m³ s⁻¹, as the most intense rainfall did not affect its main basin, whereas local rivers like the Lamone, Savio, and Reno exceeded historical flood levels (Fig. 3c). Cross-correlation analysis between accumulated precipitation and river discharge (Fig. A5) confirmed statistically significant relationships (p < 0.05), with time lags ranging from −6 to −120 h. Rapid response was observed in local rivers with short catchments, such as Rubicone (r = 0.74, lag = −6 h), while larger basins like the Adige exhibited longer delays (r = 0.71, lag = −36 h).

Figure 1 Hourly (colored lines) and daily-averaged (black dotted lines) time series of (a) conservative temperature, (b) absolute salinity, and (c) potential density anomaly recorded at 2 and 5 m depth offshore of Fano and Senigallia, respectively (product ref. 1 and 2; Table 1). Basic statistical descriptors (minimum, maximum, median, standard deviation) for each thermohaline property are reported in the top left corner of each panel. Light red shaded areas (I–V) indicate the timing of the five most significant surface freshening events identified between 2023 and 2024.

No widespread heavy precipitation was recorded in correspondence with Event II in late September 2023, but localized convective storms likely impacted the Marche region, just south of Emilia-Romagna, between 25–29 September. An increase in the discharge of the Adige and Po was reported and the freshening signal suggests a short-lived but intense localized input (Fig. 1). In accordance, two prolonged rainfall episodes occurred later during 2–5 November and 30 November–2 December (Event III). The latter was the most severe, with up to 383 mm of rainfall recorded inland, in the Parma Apennines. A third, minor event on 4–5 December followed the main storm.

During spring 2024, an unseasonable storm system between 20–24 April brought 70–80 mm of weak but steady rainfall and rivers with short response times reacted quickly to the precipitation input, leading to the abrupt salinity drop during Event IV (Fig. 1).

Autumn 2024 saw the return of exceptional rains. Between 17–20 September, Storm Boris delivered ∼ 300 mm of rainfall in eastern Emilia-Romagna, causing major river discharges in local rivers like the Savio, Fiumi Uniti and Lamone River (Fig. 3c). Less than a month later (18–19 October), another storm triggered river discharge records for all major rivers, especially for the Po River (record value = 7880 m³ s⁻¹), reflecting a marked hydrological recovery after prolonged drought.

This analysis affirms that both large-scale catchment saturation and localized storm activity contributed significantly to the observed salinity anomalies. Rivers with differing catchment sizes showed distinct response times, with immediate peaks in small rivers and delayed responses in larger ones. Figures 2 and 3, along with Appendix Figs. A2–A6, provide a spatial and temporal context for precipitation and discharge patterns associated with each event.

Figure 2 Spatial distribution of accumulated precipitation (product ref. 2; Table 1) over Emilia-Romagna and surrounding areas during key extreme rainfall events associated with surface freshening Event I and Event V (see Fig. 1). (a) Total precipitation during the combined May 2023 period (Event I). Accumulated precipitation during individual May 2023 storm episodes: (b) 1–3 May, (c) 9–10 May, (d) 12–14 May, and (e) 16–18 May. (f) Total precipitation during the two September–October 2024 events (Event V). The Emilia-Romagna region is outlined in red in each panel. (g) Time series of daily accumulated rainfall volume (m³ d⁻¹) over the domain, highlighting the timing and intensity of individual rainfall peaks between January 2023 and December 2024. Basic statistical descriptors (minimum, maximum, median, standard deviation) are reported in the upper left corner.

Figure 3 Overview of (a) the study area with the location of the 20 assessed river mouths (blue dots) and the positions of the autonomous monitoring platforms offshore of Fano and Senigallia (dark and light red star, respectively). The background corresponds to a Sentinel-2 Level 2A true-color image (https://dataspace.copernicus.eu/browser, last access: 17 August 2025) acquired on 31 October 2024, overlaid with surface current vectors derived from the numerical circulation model (product ref. 8; Table 1) for the same date. (b) Daily time series of the total accumulated river discharge (m³ s⁻¹) from all 20 rivers for the period January 2023–December 2024 (product ref.3; Table 1). Key descriptive statistics (minimum, median, maximum, standard deviation) are indicated. (c) Individual river discharge time series, grouped by magnitude, with numeric labels corresponding to the river mouth positions indicated in panel (a).

3.3 Ecological implications

Each freshening event identified in this study (Fig. 1) triggered a marked biogeochemical response in the northern Adriatic Sea, as captured by the high-frequency in-situ measurements of chlorophyll a, turbidity, and dissolved oxygen collected by the Fortunae Buoy and the Senigallia Meda Station (Fig. 4). These parameters serve as proxies for primary productivity, sediment transport, and oxygen dynamics, offering direct evidence of ecosystem perturbations following extreme hydro-meteorological events.

Figure 4 Hourly (colored lines) and daily-averaged (black dotted lines) time series of (a) chlorophyll a concentration, (b) turbidity, and (c) dissolved oxygen concentration measured at 2 and 5 m depths by the Fortunae buoy (offshore of Fano) and the Senigallia Meda station, respectively (product ref. 1 and 2; Table 1). Basic statistical descriptors (minimum, maximum, median, standard deviation) are reported in the left corner of each panel. Red shaded areas (I–V) denote the timing of the five most pronounced surface freshening events as identified in Fig. 1. Biogeochemical responses, including sharp increases in turbidity and chlorophyll a, coincide with or follow these events, reflecting ecosystem perturbations.

Following Event I (May 2023), which was preceded by a sequence of five intense rainfall episodes in the Emilia-Romagna region and culminated in catastrophic flooding (Fig. 2b–e), an extreme turbidity increase (∼ 70 NTU average) was recorded for more than a month, along with chlorophyll a concentrations exceeding 10 µg L⁻¹ (Fig. 4). These responses reflect the combined influence of excess sediment transport and nutrient enrichment, especially from local rivers like the Lamone, Reno, and Savio, which exceeded historical discharge levels (Fig. 3c). The sustained rainfall and resulting prolonged runoff seem to have led to a persistent stratified surface layer that confined biological activity and promoted bloom formation.

After the record May flood period, a particularly notable anomaly was observed in June–July 2023, with sustained extreme turbidity (> 200 NTU) and following chlorophyll a values (> 30 µg L⁻¹). Subsequent visual inspections of the buoy infrastructure revealed the presence of mucilage-forming macroaggregates, consistent with recurring gelatinous algal blooms historically observed in the Adriatic under warm, stratified, and nutrient-enriched conditions (Arrighi and Domeneghetti, 2024). These mucilage events have well-documented impacts on marine ecosystems, coastal water quality, and economic sectors such as aquaculture and tourism.

The response to Event II (late September–October 2023) presented a contrasting case. While no widespread heavy precipitation and flooding occurred in Emilia-Romagna, localized convective storms over the Marche region and rising discharges in the Adige and Po rivers triggered a short-lived but intense freshening. While turbidity remained low, the biogeochemical response included a marked increase in chlorophyll a concentration and oxygen supersaturation reaching approximately 150 %, indicating a nutrient rich but low sediment freshwater input. Moreover, dissolved oxygen rose near-synchronously with chlorophyll a concentrations (r = 0.77, lag = 0 h), consistent with well-lit bloom conditions. This event exemplifies how smaller, convective systems can fuel biological productivity without large sediment loads, while favorable light conditions likely allowed rapid phytoplankton growth.

The prolonged rainfall in early and late November 2023 preceded Event III, which produced another strong bloom and variable oxygen patterns. This response suggests that delayed runoff from larger, saturated catchments delivered substantial nutrient loads, with peak phytoplankton activity followed by oxygen variability related to organic matter remineralization and light attenuation during turbid phases.

The corresponding biogeochemical response to event IV followed an unseasonable storm system in April 2024. Rapid response from short catchment rivers like the Rubicone and Uso led to an abrupt salinity drop (Fig. 1) and sharp turbidity and chlorophyll peaks, consistent with the fast land-sea transfer typical of springtime storms over saturated terrain. These dynamics emphasize the system's sensitivity even to moderate precipitation when antecedent soil conditions favor immediate runoff.

Finally, Event V in autumn 2024 was preceded by two major storm systems: Storm Boris (17–20 September) and an intense October storm (18–19 October), the latter driving the Po River to record discharge levels (Fig. 3c). This sequence produced the most intense surface freshening and the highest chlorophyll a concentration observed during the study period, exceeding 60 µg L⁻¹. Interestingly, these extremely high chlorophyll a levels occurred in conjunction with notably low turbidity, suggesting that sediment-laden runoff from the Po River may have settled or dispersed before reaching the offshore monitoring sites at Fano and Senigallia. This decoupling implies that nutrient transport remained effective despite reduced suspended particulate matter, allowing for favorable light conditions and enhanced primary productivity which is consistent with a multi-day negative turbidity–chlorophyll relationship (r = −0.32 at −144 h lag) and near-synchronous, positive chlorophyll–oxygen coupling (r = 0.86 at 0 h lag), together indicating strong photosynthetic control under clear conditions.

Across all events, dissolved oxygen dynamics reflected a combination of photosynthetic production during bloom conditions and oxygen depletion during periods of high turbidity and organic matter degradation. Event-scale lead–lag patterns showed a biphasic turbidity–chlorophyll response with a short-lag negative branch when high turbidity limits light, and a positive branch at multi-day delays when nutrient delivery and stratification promote blooms. The chlorophyll–oxygen coupling is typically positive and near-synchronous, with longer-lag negatives appearing during potential bloom decay and post-bloom respiration.

These findings confirm that hydro-meteorological forcing, whether through widespread rainfall, localized storms, or antecedent soil saturation, directly drive ecosystem-level responses in the northern Adriatic Sea. The magnitude and timing of these responses depend not only on precipitation intensity and discharge volume, but also on the sediment load, nutrient delivery, and seasonal conditions governing stratification and light availability.

4 Discussion and conclusions

The northern Adriatic experienced a series of extreme freshening events between 2023 and 2024, underscoring its vulnerability to compound hydro-meteorological extremes. Events such as Storm Minerva and Storm Boris were marked by intense precipitation and widespread flooding, with significant hydrological and ecological repercussions. These cases highlight the increasing role of precipitation whiplash and flood-drought cycles in modulating coastal hydrography.

According to the IPCC's Sixth Assessment Report, heavy rainfall events and the corresponding overflow are projected to increase in frequency and intensity as global temperatures rise beyond 1.5 °C (IPCC, 2023). River floods in central and western Europe are already increasing, and small catchments, like those affected along the northwestern Adriatic coast, are especially vulnerable. Future storms may more frequently exceed their coping capacities, reinforcing the urgency of adaptive monitoring systems.

In the Northern Adriatic, where the Po River and smaller local rivers drive salinity and biogeochemical patterns, such shifts have significant implications. As the shallowest and northernmost basin of the Mediterranean, it is particularly vulnerable to such perturbations due to its limited flushing, strong land-sea coupling, and intensive human pressure (Coll et al., 2010). Its physical and ecological sensitivity to episodic freshwater pulses, such as those examined here, underscores the growing importance of continuous monitoring in the context of climate change.

Surface freshening alters water column structure, promotes stratification, and modulates ecological responses. In summer 2023 and 2024, sea surface temperatures exceeded 30 °C, accompanied by extensive mucilage-forming algal aggregations. These events compromise water quality, threaten aquaculture viability, and negatively impact tourism, posing socio-ecological risks.

High-resolution in-situ observations play a critical role in capturing such rapid and localized phenomena. Coastal meteo-marine platforms, such as the utilized Fortunae buoy and Senigallia Meda station, offer a cost-effective, automated system capable of operating under harsh conditions and at high sampling rates in complex nearshore zones.

We therefore recommend a collective effort that incorporates such observational datasets into Copernicus and broader databases. Such an integration would bridge existing gaps in coastal coverage, complementing the basin-wide scope of satellite and reanalysis products with finer-scale, real-time data. As climate-induced extremes accelerate, safeguarding ecological and societal resilience in shallow semi-enclosed seas, such as the Adriatic Sea, requires sustained, integrated, and transdisciplinary efforts anchored in robust, real-time observations and openly accessible data.

Appendix A

A1 Overview of the meteo-marine buoy Fortunae

The multidisciplinary buoy is equipped with a series of sensors that measure atmospheric parameters at 2 m height; wave parameters at sea level; as well as current, thermohaline and biochemical parameters at ∼ 2 m depth below the sea level.

A1.1 Atmospheric conditions at 2 m height

More precisely, the following atmospheric parameters are systematically sampled every 30 min at 2 m above sea level using a GILL MaxiMet GMX500 Compact Weather Station equipped with GPS (https://gillinstruments.com/, last access: 10 May 2026; Table A1).

Table A1 Overview of the observed atmospheric parameters with the corresponding technical details and characteristics. Additional information are available under the feature specification of the GMX500 Model: https://gillinstruments.com/wp-content/uploads/2022/08/1957-008-Maximet-gmx500-Iss-9.pdf, last access: 10 May 2026.

A1.2 Wave conditions at sea level

Note that wave measurements are available from the platform but were not utilized in the present study. The following wave parameters are systematically sampled every 60 min at sea level utilizing a Brizo-X Directional Wave Height Sensor with an integrated Global Navigation Satellite System (GNSS) from Xeos Technologies (https://xeostech.com/brizo-x, last access: 14 May 2026; Table A2).

Table A2 Overview of the observed wave characteristics with the corresponding technical details. Wave height range information stem from Wang et al. (2016). Additional information is available under the feature specification of the Brizo-X Model: https://xeostech.com/sites/default/files/2019-08/BrizoX%20Brochure%20WEB.pdf, last access: 14 May 2026.

The Brizo-X Directional Wave Height Sensor utilizes standard spectral analysis techniques and formulas to calculate wave parameters for waves of periods from 1.6 s to 33 s (Earle, 1996). The energy density spectrum of a sea state is generally designated by E(f) and the total energy is given by

$$\begin{array}{}\text{(A1)}& \stackrel{\mathrm{‾}}{E}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}E\left(f\right)\mathrm{d}f\end{array}$$

The wave frequency spectrum can, therefore, be determined from a wave record η(t) by using a Fourier transform as follows: The wave energy averaged over a period $[-T/\mathrm{2}<(t-t_{\mathrm{0}})<T/\mathrm{2}]$ is given by

$$\begin{array}{}\text{(A2)}& \stackrel{\mathrm{‾}}{E}={\displaystyle \frac{g\mathit{\rho }}{T}}\underset{-T/\mathrm{2}}{\overset{T/\mathrm{2}}{\int }}{\left(\mathit{\eta }\left(t-t_{\mathrm{0}}\right)-\langle \mathit{\eta }\rangle \right)}^{\mathrm{2}}\mathrm{d}t\end{array}$$

where 〈η〉 is the mean value.

Inserting in this expression the Fourier development of $\left(\mathit{\eta }\right(t-t_{\mathrm{0}})-\langle \mathit{\eta }\rangle )$, gives

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& \stackrel{\mathrm{‾}}{E}=\sum _{k=\mathrm{1}}^{\mathrm{\infty }}E\left(f_k\right)\mathrm{\Delta }f,f_k=k\mathrm{\Delta }f,\mathrm{\Delta }f={\displaystyle \frac{\mathrm{1}}{T}},\\ & E\left(f_k\right)\mathrm{\Delta }f={\displaystyle \frac{g\mathit{\rho }}{\mathrm{8}}}{\left|H_k\right|}^{\mathrm{2}},\\ & H_k={\displaystyle \frac{\mathrm{4}}{T}}\underset{-T/\mathrm{2}}{\overset{T/\mathrm{2}}{\int }}\left(\mathit{\eta }\left(t-t_{\mathrm{0}}\right)-\langle \mathit{\eta }\rangle \right)e^{-\mathrm{2}i\mathit{\pi }{f}_{k}t}\mathrm{d}t\end{array}\end{array}$$

where g is the gravitational acceleration, ρ is the seawater density.

A1.3 Thermohaline and biochemical conditions at 2 m depth

Additionally, the following oceanographic parameters are systematically sampled every 30 min at ∼ 2 m depth below the sea level by a submerged MicroCAT SBE-37 SIP (Serial Interface Pumped) CTD (Conductivity-Temperature-Depth) and a WET Labs ECO-FLNTU fluorescence sensor with Bio-Wiper option from Sea-Bird (https://www.seabird.com/, last access: 14 May 2026; Table A3).

Table A3 Overview of the observed oceanographic parameters^∗ with the corresponding technical details and characteristics. Additional information is available in the SBE feature specifications: https://www.seabird.com/products/sbe-37-microcat (last access: 14 May 2026).

^∗ Though not included in this study, additional parameters are available from an SBE optical oxygen sensor and a Nortek Aquadopp current meter from 14 July 2023 onwards.

The TEOS-10 GSW Matlab routines (version 3.05.5) were applied to derive Absolute Salinity (S_A) (mass fraction of dissolved material, in g kg⁻¹) and of Conservative Temperature (CT, in °C) from of Practical Salinity (S, unitless) and Potential Temperature (θ, in °C), respectively. According to McDougall and Barker (2011) these properties are calculated as follows:

$$\begin{array}{}\text{(A4)}& S_{\mathrm{A}}=\left({\displaystyle \frac{\mathrm{35.16504}}{\mathrm{35}}}\right)S+\mathit{\delta }{S}_{\mathrm{A}}(x,y,p)\end{array}$$

where $\mathit{\delta }{S}_{\mathrm{A}}=S_{\mathrm{A}}-S_{\mathrm{R}}$ is the difference between the Absolute Salinity and Reference Salinity (S_R) based on the location, with x and y being the longitude and latitude in decimal degrees and p the pressure in dbar.

$$\begin{array}{}\text{(A5)}& \mathrm{CT}={\displaystyle \frac{h^{\mathrm{0}}}{C_p^{\mathrm{0}}}}\end{array}$$

where $C_p^{\mathrm{0}}=\mathrm{3991.86795711963}$ J kg⁻¹ K⁻¹ is a reference value of the specific heat capacity, chosen to be as close as possible to the spatial average of the heat capacity over the ocean surface and h⁰ is the potential enthalpy which corresponds to the enthalpy (h, in J kg⁻¹) at a reference pressure:

$$\begin{array}{}\text{(A6)}& h^{\mathrm{0}}\left(S,\mathit{\theta },p_{\mathrm{r}}\right)=h(S,\mathit{\theta },p)-\underset{p_{\mathrm{r}}}{\overset{p}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }(S,\mathit{\theta },p)}}\mathrm{d}p\end{array}$$

Readers interested in the full observational capabilities of the buoy are referred to Table 1 and Penna et al. (2025a).

A2 Freshwater flow rate calculation

Following the salt mass balance equation (Eq. 1) and based on the in-situ data from the Fortunae buoy (Product Ref. No. 1) and the estimated impacted area from the satellite observations (Product Ref. No. 5 and 6), the freshwater input leading to the first extreme freshening can be calculated as follows.

Taking into account that

• ρ₁=1026.212 kg m⁻³ is the initial seawater density measured by the buoy;

• ρ₂=1015.725 kg m⁻³ is the final seawater density measured by the buoy;

• S₁=36.194 g kg⁻¹ is the initial seawater salinity measured by the buoy;

• S₂=23.424 g kg⁻¹ is the final seawater salinity measured by the buoy;

• $V_{\mathrm{1}}=\mathrm{22.5}\times {\mathrm{10}}^{\mathrm{3}}$ m³ is the initial volume of seawater confined by the distance from the shore (15 × 10³ m), the width (1 m) and depth (1.5 m) considering the position of the buoy sensors;

• V₂ is the final volume.

Figure A1 Indication of the considered volume of seawater.

Note that the control volume is defined by an offshore distance of 15 × 10³ m, consistent with the observed plume extent, a depth of 1.5 m corresponding to the buoy sensor level, and an alongshore width of 1 m. This width is intentionally chosen to express the freshwater requirement per meter of coastline. The resulting freshwater rate therefore represents the required freshwater inflow per meter of coastline near the buoy location, rather than the total discharge integrated along the entire coastal segment.

Equation ρ₁S₁V₁=ρ₂S₂V₂ can be solved for $V_{\mathrm{2}}=V_{\mathrm{1}}{\mathit{\rho }}_{\mathrm{1}}{S}_{\mathrm{1}}/{\mathit{\rho }}_{\mathrm{2}}{S}_{\mathrm{2}}=\mathrm{35.762}\times {\mathrm{10}}^{\mathrm{3}}$ m³ and the corresponding volume change $\mathrm{\Delta }V=V_{\mathrm{2}}-V_{\mathrm{1}}=\mathrm{13.262}\times {\mathrm{10}}^{\mathrm{3}}$ m³ provides the freshwater volume and the flow rate (Q_f) required to achieve the average salinity change of 0.982 g kg⁻¹ d⁻¹ (from 36.194 to 23.424 g kg⁻¹) over the first specified 13 d period as follows: $Q_{\mathrm{f}}=\mathrm{\Delta }V/\mathrm{\Delta }t=\mathrm{1020.222}$ m³ d⁻¹ or 42.509 m³ h⁻¹.

A3 Statistical analysis

Following Emery and Thomson (1998), the standard deviation (σ) and the corresponding standard error (S_ε) were calculated from the variance (σ²) as follows:

$$\begin{array}{}\text{(A7)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}\sum _{i=\mathrm{1}}^N{\left(x_i-\stackrel{\mathrm{‾}}{x}\right)}^{\mathrm{2}}\text{(A8)}& {\displaystyle }{S}_{\mathit{\epsilon }}={\displaystyle \frac{\mathit{\sigma }}{\sqrt{N}}}\end{array}$$

where x_i is the ith data point in the dataset and N is the number of data points.

Correspondingly, the Root Mean Square Difference (RMSD) and bias of the different datasets (x,y) were calculated based on the following equations:

$$\begin{array}{}\text{(A9)}& {\displaystyle }\mathrm{RMSD}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{\left(y_i-x_i\right)}^{\mathrm{2}}}{N}}}\text{(A10)}& {\displaystyle }\mathrm{bias}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{\left(y_i-x_i\right)}^{.}}{N}}\end{array}$$

In addition, the cross-covariance (C_xy) and cross-correlation (r_xy) function was used to assess the relationship between two time series (Emery and Thomson, 1998). Assuming that x and y are the relative variables, the function was calculated as follows:

$$\begin{array}{}\text{(A11)}& C_{xy}\left(\mathit{\tau }\right)={\displaystyle \frac{\mathrm{1}}{N-k}}\sum _{i=\mathrm{1}}^{N-k}\left[y_i-\stackrel{\mathrm{‾}}{y}\right]\left[x_{i+k}-\stackrel{\mathrm{‾}}{x}\right]\end{array}$$

where τ=kΔt is the lag time for k sampling time increments of duration Δt

$$\begin{array}{}\text{(A12)}& r_{xy}\left(\mathit{\tau }\right)={\displaystyle \frac{C_{xy}\left(\mathit{\tau }\right)}{{\mathit{\sigma }}_x{\mathit{\sigma }}_y}}\end{array}$$

where σ_x and σ_y are the standard deviations for each time series.

A4 Summary of the freshening events

Table A4 Overview of the 5 defined freshening events, representing the start and end date, the duration (in days) and the highest and lowest measured absolute salinity during each corresponding event (in g kg⁻¹).

Table A5 Overview of the most pronounced freshening period during each event, representing the duration (in days) during which the largest total decrease in surface salinity (in g kg⁻¹) was observed. The corresponding maximum and minimum salinity values and the respective timestamps are reported.

A5 Lead-lag correlations between accumulated precipitation and river discharge responses

Figure A2 Lead-lag correlation plots between accumulated precipitation over the Emilia-Romagna region and the discharge of 20 rivers during the first major flooding periods that lead to the previously-defined freshening event I. Each subplot shows the Pearson correlation coefficient (r) as a function of lag time (in hours), where negative lags indicate that precipitation leads river discharge. The maximum correlation point for each river is marked with a red dot. Corresponding correlation values and lags are annotated in bold if statistically significant (p < 0.05), and in italics otherwise.

Figure A3 Same as Fig. A2 but for the second event.

Figure A4 Same as Fig. A2 but for the third event.

Figure A5 Same as Fig. A2 but for the fourth event.

Figure A6 Same as Fig. A2 but for the fifth event.