Analyzing trend and heatwaves of 15 Years of Sea Surface Temperature Variations along the Italian Adriatic Coast

Abstract

Water temperature is a vital parameter impacting the growth and survival of aquatic life. Using satellite-derived infrared data, this study analysed the trend of sea surface temperature (SST) from 2008 to 2022 of the Adriatic coastal waters of Italian regions. The “Mediterranean Sea High Resolution and Ultra High Resolution Sea Surface Temperature Analysis" product collected from the Copernicus Marine Service of European Copernicus programme was used, as a good compromise among spatial accuracy, temporal frequency and coverage. SST were derived in 176 locations, placed in the Adriatic Sea from the southern limit of the lagoon of Venice (Veneto) to Santa Maria di Leuca (LE), at a distance from the coast between 500 m and 5000 m (0.3 - 2.7 nautical miles).

Time series analysis was applied to average value of daily SST calculated from the selected spatial locations to identify the additive model components: trend, seasonality and random effects. The trend component was isolated and assessed using a linear regression model to determine its significance and magnitude. A 0.010 °C/year increase in SST was observed. Additionally, marine heatwaves and cold spells were consistently registered throughout the entire observation period, with a north-south gradient in intensity. 

Introduction

In recent decades, the Earth system has accumulated a significant amount of heat (Collins et al. in the IPCC report of 2019). The emitted greenhouse gases, primarily carbon dioxide (CO₂), methane (CH₄), and nitrous oxides (NOₓ) prevent the diffusion of the heat into the space, while trapping it into the atmosphere. Approximately 90% of this extra heat is absorbed by the oceans, making them a massive heat reservoir. Recent years registered the highest values of ocean temperatures ever recorded (World Meteorological Organization 2022). This ocean warming affects the water cycle as well as marine ecosystems, accelerating evaporation, altering the global precipitation patterns, contributing to sea level rise (Swapna et al. 2020). Ecologically, rising ocean temperatures alter marine habitats, threatening biodiversity (Garrabou et al. 2022, Ratnarajah et al. 2023, Smith et al. 2023).

For example, climate change-induced temperature increases can complexly alter fish interactions, affecting the dynamics and stability of marine communities (Masayuki et al. 2023). Additionally, it can decrease benthic remineralization, thereby reducing the carbon export flux from the ocean surface and diminishing the size and extent of biomass (Jones et al. 2014).

Globally, there has been an accelerated warming trend in recent decades (Bâki Iz 2018, Garcia-Soto et al. 2021), with projections indicating further increases in the future (Ruela et al. 2020). Consequently, there have been changes in air-sea humidity and heat exchanges, directly impacting the global water cycle according to the "dry gets drier and wet gets wetter" paradigm (Durack et al. 2012, Yu et al. 2020).

In aquaculture, water temperature is one of the most critical environmental factors, as it directly affects the physiology, growth, metabolism, and survival of aquatic organisms, as fishes (Islam et al. 2022) and molluscs (Kroeker et al. 2014, Des et al. 2020). Elevated water temperatures affect the ability of mussels to maintain stable cell function (Anestis et al. 2010, Des et al. 2020), making it a critical factor influencing survival and mortality, particularly when prolonged exposure is observed (Masanja et al. 2023, Brachetti et al. 2024). High temperatures can reduce mussel reproduction by affecting gonad development and gamete quality, causing early or incomplete spawning and lowering reproduction efficiency (Ceccherelli and Rossi 1984, Masanja et al. 2023, Brachetti et al. 2024). These extreme environmental conditions can also increase metabolism and oxygen use in mussels: without enough energy reserves, growth slows, leading to smaller, less robust individuals (Ceccherelli and Rossi 1984, Des et al. 2020). At elevated temperatures, mussels exhibit a decline in filtration capacity, which directly impacts their ability to feed and sustain metabolic processes. Prolonged exposure to high temperatures often results in extended valve closure, reducing food intake and leading to energy deficits (Ceccherelli and Rossi 1984, Des et al. 2020). This behavior is likely a protective response to thermal stress but comes at a cost, as it limits the mussels' ability to filter suspended particles, obtain nutrients, and clear waste from their system. Over time, reduced feeding efficiency under such conditions can impair growth, reproduction, and overall resilience, making mussels more vulnerable to additional environmental stressors. In addition, water temperature can indirectly influence the health and productivity of aquaculture systems by accelerating the growth and distribution of pathogens (disease-causing organisms) (Garrabou et al. 2022, Masanja et al. 2023).

Thus, temperature was defined an Essential Climate Variable by the Global Observing Systems Information (GOSIC) and Essential Ocean Variable (Bojinski et al. 2014). These designations give the temperature a key role as effective indicator for monitoring and understanding SST variations, in the development of climatic models, in the assessment of ecosystem services (Filipponi et al. 2017).

On a regional scale, the Mediterranean SST has also experienced significant warming, as shown in recent studies using satellite data (Pastor et al. 2017, Pastor et al. 2020), estimated to be 3.7 times higher than the global ocean warming trend (Pisano et al. 2020). In the Mediterranean region, this warming has been accompanied by a consistent increase in sea level and salinity (Menna et al. 2022).

Nowadays, satellite Earth observation data are considered the premier source for environmental and climatic information, due to their ability to regularly capture high-resolution spatial data over both terrestrial and marine environments.

There are satellites constellations whose spatial resolution, revisit frequency and available sensors have been specifically designed to monitor oceans and water bodies data. SST can be estimated from infrared radiometers onboard satellites, typically operating at wavelengths 3.7-12 µm, and represents the skin temperature, which corresponds to the upper 10-20 micrometers of the water column. Measurements of skin SST are subject to a large potential diurnal cycle, while the foundation SST is defined as the SST free of diurnal temperature variability. Considering that night-time SST is more stable and closer to the foundation SST, many SST datasets are based on measurements by infrared radiometers onboard satellite platforms acquired at night-time.

SST is estimated by sensors onboard various satellites, with varying spatial resolutions (from one to a few kilometres) and hourly to daily revisit frequency at mid latitudes. Among the most important operational polar orbiting satellite platforms there are: NASA's Aqua and Terra, whose Moderate Resolution Imaging Spectroradiometer (MODIS) sensor provides SST data at spatial resolutions of 1 km, and daily revisit frequency (https://neo.gsfc.nasa.gov/view.php?%20datasetId=MYD28M); the Copernicus Sentinel-3 constellation, equipped with the Sea and Land Surface Temperature Radiometer (SLSTR) sensor, to provide SST dataset with the spatial resolution of 1 km and a daily temporal resolution; NOAA satellites equipped with the Advanced Very High Resolution Radiometer (AVHRR) sensor that has been providing global SST data with daily revisit frequency and a spatial resolution of 4 km for more than three decades. Existing services provide regular and systematic reference information on the state (physical and biogeochemical), variability, and dynamics of the global oceans and the European seas, following a free and open-access data distribution policy. For example, Copernicus Marine Service (or Copernicus Marine Environment Monitoring Service - CMEMS, https://marine.copernicus.eu/it), NASA Physical Oceanography Distributed Active Archive Center (PO.DAAC, https://podaac.jpl.nasa.gov) or Google Earth Engine Catalog (https://developers.google.com/earth-engine/datasets/catalog) provide access to EO-derived SST datasets. These datasets are distributed as regular gridded data at the global or local level. They can represent the combination of observations different acquired by various satellites and in some cases are the result of gap-filling procedures to interpolate missing SST values, that are the caused by the presence of cloud pixels in source satellite acquisitions.

Time series analysis is a statistical method used for the analysis of data collected over time. This method allows to identify, extract and forecast patterns and trends within temporal datasets, as it allows to enhance the key seasonal and inter-annual patterns observed in biological and oceanographic data (Vantrepotte and Mélin 2010, Vantrepotte and Mélin 2011, Costa Goela et al. 2016). These studies provide indicators about long-term changes in natural conditions, such as climate change, which is why such indicators are advised or even mandated for coastal water monitoring programmes (e.g. Water Framework Directive, 2000/60/EC) (Vantrepotte and Mélin 2010).

Studies on SST variations could be integrated with those of temperature extreme events as the marine heatwaves (both hot and cold), as they are based on the same dataset of SST. Marine heatwaves (MHWs) and Marine Cold Spells (MCSs) are temperature anomalies that persist for at least five consecutive days in a location (Oliver et al. 2021), whose effects are relevant on marine ecosystems (Garrabou et al. 2022), altering the presence of food web, ecosystem dynamics, biodiversity, water quality, affecting the economy and human activities linked to the sea. MHWs are caused by “anomalies” in the local interactions between the atmosphere and the sea, sea currents, atmospheric events, cloud coverage. Climatic anomalies, global warming, changes in wind patterns and their weakness, reduction in cloud cover can collectively increase the frequency and intensity of anomalously high sea temperatures (Sen Gupta et al. 2020).

Aim of this study was to describe and measure the variations (trend, marine heatwaves) in SST over the past 15 years, along the Italian coasts of the Adriatic Sea, focussing on the zones dedicated to aquaculture, to verify and measure a potential increase in surface temperature. Although water depth, mixing processes, and thermal stratification are conditions affecting the water temperature along the water column, consistent and prolonged sea surface temperature increases, it is likely to rise the subsurface temperature.

Materials and methods

Study area

The study area encompasses the Adriatic Sea from the southern limit of the Venice Lagoon (Venice province) to Santa Maria di Leuca (Lecce province): this area considers six administrative regions, namely Veneto, Emilia Romagna, Marche, Abruzzo, Molise, Puglia going from north to south (Figure 1).

Figure. 1. Study area and spatial distribution of study locations.

Along the coast, regular sampling on water and molluscs are performed by local Veterinary authorities for health monitoring and marketing of live bivalve molluscs intended for human consumption, whether bred or present in natural beds, required by EU Regulation 2019/627 (European Commission (EC) 2019). We choose the same locations of these inspections as representative of the sea surface conditions of the sea coastal area facing the selected administrative regions.

A number of 141 locations were considered in this study, having a distance from the shoreline between 500 m (0.27 nautical miles) and 5 km (2.70 nautical miles). For the Puglia region, only 7 locations were available, leaving a significantly data gap along roughly 250 km of coastline. To maintain a consistent spatial distribution along the whole study area, 35 more locations were added, latitudinally placed at a constant distance of 10 km from each other. Their distance from the coastline was a value randomly selected from among the distances of the other locations.

Sea Surface Temperature dataset

Among the available Earth observation SST datasets, the product named "Mediterranean Sea High Resolution and Ultra High Resolution Sea Surface Temperature Analysis" (SST_MED_SST_L4_NRT_OBSERVATIONS_010_004) (Buongiorno Nardelli et al. 2013) from CMEMS catalogue was identified as reference product for the data analysis. It represents the state of art SST EO product over the Mediterranean Sea, considering the number of combined satellite sensors SST estimates, the use of gap-filling methods to interpolate missing pixels, the spatial resolution about 1 km and the daily temporal resolution. The level 4 (L4) SST dataset generated from night-time satellite observations and representing the foundation SST over the Mediterranean basin were collected from CMEMS (https://data.marine.copernicus.eu/product/SST_MED_SST_L4_NRT_OBSERVATIONS_010_004/services last access on 27/06/2023), in the NetCDF-3 format.

The L4 SST dataset, is produced on a daily basis and reprocessed after one month to ensure the highest quality in the archived datasets, following the procedure developed by National Research Council-Satellite Oceanography Group (CNR-GOS) (Buongiorno Nardelli et al. 2013). These SST products are based on night-time observations collected by infrared sensors mounted on various satellite platforms. The CNR-GOS processing chain includes several modules, from data extraction and preliminary quality control, to the removal of blurred pixels and the collation/merge of satellite images. Finally, a two-step algorithm enables interpolation of SST data at high (HR 0.0625°) and ultra-high (UHR 0.01°) spatial resolution by applying statistical techniques (Buongiorno Nardelli et al. 2013), providing gap-free maps ready for use.

CMEMS data collection was done through the Copernicus Marine Toolbox Command Line Interface (CLI), that allows to download the products of interest and to optionally define a temporal subset, a geographic subset, and the variables of interest. The present study considers 15 complete years, from January 1, 2008 (first available observation date in the repository CMEMS SST_MED_SST_L4_NRT_OBSERVATIONS_010_004 product) to December 31, 2022. Downloaded dataset has a spatial resolution of 0.01° × 0.01° in the WGS84 coordinate reference system (EPSG: 4326), and a daily temporal resolution, namely each collected file represents a single daily timestep. The values of the downloaded SST dataset were converted from Kelvin to degrees Celsius (°C) for the entire 15-year period considered.

Sea Surface Temperature data analysis

From daily TIF files, SST temperature values and corresponding dates were extracted for each study locations. The daily data were exported in Text Document (.txt) format for each year and then aggregated into a single Excel file (Figure 2). The resulting dataset for subsequent statistical analysis contained information on: study location code, Region, longitude, latitude, SST and date.

Before proceeding to time series analysis, the observed SST data were evaluated through graphical representations, such as boxplots and heatmaps, and tabular representations showing median and quartiles, through different daily and monthly grouping levels, both overall (for the entire Adriatic Sea study area) and for each region.

These statistical analyses were conducted using R version 4.3.2 (R Core Team, 2023). The packages ‘terra’ (v1.7-3; Hijmans 2023), ‘ncdf4’ (v1.21; Pierce 2023) and ‘raster’ (v3.6-14; Hijmans 2023) were used to open, read, and manipulate the data in raster format. Heatmaps were created with the ‘oce’ package (v1.8-2; Kelley and Richards 2023). Finally, ‘ggplot2’ (v3.5.1; Wickham 2016) was used to generate graphical displays.

Figure. 2. Schematic representation of the entire data processing workflow, starting with downloading and processing satellite raster data, to study locations extraction, performing statistical analysis, and time series model development.

Time series analysis

Time series {Y_t}^N_t=1 is a collection of observations indexed by time t. In this study, time series were extracted spanning from January 1, 2008 to December 31, 2022, covering a time horizon of 15 years (N = 5475 daily observation) for each study location. For the leap years falling in the period under consideration (i.e. 2008, 2012, 2016 and 2020), the observations for February 29 were omitted from the analysis.

Time series modelling

Classical decomposition methods are prevalent in the marine sciences literature (e.g., Vantrepotte and Mélin 2010, Mélin et al. 2011, Loisel et al. 2014, Costa Goela et al. 2016). These methods characterize a time series using distinct components: a trend component (T_t), representing the long-term direction over the time span; a seasonal component (S_t), reflecting repetitive patterns over time; and a random component (R_t), capturing unexplained variations not attributed to trend or seasonality. An additive model

is then employed to the decompose the daily time series into these components.

The pronounced seasonality of certain time series complicates trend measurement (DeLurgio 1998), and this is particularly true for temperatures, which inherently exhibit seasonal pattern. Hence, identifying the seasonal pattern (S_t) is crucial to eliminate seasonal variations without affecting the trend. This process yields a new time series,

Through this analysis, all components (trend, seasonality, and random) were decomposed, both for the entire geographical area of the Adriatic Sea and for individual regions falling within.

Trend analysis

The classical approach for trend analysis is to estimate a linear model for seasonally adjusted time series (Y_t^*) (Lima and Wethey 2012, Astor et al. 2013, Nicastro et al. 2013), where historical data have stable behavior over the time horizon. The linear model, including t = 1,..., N as the predictor variable, is given by

where β₀ is the intercept, β₁ is the regression coefficient or slope, and ε_t is the random error.

Marine heatwaves and cold spells

Marine heatwaves (MHWs) are qualitatively defined in the literature as temperature anomalies that persist for at least 5 consecutive days at a particular location. Anomaly can be quantitatively defined as exceeding a fixed threshold (Frölicher et al. 2018), or exceeding a seasonally varying threshold (Hobday et al. 2016, 2018), although the most common definition is as exceeding the 90th percentile of values (Oliver et al. 2021, Garrabou et al. 2022). Similarly, marine cold spells (MCSs), "anomalous, discrete, and prolonged cold water events at a particular location" (Schlegel et al. 2021), were defined as anomalies below the 10th percentile of values. In our study, the anomaly in any day of the year is calculated using data from the 15 years within an 11-day window centred on the day under evaluation. When the 90th percentile is overpassed for at least 5 consecutive days, it is recorded as MHW (or MCS) event. The ‘heatwaveR’ package of the R software (v0.5.3.9001; Schlegel and Smit 2018) was used.

To statistically verify these observations, a Kruskal-Wallis test was conducted (p-value < 0.05), identifying significant differences in peak intensity between regions. A post-hoc Dunn's test with Bonferroni correction was then performed for multiple comparisons.

Furthermore, the duration of these waves was calculated and plotted alongside the maximum intensity for individual regions over the entire period under consideration.

Results

Study location characteristics

The 176 study locations are placed along the Adriatic coast, facing the Veneto, Emilia Romagna, Marche, Abruzzo, Molise and Puglia administrative regions (Figure 1). Figure 3 shows the box and whisker plots at study locations, aggregated by region, for the two main characteristics: distance to coastline and water depth. The median distance (Figure 3a) of the study locations from the coastline falls between approximately 1 and 2 kilometers (Table SM I). While statistically significant differences were found overall between regions (Kruskal-Wallis test: p-value < 0.05), Dunn's post-hoc tests with Bonferroni correction revealed only a significant difference between Abruzzo and Puglia (Bonferroni corrected: p-value < 0.0033).

Their median water depth values (Figure 3b) are around 10 m (Table SM II), and there are no statistically significant differences between regions (Kruskal-Wallis test: p-value > 0.05).

Figure. 3. Box and whisker plots depicting the distribution of distances from the coastline (a) and water depths (b) of study locations, grouped by region.

Sea Surface Temperature description

Figure 4 shows the spatial distribution of monthly SST averages, on a raster basis, across the whole Adriatic Sea, by month and based on daily data from 2008 to 2022.

Figure. 4. Monthly mean SST for the study period (2008-2022), based on daily data aggregated by month.

Figure 5 depicts the average temperature across each day during the analysed period for study locations in the Adriatic Sea, while Figure 6 illustrates these evolutions for each region. In these visualizations, cooler temperatures are depicted in deep blue hues, while warmer temperatures are represented by increasingly intense shades of red. The study area exhibits the characteristic mid-latitude pattern of SST, with lower values occurring during the winter months and peaking during the summer season.

Figure. 5. Heatmap Adriatic Sea. Mean SST values for each day of the year (shown on the vertical axis) and for each year (shown on the horizontal axis) for the Adriatic Sea area.

Figure. 6. Heatmap for specific regions. Mean SST values for each day of the year (shown on the vertical axis) and for each year (shown on the horizontal axis).

In particular, the peak of daily average of SST occurs around the beginning of August, indicating the onset of a downward phase, and the lowest point occurring towards the end of February, marking the beginning of an upward trend (figures SM 1- SM 2).

The monthly distribution of SSTs in the Adriatic Sea is shown in Figure 7, while Figure SM 3 provides the detailed view for individual regions.

Figure. 7. Box plot of the monthly distribution of SST in the 176 location for the Adriatic Sea during the study period (2008-2022).

SST time series analysis

Figures 8 and SM 4 show the time series decomposition of the average value of daily SST calculated from the selected spatial locations, for both the entire study area (Adriatic Sea) and individual regions, respectively. In the figures, moving from top to bottom, the graph shows the observed data across the study period (2008-2022), the trend component together with the regression line, the seasonal component and, lastly, the random effect. The time series decomposition employs a 12-month centered moving window. This window captures values from 6 months before and 6 months after the target month. Consequently, the first and last six months of the data series lack trend and random components by definition.

Figure. 8. Time series decomposition of the average value of daily SST calculated from the selected spatial locations for the Adriatic Sea.

The summary table (Table I) of the regression analysis reveals statistically significant positive slopes (coefficients) in all regions except Marche. These coefficients range from 0.004°C/year in Abruzzo to 0.041°C/year in Puglia. Consequently, over the 15-year observation period (2008-2022), SST increases range from 0.06°C in Abruzzo to 0.62°C in Puglia.

Table. I. Table of regression coefficients for trend analysis of SST time series data. Asterisk is reported for regions with statistically significant temperature values.

SST exhibits an overall increasing trend. Figure 9 reports the trend component for each region, along with corresponding regression lines and 95% confidence intervals.

Figure. 9. SST trend component by region, including regression lines with 95% confidence intervals. It is worth reporting that the graph shows only the trend component, not the seasonality of the SST series.

Marine heatwaves and cold spells analysis

Both MHWs and MCSs were consistently observed throughout the entire dataset for the years 2008-2022 in the study locations, both for the Adriatic Sea as a whole and for individual regions. Figure 10 shows the number of MHWs occurred in the study area by year and month, as sum of events occurred at regional level. Higher number of peaks are reported when the MHWs occurred throughout the coast (as beginning of 2014, summer 2018, autumn 2019). The distribution of events by region is reported in figure SM 5.

Figure. 10. Number of MHWs for the study area by year and month during the period 2008-2022.

The maximum intensity of MHWs and their duration (in days) were illustrated in Figure 11 for investigated regions. In May 2009, it was registered the most intense event across all regions in the study period (years 2008-2022) with more than 4°C of anomaly peak intensity; it was the most intense event for the regions Marche, Abruzzo and Puglia (figure SM 6). The MHWs duration was 9.3 days (on average across all regions). The second most intense and prolonged period with MHWs was registered in 2018 across all regions, with MHWs occurring in May (10.8 days average duration, 3.3°C average anomaly intensity), June (higher values of anomalies of 4.5°C, average duration 15.8 days), August and September (2.5°C average anomaly intensity, average duration of the events 11.25 days). May 2018 registered the highest intensity for Molise region across the study period.

14 MHWs occurred in September-October 2011, distributed across all regions, with higher intensities in northern regions and with an average duration of 10 days. In 2014, 15 MHWs occurred in January-March across all regions lasting for 11.3 days (on average); highest intensity events (anomalies of 5°C) were reported in northern regions in June, and in the month of October 2014 another wave occurred across all regions, lasting on average 7.7 days. Although with less intensity (average 2.10°C), the year 2019 registered 21 MHWs distributed across all regions from June to December, highlighting its higher temperatures than the average of the reference period; the events lasted on average 8.5 days. 21 MHWs occurred in 2022, with southernmost regions involved mainly in summer months, while northernmost regions reported MHWs with most intense anomalies late in the year; in November another wave occurred across all areas, involving all regions with an average of 2.5°C and with a duration of 20 days (average across regions, ranging 10-33), which was the longest wave in the study period.

Figure. 11. MHWs maximum intensity (°C, dark red bars, left y-axis) and length in days (blue lines, right vertical axis) for individual regions.

The number of extreme events in the northern regions (Veneto and Emilia Romagna) resulted 25% higher compared to the other regions (Figure SM 5 and 11). Furthermore, looking at the intensity of MHWs, Northern regions reported significantly higher values than the other regions, as is highlighted in Table II and supported by the application of a Kruskal-Wallis (p-value < 0.05) and Dunn's post-hoc tests in Table III.

Table. II. Summary statistics of intensity (°C) of the maximum value in MHW events for individual regions in the period 2008-2022.

Table. III. p-value of post-hoc Dunn test after Bonferroni correction (significance level: 0.0033) for intensity (°C) of the maximum value in MHW events for individual regions in the period 2008-2022.

Concerning MCSs, figure 12 illustrates the distribution of peaks intensity across the investigated regions for the 15-year period. The duration (in days) of these peak events is shown by a blue line, and the maximum intensity is indicated by light blue bars. January to May, and October, are the months where MCSs most often occurs (figure 13). The years 2010, 2014, 2021 reported the majority of number of events (26, 16, 22 respectively), although with a different distribution along the year: in 2010, cold spells occurred every month except July and August; similarly, in 2021, most of the Winter and Spring months were involved. Cold spells were concentrated from July to September in 2014.

When looking at the most intense peaks (average of 3.2-3.5°C anomaly), they were concentrated in single months, i.e. June 2009, June 2013, May 2019, when all the regions were involved. Molise region had the most intense MCS recorded in June 2009 (Figure SM 7); Marche and Abruzzo experienced their peak MCS intensity in June 2013; Veneto, Emilia Romagna and Puglia regions had their peak intensity of cold spells (MCSs) in June 2019.

Figure. 12. MCSs maximum intensity (°C, light blue bars, primary y-axis) and duration in days (blue lines, secondary vertical axis) for individual regions.

Similarly to MHWs, MCSs also show more sensitive peak intensities in northern regions than in southern regions (Table IV). Statistical analysis confirms this observation with a Kruskal-Wallis test revealing a p-value < 0.05, indicating statistically significant differences between the regions. Table V details the p-values from Dunn's post-hoc tests with Bonferroni correction for multiple comparisons. Notably, the Puglia region stands out with significantly lower values compared to all other regions, signifying less intense cold spells.

Figure 13 shows the number of MCSs occurred in the study area by year and month, as sum of events occurred at regional level. The distribution of events per region is shown in Figure SM 8.

Table. IV. Summary statistics of intensity of the MCS peaks for each region in the period 2008-2022.

Table. V. p-value of post-hoc Dunn test after Bonferroni correction (significance level: 0.0033) for intensity (°C) of the minimum value in MCS events for individual regions in the period 2008-2022.

Figure. 13. Number of peaks of MCSs for the entire study area by Year and Month during the period 2008-2022.

Discussion

In this paper, the SST along the Italian cost of the Adriatic Sea were described in 15 years, from 2008 to 2022, using satellite-derived Earth observation data. The temperature trend was analysed through time series analysis and integrated with the study of heatwaves and cold spells occurred in the same period. These different aspects of the temperatures at sea surface were evaluated in locations related to the aquaculture activities to provide insight in the knowledge/understanding of this biological communities.

The changes observed in this study, although over a relatively short period, are consistent with other research conducted in the Mediterranean region (Pastor et al. 2020, Pisano et al. 2020), highlighting a significant upward trend in SST. Additionally, prolonged increases in water temperature, referred to in the literature as marine heatwaves (MHWs), have become more evident in terms of both frequency and duration, with these events becoming more intense and widespread (Collins et al. 2019, Garrabou et al. 2022, Pastor et al. 2023, Brachetti et al. 2024).

When looking at the maps of SST averaged for the study period on a monthly basis (figure 4), they show that the Adriatic Sea follows a predictable pattern, with the warmest waters in August and the coldest in February, in accordance with the changing tendency reported in other studies (Pastor et al. 2020, 2023). Notably, nearshore temperatures fluctuate more rapidly due to shallower depths: these areas cool first in November-December and warm first in May-June, mirroring air temperature changes. This seasonal cycle influences fish migrations, driving them offshore in fall/winter and back to the coast in spring.

Along the Italian Adriatic coast, aquaculture farms and natural mollusc beds are located nearshore, approximately between 1 and 2 km from the coastline, at a median average of 10 meters depth (figure 1, figure 3). Given the relevance of the aquaculture sector in Italian economy, and considering live bivalve molluscs as “sentinel” of the quality status of the sea, either bred or present in natural beds, we focussed our analyses on these locations. Without any in-situ measurements, we characterised the evolution of surface temperature of the sea, in the areas where molluscs live.

Our analysis of the SST data in the study sites suggests a pronounced warming trend between 2008 and 2022, particularly in August (figure 5 and figure 6). Compared to previous years of the study, which exhibited alternating temperature patterns (represented by the colour changes in the graph), August has experienced a sustained shift towards higher temperatures from 2017 to 2022. September displays a similar trend, albeit with a less pronounced warming effect compared to August. This is reflected by the presence of alternating colours in the early years, followed by a consistent red colour indicating higher temperatures from 2017 onwards. When looking at individual regions, the largest SST discrepancies occur during the colder months (October to March) (figure 7 and figure SM 3). Here, northernmost regions like Veneto and Emilia Romagna exhibit temperatures that are, on average, 2-3°C lower than their southern counterparts, as following a latitude gradient (in figure 6, the blue-colored areas are with a darker colour following latitude north-south). The tendency is opposite for summer months, with differences more evident for August temperatures in southernmost regions.

We analysed also the time series in the study locations and in the fifteen years. After removing other components (seasonality, background noise), an increasing temperature trend is observed along the north-south axis (latitude) with the exception of the Veneto and Marche regions. While Figure 9 shows colder waters in Veneto and Emilia Romagna compared to other regions, Table I reveals a stronger warming effect in Veneto. The causes of this anomaly likely lie in a combination of factors: latitude, freshwater inflow from rivers, potential drought periods, seabed depth, and the influence of sea currents and coastal conformation.

For the southernmost regions (Abruzzo, Molise, and Puglia), the north-south latitudinal gradient appears to be the primary driver of the increasing trend. As shown in Figure 9, these regions have the highest temperature values, with a clear latitudinal gradient: Abruzzo (0.004) < Molise (0.019) < Puglia (0.041). Interestingly, the temperature trends for these regions overlap from 2008 to 2015, but diverge thereafter, suggesting the influence of additional factors beyond latitude.

Consistent with its central geographic position, the Marche region exhibits a temperature trend intermediate between the northern and southern regions. Notably, unlike other regions, the trend in Marche appears to be stable. This suggests that factors beyond latitude, such as coastal conformation, shoreline and seafloor morphology, and sea currents, may be exerting a significant influence. These factors could potentially have synergistic effects, leading to increased water mixing within the water column, which might be offsetting the broader effects of climate change.

Overall, the analysis of the SST data series from 2008 to 2022 revealed an annual increase in the trend, which, net of other components (seasonality, background noise), is equal to 0.010°C across all regions.

Both marine heatwaves and marine cold spells were consistently observed throughout the entire dataset for the years 2008-2022 in the study locations, both for the Adriatic Sea as a whole and for individual regions (figures 10 and 11 for MHWs, and figures 12 and 13 for MCSs). These results on MHWs are coherent with other studies (Frölicher et al. 2018, Oliver et al. 2021, Garrabou et al. 2022, Pastor et al. 2023). The most intense event across all regions occurred in May 2009, with more than 4°C anomaly in the peak intensity and an average duration of more than 9 days. Again May (but also June, August and September) of 2018 registered record values in duration (average of almost 11 days) and peak intensity (3-4°C of anomaly). While these events were globally reported along the whole area of study, others were registered locally (figure SM 5 and figure SM 6), reflecting the normal geographical differences. The number of extreme events and their intensity is higher in northern regions, and decrease going southwards (figure 10 and table II).

The MCSs normally occur during autumn and winter months in the study area, from October to March, especially in the second part of the observation period (from 2015 onwards). Notably, from 2015 onwards, MCS events are almost absent in the summer months (figure 13), coupling the results on major presence of extreme warm events. This issue is particularly concerning because, despite MCSs being less easily interpretable than MHWs, an expanding body of research indicates that low SST events can also exert significant ecological impacts and are associated with extreme weather events on land (Duchez et al. 2016).

The number of cold extreme events and their intensity is higher in northern regions, and decrease going southwards (table IV).

The use of satellite-derived SST data offers significant advantages, particularly in terms of extensive spatial and regular temporal coverage (pixels of about 1 km side and daily overpass in our study). The high detailed data (high spatio-temporal resolutions) is crucial to deeper understand the phenomena, especially when dealing with molluscs, and marine life in general, as small variations in temperature can affect the life dynamics of the animals (Dillon et al. 2016). This motivated our choice to analyze satellite data at specific locations, in order to more accurately characterize the temperature conditions at aquaculture sites on the natural seabed of the Adriatic Sea, which are primarily located in shallow waters near the coast.

Disadvantage of satellite-derived data is that temperature measured by satellite and used in this study refers to the SST, which is influenced by the atmospheric temperature and its fluctuations, solar radiation, wind, and short-term weather conditions. These factors make SST not always accurately reflecting the temperature of the underlying water column, which can differ significantly due to currents, local seabed morphology and composition influence the water mix and influence the MHWs (Schaeffer & Roughan 2017). However, in shallow waters, like our study locations, close to the coast, the variation of temperature along the water column is less pronounced and the SST can be considered, with reliable and sufficient approximation, as representative of the fish and molluscs environmental conditions.

Especially at night, may have a stronger correlation with the temperature of the entire water column due to more rapid thermal mixing.

The SST data analyzed in this study, processed using CMEMS products, from 2008 to 2022, confirm a significant increase in surface temperatures at the sampling locations along the Adriatic Sea coast. The upward trend identified in this work may pose a potential limiting factor for marine biological resources in general, and specifically for the cultivation of bivalve molluscs, directly impacting their survival and behavior.

Previous studies, utilizing climate models, have indicated that climate change may lead to a reduction in phytoplankton species richness and size, as well as the possibility of certain species outcompeting others, thereby altering the structure of the trophic network, reducing the number of potential ecological niches, and decreasing biodiversity (Henson et al. 2021, Tsirintanis et al. 2022, Ratnarajah et al. 2023).

The trend towards increasing dominance of smaller phytoplankton species has significant implications for both the ecological and biogeochemical functions of the oceans. Regions dominated by smaller phytoplankton typically support less productive food webs and sequester less organic carbon than those dominated by larger species. These shifts in species composition occur as environmental conditions exceed the tolerances of the existing phytoplankton types (Henson et al. 2021).

As a result, if primary consumers of phytoplankton, such as zooplankton or filter-feeding organisms, are unable to obtain adequate prey, significant disruptions to marine trophic networks are likely, leading to biodiversity loss, ecosystem instability, and a decrease in fishery resource biomass (Ratnarajah et al. 2023).

Moreover, climate change, characterized by strong temperature fluctuations coupled with periods of intense rainfall or drought that significantly alter sea salinity, can induce stress in bivalves, leading to changes in their physiological and behavioral responses and, in extreme cases, causing high mortality rates. The projected high temperatures and sporadic heavy rainfall in southern Europe due to climate change will have severe implications for natural and human systems and will likely be accompanied by increased ocean acidification (due to elevated CO2 levels) and decreased dissolved oxygen (due to rising temperatures). From a public health perspective, rising temperatures have also been identified as a primary trigger for the accumulation of certain biotoxins in European bivalves (Dhanji-Rapkova et al. 2023).

Continuous monitoring of sea surface temperatures is therefore essential, particularly using rapid and cost-effective methods such as satellite detection, to validate model predictions and assess the resilience of marine ecosystems (Valentini et al. 2016). Only through such efforts can we effectively implement ecological restoration policies and adapt fishing and aquaculture strategies accordingly.

In the future, we aim to utilize higher spatial resolution data over longer time periods to conduct pixel-based analyses of trends and MHWs. This approach could enable a spatial characterization of anomalies, helping to identify areas more prone to extreme events or others more suitable for aquaculture due to their stable conditions. Temporally, this would allow us to confirm and further quantify the warming trend, providing a robust foundation for resilience policies. Additionally, higher spatial resolution would enable us to better measure and characterize the effects of freshwater input, particularly in coastal regions.

A more comprehensive study should incorporate additional variables that are influenced by and interconnected with temperature variations, such as bottom depth, winds, tides, seawater density, evaporation, freshwater input, and latitude—all of which significantly affect the intensity and direction of currents and, consequently, water temperature (Lipizer et al. 2014).

Integrating temperature data with other abiotic variables derived from Earth observation (e.g., chlorophyll, suspended solids) (Ippoliti et al. 2020, Filipponi et al. 2021) and in situ measurements (e.g., pH, salinity) would further enhance our understanding of water quality.

Acknowledgements

The activities were partially funded by Abruzzo region, Servizio Sviluppo Locale ed Economia Ittica, Dipartimento Agricoltura, Italy. FEAMP 2013-2020 funds.

Conflict of interest statement

No potential conflict of interest was reported by the authors.